diff --git a/QP_main.m b/QP_main.m
new file mode 100644
index 0000000..8426e60
--- /dev/null
+++ b/QP_main.m
@@ -0,0 +1,37 @@
+%% This file contains a basic pipeline for QP retrieval from brightfield stacks
+% - 3D image stack loading
+% - 3D stack preprocessing
+% - processing parameters definition
+% - phase retrival
+% - display of the results
+
+addpath('QP_package')
+%% 3D image stack loading
+[stack,pname,fname] = loadData;
+[Nx,Ny,Nz] = size(stack);
+
+%% 3D stack preprocessing 
+if Nz == 8 % i.e. multiplane data : remove the coregistration mask
+    stack = cropCoregMask(stack); 
+    stack = cropXY(stack,temp_nPix-4); % extra safety crop 
+else
+    stack = cropXY(stack); %% crop the stack to have a square image
+end
+    [Nx,Ny,Nz] = size(stack);
+    
+%% Phase retrieval
+% define optics and processing parameters
+run('setup_phase.m');
+
+% set experimental parameters
+if size(stack,3) == 8 % i.e. MultiPlane data
+    s.optics.dz = 0.35; % [um]
+else
+    s.optics.dz = 0.2; % typical sampling for fixed cells
+end
+
+% compute the phase
+QP = getQP(stack,s);
+
+% custom 3D stack plotting 
+plotStack(QP)
\ No newline at end of file
diff --git a/QP_package/clamp.m b/QP_package/clamp.m
new file mode 100644
index 0000000..def9bcf
--- /dev/null
+++ b/QP_package/clamp.m
@@ -0,0 +1,17 @@
+function val = clamp(val,minval,maxval)
+
+if isempty(minval) % no lower bound
+    mapmax = val > maxval;
+    val(mapmax) = maxval;
+end
+if isempty(maxval) % no upper bound
+    mapmin = val < minval;
+    val(mapmin) = minval;
+end
+if ~isempty(minval) && ~isempty(maxval)
+    mapmax = val > maxval;
+    val(mapmax) = maxval;
+    
+    mapmin = val < minval;
+    val(mapmin) = minval;
+end
\ No newline at end of file
diff --git a/QP_package/cropCoregMask.m b/QP_package/cropCoregMask.m
new file mode 100644
index 0000000..5fbd1cc
--- /dev/null
+++ b/QP_package/cropCoregMask.m
@@ -0,0 +1,10 @@
+function stack = cropCoregMask(stack)
+
+    temp = stack ~= 0;
+    mask = sum(temp,3); mask = (mask == 8); 
+    masky = (sum(mask,2)./Nx > 0.5);
+    [~,ay] = max(masky); [~,by] = max(masky(end:-1:1)); by = length(masky)-by+1;
+    maskx = (sum(mask,1)./Nx > 0.5);
+    [~,ax] = max(maskx); [~,bx] = max(maskx(end:-1:1)); bx = length(maskx)-bx+1;
+    temp_nPix = min(bx-ax+1,by-ay+1);
+    stack = stack(ay:temp_nPix-1+ay,ax:temp_nPix+ax-1,:);
\ No newline at end of file
diff --git a/QP_package/cropXY.m b/QP_package/cropXY.m
new file mode 100644
index 0000000..bc2add0
--- /dev/null
+++ b/QP_package/cropXY.m
@@ -0,0 +1,21 @@
+% This function crop the input 3D stack to have the size [N,N,~]
+% if N > than Nx or Ny, it make sure that the resulting image is square
+function out = cropXY(input,N)
+[Nx,Ny,~] = size(input);
+if nargin < 2; N = max(Nx,Ny)+1;end
+
+if N > max(Nx,Ny) % invalid/ missing input N => square the stack
+    if Nx~=Ny
+        N = min(Nx,Ny);
+        ox = floor((Nx-N)/2); oy = floor((Ny-N)/2);
+        out = input(ox+1:ox+N,oy+1:oy+N,:);
+    else
+        out = input;
+    end
+else
+    ox = floor((Nx-N)/2); oy = floor((Ny-N)/2);
+    out = input(ox+1:ox+N,oy+1:oy+N,:);
+end
+
+end
+
diff --git a/QP_package/getMirroredStack.m b/QP_package/getMirroredStack.m
new file mode 100644
index 0000000..690c788
--- /dev/null
+++ b/QP_package/getMirroredStack.m
@@ -0,0 +1,38 @@
+function [stackM,kx,kz] = getMirroredStack(stack,s)
+
+[Nx,~,Nz] = size(stack);
+
+% compute real space
+x = linspace(-Nx*s.optics.dx/2,Nx*s.optics.dx/2,Nx);
+z = linspace(-Nz*s.optics.dz/2,Nz*s.optics.dz/2,Nz);
+
+temp = stack;
+if s.proc.mirrorZ
+    t = temp;
+    t(:,:,end+1:2*end) = temp(:,:,end:-1:1);
+    temp = t;
+    kz = 2*pi/(max(z)-min(z)).*linspace(-Nz/2,(Nz-1)/2,2*Nz);
+else
+    if mod(Nz,2) % i.e. if Nz is odd
+        kz = 2*pi/(max(z)-min(z)).*linspace(-Nz/2,Nz/2,Nz);
+    else
+        kz = 2*pi/(max(z)-min(z)).*linspace(-Nz/2,Nz/2-1,Nz);
+    end
+end
+
+if s.proc.mirrorX
+    t = temp;
+    t(end+1:2*end,:,:) = temp(end:-1:1,:,:);
+    t(:,end+1:2*end,:) = t(:,end:-1:1,:);
+    temp = t;
+    kx = 2*pi/(max(x)-min(x)).*linspace(-Nx/2,(Nx-1)/2,2*Nx);
+else
+    if mod(Nz,2) % i.e. if Nx is odd
+        kx = 2*pi/(max(x)-min(x)).*linspace(-Nx/2,Nx/2,Nx);
+    else
+        kx = 2*pi/(max(x)-min(x)).*linspace(-Nx/2,Nx/2-1,Nx);
+    end
+end
+
+stackM = temp;
+end
\ No newline at end of file
diff --git a/QP_package/getQP.m b/QP_package/getQP.m
new file mode 100644
index 0000000..c807756
--- /dev/null
+++ b/QP_package/getQP.m
@@ -0,0 +1,65 @@
+% Author : Adrien Descloux
+% Date : 15 Feb 2018
+% Version : 2.0
+
+% Recover the Cross-spectral density from the intensity stack
+
+% The theory supporting this m-file can be found in : 
+% "Descloux, A., et al. "Combined multi-plane phase retrieval and 
+%  super-resolution optical fluctuation imaging for 4D cell microscopy." 
+%  Nature Photonics 12.3 (2018): 165."
+
+% the private parameters are optimized for the MP-sofi setup
+
+% input parameters :
+    % stack : 3D intensity stack
+    % mask : use precomputed mask for high-speed reconstruction
+% output parameters :
+    % QP : Quantitative Phase
+    % mask : return the mask used for the reconstruction
+function [QP,mask] = getQP(stack,s,mask)
+
+if nargin < 3 % if no mask are provided
+% mirror the data and compute adequate Fourier space grid
+[stackM,kx,kz] = getMirroredStack(stack,s);
+
+% compute usefull stuff
+th = asin(s.optics.NA/s.optics.n);
+th_ill = asin(s.optics.NA_ill/s.optics.n);
+k0max = s.optics.n*2*pi/(s.optics.lambda - s.optics.dlambda/2);
+k0min = s.optics.n*2*pi/(s.optics.lambda - s.optics.dlambda/2);
+
+% compute Fourier space grid and the phase mask
+[Kx,Kz] = meshgrid(kx,kz);
+if isempty(s.optics.kzT)
+    mask2D = Kz >= k0max*(1-cos(th_ill));
+else
+    mask2D = Kz >= s.optics.kzT;
+end
+
+if s.proc.applyFourierMask % => compute the CTF mask for extra denoising
+% CTF theory 
+maskCTF = ((Kx-k0max*sin(th_ill)).^2 + (Kz-k0max*cos(th_ill)).^2 <= k0max^2) & ...
+          ((Kx+k0min*sin(th_ill)).^2 + (Kz-k0min).^2 >= k0min.^2) & ...
+          Kx>=0 & ...
+          Kz < k0max*(1-cos(th));
+maskCTF = maskCTF | maskCTF(:,end:-1:1);
+mask2D = mask2D & maskCTF;
+end
+
+% since we assume a circular symetric CTF, we expand the 2Dmask in 3D
+mask = map3D(mask2D);
+
+end
+
+% Cross-Spectral Density calculation
+Ik = fftshift(fftn(fftshift(stackM)));
+
+Gamma = Ik.*mask; % cross-spectral density
+
+csd = ifftshift(ifftn(ifftshift(Gamma)));
+csd = csd(1:size(stack,1),1:size(stack,2),1:size(stack,3)); % remove the mirrored input
+
+QP = angle(csd + mean(stack(:))/s.optics.alpha);
+
+end
\ No newline at end of file
diff --git a/QP_package/linmap.m b/QP_package/linmap.m
new file mode 100644
index 0000000..e98c963
--- /dev/null
+++ b/QP_package/linmap.m
@@ -0,0 +1,15 @@
+function  rsc = linmap(val,valMin,valMax,mapMin,mapMax)
+
+% convert the input value between 0 and 1
+tempVal = (val-valMin)./(valMax-valMin);
+
+% clamp the value between 0 and 1
+map0 = tempVal < 0;
+map1 = tempVal > 1;
+tempVal(map0) = 0;
+tempVal(map1) = 1;
+
+% rescale and return
+rsc = tempVal.*(mapMax-mapMin) + mapMin;
+
+end
\ No newline at end of file
diff --git a/QP_package/loadData.m b/QP_package/loadData.m
new file mode 100644
index 0000000..a2e800d
--- /dev/null
+++ b/QP_package/loadData.m
@@ -0,0 +1,20 @@
+function [im,pname,fname] = loadData(path)
+if nargin < 1
+    [fname,pname] = uigetfile('*.*');
+    path = [pname,filesep,fname];
+end
+    
+
+keepReading = 1; k = 1;
+warning('off')
+while keepReading 
+    try
+        im(:,:,k) = imread(path,k);
+        k = k+1;
+    catch
+        keepReading = 0;
+        disp('Finished reading... ')
+        disp(['Stack size : ',num2str(size(im))])
+    end
+end
+warning('on')
\ No newline at end of file
diff --git a/QP_package/map3D.m b/QP_package/map3D.m
new file mode 100644
index 0000000..3c67bc7
--- /dev/null
+++ b/QP_package/map3D.m
@@ -0,0 +1,31 @@
+% map a kz,klat 2D CTF into its kx,ky,kz 3D counterpart
+
+function out = map3D(in)
+
+[Nz,Nx] = size(in);
+
+temp_x = linspace(-1,1,Nx);temp_z = linspace(-1,1,Nz);
+[tX,tY] = meshgrid(temp_x,temp_x);
+r = sqrt(tX.^2 + tY.^2);
+r(r > max(temp_x)) = max(temp_x);
+
+mapr = (r-min(temp_x))./(max(temp_x)-min(temp_x)).*(length(temp_x)-2);
+mapf = floor(mapr)+1; mapc = ceil(mapr)+1;
+p = mapc-mapr-1;
+
+out = zeros(length(temp_x),length(temp_x),length(temp_z));
+for kk = 1:length(temp_z)
+    
+    CTF1D = in(kk,:);      % extract the 1D psf
+    tempf = CTF1D(mapf);          % map the 1D psf on a polar grid
+    tempc = CTF1D(mapc);          % map the 1D psf on a polar grid
+        
+    % local linear interpolation 
+    temp     = p.*tempf + (1-p).*tempc;
+    
+    out(:,:,kk) = temp;        % store the 2D psf
+end
+
+
+
+end
diff --git a/QP_package/plotStack.m b/QP_package/plotStack.m
new file mode 100644
index 0000000..1aca731
--- /dev/null
+++ b/QP_package/plotStack.m
@@ -0,0 +1,253 @@
+function plotStack(data,figID)
+
+% this function create a basic gui that allows to scroll through the stack
+% in the specified dimension
+dim = 2;
+
+[s1,s2,s3] = size(data);
+dat_min = min(data(:)); dat_max = max(data(:));
+if nargin == 1
+f = figure('visible','off','position',[360 500 400 285],'WindowScrollWheelFcn',@figScroll);
+else
+try close(figID);end
+f = figure(figID);
+set(f,'visible','off','position',[360 500 400 285],'WindowScrollWheelFcn',@figScroll);
+end
+
+
+set(f,'name',inputname(1))
+hslide = uicontrol('style','slider','String','z',...
+            'position',[330 30 30 200],...
+            'Callback',{@slider_Callback});
+set(hslide,'Min',1);set(hslide,'Max',s3);
+set(hslide,'value',round((s3)/2));set(hslide,'SliderStep',[1/(s3) 1/(s3)]);
+addlistener(hslide,'Value','PreSet',@slider_Callback);
+
+htxt = uicontrol('style','edit','string','0',...
+    'position',[330 235 30 15],'Callback',@edit_Callback);
+hcmax = uicontrol('style','edit','string',num2str(dat_max),...
+    'position',[335 250 50 8],'Callback',@hcmax_Callback,'enable','off');
+hcmin = uicontrol('style','edit','string',num2str(dat_min),...
+    'position',[335 260 50 8],'Callback',@hcmax_Callback,'enable','off');
+hcmapbox = uicontrol('style','checkbox','value',1,'string','Colomap auto',...
+                'Callback',@cmap_Callback,'position',[335 268 50 15]);
+            
+hpopcmap = uicontrol('style','popupmenu','string',{'parula','jet','hot','winter',...
+                    'gray','copper','morgenstemning','isolum'},'units','normalized','position',[.1 .85 0.1 0.1],...
+                    'value',1,'Callback',@popcmap_Callback);
+            
+hcbox = uicontrol('style','checkbox','value',0,'string','Scale colormap to 3D data',...
+                'Callback',@cbocx_Callback,'position',[200 260 80 20]);
+
+hpop = uicontrol('style','popupmenu','string',{'X','Y','Z'},'value',2,...
+                'Callback',@pop_Callback,'position',[280 245 30 30]);
+
+                 
+ha = axes('Units','pixels','Position',[50,60,260,185]);
+
+% by default, plot the midle of the data
+switch dim
+    case 1
+        him = imagesc(reshape(data(round(end/2),:,:),[s2 s3]));
+        
+    case 2
+        him = imagesc(reshape(data(:,round(end/2),:),[s1 s3])');
+        title('Sliding along the "Y" direction')
+    case 3
+        him = imagesc(data(:,:,round(end/2)));
+        colorbar
+        title('Sliding along the "Z" direction');
+    otherwise
+        dim  = 3;
+        disp('Display direction set to "z"')
+        him = imagesc(data(:,:,round(end/2)));
+        title('Sliding along the "Z" direction')
+end
+
+colormap(parula)
+colorbar
+caxis('auto')
+
+% define figure handles
+handles.data = data;
+handles.dim = dim;
+handles.s = [s1 s2 s3];
+handles.txt = htxt;
+handles.slide = hslide;
+handles.cbox = hcbox;
+handles.cmapbox = hcmapbox;
+handles.im = him;
+handles.ax = ha;
+handles.dmin = dat_min; handles.dmax = dat_max;
+handles.hcmin = hcmin; handles.hcmax = hcmax;
+handles.pop = hpop;
+handles.popcmap = hpopcmap;
+handles.f = f;
+
+guidata(handles.f,handles)
+
+% activation of the gui
+set(f,'toolbar','figure');
+set(f,'Units','normalized');
+set(ha,'Units','normalized');
+set(hslide,'Units','normalized');
+set(hpop,'Units','normalized');
+set(hcbox,'Units','normalized');
+set(htxt,'Units','normalized');
+set(hcmapbox,'Units','normalized');
+set(hcmax,'Units','normalized');
+set(hcmin,'Units','normalized');
+set(f,'position',[0.0891    0.3417    0.3885    0.5037]);
+set(f,'Visible','on');
+
+end
+
+function cbocx_Callback(source,~)
+%     edit_Callback(source)
+    slider_Callback(source);
+end
+
+function cmap_Callback(source,~)
+h = guidata(gcf);
+val = get(h.cmapbox,'value');
+if val
+    set(h.hcmax,'enable','off')
+    set(h.hcmin,'enable','off')
+else
+	set(h.hcmax,'enable','on')
+    set(h.hcmin,'enable','on')
+end
+slider_Callback(source);
+end
+
+function hcmax_Callback(source,~)
+h = guidata(gcf);
+cmax = str2double(get(h.hcmax,'string'));
+cmin = str2double(get(h.hcmin,'string'));
+cmax = clamp(cmax,h.dmin,h.dmax); cmin = clamp(cmin,h.dmin,h.dmax);
+if cmax <= cmin; cmax = h.dmax;cmin = h.dmin; end
+set(h.hcmax,'string',num2str(cmax))
+set(h.hcmin,'string',num2str(cmin))
+    slider_Callback(source);
+end
+
+function slider_Callback(source,~)
+
+h = guidata(gcf);
+h.dim = get(h.pop,'value');
+
+% update slide properties
+set(h.slide,'max',h.s(h.dim));
+slid_val = get(h.slide,'value');
+if slid_val > h.s(h.dim)
+	set(h.slide,'value',h.s(h.dim))
+end
+
+val = get(h.slide,'value');
+vmax = get(h.slide,'max');vmin = get(h.slide,'min');
+val = round(val);
+set(h.slide,'value',round(val));
+% set(h.txt,'string',num2str(2*(round(val)-1)/(h.s(h.dim)-1)-1));
+set(h.txt,'string',num2str(round(val)))
+
+val = round((vmin + vmax)-val);
+switch get(h.pop,'value')
+    case 1
+        temp = reshape(h.data(val,:,:),[h.s(2) h.s(3)])';
+        set(h.im,'Cdata',temp)
+        set(h.ax,'xlim',[0.5 h.s(2)+0.5],'ylim',[0.5 h.s(3)+0.5])
+        title('Sliding along the "X" direction')
+        xlabel('"Y" direction sdf');ylabel('"Z" direction')
+    case 2
+        temp = reshape(h.data(:,val,:),[h.s(1) h.s(3)])';
+        set(h.im,'Cdata',temp)
+        set(h.ax,'xlim',[0.5 h.s(1)+0.5],'ylim',[0.5 h.s(3)+0.5])
+        title('Sliding along the "Y" direction')
+        xlabel('"X" direction');ylabel('"Z" direction')
+    case 3
+        temp = h.data(:,:,val) ;
+        set(h.im,'Cdata',temp)
+        set(h.ax,'xlim',[0.5 h.s(2)+0.5],'ylim',[0.5 h.s(1)+0.5])
+        title('Sliding along the "Z" direction')
+        xlabel('"Y" direction');ylabel('"X" direction')
+end
+
+try
+    cmapS = get(h.popcmap,'string');
+    colormap(cmapS{get(h.popcmap,'value')})
+end
+
+if get(h.cbox,'value')
+    if (h.dmin-h.dmax)
+        caxis([h.dmin h.dmax])
+    else
+        caxis([h.dmin h.dmax+1])
+    end
+else
+    if get(h.cmapbox,'value')
+        caxis('auto')
+    else
+        cmax = str2double(get(h.hcmax,'string'));
+        cmin = str2double(get(h.hcmin,'string'));
+        caxis([cmin cmax])
+    end
+end
+
+% save global changes
+guidata(h.f,h);
+
+end
+
+function pop_Callback(source,event)
+    h = guidata(gcf);
+    set(h.txt,'string','0');
+    h.dim = get(h.pop,'value');
+    guidata(gcf,h);
+    edit_Callback(source,event);
+%     slider_Callback(source,event);
+end
+
+function popcmap_Callback(source,event)
+
+    edit_Callback(source,event);
+%     slider_Callback(source,event);
+end
+
+function edit_Callback(source,event)
+    h = guidata(gcf);
+    val = str2num(get(h.txt,'string'));
+    if val < 1
+%         val = floor(((val+1)/2)*(h.s(h.dim)-1)+1)
+        val = floor(linmap(val,-1,1,1,h.s(h.dim)));
+    end
+        
+    if val > h.s(h.dim)
+    	val = h.s(h.dim);
+    elseif val < 1
+        val = 1;
+    end
+    
+set(h.slide,'value',val);
+slider_Callback(source,event);
+
+end
+
+function figScroll(source,event)
+
+dir = event.VerticalScrollCount;
+
+h = guidata(gcf);
+val = get(h.slide,'value');
+
+% val = round(val - dir.*0.05.*h.s(h.dim));
+val = val -dir;
+if val > h.s(h.dim)
+    val = h.s(h.dim);
+elseif val < 1
+    val = 1;
+end
+set(h.slide,'value',val);
+slider_Callback(source,event);
+
+end
+
diff --git a/QP_package/setup_phase.m b/QP_package/setup_phase.m
new file mode 100644
index 0000000..640c652
--- /dev/null
+++ b/QP_package/setup_phase.m
@@ -0,0 +1,18 @@
+% run this file to setup the QP reconstruction
+
+% optics parameters
+s.optics.dx = 0.11 ;        % projected pixel size [um]
+s.optics.dz = 0.2 ;        % axial sampling [um]
+s.optics.NA = 1.15;          % Numerical Aperture detection
+s.optics.NA_ill = 0.35 ;     % Numerical Aperture illumination
+s.optics.n = 1.33;          % refractive index
+s.optics.lambda = 0.58 ;    % Central wavelength in vacuum
+s.optics.dlambda = 0.075 ;  % Spectrum bandwidth
+s.optics.alpha = 3.15 ;     % Experimental coeficient for QP "normalisation"
+s.optics.kzT = 0.04 ;        % Axial cutoff frequency. 
+                            % if set to [] use the theoretical one (nk0max*(1-cos(asin(NA_ill))))
+
+% processing parameters
+s.proc.mirrorX = 0 ;            % mirror the input stack along X 
+s.proc.mirrorZ = 1 ;            % mirror the input stack along Z
+s.proc.applyFourierMask = 1 ;   % apply the denoising Fourier mask