Changed super-ray computation to allow scattered data

LFtoSR.m

-function [SRSet,varargout] = LFtoSR(LFSet,Xq,Yq,Method,varargin)
+function [SRSet,SRDisp] = LFtoSR(LFSet,LFDisp,Method)
 %LFTOSR Builds a super-ray set from a light-field
 %   LFTOSR(LFSet)
 
 %% Initialize constants and interpolants
 % Lightfield properties
+ImgSize = LFSet.ImgSize(1:2);
+numChan = LFSet.ImgSize(3);
 ImgRes = LFSet.ImgRes;
 Offset = LFSet.Offset;
 Color  = LFSet.Color;
 Label  = LFSet.Label;
-ResOff = num2cell(floor(ImgRes/2)+1);
-numChan = size(Color,3);
-
-% Number of super-rays to be in the set
-numLab = numel(Xq);
+ResOff = floor(ImgRes/2)+1;
+LFSize = [ImgSize,ImgRes];
 numView = prod(ImgRes);
+numLab = max(Label(:));
 
-% Initialize super-ray set
-[SRImgSize,SRImgRes,SROff,SRCol,SRLab] = deal(cell(numLab,1));
-varargout = cell(1,nargin-4);
-varargout(:) = deal({cell(numLab,1)});
-
-% Initialize interpolants
-ValInt = griddedInterpolant(Color(:,:,1),Method,'none');
-
-sizeMsg = 0;
+% Initialize super-ray fields
+[SRImgSize,SRImgRes,SROff,SRCol,SRLab,SRDisp] = deal(cell(numLab,1));
+[srxqgv,sryqgv] = deal(cell(numLab,1));
 
+% Replace missing values by zero (to avoid interpolation errors)
 Color(isnan(Color))=0;
 Label(isnan(Label))=0;
+LFDisp(isnan(LFDisp))=0;
 
-for it = 1:nargin-4
-    varargin{it}(isnan(varargin{it}))=0;
-end
+% Compute projected coordinates of reference samples
+gv = arrayfun(@(x) 1:x,LFSize,'UniformOutput',false);
+ugv = reshape(gv{3},1,1,[],1);
+vgv = reshape(gv{4},1,1,1,[]);
+ugv = ugv-ResOff(1);
+vgv = vgv-ResOff(2);
 
-%% Build super-ray set
-fprintf('\n');
+[LFX,LFY] = ndgrid(gv{:});
+LFX = LFX-ugv.*reshape(LFDisp,LFSize);
+LFY = LFY-vgv.*reshape(LFDisp,LFSize);
+
+progress('');
+
+%% Initialize super-rays
+fprintf('Super-ray initialization...\n');
 for lab = 1:numLab
-    fprintf(repmat('\b',1,sizeMsg));
-    msg = ['Constructing super-ray ' num2str(lab) '/' num2str(numLab)];
-    sizeMsg = numel(msg);
-    fprintf(msg);
+    msg = ['Initializing super-ray ',num2str(lab),'/',num2str(numLab),'\n'];
+    progress(msg);
     
-    % Set current super-ray size
-    SRSize = size(Xq{lab});
-    SRSize = [SRSize(1:2),numChan,SRSize(3:end)];
-    SRImgSize{lab} = SRSize(1:3);
-    SRImgRes{lab}  = SRSize(4:end);
-    
-    % Set current super-ray position (offset from reference frame center)
-    Xqmin = floor(min(Xq{lab}(:,:,ResOff{:})));
-    Yqmin = floor(min(Yq{lab}(:,:,ResOff{:})));
+    % Compute super-ray grid boundary
+    SRX = LFX(Label==lab);
+    SRY = LFY(Label==lab);
     
