diff --git a/getLF.m b/getLF.m
index 5800f23f27cb392c98821357a3630d537d425e2c..2b2ebd8aa5ac5767fc74d8cbf899dc0dca1f273e 100644
--- a/getLF.m
+++ b/getLF.m
@@ -1,6 +1,6 @@
 function [LFRef,LFDisp] = getLF(SRRef,SRDisp,SRX,SRY,lfxgv,lfygv,Method)
-%GETLF Builds a super-ray set from a light-field
-%   SRTOLF(LFSet)
+%GETLF Summary of this function goes here
+%   Detailed explanation goes here
 
 %%
 SRCol = {SRRef.Color};
@@ -8,7 +8,6 @@ SRLab = {SRRef.Label};
 
 seh = strel('arbitrary',[1,0,1]);
 sev = strel('arbitrary',[1,0,1]');
-epsilon = 0;
 
 %% Initialize constants and interpolants
 numChan = size(SRCol{1},3);
@@ -34,7 +33,6 @@ warning('off','MATLAB:scatteredInterpolant:DupPtsAvValuesWarnId');
 % Compute light field super-ray by super-ray
 for lab = 1:numLab
     % Replace missing values by zero (to avoid interpolation errors)
-    
     for it = 1:2
         SRCol{lab} = utils.fillmissing(SRCol{lab},seh,1);
         SRDisp{lab} = utils.fillmissing(SRDisp{lab},seh,1);
@@ -43,15 +41,15 @@ for lab = 1:numLab
         SRDisp{lab} = utils.fillmissing(SRDisp{lab},sev,1);
     end
     
-    SRCol{lab}(isnan(SRCol{lab})) = 0;
-    SRLab{lab}(isnan(SRLab{lab})) = 0;
+    SRCol {lab}(isnan(SRCol {lab})) = 0;
+    SRLab {lab}(isnan(SRLab {lab})) = 0;
     SRDisp{lab}(isnan(SRDisp{lab})) = 0;
     
     % Compute super-ray view by view
     for v = 1:numView
-        msg = ['Constructing light field from super-ray ',num2str(lab),'/',num2str(numLab),...
+        msg = ['Computing light field from super-ray ',num2str(lab),'/',num2str(numLab),...
             '\nView ',num2str(v),'/',num2str(numView),'\n'];
-        progress(msg);
+        %progress(msg);
         
         SRXsub    = SRX   {lab}(:,:,1,v);
         SRYsub    = SRY   {lab}(:,:,1,v);
@@ -96,7 +94,7 @@ for lab = 1:numLab
         LabNaN = isnan(interp2(SRYsub,SRXsub,SRLabsub,LFYsub,LFXsub,'cubic'));
         
         % Combined depth + stencil buffer
-        depthStencilMask = tempDisp>(LFDisp(xgv,ygv,1,v)-epsilon) & tempLab==lab;
+        depthStencilMask = tempDisp>LFDisp(xgv,ygv,1,v) & tempLab==lab;
         combinedMask = ~(queryNaN|LabNaN)&depthStencilMask;
         
         buffDisp = LFDisp(xgv,ygv,1,v);
diff --git a/getLF_parfor.m b/getLF_parfor.m
new file mode 100644
index 0000000000000000000000000000000000000000..b35c271bb40993034f69102e03c8cafbdd76a8ef
--- /dev/null
+++ b/getLF_parfor.m
@@ -0,0 +1,147 @@
+function [LFRef,LFDisp] = getLF_parfor(SRRef,SRDisp,SRX,SRY,lfxgv,lfygv,Method)
+%GETLF_PARFOR Summary of this function goes here
+%   Detailed explanation goes here
+
+%%
+SRCol = {SRRef.Color};
+SRLab = {SRRef.Label};
+
+seh = strel('arbitrary',[1,0,1]);
+sev = strel('arbitrary',[1,0,1]');
+
+%% Initialize constants and interpolants
+numChan = size(SRCol{1},3);
+ImgSize = [numel(lfxgv),numel(lfygv),numChan];
+ImgRes = size(SRRef(1).Color);
+ImgRes = ImgRes(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);
+
+wgriddata = warning('query','MATLAB:griddata:DuplicateDataPoints');
+wscattered = warning('query','MATLAB:scatteredInterpolant:DupPtsAvValuesWarnId');
+warning('off','MATLAB:griddata:DuplicateDataPoints');
+warning('off','MATLAB:scatteredInterpolant:DupPtsAvValuesWarnId');
+
+%% Compute light field using interpolation
+% Compute light field super-ray by super-ray
+for lab = 1:numLab
+    msg = ['Computing light field from super-ray ',num2str(lab),'/',num2str(numLab),'\n'];
+    progress(msg);
+    
+    % Replace missing values by zero (to avoid interpolation errors)
+    for it = 1:2
+        SRCol{lab} = utils.fillmissing(SRCol{lab},seh,1);
+        SRDisp{lab} = utils.fillmissing(SRDisp{lab},seh,1);
+        
+        SRCol{lab} = utils.fillmissing(SRCol{lab},sev,1);
+        SRDisp{lab} = utils.fillmissing(SRDisp{lab},sev,1);
+    end
+    
+    SRCol{lab}(isnan(SRCol{lab})) = 0;
+    SRLab{lab}(isnan(SRLab{lab})) = 0;
+    SRDisp{lab}(isnan(SRDisp{lab})) = 0;
+    
+    SRCol_lab = SRCol{lab};
+    SRLab_lab = SRLab{lab};
+    SRDisp_lab = SRDisp{lab};
+    SRX_lab = SRX{lab};
+    SRY_lab = SRY{lab};
+    
+    % Compute super-ray view by view
+    parfor v = 1:numView
+        SRXsub    = SRX_lab   (:,:,1,v);
+        SRYsub    = SRY_lab   (:,:,1,v);
+        SRColsub  = SRCol_lab (:,:,:,v);
+        SRLabsub  = SRLab_lab (:,:,1,v);
+        SRDispsub = SRDisp_lab(:,:,1,v);
+        
+        SRColsub  = utils.fillmissing(SRColsub,seh,1);
+        SRLabsub  = utils.fillmissing(SRLabsub,seh,1);
+        SRDispsub = utils.fillmissing(SRDispsub,seh,1);
+        
+        SRColsub  = utils.fillmissing(SRColsub,sev,1);
+        SRLabsub  = utils.fillmissing(SRLabsub,sev,1);
+        SRDispsub = utils.fillmissing(SRDispsub,sev,1);
+        
+        % Compute super-ray boundary
+        srxqmin = ceil (min(SRXsub(:)));
+        srxqmax = floor(max(SRXsub(:)));
+        sryqmin = ceil (min(SRYsub(:)));
+        sryqmax = floor(max(SRYsub(:)));
+        
+        % Reduce number of query points for speed
+        xgv = lfxgv>=srxqmin & lfxgv<=srxqmax;
+        ygv = lfygv>=sryqmin & lfygv<=sryqmax;
+        
+        [LFXsub,LFYsub] = ndgrid(lfxgv(xgv),lfygv(ygv));
+        
+        [SRXsub,SRYsub,LFXsub,LFYsub] = utils.gridCoords(SRXsub,SRYsub,LFXsub,LFYsub);
+        
+        queryNaN = isnan(LFXsub)|isnan(LFYsub);
+        
+        LFXsub(queryNaN) = 1;
+        LFYsub(queryNaN) = 1;
+        
+        % Interpolate disparity
+        tempDisp = interp2(SRYsub,SRXsub,SRDispsub,LFYsub,LFXsub,Method);
+        
+        % Interpolate label
+        tempLab = interp2(SRYsub,SRXsub,SRLabsub,LFYsub,LFXsub,'nearest');
+        
+        % Find out of boundary values in label using cubic interpolation
+        LabNaN = isnan(interp2(SRYsub,SRXsub,SRLabsub,LFYsub,LFXsub,'cubic'));
+        
+        LFDisp_ = LFDisp(:,:,:,v);
+        Label_  = Label (:,:,:,v);
+        Color_  = Color (:,:,:,v);
+        
+        % Combined depth + stencil buffer
+        depthStencilMask = tempDisp>LFDisp_(xgv,ygv,:,:) & tempLab==lab;
+        combinedMask = ~(queryNaN|LabNaN)&depthStencilMask;
+        
+        buffDisp = LFDisp_(xgv,ygv,:,:);
+        buffLab  = Label_ (xgv,ygv,:,:);
+        
+        buffDisp(combinedMask) = tempDisp(combinedMask);
+        buffLab (combinedMask) = tempLab (combinedMask);
+        
+        LFDisp_(xgv,ygv,:,:) = buffDisp;
+        Label_ (xgv,ygv,:,:) = buffLab;
+        
+        % Interpolate color
+        for c = 1:numChan
+            buffCol = Color_(xgv,ygv,c,:);
+            tempCol = interp2(SRYsub,SRXsub,SRColsub(:,:,c),LFYsub,LFXsub,Method);
+            buffCol(combinedMask) = tempCol(combinedMask);
+            Color_(xgv,ygv,c,:) = buffCol;
+        end
+        
+        LFDisp(:,:,:,v) = LFDisp_;
+        Label (:,:,:,v) = Label_;
+        Color (:,:,:,v) = Color_;
+    end
+end
+fprintf('\n\n');
+
+warning(wgriddata.state,'MATLAB:griddata:DuplicateDataPoints');
+warning(wscattered.state,'MATLAB:scatteredInterpolant:DupPtsAvValuesWarnId');
+
+LFRef = SR.FieldsToSet(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
diff --git a/getSR.m b/getSR.m
index 63ab5c5f0250536dcab6673d332560f3b71a3b19..f578603bb2c21f07adefc81552d55713874b20fd 100644
--- a/getSR.m
+++ b/getSR.m
@@ -1,10 +1,11 @@
 function [SRRef,SRDisp] = getSR(LFRef,LFDisp,SRXq,SRYq,Method)
 %GETSR Summary of this function goes here
 %   Detailed explanation goes here
+
 %% Initialize constants and interpolants
 % Lightfield properties
 LFSize = size(LFRef.Color);
-ImgSize = LFSize(1:3);
+ImgSize = LFSize(1:2);
 numChan = LFSize(3);
 ImgRes = LFSize(4:end);
 Offset = LFRef.Offset;
@@ -21,9 +22,8 @@ numLab = numel(SRXq);
 Color = padarray(Color,[1,1,0,0,0],nan,'both');
 Label = padarray(Label,[1,1,0,0,0],nan,'both');
 LFDisp = padarray(LFDisp,[1,1,0,0,0],nan,'both');
-ImgSize(1:2) = ImgSize(1:2)+2;
+ImgSize = ImgSize+2;
 Offset(1:2) = Offset(1:2)-1;
-LFSize = [ImgSize,ImgRes];
 
 seh = strel('arbitrary',[1,0,1]);
 sev = strel('arbitrary',[1,0,1]');
@@ -42,7 +42,7 @@ Label(isnan(Label))=0;
 LFDisp(isnan(LFDisp))=0;
 
 % Compute projected coordinates of reference samples
-gv = arrayfun(@(x) 1:x,LFSize,'UniformOutput',false);
+gv = arrayfun(@(x) 1:x,ImgSize,'UniformOutput',false);
 
 [LFX,LFY] = ndgrid(gv{1:2});
 LFX = LFX + Offset(1);
@@ -58,9 +58,6 @@ warning('off','MATLAB:scatteredInterpolant:DupPtsAvValuesWarnId');
 %% Initialize super-rays
 fprintf('Super-ray computation...\n');
 for lab = 1:numLab
-    msg = ['Initializing super-ray ',num2str(lab),'/',num2str(numLab),'\n'];
-    progress(msg);
-    
     % Compute super-ray grid boundary
     SRXc = SRXq{lab}(:,:,1,centerView(1),centerView(2));
     SRYc = SRYq{lab}(:,:,1,centerView(1),centerView(2));
@@ -74,19 +71,22 @@ for lab = 1:numLab
     SROff{lab} = [srxqmin,sryqmin,0,0,0];
     
     % Initialize super-ray color, label and disparities
-    SRCol {lab} = nan(size(SRXq{lab}));
-    SRCol {lab} = repmat(SRCol{lab},[1,1,numChan]);
-    SRLab {lab} = nan(size(SRXq{lab}));
-    SRDisp{lab} = nan(size(SRXq{lab}));
+    size_lab = size(SRXq{lab});
+    size_lab_col = size_lab;
+    size_lab_col(3) = numChan;
+    
+    SRCol {lab} = nan(size_lab_col);
+    SRLab {lab} = nan(size_lab);
+    SRDisp{lab} = nan(size_lab);
     
     % Compute super-ray view by view
     for v = 1:numView
-        msg = ['Interpolating super-ray ',num2str(lab),'/',num2str(numLab),...
-            '\nView ',num2str(v),'/',num2str(numView),'\n'];
+        msg = ['Computing super-ray ',num2str(lab),'/',num2str(numLab),...
+            ' from light field\nView ',num2str(v),'/',num2str(numView),'\n'];
         progress(msg);
         
-        Xq = SRXq{lab}(:,:,1,v);
-        Yq = SRYq{lab}(:,:,1,v);
+        Xq = SRXq{lab}(:,:,:,v);
+        Yq = SRYq{lab}(:,:,:,v);
         
         % Reduce number of reference points for speed
         Mask = ...true(size(LFX(:,:,v)));
@@ -115,7 +115,7 @@ for lab = 1:numLab
         LFDispsub = LFDisp(mgv{:},v);
         temp = interp2(LFYsub,LFXsub,LFDispsub,Yq,Xq,Method);
         temp(M) = nan;
-        SRDisp{lab}(:,:,1,v) = temp;
+        SRDisp{lab}(:,:,:,v) = temp;
         
         % Interpolate label
         Labelsub = Label(mgv{:},v);
@@ -126,7 +126,7 @@ for lab = 1:numLab
         % Replace out of boundary values in label using cubic interpolation
         NaNLab = isnan(interp2(LFYsub,LFXsub,Labelsub,Yq,Xq,'cubic'));
         Lab(NaNLab) = nan;
-        SRLab{lab}(:,:,1,v) = Lab;
+        SRLab{lab}(:,:,:,v) = Lab;
     end
 end
 fprintf('\n\n');
diff --git a/getSR_parfor.m b/getSR_parfor.m
new file mode 100644
index 0000000000000000000000000000000000000000..03711c20d7569f48e703f069a5f1617196c2ceaa
--- /dev/null
+++ b/getSR_parfor.m
@@ -0,0 +1,157 @@
+function [SRRef,SRDisp] = getSR_parfor(LFRef,LFDisp,SRXq,SRYq,Method)
+%GETSR_PARFOR Summary of this function goes here
+%   Detailed explanation goes here
+
+%% Initialize constants and interpolants
+% Lightfield properties
+LFSize = size(LFRef.Color);
+ImgSize = LFSize(1:2);
+numChan = LFSize(3);
+ImgRes = LFSize(4:end);
+Offset = LFRef.Offset;
+Color  = LFRef.Color;
+Label  = LFRef.Label;
+centerView = floor(ImgRes/2)+1;
+numView = prod(ImgRes);
+numLab = numel(SRXq);
+
+% Initialize super-ray fields
+[SROff,SRCol,SRLab,SRDisp] = deal(cell(numLab,1));
+
+% Pad values
+Color = padarray(Color,[1,1,0,0,0],nan,'both');
+Label = padarray(Label,[1,1,0,0,0],nan,'both');
+LFDisp = padarray(LFDisp,[1,1,0,0,0],nan,'both');
+ImgSize = ImgSize+2;
+Offset(1:2) = Offset(1:2)-1;
+
+seh = strel('arbitrary',[1,0,1]);
+sev = strel('arbitrary',[1,0,1]');
+
+Color = utils.fillmissing(Color,seh,1);
+Label = utils.fillmissing(Label,seh,1);
+LFDisp = utils.fillmissing(LFDisp,seh,1);
+
+Color = utils.fillmissing(Color,sev,1);
+Label = utils.fillmissing(Label,sev,1);
+LFDisp = utils.fillmissing(LFDisp,sev,1);
+
+% Replace missing values by zero (to avoid interpolation errors)
+Color(isnan(Color))=0;
+Label(isnan(Label))=0;
+LFDisp(isnan(LFDisp))=0;
+
+% Compute projected coordinates of reference samples
+gv = arrayfun(@(x) 1:x,ImgSize,'UniformOutput',false);
+
+[LFX,LFY] = ndgrid(gv{1:2});
+LFX = LFX + Offset(1);
+LFY = LFY + Offset(2);
+
+progress('');
+
+wgriddata = warning('query','MATLAB:griddata:DuplicateDataPoints');
+wscattered = warning('query','MATLAB:scatteredInterpolant:DupPtsAvValuesWarnId');
+warning('off','MATLAB:griddata:DuplicateDataPoints');
+warning('off','MATLAB:scatteredInterpolant:DupPtsAvValuesWarnId');
+
+%% Initialize super-rays
+fprintf('Super-ray computation...\n');
+for lab = 1:numLab
+    msg = ['Computing super-ray ',num2str(lab),'/',num2str(numLab),' from light field\n'];
+    progress(msg);
+    
+    % Compute super-ray grid boundary
+    SRXc = SRXq{lab}(:,:,1,centerView(1),centerView(2));
+    SRYc = SRYq{lab}(:,:,1,centerView(1),centerView(2));
+    
+    Mask = ~(isnan(SRXc)|isnan(SRYc));
+    [~,mgv] = utils.tighten(Mask);
+    
+    srxqmin = SRXc(mgv{1}(1),mgv{2}(1)) - mgv{1}(1);
+    sryqmin = SRYc(mgv{1}(1),mgv{2}(1)) - mgv{2}(1);
+    
+    SROff{lab} = [srxqmin,sryqmin,0,0,0];
+    
+    % Initialize super-ray color, label and disparities
+    size_lab = size(SRXq{lab});
+    size_lab_col = size_lab;
+    size_lab_col(3) = numChan;
+    
+    SRCol_lab = nan(size_lab_col);
+    SRLab_lab = nan(size_lab);
+    SRDisp_lab = nan(size_lab);
+    SRXq_lab = SRXq{lab};
+    SRYq_lab = SRYq{lab};
+    
+    % Compute super-ray view by view
+    parfor v = 1:numView
+        Xq = SRXq_lab(:,:,:,v);
+        Yq = SRYq_lab(:,:,:,v);
+        
+        % Reduce number of reference points for speed
+        Mask = ...true(size(LFX(:,:,v)));
+            LFX>=min(Xq(:))-1 & LFX<=max(Xq(:))+1 & ...
+            LFY>=min(Yq(:))-1 & LFY<=max(Yq(:))+1;
+        
+        [~,mgv] = utils.tighten(Mask);
+        
+        LFXsub = LFX(mgv{:});
+        LFYsub = LFY(mgv{:});
+        
+        M = isnan(Xq)|isnan(Yq);
+        
+        Xq(M) = 1;
+        Yq(M) = 1;
+        
+        % Interpolate color
+        for c = 1:numChan
+            Colsub = Color(:,:,c,v);
+            Colsub = Colsub(mgv{:},:,:);
+            temp = interp2(LFYsub,LFXsub,Colsub,Yq,Xq,Method);
+            temp(M) = nan;
+            SRCol_lab(:,:,c,v) = temp;
+        end
+        
+        % Interpolate disparity
+        LFDispsub = LFDisp(:,:,1,v);
+        LFDispsub = LFDispsub(mgv{:},:,:);
+        temp = interp2(LFYsub,LFXsub,LFDispsub,Yq,Xq,Method);
+        temp(M) = nan;
+        SRDisp_lab(:,:,1,v) = temp;
+        
+        % Interpolate label
+        Labelsub = Label(:,:,1,v);
+        Labelsub = Labelsub(mgv{:},:,:);
+        temp = interp2(LFYsub,LFXsub,Labelsub,Yq,Xq,'nearest');
+        temp(M) = nan;
+        Lab = temp;
+        
+        % Replace out of boundary values in label using cubic interpolation
+        NaNLab = isnan(interp2(LFYsub,LFXsub,Labelsub,Yq,Xq,'cubic'));
+        Lab(NaNLab) = nan;
+        SRLab_lab(:,:,1,v) = Lab;
+    end
+    
+    SRCol {lab} = SRCol_lab;
+    SRLab {lab} = SRLab_lab;
+    SRDisp{lab} = SRDisp_lab;
+end
+fprintf('\n\n');
+
+warning(wgriddata.state,'MATLAB:griddata:DuplicateDataPoints');
+warning(wscattered.state,'MATLAB:scatteredInterpolant:DupPtsAvValuesWarnId');
+
+% Create super-ray structure from fields
+SRRef = SR.FieldsToSet(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