Dual regularization-based image resolution enhancement for asymmetric stereoscopic images

by

 

17 Apr 2012 (Updated )

Dual regularization-based image resolution enhancement for asymmetric stereoscopic images

SP_2012_SR_main
function SP_2012_SR_main

% This is a demo program of the paper J. Tian, L. Chen, and Z. Liu,
% "Dual regularization-based image resolution enhancement for asymmetric 
% stereoscopic images," Signal Processing, Vol. 92, No. 2, Feb. 2012, pp. 490-497. 
clear all; close all; clc; warning off;

img_low = double(imread('image_low.bmp'));
img_truth = double(imread('image_truth.bmp'));
img_reference = double(imread('image_reference.bmp'));

param.NumberOfFrames = 1;
param.nEnlargeRatio = 2;    
param.lambda = 0.02;
param.beta = 2*param.lambda;
param.iter = 15;
[param.nRowLR, param.nColLR, param.nChannel] = size(img_low);
param.nRowHR = param.nRowLR*param.nEnlargeRatio;
param.nColHR = param.nColLR*param.nEnlargeRatio; 
%Point spread function
Hpsf = fspecial('average', 3);               
param.step = 1;

% Initial estimation of HR
HR = func_single_image_interp(img_low,param.nEnlargeRatio,'cubic');

param.BMESIZE = 16;
param.RangeGlobal = 20;    
HR_reference = func_disparity_estimation_compensation(img_low, img_reference, param.nEnlargeRatio, param.BMESIZE, param.RangeGlobal);

iter = 1;
while iter<=param.iter
    delta_likelihood = func_SR_likelihood_term(HR, img_low, Hpsf, param);
    delta_reg = func_SR_regularization_term(HR);
    delta_view = func_SR_view_term(HR, HR_reference);

    view_thre = 30;
    mask_view = (abs(HR-HR_reference)<=view_thre);
    map_phase = func_phasecong(HR);
    map_phase = map_phase.^2;
    mask_phase = (map_phase<=mean(map_phase(:)));

    HR = HR - param.step.*(delta_likelihood + param.lambda.*mask_phase.*delta_reg + param.beta.*mask_view.*delta_view);

    fprintf('%d/%d iteration is finished\n',iter, param.iter);
    iter = iter+1;      
end

imwrite(uint8(HR),['image_rec.bmp'],'bmp');
fprintf('PSNR=%.2f dB\n',func_psnr_gray(img_truth, HR));

%--------------------------------------------------------------------------
%----------------------------Inner Function -------------------------------
%--------------------------------------------------------------------------

function result = func_single_image_interp(x, factor, mode)
%perform image interpolation for single image only
[nRowLR, nColumnLR] = size(x);
nRowHR = nRowLR * factor;
nColumnHR = nColumnLR * factor;
[X_LR, Y_LR] = meshgrid(1:factor:nColumnHR, 1:factor:nRowHR);
[X_HR, Y_HR] = meshgrid(1:nColumnHR, 1:nRowHR);
result = interp2(X_LR, Y_LR, x, X_HR, Y_HR, mode);
result = inpaint_nans(result);

%--------------------------------------------------------------------------
%----------------------------Inner Function -------------------------------
%--------------------------------------------------------------------------

function img_MC = func_disparity_estimation_compensation(img_low, img_reference, nENLARGERATIO, block_size, range_global)
range_local = 5;
img_current = imresize(img_low, nENLARGERATIO, 'bicubic');

% first step
error_min = size(img_reference,1)*size(img_reference,2)*255;
disparity_global = 0;
for i=0:2:range_global
    temp = img_reference;
    temp(:,i+1:end) = img_current(:,1:end-i);
    if sum(sum(abs(temp-img_reference)))<error_min;
        disparity_global = i;
        error_min=sum(sum(abs(temp-img_reference)));
    end
end

