diff --git a/warp.m b/warp.m
index 8d0b3edd29af56d17aab949514755ff98765364c..145cd8f1aca8492649e583f066c0e1c28559409f 100644
--- a/warp.m
+++ b/warp.m
@@ -1,213 +1,200 @@
-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
diff --git a/warpSet.m b/warpSet.m
index 4f57cbe28e257809455e2263017fb3a67786ba61..d67a1221c73ce87d11520388a7c5de3cc3ca4f12 100644
--- a/warpSet.m
+++ b/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);