-    SROff{lab} = Offset+[min(Xqmin(:)),min(Yqmin(:)),1,1,1]-1;
-    
-    % Set current super-ray color, label and disparities
-    SRCol{lab} = nan(SRSize);
-    SRLab{lab} = nan(SRSize);
-    
-    for it = 1:nargin-4
-        varargout{it}{lab} = nan(SRSize);
-    end
+    srxqmin = floor(min(SRX(:)));
+    srxqmax = ceil (max(SRX(:)));
+    sryqmin = floor(min(SRY(:)));
+    sryqmax = ceil (max(SRY(:)));
     
-    ValInt.Method = Method;
-    ValInt.ExtrapolationMethod = 'none';
+    srxqgv{lab} = srxqmin:srxqmax;
+    sryqgv{lab} = sryqmin:sryqmax;
     
-    for c = 1:numChan
-        for v = 1:numView
-            ValInt.Values = Color(:,:,c,v);
-            SRCol{lab}(:,:,c,v) = ValInt(Xq{lab}(:,:,v),Yq{lab}(:,:,v));
-            %SRCol{lab}(:,:,c,v) = interp2(Color(:,:,c,v),Yq(:,:,v),Xq(:,:,v),Method,nan);
-        end
-    end
+    % Compute super-ray size and offset wrt reference origin
+    SRSize = [numel(srxqgv{lab}),numel(sryqgv{lab}),numChan,ImgRes];
+    SRImgSize{lab} = SRSize(1:3);
+    SRImgRes{lab}  = SRSize(4:end);
     
-    for it = 1:nargin-4
-        for c = 1:size(varargin{it},3)
-            for v = 1:numView
-                ValInt.Values = varargin{it}(:,:,c,v);
-                varargout{it}{lab}(:,:,c,v) = ValInt(Xq{lab}(:,:,v),Yq{lab}(:,:,v));
-                %varargout{it}{lab}(:,:,c,v) = interp2(Color(:,:,c,v),Yq(:,:,v),Xq(:,:,v),Method,nan);
-            end
-        end
-    end
+    SROff{lab} = Offset+[srxqmin,sryqmin,1,1,1]-1;
     
-    ValInt.Method = 'nearest';
-    ValInt.ExtrapolationMethod = 'none';
+    % Initialize super-ray color, label and disparities
+    SRCol {lab} = nan([numel(srxqgv{lab}),numel(sryqgv{lab}),numChan,ImgRes]);
+    SRLab {lab} = nan([numel(srxqgv{lab}),numel(sryqgv{lab}),1      ,ImgRes]);
+    SRDisp{lab} = nan([numel(srxqgv{lab}),numel(sryqgv{lab}),1      ,ImgRes]);
+end
+fprintf('\n');
+
+progress('',0);
+
+%% Compute super-rays using interpolation
+fprintf('Super-ray interpolation...\n');
+for lab = 1:numLab
+    % Compute query grid for current super-ray
+    [SRXq,SRYq] = ndgrid(srxqgv{lab},sryqgv{lab});
     
+    % Compute super-ray view by view
     for v = 1:numView