% second step
img_MC = img_current;
for i=1:block_size:size(img_current,1)
    for j=1:block_size:size(img_current,2)
        temp_current = img_current(i:i+block_size-1,j:j+block_size-1);
        temp_MC = temp_current;
        error_min = block_size*block_size*255;
        for k=-range_local:range_local
            if (j+disparity_global+k>=1) && (j+disparity_global+k<=size(img_current,2)) && (j+block_size-1+disparity_global+k>=1) & (j+block_size-1+disparity_global+k<=size(img_current,2))
                temp_reference = img_current(i:i+block_size-1,j+disparity_global+k:j+block_size-1+disparity_global+k);

                if sum(sum(abs(temp_current-temp_reference)))<error_min;
                    error_min=sum(sum(abs(temp_current-temp_reference)));
                    temp_MC = temp_reference;
                end
            end
        end
        img_MC(i:i+block_size-1,j:j+block_size-1) = temp_MC;
    end
end


%--------------------------------------------------------------------------
%----------------------------Inner Function -------------------------------
%--------------------------------------------------------------------------

function G = func_SR_likelihood_term(HR, LR, Hpsf, param)

nEnlargeRatio = param.nEnlargeRatio;
Zn = imfilter(HR, Hpsf, 'symmetric');
HR_shift_down = Zn(1:nEnlargeRatio:end,1:nEnlargeRatio:end);
             
Gsign = sign(HR_shift_down-LR);
Gsign_HR = func_single_image_interp(Gsign,param.nEnlargeRatio,'cubic');
G = imfilter(Gsign_HR, flipud(fliplr(Hpsf)), 'symmetric');

G = sum(G, 3);

%--------------------------------------------------------------------------
%----------------------------Inner Function -------------------------------
%--------------------------------------------------------------------------

function G = func_SR_regularization_term(Xn)

G = zeros(size(Xn));

Xpad = padarray(Xn, [1 1], 'symmetric');

for l=-1:1
    for m=-1:1
        Xshift = Xpad(1+1-l:end-1-l, 1+1-m:end-1-m);
        Xsign = sign(Xn-Xshift);
        Xsignpad = padarray(Xsign, [1 1], 0);
        Xshift = Xsignpad(1+1+l:end-1+l, 1+1+m:end-1+m);
        G = G + Xsign-Xshift;
    end
end

%--------------------------------------------------------------------------
%----------------------------Inner Function -------------------------------
%--------------------------------------------------------------------------

function Gsign = func_SR_view_term(HR, HR_reference)
             
Gsign = sign(HR-HR_reference);

%--------------------------------------------------------------------------
%----------------------------Inner Function -------------------------------
%--------------------------------------------------------------------------

function result = func_psnr_gray(f, g)

f = double(f);
g = double(g);
Q=255; MSE=0;
[M,N]=size(f);
h = f - g;
MSE = sum(sum(h.*h));
MSE=MSE/M/N;
result=10*log10(Q*Q/MSE);


%-------------------------------------------------------------------------
%------------------------------Inner Function ----------------------------
%-------------------------------------------------------------------------
% Calculate local phase coherence at each pixel position
function phaseCongruency = func_phasecong(im)                    

sze = size(im);

if nargin < 2
    nscale          = 3;     % Number of wavelet scales.
end
if nargin < 3
    norient         = 6;     % Number of filter orientations.
end
if nargin < 4
    minWaveLength   = 3;     % Wavelength of smallest scale filter.
end
if nargin < 5
    mult            = 2;     % Scaling factor between successive filters.
end
if nargin < 6
    sigmaOnf        = 0.55;  % Ratio of the standard deviation of the
                             % Gaussian describing the log Gabor filter's transfer function 
			     % in the frequency domain to the filter center frequency.
end
if nargin < 7
    dThetaOnSigma   = 1.7;   % Ratio of angular interval between filter orientations
			     % and the standard deviation of the angular Gaussian
			     % function used to construct filters in the
                             % freq. plane.
