Reworked warping entirely

warp.m

-function [Value,Mask,Pos,DispX,DispY] = warp(Value,Mask,Pos,DispX,DispY,varargin)
+function [Value,Mask,Pos,dX,dY] = warp(Value,Mask,Pos,dX,dY,varargin)
 %WARP Summary of this function goes here
 %   Detailed explanation goes here
 
 %% Set parameter values according to given arguments
 interps = {'nearest','linear','cubic','spline'};
-boundaryTypes = {'current','any','all'};
 directions = {'forward','backward'};
 
 p = inputParser; p.KeepUnmatched = true; p.StructExpand = true;
-p.addParameter('interpMethod', 'nearest' , @(x) any(validatestring(x,interps)));
-p.addParameter('SRBoundary'  , 'current' , @(x) any(validatestring(x,boundaryTypes)));
-p.addParameter('LFBoundary'  , false     , @islogical);
-p.addParameter('direction'   , 'backward', @(x) any(validatestring(x,directions)));
-p.addParameter('relative'    , true      , @islogical);
+p.addParameter('interpMethod', 'cubic'  , @(x) any(validatestring(x,interps   )));
+p.addParameter('direction'   , 'forward', @(x) any(validatestring(x,directions)));
+p.addParameter('backproject' , true     , @islogical);
 
 p.parse(varargin{:});
-interpMethod = p.Results.interpMethod; % interpolation method of Gridded Interpolant
-SRBoundary = p.Results.SRBoundary; % Query values to retain
-LFBoundary = p.Results.LFBoundary; % Retain values inside original boundaries
-direction = p.Results.direction; % Forward or backward projection
-relative = p.Results.relative; % Scale disparity relative to reference view
+interpMethod = p.Results.interpMethod;
+direction = p.Results.direction;
+backproject = p.Results.backproject;
 
-%%
+%% Compute reference grid
+% Reference size
 Size = size(Value);
-ValMask = ~isnan(Value);
 
-% Lightfield grid vectors representing super-ray bounding box
-ugv = Pos(1)+reshape(1:Size(1),[],1);
-vgv = Pos(2)+reshape(1:Size(2),1,[]);
-xgv = Pos(3)+reshape(1:Size(3),1,1,[]);
-ygv = Pos(4)+reshape(1:Size(4),1,1,1,[]);
+% Mask bounding box
+gv = subarrayIdx(Mask);
 
-gv = {ugv,vgv,xgv,ygv};
+% Mask grid coordinates
+[~,~,X,Y] = ndgrid(gv{:});
 
-% Lightfield coordinates (4D)
-[U,V,X,Y] = ndgrid(ugv,vgv,xgv,ygv);
+% Mask grid vectors
+[ugv,vgv,~,~] = deal(gv{:});
+
+ugv = reshape(ugv,[],1);
+vgv = reshape(vgv,1,[]);
 
 % Reference view
-uref = median(ugv);
-vref = median(vgv);
+u_ref = floor(numel(ugv)/2)+1;
+v_ref = floor(numel(vgv)/2)+1;
+
+% Disparity
+DispX = dX*(ugv-u_ref);
+DispY = dY*(vgv-v_ref);
+
+if backproject
+    % Temporary mask
+    tempMask = Mask(gv{:});
+    
+    % Backproject coordinates
+    SRX = X + DispX;
+    SRY = Y + DispY;
+    
+    % Backprojected boundaries
+    minSRX = floor(min(SRX(tempMask))); maxSRX = ceil(max(SRX(tempMask)));
+    minSRY = floor(min(SRY(tempMask))); maxSRY = ceil(max(SRY(tempMask)));
+    
+    minSRX = max(minSRX,1); maxSRX = min(maxSRX,Size(3));
+    minSRY = max(minSRY,1); maxSRY = min(maxSRY,Size(4));
+    
+    % Reference grid vectors
+    xgv = minSRX:maxSRX;
+    ygv = minSRY:maxSRY;
+    
+    % Reference bounding box
+    gv = {ugv,vgv,xgv,ygv};
+end
 