-        ValInt.Values = Label(:,:,1,v);
-        SRLab{lab}(:,:,1,v) = ValInt(Xq{lab}(:,:,v),Yq{lab}(:,:,v));
-        %SRLab{lab}(:,:,1,v) = interp2(Label(:,:,c,v),Yq(:,:,v),Xq(:,:,v),Method,nan);
+        msg = ['Interpolating super-ray ',num2str(lab),'/',num2str(numLab),...
+            '\nView ',num2str(v),'/',num2str(numView),'\n'];
+        progress(msg);
+        
+        % Reduce number of reference points for speed
+        Mask = ...true(size(LFX(:,:,v)));
+            LFX(:,:,v)>=min(srxqgv{lab})-1 & LFX(:,:,v)<=max(srxqgv{lab})+1 & ...
+            LFY(:,:,v)>=min(sryqgv{lab})-1 & LFY(:,:,v)<=max(sryqgv{lab})+1;
+        
+        LFXsub = reshape(LFX(:,:,v),[],1);
+        LFYsub = reshape(LFY(:,:,v),[],1);
+        
+        LFXsub = LFXsub(Mask,:);
+        LFYsub = LFYsub(Mask,:);
+        
+        % Interpolate color
+        for c = 1:numChan
+            Colsub = reshape(Color(:,:,c,v),[],1);
+            Colsub = Colsub(Mask,:);
+            SRCol{lab}(:,:,c,v) = griddata(LFXsub,LFYsub,Colsub,SRXq,SRYq,Method);
+        end
+        
+        % Interpolate disparity
+        LFDispsub = reshape(LFDisp(:,:,1,v),[],1);
+        LFDispsub = LFDispsub(Mask,:);
+        SRDisp{lab}(:,:,1,v) = griddata(LFXsub,LFYsub,LFDispsub,SRXq,SRYq,Method);
+        
+        % Interpolate label
+        Labelsub = reshape(Label(:,:,1,v),[],1);
+        Labelsub = Labelsub(Mask,:);
+        Lab = griddata(LFXsub,LFYsub,Labelsub,SRXq,SRYq,'nearest');
+        
+        % Replace out of boundary values in label using cubic interpolation
+        NaNLab = isnan(griddata(LFXsub,LFYsub,Labelsub,SRXq,SRYq,'cubic'));
+        Lab(NaNLab) = nan;
+        SRLab{lab}(:,:,1,v) = Lab;
     end
 end
 fprintf('\n\n');
 
+% Create super-ray structure from fields
 SRSet = SR.FieldsToSet(SRImgSize,SRImgRes,SROff,SRCol,SRLab);
 
+    function progress(msg,varargin)
+        persistent sz
+        if isempty(sz); sz = 0; end
+        if nargin>1; sz = varargin{1}; end
+        
+        fprintf(repmat('\b',1,sz));
+        fprintf(msg);
+        sz = numel(msg)-count(msg,'\n');
+    end
 end
\ No newline at end of file

SRtoLF.m

-function [LFSet,LFDisp] = SRtoLF(SRSet,SRXq,SRYq,SRDisp,LFXq,LFYq,Method)
+function [LFSet,LFDisp] = SRtoLF(SRSet,SRDisp,lfxgv,lfygv,Method)
 %LFTOSR Builds a super-ray set from a light-field
-%   LFTOSR(LFSet)
+%   SRTOLF(LFSet)
 
 %%
-SRCol  = {SRSet.Color};
-SRLab  = {SRSet.Label};
+SRCol = {SRSet.Color};
+SRLab = {SRSet.Label};
+SROffset = {SRSet.Offset};
+SRImgSize = {SRSet.ImgSize};
+SRImgRes  = {SRSet.ImgRes};
 
 %% Initialize constants and interpolants
-% Lightfield properties
 numChan = size(SRCol{1},3);
+ImgSize = [numel(lfxgv),numel(lfygv),numChan];
+ImgRes = SRImgRes{1};
+Offset = [min(lfxgv),min(lfygv),1,1,1]-1;
 numLab = numel(SRCol);
-se = strel('square',3);
-
-LFSize = size(LFXq);
-LFSize = [LFSize(1:2),numChan,LFSize(3:end)];
-ImgSize = LFSize(1:3);
-ImgRes = LFSize(4:end);
 numView = prod(ImgRes);
-ResOff = num2cell(floor(ImgRes/2)+1);
 
-lfxqmin = floor(min(LFXq(:,:,ResOff{:}))); lfxqmin = min(lfxqmin(:));
-lfxqmax = ceil (max(LFXq(:,:,ResOff{:}))); lfxqmax = max(lfxqmax(:));
-lfyqmin = floor(min(LFYq(:,:,ResOff{:}))); lfyqmin = min(lfyqmin(:));
-lfyqmax = ceil (max(LFYq(:,:,ResOff{:}))); lfyqmax = max(lfyqmax(:));
+[LFXq,LFYq] = ndgrid(lfxgv,lfygv);
+[SRX,SRY] = deal(cell(numLab,1));
 