end
if nargin < 8
    k               = 3.0;   % No of standard deviations of the noise energy beyond the
			     % mean at which we set the noise threshold point.
			     % standard deviation to its maximum effect
                             % on Energy.
end
if nargin < 9
    cutOff          = 0.4;   % The fractional measure of frequency spread
                             % below which phase congruency values get penalized.
end
   
g               = 10;    % Controls the sharpness of the transition in the sigmoid
                         % function used to weight phase congruency for frequency
                         % spread.
epsilon         = .0001; % Used to prevent division by zero.


thetaSigma = pi/norient/dThetaOnSigma;  % Calculate the standard deviation of the
                                        % angular Gaussian function used to
                                        % construct filters in the freq. plane.

imagefft = fft2(im);                    % Fourier transform of image
sze = size(imagefft);
rows = sze(1);
cols = sze(2);
zero = zeros(sze);

totalEnergy = zero;                     % Matrix for accumulating weighted phase 
                                        % congruency values (energy).
totalSumAn  = zero;                     % Matrix for accumulating filter response
                                        % amplitude values.
orientation = zero;                     % Matrix storing orientation with greatest
                                        % energy for each pixel.
estMeanE2n = [];

% Pre-compute some stuff to speed up filter construction

x = ones(rows,1) * (-cols/2 : (cols/2 - 1))/(cols/2);  
y = (-rows/2 : (rows/2 - 1))' * ones(1,cols)/(rows/2);
radius = sqrt(x.^2 + y.^2);       % Matrix values contain *normalised* radius from centre.
radius(round(rows/2+1),round(cols/2+1)) = 1; % Get rid of the 0 radius value in the middle 
                                             % so that taking the log of the radius will 
                                             % not cause trouble.
theta = atan2(-y,x);              % Matrix values contain polar angle.
                                  % (note -ve y is used to give +ve
                                  % anti-clockwise angles)
sintheta = sin(theta);
costheta = cos(theta);
clear x; clear y; clear theta;      % save a little memory

% The main loop...