-% Fill NaN values to allow cubic interpolation (way faster than fillmissing)
-seAng = strel('arbitrary',[0,1,0]|[0,1,0]');
-seSpa = strel('arbitrary',reshape([0,1,0]|[0,1,0]',1,1,3,3));
+% Reference grid
+[~,~,X,Y] = ndgrid(gv{:});
 
-if strcmp(interpMethod,'cubic')&&any(isnan(Value(:)))
-    Value = extend(Value,seSpa,6);
-    Value = extend(Value,seAng,6);
-end
-if strcmp(interpMethod,'cubic')&&any(isnan(DispX(:)))
-    DispX = extend(DispX,seSpa,6);
-    DispX = extend(DispX,seAng,6);
-end
-if strcmp(interpMethod,'cubic')&&any(isnan(DispY(:)))
-    DispY = extend(DispY,seSpa,6);
-    DispY = extend(DispY,seAng,6);
-end
+%% Crop reference
+Value = Value(gv{:});
+Mask  = Mask (gv{:});
+Pos_ = cellfun(@(x) x(1),gv)-1;
+End_ = cellfun(@(x) x(end),gv);
 
-% If view-wise disparity, scale accordingly
-if relative
-    DispX_ = (ugv-uref).*DispX;
-    DispY_ = (vgv-vref).*DispY;
-else
-    DispX_ = DispX;
-    DispY_ = DispY;
-end
+%% Compute query grid
+% Project coordinates
+SRX = X - DispX;
+SRY = Y - DispY;
 
-% Check size of interpolated data (for GriddedInterpolant)
-if numel(Value)==1, szVal=1; else, szVal = size(Value); end
-singValue = szVal==1;
-if numel(DispX)==1, szDispX=1; else, szDispX = size(DispX); end
-singDispX = szDispX==1;
-if numel(DispY)==1, szDispY=1; else, szDispY = size(DispX); end
-singDispY = szDispY==1;
-
-% expand any singleton dimension (need at least 2)
-if any(singValue),Value = repmat(Value,2*singValue); end
-szVal = size(Value);
-
-singDispX(end+1:numel(szVal)) = 1;
-szDispX(end+1:numel(szVal)) = 1;
-singDispY(end+1:numel(szVal)) = 1;
-szDispY(end+1:numel(szVal)) = 1;
-if any(singDispX),DispX = repmat(DispX,szVal./szDispX); end
-if any(singDispY),DispY = repmat(DispY,szVal./szDispY); end
+% Projected boundaries
+minPosX = floor(min(SRX(Mask))); maxPosX = ceil(max(SRX(Mask)));
+minPosY = floor(min(SRY(Mask))); maxPosY = ceil(max(SRY(Mask)));
 
-% Initialize interpolants
-ValueInt = griddedInterpolant(gv,Value,interpMethod,'nearest'); % Interpolant on pixel values
-DispXInt = griddedInterpolant(gv,DispX,interpMethod,'nearest'); % Interpolant on pixel values
-DispYInt = griddedInterpolant(gv,DispY,interpMethod,'nearest'); % Interpolant on pixel values
-MaskInt = griddedInterpolant(gv,double(Mask),'nearest','nearest'); % Interpolant on masked pixels
+% Query grid vectors
+uqgv = ugv;
+vqgv = vgv;
+xqgv = minPosX:maxPosX;
+yqgv = minPosY:maxPosY;
 
-disp(['Value:  "' ValueInt.Method '" interp., "' ValueInt.ExtrapolationMethod '" extrap.']);
-disp(['Mask:   "'  MaskInt.Method '" interp., "'  MaskInt.ExtrapolationMethod '" extrap.']);
+% Query bounding box
+qgv = {uqgv,vqgv,xqgv,yqgv};
 
-% Projected values
-projX = X-DispX_;
-projY = Y-DispY_;
+% Query grid
+[Uq,Vq,Xq,Yq] = ndgrid(qgv{:});
 
-% Rounded projected values
+Xq = Xq + DispX;
+Yq = Yq + DispY;
+
+% Use invertible grid coordinates to compute mask
 switch direction
     case 'forward'
-        roundProjX = floor(projX+0.5);
-        roundProjY = floor(projY+0.5);
+        Xq_Mask = floor(Xq+0.5); 
+        Yq_Mask = floor(Yq+0.5); 
     case 'backward'
-        roundProjX = ceil(projX-0.5);
-        roundProjY = ceil(projY-0.5);
+        Xq_Mask = ceil(Xq-0.5); 
+        Yq_Mask = ceil(Yq-0.5); 
 end
 
-% Query coordinates
-Uq = U;
-Vq = V;
-Xq = X+roundProjX-projX;
-Yq = Y+roundProjY-projY;
+%% Interpolate
+% Extend signal to allow cubic interpolation (way faster than fillmissing)
+seAng = strel('arbitrary',[0,1,0]|[0,1,0]');
+seSpa = strel('arbitrary',reshape([0,1,0]|[0,1,0]',1,1,3,3));
 
-% Interpolate at query coordinates
-Value = ValueInt(Uq,Vq,Xq,Yq);
-DispX = DispXInt(Uq,Vq,Xq,Yq);
-DispY = DispYInt(Uq,Vq,Xq,Yq);
-Mask  = MaskInt(Uq,Vq,Xq,Yq);
-
-Mask(isnan(Mask))=0;
-Mask = logical(Mask);
-
-% Create mask corresponding to query points
-switch SRBoundary
-    case 'current'
-        % inside boundary in current view
-        % Mask = Mask;
-    case 'any'
-        % inside boundary in at least one view
-        Mask = repmat(any(any(Mask,1),2),size(Mask,1),size(Mask,2),1,1);
-    case 'all'
-        % keep all query points
-        Mask = true(Size);
+if strcmp(interpMethod,'cubic')&&any(isnan(Value(:)))
+    Value = extend(Value,seSpa,6);
+    Value = extend(Value,seAng,6);
 end
 
-if LFBoundary
-    % Mask corresponding to query points inside lightfield boundary
-    Mask = Mask & ...
-        Xq-Pos(3)>=0.5 & Xq-Pos(3)<=Size(3)+0.5 & ... % <==> round(Xq)>=1 & round(Xq)<=LFSize(3)
-        Yq-Pos(4)>=0.5 & Yq-Pos(4)<=Size(4)+0.5;      %    & round(Yq)>=1 & round(Xq)<=LFSize(4)
+% Expand any singleton dimension for interpolation (at least 3 for cubic)
+szVal = size(Value); szVal(end+1:4)=1;
+
+singDims = double(szVal==1);
+if any(singDims)
+    Value = repmat(Value,2*double(singDims)+1);
+    for it = find(singDims)
+        gv{it} = gv{it}-1:gv{it}+1;
+    end
 end
 
-% Min and Max x and y coordinates of the pixels of the super-ray
-minProjX = floor(min(projX(:))); maxProjX = ceil(max(projX(:)));
-minProjY = floor(min(projY(:))); maxProjY = ceil(max(projY(:)));
+dualDims = double(szVal==2);
+if any(dualDims)
+    Value = repmat(Value,double(dualDims)+1);
+    for it = find(dualDims)
+        gv{it} = [gv{it},gv{it}(end)+1];
+    end
+    
+end
 
-% Query grid vectors
-projxgv = reshape(minProjX:maxProjX,1,1,[]);
-projygv = reshape(minProjY:maxProjY,1,1,1,[]);
 
-% New position
-projPos = Pos;
-projPos(3) = minProjX-1;
-projPos(4) = minProjY-1;
-
-% New size
-projSize = Size;    
-projSize(3) = numel(projxgv);
-projSize(4) = numel(projygv);
-
-% Determine indices corresponding to original and warped pixels
-ind = sub2ind_(Size,Pos,U(ValMask),V(ValMask),X(ValMask),Y(ValMask));
-roundInd = sub2ind_(projSize,projPos,U(ValMask),V(ValMask),roundProjX(ValMask),roundProjY(ValMask));
-roundMaskInd = sub2ind_(projSize,projPos,U(Mask),V(Mask),roundProjX(Mask),roundProjY(Mask));
-
-% Performing warping 
-projValue = nan(projSize);
-projMask  = false(projSize);
-projDispX = nan(projSize);
-projDispY = nan(projSize);
-
-projValue(roundInd) = Value(ind);
-projMask (roundMaskInd) = true;
-projDispX(roundInd) = DispX(ind);
-projDispY(roundInd) = DispY(ind);
-
-Value = projValue;
-Mask  = projMask;
-DispX = -projDispX;
-DispY = -projDispY;
-Pos = projPos;
-end
+% Cast mask as double for interpolation
+Mask = double(Mask);
+
+% Initialize interpolants
+ValueInt = griddedInterpolant(gv,Value,interpMethod,'nearest');
+MaskInt  = griddedInterpolant(gv,Mask ,'nearest'   ,'nearest');
 
-function ind = sub2ind_(Size,Pos,varargin)
-subs = varargin;
-Pos = num2cell(Pos);
+% Interpolation
+Value = ValueInt(Uq,Vq,Xq,Yq);
+Mask  = MaskInt(Uq,Vq,Xq_Mask,Yq_Mask); Mask = logical(Mask);
 
-subs = cellfun(@minus,subs,Pos,'UniformOutput',false);
-subs = cellfun(@(x) x(:),subs,'UniformOutput',false);
-ind = sub2ind(Size,subs{:});
-end
+% Mask of query coordinates inside reference boundary
+refMask = Xq_Mask>Pos_(3) & Xq_Mask<End_(3)+1 & Yq_Mask>Pos_(4) & Yq_Mask<End_(4)+1;
 
-function varargout = ind2sub_(Size,Pos,ind)
-subs = cell(1,numel(Size));
-[subs{:}] = ind2sub(Size,ind);
-Pos = num2cell(Pos);
+% New mask
+Mask = Mask & refMask;
+
+% New position
+Pos = Pos + cellfun(@(x) x(1),qgv)-1;
 
-varargout = cellfun(@plus,subs,Pos,'UniformOutput',false);
+% New disparities
+dX = -dX;
+dY = -dY;
+end
+
+%% Auxiliary functions
+% Compute bounding box of mask
+function gv = subarrayIdx(Mask)
+s = regionprops(Mask,'SubarrayIdx');
+if numel(s)==1
+    gv = s.SubarrayIdx;
+else
+    gv = {s.SubarrayIdx};
+    [ugv,vgv,xgv,ygv] = cellfun(@(x)deal(x{:}),gv,'UniformOutput',false);
+    
+    for it = 1:numel(s)-1
+        ugv = [union(ugv{1},ugv{2}),ugv(3:end)];
+        vgv = [union(vgv{1},vgv{2}),vgv(3:end)];
+        xgv = [union(xgv{1},xgv{2}),xgv(3:end)];
+        ygv = [union(ygv{1},ygv{2}),ygv(3:end)];
+    end
+    
+    ugv = min(ugv{:}):max(ugv{:});
+    vgv = min(vgv{:}):max(vgv{:});
+    xgv = min(xgv{:}):max(xgv{:});
+    ygv = min(ygv{:}):max(ygv{:});
+    gv = {ugv,vgv,xgv,ygv};
+end
 end
 
+% Extend signal
 function Value = extend(Value,se,numIter)
 for i=1:numIter
-        Mask = ~isnan(Value);
-        Temp = Value;
-        Temp(~Mask)=-Inf;
-        Temp = imdilate(Temp,se);
-        Temp(isinf(Temp))=NaN;
-        Value(~Mask) = Temp(~Mask);
+    Mask = isnan(Value);
+    Temp = Value;
+    Temp(Mask)=-Inf;
+    Temp = imdilate(Temp,se);
+    Temp(isinf(Temp))=NaN;
+    Value(Mask) = Temp(Mask);
 end
 end
\ No newline at end of file

warpSet.m

@@ -4,19 +4,17 @@ function SRSet = warpSet(SRSet,varargin)
 
 [Val,Lab,Pos,DispX,DispY] = SR.SetToFields(SRSet);
 
-Mask = Lab;
-
 for it = 1:numel(Lab)
-    Mask{it} = Lab{it}==it;
+    Lab{it} = Lab{it}==it;
 end
 
 for it = 1:numel(Val)
-    [Val{it},Mask{it},Pos{it},DispX{it},DispY{it}] = SR.warp(...
-        Val{it},Mask{it},Pos{it},DispX{it},DispY{it},varargin{:});
+    [Val{it},Lab{it},Pos{it},DispX{it},DispY{it}] = SR.warp(...
+        Val{it},Lab{it},Pos{it},DispX{it},DispY{it},varargin{:});
 end
 
 for it = 1:numel(Lab)
-    Lab{it} = it*double(Mask{it});
+    Lab{it} = it*double(Lab{it});
 end
 
 SRSet = SR.FieldsToSet(Val,Lab,Pos,DispX,DispY);