diff --git a/getLF.m b/getLF.m
index 559b15a3bc9c722d089af053d4aa832d0d3e16b9..578179e99e83b185a178beef3a104806b9601852 100644
--- a/getLF.m
+++ b/getLF.m
@@ -1,29 +1,30 @@
-function [LFRef,LFDisp] = getLF(SRRef,SRDisp,SRX,SRY,lfxgv,lfygv,Method)
+function [LFRef,LFDisp] = getLF(SRRef,SRDisp,SRX,SRY,lfxgv,lfygv,Method,numLay)
 %GETLF Summary of this function goes here
 %   Detailed explanation goes here
 
-%%
+%% Initialize constants and interpolants
 SRCol = {SRRef.Color};
 SRLab = {SRRef.Label};
+SRSize = size(SRRef(1).Color);
+SRSize(end+1:5) = 1;
 
 seh = strel('arbitrary',[1,0,1]);
 sev = strel('arbitrary',[1,0,1]');
 
-%% Initialize constants and interpolants
+if ~exist('numLay','var'); numLay = 1; end
 numChan = size(SRCol{1},3);
-ImgSize = [numel(lfxgv),numel(lfygv),numChan];
-ImgRes = size(SRRef(1).Color);
-ImgRes = ImgRes(4:end);
+ImgSize = [numel(lfxgv),numel(lfygv)];
+ImgRes = SRSize(4:end);
 Offset = [min(lfxgv),min(lfygv),1,1,1]-1;
 numLab = numel(SRCol);
 numView = prod(ImgRes);
 
-Color = nan([ImgSize(1:2),numChan,ImgRes]);
-Label = zeros([ImgSize(1:2),1,ImgRes]);
-LFDisp = -inf([ImgSize(1:2),1,ImgRes]);
-
-progress('',0);
+%% Reshape variables
+Color  =  nan  ([ImgSize,numChan,numView,numLay]);
+Label  =  zeros([ImgSize,1      ,numView,numLay]);
+LFDisp = -inf  ([ImgSize,1      ,numView,numLay]);
 
+%% Remove warnings
 wgriddata = warning('query','MATLAB:griddata:DuplicateDataPoints');
 wscattered = warning('query','MATLAB:scatteredInterpolant:DupPtsAvValuesWarnId');
 warning('off','MATLAB:griddata:DuplicateDataPoints');
@@ -31,6 +32,7 @@ warning('off','MATLAB:scatteredInterpolant:DupPtsAvValuesWarnId');
 
 %% Compute light field using interpolation
 % Compute light field super-ray by super-ray
+progress('',0);
 for lab = 1:numLab
     % Replace missing values by zero (to avoid interpolation errors)
     for it = 1:2
@@ -93,25 +95,60 @@ for lab = 1:numLab
         % Find out of boundary values in label using cubic interpolation
         LabNaN = isnan(interp2(SRYsub,SRXsub,SRLabsub,LFYsub,LFXsub,'cubic'));
         
-        % Combined depth + stencil buffer
-        depthStencilMask = tempDisp>LFDisp(xgv,ygv,1,v) & tempLab==lab;
-        combinedMask = ~(queryNaN|LabNaN)&depthStencilMask;
+        % Stencil buffer
+        Mask = ~(queryNaN|LabNaN) & tempLab==lab;
+        
+        Mask = reshape(Mask,[],1);
+        
+        buffDisp = LFDisp(xgv,ygv,1,v,:);
+        buffLab  = Label (xgv,ygv,1,v,:);
+        
+        szBuff = size(buffDisp);
         
-        buffDisp = LFDisp(xgv,ygv,1,v);
-        buffLab  = Label (xgv,ygv,1,v);
+        buffDisp = reshape(buffDisp,[],numLay);
+        buffLab  = reshape(buffLab ,[],numLay);
         
-        buffDisp(combinedMask) = tempDisp(combinedMask);
-        buffLab (combinedMask) = tempLab (combinedMask);
+        tempDisp = reshape(tempDisp,[],1);
+        tempLab  = reshape(tempLab ,[],1);
         
-        LFDisp(xgv,ygv,1,v) = buffDisp;
-        Label (xgv,ygv,1,v) = buffLab;
+        jointDisp = [buffDisp,tempDisp];
+        jointLab  = [buffLab ,tempLab ];
+        
+        [~,sig] = sort(jointDisp,2,'descend');
+        
+        sigsz = size(sig);
+        x = ndgrid(1:sigsz(1),1:sigsz(2));
+        sig = sub2ind(sigsz,x,sig);
+        
+        jointDisp = jointDisp(sig);
+        jointLab  = jointLab (sig);
+        
+        buffDisp(Mask,:) = jointDisp(Mask,1:numLay);
+        buffLab (Mask,:) = jointLab (Mask,1:numLay);
+        
+        buffDisp = reshape(buffDisp,szBuff);
+        buffLab  = reshape(buffLab ,szBuff);
+        
+        LFDisp(xgv,ygv,1,v,:) = buffDisp;
+        Label (xgv,ygv,1,v,:) = buffLab;
         
         % Interpolate color
         for c = 1:numChan
-            buffCol = Color(xgv,ygv,c,v);
             tempCol = interp2(SRYsub,SRXsub,SRColsub(:,:,c),LFYsub,LFXsub,Method);
-            buffCol(combinedMask) = tempCol(combinedMask);
-            Color(xgv,ygv,c,v) = buffCol;
+            
+            buffCol = Color(xgv,ygv,c,v,:);
+            
+            buffCol = reshape(buffCol,[],numLay);
+            tempCol = reshape(tempCol,[],1);
+            
+            jointCol = [buffCol,tempCol];
+            jointCol = jointCol(sig);
+            
+            buffCol(Mask,:) = jointCol(Mask,1:numLay);
+            
+            buffCol = reshape(buffCol,szBuff);
+            
+            Color(xgv,ygv,c,v,:) = buffCol;
         end
     end
 end
@@ -119,11 +156,19 @@ end
 progress('',0);
 fprintf('\n');
 
+%% Reset warnings
 warning(wgriddata.state,'MATLAB:griddata:DuplicateDataPoints');
 warning(wscattered.state,'MATLAB:scatteredInterpolant:DupPtsAvValuesWarnId');
 
+%% Reshape variables
+Color  = reshape(Color ,[ImgSize,numChan,ImgRes,numLay]);
+Label  = reshape(Label ,[ImgSize,1      ,ImgRes,numLay]);
+LFDisp = reshape(LFDisp,[ImgSize,1      ,ImgRes,numLay]);
+
+%%
 LFRef = SR.FieldsToSet(Offset,Color,Label);
 
+%%
     function progress(msg,varargin)
         persistent sz
         if isempty(sz); sz = 0; end