for o = 1:norient,                   % For each orientation.
%   disp(['Processing orientation ' num2str(o)]);
  angl = (o-1)*pi/norient;           % Calculate filter angle.
  wavelength = minWaveLength;        % Initialize filter wavelength.
  sumE_ThisOrient   = zero;          % Initialize accumulator matrices.
  sumO_ThisOrient   = zero;       
  sumAn_ThisOrient  = zero;      
  Energy_ThisOrient = zero;      
  EOArray = [];          % Array of complex convolution images - one for each scale.
  ifftFilterArray = [];  % Array of inverse FFTs of filters

  ds = sintheta * cos(angl) - costheta * sin(angl); % Difference in sine.
  dc = costheta * cos(angl) + sintheta * sin(angl); % Difference in cosine.
  dtheta = abs(atan2(ds,dc));                           % Absolute angular distance.
  spread = exp((-dtheta.^2) / (2 * thetaSigma^2));      % Calculate the angular filter component.

  for s = 1:nscale,                  % For each scale.

    % Construct the filter - first calculate the radial filter component.
    fo = 1.0/wavelength;                  % Centre frequency of filter.
    rfo = fo/0.5;                         % Normalised radius from centre of frequency plane 
                                          % corresponding to fo.
    logGabor = exp((-(log(radius/rfo)).^2) / (2 * log(sigmaOnf)^2));  
    logGabor(round(rows/2+1),round(cols/2+1)) = 0; % Set the value at the center of the filter
                                                   % back to zero (undo the radius fudge).

    filter = logGabor .* spread;          % Multiply by the angular spread to get the filter.
    filter = fftshift(filter);            % Swap quadrants to move zero frequency 
                                          % to the corners.

    ifftFilt = real(ifft2(filter))*sqrt(rows*cols);  % Note rescaling to match power
    ifftFilterArray = [ifftFilterArray ifftFilt];    % record ifft2 of filter

    % Convolve image with even and odd filters returning the result in EO
    EOfft = imagefft .* filter;           % Do the convolution.
    EO = ifft2(EOfft);                    % Back transform.

    EOArray = [EOArray, EO];              % Record convolution result
    An = abs(EO);                         % Amplitude of even & odd filter response.

    sumAn_ThisOrient = sumAn_ThisOrient + An;     % Sum of amplitude responses.
    sumE_ThisOrient = sumE_ThisOrient + real(EO); % Sum of even filter convolution results.
    sumO_ThisOrient = sumO_ThisOrient + imag(EO); % Sum of odd filter convolution results.

    if s == 1                             % Record the maximum An over all scales
      maxAn = An;
    else
      maxAn = max(maxAn, An);
    end
    
    if s==1
      EM_n = sum(sum(filter.^2));           % Record mean squared filter value at smallest
    end                                     % scale. This is used for noise estimation.

    wavelength = wavelength * mult;         % Finally calculate Wavelength of next filter
  end                                       % ... and process the next scale

  % Get weighted mean filter response vector, this gives the weighted mean phase angle.

  XEnergy = sqrt(sumE_ThisOrient.^2 + sumO_ThisOrient.^2) + epsilon;   
  MeanE = sumE_ThisOrient ./ XEnergy; 
  MeanO = sumO_ThisOrient ./ XEnergy; 

  for s = 1:nscale,       
      EO = submat(EOArray,s,cols);  % Extract even and odd filter 
      E = real(EO); O = imag(EO);
      Energy_ThisOrient = Energy_ThisOrient ...
        + E.*MeanE + O.*MeanO - abs(E.*MeanO - O.*MeanE);
  end

  medianE2n = median(reshape(abs(submat(EOArray,1,cols)).^2,1,rows*cols));
  meanE2n = -medianE2n/log(0.5);
  estMeanE2n = [estMeanE2n meanE2n];

  noisePower = meanE2n/EM_n;                       % Estimate of noise power.

  % Now estimate the total energy^2 due to noise
  % Estimate for sum(An^2) + sum(Ai.*Aj.*(cphi.*cphj + sphi.*sphj))

  EstSumAn2 = zero;
  for s = 1:nscale
    EstSumAn2 = EstSumAn2+submat(ifftFilterArray,s,cols).^2;
  end

  EstSumAiAj = zero;
  for si = 1:(nscale-1)
    for sj = (si+1):nscale
      EstSumAiAj = EstSumAiAj + submat(ifftFilterArray,si,cols).*submat(ifftFilterArray,sj,cols);
    end
  end

  EstNoiseEnergy2 = 2*noisePower*sum(sum(EstSumAn2)) + 4*noisePower*sum(sum(EstSumAiAj));

  tau = sqrt(EstNoiseEnergy2/2);                     % Rayleigh parameter
  EstNoiseEnergy = tau*sqrt(pi/2);                   % Expected value of noise energy
  EstNoiseEnergySigma = sqrt( (2-pi/2)*tau^2 );

  T =  EstNoiseEnergy + k*EstNoiseEnergySigma;       % Noise threshold

  T = T/1.7;        % Empirical rescaling of the estimated noise effect to 
                    % suit the PC_2 phase congruency measure

  Energy_ThisOrient = max(Energy_ThisOrient - T, zero);  % Apply noise threshold

  width = sumAn_ThisOrient ./ (maxAn + epsilon) / nscale;    

  % Now calculate the sigmoidal weighting function for this orientation.

  weight = 1.0 ./ (1 + exp( (cutOff - width)*g)); 

  % Apply weighting

  Energy_ThisOrient =   weight.*Energy_ThisOrient;

  % Update accumulator matrix for sumAn and totalEnergy

  totalSumAn  = totalSumAn + sumAn_ThisOrient;
  totalEnergy = totalEnergy + Energy_ThisOrient;

  if(o == 1),
    maxEnergy = Energy_ThisOrient;
    featType = E + i*O;
  else
    change = Energy_ThisOrient > maxEnergy;
    orientation = (o - 1).*change + orientation.*(~change);
    featType = (E+i*O).*change + featType.*(~change);
    maxEnergy = max(maxEnergy, Energy_ThisOrient);
  end

