source: proiecte/pmake3d/make3d_original/Make3dSingleImageStanford_version0.1/third_party/missing-data/create_nullspace.m @ 37

function NULLSPACE  = create_nullspace(M,INC,r,n)

% Fit a rank r matrix to a matrix with missing data.  This uses my
% algorithm of taking the most orthogonal space to the null space generated
% by r-tuples of the columns of M.  In principle, it would be nice
% to use all r-tuples, but instead we'll pick n random ones.

%global NULLSPACEexternal

% 10 is the default for the simulated motion experiments.
% however, for the real intensity images, may want to change
% since matrix dimensions are very different, eg, 9 rows
% instead of 40 in M.

%MAXNULLSIZE = 100;  
MAXNULLSIZE = 10;
% if the nullspace gets MAXNULLSIZE times the number of columns in M, we say
% that that's enough, and stop generating elements.

Mdim = size(M);
numrows = Mdim(1);
numcols = Mdim(2);
incsum = sum(INC);
NULLSPACE = [];
current_row = 1;
for i = 1:n
  if size(NULLSPACE,2) >= numrows*MAXNULLSIZE, break, end

  col_nums = find(INC(current_row,:) ~= 0);

% first check that there actually is enough data to recover that row.
% note that although we simply choose columns at random if the check fails,
% it will be impossible to recover the row at any point in the process,
% since it is under-constrained.

  if (length(col_nums) < r)
    col_nums = diff_rand_ints([],r,1,num_cols);
  else
    col_nums = col_nums(diff_rand_ints([],r,1,length(col_nums)));
  end

  incsub = INC(:,col_nums);
  rowsum = sum(incsub')';
  fullrows = rowsum==r;
  if sum(fullrows) > r

    submatrix = M(find(fullrows),col_nums);
    subnull = nulleps(submatrix,.001,-1);
    % .1 is a total hack.  Also, I'm not sure if its good to check sum(fullrows) > r.
    if size(submatrix,1) == size(submatrix,2) + size(subnull,2)
    %% The null space of the submatrix, combined with the submatrix should be a full
    %% rank square matrix.  However, if not, it means that the submatrix was rank
    %% deficient (to within tolerance) and the null matrix is too big.  We can't know
    %% which of these components are the correct ones, and which just arise because the
    %% submatrix was rank deficient.  This test will make the whole method not work if
    %% M really does have rank less than r (in which case the whole null space is right.
      nullTEMP = zeros(numrows,size(subnull,2));
      nullTEMP(fullrows,:) = subnull;
      NULLSPACE = [NULLSPACE, nullTEMP];
    end
  end

  %% Below is a hack added to try to focus attention on the rows for which
  %% we have no information.  Half the time, it makes the current row to be
  %% one for which we don't have any information yet in the nullspace.  The
  %% other half of the time it's just incremented.  I'm not sure whether this
  %% helps significantly.
  if isempty(NULLSPACE)
    sumnonzeros = [];
    nonsamprows = [];
  else
    sumnonzeros = (sum(NULLSPACE' ~= 0)')';
    nonsamprows = find(0 == sumnonzeros);
  end
  if (rand(1) > .5) & (length(nonsamprows) > 0)
    current_row = nonsamprows(random_int(1,length(nonsamprows)));
  else
    current_row = 1+rem(current_row,numrows);
  end
end
%NULLSPACEexternal = NULLSPACE;