-lfxgv = lfxqmin:lfxqmax;
-lfygv = lfyqmin:lfyqmax;
+Color = nan([ImgSize(1:2),numChan,ImgRes]);
+Label = zeros([ImgSize(1:2),1,ImgRes]);
+LFDisp = -inf([ImgSize(1:2),1,ImgRes]);
 
-Offset = [lfxqmin,lfyqmin,1,1,1]-1;
+progress('',0);
 
-Color = nan(LFSize);
-Label = zeros(LFSize);
-LFDisp = -inf(LFSize);
+%% Compute super-ray coordinates
+for lab = 1:numLab
+    msg = ['Initializing super-ray ',num2str(lab),'/',num2str(numLab),'\n'];
+    progress(msg);
+    
+    % Compute grid coordinates
+    xgv = SROffset{lab}(1)+(1:SRImgSize{lab}(1));
+    ygv = SROffset{lab}(2)+(1:SRImgSize{lab}(2));
+    ugv = 1:SRImgRes{lab}(1);
+    ugv = ugv-floor(ugv(end)/2)-1;
+    vgv = 1:SRImgRes{lab}(2);
+    vgv = vgv-floor(vgv(end)/2)-1;
+    
+    % Reshape grid coordinate along corresponding dimension
+    xgv = reshape(xgv,[],1,1,1);
+    ygv = reshape(ygv,1,[],1,1);
+    ugv = reshape(ugv,1,1,[],1);
+    vgv = reshape(vgv,1,1,1,[]);
+    
+    % Compute projected coordinates of reference samples
+    Disp = reshape(SRDisp{lab},[SRImgSize{lab}(1:2),SRImgRes{lab}]);
+    SRX{lab} = xgv+ugv.*Disp.*ones(size(ygv)).*ones(size(vgv));
+    SRY{lab} = ygv+vgv.*Disp.*ones(size(xgv)).*ones(size(ugv));
+end
+fprintf('\n');
 
-sizeMsg = 0;
+progress('',0);
 
-fprintf('\n');
+%% Compute light field using interpolation
+% Compute light field super-ray by super-ray
 for lab = 1:numLab
-    srxqmin = floor(min(SRXq{lab}(:)));
-    srxqmax = ceil (max(SRXq{lab}(:)));
-    sryqmin = floor(min(SRYq{lab}(:)));
-    sryqmax = ceil (max(SRYq{lab}(:)));
+    % Compute super-ray boundary
+    srxqmin = floor(min(SRX{lab}(:)));
+    srxqmax = ceil (max(SRX{lab}(:)));
+    sryqmin = floor(min(SRY{lab}(:)));
+    sryqmax = ceil (max(SRY{lab}(:)));
     
+    % Reduce number of query points for speed
     xgv = lfxgv>=srxqmin & lfxgv<=srxqmax;
     ygv = lfygv>=sryqmin & lfygv<=sryqmax;
     
+    % Replace missing values by zero (to avoid interpolation errors)
     SRLab{lab}(isnan(SRLab{lab})) = 0;
     SRDisp{lab}(isnan(SRDisp{lab})) = 0;
     SRCol{lab}(isnan(SRCol{lab})) = 0;
-        
+    
+    % Compute super-ray view by view
     for v = 1:numView
-        fprintf(repmat('\b',1,sizeMsg-2));
-        msg = ['Constructing super-ray ',num2str(lab),'/',num2str(numLab)...
-           '\nView ',num2str(v),'/',num2str(numView),'\n'];
-        sizeMsg = numel(msg);
-        fprintf(msg);
+        msg = ['Constructing super-ray ',num2str(lab),'/',num2str(numLab),...
+            '\nView ',num2str(v),'/',num2str(numView),'\n'];
+        progress(msg);
+        
+        % Compute missing reference data mask
+        NaNMask = ~(isnan(SRX{lab}(:,:,v)) | isnan(SRY{lab}(:,:,v)));
+        
+        % Remove missing reference data
+        SRXSub = SRX{lab}(:,:,v);
+        SRXSub = reshape(SRXSub(NaNMask),[],1);
+        SRYSub = SRY{lab}(:,:,v);
+        SRYSub = reshape(SRYSub(NaNMask),[],1);
+        SRLabSub = SRLab{lab}(:,:,v);
+        SRLabSub = reshape(SRLabSub(NaNMask),[],1);
+        SRDispSub = SRDisp{lab}(:,:,v);
+        SRDispSub = reshape(SRDispSub(NaNMask),[],1);
         