end  % For each orientation


phaseCongruency = totalEnergy ./ (totalSumAn + epsilon);

orientation = orientation * (180 / norient);

featType = featType*i;   % Rotate feature phase angles by 90deg so that 0
                         % phase corresponds to a step edge (this is a
                         % fudge I must have something the wrong way
                         % around somewhere)

phaseCongruency = exp(phaseCongruency);
phaseCongruency = phaseCongruency./(sum(sum(phaseCongruency)));

%-------------------------------------------------------------------------
%------------------------------Inner Function ----------------------------
%-------------------------------------------------------------------------
function a = submat(big,i,cols)

a = big(:,((i-1)*cols+1):(i*cols));

%-------------------------------------------------------------------------
%------------------------------Inner Function ----------------------------
%-------------------------------------------------------------------------

function B=inpaint_nans(A,method)
% inpaint_nans: in-paints over nans in an array
% usage: B=inpaint_nans(A)
%
% solves approximation to one of several pdes to
% interpolate and extrapolate holes
%
% arguments (input):
%   A - nxm array with some NaNs to be filled in
%
%   method - (OPTIONAL) scalar numeric flag - specifies
%       which approach (or physical metaphor to use
%       for the interpolation.) All methods are capable
%       of extrapolation, some are better than others.
%       There are also speed differences, as well as
%       accuracy differences for smooth surfaces.
%
%       methods {0,1,2} use a simple plate metaphor
%       methods {3} use a better plate equation,
%                   but will be slower
%       methods 4 use a spring metaphor
%
%       method == 0 --> (DEFAULT) see method 1, but
%         this method does not build as large of a
%         linear system in the case of only a few
%         NaNs in a large array.
%         Extrapolation behavior is linear.
%         
%       method == 1 --> simple approach, applies del^2
%         over the entire array, then drops those parts
%         of the array which do not have any contact with
%         NaNs. Uses a least squares approach, but it
%         does not touch existing points.
%         In the case of small arrays, this method is
%         quite fast as it does very little extra work.
%         Extrapolation behavior is linear.
%         
%       method == 2 --> uses del^2, but solving a direct
%         linear system of equations for nan elements.
%         This method will be the fastest possible for
%         large systems since it uses the sparsest
%         possible system of equations. Not a least
%         squares approach, so it may be least robust
%         to noise on the boundaries of any holes.
%         This method will also be least able to
%         interpolate accurately for smooth surfaces.
%         Extrapolation behavior is linear.
%         
%       method == 3 --+ See method 0, but uses del^4 for
%         the interpolating operator. This may result
%         in more accurate interpolations, at some cost
%         in speed.
%         
%       method == 4 --+ Uses a spring metaphor. Assumes
%         springs (with a nominal length of zero)
%         connect each node with every neighbor
%         (horizontally, vertically and diagonally)
%         Since each node tries to be like its neighbors,
%         extrapolation is as a constant function where
%         this is consistent with the neighboring nodes.
%
%
% arguments (output):
%   B - nxm array with NaNs replaced

% I always need to know which elements are NaN,
% and what size the array is for any method
[n,m]=size(A);
nm=n*m;
k=isnan(A(:));

% list the nodes which are known, and which will
% be interpolated
nan_list=find(k);
known_list=find(~k);

% how many nans overall
nan_count=length(nan_list);

% convert NaN indices to (r,c) form
% nan_list==find(k) are the unrolled (linear) indices
% (row,column) form
[nr,nc]=ind2sub([n,m],nan_list);

% both forms of index in one array:
% column 1 == unrolled index
% column 2 == row index
% column 3 == column index
nan_list=[nan_list,nr,nc];

% supply default method
if (nargin<2)|isempty(method)
  method = 0;
elseif ~ismember(method,0:4)
  error 'If supplied, method must be one of: {0,1,2,3,4}.'
end

% for different methods
switch method
 case 0
  % The same as method == 1, except only work on those
  % elements which are NaN, or at least touch a NaN.
  
  % horizontal and vertical neighbors only
  talks_to = [-1 0;0 -1;1 0;0 1];
  neighbors_list=identify_neighbors(n,m,nan_list,talks_to);
  
  % list of all nodes we have identified
  all_list=[nan_list;neighbors_list];
  
  % generate sparse array with second partials on row
  % variable for each element in either list, but only
  % for those nodes which have a row index > 1 or < n
  L = find((all_list(:,2) > 1) & (all_list(:,2) < n)); 
  nl=length(L);
  if nl>0
    fda=sparse(repmat(all_list(L,1),1,3), ...
      repmat(all_list(L,1),1,3)+repmat([-1 0 1],nl,1), ...
      repmat([1 -2 1],nl,1),nm,nm);
  else
    fda=spalloc(n*m,n*m,size(all_list,1)*5);
  end
  
  % 2nd partials on column index
  L = find((all_list(:,3) > 1) & (all_list(:,3) < m)); 
  nl=length(L);
  if nl>0
    fda=fda+sparse(repmat(all_list(L,1),1,3), ...
      repmat(all_list(L,1),1,3)+repmat([-n 0 n],nl,1), ...
      repmat([1 -2 1],nl,1),nm,nm);
  end
  
  % eliminate knowns
  rhs=-fda(:,known_list)*A(known_list);
  k=find(any(fda(:,nan_list(:,1)),2));
  
  % and solve...
  B=A;
  B(nan_list(:,1))=fda(k,nan_list(:,1))\rhs(k);
  
 case 1
  % least squares approach with del^2. Build system
  % for every array element as an unknown, and then
  % eliminate those which are knowns.

  % Build sparse matrix approximating del^2 for
  % every element in A.
  % Compute finite difference for second partials
  % on row variable first
  [i,j]=ndgrid(2:(n-1),1:m);
  ind=i(:)+(j(:)-1)*n;
  np=(n-2)*m;
  fda=sparse(repmat(ind,1,3),[ind-1,ind,ind+1], ...
      repmat([1 -2 1],np,1),n*m,n*m);
  
  % now second partials on column variable
  [i,j]=ndgrid(1:n,2:(m-1));
  ind=i(:)+(j(:)-1)*n;
  np=n*(m-2);
  fda=fda+sparse(repmat(ind,1,3),[ind-n,ind,ind+n], ...
      repmat([1 -2 1],np,1),nm,nm);
  
  % eliminate knowns
  rhs=-fda(:,known_list)*A(known_list);
  k=find(any(fda(:,nan_list),2));
  
  % and solve...
  B=A;
  B(nan_list(:,1))=fda(k,nan_list(:,1))\rhs(k);
  
 case 2
  % Direct solve for del^2 BVP across holes

  % generate sparse array with second partials on row
  % variable for each nan element, only for those nodes
  % which have a row index > 1 or < n
  L = find((nan_list(:,2) > 1) & (nan_list(:,2) < n)); 
  nl=length(L);
  if nl>0
    fda=sparse(repmat(nan_list(L,1),1,3), ...
      repmat(nan_list(L,1),1,3)+repmat([-1 0 1],nl,1), ...
      repmat([1 -2 1],nl,1),n*m,n*m);
  else
    fda=spalloc(n*m,n*m,size(nan_list,1)*5);
  end
  
  % 2nd partials on column index
  L = find((nan_list(:,3) > 1) & (nan_list(:,3) < m)); 
  nl=length(L);
  if nl>0
    fda=fda+sparse(repmat(nan_list(L,1),1,3), ...
      repmat(nan_list(L,1),1,3)+repmat([-n 0 n],nl,1), ...
      repmat([1 -2 1],nl,1),n*m,n*m);
  end
  
  % fix boundary conditions at extreme corners
  % of the array in case there were nans there
  if ismember(1,nan_list(:,1))
    fda(1,[1 2 n+1])=[-2 1 1];
  end
  if ismember(n,nan_list(:,1))
    fda(n,[n, n-1,n+n])=[-2 1 1];
  end
  if ismember(nm-n+1,nan_list(:,1))
    fda(nm-n+1,[nm-n+1,nm-n+2,nm-n])=[-2 1 1];
  end
  if ismember(nm,nan_list(:,1))
    fda(nm,[nm,nm-1,nm-n])=[-2 1 1];
  end
  
  % eliminate knowns
  rhs=-fda(:,known_list)*A(known_list);
  
  % and solve...
  B=A;
  k=nan_list(:,1);
  B(k)=fda(k,k)\rhs(k);
  
 case 3
  % The same as method == 0, except uses del^4 as the
  % interpolating operator.
  
  % del^4 template of neighbors
  talks_to = [-2 0;-1 -1;-1 0;-1 1;0 -2;0 -1; ...
      0 1;0 2;1 -1;1 0;1 1;2 0];
  neighbors_list=identify_neighbors(n,m,nan_list,talks_to);
  
  % list of all nodes we have identified
  all_list=[nan_list;neighbors_list];
  
  % generate sparse array with del^4, but only
  % for those nodes which have a row & column index
  % >= 3 or <= n-2
  L = find( (all_list(:,2) >= 3) & ...
            (all_list(:,2) <= (n-2)) & ...
            (all_list(:,3) >= 3) & ...
            (all_list(:,3) <= (m-2)));
  nl=length(L);
  if nl>0
    % do the entire template at once
    fda=sparse(repmat(all_list(L,1),1,13), ...
        repmat(all_list(L,1),1,13) + ...
        repmat([-2*n,-n-1,-n,-n+1,-2,-1,0,1,2,n-1,n,n+1,2*n],nl,1), ...
        repmat([1 2 -8 2 1 -8 20 -8 1 2 -8 2 1],nl,1),nm,nm);
  else
    fda=spalloc(n*m,n*m,size(all_list,1)*5);
  end
  
  % on the boundaries, reduce the order around the edges
  L = find((((all_list(:,2) == 2) | ...
             (all_list(:,2) == (n-1))) & ...
            (all_list(:,3) >= 2) & ...
            (all_list(:,3) <= (m-1))) | ...
           (((all_list(:,3) == 2) | ...
             (all_list(:,3) == (m-1))) & ...
            (all_list(:,2) >= 2) & ...
            (all_list(:,2) <= (n-1))));
  nl=length(L);
  if nl>0
    fda=fda+sparse(repmat(all_list(L,1),1,5), ...
      repmat(all_list(L,1),1,5) + ...
        repmat([-n,-1,0,+1,n],nl,1), ...
      repmat([1 1 -4 1 1],nl,1),nm,nm);
  end
  
  L = find( ((all_list(:,2) == 1) | ...
             (all_list(:,2) == n)) & ...
            (all_list(:,3) >= 2) & ...
            (all_list(:,3) <= (m-1)));
  nl=length(L);
  if nl>0
    fda=fda+sparse(repmat(all_list(L,1),1,3), ...
      repmat(all_list(L,1),1,3) + ...
        repmat([-n,0,n],nl,1), ...
      repmat([1 -2 1],nl,1),nm,nm);
  end
  
  L = find( ((all_list(:,3) == 1) | ...
             (all_list(:,3) == m)) & ...
            (all_list(:,2) >= 2) & ...
            (all_list(:,2) <= (n-1)));
  nl=length(L);
  if nl>0
    fda=fda+sparse(repmat(all_list(L,1),1,3), ...
      repmat(all_list(L,1),1,3) + ...
        repmat([-1,0,1],nl,1), ...
      repmat([1 -2 1],nl,1),nm,nm);
  end
  
  % eliminate knowns
  rhs=-fda(:,known_list)*A(known_list);
  k=find(any(fda(:,nan_list(:,1)),2));
  
  % and solve...
  B=A;
  B(nan_list(:,1))=fda(k,nan_list(:,1))\rhs(k);
  
 case 4
  % Spring analogy
  % interpolating operator.
  
  % list of all springs between a node and a horizontal
  % or vertical neighbor
  hv_list=[-1 -1 0;1 1 0;-n 0 -1;n 0 1];
  hv_springs=[];
  for i=1:4
    hvs=nan_list+repmat(hv_list(i,:),nan_count,1);
    k=(hvs(:,2)>=1) & (hvs(:,2)<=n) & (hvs(:,3)>=1) & (hvs(:,3)<=m);
    hv_springs=[hv_springs;[nan_list(k,1),hvs(k,1)]];
  end

  % delete replicate springs
  hv_springs=unique(sort(hv_springs,2),'rows');
  
  % build sparse matrix of connections, springs
  % connecting diagonal neighbors are weaker than
  % the horizontal and vertical springs
  nhv=size(hv_springs,1);
  springs=sparse(repmat((1:nhv)',1,2),hv_springs, ...
     repmat([1 -1],nhv,1),nhv,nm);
  
  % eliminate knowns
  rhs=-springs(:,known_list)*A(known_list);
  
  % and solve...
  B=A;
  B(nan_list(:,1))=springs(:,nan_list(:,1))\rhs;
  
end

% ====================================================
%      end of main function
% ====================================================
% ====================================================
%      begin subfunctions
% ====================================================
function neighbors_list=identify_neighbors(n,m,nan_list,talks_to)
% identify_neighbors: identifies all the neighbors of
%   those nodes in nan_list, not including the nans
%   themselves
%
% arguments (input):
%  n,m - scalar - [n,m]=size(A), where A is the
%      array to be interpolated
%  nan_list - array - list of every nan element in A
%      nan_list(i,1) == linear index of i'th nan element
%      nan_list(i,2) == row index of i'th nan element
%      nan_list(i,3) == column index of i'th nan element
%  talks_to - px2 array - defines which nodes communicate
%      with each other, i.e., which nodes are neighbors.
%
%      talks_to(i,1) - defines the offset in the row
%                      dimension of a neighbor
%      talks_to(i,2) - defines the offset in the column
%                      dimension of a neighbor
%      
%      For example, talks_to = [-1 0;0 -1;1 0;0 1]
%      means that each node talks only to its immediate
%      neighbors horizontally and vertically.
% 
% arguments(output):
%  neighbors_list - array - list of all neighbors of
%      all the nodes in nan_list

if ~isempty(nan_list)
  % use the definition of a neighbor in talks_to
  nan_count=size(nan_list,1);
  talk_count=size(talks_to,1);
  
  nn=zeros(nan_count*talk_count,2);
  j=[1,nan_count];
  for i=1:talk_count
    nn(j(1):j(2),:)=nan_list(:,2:3) + ...
        repmat(talks_to(i,:),nan_count,1);
    j=j+nan_count;
  end
  
  % drop those nodes which fall outside the bounds of the
  % original array
  L = (nn(:,1)<1)|(nn(:,1)>n)|(nn(:,2)<1)|(nn(:,2)>m); 
  nn(L,:)=[];
  
  % form the same format 3 column array as nan_list
  neighbors_list=[sub2ind([n,m],nn(:,1),nn(:,2)),nn];
  
  % delete replicates in the neighbors list
  neighbors_list=unique(neighbors_list,'rows');
  
  % and delete those which are also in the list of NaNs.
  neighbors_list=setdiff(neighbors_list,nan_list,'rows');
  
else
  neighbors_list=[];
end


Contact us