-        LFLabInt  = griddata(SRYq{lab}(:,:,v),SRXq{lab}(:,:,v),...
-            SRLab{lab}(:,:,1,v),LFYq(xgv,ygv,v),LFXq(xgv,ygv,v),'nearest');
+        % Interpolate label
+        LabInt = griddata(SRXSub,SRYSub,SRLabSub,LFXq(xgv,ygv),LFYq(xgv,ygv),'nearest');
         
-        LFDispInt = griddata(SRYq{lab}(:,:,v),SRXq{lab}(:,:,v),...
-            SRDisp{lab}(:,:,1,v),LFYq(xgv,ygv,v),LFXq(xgv,ygv,v),Method);
+        %Interpolate disparity
+        DispInt = griddata(SRXSub,SRYSub,SRDispSub,LFXq(xgv,ygv),LFYq(xgv,ygv),Method);
         
-        Mask = LFDispInt>LFDisp(xgv,ygv,1,v) & LFLabInt==lab;
+        % Combined depth + stencil buffer
+        Mask = DispInt>LFDisp(xgv,ygv,1,v) & LabInt==lab;
         
-        LFLabTemp = Label(xgv,ygv,1,v);
-        LFLabTemp(Mask) = lab;
-        Label(xgv,ygv,1,v) = LFLabTemp;
+        % Update light field label
+        LabTemp = Label(xgv,ygv,1,v);
+        LabTemp(Mask) = lab;
+        Label(xgv,ygv,1,v) = LabTemp;
         
-        LFDispTemp = LFDisp(xgv,ygv,1,v);
-        LFDispTemp(Mask) = LFDispInt(Mask);
-        LFDisp(xgv,ygv,1,v) = LFDispTemp;
+        % Update light field disparity
+        DispTemp = LFDisp(xgv,ygv,1,v);
+        DispTemp(Mask) = DispInt(Mask);
+        LFDisp(xgv,ygv,1,v) = DispTemp;
         
         for c = 1:numChan
-            LFColInt  = griddata(SRYq{lab}(:,:,v),SRXq{lab}(:,:,v),...
-                SRCol{lab}(:,:,c,v),LFYq(xgv,ygv,v),LFXq(xgv,ygv,v),Method);
+            % Remove missing reference data
+            SRColSub = SRCol{lab}(:,:,c,v);
+            SRColSub = reshape(SRColSub(NaNMask),[],1);
             
+            %Interpolate disparity
+            LFColInt = griddata(SRXSub,SRYSub,SRColSub,LFXq(xgv,ygv),LFYq(xgv,ygv),Method);
+            
+            % Update light field color
             LFColTemp = Color(xgv,ygv,c,v);
             LFColTemp(Mask) = LFColInt(Mask);
             Color(xgv,ygv,c,v) = LFColTemp;
@@ -86,4 +129,13 @@ fprintf('\n\n');
 
 LFSet = SR.FieldsToSet(ImgSize,ImgRes,Offset,Color,Label);
 
+    function progress(msg,varargin)
+        persistent sz
+        if isempty(sz); sz = 0; end
+        if nargin>1; sz = varargin{1}; end
+        
+        fprintf(repmat('\b',1,sz));
+        fprintf(msg);
+        sz = numel(msg)-count(msg,'\n');
+    end
 end
\ No newline at end of file