diff --git a/pom.xml b/pom.xml
old mode 100644
new mode 100755
index 732bee8..2760901
--- a/pom.xml
+++ b/pom.xml
@@ -14,7 +14,7 @@
 
 	<groupId>org.pasteur</groupId>
 	<artifactId>ZOLA_</artifactId>
-	<version>0.2.7-SNAPSHOT</version>
+	<version>0.2.8-SNAPSHOT</version>
 
 	<name>ZOLA</name>
 	<description>An ImageJ plugin for 3D localization microscopy.</description>
@@ -78,43 +78,43 @@
                 <dependency>
                     <groupId>org.jcuda</groupId>
                     <artifactId>jcuda</artifactId>
-                    <version>0.8.0</version>
+                    <version>11.1.1</version>                
                 </dependency>
                 <dependency>
                         <groupId>org.jcuda</groupId>
                         <artifactId>jcublas</artifactId>
-                        <version>0.8.0</version>
+                        <version>11.1.1</version>
                 </dependency>
                 <dependency>
                         <groupId>org.jcuda</groupId>
                         <artifactId>jcufft</artifactId>
-                        <version>0.8.0</version>
+                        <version>11.1.1</version>
                 </dependency>
                 <dependency>
                         <groupId>org.jcuda</groupId>
                         <artifactId>jcusparse</artifactId>
-                        <version>0.8.0</version>
+                        <version>11.1.1</version>
                 </dependency>
                 <dependency>
                         <groupId>org.jcuda</groupId>
                         <artifactId>jcusolver</artifactId>
-                        <version>0.8.0</version>
+                        <version>11.1.1</version>
                 </dependency>
                 <dependency>
                         <groupId>org.jcuda</groupId>
                         <artifactId>jcurand</artifactId>
-                        <version>0.8.0</version>
-                </dependency>
-                <dependency>
-                        <groupId>org.jcuda</groupId>
-                        <artifactId>jnvgraph</artifactId>
-                        <version>0.8.0</version>
+                        <version>11.1.1</version>
                 </dependency>
                 <dependency>
                         <groupId>org.jcuda</groupId>
                         <artifactId>jcudnn</artifactId>
-                        <version>0.8.0</version>
+                        <version>11.1.1</version>
                 </dependency>
+            <dependency>
+                <groupId>org.apache.commons</groupId>
+                <artifactId>commons-math3</artifactId>
+                <version>3.6.1</version>
+            </dependency>
                 <!--<dependency>
                     <groupId>org.pasteur</groupId>
                     <artifactId>SMLM_IO</artifactId>
diff --git a/src/main/java/org/pasteur/imagej/ZOLA.java b/src/main/java/org/pasteur/imagej/ZOLA.java
old mode 100644
new mode 100755
index 4f62572..48566ee
--- a/src/main/java/org/pasteur/imagej/ZOLA.java
+++ b/src/main/java/org/pasteur/imagej/ZOLA.java
@@ -375,7 +375,7 @@ else if (command.startsWith("zola_filtering")){
         
         
         else if (command.startsWith("zola_frc_fsc")){
-            
+            IJ.log("sorry FRC is not working yet");
             frc_fsc();
             
         }
@@ -456,7 +456,7 @@ else if (command.startsWith("zola_test3")){
         }
         else if (command.startsWith("zola_test4")){
             
-            IJ.log("no code for this function");
+            test4();
             
             
         }
@@ -609,6 +609,9 @@ else if (command.startsWith("zola_phaseoptim")){
         else if (command.startsWith("zola_imageFromTable")){
             simulationFromLocTable();
         }
+        else if (command.startsWith("zola_PSF_Generator")){
+            psfGenerator();
+        }
         else if (command.startsWith("zola_getModels_GPUalloc")){
             getModels_GPUalloc();
         }
@@ -673,6 +676,17 @@ else if (command.startsWith("zola_photonCountEstimation")){
             photonCountEstimation();
         }
         
+        else if (command.startsWith("zola_SCMOS_offset_variance_map")){
+            
+            scmosComputeOffsetAndVariance();
+        }
+        else if (command.startsWith("zola_SCMOS_gain_map")){
+            
+            scmosComputeGain();
+        }
+        else if (command.startsWith("zola_convertImageToPhotonCount")){
+            convertImageToPhotonCount();
+        }
         else if (command.startsWith("zola_crlb_FromFileDualObj")){
             
             crlbfromfileDualObj();
@@ -697,6 +711,7 @@ else if (command.startsWith("zola_crlb")){
             
             
         }
+        
         else if (command.startsWith("zola_photonConversion")){
             
             IJ.log("not implemented");
@@ -1100,7 +1115,7 @@ void localization(){
                 
                 
                 
-                if (gain<=0||adu<=0){
+                if (!isSCMOS && (gain<=0||adu<=0)){
                     new WaitForUserDialog("Error message", "gain has to be positive (default value = 1)").show();
                 }
                 else{
@@ -1443,7 +1458,7 @@ void localizationHD(){
                 
                 
                 
-                if (gain<=0||adu<=0){
+                if (!isSCMOS && (gain<=0||adu<=0)){
                     new WaitForUserDialog("Error message", "gain has to be positive (default value = 1)").show();
                 }
                 else{
@@ -1714,7 +1729,7 @@ void localizationtest(){
                 
                 
                 
-                if (gain<=0||adu<=0){
+                if (!isSCMOS && (gain<=0||adu<=0)){
                     new WaitForUserDialog("Error message", "gain has to be positive (default value = 1)").show();
                 }
                 else{
@@ -1765,6 +1780,9 @@ void localizationtest(){
     
     
     
+    
+    
+    
     void crlb(){
         
         if ((processing !=null)&&(processing.booleanValue()==true)){
@@ -2466,7 +2484,7 @@ void calibration(){
                 return;
             }
             
-            sigma=0.9;
+            sigma=1;
             
             double [][][] image =new double [nbImage][width][height];
             
@@ -2642,7 +2660,7 @@ void calibration(){
             
             
             
-            if (gain<=0||adu<=0){
+            if (!isSCMOS && (gain<=0||adu<=0)){
                 new WaitForUserDialog("Error message", "gain has to be positive (default value = 1)").show();
             }
             else{
@@ -2715,6 +2733,214 @@ void calibration(){
     
     
     
+    
+    
+        
+    void psfGenerator(){
+        
+        
+        
+        
+        if ((processing !=null)&&(processing.booleanValue()==true)){
+            printErrorMessage();
+            
+            return;
+        }
+        
+        processing=new Boolean(true);
+        
+        
+        
+        
+        {
+            
+            String pathdir=IJ.getDirectory("current");
+            String path=pathdir;
+            
+//            
+            path=path+"ZOLA_generated_PSF.json";
+            
+            
+            
+            
+            
+            
+            //classes.ImageShow.imshow(image,"image input");
+            //String [] fieldZernike = {"6","10","15","21","28","36","45"};
+            //String [] fieldZernike = {"pixel-based (prototype)","Zernike 15 coefs","Zernike 28 coefs","Zernike 45 coefs","Zernike 66 coefs"};
+            String [] fieldShape = {"tophat","donut"};
+            int [] correspondingFieldShape = {0,1};
+            
+            
+            GenericDialog gd = new GenericDialog("ZOLA: PSF generator");
+
+            
+            
+                
+            Font font = gd.getFont();
+            Font fontBold= gd.getFont();
+            try{
+                fontBold=new Font(font.getName(),Font.BOLD,font.getSize());
+            }catch(Exception e){}
+            
+            gd.addCheckbox("Run_on_GPU", GPU_computation);
+            
+            
+            gd.addNumericField("Pixel_size:", xystep*1000, 1,6,"(nm)");
+            
+            
+            gd.addMessage("Sample parameters",fontBold);
+            
+            gd.addNumericField("Z_step: ", zstep*1000, 1,6,"(nm)");
+            
+            gd.addNumericField("Axial_range: ", axialRange,2,6,"(µm)");
+            
+            gd.addNumericField("Numerical_aperture: ", na, 3);
+            gd.addNumericField("immersion_refractive index: ", noil, 3,6,"~1.518 for oil, ~1.33 for water");
+            
+            gd.addNumericField("Mounting_medium_refractive index: ", nwat, 3,6,"1.33 for water");
+            
+            
+            gd.addNumericField("Distance_focus_to_coverslip (µm)", zfocus, 3);
+            
+            gd.addNumericField("wavelength: ", wavelength, 3,6,"(µm)");
+            
+            
+            
+        
+        
+        
+        
+        
+        
+
+            gd.addMessage("Advanced parameters",fontBold);
+            
+            
+            gd.addNumericField("Patch_size: ", sizePatchPhaseRet, 0,6,"(pixels)");
+            //gd.addSlider("Zernike_coefficient #", 0, 2*((int)Math.sqrt(nbImage)), (int)Math.sqrt(nbImage));
+            gd.addCheckbox("Use_pupil_apodization factor", withApoFactor);
+            gd.addChoice("Zernike_coefficient number", fieldShape,  fieldShape[0]);
+            
+
+            gd.addStringField("Result_calibration_file:",path,sizeTextString); 
+            
+        int shape=0;
+        TextField textField2; 
+        // Add a mouse listener to the config file field 
+        if (!(java.awt.GraphicsEnvironment.isHeadless() || IJ.isMacro())) 
+        {
+            int t=0;
+            
+            Vector<TextField> texts = (Vector<TextField>) gd.getStringFields(); 
+            
+            
+            textField2 = texts.get(t++); 
+            MouseOptionSave mos=new MouseOptionSave(textField2,pathdir,"Export calibration file");
+            textField2.addMouseListener(mos); 
+            
+            
+        }
+            
+            gd.showDialog();
+            
+            
+            if (!gd.wasCanceled()){
+                
+                GPU_computation = (boolean)gd.getNextBoolean();
+                
+                xystep = (double)gd.getNextNumber()/1000;
+                zstep = (double)gd.getNextNumber()/1000;
+                axialRange = (double)gd.getNextNumber();
+                
+                
+                na = (double)gd.getNextNumber();
+                noil = (double)gd.getNextNumber();
+                nwat = (double)gd.getNextNumber();
+                zfocus = (double)gd.getNextNumber();
+                        
+                wavelength = (double)gd.getNextNumber();
+                
+                sizePatchPhaseRet = (int)gd.getNextNumber();
+                withApoFactor=(boolean)gd.getNextBoolean();
+                shape=gd.getNextChoiceIndex();
+                
+                //sigma = (double)gd.getNextNumber();
+                
+                
+                path_calibration = gd.getNextString(); 
+            }
+            else{
+                processing=new Boolean(false);
+                return;
+            }
+            
+            int nnn=(path_calibration.lastIndexOf("."));
+            //IJ.log("path_calibration "+path_calibration+"  "+nnn+"  "+path_calibration.length());
+            if ((nnn>=0)&&(nnn>path_calibration.length()-6)){
+                path_calibration=path_calibration.substring(0,nnn);
+            }
+            
+            path_calibration+=".json";
+            //IJ.log("path_calibration "+path_calibration);
+            if (sizePatchPhaseRet%2==1){
+                sizePatchPhaseRet++;
+            }
+            
+            int sizeFFT=180;
+            sizeFFT=Math.max(sizeFFT,sizePatchPhaseRet*2);
+            sizeFFT=Math.min(sizeFFT,maxFFT);
+            
+            
+            
+            {
+                IJ.log("ZOLA PSF generator started");
+                
+                
+                long t=System.currentTimeMillis();
+                
+                if (GPU_computation){
+                    //maybe test if image is a square
+                    //IJ.log("init...");
+                    MyCudaStream.init(1);
+                    //IJ.log("init cuda ok "+xp+"  "+yp+"  "+sizePatchPhaseRet+"  "+rescale_slope+"  "+rescale_intercept+"  "+zstep+"  "+image);
+                    IJ.log("todo");
+                    //IJ.log("param ok");
+                    
+                    org.pasteur.imagej.process.gpu.PSFgenerator_ psfgen = new org.pasteur.imagej.process.gpu.PSFgenerator_(sizeFFT,xystep,zstep,wavelength,noil,na,path_calibration,sigma,axialside,withApoFactor,axialRange,zstep);
+                    
+                    psfgen.run();
+                    
+                    psfgen.free();
+                            
+                    MyCudaStream.destroy();
+                }
+                else{
+                    
+                    
+                    IJ.log("todo");
+                    
+                }
+                long t2=System.currentTimeMillis();
+                IJ.log("elapsed time "+Math.floor((double)(t2-t)/(double)(1000*60))+" min "+Math.floor(60*((((double)(t2-t)/(double)(1000*60)))-Math.floor((double)(t2-t)/(double)(1000*60))))+" sec");
+            }
+
+            
+            
+        }
+        
+        processing=new Boolean(false);
+        
+        
+        
+        
+        
+    }
+    
+    
+    
+    
+    
     void wobblecorrection(){
         if ((processing !=null)&&(processing.booleanValue()==true)){
             printErrorMessage();
@@ -3025,7 +3251,7 @@ void wobblecalibration(){
                 
                 
                 
-                if (gain<=0||adu<=0){
+                if (!isSCMOS && (gain<=0||adu<=0)){
                     new WaitForUserDialog("Error message", "gain has to be positive (default value = 1)").show();
                 }
                 else{
@@ -5447,6 +5673,8 @@ void simulationFromLocTable(){
             
         double bckg=0;
         
+        int imagesize=-1;
+        
         gd.addStringField("Calibration_path:",path_calibration ,sizeTextString);
         
         gd.addStringField("Localization_path:",path_localization ,sizeTextString); 
@@ -5462,6 +5690,11 @@ void simulationFromLocTable(){
         gd.addNumericField("Background: ", bckg, 3,6,"");
         
         
+        gd.addNumericField("Image size: ", imagesize, 3,6,"-1 for automatic sizing");
+        
+        
+        
+        
         gd.addStringField("Result_image_path:",path_result ,sizeTextString); 
         
         // Add a mouse listener to the config file field 
@@ -5495,6 +5728,7 @@ void simulationFromLocTable(){
             nwat = gd.getNextNumber();
             this.zfocus=gd.getNextNumber();
             bckg = gd.getNextNumber();
+            imagesize = (int)gd.getNextNumber();
             path_result = gd.getNextString(); 
         }
         else{
@@ -5540,7 +5774,7 @@ void simulationFromLocTable(){
         
         mg=new org.pasteur.imagej.process.gpu.Model_generator_(dp);
         mg.setBackground(bckg);
-        mg.computePSFfromTable(sl, path_result,startAt0);
+        mg.computePSFfromTable(sl, path_result,startAt0,imagesize);
         mg.free();
         
         processing=new Boolean(false);
@@ -8766,7 +9000,7 @@ void multiEmitterLocalization(){
                 
                 
                 
-                if (gain<=0||adu<=0){
+                if (!isSCMOS && (gain<=0||adu<=0)){
                     new WaitForUserDialog("Error message", "gain has to be positive (default value = 1)").show();
                 }
                 else{
@@ -9566,16 +9800,18 @@ else if (i>merge){
     
     
     
-    
-    
-    
-    void photonCountEstimation(){
+        
+                
+                
+                
+                
+    void scmosComputeGain(){
         
         
         
         
             
-        ImagePlus imp=IJ.getImage();
+        String [] images=WindowManager.getImageTitles();
 
 
         if ((processing !=null)&&(processing.booleanValue()==true)){
@@ -9583,41 +9819,752 @@ void photonCountEstimation(){
             return;
         }
         
-        processing=new Boolean(true);
-
-
-        if (imp==null){
-            IJ.log("no image opened");
-            processing=new Boolean(false);
-        }
-        /*else if ((imp.getWidth()>180)||(imp.getHeight()>180)){
-            IJ.log("image should be < 180 pixels for the moment");
-        }*/
-        else{
-
-            GenericDialog gd=null;
-
-            //since NonBlockingGenericDialog incompatible with headless mode:
-            if (!(java.awt.GraphicsEnvironment.isHeadless() || IJ.isMacro())){
-                gd = new NonBlockingGenericDialog("ZOLA: photon count estimation");
+            processing=new Boolean(true);
+            
+            
+            if (images.length==0){
+                IJ.log("no image opened");
+                processing=new Boolean(false);
             }
+            /*else if ((imp.getWidth()>180)||(imp.getHeight()>180)){
+                IJ.log("image should be < 180 pixels for the moment");
+            }*/
             else{
-                gd = new GenericDialog("ZOLA: photon count estimation");
-            }
-
-                Font font = gd.getFont();
-                Font fontBold= gd.getFont();
-                try{
-                    fontBold=new Font(font.getName(),Font.BOLD,font.getSize());
-                }catch(Exception e){}
-
-                gd.addMessage("This function estimate ADU and offset to convert pixels in photon count",fontBold);
-                gd.addMessage("Please, import a stack of images (~100 images) containing non moving objects (such as beads)",fontBold);
                 
-                gd.addNumericField("Gain: ", this.gain, 3,6,"");
-
                 
-                gd.showDialog();
+                
+                
+                
+                
+                
+                String pathdir=IJ.getDirectory("image");
+                if (pathdir==null){
+                    pathdir=IJ.getDirectory("current");
+                }
+                String path=pathdir;
+                
+                path=path+"ZOLA_SCMOS_gain.csv";
+                
+                
+                
+                
+                GenericDialog gd=null;
+                
+                //since NonBlockingGenericDialog incompatible with headless mode:
+                if (!(java.awt.GraphicsEnvironment.isHeadless() || IJ.isMacro())){
+                    gd = new NonBlockingGenericDialog("ZOLA: SCMOS calibration 2/2");
+                }
+                else{
+                    gd = new GenericDialog("ZOLA: SCMOS calibration 2/2");
+                }
+
+                    Font font = gd.getFont();
+                    Font fontBold= gd.getFont();
+                try{
+                    fontBold=new Font(font.getName(),Font.BOLD,font.getSize());
+                }catch(Exception e){}
+
+
+                    for (int i=0;i<images.length;i++){
+                        gd.addMessage("input image # "+i+" : "+images[i]);
+                    }
+
+
+                    gd.addMessage("Already computed maps:", fontBold);
+                    
+                    gd.addStringField("SCMOS offset map:",path_SCMOSoffset,sizeTextString); 
+                    
+                    gd.addStringField("SCMOS variance map:",path_SCMOSvariance,sizeTextString); 
+                    
+                    gd.addMessage("Remaining map to compute:", fontBold);
+                    
+                    gd.addStringField("SCMOS gain map:",path,sizeTextString); 
+
+
+
+
+                    TextField textField; 
+                    TextField textField2; 
+                    TextField textField3; 
+                    // Add a mouse listener to the config file field 
+                    if (!(java.awt.GraphicsEnvironment.isHeadless() || IJ.isMacro())) 
+                    {
+                        int t=0;
+                        //textField
+                        Vector<TextField> texts = (Vector<TextField>) gd.getStringFields(); 
+
+                        textField = texts.get(t++); 
+                        MouseOptionLoad mol=new MouseOptionLoad(textField,path_SCMOSoffset,"Import SCMOS offset map");
+                        textField.addMouseListener(mol); 
+                        
+                        textField2 = texts.get(t++); 
+                         mol=new MouseOptionLoad(textField2,path_SCMOSvariance,"Import SCMOS variance map");
+                        textField2.addMouseListener(mol); 
+
+                        textField3 = texts.get(t++); 
+                        MouseOptionSave mos=new MouseOptionSave(textField3,path,"Export localization result table");
+                        textField3.addMouseListener(mos); 
+
+
+                    }
+
+
+                    gd.showDialog();
+
+
+
+
+
+                    if (!gd.wasCanceled()){
+                        
+                        path_SCMOSoffset = gd.getNextString(); 
+                        path_SCMOSvariance = gd.getNextString(); 
+
+                        
+                        path_SCMOSgain = gd.getNextString(); 
+                    }
+                    else{
+                        processing=new Boolean(false);
+                        return;
+                    }
+                
+                
+                
+                int nnn=(path_SCMOSgain.lastIndexOf("."));
+                if ((nnn>=0)&&(nnn>path_SCMOSgain.length()-6)){
+                    path_SCMOSgain=path_SCMOSgain.substring(0,nnn);
+                }
+                path_SCMOSgain+=".tif";
+            
+                
+                ImagePlus [] imp = new ImagePlus[images.length];
+                
+                for (int i=0;i<images.length;i++){
+                    IJ.selectWindow(images[i]);
+                    imp[i]=IJ.getImage();
+                    
+                }
+                
+                
+                int width=imp[0].getWidth();
+                int height=imp[0].getHeight();
+                
+                for (int i=1;i<images.length;i++){
+                    if ((imp[i].getWidth()!=width)||(imp[i].getHeight()!=height)){
+                        IJ.log("ERROR: sizes in opened image stacks are not the same");processing=new Boolean(false);return;
+                    }
+                }
+                
+            
+                int xstart=0;
+                int ystart=0;
+                
+                
+                
+                
+                
+                
+            
+            
+                
+                
+                
+                
+                
+                
+                ImagePlus impoffset=null;
+                try{
+                    impoffset=new ImagePlus(path_SCMOSoffset);
+                }catch(Exception ee){IJ.log("impossible to load offset map "+ee);processing=new Boolean(false);return;}
+                
+                ImagePlus impvariance=null;
+                try{
+                    impvariance=new ImagePlus(path_SCMOSvariance);
+                }catch(Exception eee){IJ.log("impossible to load variance map "+eee);processing=new Boolean(false);return;}
+                
+                
+                if ((impoffset.getWidth()!=width)||(impoffset.getHeight()!=height)){
+                    IJ.log("ERROR: input offset map size is different than opened image size");processing=new Boolean(false);return;
+                }
+                if ((impvariance.getWidth()!=width)||(impvariance.getHeight()!=height)){
+                    IJ.log("ERROR: input variance map size is different than opened image size");processing=new Boolean(false);return;
+                }
+                
+                
+                
+                
+                
+                
+                
+                
+                
+                ImageProcessor ipgain= new FloatProcessor(width,height);
+                ImageProcessor ipoffset= impoffset.getProcessor();
+                ImageProcessor ipvariance= impvariance.getProcessor();
+                float [][][][] pix= new float[images.length][][][];
+                float [][][] var = new float [width][height][images.length];
+                float [][][] mean = new float [width][height][images.length];
+                
+                for (int k=0;k<images.length;k++){
+                    IJ.showProgress(0.4*(double)k/(double)(images.length));
+                    ImageProcessor ip;
+                    boolean stacked=true;
+                    int nbImage=imp[k].getNSlices();
+                    if (imp[k].getNSlices()==1){
+                        nbImage=imp[k].getNFrames();
+                        stacked=false;
+                    }
+                    ImageStack [] ims=null;
+                    if (stacked){
+                        ims = new ImageStack[1];
+                        ims[0]=imp[k].getStack();
+
+                    }
+                    ArrayList<float [][]> al = new ArrayList<float [][]>();
+                    
+                    pix[k]= new float[width][height][nbImage];
+                    //ImageStack imsres = new ImageStack(width,height);
+                    int index=0;
+                    
+                    //copy image stack in pix vector
+                    for (int i=0;i<nbImage;i++){
+                        
+                        if (stacked){
+                            ip = ims[0].getProcessor(i+1);
+                        }
+                        else{
+                            imp[k].setSlice(i+1);
+                            ip = imp[k].getProcessor();
+                        }
+
+                        //fill tmp stack at the beginning
+
+                        for (int u=0;u<width;u++){
+                            for (int uu=0;uu<height;uu++){
+                                float v=ip.getPixelValue(u+xstart, uu+ystart);
+                                pix[k][u][uu][i]=v;
+                            }
+                        }
+                    }
+                }
+                
+                for (int k=0;k<images.length;k++){
+                    IJ.showProgress(.4+0.4*(double)k/(double)(images.length));
+                    for (int i=0;i<width;i++){
+                        for (int ii=0;ii<height;ii++){
+                            for (int j=0;j<pix[k][i][ii].length;j++){
+                                mean[i][ii][k]+=pix[k][i][ii][j];
+                            }
+                            mean[i][ii][k]/=pix[k][i][ii].length;
+                            
+                            for (int j=0;j<pix[k][i][ii].length;j++){
+                                var[i][ii][k]+=(pix[k][i][ii][j]-mean[i][ii][k])*(pix[k][i][ii][j]-mean[i][ii][k]);
+                            }
+                            var[i][ii][k]/=pix[k][i][ii].length;
+                            
+                        }
+                    }  
+                }
+                
+                
+                for (int i=0;i<width;i++){
+                    IJ.showProgress(.8+0.2*(double)i/(double)(width));
+                    for (int ii=0;ii<height;ii++){
+                        float avgGain=0;
+                        for (int k=0;k<images.length;k++){
+                            avgGain+=(var[i][ii][k]-ipvariance.getPixelValue(i, ii))/(mean[i][ii][k]-ipoffset.getPixelValue(i, ii));
+                        }
+                        avgGain/=images.length;
+                        ipgain.putPixelValue(i, ii, avgGain);
+                    }
+                }
+                
+                
+                
+                ImagePlus impgain = new ImagePlus("SCMOS_gain",ipgain);
+                FileSaver fs_gain = new FileSaver(impgain);
+                fs_gain.saveAsTiff(path_SCMOSgain);
+                
+                
+                IJ.log("SCMOS calibration part 2/2 finished.");
+                IJ.showProgress(1);
+                
+                isSCMOS=true;
+            }
+        
+        processing=new Boolean(false);
+        
+        
+    }
+    
+    
+                
+    void convertImageToPhotonCount(){
+        
+        ImagePlus imp = IJ.getImage();
+        
+        if ((processing !=null)&&(processing.booleanValue()==true)){
+            printErrorMessage();
+            
+            return;
+        }
+        
+        processing=new Boolean(true);
+        
+        
+        
+        
+        if (imp==null){
+            IJ.log("no image opened");
+        }
+        else{
+            
+            
+            
+            boolean stacked=true;
+            int nbImage=imp.getNSlices();
+            if (nbImage==1){
+                nbImage=imp.getNFrames();
+                stacked=false;
+            }
+            if (nbImage==1){
+                IJ.log("WARNING, you should provide a stack with more than 1 image");
+            }
+            ImageStack ims= imp.getStack();
+            
+            int width=imp.getWidth();
+            int height=imp.getHeight();
+            
+            boolean scmosLoad=false;
+            SCMOScamera scmoscam=null;
+            if (isSCMOS){
+                scmoscam = new SCMOScamera(path_SCMOSvariance,path_SCMOSoffset,path_SCMOSgain);
+                scmosLoad=scmoscam.loadSCMOScameraFiles(width, height, null);
+            }
+                
+            
+            
+            
+            double [][][] image =new double [nbImage][width][height];
+            
+            for (int i=0;i<nbImage;i++){
+                ImageProcessor ip;
+                if (stacked){
+                    ip = ims.getProcessor(i+1);
+                }
+                else{
+                    imp.setSlice(i+1);
+                    ip = imp.getProcessor();
+                }
+                for (int ii=0,jj=0;ii<width;ii++,jj++){
+                    for (int iii=0,jjj=0;iii<height;iii++,jjj++){
+                        image[i][ii][iii]=ip.getPixelValue(jj, jjj);
+                    }
+                }
+            }
+            
+            
+            GenericDialog gd = new GenericDialog("ZOLA: Convert to photon count");
+
+            
+            
+                
+            Font font = gd.getFont();
+            Font fontBold= gd.getFont();
+            try{
+                fontBold=new Font(font.getName(),Font.BOLD,font.getSize());
+            }catch(Exception e){}
+            
+            
+            
+            
+                if (!isSCMOS){
+                    gd.addMessage("Camera is EMCCD",fontBold);
+                    gd.addMessage("ADU = "+adu);
+                    gd.addNumericField("Gain = ", gain, 1,6,"(1 if camera provides photon counts)");
+                    gd.addMessage("Offset = "+offset);
+                    
+                }
+                else{
+                    
+                    if (scmosLoad){
+                        gd.addMessage("Camera is SCMOS",fontBold);
+                    }
+                    else{
+                        gd.addMessage("Camera is SCMOS / problem with image sizes / ERROR will occur",fontBold);
+                    }
+//                    String sep=File.separator;
+//                    String [] var=this.path_SCMOSvariance.split(sep);
+//                    String [] off=this.path_SCMOSoffset.split(sep);
+//                    String [] gain=this.path_SCMOSgain.split(sep);
+//                    gd.addMessage("Offset file = "+off[var.length-1]);
+//                    gd.addMessage("Variance file = "+var[var.length-1]);
+//                    gd.addMessage("Gain file = "+gain[var.length-1]);
+                    
+                    gd.addMessage("Offset file = "+path_SCMOSoffset);
+                    gd.addMessage("Variance file = "+path_SCMOSvariance);
+                    gd.addMessage("Gain file = "+path_SCMOSgain);
+                }
+                
+            
+            
+            gd.showDialog();
+            
+            
+            if (!gd.wasCanceled()){
+                
+                if (!isSCMOS){
+                    gain = (double)gd.getNextNumber();
+                }
+                
+            }
+            else{
+                processing=new Boolean(false);
+                return;
+            }
+            
+            
+            
+            
+            
+            if (!isSCMOS && (gain<=0||adu<=0)){
+                new WaitForUserDialog("Error message", "gain has to be positive (default value = 1)").show();
+            }
+            else{
+                
+                
+                InitBackgroundAndPhotonNumber paramImage ;
+                if (isSCMOS){
+                    
+                    double [][][] imageH =new double [nbImage][width][height];
+
+                    for (int z=0;z<image.length;z++){
+                        for (int i=0;i<image[0].length;i++){
+                            for (int ii=0;ii<image[0][0].length;ii++){
+                                
+                                imageH[z][i][ii]=((image[z][i][ii]-scmoscam.scmosoffset[i][ii])/scmoscam.scmosgain[i][ii]);
+                                image[z][i][ii]=((image[z][i][ii]-scmoscam.scmosoffset[i][ii])/scmoscam.scmosgain[i][ii])+scmoscam.scmosvargain[i][ii];
+
+                                if (image[z][i][ii]<0){
+                                    image[z][i][ii]=0;
+                                }
+                            }
+                        }
+
+                    }
+                    IJ.log("homogeneous image should have homogeneous aspect but does not follow Poisson statistic (mean=var) ; the difference is var/gain^2");
+                    ImageShow.imshow(imageH,"homogeneous image");
+                    
+                }
+                else{
+        
+        
+                    width=image[0].length;
+                    height=image[0][0].length;
+
+                    double slope=adu/gain;
+                    double intercept=-offset*adu/gain;
+                    for (int z=0;z<image.length;z++){
+                        for (int i=0;i<image[0].length;i++){
+                            for (int ii=0;ii<image[0][0].length;ii++){
+
+                                image[z][i][ii]=(image[z][i][ii]*slope+intercept)+1;
+
+                                if (image[z][i][ii]<0){
+                                    image[z][i][ii]=0;
+                                }
+
+
+
+                            }
+                        }
+
+
+                    }
+
+                
+
+                }
+                ImageShow.imshow(image,"photon count image");
+            }
+        }
+        processing=new Boolean(false);
+    }
+    
+            
+    void scmosComputeOffsetAndVariance(){
+        
+        
+        
+        
+            
+        ImagePlus imp=IJ.getImage();
+
+
+        if ((processing !=null)&&(processing.booleanValue()==true)){
+            printErrorMessage();
+            return;
+        }
+        
+        processing=new Boolean(true);
+
+
+        if (imp==null){
+            IJ.log("no image opened");
+            processing=new Boolean(false);
+        }
+        /*else if ((imp.getWidth()>180)||(imp.getHeight()>180)){
+            IJ.log("image should be < 180 pixels for the moment");
+        }*/
+        else{
+
+            GenericDialog gd=null;
+
+            //since NonBlockingGenericDialog incompatible with headless mode:
+            if (!(java.awt.GraphicsEnvironment.isHeadless() || IJ.isMacro())){
+                gd = new NonBlockingGenericDialog("ZOLA: SCMOS calibration 1/2");
+            }
+            else{
+                gd = new GenericDialog("ZOLA: SCMOS calibration 1/2");
+            }
+
+            Font font = gd.getFont();
+            Font fontBold= gd.getFont();
+            try{
+                fontBold=new Font(font.getName(),Font.BOLD,font.getSize());
+            }catch(Exception e){}
+            gd.addMessage("This function convert pixels in photon count for SCMOS cameras",fontBold);
+            gd.addMessage("it estimates variance and offset for each pixel from a stack of black images",fontBold);
+            
+            
+            
+            String pathdir=IJ.getDirectory("image");
+            if (pathdir==null){
+                pathdir=IJ.getDirectory("current");
+            }
+
+            String pathoffset=pathdir+"ZOLA_SCMOS_offset.tif";
+            String pathvariance=pathdir+"ZOLA_SCMOS_variance.tif";
+            
+            
+            gd.addStringField("Offset path to save:",pathoffset,sizeTextString); 
+            
+            gd.addStringField("Variance path to save:",pathvariance,sizeTextString); 
+            
+
+            TextField textField; 
+            TextField textField2; 
+            // Add a mouse listener to the config file field 
+            if (!(java.awt.GraphicsEnvironment.isHeadless() || IJ.isMacro())) 
+            {
+                int t=0;
+                //textField
+                Vector<TextField> texts = (Vector<TextField>) gd.getStringFields(); 
+
+
+
+                textField = texts.get(t++); 
+                MouseOptionSave mos=new MouseOptionSave(textField,pathoffset,"Path offset image");
+                textField.addMouseListener(mos); 
+                
+                textField2 = texts.get(t++); 
+                mos=new MouseOptionSave(textField2,pathvariance,"Path variance image");
+                textField2.addMouseListener(mos); 
+
+
+
+                //NumericField
+                Vector<TextField> nums = (Vector<TextField>) gd.getNumericFields(); 
+
+
+
+                
+
+            }
+
+
+            gd.showDialog();
+
+
+
+
+
+            if (!gd.wasCanceled()){
+
+                path_SCMOSoffset = gd.getNextString(); 
+                path_SCMOSvariance = gd.getNextString(); 
+            }
+            
+            int nnn=(path_SCMOSoffset.lastIndexOf("."));
+            if ((nnn>=0)&&(nnn>path_SCMOSoffset.length()-6)){
+                path_SCMOSoffset=path_SCMOSoffset.substring(0,nnn);
+            }
+            path_SCMOSoffset+=".tif";
+            
+            nnn=(path_SCMOSvariance.lastIndexOf("."));
+            if ((nnn>=0)&&(nnn>path_SCMOSvariance.length()-6)){
+                path_SCMOSvariance=path_SCMOSvariance.substring(0,nnn);
+            }
+            path_SCMOSvariance+=".tif";
+                
+            int width=imp.getWidth();
+            int height=imp.getHeight();
+
+            
+                int xstart=0;
+                int ystart=0;
+                
+                
+                
+                
+                
+                boolean stacked=true;
+                int nbImage=imp.getNSlices();
+                if (imp.getNSlices()==1){
+                    nbImage=imp.getNFrames();
+                    stacked=false;
+                }
+                ImageStack [] ims=null;
+                if (stacked){
+                    ims = new ImageStack[1];
+                    ims[0]=imp.getStack();
+
+                }
+            
+            
+            
+                
+                ImageProcessor ip;
+                
+                ArrayList<float [][]> al = new ArrayList<float [][]>();
+                
+                //double[][] offset = new double[width][height];
+                //double[][] var = new double[width][height];
+                float [][][] pix= new float[width][height][nbImage];
+                //ImageStack imsres = new ImageStack(width,height);
+                int index=0;
+                
+                //copy image stack in pix vector
+                for (int i=0;i<nbImage;i++){
+                    IJ.showProgress(0.25*(double)i/(double)(nbImage));
+                    if (stacked){
+                        ip = ims[0].getProcessor(i+1);
+                    }
+                    else{
+                        imp.setSlice(i+1);
+                        ip = imp.getProcessor();
+                    }
+                    
+                    //fill tmp stack at the beginning
+                    
+                    for (int u=0;u<width;u++){
+                        for (int uu=0;uu<height;uu++){
+                            float v=ip.getPixelValue(u+xstart, uu+ystart);
+                            pix[u][uu][i]=v;
+                        }
+                    }
+                }
+                
+                ImageProcessor ipoffset= new FloatProcessor(width,height);
+                ImageProcessor ipvariance= new FloatProcessor(width,height);
+                //compute mean 
+                float offset,var;
+                for (int i=0;i<width;i++){
+                    for (int ii=0;ii<height;ii++){
+                        offset=0;
+                        for (int j=0;j<nbImage;j++){
+                            offset+=pix[i][ii][j];
+                        }
+                        offset/=nbImage;
+                        ipoffset.putPixelValue(i, ii, offset);
+                        
+                        var=0;
+                        for (int j=0;j<nbImage;j++){
+                            var+=(pix[i][ii][j]-offset)*(pix[i][ii][j]-offset);
+                        }
+                        var/=nbImage;
+                        ipvariance.putPixelValue(i, ii, var);
+                    }
+                }  
+                
+                
+                
+                ImagePlus impoffset = new ImagePlus("SCMOS_offset",ipoffset);
+                ImagePlus impvariance = new ImagePlus("SCMOS_variance",ipvariance);
+                FileSaver fs_offset = new FileSaver(impoffset);
+                fs_offset.saveAsTiff(path_SCMOSoffset);
+                FileSaver fs_variance = new FileSaver(impvariance);
+                fs_variance.saveAsTiff(path_SCMOSvariance);
+                
+                
+                
+                //ImagePlus impres = new ImagePlus("result",imsres);
+                //impres.show();
+                IJ.log("SCMOS calibration part 1/2 finished.");
+                IJ.showProgress(1);
+                
+                isSCMOS=true;
+                
+                
+            }
+        
+        processing=new Boolean(false);
+        
+        
+    }
+    
+    
+    
+    
+    void photonCountEstimation(){
+        
+        
+        
+        
+            
+        ImagePlus imp=IJ.getImage();
+
+
+        if ((processing !=null)&&(processing.booleanValue()==true)){
+            printErrorMessage();
+            return;
+        }
+        
+        processing=new Boolean(true);
+
+
+        if (imp==null){
+            IJ.log("no image opened");
+            processing=new Boolean(false);
+        }
+        /*else if ((imp.getWidth()>180)||(imp.getHeight()>180)){
+            IJ.log("image should be < 180 pixels for the moment");
+        }*/
+        else{
+
+            GenericDialog gd=null;
+
+            //since NonBlockingGenericDialog incompatible with headless mode:
+            if (!(java.awt.GraphicsEnvironment.isHeadless() || IJ.isMacro())){
+                gd = new NonBlockingGenericDialog("ZOLA: photon count estimation");
+            }
+            else{
+                gd = new GenericDialog("ZOLA: photon count estimation");
+            }
+
+                Font font = gd.getFont();
+                Font fontBold= gd.getFont();
+                try{
+                    fontBold=new Font(font.getName(),Font.BOLD,font.getSize());
+                }catch(Exception e){}
+
+                gd.addMessage("This function estimate ADU and offset to convert pixels in photon count",fontBold);
+                gd.addMessage("Please, import a stack of images (~100 images) containing non moving objects (such as beads)",fontBold);
+                
+                gd.addNumericField("Gain: ", this.gain, 3,6,"");
+
+                
+                gd.showDialog();
 
 
 
@@ -9896,6 +10843,9 @@ void frc_fsc(){
         if (imageInput1.length==1){
             FRC frc = new FRC(imageInput1[0],imageInput2[0],sizeRendering);
             frc.run();
+            
+            FRCbiop bp = new FRCbiop();
+            bp.test(imp1.getProcessor(),imp2.getProcessor(),sizeRendering);
         }
         else{
             FSC fsc = new FSC(imageInput1,imageInput2,sizeRendering,sizeRenderingZ);
@@ -10067,78 +11017,277 @@ public int compare(int[] o1, int[] o2) {
                 boolean isstack1=true;
                 for (int i=0;i<sl.fl.size();i++){
 
-                    FrameLocalization fl = new FrameLocalization(sl.fl.get(frameList[i][0]).numFrame);
-                    for (int ii=0;ii<sl.fl.get(frameList[i][0]).loc.size();ii++){
-                        fl.loc.add(sl.fl.get(frameList[i][0]).loc.get(ii).copy());
-                        locNumber++;
-                        if ((isstack1)&&(locNumber>totLocNumber/2)){
-                            sl1.fl.add(fl);
-                            isstack1=false;
-                            fl = new FrameLocalization(sl.fl.get(frameList[i][0]).numFrame);
-                        }
+                    FrameLocalization fl = new FrameLocalization(sl.fl.get(frameList[i][0]).numFrame);
+                    for (int ii=0;ii<sl.fl.get(frameList[i][0]).loc.size();ii++){
+                        fl.loc.add(sl.fl.get(frameList[i][0]).loc.get(ii).copy());
+                        locNumber++;
+                        if ((isstack1)&&(locNumber>totLocNumber/2)){
+                            sl1.fl.add(fl);
+                            isstack1=false;
+                            fl = new FrameLocalization(sl.fl.get(frameList[i][0]).numFrame);
+                        }
+
+                    }
+                    if (isstack1){
+                        sl1.fl.add(fl);
+                    }
+                    else{
+                        sl2.fl.add(fl);
+                    }
+                }
+            }
+            else{//here, we split randomly
+                IJ.log("localizations splitted in 2 randomly");
+                for (int i=0;i<sl.fl.size();i++){
+
+                    FrameLocalization fl1 = new FrameLocalization(sl.fl.get(i).numFrame);
+                    FrameLocalization fl2 = new FrameLocalization(sl.fl.get(i).numFrame);
+                    for (int ii=0;ii<sl.fl.get(i).loc.size();ii++){
+                        if (Math.random()<0.5){
+                            fl1.loc.add(sl.fl.get(i).loc.get(ii).copy());
+                        }
+                        else{
+                            fl2.loc.add(sl.fl.get(i).loc.get(ii).copy());
+                        }
+
+                    }
+                    sl1.fl.add(fl1);
+                    sl2.fl.add(fl2);
+                    
+                }
+            }
+            if (is3Drendering){
+                ZRendering.nameHistPlot="part1";
+                ZRendering.hist3Dgaussian(sl1,sizeRendering,sizeRenderingZ,minX,maxX,minY,maxY,minZ,maxZ);
+                //ZRendering.hist3D(sl1,sizeRendering,sizeRenderingZ,minX,maxX,minY,maxY,minZ,maxZ,0);
+                ZRendering.nameHistPlot="part2";
+                ZRendering.hist3Dgaussian(sl2,sizeRendering,sizeRenderingZ,minX,maxX,minY,maxY,minZ,maxZ);
+                //ZRendering.hist3D(sl2,sizeRendering,sizeRenderingZ,minX,maxX,minY,maxY,minZ,maxZ,0);
+            }
+            else{
+                ZRendering.nameHistPlot="part1";
+                ZRendering.hist2Dgaussian(sl1,sizeRendering,minX,maxX,minY,maxY,minZ,maxZ);
+                //ZRendering.hist2D(sl1,sizeRendering,minX,maxX,minY,maxY,minZ,maxZ,0);
+                ZRendering.nameHistPlot="part2";
+                ZRendering.hist2Dgaussian(sl2,sizeRendering,minX,maxX,minY,maxY,minZ,maxZ);
+                //ZRendering.hist2D(sl2,sizeRendering,minX,maxX,minY,maxY,minZ,maxZ,0);
+            }
+            
+        }
+        else{
+            IJ.log("please import a localization table first");
+        }
+        
+    }
+    
+    
+    
+    
+    
+    
+    void test2(){
+        
+        
+        
+        
+        
+        if ((processing !=null)&&(processing.booleanValue()==true)){
+            printErrorMessage();
+            return;
+        }
+        
+        processing=new Boolean(true);
+            
+        
+        
+        double maxValuePlot=-1;
+        
+        String pathResult="";
+        
+        double photonNumber=3000;
+        double background=50;
+        
+        String [] fieldSimul = {"molecule moves","objective moves"};
+        
+        GenericDialog gd = new GenericDialog("ZOLA: test 2 PSF cross correl");
+        
+        
+        Font font = gd.getFont();
+        Font fontBold= gd.getFont();
+            try{
+                fontBold=new Font(font.getName(),Font.BOLD,font.getSize());
+            }catch(Exception e){}        
+                
+        gd.addCheckbox("Run_on_GPU", GPU_computation);
+        
+        //String path_calibration2="";
+        
+        gd.addStringField("Calibration_file:",path_calibration,sizeTextString); 
+        
+        //gd.addStringField("Calibration_file_cam2 [optional]:",path_calibration2,sizeTextString); 
+        
+        
+        
+        
+        gd.addMessage("Parameters",fontBold);
+
+        gd.addNumericField("patch_size: ", sizePatchPhaseRet,0,6,"(pixels)");
+        
+        
+        
+        
+        gd.addNumericField("Axial_range: ", axialRange,2,6,"(µm)");
+        
+        gd.addNumericField("Z_step: ", stepZloc*1000,1,6,"(nm)");
+        
+        gd.addNumericField("Mounting_medium_refractive index: ", nwat, 3,6,"1.33 for water");
+        
+        gd.addNumericField("Distance_focus_to_coverslip (µm)", zfocus, 3);
+        
+        
+        gd.addNumericField("Photon_number: ", photonNumber,0);
+        gd.addNumericField("Background_intensity: ", background,0);
+        
+        //gd.addStringField("Plot_max_value (µm) [optional]: ", "",6);
+        
+        //gd.addStringField("Result_table [optional]:",pathResult,sizeTextString); 
+        
+        TextField textField; 
+        TextField textField2; 
+        TextField textField3; 
+        // Add a mouse listener to the config file field 
+        if (!(java.awt.GraphicsEnvironment.isHeadless() || IJ.isMacro())) 
+        {
+            int t=0;
+            
+            Vector<TextField> texts = (Vector<TextField>) gd.getStringFields(); 
+            textField = texts.get(t++); 
+            MouseOptionLoad mol=new MouseOptionLoad(textField,path_calibration,"Import calibration file");
+            textField.addMouseListener(mol); 
+            
+            //textField3 = texts.get(t++); 
+            //mol=new MouseOptionLoad(textField3,path_calibration,"Import cal. file cam 2 (optional)");
+            //textField3.addMouseListener(mol); 
+//            textField2 = texts.get(t++); 
+//            textField3 = texts.get(t++); 
+//            MouseOptionSave mos=new MouseOptionSave(textField3,path_calibration,"Export CRLB result file");
+//            textField3.addMouseListener(mos); 
+            
+            
+        }
+        
+        
+        gd.showDialog();
+
+
+
+        if (!gd.wasCanceled()){
+            
+            GPU_computation = (boolean)gd.getNextBoolean();
+            
+            path_calibration = gd.getNextString(); 
+            //path_calibration2 = gd.getNextString(); 
+            
+            sizePatchPhaseRet = (int)gd.getNextNumber();
+            
+            
+            axialRange= (double)gd.getNextNumber();
+            stepZloc=((double)gd.getNextNumber())/1000.;
+            nwat=(double)gd.getNextNumber();
+            zfocus=(double)gd.getNextNumber();
+            
+            
+            photonNumber=(double)gd.getNextNumber();
+            background=(double)gd.getNextNumber();
+            
+//            try{
+//                String str_=gd.getNextString();
+//                if (str_.length()>0){
+//                    maxValuePlot=(double)Double.parseDouble(str_);
+//                }
+//                else{
+//                    maxValuePlot=-1;
+//                }
+//            }catch(Exception ee){maxValuePlot=-1;}
+//            pathResult = gd.getNextString(); 
+            
+
+        }
+        else{
+            processing=new Boolean(false);
+            return;
+        }
+        
+        
+        if (pathResult.equals(path_calibration)){
+            IJ.log("Process stopped: the result crlb path has to be different than calibration path");
+            processing=new Boolean(false);
+            return;
+        }
+        
+//        if (pathResult.equals(path_calibration2)){
+//            if (pathResult.length()>0){//only if path not null
+//                IJ.log("Process stopped: the result crlb path has to be different than calibration path");
+//                processing=new Boolean(false);
+//                return;
+//            }
+//        }
+        
+        
+        
+        
+        int sizeFFT=128;
+        sizeFFT=Math.max(sizeFFT,sizePatchPhaseRet*2);
+        sizeFFT=Math.min(sizeFFT,maxFFT);
+        
+        if (GPU_computation){
+            MyCudaStream.init(1);
+            org.pasteur.imagej.process.gpu.DataPhase_ dp = new org.pasteur.imagej.process.gpu.DataPhase_(sizeFFT,path_calibration);
+            if (!dp.loading){
+                IJ.log("impossible to load "+path_calibration);
+                processing=new Boolean(false);
+                return;
+            }
+
+            if (dp.param!=null){
 
-                    }
-                    if (isstack1){
-                        sl1.fl.add(fl);
-                    }
-                    else{
-                        sl2.fl.add(fl);
-                    }
-                }
-            }
-            else{//here, we split randomly
-                IJ.log("localizations splitted in 2 randomly");
-                for (int i=0;i<sl.fl.size();i++){
 
-                    FrameLocalization fl1 = new FrameLocalization(sl.fl.get(i).numFrame);
-                    FrameLocalization fl2 = new FrameLocalization(sl.fl.get(i).numFrame);
-                    for (int ii=0;ii<sl.fl.get(i).loc.size();ii++){
-                        if (Math.random()<0.5){
-                            fl1.loc.add(sl.fl.get(i).loc.get(ii).copy());
-                        }
-                        else{
-                            fl2.loc.add(sl.fl.get(i).loc.get(ii).copy());
-                        }
+                //dp.
+                dp.setSizeoutput(sizePatchPhaseRet);
+                dp.setNwat(nwat);
+                dp.param.Zfocus=zfocus;
+                
+                
+                org.pasteur.imagej.process.gpu.PSFCrossCorrel psfcc = new org.pasteur.imagej.process.gpu.PSFCrossCorrel(dp,axialRange, stepZloc);
+                psfcc.run();
+                
+                
+                dp.free();
 
-                    }
-                    sl1.fl.add(fl1);
-                    sl2.fl.add(fl2);
-                    
-                }
-            }
-            if (is3Drendering){
-                ZRendering.nameHistPlot="part1";
-                ZRendering.hist3Dgaussian(sl1,sizeRendering,sizeRenderingZ,minX,maxX,minY,maxY,minZ,maxZ);
-                //ZRendering.hist3D(sl1,sizeRendering,sizeRenderingZ,minX,maxX,minY,maxY,minZ,maxZ,0);
-                ZRendering.nameHistPlot="part2";
-                ZRendering.hist3Dgaussian(sl2,sizeRendering,sizeRenderingZ,minX,maxX,minY,maxY,minZ,maxZ);
-                //ZRendering.hist3D(sl2,sizeRendering,sizeRenderingZ,minX,maxX,minY,maxY,minZ,maxZ,0);
-            }
-            else{
-                ZRendering.nameHistPlot="part1";
-                ZRendering.hist2Dgaussian(sl1,sizeRendering,minX,maxX,minY,maxY,minZ,maxZ);
-                //ZRendering.hist2D(sl1,sizeRendering,minX,maxX,minY,maxY,minZ,maxZ,0);
-                ZRendering.nameHistPlot="part2";
-                ZRendering.hist2Dgaussian(sl2,sizeRendering,minX,maxX,minY,maxY,minZ,maxZ);
-                //ZRendering.hist2D(sl2,sizeRendering,minX,maxX,minY,maxY,minZ,maxZ,0);
             }
-            
+
+            MyCudaStream.destroy();
         }
         else{
-            IJ.log("please import a localization table first");
+            IJ.log("nothing");
+
         }
         
-    }
-    
-    
-    
-    
-    
-    
-    void test2(){
+//        classes.cudaProcess.DataPhase dp2=null;
+//        if (path_calibration2.length()>3){
+//            dp2 = new classes.cudaProcess.DataPhase(sizeFFT,path_calibration2);
+//        }
+        
+
+        processing=new Boolean(false);
         
         
+    }
+     
+    void test3(){
+        
         
+        //try wiener deconvolution
         
         
         if ((processing !=null)&&(processing.booleanValue()==true)){
@@ -10148,6 +11297,9 @@ void test2(){
         
         processing=new Boolean(true);
             
+        ImagePlus imp=IJ.getImage();
+        
+        
         
         
         double maxValuePlot=-1;
@@ -10300,16 +11452,54 @@ void test2(){
             }
 
             if (dp.param!=null){
-
+                
+                
+                
+                Roi r=imp.getRoi();
+                int width=imp.getWidth();
+                int height=imp.getHeight();
+                int startX=0;
+                int startY=0;
+                if (r!=null){
+                    Rectangle bound=imp.getRoi().getBounds();
+                    startX=bound.x;
+                    startY=bound.y;
+                    width=bound.width;
+                    height=bound.height;
+                }
+                
+                int wh=Math.min(width, height);
+                double [][][] imageInput = new double[imp.getNFrames()*imp.getNSlices()][wh][wh];
+                
+                for (int i=0;i<imp.getNSlices();i++){
+                    for (int ii=0;ii<imp.getNFrames();ii++){
+                        imp.setSlice(i*imp.getNFrames()+ii+1);
+                        ImageProcessor p=imp.getProcessor();
+                        for (int u=0;u<wh;u++){
+                            for (int uu=0;uu<wh;uu++){
+                                double pixVal = p.getPixelValue(u+startX, uu+startY);
+                                imageInput[i*imp.getNFrames()+ii][u][uu]=pixVal;
+                            }
+                        }
+                    }
+                }
+                
+                IJ.log("sizFFT "+sizeFFT+"  "+dp.param.sizeoutput);
+                
 
                 //dp.
-                dp.setSizeoutput(sizePatchPhaseRet);
+                //dp.setSizeoutput(128);
                 dp.setNwat(nwat);
                 dp.param.Zfocus=zfocus;
                 
                 
-                org.pasteur.imagej.process.gpu.PSFCrossCorrel psfcc = new org.pasteur.imagej.process.gpu.PSFCrossCorrel(dp,axialRange, stepZloc);
-                psfcc.run();
+                org.pasteur.imagej.process.gpu.InverseFiltering invfil = new org.pasteur.imagej.process.gpu.InverseFiltering(dp,axialRange, stepZloc);
+                
+                
+                
+                
+                
+                invfil.run(imageInput);
                 
                 
                 dp.free();
@@ -10333,8 +11523,11 @@ void test2(){
         
         
     }
-     
-    void test3(){
+    
+    
+    
+    
+    void test4(){
         
         
         //try wiener deconvolution
@@ -10349,15 +11542,13 @@ void test3(){
             
         ImagePlus imp=IJ.getImage();
         
-        
+        double minPosizionZ=-axialRange/2;
         
         
         double maxValuePlot=-1;
         
         String pathResult="";
         
-        double photonNumber=3000;
-        double background=50;
         
         String [] fieldSimul = {"molecule moves","objective moves"};
         
@@ -10386,7 +11577,7 @@ void test3(){
         gd.addNumericField("patch_size: ", sizePatchPhaseRet,0,6,"(pixels)");
         
         
-        
+        gd.addNumericField("min_Z_position: ", minPosizionZ,2,6,"(µm)");
         
         gd.addNumericField("Axial_range: ", axialRange,2,6,"(µm)");
         
@@ -10397,9 +11588,6 @@ void test3(){
         gd.addNumericField("Distance_focus_to_coverslip (µm)", zfocus, 3);
         
         
-        gd.addNumericField("Photon_number: ", photonNumber,0);
-        gd.addNumericField("Background_intensity: ", background,0);
-        
         //gd.addStringField("Plot_max_value (µm) [optional]: ", "",6);
         
         //gd.addStringField("Result_table [optional]:",pathResult,sizeTextString); 
@@ -10441,7 +11629,7 @@ void test3(){
             //path_calibration2 = gd.getNextString(); 
             
             sizePatchPhaseRet = (int)gd.getNextNumber();
-            
+            minPosizionZ= (double)gd.getNextNumber();
             
             axialRange= (double)gd.getNextNumber();
             stepZloc=((double)gd.getNextNumber())/1000.;
@@ -10449,9 +11637,6 @@ void test3(){
             zfocus=(double)gd.getNextNumber();
             
             
-            photonNumber=(double)gd.getNextNumber();
-            background=(double)gd.getNextNumber();
-            
 //            try{
 //                String str_=gd.getNextString();
 //                if (str_.length()>0){
@@ -10488,9 +11673,7 @@ void test3(){
         
         
         
-        int sizeFFT=128;
-        sizeFFT=Math.max(sizeFFT,sizePatchPhaseRet*2);
-        sizeFFT=Math.min(sizeFFT,maxFFT);
+        int sizeFFT=512;
         
         if (GPU_computation){
             MyCudaStream.init(1);
@@ -10518,38 +11701,39 @@ void test3(){
                     height=bound.height;
                 }
                 
-                int wh=Math.min(width, height);
-                double [][][] imageInput = new double[imp.getNFrames()*imp.getNSlices()][wh][wh];
+                double [][][] imageInput = new double[imp.getNFrames()*imp.getNSlices()][sizeFFT][sizeFFT];
                 
                 for (int i=0;i<imp.getNSlices();i++){
                     for (int ii=0;ii<imp.getNFrames();ii++){
                         imp.setSlice(i*imp.getNFrames()+ii+1);
                         ImageProcessor p=imp.getProcessor();
-                        for (int u=startX;u<wh;u++){
-                            for (int uu=startY;uu<wh;uu++){
-                                double pixVal = p.getPixelValue(u, uu);
+                        for (int u=0;u<Math.min(width,sizeFFT);u++){
+                            for (int uu=0;uu<Math.min(height,sizeFFT);uu++){
+                                double pixVal = p.getPixelValue(u+startX, uu+startY);
                                 imageInput[i*imp.getNFrames()+ii][u][uu]=pixVal;
                             }
                         }
                     }
                 }
                 
-                
+                ImageShow.imshow(imageInput);
                 
 
-                //dp.
-                dp.setSizeoutput(wh);
+                IJ.log("sizFFT "+sizeFFT+"  "+dp.param.sizeoutput);
+                /*if (wh<sizeFFT){
+                    dp.setSizeoutput(wh);
+                }*/
                 dp.setNwat(nwat);
                 dp.param.Zfocus=zfocus;
                 
                 
-                org.pasteur.imagej.process.gpu.InverseFiltering invfil = new org.pasteur.imagej.process.gpu.InverseFiltering(dp,axialRange, stepZloc);
+                org.pasteur.imagej.process.gpu.Filtering fil = new org.pasteur.imagej.process.gpu.Filtering(dp,minPosizionZ,axialRange, stepZloc);
                 
                 
                 
                 
                 
-                invfil.run(imageInput);
+                fil.run(imageInput);
                 
                 
                 dp.free();
@@ -10703,7 +11887,7 @@ void test(){
         
         //simulation different PSF
         //rose
-        {
+        if (false){
             int num=60;
             double x0=0;
             double y0=0;
@@ -10746,8 +11930,54 @@ void test(){
             ImageShow.imshow(image,"rose");
         }
         
+        
+        
+        //fermat
+        if (true){
+            int num=60;
+            double x0=0;
+            double y0=0;
+            double [] x1=new double[num];
+            double [] y1=new double[num];
+
+            int size=128;
+            double [][][] image = new double[num][size][size];
+
+            for (int i=0;i<num;i++){
+                double angle=0+((double)(i)/(double)(num))*Math.PI;
+                double r=16*Math.cos(angle);
+                x1[i]=r*Math.cos(angle);
+                y1[i]=r*Math.sin(angle);
+            }
+
+            double amplitude0=10000;
+            double amplitude1=10000;
+
+            double sigma=1.5;
+            double photonNumber=0;
+            for (int z=0;z<num;z++){
+                double centerX=(x0+x1[z])/2.;
+                double centerY=(y0+y1[z])/2.;
+                double cx0=x0+(double)size/2.-centerX;
+                double cy0=y0+(double)size/2.-centerY;
+                double cx1=x1[z]+(double)size/2.-centerX;
+                double cy1=y1[z]+(double)size/2.-centerY;
+                for (int i=0;i<size;i++){
+                    for (int ii=0;ii<size;ii++){
+                        double amp0=amplitude0*Math.exp(-0.5*((i-cx0)*(i-cx0)+(ii-cy0)*(ii-cy0))/(sigma*sigma));
+                        double amp1=amplitude1*Math.exp(-0.5*((i-cx1)*(i-cx1)+(ii-cy1)*(ii-cy1))/(sigma*sigma));
+                        image[z][i][ii]=amp0+amp1+1;
+                        photonNumber+=amp0+amp1;
+                    }
+                }
+            }
+
+            IJ.log("photonNumber : "+photonNumber);
+            ImageShow.imshow(image,"fermat");
+        }
+        
         //DTET
-        {
+        if (false){
             int num=100;
             double x0=0;
             double y0=0;
@@ -10813,7 +12043,7 @@ void test(){
         
         
         //SPIRAL
-        {
+        if (false){
             int num=90;
             double x0=0;
             double y0=0;
@@ -10858,7 +12088,7 @@ void test(){
         
         
         //LINE
-        {
+        if (false){
             int num=90;
             double x0=0;
             double y0=0;
@@ -10897,7 +12127,7 @@ void test(){
         }
         
         //RING
-        {
+        if (false){
             int num=60;
             double x0=0;
             double y0=0;
@@ -10939,6 +12169,48 @@ void test(){
             IJ.log("photonNumber : "+photonNumber);
             ImageShow.imshow(image,"RING");
         }
+        
+        
+        
+        //SPHERE
+        {
+            int num=32;
+
+            int size=32;
+            double [][][] image = new double[num][size][size];
+
+            
+
+            double amplitude0=10000;
+            double amplitude1=10000;
+            double centerX=size/2;
+            double centerY=size/2;
+            double centerZ=num/2;
+            double diameter=8;
+            double radius=diameter/2;
+            double radius2=radius*radius;
+            double photonNumber=0;
+            IJ.log(""+centerX+"  "+centerY+"  "+centerZ);
+            
+            for (int z=0;z<num;z++){
+                
+                
+                for (int i=0;i<size;i++){
+                    for (int ii=0;ii<size;ii++){
+                        double dist=Math.sqrt((i-centerX)*(i-centerX)+(ii-centerY)*(ii-centerY)+(z-centerZ)*(z-centerZ));
+                        
+                        dist=Math.exp(-Math.pow(Math.abs(dist-radius),1));
+                        
+                        
+                        image[z][i][ii]=dist;
+                        photonNumber+=dist;
+                    }
+                }
+            }
+
+            IJ.log("photonNumber : "+photonNumber);
+            ImageShow.imshow(image,"SPHERE");
+        }
     }
     
     
@@ -11419,7 +12691,7 @@ public void run(){
                 
                 ImagePlus impTest=new ImagePlus("test",ip2);
                 
-                if (gain<=0||adu<=0){
+                if (!isSCMOS && (gain<=0||adu<=0)){
                     new WaitForUserDialog("Error message", "gain has to be positive (default value = 1)").show();
                 }
                 else{
@@ -11567,7 +12839,7 @@ public void run(){
                 }
                 
                 
-                if (gain<=0||adu<=0){
+                if (!isSCMOS && (gain<=0||adu<=0)){
                     new WaitForUserDialog("Error message", "gain has to be positive (default value = 1)").show();
                 }
                 else{
diff --git a/src/main/java/org/pasteur/imagej/ZOLA_.java b/src/main/java/org/pasteur/imagej/ZOLA_.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/configuration/ConfigFile.java b/src/main/java/org/pasteur/imagej/configuration/ConfigFile.java
old mode 100644
new mode 100755
index a505648..84cfd4d
--- a/src/main/java/org/pasteur/imagej/configuration/ConfigFile.java
+++ b/src/main/java/org/pasteur/imagej/configuration/ConfigFile.java
@@ -32,7 +32,11 @@ public class ConfigFile {
     
     static String [] option = {
         
-        "Plugins>ZOLA>  Camera setup>, \"photon count estimation\", org.pasteur.imagej.ZOLA(\"zola_photonCountEstimation\")\n",
+        "Plugins>ZOLA>  Camera setup>, \"EMCCD calibration photon count\", org.pasteur.imagej.ZOLA(\"zola_photonCountEstimation\")\n"+
+        
+        "Plugins>ZOLA>  Camera setup>SCMOS calibration photon count, \"compute offset&variance maps (1/2)\", org.pasteur.imagej.ZOLA(\"zola_SCMOS_offset_variance_map\")\n"+
+        "Plugins>ZOLA>  Camera setup>SCMOS calibration photon count, \"compute gain map (2/2)\", org.pasteur.imagej.ZOLA(\"zola_SCMOS_gain_map\")\n"+
+        "Plugins>ZOLA>  Camera setup>, \"test photon count conversion\", org.pasteur.imagej.ZOLA(\"zola_convertImageToPhotonCount\")\n",
         
         "Plugins>ZOLA>additional tools, \"Best PSF model optimization\", org.pasteur.imagej.ZOLA(\"zola_phaseoptim\")\n",
         
@@ -91,6 +95,7 @@ public class ConfigFile {
 "Plugins>ZOLA>dev>, \"RemoveLocalizations\", org.pasteur.imagej.ZOLA(\"zola_removeLocalizations\")\n" +
 "Plugins>ZOLA>dev>, \"NPC detection\", org.pasteur.imagej.ZOLA(\"zola_NPCdetection\")\n" +
 "Plugins>ZOLA>dev>, \"CRLB_3D\", org.pasteur.imagej.ZOLA(\"zola_crlb_3D\")\n" +
+"Plugins>ZOLA>dev>, \"PSF_Generator\", org.pasteur.imagej.ZOLA(\"zola_PSF_Generator\")\n" +
 "Plugins>ZOLA>dev>, \"Dec\", org.pasteur.imagej.ZOLA(\"zola_deconvolution\") "};
     
     
diff --git a/src/main/java/org/pasteur/imagej/cuda/MyCudaStream.java b/src/main/java/org/pasteur/imagej/cuda/MyCudaStream.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/cuda/MyDefaultVecKernels.java b/src/main/java/org/pasteur/imagej/cuda/MyDefaultVecKernels.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/cuda/MyVecDouble.java b/src/main/java/org/pasteur/imagej/cuda/MyVecDouble.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/cuda/MyVecKernels.java b/src/main/java/org/pasteur/imagej/cuda/MyVecKernels.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/data/FrameLocalization.java b/src/main/java/org/pasteur/imagej/data/FrameLocalization.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/data/PLocalization.java b/src/main/java/org/pasteur/imagej/data/PLocalization.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/data/StackLocalization.java b/src/main/java/org/pasteur/imagej/data/StackLocalization.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/development/FilamentResolution.java b/src/main/java/org/pasteur/imagej/development/FilamentResolution.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/development/PointResolution.java b/src/main/java/org/pasteur/imagej/development/PointResolution.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/development/Simulation.java b/src/main/java/org/pasteur/imagej/development/Simulation.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/development/Simulation2Obj.java b/src/main/java/org/pasteur/imagej/development/Simulation2Obj.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/development/SimulationOneBead.java b/src/main/java/org/pasteur/imagej/development/SimulationOneBead.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/development/SimulationOverlappingBeads.java b/src/main/java/org/pasteur/imagej/development/SimulationOverlappingBeads.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/postprocess/ChainList.java b/src/main/java/org/pasteur/imagej/postprocess/ChainList.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/postprocess/Condensation.java b/src/main/java/org/pasteur/imagej/postprocess/Condensation.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/postprocess/Deconvolution.java b/src/main/java/org/pasteur/imagej/postprocess/Deconvolution.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/postprocess/DriftCorrectionFiducial.java b/src/main/java/org/pasteur/imagej/postprocess/DriftCorrectionFiducial.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/postprocess/DualCamRegistration.java b/src/main/java/org/pasteur/imagej/postprocess/DualCamRegistration.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/postprocess/DualColorFusion.java b/src/main/java/org/pasteur/imagej/postprocess/DualColorFusion.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/postprocess/RegistrationCrossCorrel.java b/src/main/java/org/pasteur/imagej/postprocess/RegistrationCrossCorrel.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/postprocess/SMLM_CrossCorrel2Model.java b/src/main/java/org/pasteur/imagej/postprocess/SMLM_CrossCorrel2Model.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/postprocess/Test_deconv.java b/src/main/java/org/pasteur/imagej/postprocess/Test_deconv.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/postprocess/WobbleCorrection.java b/src/main/java/org/pasteur/imagej/postprocess/WobbleCorrection.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/postprocess/ZRendering.java b/src/main/java/org/pasteur/imagej/postprocess/ZRendering.java
old mode 100644
new mode 100755
index dd015d7..03505fe
--- a/src/main/java/org/pasteur/imagej/postprocess/ZRendering.java
+++ b/src/main/java/org/pasteur/imagej/postprocess/ZRendering.java
@@ -1035,7 +1035,9 @@ public static ImagePlus colorRendering(ImagePlus imp,StackLocalization sl,double
         
         
         if (printLUT){
-            addLUT(imp,minZsmart, maxZsmart,lut);
+            if (!(java.awt.GraphicsEnvironment.isHeadless() || IJ.isMacro())){
+                addLUT(imp,minZsmart, maxZsmart,lut);
+            }
         }
         
         
@@ -1500,7 +1502,7 @@ public static ImagePlus hist2D(StackLocalization sl,double pixelsizeNM,double mi
                     xx=(int)((x-minX)/pixelsizeNM)+shift;
                     yy=(int)((y-minY)/pixelsizeNM)+shift;
                     z=sl.fl.get(i).loc.get(j).Z;
-                    if ((z>minZ)&&(z<maxZ)){
+                    if ((z>=minZ)&&(z<maxZ)){
                         //IJ.log("z "+z);
                         if ((xx-shift>=0)&&(yy-shift>=0)&&(xx+shift<width)&&(yy+shift<height)){
                             for (int a=-shift,aa=0;a<=shift;a++,aa++){
@@ -1642,13 +1644,17 @@ public static ImagePlus hist2Dgaussian(StackLocalization sl,double pixelsizeNM,d
                     xx=(int)x;
                     yy=(int)y;
                     z=sl.fl.get(i).loc.get(j).Z;
+                    
                     sx=sl.fl.get(i).loc.get(j).crlb_X;
                     sy=sl.fl.get(i).loc.get(j).crlb_Y;
+                    if (Math.random()<.001){
+                        IJ.log("z"+z+"  "+sx+"  "+sy);
+                    }
                     if (sx<=0){
-                        sx=.1;
+                        sx=pixelsizeNM;
                     }
                     if (sy<=0){
-                        sy=.1;
+                        sy=pixelsizeNM;
                     }
                     double leftGauss=cste/(sx*sy);
                     decx=(x-xx)*pixelsizeNM;
@@ -1656,7 +1662,7 @@ public static ImagePlus hist2Dgaussian(StackLocalization sl,double pixelsizeNM,d
                     
                     shiftA=(int)Math.ceil(2*sx/pixelsizeNM);
                     shiftAA=(int)Math.ceil(2*sy/pixelsizeNM);
-                    if ((z>minZ)&&(z<maxZ)){
+                    if ((z>=minZ)&&(z<maxZ)){
                         //IJ.log("z "+z);
                         for (int a=-shiftA;a<=shiftA;a++){
                             for (int aa=-shiftAA;aa<=shiftAA;aa++){
@@ -2001,19 +2007,21 @@ public static ImagePlus hist3Dgaussian(StackLocalization sl,double pixelsizeXYNM
                     sz=sl.fl.get(i).loc.get(j).crlb_Z;
                     
                     if (sx<=0){
-                        sx=.1;
+                        sx=pixelsizeXYNM;
                     }
                     if (sy<=0){
-                        sy=.1;
+                        sy=pixelsizeXYNM;
                     }
                     if (sz<=0){
-                        sz=.1;
+                        sz=pixelsizeZNM;
                     }
                     double leftGauss=cste/(sx*sy*sz);
                     decx=(x-xx)*pixelsizeXYNM;
                     decy=(y-yy)*pixelsizeXYNM;
                     decz=(z-zz)*pixelsizeZNM;
-                    //IJ.log("z "+z);
+                    //if (Math.random()<.01){
+                        //IJ.log("z "+z+"  "+x+"  "+y+"  "+sx+"  "+sy+"  "+sz);
+                    //}
                     shiftA=(int)Math.ceil(2*sx/pixelsizeXYNM);
                     shiftAA=(int)Math.ceil(2*sy/pixelsizeXYNM);
                     shiftAAA=(int)Math.ceil(2*sz/pixelsizeZNM);
diff --git a/src/main/java/org/pasteur/imagej/process/InitBackgroundAndPhotonNumber.java b/src/main/java/org/pasteur/imagej/process/InitBackgroundAndPhotonNumber.java
old mode 100644
new mode 100755
index 4475c87..4220aa9
--- a/src/main/java/org/pasteur/imagej/process/InitBackgroundAndPhotonNumber.java
+++ b/src/main/java/org/pasteur/imagej/process/InitBackgroundAndPhotonNumber.java
@@ -86,7 +86,7 @@ public InitBackgroundAndPhotonNumber(double [][][] image,float [] xpoints, float
             for (int i=0;i<width;i++){
                 for (int ii=0;ii<height;ii++){
                     
-                    image[z][i][ii]=((image[z][i][ii]-scmoscam.scmosoffset[i][ii])/scmoscam.scmosgain[i][ii]);
+                    image[z][i][ii]=((image[z][i][ii]-scmoscam.scmosoffset[i][ii])/scmoscam.scmosgain[i][ii]);// here we avoid var/gain2 term that simplify computation
                     
                     if (image[z][i][ii]<0){
                         image[z][i][ii]=0;
@@ -677,8 +677,9 @@ public void computeBackground(double [][][] image){
             }
         }
         
-        
-        mean/=countMean0;
+        if (countMean0>0){
+            mean/=countMean0;
+        }
         for (i=0;i<image.length;i++){
             for (ii=0;ii<image[0].length;ii++){
                 for (iii=0;iii<image[0][0].length;iii++){
@@ -688,8 +689,11 @@ public void computeBackground(double [][][] image){
                 }
             }
         }
-        std/=countMean0;
-        std=Math.sqrt(std);
+        if (countMean0>0){
+            std/=countMean0;
+            std=Math.sqrt(std);
+        }
+        
         
         
         for (i=0;i<image.length;i++){
@@ -751,7 +755,9 @@ public void computeBackground(double [][][] image){
                     count++;
                 }
             }
-            filtMean[i]/=count;
+            if (count!=0){
+                filtMean[i]/=count;
+            }
             if (filtMean[i]<filtMean[miniposition]){
                 miniposition=i;
             }
@@ -844,8 +850,9 @@ else if (i>miniposition+(numberImageNeighbor)){
                         }
                     }
                     //Arrays.sort(vect);
-                    
-                    meantmp/=count;
+                    if (count!=0){
+                        meantmp/=count;
+                    }
                     double [] data=new double [3];
                     data[0]=ii;
                     data[1]=iii;
diff --git a/src/main/java/org/pasteur/imagej/process/PhaseParameters.java b/src/main/java/org/pasteur/imagej/process/PhaseParameters.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/process/cpu/CRLB.java b/src/main/java/org/pasteur/imagej/process/cpu/CRLB.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/process/cpu/CRLBdualCam.java b/src/main/java/org/pasteur/imagej/process/cpu/CRLBdualCam.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/process/cpu/DataPhase.java b/src/main/java/org/pasteur/imagej/process/cpu/DataPhase.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/process/cpu/DriftCorrection.java b/src/main/java/org/pasteur/imagej/process/cpu/DriftCorrection.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/process/cpu/GenericPhaseRetrieval.java b/src/main/java/org/pasteur/imagej/process/cpu/GenericPhaseRetrieval.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/process/cpu/Localization.java b/src/main/java/org/pasteur/imagej/process/cpu/Localization.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/process/cpu/LocalizationPipeline.java b/src/main/java/org/pasteur/imagej/process/cpu/LocalizationPipeline.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/process/cpu/PSFPhase.java b/src/main/java/org/pasteur/imagej/process/cpu/PSFPhase.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/process/cpu/PSFPhaseMany.java b/src/main/java/org/pasteur/imagej/process/cpu/PSFPhaseMany.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/process/cpu/Phase.java b/src/main/java/org/pasteur/imagej/process/cpu/Phase.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/process/cpu/Predetection.java b/src/main/java/org/pasteur/imagej/process/cpu/Predetection.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/process/cpu/SearchPSFcenter.java b/src/main/java/org/pasteur/imagej/process/cpu/SearchPSFcenter.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/process/cpu/ZernikePhase.java b/src/main/java/org/pasteur/imagej/process/cpu/ZernikePhase.java
old mode 100644
new mode 100755
index 9a0d13d..a77e8fd
--- a/src/main/java/org/pasteur/imagej/process/cpu/ZernikePhase.java
+++ b/src/main/java/org/pasteur/imagej/process/cpu/ZernikePhase.java
@@ -53,7 +53,7 @@ public ZernikePhase(PhaseParameters param,int zernikePolyNumber,int method){
         this.multiTrainingWithParabola=false;
         this.param=param;
         Zernike z = new Zernike(param.size_cpu,param.sizeRadiusRingPixel_cpu,zernikePolyNumber);
-        //ImageShow.imshow(z.Z,"Z");
+        
         
         numberOfCoef=zernikePolyNumber;
         
@@ -157,7 +157,7 @@ public ZernikePhase(PhaseParameters param,int [] coef){
         
         this.param=param;
         Zernike z = new Zernike(param.size_cpu,param.sizeRadiusRingPixel_cpu,maxiCoef+1);
-        
+        ImageShow.imshow(z.Z,"Z");
         numberOfCoef=coef.length;
         
             
diff --git a/src/main/java/org/pasteur/imagej/process/cpu/ZernikePhaseRetrieval.java b/src/main/java/org/pasteur/imagej/process/cpu/ZernikePhaseRetrieval.java
old mode 100644
new mode 100755
index 2936706..12af6da
--- a/src/main/java/org/pasteur/imagej/process/cpu/ZernikePhaseRetrieval.java
+++ b/src/main/java/org/pasteur/imagej/process/cpu/ZernikePhaseRetrieval.java
@@ -173,7 +173,6 @@ public ZernikePhaseRetrieval(int sizeFFT,double xystep,double zstep,double wavel
             //init background fixed value
             int [] ind=pbg[z].getIndexOfOrder(0);
             pbg[z].a[ind[0]]=paramImage.B[z];
-            
             computeBackground2D(z);
             
             
@@ -402,7 +401,11 @@ void showPhase(){
         for (int ii=0;ii<image[image_id].length;ii++){
             for (int iii=0;iii<image[image_id][ii].length;iii++){
                 for (int iiii=0;iiii<image[image_id][ii][iii].length;iiii++){
-                    som+=(image[image_id][ii][iii][iiii]-mod[ii][iii][iiii])*(image[image_id][ii][iii][iiii]-mod[ii][iii][iiii])/mod[ii][iii][iiii];
+                    double denom=Math.abs(mod[ii][iii][iiii]);
+                    if (denom==0){
+                        denom=0.00001;
+                    }
+                    som+=(image[image_id][ii][iii][iiii]-mod[ii][iii][iiii])*(image[image_id][ii][iii][iiii]-mod[ii][iii][iiii])/denom;
                     nbData++;
                 }
             }
@@ -1515,17 +1518,18 @@ private double init(double likelihood,int stackNumber,int startcoefNumber,int en
     void computeBackground2D(int idZ){
         
         for (int i=0;i<background[idZ].length;i++){
-            for (int ii=0;ii<background[idZ][0].length;ii++){
-
+            for (int ii=0;ii<background[idZ][0].length;ii++){   
+                    
                     double [] v = new double [2];
                     v[0]=i;
                     v[1]=ii;
                     background[idZ][i][ii]=pbg[idZ].transform(v);
-
                     
             }
         }
         
+        //ImageShow.imshow(background);
+        
     }
     
     
@@ -1617,6 +1621,7 @@ public double transform(double [] X){//transform one position of N dimensions
     //                        prod*=X[ii];
     //                    }
     //                    else{
+                            
                             prod*=Math.pow(X[ii], indiceLigne[ii]);
     //                    }
                         //IJ.log("* "+X[ii]+" ^ "+indiceLigne[ii]);
@@ -1857,6 +1862,9 @@ void computeLikelihood(){
                     else{
                         mod[i][ii]=psf[i][ii]*(paramImage.A[idStack][idSlice]+fit_a[idStack]+fit_a_each[idStack][idSlice])+background[idStack][i][ii];
                     }
+                    if (mod[i][ii]<=0){
+                        mod[i][ii]=0.00000001;
+                    }
                     likelihood+=mod[i][ii]-image[idStack][idSlice][i][ii]*Math.log(mod[i][ii]);
                     
                     
diff --git a/src/main/java/org/pasteur/imagej/process/gpu/CRLB_.java b/src/main/java/org/pasteur/imagej/process/gpu/CRLB_.java
old mode 100644
new mode 100755
index 3ea5e42..8a2c001
--- a/src/main/java/org/pasteur/imagej/process/gpu/CRLB_.java
+++ b/src/main/java/org/pasteur/imagej/process/gpu/CRLB_.java
@@ -85,7 +85,7 @@ public void run(String path){
         SearchPSFcenter_ spsfc= new SearchPSFcenter_(dp,axialRange);
         double position=spsfc.getPosition();
         
-        //IJ.log("position "+position);
+        IJ.log("position "+position);
         
         boolean computeSTD=false;
         
@@ -466,6 +466,9 @@ public void run(String path){
     
     
     
+    
+    
+    
     //perform 3D localization precision (in case of 3D localization: stack of images)
     public void run3D(String path){
         
diff --git a/src/main/java/org/pasteur/imagej/process/gpu/Cross.java b/src/main/java/org/pasteur/imagej/process/gpu/Cross.java
old mode 100644
new mode 100755
index 6a81033..d3a9c61
--- a/src/main/java/org/pasteur/imagej/process/gpu/Cross.java
+++ b/src/main/java/org/pasteur/imagej/process/gpu/Cross.java
@@ -542,6 +542,7 @@ public void runCrossCorrelationOnGPU(){
         
         cudaResult =cudaMemcpyAsync(this.host_resultCorrelation, device_correlImage, totalSizeUnpadded*range.length*Sizeof.FLOAT, cudaMemcpyDeviceToHost,MyCudaStream.getCudaStream_t(streamId));
         
+        //IJ.log("siz "+range.length+"  "+width+"  "+height);
         
         for (int z=0;z<range.length;z++){
             for (int i=0;i<this.width;i++){
diff --git a/src/main/java/org/pasteur/imagej/process/gpu/DataPhase_.java b/src/main/java/org/pasteur/imagej/process/gpu/DataPhase_.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/process/gpu/DriftCorrection_.java b/src/main/java/org/pasteur/imagej/process/gpu/DriftCorrection_.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/process/gpu/DualCamLocalizationPipeline_.java b/src/main/java/org/pasteur/imagej/process/gpu/DualCamLocalizationPipeline_.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/process/gpu/DualCamLocalizationPipeline_save.java b/src/main/java/org/pasteur/imagej/process/gpu/DualCamLocalizationPipeline_save.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/process/gpu/Filtering.java b/src/main/java/org/pasteur/imagej/process/gpu/Filtering.java
new file mode 100755
index 0000000..b2fe5d8
--- /dev/null
+++ b/src/main/java/org/pasteur/imagej/process/gpu/Filtering.java
@@ -0,0 +1,190 @@
+/*
+ * To change this license header, choose License Headers in Project Properties.
+ * To change this template file, choose Tools | Templates
+ * and open the template in the editor.
+ */
+package org.pasteur.imagej.process.gpu;
+
+import ij.IJ;
+import org.pasteur.imagej.utils.*;
+import java.util.Random;
+
+/**
+ *
+ * @author benoit
+ */
+public class Filtering {
+    
+    
+    
+    
+    double axialRange; 
+    double stepZ;
+    double minPosizionZ;
+    DataPhase_ dp;
+    
+    
+    public double [] range;
+    
+    
+    public Filtering(DataPhase_ dp,double minPosizionZ,double axialRange, double stepZ){
+        this.dp=dp;
+        this.axialRange=axialRange;
+        this.stepZ=stepZ;
+        this.minPosizionZ=minPosizionZ;
+        
+    }
+    
+    
+    double [][] pad(double [][] im){
+        Random rand = new Random();
+        double mean=mean(im);
+        double var=var(im);
+        double epsilon=0;
+        
+        double [][] r = new double [im.length*2][im[0].length*2];
+        if((mean>0)&&(var>0)){
+            for (int i=0;i<im.length*2;i++){
+                for (int ii=0;ii<im[0].length*2;ii++){
+                    do{
+                        //r[i][ii]=rand.nextGaussian()*Math.sqrt(var)+mean;
+                        r[i][ii]=mean+rand.nextGaussian()*.00001;
+                    }while(r[i][ii]<0);
+                }
+            }
+        }
+        else{
+            IJ.log("error in padding");
+        }
+        for (int i=0;i<im.length;i++){
+            for (int ii=0;ii<im[0].length;ii++){
+                r[i+im.length/2][ii+im[i].length/2]=im[i][ii];
+            }
+        }
+        
+        return r;
+    }
+    
+    
+    
+    
+    public double mean(double [][] x){
+        double mean=0;
+        
+        double num=(x.length*x[0].length);
+        for (int i=0;i<x.length;i++){
+            for (int ii=0;ii<x[0].length;ii++){
+                mean+=x[i][ii];
+            }
+        }
+        mean/=num;
+        
+        return mean;
+    }
+    
+    
+    public double var(double [][] x){
+        double mean=mean(x);
+        double var=0;
+        double num=(x.length*x[0].length);
+        for (int i=0;i<x.length;i++){
+            for (int ii=0;ii<x[0].length;ii++){
+                var+=(x[i][ii]-mean)*(x[i][ii]-mean);
+            }
+        }
+        var/=num;
+        return var;
+    }
+    
+    
+    public double [][] addnoise(double [][] x,double photon,double bckg,double sigma){
+        double [][] y = new double[x.length][x[0].length];
+        Random r = new Random();
+        for (int i=0;i<x.length;i++){
+            for (int ii=0;ii<x[0].length;ii++){
+                y[i][ii]=x[i][ii]*photon+r.nextGaussian()*sigma+bckg;
+            }
+        }
+        
+        return y;
+    }
+    
+    
+    
+    public void run(double [][][] image){
+        
+        int sizeFilt=dp.param.sizeoutput;
+        
+        dp.psf.resetKz();
+        
+        //SearchPSFcenter_ spsfc= new SearchPSFcenter_(dp,axialRange);
+        //double position=spsfc.getPosition();
+        
+        double minZ=minPosizionZ;
+        double maxZ=minPosizionZ+axialRange;
+        
+        int w,h;
+        IJ.log(""+minZ+"  "+maxZ+"  "+stepZ);
+        Cross predGPU=new Cross(dp.param.sizeoutput,dp,minZ, maxZ, stepZ,.3);
+        
+        double []range=predGPU.getRange();
+        //double [][] padded=pad(psf[0]);
+        
+        
+        
+        w=image[0].length;
+        h=image[0][0].length;
+        float [][][] ResultFilt = new float[range.length][sizeFilt][sizeFilt];
+        
+        float [][][] ResultFinal = new float[range.length][w][h];
+        
+        
+        double [][] tmp = new double[sizeFilt][sizeFilt];
+        
+        for (int i=0;i<Math.ceil((double)w/(double)sizeFilt);i+=1){
+            for (int ii=0;ii<Math.ceil((double)h/(double)sizeFilt);ii+=1){
+                for (int u=0;u<sizeFilt;u++){
+                    for (int uu=0;uu<sizeFilt;uu++){
+                        if ((i*sizeFilt+u)<w && (ii*sizeFilt+uu)<h){
+                            tmp[u][uu]=(double)image[0][(int)(i*sizeFilt)+u][(int)(ii*sizeFilt)+uu];
+                        }
+                    }
+                }
+                predGPU.setImage(tmp);
+                predGPU.convolveNormalizedInFourierDomain(ResultFilt);
+                for (int y=0;y<range.length;y++){
+                    for (int u=0;u<sizeFilt;u++){
+                        for (int uu=0;uu<sizeFilt;uu++){
+                            if ((i*sizeFilt+u)<w && (ii*sizeFilt+uu)<h){
+                                ResultFinal[y][(int)(i*sizeFilt)+u][(int)(ii*sizeFilt)+uu]=ResultFilt[y][u][uu];
+                            }
+                        }
+                    }
+                }
+                ImageShow.imshow(ResultFilt,"rf "+i+"  "+ii);
+                
+            }
+            
+        }
+        
+        
+        double [][][] psf=predGPU.getPSFNonNormalized3D();
+        ImageShow.imshow(psf,"psf");
+        
+        
+        ImageShow.imshow(ResultFinal,"result filtering");
+        
+        
+        predGPU.free();
+        
+        
+        
+    }
+    
+    
+    
+    
+    
+    
+    
+}
diff --git a/src/main/java/org/pasteur/imagej/process/gpu/GS_.java b/src/main/java/org/pasteur/imagej/process/gpu/GS_.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/process/gpu/GaussianKernel_.java b/src/main/java/org/pasteur/imagej/process/gpu/GaussianKernel_.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/process/gpu/GenericPhaseOptimization_.java b/src/main/java/org/pasteur/imagej/process/gpu/GenericPhaseOptimization_.java
old mode 100644
new mode 100755
index 983f35e..61f511a
--- a/src/main/java/org/pasteur/imagej/process/gpu/GenericPhaseOptimization_.java
+++ b/src/main/java/org/pasteur/imagej/process/gpu/GenericPhaseOptimization_.java
@@ -187,73 +187,94 @@ private void phase_optim(int iterations,boolean init){
             
             
         double likelihood=Double.POSITIVE_INFINITY;
-        for (int iter=0;iter<iterations;iter++){
-            
-            IJ.showProgress(((double)iter/(double)iterations));
-            
-            
-            
-            int splitNumber=angleSplit*distSplit;
-            
-            dparam.phaseNonZer.createSplits(angleSplit, distSplit);
-            
-            double [][] tmpPhase=new double [splitNumber][];
-            for (int i=0;i<splitNumber;i++){
-                tmpPhase[i]=dparam.phaseNonZer.getValuesPhase(i);
-            }
-            
-            for (int i=0;i<splitNumber;i++){
-                double random = Math.random()*2.*Math.PI;
-                for (int ii=0;ii<tmpPhase[i].length;ii++){
-                    tmpPhase[i][ii]=random;
-                }
-            }
+        
+
+
+
+        int splitNumber=angleSplit*distSplit;
 
-            for (int i=0;i<splitNumber;i++){
-                dparam.phaseNonZer.setValuesPhase(i, tmpPhase[i]);
+        dparam.phaseNonZer.createSplits(angleSplit, distSplit);
+
+        double [][] tmpPhase=new double [splitNumber][];
+        for (int i=0;i<splitNumber;i++){
+            tmpPhase[i]=dparam.phaseNonZer.getValuesPhase(i);
+        }
+
+        for (int i=0;i<splitNumber;i++){
+            double random = Math.random()*2.*Math.PI;
+            for (int ii=0;ii<tmpPhase[i].length;ii++){
+                //tmpPhase[i][ii]=random;
+                tmpPhase[i][ii]=0;
             }
+        }
 
-            dparam.psf_fMany.updatePhase(dparam.phaseNonZer.getPointerPhase());
+        for (int i=0;i<splitNumber;i++){
+            dparam.phaseNonZer.setValuesPhase(i, tmpPhase[i]);
+        }
 
-            lik=compute_sum_CRLB((float)hdec);
 
-            //bestPhasetmp=dparam.phaseNonZer.getValuesPhase(); 
-            
+        double [] ph = dparam.phaseNonZer.getValuesPhase();
+        this.wrapPhase(ph);
+        dparam.phaseNonZer.setValuesPhase(ph);
+        dparam.psf_fMany.updatePhase(dparam.phaseNonZer.getPointerPhase());
 
+        lik=compute_sum_CRLB((float)hdec);
 
-            //dparam.phaseNonZer.setValuesPhase(bestPhasetmp);
-            //dparam.psf_fMany.updatePhase(dparam.phaseNonZer.getPointerPhase());
-            IJ.showProgress(1);
+        //bestPhasetmp=dparam.phaseNonZer.getValuesPhase(); 
 
 
-            //pattern based
-            for (int t=0;t<50;t++){
-                angleSplit=12+(int)(Math.random()*3);
-                distSplit=5+(int)(Math.random()*3);
-                //dparam.phaseNonZer.createSplits(angleSplit, distSplit);
-                lik=this.updatePhaseGeneric(.05,angleSplit,distSplit);
+
+        //dparam.phaseNonZer.setValuesPhase(bestPhasetmp);
+        //dparam.psf_fMany.updatePhase(dparam.phaseNonZer.getPointerPhase());
+        IJ.showProgress(1);
+        
+        initPhaseGeneric(0.05);
+
+        //pattern based
+        for (int t=0;t<500;t++){
+            angleSplit=16+(int)(Math.random()*5);
+            distSplit=10+(int)(Math.random()*5);
+            dparam.phaseNonZer.createSplits(angleSplit, distSplit);
+            lik=this.updatePhaseGeneric(.05,angleSplit,distSplit);
+            if (t%10==0){
                 
-                /*ImageShow.imshow(dparam.psf_fMany.getPhase(),"phase");
+            
+                ph = dparam.phaseNonZer.getValuesPhase();
+                this.wrapPhase(ph);
+                dparam.phaseNonZer.setValuesPhase(ph);
+                dparam.psf_fMany.updatePhase(dparam.phaseNonZer.getPointerPhase());
+
+                smoothPhase();
+                        
+                ImageShow.imshow(dparam.psf_fMany.getPhase(),"phaseFinal");
                 float [][][] im = new float [deltaZ.length][][];
                 for (int i=0;i<deltaZ.length;i++){
                     im[i]=dparam.psf_fMany.getPSF(i*7);
                 }
-                ImageShow.imshow(im,"PSF"+iter);*/
-
+                ImageShow.imshow(im,"PSF_tmp"+t);
+                dparam.save(path_result+"_"+t+"_.json");
             }
+        }
+        
+        
+            
+        for (int iter=0;iter<iterations;iter++){
             //pixel based
-            for (int t=0;t<5;t++){
 
 
-                IJ.showProgress(((double)t/(double)iterations));
+            IJ.showProgress(((double)iter/(double)iterations));
 
 
-                lik=this.updatePhaseGeneric(hdec);
+            lik=this.updatePhaseGeneric(hdec);
 
 
-            }
             
             
+            ph = dparam.phaseNonZer.getValuesPhase();
+            this.wrapPhase(ph);
+            dparam.phaseNonZer.setValuesPhase(ph);
+            dparam.psf_fMany.updatePhase(dparam.phaseNonZer.getPointerPhase());
+            
             ImageShow.imshow(dparam.psf_fMany.getPhase(),"phaseFinal");
             float [][][] im = new float [deltaZ.length][][];
             for (int i=0;i<deltaZ.length;i++){
@@ -261,6 +282,7 @@ private void phase_optim(int iterations,boolean init){
             }
             ImageShow.imshow(im,"PSF_tmp"+iter);
                 
+            dparam.save(path_result+"_optim_"+iter+"_.json");
                 
             if (lik<likelihood){
                 likelihood=lik;
@@ -332,6 +354,53 @@ private void phase_optim(int iterations,boolean init){
         
     }
     
+    private void smoothPhase(){
+        double [] phase= dparam.phaseNonZer.getValuesPhase();
+        double [] phaseMean= dparam.phaseNonZer.getValuesPhase();
+        int size=dparam.param.size;
+        int [][] ph=new int[size][size];
+        for (int i=0;i<size;i++){
+            for (int ii=0;ii<size;ii++){
+                ph[i][ii]=-1;
+            }
+        }
+        for (int p=0;p<dparam.param.sizeDisk;p++){
+            ph[dparam.param.disk2D[p][0]][dparam.param.disk2D[p][1]]=p;
+        }
+        for (int i=0;i<size;i++){
+            for (int ii=0;ii<size;ii++){
+                if (ph[i][ii]>=0){
+                    phaseMean[ph[i][ii]]=0;
+                    double count=0;
+                    for (int j=-1;j<=1;j++){
+                        for (int jj=-1;jj<=1;jj++){
+                            if (((i+j)>=0)&&((i+j)<size)&&((ii+jj)>=0)&&((ii+jj)<size)){
+                                if (ph[i+j][ii+jj]>=0){
+                                    double shift=0;
+                                    double minDist1=Math.abs(phase[ph[i][ii]]-phase[ph[i+j][ii+jj]]);
+                                    double minDist2=Math.abs(phase[ph[i][ii]]-phase[ph[i+j][ii+jj]]-Math.PI*2);
+                                    double minDist3=Math.abs(phase[ph[i][ii]]-phase[ph[i+j][ii+jj]]+Math.PI*2);
+                                    if (minDist2<minDist1){
+                                        shift=Math.PI*2;
+                                    }
+                                    if (minDist3<minDist1){
+                                        shift=-Math.PI*2;
+                                    }
+                                    
+                                    phaseMean[ph[i][ii]]+=phase[ph[i+j][ii+jj]]+shift;
+                                    count++;
+                                }
+                            }
+                        }
+                    }
+                    phaseMean[ph[i][ii]]/=count;
+                }
+            }
+        }
+        dparam.phaseNonZer.setValuesPhase(phaseMean);
+    }
+    
+    
     private void phase_optim_old(int iterations,boolean init){
         
         
@@ -812,7 +881,25 @@ private double updatePhaseGeneric(double h){
     
     
     
-    
+    void wrapPhase(double [] phase){
+        for (int i=0;i<phase.length;i++){
+            if (phase[i]>=0){
+                int v = (int)(phase[i]/(2*Math.PI));
+                phase[i]-=v*2*Math.PI;
+                if (phase[i]>Math.PI){
+                    phase[i]-=2*Math.PI;
+                }
+            }
+            else{
+                int v = (int)(Math.abs(phase[i])/(2*Math.PI));
+                phase[i]+=v*2*Math.PI;
+                if (phase[i]<-Math.PI){
+                    phase[i]+=2*Math.PI;
+                }
+            }
+        }
+        
+    }
     
     
     
@@ -916,6 +1003,7 @@ private double updatePhaseGeneric(double h,int angleSplit,int distSplit){
                         }
                     }
                     if (!isNan){
+                        wrapPhase(tmp);
                         dparam.phaseNonZer.setValuesPhase(p, tmp);
                         dparam.psf_fMany.updatePhase(dparam.phaseNonZer.getPointerPhase());
 
@@ -936,6 +1024,7 @@ private double updatePhaseGeneric(double h,int angleSplit,int distSplit){
                     
                 }
                 if (!found){
+                    wrapPhase(save);
                     dparam.phaseNonZer.setValuesPhase(p, save);
                     dparam.psf_fMany.updatePhase(dparam.phaseNonZer.getPointerPhase());
                     
diff --git a/src/main/java/org/pasteur/imagej/process/gpu/GenericPhaseOptimization_pi.java b/src/main/java/org/pasteur/imagej/process/gpu/GenericPhaseOptimization_pi.java
new file mode 100644
index 0000000..05f61fe
--- /dev/null
+++ b/src/main/java/org/pasteur/imagej/process/gpu/GenericPhaseOptimization_pi.java
@@ -0,0 +1,1835 @@
+/*
+ * To change this license header, choose License Headers in Project Properties.
+ * To change this template file, choose Tools | Templates
+ * and open the template in the editor.
+ */
+package org.pasteur.imagej.process.gpu;
+
+
+
+import org.pasteur.imagej.utils.*;
+import jcuda.runtime.JCuda;
+
+import jcuda.runtime.cudaError;
+import ij.IJ;
+import ij.gui.Plot;
+import java.awt.Color;
+import jcuda.Pointer;
+
+
+/**
+ *
+ * @author benoit
+ */
+public class GenericPhaseOptimization_pi {
+    
+    //int stream=0;//not parallel process: iterative
+    double [] deltaZ;
+    
+    int methodLikelihood=0;
+    
+    DataPhase_ dparam;
+    
+    
+    
+    double axialRange; 
+    double stepZ;
+    
+    double [] zCRLB;
+    
+    
+    String path_result;
+    
+    float photonNumber; 
+    float background;
+    double nwat;
+    double Zfocus;
+    public GenericPhaseOptimization_pi(int sizeFFT,int sizePatch,double xystep,double axialRange, double stepZ,double photonNumber, double background,double wavelength,double noil,double nwat,double zfocus,double na,String path_result,boolean withApoFactor){
+        
+        this.photonNumber=(float)photonNumber;
+        this.background=(float)background;
+        
+        this.axialRange=axialRange;
+        this.stepZ=stepZ;
+        
+        
+        dparam = new DataPhase_(sizeFFT,sizePatch,0,xystep,stepZ,wavelength,noil,na,1,withApoFactor);
+        dparam.param.zernikedPSF=false;
+        this.nwat=nwat;
+        this.Zfocus=zfocus;
+        
+        int nb=(int)((axialRange)/stepZ)+1;
+        
+        
+        
+        
+        deltaZ=new double [nb];
+        //IJ.log("center is : "+center+"  !!!!!!!!!!§§§§§§§§§§§§§§§!!!!!!!");
+        for (int s=0;s<nb;s++){
+            deltaZ[s]=0-(axialRange/2.)+dparam.param.zstep*(double)s;
+            //IJ.log("deltaZ "+deltaZ[s]);
+        }
+        
+        
+        
+        dparam.setMany(nb*7);//images of 1 stack are computed in parallel * 7 for fisher matrix
+        
+        
+        
+        
+        
+        this.path_result=path_result;
+        
+        
+    }
+    
+    
+    
+    
+    
+    public void run(int nbIter){
+        
+        
+        
+        dparam.setNwat(dparam.param.noil);
+        dparam.param.Zfocus=0;
+        //IJ.log("zstep "+dparam.param.zstep);
+        
+        this.phase_optim(nbIter,true);
+        
+        
+        dparam.setNwat(nwat);
+        dparam.param.Zfocus=this.Zfocus;
+                
+        
+        if (dparam.param.nwat!=dparam.param.noil){
+            IJ.log("process ran a second time to deal with refractive index mismatch");
+            SearchPSFcenter_ spsfc= new SearchPSFcenter_(dparam,axialRange);
+            double position=spsfc.getPosition();
+            for (int s=0;s<deltaZ.length;s++){
+                deltaZ[s]=position-(axialRange/2.)+dparam.param.zstep*(double)s;
+                //IJ.log("deltaZ "+deltaZ[s]);
+            }
+
+            this.phase_optim(nbIter/2,false);
+        }
+        
+        
+        float [][][] im = new float [deltaZ.length][][];
+        
+        for (int i=0;i<deltaZ.length;i++){
+            im[i]=dparam.psf_fMany.getPSF(i*7);
+        }
+        
+        //ImageShow.imshow(im,"PSF");
+        
+        
+        /*double [][] rescrlb=compute_CRLB((float)0.05);
+        for (int i=0;i<rescrlb.length;i++){
+            String s="";
+            for (int ii=0;ii<3;ii++){
+                s+=rescrlb[i][ii]+",";
+            }
+            IJ.write(""+s);
+        }*/
+        
+        if (path_result.length()>2){
+            if (path_result.endsWith("json")||path_result.endsWith("csv")){
+                dparam.save(path_result);
+            }
+            else{
+                dparam.save(path_result+".json");
+            }
+        }
+        else{
+            IJ.log("calibration file not saved...you should select a path to save it");
+        }
+        
+    }
+    
+    private void phase_optim(int iterations,boolean init){
+        
+        
+        
+        double hdec=.05;
+        
+        
+        
+        
+        double lik=Double.POSITIVE_INFINITY;
+        
+        
+        
+        IJ.showProgress(.02);
+        
+        
+        
+
+
+        
+        IJ.log("optimization started");
+        
+        
+        
+        
+        
+        IJ.log("optimization started");
+        double [] bestPhase=new double [dparam.phaseNonZer.nbDataPerImage];
+        double [] bestPhasetmp=new double [dparam.phaseNonZer.nbDataPerImage];
+        
+        
+        init_pi(hdec,6,4);
+        
+            
+            
+        double likelihood=Double.POSITIVE_INFINITY;
+        
+            
+            
+            
+
+        lik=compute_sum_CRLB((float)hdec);
+
+
+
+        ImageShow.imshow(dparam.psf_fMany.getPhase(),"phaseFinal");
+        float [][][] im = new float [deltaZ.length][][];
+        for (int i=0;i<deltaZ.length;i++){
+            im[i]=dparam.psf_fMany.getPSF(i*7);
+        }
+        ImageShow.imshow(im,"PSF_tmp");
+
+
+        if (lik<likelihood){
+            likelihood=lik;
+
+            bestPhase=dparam.phaseNonZer.getValuesPhase(); 
+
+        }
+        
+        dparam.save(path_result+"_"+lik+"_.json");
+
+            
+        
+        
+        dparam.phaseNonZer.setValuesPhase(bestPhase);
+        dparam.psf_fMany.updatePhase(dparam.phaseNonZer.getPointerPhase());
+        
+        ImageShow.imshow(dparam.psf_fMany.getPhase(),"phaseFinal");
+        for (int i=0;i<deltaZ.length;i++){
+            im[i]=dparam.psf_fMany.getPSF(i*7);
+        }
+        ImageShow.imshow(im,"PSF_Final");
+        
+        
+        IJ.log("finished");
+        
+        IJ.showProgress(1);
+        
+        
+        
+        
+    }
+    
+    private void phase_optim_old(int iterations,boolean init){
+        
+        
+        
+        double hdec=.05;
+        
+        
+        
+        
+        double lik=Double.MAX_VALUE;
+        int [] compl=this.dparam.phaseZer.complexity;
+        
+        int maxCompl=0;
+        for (int i=0;i<compl.length;i++){
+            if (maxCompl<compl[i]){
+                maxCompl=compl[i];
+            }
+        }
+        
+        
+        IJ.showProgress(.02);
+        if (init){
+            IJ.log("initialization started");
+            initPhaseGeneric(hdec);
+        }
+        
+        
+
+
+        
+        IJ.log("optimization started");
+        
+        
+        dparam.psf_fMany.updatePhase(dparam.phaseNonZer.getPointerPhase());
+        //ImageShow.imshow(dparam.psf_fMany.getPhase(),"phase "+0);
+        
+        for (int t=0;t<iterations;t++){
+            
+            
+            IJ.log("remaining iterations: "+(iterations-t));
+            IJ.showProgress(((double)t/(double)iterations));
+            
+            
+            lik=this.updatePhaseGeneric(hdec);
+                
+            
+            
+            
+            
+        
+            
+        }
+        
+        
+        IJ.showProgress(1);
+        
+        ImageShow.imshow(dparam.psf_fMany.getPhase(),"phase");
+        
+        //initPhasePistonOneValue(hdec);
+        /*for (int i=0;i<iterations;i++){
+            IJ.log("iter..."+(iterations-i));
+            this.updatePhasePistonOneValue(hdec);
+        }*/
+        
+        //ImageShow.imshow(dparam.psf_fMany.getPhase(),"phase after modif");
+        
+        
+        
+        double [][] crlb=this.compute_CRLB((float)hdec);
+        double [] xcrlb= new double [crlb.length];
+        double [] ycrlb= new double [crlb.length];
+        double [] zcrlb= new double [crlb.length];
+        
+        double meanX=0;
+        double meanY=0;
+        double meanZ=0;
+        
+        for (int i=0;i<crlb.length;i++){
+            xcrlb[i]=crlb[i][0];
+            ycrlb[i]=crlb[i][1];
+            zcrlb[i]=crlb[i][2];
+            meanX+=xcrlb[i];
+            meanY+=ycrlb[i];
+            meanZ+=zcrlb[i];
+        }
+        meanX/=(double)crlb.length;
+        meanY/=(double)crlb.length;
+        meanZ/=(double)crlb.length;
+        
+        //IJ.write(""+(deltaZ[deltaZ.length-1]-deltaZ[0])+" "+(deltaZ[1]-deltaZ[0])+" "+dparam.phaseZer.numCoef+" "+this.photonNumber+" "+this.background+" "+meanX+" "+meanY+" "+meanZ);
+        
+        //this.plot(deltaZ, xcrlb,"CRLB (X)","Z (µm)","CRLB (µm)");
+        //this.plot(deltaZ, ycrlb,"CRLB (Y)","Z (µm)","CRLB (µm)");
+        //this.plot(deltaZ, zcrlb,"CRLB (Z)","Z (µm)","CRLB (µm)");
+        
+        
+        
+    }
+    
+    
+    public void plot(double [] x, double [] y,String title,String xlabel,String ylabel){
+        Plot p = new Plot(""+title,xlabel,ylabel);
+        p.setFont(0, 18);
+        double xmin=Double.MAX_VALUE;
+        double xmax=Double.NEGATIVE_INFINITY;
+        double ymin=Double.MAX_VALUE;
+        double ymax=Double.NEGATIVE_INFINITY;
+        for (int i=0;i<x.length;i++){
+            if (xmin>x[i]){
+                xmin=x[i];
+            }
+            if (xmax<x[i]){
+                xmax=x[i];
+            }
+            if (ymin>y[i]){
+                ymin=y[i];
+            }
+            if (ymax<y[i]){
+                ymax=y[i];
+            }
+        }
+        if (ymin<.025){
+            ymin=0;
+            ymax=.025;
+        }
+        p.setLimits(xmin, xmax, ymin, ymax);
+        for (int ii=0;ii<x.length;ii++){
+            
+            p.setColor(Color.red);
+            p.add("CIRCLE",  x,y);
+            //p.show();
+            p.setLineWidth(2);
+        }
+        p.show();
+        
+    }
+    
+    
+    
+    private double updatePhasePistonOneValue(double h){
+        
+        
+        
+        dparam.phaseZer.setA(0, 1);
+        
+        
+        double lastLikelihood=-1;
+        
+        for (int i=0;i<dparam.phaseZer.nbDataPerImage;i++){
+        //for (int i=0;i<10;i++){
+            
+            double value=dparam.phaseZer.getPistonValue(i);
+            
+            
+            
+            dparam.psf_fMany.updatePhase(dparam.phaseZer.computeCombination());
+
+            double lik2=this.compute_sumstd_CRLB((float)h);
+            
+            
+            
+            
+
+            dparam.phaseZer.modifyOnePixelOfPiston(i, value-h);
+            
+            dparam.psf_fMany.updatePhase(dparam.phaseZer.computeCombination());
+            
+            double lik1=this.compute_sumstd_CRLB((float)h);
+
+
+
+
+
+
+            
+
+            dparam.phaseZer.modifyOnePixelOfPiston(i, value+h);
+            
+            dparam.psf_fMany.updatePhase(dparam.phaseZer.computeCombination());
+            
+            double lik3=this.compute_sumstd_CRLB((float)h);
+
+
+
+
+
+
+
+            double grad=0;
+            if (Math.abs((lik3+lik1-2*lik2)/(h*h))==0){
+                grad=((lik3-lik1)/(2*h));
+            }
+            else{
+                grad=((lik3-lik1)/(2*h))/Math.abs((lik3+lik1-2*lik2)/(h*h));
+            }
+
+
+            //IJ.log("lik("+i+") "+lik1+"  "+lik2+"  "+lik3+"  "+grad);
+
+
+            boolean found=false;
+            loop:for (double gamma=1;gamma>.02;gamma/=10){
+                
+                
+                dparam.phaseZer.modifyOnePixelOfPiston(i, value-gamma*grad);
+                dparam.psf_fMany.updatePhase(dparam.phaseZer.computeCombination());
+                double lik=this.compute_sumstd_CRLB((float)h);
+                
+                
+                
+                
+                if (lik<lik2){
+                    lastLikelihood=lik;
+                    found=true;
+                    break loop;
+                }
+                else{
+                    dparam.phaseZer.modifyOnePixelOfPiston(i, value);
+                }
+            }
+            
+            if (!found){
+                dparam.psf_fMany.updatePhase(dparam.phaseZer.computeCombination());
+            }
+            
+
+        }
+        
+        return lastLikelihood;
+        
+    }
+    
+    
+    
+    
+    
+    
+    private double initPhasePistonOneValue(double h){
+        
+        
+        
+        dparam.phaseZer.setA(0, 1);
+        
+        
+        double lastLikelihood=-1;
+        
+        for (int i=0;i<dparam.phaseZer.nbDataPerImage;i+=1){
+        //for (int i=0;i<10;i++){
+            IJ.log(""+(dparam.phaseZer.nbDataPerImage-i));
+            double lik=Double.MAX_VALUE;
+            
+            double v=dparam.phaseZer.getA(i);
+            
+            double pi=Math.PI;
+            
+            double finVal=-pi;
+            for (double value=v-pi;value<v+pi;value+=.1){
+                dparam.phaseZer.modifyOnePixelOfPiston(i, value);
+            
+                dparam.psf_fMany.updatePhase(dparam.phaseZer.computeCombination());
+
+                double liktmp=this.compute_sumstd_CRLB((float)h);
+                if (liktmp<lik){
+                    lik=liktmp;
+                    finVal=value;
+                }
+            }
+            
+            dparam.phaseZer.modifyOnePixelOfPiston(i, finVal);
+            
+
+        }
+        
+        return lastLikelihood;
+        
+    }
+    
+    
+    
+    
+    
+    private double updatePhaseZernike_std(double h,int numbCoef){
+        
+        
+        if (numbCoef>dparam.phaseZer.numCoef){
+            numbCoef=dparam.phaseZer.numCoef;
+        }
+        
+        double lastLikelihood=-1;
+        
+       
+        for (int p=3;p<numbCoef;p++){
+            
+            if (true){
+            
+                
+                dparam.psf_fMany.updatePhase(dparam.phaseZer.computeCombination(p, -h));
+                
+                double lik1=this.compute_sumstd_CRLB((float)h);
+                
+                
+                
+
+                
+                
+                
+                dparam.psf_fMany.updatePhase(dparam.phaseZer.computeCombination());
+                
+                double lik2=this.compute_sumstd_CRLB((float)h);
+                
+                
+                
+                
+                
+                
+                
+                dparam.psf_fMany.updatePhase(dparam.phaseZer.computeCombination(p, h));
+                
+                double lik3=this.compute_sumstd_CRLB((float)h);
+                
+                
+                
+                
+                
+                
+                
+                double save=dparam.phaseZer.getA(p);
+                double grad=0;
+                if (Math.abs((lik3+lik1-2*lik2)/(h*h))==0){
+                    grad=((lik3-lik1)/(2*h));
+                }
+                else{
+                    grad=((lik3-lik1)/(2*h))/Math.abs((lik3+lik1-2*lik2)/(h*h));
+                }
+                
+                
+                //IJ.log("lik("+p+") "+lik1+"  "+lik2+"  "+lik3+"  "+grad);
+                
+                
+                boolean found=false;
+                loop:for (double gamma=1;gamma>.02;gamma/=10){
+                    
+                    
+                    
+                    dparam.phaseZer.setA(p, save-gamma*grad);
+                    dparam.psf_fMany.updatePhase(dparam.phaseZer.computeCombination());
+                    double lik=this.compute_sumstd_CRLB((float)h);
+                
+                    
+                    
+                    
+                    if (lik<lik2){
+                        lastLikelihood=lik;
+                        found=true;
+                        break loop;
+                    }
+                    else{
+                        dparam.phaseZer.setA(p, save);
+                    }
+                }
+                if (!found){
+                    dparam.psf_fMany.updatePhase(dparam.phaseZer.computeCombination());
+                    
+                    
+                }
+            }
+        }
+        
+        return lastLikelihood;
+        
+    }
+    
+    
+    
+    
+    
+    
+    
+    private double updatePhaseGeneric(double h){
+        
+        
+        
+        double lastLikelihood=-1;
+        
+       
+        for (int p=0;p<dparam.param.sizeDisk;p++){
+            
+            if (true){
+                
+                double save=dparam.phaseNonZer.getValuePixel(p);
+                
+                
+                dparam.phaseNonZer.setValuePixel(p, save-h);
+                
+                dparam.psf_fMany.updatePhase(dparam.phaseNonZer.getPointerPhase());
+                
+                double lik1=this.compute_sum_CRLB((float)h);
+                
+                
+                
+
+                
+                
+                dparam.phaseNonZer.setValuePixel(p, save);
+                dparam.psf_fMany.updatePhase(dparam.phaseNonZer.getPointerPhase());
+                
+                double lik2=this.compute_sum_CRLB((float)h);
+                
+                
+                
+                
+                
+                
+                dparam.phaseNonZer.setValuePixel(p, save+h);
+                dparam.psf_fMany.updatePhase(dparam.phaseNonZer.getPointerPhase());
+                
+                double lik3=this.compute_sum_CRLB((float)h);
+                
+                
+                
+                
+                
+                
+                
+                double grad=0;
+                if (Math.abs((lik3+lik1-2*lik2)/(h*h))==0){
+                    grad=((lik3-lik1)/(2*h));
+                }
+                else{
+                    grad=((lik3-lik1)/(2*h))/Math.abs((lik3+lik1-2*lik2)/(h*h));
+                }
+                
+                
+                if (grad>Math.PI/10){
+                    grad=Math.PI/10;
+                }
+                if (grad<-Math.PI/10.){
+                    grad=-Math.PI/10.;
+                }
+                
+                //IJ.log("lik("+p+") "+lik1+"  "+lik2+"  "+lik3+"  "+grad);
+                
+                
+                boolean found=false;
+                loop:for (double gamma=1;gamma>.02;gamma/=10){
+                    
+                    
+                    dparam.phaseNonZer.setValuePixel(p, save-gamma*grad);
+                    dparam.psf_fMany.updatePhase(dparam.phaseNonZer.getPointerPhase());
+                    
+                    double lik=this.compute_sum_CRLB((float)h);
+                
+                    
+                    
+                    
+                    if (lik<lik2){
+                        lastLikelihood=lik;
+                        found=true;
+                        break loop;
+                    }
+                    else{
+                        dparam.phaseNonZer.setValuePixel(p, save);
+                        dparam.psf_fMany.updatePhase(dparam.phaseNonZer.getPointerPhase());
+                    }
+                }
+                if (!found){
+                    dparam.psf_fMany.updatePhase(dparam.phaseNonZer.getPointerPhase());
+                    
+                    
+                }
+            }
+        }
+        
+        return lastLikelihood;
+        
+    }
+    
+    
+    
+    
+    
+    
+    
+    
+    
+    
+    private double updatePhaseGeneric(double h,int angleSplit,int distSplit){
+        
+        h=.1;
+        
+        int splitNumber=angleSplit*distSplit;
+        
+        if ((dparam.phaseNonZer.angleNumber!=angleSplit)||(dparam.phaseNonZer.distNumber!=distSplit)){
+            dparam.phaseNonZer.createSplits(angleSplit, distSplit);
+        }
+        
+        
+        
+        double lastLikelihood=-1;
+        
+       
+        for (int p=0;p<splitNumber;p++){
+        
+            double [] save=dparam.phaseNonZer.getValuesPhase(p);
+            
+            if (save!=null){
+                
+                
+                
+                double [] tmp=new double [save.length];
+                
+                for (int i=0;i<tmp.length;i++){
+                    tmp[i]=save[i]-h;
+                    
+                }
+                dparam.phaseNonZer.setValuesPhase(p, tmp);
+                dparam.psf_fMany.updatePhase(dparam.phaseNonZer.getPointerPhase());
+                
+                
+                double lik1=this.compute_sum_CRLB((float)h);
+                
+                
+                
+
+                
+                
+                for (int i=0;i<tmp.length;i++){
+                    tmp[i]=save[i];
+                    
+                }
+                dparam.phaseNonZer.setValuesPhase(p, tmp);
+                dparam.psf_fMany.updatePhase(dparam.phaseNonZer.getPointerPhase());
+                
+                double lik2=this.compute_sum_CRLB((float)h);
+                
+                
+                
+                
+                
+                
+                for (int i=0;i<tmp.length;i++){
+                    tmp[i]=save[i]+h;
+                    
+                }
+                dparam.phaseNonZer.setValuesPhase(p, tmp);
+                dparam.psf_fMany.updatePhase(dparam.phaseNonZer.getPointerPhase());
+                
+                double lik3=this.compute_sum_CRLB((float)h);
+                
+                
+                
+                
+                
+                
+                
+                double grad=0;
+                if (Math.abs((lik3+lik1-2*lik2)/(h*h))==0){
+                    grad=((lik3-lik1)/(2*h));
+                }
+                else{
+                    grad=((lik3-lik1)/(2*h))/Math.abs((lik3+lik1-2*lik2)/(h*h));
+                }
+                
+                if (grad>Math.PI/10){
+                    grad=Math.PI/10;
+                }
+                if (grad<-Math.PI/10.){
+                    grad=-Math.PI/10.;
+                }
+                
+                
+                
+                
+                boolean found=false;
+                loop:for (double gamma=1;gamma>.02;gamma/=10){
+                    
+                    boolean isNan=false;
+                    for (int i=0;i<tmp.length;i++){
+                        tmp[i]=save[i]-gamma*grad;
+                        if (Double.isNaN(tmp[i])){
+                            isNan=true;
+                        }
+                    }
+                    if (!isNan){
+                        dparam.phaseNonZer.setValuesPhase(p, tmp);
+                        dparam.psf_fMany.updatePhase(dparam.phaseNonZer.getPointerPhase());
+
+                        double lik=this.compute_sum_CRLB((float)h);
+
+
+
+
+                        if (lik<lik2){
+                            lastLikelihood=lik;
+                            found=true;
+                            break loop;
+                        }
+                        else{
+                            dparam.phaseNonZer.setValuesPhase(p, save);
+                        }
+                    }
+                    
+                }
+                if (!found){
+                    dparam.phaseNonZer.setValuesPhase(p, save);
+                    dparam.psf_fMany.updatePhase(dparam.phaseNonZer.getPointerPhase());
+                    
+                    
+                }
+            }
+        }
+        
+        return lastLikelihood;
+        
+    }
+    
+    
+    
+    
+    
+    
+    
+    
+    private void initPhaseGeneric(double h){
+        
+        float [][][] im;
+        
+        IJ.log("initialization 1/6");
+        
+        
+        
+        
+        double lastLikelihood=-1;
+        
+        
+        /*for (int u=0;u<3;u++){
+            double v=Math.random()*3;
+            this.init(h, 4, 1, -v, v, v);
+            ImageShow.imshow(dparam.psf_fMany.getPhase(),"phase");
+            im = new float [deltaZ.length][][];
+            for (int i=0;i<deltaZ.length;i++){
+                im[i]=dparam.psf_fMany.getPSF(i*7);
+            }
+            ImageShow.imshow(im,"PSF1_"+u);
+            
+        }
+        
+        
+        //for (int u=0;u<1;u++){
+        //    double v=Math.random()*5;
+        //    double angle=Math.random()*2*Math.PI;
+        //    double radius=Math.random();
+        //    double dist=Math.max(1,Math.random()+radius);
+        //    this.initRing(h, angle, dist, radius, -v, v, v);
+        //    ImageShow.imshow(dparam.psf_fMany.getPhase(),"phase");
+        //    im = new float [deltaZ.length][][];
+        //    for (int i=0;i<deltaZ.length;i++){
+        //        im[i]=dparam.psf_fMany.getPSF(i*7);
+        //    }
+        //    ImageShow.imshow(im,"PSF2_"+u);
+        //}
+        
+        
+        
+        for (int u=0;u<3;u++){
+            double v=Math.random()*5;
+            this.init(h, 4, 1, -v, v, v);
+            ImageShow.imshow(dparam.psf_fMany.getPhase(),"phase");
+            im = new float [deltaZ.length][][];
+            for (int i=0;i<deltaZ.length;i++){
+                im[i]=dparam.psf_fMany.getPSF(i*7);
+            }
+            ImageShow.imshow(im,"PSF2_");
+            
+        }*/
+        
+        
+        //INITIALISE checking differents values for each split
+        int a=8;//angleSplit
+        int b=3;//distSplit
+        int splitNumber=a*b;
+        for (int u=0;u<5;u++){
+            IJ.log("WARNING... HERE, we should check all combinations at the same time");
+            dparam.phaseNonZer.createSplits(a, b);
+            for (int p=0;p<splitNumber;p++){
+                double [] save=dparam.phaseNonZer.getValuesPhase(p);
+                dparam.psf_fMany.updatePhase(dparam.phaseNonZer.getPointerPhase());
+                lastLikelihood=this.compute_sum_CRLB((float)h);
+
+                if (save!=null){
+
+                    double [] tmp=new double [save.length];
+                    double [] remember=new double [save.length];
+                    for (double t=-Math.PI;t<=Math.PI;t+=.5){
+                        for (int i=0;i<tmp.length;i++){
+                            tmp[i]=save[i]+t;
+
+                        }
+                        dparam.phaseNonZer.setValuesPhase(p, tmp);
+                        dparam.psf_fMany.updatePhase(dparam.phaseNonZer.getPointerPhase());
+
+
+                        double lik1=this.compute_sum_CRLB((float)h);
+
+                        if (lik1<lastLikelihood){
+                            for (int i=0;i<tmp.length;i++){
+                                remember[i]=tmp[i];
+
+                            }
+                            lastLikelihood=lik1;
+                        }
+                    }
+
+                    dparam.phaseNonZer.setValuesPhase(p, remember);
+                    dparam.psf_fMany.updatePhase(dparam.phaseNonZer.getPointerPhase());
+
+                }
+
+            }
+        }
+        
+        
+        ImageShow.imshow(dparam.psf_fMany.getPhase(),"phase");
+        im = new float [deltaZ.length][][];
+        for (int i=0;i<deltaZ.length;i++){
+            im[i]=dparam.psf_fMany.getPSF(i*7);
+        }
+        ImageShow.imshow(im,"PSF");
+        
+        
+        this.updatePhaseGeneric(.001,a,b);
+        
+        
+        
+        for (int u=0;u<20;u++){
+            int angleSplit=4+(int)(Math.random()*3);
+            int distSplit=2+(int)(Math.random()*3);
+            this.updatePhaseGeneric(.001,angleSplit,distSplit);
+            
+        }
+        
+        
+        ImageShow.imshow(dparam.psf_fMany.getPhase(),"phase");
+        im = new float [deltaZ.length][][];
+        for (int i=0;i<deltaZ.length;i++){
+            im[i]=dparam.psf_fMany.getPSF(i*7);
+        }
+        ImageShow.imshow(im,"PSF");
+        
+        
+        IJ.log("initialization 2/6");
+        
+        
+        for (int u=0;u<20;u++){
+            int angleSplit=4+(int)(Math.random()*12);
+            int distSplit=2+(int)(Math.random()*12);
+            this.updatePhaseGeneric(.001,angleSplit,distSplit);
+        }
+        
+        
+        ImageShow.imshow(dparam.psf_fMany.getPhase(),"phase");
+        im = new float [deltaZ.length][][];
+        for (int i=0;i<deltaZ.length;i++){
+            im[i]=dparam.psf_fMany.getPSF(i*7);
+        }
+        ImageShow.imshow(im,"PSF");
+        
+        
+        IJ.log("initialization 3/6");
+        
+        for (int u=0;u<20;u++){
+            int angleSplit=4+(int)(Math.random()*15);
+            int distSplit=2+(int)(Math.random()*15);
+            this.updatePhaseGeneric(.001,angleSplit,distSplit);
+        }
+        
+        
+        
+        ImageShow.imshow(dparam.psf_fMany.getPhase(),"phase");
+        im = new float [deltaZ.length][][];
+        for (int i=0;i<deltaZ.length;i++){
+            im[i]=dparam.psf_fMany.getPSF(i*7);
+        }
+        ImageShow.imshow(im,"PSF");
+        
+        
+        IJ.log("initialization 4/6");
+        
+        for (int u=0;u<15;u++){
+            int angleSplit=4+(int)(Math.random()*20);
+            int distSplit=2+(int)(Math.random()*20);
+            this.updatePhaseGeneric(.001,angleSplit,distSplit);
+        }
+        
+        
+        
+        
+        IJ.log("initialization 5/6");
+        
+        
+        for (int u=0;u<10;u++){
+            int angleSplit=10+(int)(Math.random()*25);
+            int distSplit=8+(int)(Math.random()*14);
+            this.updatePhaseGeneric(.001,angleSplit,distSplit);
+        }
+        
+        
+        
+        
+        
+        
+        
+        
+        
+        IJ.log("initialization 6/6");
+        
+        for (int u=0;u<5;u++){
+            int angleSplit=4+(int)(Math.random()*58);
+            int distSplit=20+(int)(Math.random()*12);
+            this.updatePhaseGeneric(.001,angleSplit,distSplit);
+        }
+        
+        
+        
+        
+    }
+    
+    
+    
+    private double init(double h,int angleSplit,int distSplit,double min,double max,double step){
+        
+        
+        
+        dparam.psf_fMany.updatePhase(dparam.phaseNonZer.getPointerPhase());
+        double likelihood=this.compute_sum_CRLB((float)h);
+        
+        
+        
+        
+        
+        int splitNumber=angleSplit*distSplit;
+        
+        
+        
+        //INITIALISE checking differents values for each split
+        dparam.phaseNonZer.createSplits(angleSplit, distSplit);
+        
+        
+        double [][] init=new double [splitNumber][];
+        for (int i=0;i<splitNumber;i++){
+            init[i]=dparam.phaseNonZer.getValuesPhase(i);
+        }
+        
+        
+        
+        int t=0;
+        for (double val=min;val<=max;val+=step){
+            t++;
+        }
+        double [] valstep = new double [t];
+        t=0;
+        for (double val=min;val<=max;val+=step){
+            //IJ.log("val "+val+"  "+min+"  "+max+"  "+step);
+            valstep[t]=val;
+            t++;
+        }
+        int numberOfCoef=splitNumber;
+        int numberOfVal=valstep.length;
+        
+        IJ.log("size   "+dparam.param.sizeDisk+" * "+Math.pow(numberOfVal,numberOfCoef));
+        double [][] coco = new double [(int)Math.round(Math.pow(numberOfVal,numberOfCoef))][dparam.param.sizeDisk];
+        for (int j=0;j<coco.length;j++){
+            coco[j]=dparam.phaseNonZer.getValuesPhase();
+        }
+        
+        int [] id = new int [numberOfCoef];
+        for (int i=0;i<id.length;i++){
+            id[i]=0;
+        }
+        for (int i=0;i<coco.length;i++){
+            
+            
+            //on compte en base numberOfVal
+            int number=i;
+            for (int b=0;b<id.length;b++){
+                id[b]=number%numberOfVal;
+                number=number/numberOfVal;
+            }
+            
+            for (int b=0,bb=0;b<splitNumber;b++,bb++){
+                int [] posit=dparam.phaseNonZer.getPositSplitPointer(b);
+                for (int u=0;u<posit.length;u++){
+                        coco[i][posit[u]]+=valstep[id[bb]];
+                }
+                
+            }
+        }
+        
+        int bestId=0;
+        
+        for (int p=0;p<coco.length;p++){
+            
+            
+            if (p%10==0){
+                IJ.log(""+(coco.length-p)/10);
+            }
+            
+            dparam.phaseNonZer.setValuesPhase(coco[p]);
+            dparam.psf_fMany.updatePhase(dparam.phaseNonZer.getPointerPhase());
+                
+            double lik=this.compute_sum_CRLB((float)h);
+            
+
+            if (lik<=likelihood){
+                likelihood=lik;
+                bestId=p;
+            }
+
+        }
+        
+        dparam.phaseNonZer.setValuesPhase(coco[bestId]);
+        dparam.psf_fMany.updatePhase(dparam.phaseNonZer.getPointerPhase());
+        
+        
+        return likelihood;
+    }
+    
+    
+    
+    
+    private double init_pi(double h,int angleSplit,int distSplit){
+        
+        IJ.log("init");
+        
+        dparam.psf_fMany.updatePhase(dparam.phaseNonZer.getPointerPhase());
+        double likelihood=this.compute_sum_CRLB((float)h);
+        
+        
+        
+        
+        
+        int splitNumber=angleSplit*distSplit;
+        
+        
+        
+        //INITIALISE checking differents values for each split
+        dparam.phaseNonZer.createSplits(angleSplit, distSplit);
+        
+        
+        IJ.log("a");
+        
+        double [][] phase=new double [splitNumber][];
+        double [][] bestPhase=new double[splitNumber][];
+        for (int i=0;i<splitNumber;i++){
+            phase[i]=dparam.phaseNonZer.getValuesPhase(i);
+            bestPhase[i]=dparam.phaseNonZer.getValuesPhase(i);
+        }
+        
+        
+        
+        IJ.log("lik init  "+likelihood);
+        
+        double pi=Math.PI;
+        
+        
+        
+        
+        fullLoop:while (true){
+            
+            loop:for (int i=0;i<splitNumber;i++){
+                if (phase[i][0]==0){
+                    for (int k=0;k<phase[i].length;k++){
+                        phase[i][k]=pi;
+                    }
+                    dparam.phaseNonZer.setValuesPhase(i,phase[i]);
+                    break loop;
+                }
+                else{
+                    for (int k=0;k<phase[i].length;k++){
+                        phase[i][k]=0;
+                        
+                    }
+                    dparam.phaseNonZer.setValuesPhase(i,phase[i]);
+                    if (i==splitNumber-1){
+                        break fullLoop;
+                    }
+                }
+                if (i==splitNumber-Math.min(splitNumber,9)){
+                    String s="";
+                    for (int j=i+1;j<splitNumber;j++){
+                        s+=" "+(phase[j][0]/pi);
+                        
+                    }
+                    IJ.log("wait...  "+s);
+                }
+            }
+            
+            
+            //dparam.phaseNonZer.setValuesPhase(phase);
+            dparam.psf_fMany.updatePhase(dparam.phaseNonZer.getPointerPhase());
+                
+            double lik=this.compute_sum_CRLB((float)h);
+            
+            //IJ.log("lik "+p+"  "+lik);
+
+            if ((lik<=likelihood)){
+                likelihood=lik;
+                IJ.log("better phase found");
+                for (int i=0;i<splitNumber;i++){
+                    for (int k=0;k<phase[i].length;k++){
+                        bestPhase[i][k]=phase[i][k];
+                    }
+                }
+                
+                dparam.save(path_result+"_"+lik+"_.json");
+            }
+            if (false){
+                break;
+            }
+
+        }
+        
+        for (int i=0;i<splitNumber;i++){
+            dparam.phaseNonZer.setValuesPhase(i,bestPhase[i]);
+        }
+        dparam.psf_fMany.updatePhase(dparam.phaseNonZer.getPointerPhase());
+        
+        
+        return likelihood;
+        
+        
+    }
+    
+    
+    
+    
+    private double init_pi(double h){
+        
+        IJ.log("size disk :"+dparam.param.sizeDisk);
+        int size=dparam.param.sizeDisk;
+        
+        dparam.psf_fMany.updatePhase(dparam.phaseNonZer.getPointerPhase());
+        
+        double likelihood=this.compute_sum_CRLB((float)h);
+        
+        
+        
+        
+        IJ.log("lik init  "+likelihood);
+        
+        double pi=Math.PI;
+        
+        
+        double [] phase=dparam.phaseNonZer.getValuesPhase();
+                
+        
+        
+        
+        double [] bestPhase=new double[size];
+        for (int i=0;i<size;i++){
+            bestPhase[i]=phase[i];
+        }
+        
+        
+        fullLoop:while (true){
+            
+            loop:for (int i=0;i<size;i++){
+                if (phase[i]==0){
+                    phase[i]=pi;
+                    break loop;
+                }
+                else{
+                    phase[i]=0;
+                    if (i==size-1){
+                        break fullLoop;
+                    }
+                }
+                if (i==size-9){
+                    String s="";
+                    for (int j=i+1;j<size;j++){
+                        s+=" "+(phase[j]/pi);
+                        
+                    }
+                    IJ.log("wait...  "+s);
+                }
+            }
+            
+            
+            dparam.phaseNonZer.setValuesPhase(phase);
+            dparam.psf_fMany.updatePhase(dparam.phaseNonZer.getPointerPhase());
+                
+            double lik=this.compute_sum_CRLB((float)h);
+            
+            //IJ.log("lik "+p+"  "+lik);
+
+            if ((lik<=likelihood)){
+                likelihood=lik;
+                IJ.log("better phase found");
+                for (int i=0;i<size;i++){
+                    bestPhase[i]=phase[i];
+                    
+                }
+                
+                dparam.save(path_result+"_"+lik+"_.json");
+            }
+            if (false){
+                break;
+            }
+
+        }
+        
+        dparam.phaseNonZer.setValuesPhase(bestPhase);
+        dparam.psf_fMany.updatePhase(dparam.phaseNonZer.getPointerPhase());
+        
+        
+        return likelihood;
+        
+        
+    }
+    
+    
+    
+    
+    
+    
+    private double initRing(double h,double angle,double dist,double radius,double min,double max,double step){
+        
+        
+        
+        dparam.psf_fMany.updatePhase(dparam.phaseNonZer.getPointerPhase());
+        double likelihood=this.compute_sum_CRLB((float)h);
+        
+        
+        
+        
+        
+        int splitNumber=2;
+        
+        
+        
+        //INITIALISE checking differents values for each split
+        dparam.phaseNonZer.createSplitsAs2Rings(angle, dist, radius);
+        
+        
+        double [][] init=new double [splitNumber][];
+        for (int i=0;i<splitNumber;i++){
+            init[i]=dparam.phaseNonZer.getValuesPhase(i);
+        }
+        
+        
+        
+        int t=0;
+        for (double val=min;val<=max;val+=step){
+            t++;
+        }
+        double [] valstep = new double [t];
+        t=0;
+        for (double val=min;val<=max;val+=step){
+            //IJ.log("val "+val+"  "+min+"  "+max+"  "+step);
+            valstep[t]=val;
+            t++;
+        }
+        int numberOfCoef=splitNumber;
+        int numberOfVal=valstep.length;
+        
+        
+        double [][] coco = new double [(int)Math.round(Math.pow(numberOfVal,numberOfCoef))][dparam.param.sizeDisk];
+        for (int j=0;j<coco.length;j++){
+            coco[j]=dparam.phaseNonZer.getValuesPhase();
+        }
+        
+        int [] id = new int [numberOfCoef];
+        for (int i=0;i<id.length;i++){
+            id[i]=0;
+        }
+        for (int i=0;i<coco.length;i++){
+            
+            
+            //on compte en base numberOfVal
+            int number=i;
+            for (int b=0;b<id.length;b++){
+                id[b]=number%numberOfVal;
+                number=number/numberOfVal;
+            }
+            
+            for (int b=0,bb=0;b<splitNumber;b++,bb++){
+                int [] posit=dparam.phaseNonZer.getPositSplitPointer(b);
+                IJ.log("error here "+posit+"  "+b);
+                for (int u=0;u<posit.length;u++){
+                    coco[i][posit[u]]+=valstep[id[bb]];
+                }
+                
+            }
+        }
+        
+        int bestId=0;
+        
+        for (int p=0;p<coco.length;p++){
+            
+            
+            
+            
+            dparam.phaseNonZer.setValuesPhase(coco[p]);
+            dparam.psf_fMany.updatePhase(dparam.phaseNonZer.getPointerPhase());
+                
+            double lik=this.compute_sum_CRLB((float)h);
+            
+
+            if (lik<=likelihood){
+                likelihood=lik;
+                bestId=p;
+            }
+
+        }
+        
+        dparam.phaseNonZer.setValuesPhase(coco[bestId]);
+        dparam.psf_fMany.updatePhase(dparam.phaseNonZer.getPointerPhase());
+        
+        
+        return likelihood;
+    }
+    
+    
+    
+    
+    
+    
+    double compute_sumstd_CRLB(float hdec){
+        double [][] crlb=compute_CRLB( hdec);
+        double som=0;
+        double somX=0;
+        double somY=0;
+        for (int i=0;i<crlb.length;i++){
+            som+=crlb[i][2];
+            somX+=crlb[i][0];
+            somY+=crlb[i][1];
+        }
+        
+        /////////////////begin modif
+        double som1=som/(double)crlb.length;
+        double som1X=somX/(double)crlb.length;
+        double som1Y=somY/(double)crlb.length;
+        double mean=som1;
+        double meanX=som1X;
+        double meanY=som1Y;
+        double som2=0;
+        double som2X=0;
+        double som2Y=0;
+        for (int i=0;i<crlb.length;i++){
+            som2+=(crlb[i][2]-mean)*(crlb[i][2]-mean);
+            som2X+=(crlb[i][0]-meanX)*(crlb[i][0]-meanX);
+            som2Y+=(crlb[i][1]-meanY)*(crlb[i][1]-meanY);
+        }
+        som2/=(double)crlb.length;
+        som2X/=(double)crlb.length;
+        som2Y/=(double)crlb.length;
+        
+        
+        som=som1+Math.sqrt(som2);
+        
+        ////////////////end modif
+        
+        
+        //add X & Y
+        //som+=som1X+Math.sqrt(som2X);
+        //som+=som1Y+Math.sqrt(som2Y);
+        
+        
+        //IJ.log("X Y Z");
+        
+        return som;
+    }
+    
+    
+    
+    
+    
+    double compute_sum_CRLB(float hdec){
+        double [][] crlb=compute_CRLB( hdec);
+        double som=0;
+        double somX=0;
+        double somY=0;
+        for (int i=0;i<crlb.length;i++){
+            som+=crlb[i][2];
+            somX+=crlb[i][0];
+            somY+=crlb[i][1];
+        }
+        
+        som/=(double)crlb.length;
+        somX/=(double)crlb.length;
+        somY/=(double)crlb.length;
+        
+        
+        
+        ////////////////begin modif to delete... max optim
+        double std=0;
+        float [][][] im = new float [deltaZ.length][][];
+        for (int i=0;i<deltaZ.length;i++){
+            im[i]=dparam.psf_fMany.getPSF(i*7);
+        }
+        
+        double stdcount=0;
+        for (int i=0;i<deltaZ.length;i++){
+            double count=0;
+            double moyX=0;
+            double moyY=0;
+            for (int ii=0;ii<im[i].length;ii++){
+                for (int iii=0;iii<im[i][ii].length;iii++){
+                    moyX+=im[i][ii][iii]*((double)(ii));
+                    moyY+=im[i][ii][iii]*((double)(iii));
+                    count+=im[i][ii][iii];
+                }
+            }
+            moyX/=count;
+            moyY/=count;
+            for (int ii=0;ii<im[i].length;ii++){
+                int sizeY=im[i][ii].length;
+                for (int iii=0;iii<im[i][ii].length;iii++){
+                    std+=im[i][ii][iii]*((double)(ii-moyX)*(double)(ii-moyX)+(double)(iii-moyY)*(double)(iii-moyY));
+                    stdcount+=im[i][ii][iii];
+                }
+            }
+        }
+        std/=stdcount;
+        std=Math.sqrt(std);
+        double lambda=0.001;
+        ////////////////end modif to delete
+        
+        
+        
+        //add X & Y
+        //som+=somX+somY;
+        
+        return som+std*lambda;
+    }
+    
+    
+    
+    double compute_sum_CRLB_save(float hdec){
+        double [][] crlb=compute_CRLB( hdec);
+        double som=0;
+        double somX=0;
+        double somY=0;
+        for (int i=0;i<crlb.length;i++){
+            som+=crlb[i][2];
+            somX+=crlb[i][0];
+            somY+=crlb[i][1];
+        }
+        
+        som/=(double)crlb.length;
+        somX/=(double)crlb.length;
+        somY/=(double)crlb.length;
+        
+        
+        
+        ////////////////begin modif to delete... max optim
+        /*double max=Double.NEGATIVE_INFINITY;
+        for (int i=0;i<crlbZ.length;i++){
+            if (crlbZ[i][2]>max){
+                max=crlbZ[i][2];
+            }
+        }
+        som=max;*/
+        ////////////////end modif to delete
+        
+        
+        
+        //add X & Y
+        //som+=somX+somY;
+        
+        return som;
+    }
+    
+    
+    
+    
+    
+    
+    
+    
+    
+    
+    double [][] compute_CRLB(float hdec){
+        
+        
+        int nbParam=5;
+        
+        //dparam.psf.computePSF(dparam.param.centerX,dparam.param.centerY,z);.
+        
+        double [][] crlb=new double[deltaZ.length][nbParam];
+        
+        double [] x = new double [deltaZ.length*7];
+        double [] y = new double [deltaZ.length*7];
+        double [] z = new double [deltaZ.length*7];
+        double [] zoil = new double [deltaZ.length*7];
+        
+        for (int i=0;i<deltaZ.length;i++){
+            x[i*7+0]=0;
+            y[i*7+0]=0;
+            z[i*7+0]=deltaZ[i];
+            
+            
+            
+            x[i*7+1]=+hdec;
+            y[i*7+1]=0;
+            z[i*7+1]=deltaZ[i];
+            
+            x[i*7+2]=-hdec;
+            y[i*7+2]=0;
+            z[i*7+2]=deltaZ[i];
+            
+            
+            
+            x[i*7+3]=0;
+            y[i*7+3]=+hdec;
+            z[i*7+3]=deltaZ[i];
+            
+            x[i*7+4]=0;
+            y[i*7+4]=-hdec;
+            z[i*7+4]=deltaZ[i];
+            
+            
+            x[i*7+5]=0;
+            y[i*7+5]=0;
+            z[i*7+5]=deltaZ[i]+hdec;
+            
+            x[i*7+6]=0;
+            y[i*7+6]=0;
+            z[i*7+6]=deltaZ[i]-hdec;
+            
+            
+        }
+        
+        for (int i=0;i<deltaZ.length;i++){
+            for (int ii=0;ii<7;ii++){
+                zoil[i*7+ii]=dparam.param.Zfocus;
+            }
+        }
+        
+        dparam.psf_fMany.computePSF(x, y, zoil,z);
+        
+        //dparam.psf_fMany.imshowFloat(dparam.psf_fMany.device_phase, "phase");
+        
+        
+        
+        double [] u=new double [nbParam];
+        double [] v=new double [nbParam];
+        
+        for (int t=0;t<deltaZ.length;t++){
+        //for (int t=0;t<1;t++){
+            
+            u[0]=0;
+            u[1]=0;
+            u[2]=deltaZ[t];
+            u[3]=photonNumber;
+            u[4]=background;
+            v[0]=0;
+            v[1]=0;
+            v[2]=deltaZ[t];
+            v[3]=photonNumber;
+            v[4]=background;
+            
+            
+            
+            float [][] f1 =dparam.psf_fMany.getPSF(t*7+0);
+            
+            
+            for (int i=0;i<f1.length;i++){
+                for (int ii=0;ii<f1[0].length;ii++){
+                    f1[i][ii]=f1[i][ii]*this.photonNumber+this.background;
+                }
+            }
+            
+            
+            float [][][] f0=new float [nbParam][f1.length][f1[0].length];
+            float [][][] f2=new float [nbParam][f1.length][f1[0].length];
+
+            
+                        
+            float [][] x2 =dparam.psf_fMany.getPSF(t*7+1);
+            
+            float [][] x0 =dparam.psf_fMany.getPSF(t*7+2);
+            for (int i=0;i<x0.length;i++){
+                for (int ii=0;ii<x0[0].length;ii++){
+                    f0[0][i][ii]=x0[i][ii]*this.photonNumber+this.background;
+                    f2[0][i][ii]=x2[i][ii]*this.photonNumber+this.background;
+                }
+            }
+            
+            
+            
+            float [][] y2 =dparam.psf_fMany.getPSF(t*7+3);
+            
+            float [][] y0 =dparam.psf_fMany.getPSF(t*7+4);
+            for (int i=0;i<y0.length;i++){
+                for (int ii=0;ii<y0[0].length;ii++){
+                    f0[1][i][ii]=y0[i][ii]*this.photonNumber+this.background;
+                    f2[1][i][ii]=y2[i][ii]*this.photonNumber+this.background;
+                }
+            }
+            
+            
+            float [][] z2 =dparam.psf_fMany.getPSF(t*7+5);
+            
+            float [][] z0 =dparam.psf_fMany.getPSF(t*7+6);
+            for (int i=0;i<z0.length;i++){
+                for (int ii=0;ii<z0[0].length;ii++){
+                    f0[2][i][ii]=z0[i][ii]*this.photonNumber+this.background;
+                    f2[2][i][ii]=z2[i][ii]*this.photonNumber+this.background;
+                }
+            }
+            
+            //double  [][] a2 =new double[f1.length][f1[0].length];
+            //double  [][] a0 =new double[f1.length][f1[0].length];
+            for (int i=0;i<f1.length;i++){
+                for (int ii=0;ii<f1[0].length;ii++){
+                    f0[3][i][ii]=f1[i][ii]*(this.photonNumber-hdec)+this.background;
+                    f2[3][i][ii]=f1[i][ii]*(this.photonNumber+hdec)+this.background;
+                }
+            }
+
+            
+            //double  [][] b2 =new double[f1.length][f1[0].length];
+            //double  [][] b0 =new double[f1.length][f1[0].length];
+            for (int i=0;i<f1.length;i++){
+                for (int ii=0;ii<f1[0].length;ii++){
+                    f0[4][i][ii]=f1[i][ii]*this.photonNumber+(this.background-hdec);
+                    f2[4][i][ii]=f1[i][ii]*this.photonNumber+(this.background+hdec);
+                }
+            }
+            
+            
+            //compute fisher
+            double d1,d2;
+            double [][] I=new double [nbParam][nbParam];
+            
+            for (int p=0;p<nbParam;p++){
+                //String s="";
+                for (int pp=0;pp<nbParam;pp++){
+                    I[p][pp]=0;
+                    for (int i=0;i<f1.length;i++){
+                        for (int ii=0;ii<f1[0].length;ii++){
+                            d1=(f2[p][i][ii]-f0[p][i][ii])/(2*Math.abs(hdec));
+                            d2=(f2[pp][i][ii]-f0[pp][i][ii])/(2*Math.abs(hdec));
+                            I[p][pp]+=(1/f1[i][ii])*d1*d2;
+                        }
+                    }
+                    //s+=" "+I[p][pp];
+                }
+                //IJ.log("r "+s);
+            }
+
+            double [] std =new double[nbParam];
+            boolean inversible=true;
+            Matrixe mat = new Matrixe(I);
+            try{
+                mat=Matrixe.inverse(mat);
+                I=mat.getMatrixe();
+            }catch(Exception ee){
+                
+                //dont take into account covar if non inversible
+                IJ.log("fisher matrix non inversible at z="+z);
+                for (int i=0;i<nbParam;i++){
+                    if (I[i][i]!=0){
+                        std[i]=Math.sqrt(1/I[i][i]);
+                    }
+                    else{
+                        std[i]=Double.MAX_VALUE;
+                    }
+                }
+                inversible=false;
+
+            }
+            
+            
+            
+            
+            if (inversible){
+                for (int i=0;i<nbParam;i++){
+                    std[i]=Math.sqrt(I[i][i]);
+                }
+            }
+            
+            crlb[t]=std;
+            
+        }
+        
+        if (Double.isNaN(crlb[0][0])){
+            IJ.log("ERROR: NanValue Phase optimization... maybe h is set too low -> no PSF difference ; h="+hdec);
+        }
+        return crlb;
+        
+    }
+    
+    
+    public void free(){ 
+        dparam.free();
+    }
+    
+    
+    
+}
diff --git a/src/main/java/org/pasteur/imagej/process/gpu/GenericPhaseOptimization_test.java b/src/main/java/org/pasteur/imagej/process/gpu/GenericPhaseOptimization_test.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/process/gpu/GenericPhaseRetrieval_.java b/src/main/java/org/pasteur/imagej/process/gpu/GenericPhaseRetrieval_.java
old mode 100644
new mode 100755
index 264e76c..58b3865
--- a/src/main/java/org/pasteur/imagej/process/gpu/GenericPhaseRetrieval_.java
+++ b/src/main/java/org/pasteur/imagej/process/gpu/GenericPhaseRetrieval_.java
@@ -375,7 +375,6 @@ private void phase_retrieve(int iterations){
                 IJ.showProgress((1.-(double)(iterations-t)/iterations));
                 
 
-                
 
                 this.updateRegistrationStacks(nbstack);
                 
@@ -429,8 +428,11 @@ private void phase_retrieve(int iterations){
                 
                 likelihood=lik;
                 
-                
-
+                if (t>0&&(t%10==0||t==iterations-1)){
+                    smoothPhase();
+                    showPhase("phase smoothed iter "+t);
+                    showImageAndModel("psf with smoothed phase iter "+t);
+                }
                 //IJ.log("sig "+dparam.param.sigmaGaussianKernel+"   wz "+dparam.param.getweightZ());
                 
 
@@ -829,7 +831,51 @@ public void plot(double [] x, double [] y,String title,String xlabel,String ylab
     
     
     
-    
+    private void smoothPhase(){
+        double [] phase= dparam.phaseNonZer.getValuesPhase();
+        double [] phaseMean= dparam.phaseNonZer.getValuesPhase();
+        int size=dparam.param.size;
+        int [][] ph=new int[size][size];
+        for (int i=0;i<size;i++){
+            for (int ii=0;ii<size;ii++){
+                ph[i][ii]=-1;
+            }
+        }
+        for (int p=0;p<dparam.param.sizeDisk;p++){
+            ph[dparam.param.disk2D[p][0]][dparam.param.disk2D[p][1]]=p;
+        }
+        for (int i=0;i<size;i++){
+            for (int ii=0;ii<size;ii++){
+                if (ph[i][ii]>=0){
+                    phaseMean[ph[i][ii]]=0;
+                    double count=0;
+                    for (int j=-1;j<=1;j++){
+                        for (int jj=-1;jj<=1;jj++){
+                            if (((i+j)>=0)&&((i+j)<size)&&((ii+jj)>=0)&&((ii+jj)<size)){
+                                if (ph[i+j][ii+jj]>=0){
+                                    double shift=0;
+                                    double minDist1=Math.abs(phase[ph[i][ii]]-phase[ph[i+j][ii+jj]]);
+                                    double minDist2=Math.abs(phase[ph[i][ii]]-phase[ph[i+j][ii+jj]]-Math.PI*2);
+                                    double minDist3=Math.abs(phase[ph[i][ii]]-phase[ph[i+j][ii+jj]]+Math.PI*2);
+                                    if (minDist2<minDist1){
+                                        shift=Math.PI*2;
+                                    }
+                                    if (minDist3<minDist1){
+                                        shift=-Math.PI*2;
+                                    }
+                                    
+                                    phaseMean[ph[i][ii]]+=phase[ph[i+j][ii+jj]]+shift;
+                                    count++;
+                                }
+                            }
+                        }
+                    }
+                    phaseMean[ph[i][ii]]/=count;
+                }
+            }
+        }
+        dparam.phaseNonZer.setValuesPhase(phaseMean);
+    }
     
     
     
@@ -1898,7 +1944,13 @@ private void updateSigmaGaussianKernel(int stackNumber){
 
         boolean found=false;
         loop:for (double gamma=1;gamma>.02;gamma/=10){
-            
+            double shift=0;
+            if (gamma*grad>0){
+                shift=Math.min(gamma*grad,0.05);
+            }
+            else{
+                shift=Math.max(gamma*grad,-0.05);
+            }
             dparam.param.sigmaGaussianKernel=save-gamma*grad;
             dparam.psf_fMany.updateSigmaGaussianKernel(dparam.param.sigmaGaussianKernel);
             double lik=0;
diff --git a/src/main/java/org/pasteur/imagej/process/gpu/GenericPhase_.java b/src/main/java/org/pasteur/imagej/process/gpu/GenericPhase_.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/process/gpu/InverseFiltering.java b/src/main/java/org/pasteur/imagej/process/gpu/InverseFiltering.java
old mode 100644
new mode 100755
index ec162a7..9a508ed
--- a/src/main/java/org/pasteur/imagej/process/gpu/InverseFiltering.java
+++ b/src/main/java/org/pasteur/imagej/process/gpu/InverseFiltering.java
@@ -148,7 +148,7 @@ public void run(double [][][] image){
 
             FastFourierTransform fftFilter = new FastFourierTransform(w,h);
             FastFourierTransform fftImage = new FastFourierTransform(w,h);
-
+            
             fftFilter.setReal(pad(psf[pp]));
             fftFilter.setImag(a);
             fftFilter.fft2D();
@@ -163,7 +163,8 @@ public void run(double [][][] image){
             
             
 
-            fftImage.setReal(addnoise(pad(image[0]),2000,10,2));
+            //fftImage.setReal(addnoise(pad(image[0]),2000,10,2));
+            fftImage.setReal(pad(image[0]));
             fftImage.setImag(imageB);
             fftImage.fft2D();
             imageA=fftImage.getPointerRealOut2D();
@@ -177,7 +178,7 @@ public void run(double [][][] image){
             for (int i=0;i<w;i++){
                 for (int ii=0;ii<h;ii++){
                     double tmp=(a[i][ii]*a[i][ii]+b[i][ii]*b[i][ii]);
-                    int method=2;
+                    int method=1;
                     if (method==1){
                         if (tmp!=0){
                             if (Math.sqrt((a[i][ii]*a[i][ii]+b[i][ii]*b[i][ii])/(tmp*tmp))<gamma){;
diff --git a/src/main/java/org/pasteur/imagej/process/gpu/LocalizationAutoFocusPipeline_.java b/src/main/java/org/pasteur/imagej/process/gpu/LocalizationAutoFocusPipeline_.java
old mode 100644
new mode 100755
index 36aad03..79affbd
--- a/src/main/java/org/pasteur/imagej/process/gpu/LocalizationAutoFocusPipeline_.java
+++ b/src/main/java/org/pasteur/imagej/process/gpu/LocalizationAutoFocusPipeline_.java
@@ -1153,7 +1153,7 @@ public void run(){
                                                             double my_crlbz=(1000*loc[partNumber][threadNumber].getCRLBZ());
                                                             
                                                             
-                                                            if (my_Score<4&&(my_crlbx<25)&&(my_crlby<25)&&(my_crlbz<50)){
+                                                            if ((my_crlbx<50)&&(my_crlby<50)&&(my_crlbz<100)){
                                                                 if (my_z>0){
                                                                     focus_likelihood_list[foc].add(my_Score);
                                                                 }
diff --git a/src/main/java/org/pasteur/imagej/process/gpu/LocalizationHD_.java b/src/main/java/org/pasteur/imagej/process/gpu/LocalizationHD_.java
old mode 100644
new mode 100755
index 98e40c3..f3bbbd1
--- a/src/main/java/org/pasteur/imagej/process/gpu/LocalizationHD_.java
+++ b/src/main/java/org/pasteur/imagej/process/gpu/LocalizationHD_.java
@@ -462,7 +462,7 @@ public FrameLocalization optimize(int nb_iter){
             JCuda.cudaStreamSynchronize(MyCudaStream.getCudaStream_t(streamId));
             int count=0;
             
-            //dp.psf.imshow(wh,height, device_model, "model ML","DOUBLE");//show first image
+            dp.psf.imshow(wh,height, device_model, "model ML","DOUBLE");//show first image
             
             for (int i=0;i<numberToCompute;i++){
                 if (true){//((Math.abs(position[i+index_x*numberToCompute])<dp.param.xystep*5)&&(Math.abs(position[i+index_y*numberToCompute])<dp.param.xystep*5)&&(position[i+index_z*numberToCompute]>minZ)&&(position[i+index_z*numberToCompute]<maxZ)){
diff --git a/src/main/java/org/pasteur/imagej/process/gpu/LocalizationMotionBlurPipeline_.java b/src/main/java/org/pasteur/imagej/process/gpu/LocalizationMotionBlurPipeline_.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/process/gpu/LocalizationMultiEmitterPipeline_.java b/src/main/java/org/pasteur/imagej/process/gpu/LocalizationMultiEmitterPipeline_.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/process/gpu/LocalizationMultiEmitter_.java b/src/main/java/org/pasteur/imagej/process/gpu/LocalizationMultiEmitter_.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/process/gpu/LocalizationPipelineHD_.java b/src/main/java/org/pasteur/imagej/process/gpu/LocalizationPipelineHD_.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/process/gpu/LocalizationPipelineHDsaveOld_.java b/src/main/java/org/pasteur/imagej/process/gpu/LocalizationPipelineHDsaveOld_.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/process/gpu/LocalizationPipeline_.java b/src/main/java/org/pasteur/imagej/process/gpu/LocalizationPipeline_.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/process/gpu/Localization_.java b/src/main/java/org/pasteur/imagej/process/gpu/Localization_.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/process/gpu/Localization_old.java b/src/main/java/org/pasteur/imagej/process/gpu/Localization_old.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/process/gpu/Model3D_.java b/src/main/java/org/pasteur/imagej/process/gpu/Model3D_.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/process/gpu/Model_generator_.java b/src/main/java/org/pasteur/imagej/process/gpu/Model_generator_.java
old mode 100644
new mode 100755
index 2006dd8..2f5850f
--- a/src/main/java/org/pasteur/imagej/process/gpu/Model_generator_.java
+++ b/src/main/java/org/pasteur/imagej/process/gpu/Model_generator_.java
@@ -145,10 +145,11 @@ public void computePSF(double [][] list,String path_result){
     
     
     
-    
-    
-    
     public void computePSFfromTable(StackLocalization sl,String path_result,boolean startAt0){
+        computePSFfromTable(sl,path_result,startAt0,-1);
+    }
+    
+    public void computePSFfromTable(StackLocalization sl,String path_result,boolean startAt0,int fixedimagesize){
         
         
         
@@ -186,10 +187,14 @@ public void computePSFfromTable(StackLocalization sl,String path_result,boolean
             minX=0;
             minY=0;
         }
-        
-        width=(int)Math.ceil(((maxX-minX)/1000)/dp.param.xystep)+dp.param.sizeoutput;
-        height=(int)Math.ceil(((maxY-minY)/1000)/dp.param.xystep)+dp.param.sizeoutput;
-        
+        if (fixedimagesize==-1){
+            width=(int)Math.ceil(((maxX-minX)/1000)/dp.param.xystep)+dp.param.sizeoutput;
+            height=(int)Math.ceil(((maxY-minY)/1000)/dp.param.xystep)+dp.param.sizeoutput;
+        }
+        else{
+            width=fixedimagesize;
+            height=fixedimagesize;
+        }
         
         IJ.log("min max "+minX+"  "+maxX+"  "+Math.ceil(((maxX-minX)/1000)/dp.param.xystep)+"  "+dp.param.xystep);
         IJ.log("wh "+width+"  "+height+"  "+dp.param.sizeoutput);
@@ -235,6 +240,13 @@ public void computePSFfromTable(StackLocalization sl,String path_result,boolean
         int [] frame = new int [model_number];
         double [] photon = new double [model_number];
         
+        if (fixedimagesize>=0){
+            if (dp.param.sizeoutput>fixedimagesize){
+
+                IJ.log("WARNING: image size extended to psf output size");
+
+            }
+        }
         
         
         if (dp.param.sizeoutput>width){
@@ -299,7 +311,7 @@ public void computePSFfromTable(StackLocalization sl,String path_result,boolean
         
         //path_result
         
-        org.pasteur.imagej.utils.ImageShow.imshow(res,"PSF generated");
+        //org.pasteur.imagej.utils.ImageShow.imshow(res,"PSF generated");
          
         
         ImageStack ims = new ImageStack(res[0].length,res[0][0].length);
diff --git a/src/main/java/org/pasteur/imagej/process/gpu/Modelmany_double_.java b/src/main/java/org/pasteur/imagej/process/gpu/Modelmany_double_.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/process/gpu/Modelmany_double_old.java b/src/main/java/org/pasteur/imagej/process/gpu/Modelmany_double_old.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/process/gpu/PSFCrossCorrel.java b/src/main/java/org/pasteur/imagej/process/gpu/PSFCrossCorrel.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/process/gpu/PSF_double_.java b/src/main/java/org/pasteur/imagej/process/gpu/PSF_double_.java
old mode 100644
new mode 100755
index 1a18e10..b0b9fe1
--- a/src/main/java/org/pasteur/imagej/process/gpu/PSF_double_.java
+++ b/src/main/java/org/pasteur/imagej/process/gpu/PSF_double_.java
@@ -245,14 +245,19 @@ public PSF_double_(PhaseParameters param){
                     
                     if (param.withApoFactor){
                         if (left_oil>right){
+                            
+                            
                             pupil[id]=(double)(1/Math.pow(1-(right/left_oil),.25));//with apodization factor
                         }
                         else{
                             pupil[id]=0;
+                            
+                            
                         }
                     }
                     else{
                         pupil[id]=1/(double)Math.sqrt(param.sizeDisk);//like that -> final sum=1 (unuseful actually because we normalize at the end)
+                        //pupil[id]=(double)-(i-center)/(0.5*(Math.sqrt(param.ringsize)*(param.size*param.xystep)))+.5;IJ.log("WARNING: apod remov for testing polarization");
                     }
                     
                     
@@ -973,6 +978,7 @@ public void computeGaussianKernel(double sigmaInPixels){
 
         for (int i=0;i<param.sizeDisk;i++){
             mat[param.disk2D[i][0]][param.disk2D[i][1]]=phase[i];
+            //IJ.write(""+i+","+param.disk2D[i][0]+","+param.disk2D[i][1]);
         }
             
         return mat;
diff --git a/src/main/java/org/pasteur/imagej/process/gpu/PSFgenerator_.java b/src/main/java/org/pasteur/imagej/process/gpu/PSFgenerator_.java
new file mode 100644
index 0000000..bc39fd8
--- /dev/null
+++ b/src/main/java/org/pasteur/imagej/process/gpu/PSFgenerator_.java
@@ -0,0 +1,338 @@
+/*
+ * To change this license header, choose License Headers in Project Properties.
+ * To change this template file, choose Tools | Templates
+ * and open the template in the editor.
+ */
+package org.pasteur.imagej.process.gpu;
+
+import org.pasteur.imagej.process.InitBackgroundAndPhotonNumber;
+import org.pasteur.imagej.process.gpu.DataPhase_;
+import org.pasteur.imagej.utils.PolynomialFit;
+import org.pasteur.imagej.utils.ImageShow;
+import org.pasteur.imagej.cuda.*;
+import jcuda.runtime.JCuda;
+
+import jcuda.runtime.cudaError;
+import ij.IJ;
+import java.awt.Color;
+import ij.gui.Plot;
+import jcuda.Pointer;
+import org.pasteur.imagej.utils.Zernike;
+
+/**
+ *
+ * @author benoit
+ */
+
+public class PSFgenerator_ {
+    //int stream=0;//not parallel process: iterative
+    double [] deltaZ;
+    //double [][][][] imageBlur;
+    int methodLikelihood=0;
+    int center;
+    DataPhase_ dparam;
+    
+    double [][] x;
+    double [][] y;
+    double [][] z;
+    int nbSlice;
+    int nbstack;
+    
+    double imageLength;
+    
+    //double [] fit_b_0;
+    //double [] fit_b_1;
+    double [] fit_a;
+    double [][] fit_a_each;
+    double [][][] background;
+    
+    double epsilon=.00001;//stop criterion
+    
+    double [][] registrationStack;//z,x,y registration for each stack
+    
+    String path_calibration;
+    
+    int fitOrderBackground=2;
+    double sigma;
+    double [] zwat ;
+    
+    double dx=0;//drift X supposed linear
+    double dy=0;//drift Y supposed linear
+    double dz=0;//drift Z supposed linear
+    
+    double maxDrift=10;//max drift between each frame (nm)
+    double axialRange,stepZ;
+    
+    public PSFgenerator_(int sizeFFT,double xystep,double zstep,double wavelength,double noil,double na ,String path_calibration,double sigma,int axialside,boolean withApoFactor,double axialRange, double stepZ){
+        this.sigma=sigma;
+        
+        
+        dparam = new DataPhase_(sizeFFT,sizeFFT,0,xystep,zstep,wavelength,noil,na,1.0,withApoFactor);
+        dparam.setNwat(dparam.param.noil);//ca ne change rien normalement car bille collée à lamelle
+        
+        dparam.param.Zfocus=0;
+        //astuce for registration
+        
+        this.axialRange=axialRange;
+        this.stepZ=stepZ;
+        
+        registrationStack=new double [3][nbstack];
+        
+        center=sizeFFT/2;
+        
+        
+        
+        deltaZ=new double [nbSlice];
+        //IJ.log("center is : "+center+"  !!!!!!!!!!§§§§§§§§§§§§§§§!!!!!!!");
+        for (int s=0;s<nbSlice;s++){
+            if (axialside==0){
+                deltaZ[s]=(((double)s)-center)*dparam.param.zstep;
+            }
+            else{
+                deltaZ[s]=-(((double)s)-center)*dparam.param.zstep;
+            }
+            
+        }
+        
+        zwat = new double [nbSlice];
+        for (int i=0;i<nbSlice;i++){
+            zwat[i]=0;
+        }
+        
+        x=new double[nbstack][nbSlice];//registration of each stack
+        y=new double[nbstack][nbSlice];//registration of each stack
+        z=new double[nbstack][nbSlice];//registration of each stack
+        
+        
+        
+        this.path_calibration=path_calibration;
+    }
+    
+    
+    
+    private void tophat(){
+        
+        //compute phase:
+        int nbDataPerImage=dparam.param.sizeDisk;
+
+        double maxDist=dparam.param.sizeRadiusRingPixel;
+
+        double center=dparam.param.centerFourierImage;
+        double dist;
+        
+        double [] phase = new double [nbDataPerImage];
+        
+        for (int i=0;i<nbDataPerImage;i++){
+            int x=dparam.param.disk2D[i][0];
+            int y=dparam.param.disk2D[i][1];
+
+
+            dist=Math.sqrt(((double)x-center)*((double)x-center)+((double)y-center)*((double)y-center));
+            
+
+            if (dist<0.8*maxDist){
+                phase[i]=0;
+            }
+            else {
+                phase[i]=3.141592;
+            }
+            
+        }
+        
+        dparam.phaseNonZer.setValuesPhase(phase);
+
+        Pointer devicePhase=dparam.phaseNonZer.getPointerPhase();
+        
+        dparam.psf.updatePhase(devicePhase);
+    }
+    
+    
+    
+    
+    
+    private void fermat(){
+        
+        
+        Zernike z = new Zernike(dparam.param.size,dparam.param.sizeRadiusRingPixel,28);
+        double [][][] zer=z.Z;
+        ImageShow.imshow(zer,"zer");
+        
+        ImageShow.imshow(zer[5],"zer");
+        //compute phase:
+        int nbDataPerImage=dparam.param.sizeDisk;
+
+        double a=dparam.param.sizeRadiusRingPixel;
+
+        double center=dparam.param.centerFourierImage;
+        double dist;
+        
+        double [] phase = new double [nbDataPerImage];
+        double order=1;
+        for (int i=0;i<nbDataPerImage;i++){
+            int x=dparam.param.disk2D[i][0];
+            int y=dparam.param.disk2D[i][1];
+            double xx=((double)x-center);
+            double yy=((double)y-center);
+            double angle=Math.atan2(yy, xx);
+            double r=xx*xx+yy*yy;
+            
+            
+            
+            // fermat
+            double dn=Double.POSITIVE_INFINITY;
+            double dp=Double.POSITIVE_INFINITY;
+            
+            for (double o=order%(int)order;o<=order;o++){
+                
+                double fp=a*a*(angle);
+                double fn=a*a*(angle+Math.PI);
+                if (Double.isNaN(fp)){
+                    IJ.log(""+angle);
+                }
+                if (r<fp){
+                    dp=fp-r;
+                }
+
+                if (r<fn){
+                    dn=fn-r;
+
+                }
+                //IJ.log("fn fp "+fn+"  "+fp);
+                if (dn<dp){
+                    phase[i]=0;//-Math.PI*Math.sqrt(r)/a;
+                }
+                else{
+                    phase[i]=1.5*Math.PI/2;
+                }
+                
+                
+            
+            }
+            //phase[i]=angle;
+            
+        }
+        
+        dparam.phaseNonZer.setValuesPhase(phase);
+
+        Pointer devicePhase=dparam.phaseNonZer.getPointerPhase();
+        
+        dparam.psf.updatePhase(devicePhase);
+    }
+    
+    
+    
+    
+    
+    
+    
+    private void split(){
+        
+        //compute phase:
+        int nbDataPerImage=dparam.param.sizeDisk;
+
+        double a=dparam.param.sizeRadiusRingPixel;
+
+        double center=dparam.param.centerFourierImage;
+        double dist;
+        
+        double [] phase = new double [nbDataPerImage];
+        double order=1;
+        for (int i=0;i<nbDataPerImage;i++){
+            int x=dparam.param.disk2D[i][0];
+            int y=dparam.param.disk2D[i][1];
+            double xx=((double)x-center);
+            double yy=((double)y-center);
+            double angle=Math.atan2(yy, xx);
+            double r=xx*xx+yy*yy;
+            
+            
+            
+            
+            if ((xx>0)&&(yy>0)){
+                phase[i]=angle;//-Math.PI*Math.sqrt(r)/a;
+            }
+            else if ((xx<0)&&(yy<0)){
+                phase[i]=angle;//-Math.PI*Math.sqrt(r)/a;
+            }
+            else{
+                phase[i]=Math.PI/2+angle;
+            }
+                
+                
+            
+            
+            //phase[i]=angle;
+            
+        }
+        
+        dparam.phaseNonZer.setValuesPhase(phase);
+
+        Pointer devicePhase=dparam.phaseNonZer.getPointerPhase();
+        
+        dparam.psf.updatePhase(devicePhase);
+    }
+    
+    
+    
+    
+    
+    
+    public void run(){
+        
+        
+        dparam.param.sigmaGaussianKernel=0;
+        
+        
+        
+        
+        //this.tophat();
+        //this.fermat();
+        
+        this.split();
+        // compute psf
+        
+        
+        
+        dparam.psf.resetKz();
+        
+        SearchPSFcenter_ spsfc= new SearchPSFcenter_(dparam,axialRange);
+        double position=spsfc.getPosition();
+        IJ.log("position"+position);
+        int nb=(int)((axialRange)/stepZ)+1;
+        
+        
+        double [][][] psf = new double [nb][][];
+        
+        
+        int k=0;
+        
+        for (double u=-axialRange/2+position;k<nb;u+=stepZ,k++){
+            dparam.psf.computePSF(0, 0,this.dparam.param.Zfocus,u);
+            psf[k]=dparam.psf.getPSF();
+        }
+        
+        ImageShow.imshow(psf,"psf");
+        
+        
+        
+        double [][] phaseres = dparam.psf.getPhase();
+        
+        
+        ImageShow.imshow(phaseres,"phaseres");
+        
+        
+    }
+
+
+    public void free(){
+        dparam.free();
+    }
+
+
+
+}
+    
+    
+    
+
+ 
\ No newline at end of file
diff --git a/src/main/java/org/pasteur/imagej/process/gpu/PSFmany_float_.java b/src/main/java/org/pasteur/imagej/process/gpu/PSFmany_float_.java
old mode 100644
new mode 100755
index 4a4f4ba..dc73ed6
--- a/src/main/java/org/pasteur/imagej/process/gpu/PSFmany_float_.java
+++ b/src/main/java/org/pasteur/imagej/process/gpu/PSFmany_float_.java
@@ -1344,6 +1344,30 @@ public void updatePhase(Pointer device_p){
     }
     
     
+    public void updatePhaseAndAdd(Pointer device_p,float [] adder){
+        
+        //this.imshowDouble(device_p, "phaseDouble");
+        
+        //double -> float
+        int cudaResult=JCuda.cudaStreamSynchronize(MyCudaStream.getCudaStream_t(param.stream));if (cudaResult != cudaError.cudaSuccess){IJ.log("ERROR synchro cuda gaussianKernel Set Image_0 "+cudaResult);}
+        MyVecDouble.double2float(custream,param.sizeDisk, device_phase, device_p);
+        cudaResult=JCuda.cudaStreamSynchronize(MyCudaStream.getCudaStream_t(param.stream));if (cudaResult != cudaError.cudaSuccess){IJ.log("ERROR synchro cuda gaussianKernel Set Image_1 "+cudaResult);}
+        
+        cudaMemcpyAsync(device_tmpFull, Pointer.to(adder), param.sizeDisk*Sizeof.FLOAT, cudaMemcpyHostToDevice,MyCudaStream.getCudaStream_t(param.stream));
+        
+        cudaResult=JCuda.cudaStreamSynchronize(MyCudaStream.getCudaStream_t(param.stream));if (cudaResult != cudaError.cudaSuccess){IJ.log("ERROR synchro cuda gaussianKernel Set Image_1 "+cudaResult);}
+        
+        MyVecDouble.add(custream, param.sizeDisk, device_phase, device_phase, device_tmpFull);
+        
+        cudaResult=JCuda.cudaStreamSynchronize(MyCudaStream.getCudaStream_t(param.stream));if (cudaResult != cudaError.cudaSuccess){IJ.log("ERROR synchro cuda gaussianKernel Set Image_1 "+cudaResult);}
+        
+        
+        //this.imshowFloat(device_phase, "phaseFloat");
+        
+        //this.imshowFloat(device_phase, "phaseInit");
+    }
+    
+    
 //    public void updatePupil(Pointer device_p){
 //        this.device_pupil=device_p;
 //        
diff --git a/src/main/java/org/pasteur/imagej/process/gpu/Predetection_cross_.java b/src/main/java/org/pasteur/imagej/process/gpu/Predetection_cross_.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/process/gpu/Predetection_deconvolutionFFT_.java b/src/main/java/org/pasteur/imagej/process/gpu/Predetection_deconvolutionFFT_.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/process/gpu/Predetection_deconvolutionFastSkellamBckgAccelerated_.java b/src/main/java/org/pasteur/imagej/process/gpu/Predetection_deconvolutionFastSkellamBckgAccelerated_.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/process/gpu/Predetection_deconvolutionFastSkellamBckg_.java b/src/main/java/org/pasteur/imagej/process/gpu/Predetection_deconvolutionFastSkellamBckg_.java
old mode 100644
new mode 100755
index 19e7a13..3411002
--- a/src/main/java/org/pasteur/imagej/process/gpu/Predetection_deconvolutionFastSkellamBckg_.java
+++ b/src/main/java/org/pasteur/imagej/process/gpu/Predetection_deconvolutionFastSkellamBckg_.java
@@ -1324,7 +1324,7 @@ public boolean setOinit(){
                if (withL1Norm){
                 for (int i=0;i<this.nbFrame;i++){
                     //minimise localizations
-                    MyVecDouble.addScalarFloat(MyCudaStream.getCUstream(streamId), (long)wh*(long)(range.length-1), device_tmpFFT2.withByteOffset((long)wh*(long)range.length*i*Sizeof.FLOAT), device_tmpFFT2.withByteOffset((long)wh*(long)range.length*i*Sizeof.FLOAT), (float)-.01);
+                    MyVecDouble.addScalarFloat(MyCudaStream.getCUstream(streamId), (long)wh*(long)(range.length-1), device_tmpFFT2.withByteOffset((long)wh*(long)range.length*i*Sizeof.FLOAT), device_tmpFFT2.withByteOffset((long)wh*(long)range.length*i*Sizeof.FLOAT), (float)-.05);
                     //maximize background
                     //MyVecDouble.addScalarFloat(MyCudaStream.getCUstream(streamId), (long)wh, device_tmpFFT2.withByteOffset((long)wh*(long)range.length*i*Sizeof.FLOAT+wh*(depth-1)*Sizeof.FLOAT), device_tmpFFT2.withByteOffset((long)wh*(long)range.length*i*Sizeof.FLOAT+wh*(depth-1)*Sizeof.FLOAT), (float)+.05);
                 
diff --git a/src/main/java/org/pasteur/imagej/process/gpu/Predetection_deconvolutionFastSkellam_.java b/src/main/java/org/pasteur/imagej/process/gpu/Predetection_deconvolutionFastSkellam_.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/process/gpu/Predetection_deconvolutionFast_.java b/src/main/java/org/pasteur/imagej/process/gpu/Predetection_deconvolutionFast_.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/process/gpu/Predetection_deconvolutionTest_.java b/src/main/java/org/pasteur/imagej/process/gpu/Predetection_deconvolutionTest_.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/process/gpu/SearchPSFcenter_.java b/src/main/java/org/pasteur/imagej/process/gpu/SearchPSFcenter_.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/process/gpu/SortTest.java b/src/main/java/org/pasteur/imagej/process/gpu/SortTest.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/process/gpu/Test.java b/src/main/java/org/pasteur/imagej/process/gpu/Test.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/process/gpu/ZernikePhaseOptimization_.java b/src/main/java/org/pasteur/imagej/process/gpu/ZernikePhaseOptimization_.java
old mode 100644
new mode 100755
index ef60ee4..662b090
--- a/src/main/java/org/pasteur/imagej/process/gpu/ZernikePhaseOptimization_.java
+++ b/src/main/java/org/pasteur/imagej/process/gpu/ZernikePhaseOptimization_.java
@@ -146,6 +146,67 @@ public void run(int nbIter){
         
     }
     
+    
+    private double [] fermat(){
+        
+        
+        
+        //compute phase:
+        int nbDataPerImage=dparam.param.sizeDisk;
+
+        double a=dparam.param.sizeRadiusRingPixel;
+
+        double center=dparam.param.centerFourierImage;
+        double dist;
+        
+        double [] phase = new double [nbDataPerImage];
+        double order=1;
+        for (int i=0;i<nbDataPerImage;i++){
+            int x=dparam.param.disk2D[i][0];
+            int y=dparam.param.disk2D[i][1];
+            double xx=((double)x-center);
+            double yy=((double)y-center);
+            double angle=Math.atan2(yy, xx);
+            double r=xx*xx+yy*yy;
+            
+            
+            
+            // fermat
+            double dn=Double.POSITIVE_INFINITY;
+            double dp=Double.POSITIVE_INFINITY;
+            
+            for (double o=order%(int)order;o<=order;o++){
+                
+                double fp=a*a*(angle);
+                double fn=a*a*(angle+Math.PI);
+                if (Double.isNaN(fp)){
+                    IJ.log(""+angle);
+                }
+                if (r<fp){
+                    dp=fp-r;
+                }
+
+                if (r<fn){
+                    dn=fn-r;
+
+                }
+                //IJ.log("fn fp "+fn+"  "+fp);
+                if (dn<dp){
+                    phase[i]=0;//-Math.PI*Math.sqrt(r)/a;
+                }
+                else{
+                    phase[i]=Math.PI;
+                }
+                
+                
+            
+            }
+            //phase[i]=angle;
+            
+        }
+        return phase;
+    }
+    
     private void phase_optim_zernike(int iterations,boolean init){
         
         
diff --git a/src/main/java/org/pasteur/imagej/process/gpu/ZernikePhaseRetrieval_.java b/src/main/java/org/pasteur/imagej/process/gpu/ZernikePhaseRetrieval_.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/process/gpu/ZernikePhase_.java b/src/main/java/org/pasteur/imagej/process/gpu/ZernikePhase_.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/utils/FRC.java b/src/main/java/org/pasteur/imagej/utils/FRC.java
old mode 100644
new mode 100755
index cc7c692..3e04355
--- a/src/main/java/org/pasteur/imagej/utils/FRC.java
+++ b/src/main/java/org/pasteur/imagej/utils/FRC.java
@@ -6,6 +6,7 @@
 package org.pasteur.imagej.utils;
 
 import ij.IJ;
+import ij.ImagePlus;
 import ij.gui.Plot;
 import java.awt.Color;
 import java.util.Arrays;
@@ -101,10 +102,26 @@ public int compare(double[] o1, double[] o2) {
         
     }
     
+    private double [][] smooth(double [][] x){
+        double [][] y = new double [x.length][x[0].length];
+        for (int i=0;i<x.length;i++){
+            for (int ii=0;ii<x[0].length;ii++){
+                for (int a=-1;a<=1;a++){
+                    for (int aa=-1;aa<=1;aa++){
+                        if ((i+a>=0)&&(i+a<x.length)&&(ii+aa>=0)&&(ii+aa<x[0].length)){
+                            y[i][ii]+=x[i+a][ii+aa];
+                        }
+                    }
+                }
+            }
+        }
+        return y;
+    }
+    
     
     public void run(){
         
-        int size=1024;
+        int size=512;
         IJ.log("SIZE "+size);
         double [][] weight = new double [w][h];
         double xtmp,ytmp;
@@ -125,7 +142,7 @@ public void run(){
                     ytmp=1;
                 }
                 
-                weight[i][ii]=1;//xtmp*ytmp*ztmp;
+                weight[i][ii]=xtmp*ytmp;
                 
             }
         }
@@ -142,6 +159,7 @@ public void run(){
         double [][] F = new double [size][size];
         double [][] G = new double [size][size];
         double [][] H = new double [size][size];
+        double [][] numeratorI = new double [size][size];
         double [][] tmpointer;
         FastFourierTransform fft1 = new FastFourierTransform(size,size);
         FastFourierTransform fft2 = new FastFourierTransform(size,size);
@@ -160,6 +178,8 @@ public void run(){
                         
                     }
                 }
+                C = new double [size][size];
+                D = new double [size][size];
 
                 fft1.setReal(A);
                 fft1.setImag(C);
@@ -173,12 +193,12 @@ public void run(){
 
                 tmpointer=fft1.getPointerRealOut2D();
                 addMatriceTo(E,tmpointer);
-                tmpointer=fft1.getPointerRealOut2D();
+                tmpointer=fft1.getPointerImagOut2D();
                 addMatriceTo(F,tmpointer);
 
                 tmpointer=fft2.getPointerRealOut2D();
                 addMatriceTo(G,tmpointer);
-                tmpointer=fft2.getPointerRealOut2D();
+                tmpointer=fft2.getPointerImagOut2D();
                 addMatriceTo(H,tmpointer);
                     
                 
@@ -186,14 +206,16 @@ public void run(){
         }
         
         
-        
+        //ImageShow.imshow(A,"inA");
+        //ImageShow.imshow(B,"inB");
         
         IJ.log("FFT OK");
         
         
         shift2D(A,E);
         shift2D(B,F);
-        
+//        ImageShow.imshow(A,"A");
+//        ImageShow.imshow(B,"B");
         IJ.log("SHIFT image 1 OK");
         
         shift2D(C,G);
@@ -207,19 +229,30 @@ public void run(){
 
         for (int i=0;i<size;i++){
             for (int ii=0;ii<size;ii++){
-                for (int iii=0;iii<size;iii++){
-                    E[i][ii]=(A[i][ii]*A[i][ii]+B[i][ii]*B[i][ii]);
-                    F[i][ii]=(C[i][ii]*C[i][ii]+D[i][ii]*D[i][ii]);
-                    G[i][ii]=(A[i][ii]*C[i][ii]+B[i][ii]*D[i][ii]);
-                    H[i][ii]=(B[i][ii]*C[i][ii]-A[i][ii]*D[i][ii]);
-                }
+                
+                E[i][ii]=(A[i][ii]*A[i][ii]+B[i][ii]*B[i][ii]);
+                F[i][ii]=(C[i][ii]*C[i][ii]+D[i][ii]*D[i][ii]);
+                G[i][ii]=(A[i][ii]*C[i][ii]+B[i][ii]*D[i][ii]);
+                H[i][ii]=(B[i][ii]*C[i][ii]-A[i][ii]*D[i][ii]);
+                numeratorI[i][ii]=(G[i][ii]*G[i][ii])+(H[i][ii]*H[i][ii]);
+                    
+                
             }
         }
         
+        /*ImagePlus imptmp = new ImagePlus("/home/benoit/atlas/gaia/Mickael/Dual_Objective/20191128_NUP96SNAP_65ms_200_400/3D_FoV__1/analysis1/FSC_benoit/test/numerator.tif");
+        
+        double [][] imptmpp = new double [imptmp.getWidth()][imptmp.getHeight()];
+        for (int i=0;i<imptmp.getWidth();i++){
+            for (int ii=0;ii<imptmp.getHeight();ii++){
+                imptmpp[i][ii]=imptmp.getProcessor().getPixelValue(i, ii);
+            }
+        }*/
+        
         IJ.showProgress(.99);
         
         IJ.log("NORMALIZED CORRELATION OK");
-//        ImageShow.imshow(A,"A");
+        //ImageShow.imshow(numeratorI,"numeratorI");
 //        ImageShow.imshow(B,"B");
 //        ImageShow.imshow(C,"C");
 //        ImageShow.imshow(D,"D");
@@ -234,10 +267,34 @@ public void run(){
         //A contains distance matrix
         IJ.log("DIST MATRIX OK");
         
+        E=smooth(E);
+        F=smooth(F);
+        G=smooth(G);
+        H=smooth(H);
+        E=smooth(E);
+        F=smooth(F);
+        G=smooth(G);
+        H=smooth(H);
+        E=smooth(E);
+        F=smooth(F);
+        G=smooth(G);
+        H=smooth(H);
+        E=smooth(E);
+        F=smooth(F);
+        G=smooth(G);
+        H=smooth(H);
+        E=smooth(E);
+        F=smooth(F);
+        G=smooth(G);
+        H=smooth(H);
+        
+        ImageShow.imshow(E,"E");
+        
+        double step=1;//pixel step to convert in vector//
+        int sizeVect=(int)Math.ceil(size*Math.sqrt(2)/(2*step));//sqrt(2) because 2D
         
+        double [] numerator = new double[sizeVect];
         
-        double step=1;//pixel step to convert in vector
-        int sizeVect=(int)Math.ceil(size*Math.sqrt(2)/(2*step));//sqrt(2) because 2D
         double [] vectE = new double[sizeVect];
         double [] vectF = new double[sizeVect];
         double [] vectG = new double[sizeVect];
@@ -247,6 +304,8 @@ public void run(){
         
         //sum over distances to convert 2D images to vectors
         for (int i=0;i<v.length;i++){
+            
+            
             int bin=(int)(v[i][0]/step);
             
             
@@ -256,26 +315,35 @@ public void run(){
             vectH[bin]+=H[(int)v[i][1]][(int)v[i][2]];
             vectN[bin]++;
             
+            numerator[bin]+=G[(int)v[i][1]][(int)v[i][2]]*G[(int)v[i][1]][(int)v[i][2]]+H[(int)v[i][1]][(int)v[i][2]]*H[(int)v[i][1]][(int)v[i][2]];
+            
         }
         
-        
-        
-        
+        plot(numerator,"numerator");
+        plot(vectE,"vectE");
+        plot(vectF,"vectF");
         
         double [] vectFinal = new double[sizeVect];
         double [] xaxisXY = new double[sizeVect];
-        for (int i=0;i<sizeVect;i++){
-            vectFinal[i]=(vectG[i]/Math.sqrt(vectE[i]*vectF[i]));// - 2/Math.sqrt(vectN[i]/2);
-            xaxisXY[i]=i*step/(size*this.pixSizeXY);
+        for (int i=1;i<sizeVect;i++){
+            vectFinal[i]=(Math.sqrt(numerator[i])/Math.sqrt(vectE[i]*vectF[i]));// - 2/Math.sqrt(vectN[i]/2);
+            xaxisXY[i]=i*step/(size);//*this.pixSizeXY);
+            IJ.log("xax "+xaxisXY[i]+"  "+size+"  "+pixSizeXY);
         }
         
+        vectFinal[0]=1;
+        xaxisXY[0]=0;
+        
+        IJ.log("PIXSIZ="+pixSizeXY+"  "+xaxisXY[sizeVect-1]);
+        
+        
         double [] vectFinalSmooth = new double[sizeVect];
         
         double m,c;
         for (int i=0;i<sizeVect;i++){
             m=0;
             c=0;
-            for (int u=-5;u<=5;u++){
+            for (int u=-10;u<=10;u++){
                 if ((u+i>=0)&&(u+i<sizeVect)){
                     m+=vectFinal[i+u];
                     c++;
@@ -287,13 +355,19 @@ public void run(){
         double threshold=1./7.;
         
         int xresol=0;
-        while (vectFinalSmooth[xresol]>threshold){
+        looper:while (vectFinalSmooth[xresol]>threshold){
             xresol++;
+            if (xresol>=vectFinalSmooth.length-1){
+                break looper;
+            }
         }
-        IJ.log("r "+xaxisXY[xresol]+"  "+xresol);
+        //IJ.log("r "+xaxisXY[xresol]+"  "+xresol);
         double resolutionXY=xaxisXY[xresol];
         
         
+        
+        //plot(xaxisXY,vectN,vectFinalSmooth,xaxisXY[xresol],threshold,"N","1/nm","N");
+        
         plot(xaxisXY,vectFinal,vectFinalSmooth,xaxisXY[xresol],threshold,"FRC","1/nm","FRC");
         
         
@@ -358,4 +432,35 @@ public void plot(double [] x, double [] y, double [] y2,double lineX,double line
     
     
     
+    
+    public void plot(double [] y,String title){
+        
+        double yMin=y[0]; double yMax=y[0];
+        double [] x = new double [y.length];
+        for (int i=0;i<y.length;i++){
+            x[i]=i;
+            if (y[i]<yMin){
+                yMin=y[i];
+            }
+            if (y[i]>yMax){
+                yMax=y[i];
+            }
+            
+        }
+        
+        
+        Plot p = new Plot(""+title,"","");
+        p.setLimits(0, x.length, yMin-.01*(yMax-yMin), yMax+.01*(yMax-yMin));
+        p.setLineWidth(2);
+        p.setFont(0, 18);
+        
+        p.addPoints(x, y,Plot.LINE);
+        p.show();
+        
+    }
+    
+    
+    
+    
+    
 }
diff --git a/src/main/java/org/pasteur/imagej/utils/FRCbiop.java b/src/main/java/org/pasteur/imagej/utils/FRCbiop.java
new file mode 100755
index 0000000..2f93ad5
--- /dev/null
+++ b/src/main/java/org/pasteur/imagej/utils/FRCbiop.java
@@ -0,0 +1,851 @@
+/*
+ * To change this license header, choose License Headers in Project Properties.
+ * To change this template file, choose Tools | Templates
+ * and open the template in the editor.
+ */
+package org.pasteur.imagej.utils;
+
+
+
+import ij.IJ;
+import ij.ImagePlus;
+import ij.ImageStack;
+import ij.gui.Plot;
+import ij.gui.PlotWindow;
+import ij.measure.ResultsTable;
+import ij.process.FHT;
+import ij.process.FloatProcessor;
+import ij.process.ImageProcessor;
+
+import java.awt.Color;
+import java.io.File;
+import java.io.FilenameFilter;
+import java.util.ArrayList;
+import java.util.Arrays;
+
+
+import org.apache.commons.math3.analysis.interpolation.LoessInterpolator;
+import org.apache.commons.math3.util.FastMath;
+
+/**
+ * @author Alex Herbert
+ * @author Olivier Burri
+ * 
+ * Copyright (C) 2016  Alex Herbert
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/**
+ * Compute the Fourier Ring Correlation, a measure of the resolution of a microscopy image.
+ * <p>
+ * Adapted from the FIRE (Fourier Image REsolution) plugin produced as part of the paper:
+ * <a href="http://www.nature.com/nmeth/journal/v10/n6/full/nmeth.2448.html">
+ * Niewenhuizen, et al (2013). Measuring image resolution in optical nanoscopy. Nature Methods, 10, 557</a>
+ * 
+ */
+public class FRCbiop {
+    
+
+
+
+	public enum ThresholdMethod
+	{
+		FIXED_1_OVER_7("Fixed 1/7"), HALF_BIT("Half-bit"), THREE_SIGMA("Three sigma");
+
+		private String name;
+
+		private ThresholdMethod(String name)
+		{
+			this.name = name;
+		}
+
+		@Override
+		public String toString()
+		{
+			return name;
+		}
+	};
+
+	// Properties controlling the algorithm
+
+	/**
+	 * The correlation is computed using intervals around the circle circumference. The number of samples for half the
+	 * circle is computed as: Pi * radius * sampling factor
+	 */
+	public double perimeterSamplingFactor = 1;
+
+	/**
+	 * The correlation is computed using intervals around the circle circumference of the Fourier transform. The Fourier
+	 * image is 2-fold radially symmetric and so the calculation can use only half the circle for speed. Note: The
+	 * results will differ slightly due to the implementation of the Fourier image not being exactly symmetric and the
+	 * sample points used on the perimeter not matching between the two semi-circles.
+	 */
+	public boolean useHalfCircle = true;
+
+        
+        public void test(ImageProcessor ip1, ImageProcessor ip2,double pixSize){
+            
+            
+            double[][] frc_curve = calculateFrcCurve(ip1, ip2);
+
+
+                // Smooth Curve using LOESS
+                double[][] smooth_frc = getSmoothedCurve(frc_curve);
+                
+                // Fourier Image REsolution number ("FIRE")
+                double fire = calculateFireNumber(smooth_frc, ThresholdMethod.FIXED_1_OVER_7);
+                
+                
+                IJ.log("FRC "+ fire);
+                IJ.log("FRC (nm) "+ fire*pixSize);
+                // Plot if requested
+                Plot p = doPlot(frc_curve, smooth_frc, ThresholdMethod.FIXED_1_OVER_7, fire, "result FRC biop");
+                p.show();
+        }
+        
+	/**
+	 * Calculate the Fourier Ring Correlation curve for two images.
+	 * 
+	 * @param ip1
+	 * @param ip2
+	 * @return An array of triplets representing [][radius,correlation,N] where correlation is the FRC at the given
+	 *         radius from the centre of the Fourier transform (i.e. 1/spatial frequency) and N is the number of samples
+	 *         used to compute the correlation
+	 */
+	public double[][] calculateFrcCurve(ImageProcessor ip1, ImageProcessor ip2)
+	{
+		IJ.showStatus("Calculating complex FFT images...");
+                
+                float[] input1 = (float[]) ip1.getPixels();
+                //ImageShow.imshow(input1,ip1.getWidth(),ip1.getHeight(),"input1");
+
+		// Pad images to the same size
+		final int maxWidth = FastMath.max(ip1.getWidth(), ip2.getWidth());
+		final int maxHeight = FastMath.max(ip1.getHeight(), ip2.getHeight());
+		ip1 = mypad(ip1, maxWidth, maxHeight);IJ.log("WARNING: I use fixed padding");
+		ip2 = mypad(ip2, maxWidth, maxHeight);
+                
+                float[] input2 = (float[]) ip1.getPixels();
+                //ImageShow.imshow(input2,ip1.getWidth(),ip1.getHeight(),"input2");
+                
+		// TODO - Can this be sped up by computing both FFT transforms together?
+		FloatProcessor[] fft1 = getComplexFFT(ip1);
+		FloatProcessor[] fft2 = getComplexFFT(ip2);
+
+		IJ.showStatus("Preparing FRC curve calculation...");
+
+		final int size = fft1[0].getWidth();
+
+		// In-line for speed
+		float[] numerator = new float[size * size];
+		float[] absFFT1 = new float[numerator.length];
+		float[] absFFT2 = new float[numerator.length];
+
+		float[] dataA1 = (float[]) fft1[0].getPixels();
+		float[] dataB1 = (float[]) fft1[1].getPixels();
+		float[] dataA2 = (float[]) fft2[0].getPixels();
+		float[] dataB2 = (float[]) fft2[1].getPixels();
+                
+                
+		for (int y = size, i = size * size - 1; y-- > 0;)
+		{
+			for (int x = size; x-- > 0; i--)
+			{
+				// A remark from Brian Patton, University of Strathclyde, Glasgow says that the original
+				// computation for the numerator was wrong, namely:
+				// numerator[i] = dataA1[i] * dataA2[i] + dataB1[i] * dataB2[i];
+
+				// And should be instead:
+				// abs((FFT(Image1)*Complex_Conj(FFT(Image2)))
+				// Code is modified as follows.
+				float numrel =  dataA1[i] * dataA2[i] +  dataB1[i] * dataB2[i];
+				float numim  =  dataB1[i] * dataA2[i] -  dataA1[i] * dataB2[i];
+				// New numerator
+				numerator[i] = numrel * numrel + numim * numim;
+
+				absFFT1[i] = dataA1[i] * dataA1[i] + dataB1[i] * dataB1[i];
+				absFFT2[i] = dataA2[i] * dataA2[i] + dataB2[i] * dataB2[i];
+			}
+		}
+                
+                ImageShow.imshow(dataA1,size,size,"A1");
+//                ImageShow.imshow(dataA2,size,size,"A2");
+                ImageShow.imshow(numerator,size,size,"numerator");
+		int radius = 1;
+		final double centre = size / 2;
+                IJ.log("WARNING: radius is fixed");
+	    final double max = (Math.max(maxWidth, maxHeight)/2) - 1;
+		IJ.showStatus("Calculating FRC curve... radiusMax="+max);
+
+		double[][] frcCurve = new double[(int) max][3];
+
+		// Radius zero is always 1
+		frcCurve[0][0] = 0;
+		frcCurve[0][1] = 1;
+		frcCurve[0][2] = 1; // Avoid divide by zero errors. Not sure if this is OK
+                double [] numerator_ = new double[(int)max];
+                double [] denominator1_ = new double[(int)max];
+                double [] denominator2_ = new double[(int)max];
+                double [] number_ = new double[(int)max];
+		float[][] images = new float[][] { numerator, absFFT1, absFFT2 };
+                
+		while (radius < max)
+		{
+			final double progress = (1.0 * radius) / max;
+			IJ.showProgress(progress);
+			IJ.log("Calculating FRC curve...[Radius = " + radius + "px]");
+
+			// Inline the calculation for speed
+			double sum1 = 0;
+			double sum2 = 0;
+			double sum3 = 0;
+
+			// Circumference is 2*pi*r
+		    double angleStep = 1 / (perimeterSamplingFactor * radius);
+
+			double angle = 0D;
+			int numSum = 0;
+
+			final double limit = (useHalfCircle) ?  Math.PI :  2 * Math.PI;
+			while (angle < limit)
+			{
+				double x = centre + radius * Math.cos(angle);
+				double y = centre + radius * Math.sin(angle);
+
+				double[] values = getInterpolatedValues(x, y, images, size);
+				sum1 += values[0];
+                                
+				sum2 += values[1];
+				sum3 += values[2];
+
+				numSum++;
+				angle += angleStep;
+			}
+                        
+			double val = Math.sqrt(sum1) / Math.sqrt(sum2 * sum3);
+                        numerator_[radius]=sum1;
+                        denominator1_[radius]=sum2;
+                        denominator2_[radius]=sum3;
+                        number_[radius]=numSum;
+			frcCurve[radius][0] = radius;
+			frcCurve[radius][1] = val;
+			frcCurve[radius][2] = numSum;
+
+			radius++;
+		}
+                plot(numerator_,"numerator_");
+                plot(denominator1_,"denominator1_");
+                plot(denominator2_,"denominator2_");
+                plot(number_,"number_");
+		IJ.showProgress(1);
+		IJ.showStatus("Finished calculating FRC curve...");
+		return frcCurve;
+	}
+
+	private ImageProcessor pad(ImageProcessor ip, int width, int height)
+	{
+		// Make sure width and height are even
+		width = width%2 == 0 ? width : width+1;
+		height = height%2 == 0 ? height : height+1;
+		if (ip.getWidth() != width || ip.getHeight() != height)
+		{
+			ImageProcessor ip2 = ip.createProcessor(width, height);
+			ip2.insert(ip, 0, 0);
+			ip = ip2;
+		}
+		return ip;
+	}
+        
+        
+        private ImageProcessor mypad(ImageProcessor ip, int width, int height)
+	{
+		// Make sure width and height are even
+		width = width%2 == 0 ? width : width+1;
+		height = height%2 == 0 ? height : height+1;
+                width=512;
+                height=512;
+		if (ip.getWidth() != width || ip.getHeight() != height)
+		{
+			ImageProcessor ip2 = ip.createProcessor(width, height);
+			ip2.insert(ip, 128 , 128);
+			ip = ip2;
+		}
+		return ip;
+	}
+
+
+	/**
+	 * Convert an image into a Fourier image with real and imaginary parts
+	 * 
+	 * @param ip
+	 * @return the real and imaginary parts
+	 */
+	public FloatProcessor[] getComplexFFT(ImageProcessor ip)
+	{
+		FloatProcessor taperedDataImage = getSquareTaperedImage(ip);
+                //ImagePlus imp = new ImagePlus("toto",taperedDataImage);
+                //imp.show();
+                IJ.log("WARNING: taperedDataImage looks weird");
+		FHT fht = new FHT(taperedDataImage);//FHT fht = new FHT(ip);
+		fht.transform();
+
+		FloatProcessor[] ret = new FloatProcessor[2];
+
+		ImageStack stack1 = fht.getComplexTransform();
+		ret[0] = ((FloatProcessor) stack1.getProcessor(1));
+		ret[1] = ((FloatProcessor) stack1.getProcessor(2));
+
+		return ret;
+	}
+
+	// Cache the Tukey window function
+	private float[] taperX = new float[0];
+	private float[] taperY = new float[0];
+
+	/**
+	 * Applies a Tukey window function to the image and then pads it to the next square size power of two.
+	 * 
+	 * @param dataImage
+	 * @return The square tapered image
+	 */
+	public FloatProcessor getSquareTaperedImage(ImageProcessor dataImage)
+	{
+		// Use a Tukey window function
+		taperX = getWindowFunction(taperX, dataImage.getWidth());
+		taperY = getWindowFunction(taperY, dataImage.getHeight());
+
+		final int size = FastMath.max(dataImage.getWidth(), dataImage.getHeight());
+
+		// Pad to a power of 2
+		int newSize = 0;
+		for (int i = 4; i < 15; i++)
+		{
+			newSize = (int) Math.pow(2.0, i);
+			if (size <= newSize)
+			{
+				break;
+			}
+		}
+
+		if (size > newSize)
+			return null; // Error
+
+		dataImage = dataImage.toFloat(0, null);
+		float[] data = (float[]) dataImage.getPixels();
+		float[] pixels = new float[newSize * newSize];
+		// Note that the limits at 0 and size-1 the taper is zero so this can be ignored
+		final int maxy_1 = dataImage.getHeight() - 1;
+		final int maxx_1 = dataImage.getWidth() - 1;
+		final int oldWidth = dataImage.getWidth();
+		for (int y = 1; y < maxy_1; y++)
+		{
+			final float yTmp = taperY[y];
+			for (int x = 1, i = y * oldWidth + 1, ii = y * newSize + 1; x < maxx_1; x++, i++, ii++)
+			{
+				pixels[ii] = data[i] * taperX[x] * yTmp;
+			}
+		}
+
+		return new FloatProcessor(newSize, newSize, pixels, null);
+	}
+
+	private float[] getWindowFunction(float[] taper, int size)
+	{
+		if (taper.length != size)
+		{
+			// Re-use cached values
+			if (taperX.length == size)
+				return taperX;
+			if (taperY.length == size)
+				return taperY;
+
+			final int boundary = size / 8;
+			final int upperBoundary = size - boundary;
+			taper = new float[size];
+			for (int i = 0; i < size; i++)
+			{
+				if ((i < boundary) || (i > size - upperBoundary))
+				{
+					taper[i] = (float) Math.pow(Math.sin(12.566370614359172D * i / size), 2.0D);
+				}
+				else
+				{
+					taper[i] = 1;
+				}
+			}
+		}
+		return taper;
+	}
+
+	/**
+	 * Adapted from ij.process.ImageProcessor.getInterpolatedValue(int,int).
+	 * <p>
+	 * Removed bounds checking and compute multiple values at the same time for multiple images.
+	 * 
+	 * @param x
+	 * @param y
+	 * @return
+	 */
+	private double[] getInterpolatedValues(final double x, final double y, float[][] images, final int maxx)
+	{
+		final int xbase = (int) x;
+		final int ybase = (int) y;
+		double xFraction = x - xbase;
+		double yFraction = y - ybase;
+		if (xFraction < 0.0)
+			xFraction = 0.0;
+		if (yFraction < 0.0)
+			yFraction = 0.0;
+
+		final int lowerLeftIndex = ybase * maxx + xbase;
+		final int lowerRightIndex = lowerLeftIndex + 1;
+		final int upperLeftIndex = lowerLeftIndex + maxx;
+		final int upperRightIndex = upperLeftIndex + 1;
+
+		final int noImages = 3; //images.length;
+		double[] values = new double[noImages];
+		for (int i = 0; i < noImages; i++)
+		{
+			final float[] image = images[i];
+			final double lowerLeft = image[lowerLeftIndex];
+			final double lowerRight = image[lowerRightIndex];
+			final double upperRight = image[upperLeftIndex];
+			final double upperLeft = image[upperRightIndex];
+
+			final double upperAverage = upperLeft + xFraction * (upperRight - upperLeft);
+			final double lowerAverage = lowerLeft + xFraction * (lowerRight - lowerLeft);
+			values[i] = lowerAverage + yFraction * (upperAverage - lowerAverage);
+		}
+		return values;
+	}
+
+	/**
+	 * Perform LOESS smoothing on the FRC curve.
+	 * <p>
+	 * The input curve is copied and then the correlation values are smoothed using a LOESS interpolation with the given
+	 * parameters. If smoothing fails the original curve values are returned.
+	 * 
+	 * @param frcCurve
+	 * @param bandwidth
+	 * @param robustness
+	 * @return A new FRC curve
+	 */
+	public double[][] getSmoothedCurve(double[][] frcCurve, double bandwidth, int robustness)
+	{
+		double[][] sCurve = new double[frcCurve.length][3];
+
+		double[] xVals = new double[frcCurve.length];
+		double[] yVals = new double[frcCurve.length];
+
+		for (int i = 0; i < frcCurve.length; i++)
+		{
+			xVals[i] = frcCurve[i][0];
+			yVals[i] = frcCurve[i][1];
+
+			sCurve[i][0] = frcCurve[i][0];
+			sCurve[i][1] = frcCurve[i][1];
+			sCurve[i][2] = frcCurve[i][2];
+		}
+
+		double[] ySmoothed = new double[frcCurve.length];
+
+		try
+		{
+			LoessInterpolator loess = new LoessInterpolator(bandwidth, robustness);
+			ySmoothed = loess.smooth(xVals, yVals);
+
+			for (int i = 0; i < frcCurve.length; i++)
+				sCurve[i][1] = ySmoothed[i];
+		}
+		catch (Exception e)
+		{
+			e.printStackTrace();
+		}
+
+		return sCurve;
+	}
+
+	/**
+	 * Perform LOESS smoothing on the FRC curve
+	 * <p>
+	 * The input curve is copied and then the correlation values are smoothed using a LOESS interpolation with a
+	 * bandwidth of 0.0707 and robustness of 0. If smoothing fails the original curve values are returned.
+	 * 
+	 * @param frcCurve
+	 * @return A new FRC curve
+	 */
+	public double[][] getSmoothedCurve(double[][] frcCurve)
+	{
+		double bandwidth = 0.0707;
+		int robustness = 0;
+		return getSmoothedCurve(frcCurve, bandwidth, robustness);
+	}
+
+	/**
+	 * Calculate the curve representing the minimum correlation required to distinguish two images for each resolution
+	 * in the input FRC curve.
+	 * 
+	 * @param frcCurve
+	 * @param method
+	 * @return The threshold curve representing the threshold for each input spatial frequency
+	 */
+	public double[] calculateThresholdCurve(double[][] frcCurve, ThresholdMethod method)
+	{
+		double[] threshold = new double[frcCurve.length];
+
+		for (int i = 0; i < threshold.length; i++)
+		{
+			switch (method)
+			{
+				case HALF_BIT:
+					threshold[i] = ((0.2071 * Math.sqrt(frcCurve[i][2]) + 1.9102) / (1.2071 * Math.sqrt(frcCurve[i][2]) + 0.9102));
+					break;
+				case THREE_SIGMA:
+					threshold[i] = (3.0 / Math.sqrt(frcCurve[i][2] / 2.0));
+					break;
+				case FIXED_1_OVER_7:
+				default:
+					threshold[i] = 0.1428;
+			}
+		}
+
+		return threshold;
+	}
+
+	/**
+	 * Computes the crossing points of the FRC curve and the threshold curve. The intersections can be used to determine
+	 * the image resolution using {@link #getCorrectIntersection(ArrayList, ThresholdMethod)}
+	 * 
+	 * @param frcCurve
+	 * @param thresholdCurve
+	 * @return The crossing points
+	 */
+	public double[] getIntersections(double[][] frcCurve, double[] thresholdCurve)
+	{
+		if (frcCurve.length != thresholdCurve.length)
+		{
+			IJ.error("Error", "Unable to calculate FRC curve intersections due to input length mismatch.");
+			return null;
+		}
+
+		double[] intersections = new double[frcCurve.length - 1];
+		int count = 0;
+
+		for (int i = 1; i < frcCurve.length; i++)
+		{
+			// http://en.wikipedia.org/wiki/Line-line_intersection
+			//
+			//     x1,y1            x4,y4      
+			//         **        ++ 
+			//           **    ++
+			//             **++ P(x,y)
+			//            ++ **
+			//          ++     **
+			//        ++         **
+			//    x3,y3            ** 
+			//                       x2,y2  
+
+			final double y1 = frcCurve[i - 1][1];
+			final double y2 = frcCurve[i][1];
+			final double y3 = thresholdCurve[i - 1];
+			final double y4 = thresholdCurve[i];
+
+			// Check if they cross
+			if (!((y3 >= y1 && y4 < y2) || (y1 >= y3 && y2 < y4)))
+			{
+				continue;
+			}			
+
+			final double x1 = frcCurve[i - 1][0];
+			final double x2 = frcCurve[i][0];
+			final double x3 = x1;
+			final double x4 = x2;
+			
+			final double x1_x2 = x1 - x2;
+			final double x3_x4 = x3 - x4;
+			final double y1_y2 = y1 - y2;
+			final double y3_y4 = y3 - y4;
+
+			// Check if lines are parallel
+			if (x1_x2 * y3_y4 - y1_y2 * x3_x4 == 0)
+			{
+				if (y1 == y3)
+					// The lines are the same
+					intersections[count++] = x1;
+			}
+			else
+			{
+				// Find intersection
+				double px = ((x1 * y2 - y1 * x2) * x3_x4 - x1_x2 * (x3 * y4 - y3 * x4)) /
+						(x1_x2 * y3_y4 - y1_y2 * x3_x4);
+				//double px = ((x1 * y2 - y1 * x2) * (x3 - x4) - (x1 - x2) * (x3 * y4 - y3 * x4)) /
+				//		((x1 - x2) * (y3 - y4) - (y1 - y2) * (x3 - x4));
+				//double py = ((x1 * y2 - y1 * x2) * (y3 - y4) - (y1 - y2) * (x3 * y4 - y3 * x4)) /
+				//		((x1 - x2) * (y3 - y4) - (y1 - y2) * (x3 - x4));
+
+				// Check if the intersection is within the two points
+				// Q. Is this necessary given the intersection check above?
+				if (px >= x1 && px < x2)
+				{
+					intersections[count++] = px;
+				}
+			}
+		}
+
+		return Arrays.copyOf(intersections, count);
+	}
+
+	/**
+	 * Get the correction intersection representing the image resolution. The intersection chosen depends on the method
+	 * used to calculate the threshold curve using {@link #calculateThresholdCurve(double[][], ThresholdMethod)}
+	 * <p>
+	 * The intersection corresponds the lowest spatial frequency at which there is no significant correlation between
+	 * the images.
+	 * 
+	 * @param intersections
+	 * @param method
+	 * @return The intersection (or zero if no crossings)
+	 */
+	public double getCorrectIntersection(double[] intersections, ThresholdMethod method)
+	{
+		if (intersections == null || intersections.length == 0)
+			return 0;
+
+		switch (method)
+		{
+		// The half-bit and 3-sigma curves are often above 1 at close to zero spatial frequency.
+		// This means that any FRC curve starting at 1 may cross the line twice. 
+		// If so the second crossing is the one that is desired.
+			case HALF_BIT:
+			case THREE_SIGMA:
+				return (intersections.length > 1) ? intersections[1] : intersections[0];
+
+			case FIXED_1_OVER_7:
+			default:
+				return intersections[0];
+		}
+	}
+
+	/**
+	 * Utility function that calculates the Fourier Image Resolution (FIRE) number using the provided images.
+	 * 
+	 * @param ip1
+	 * @param ip2
+	 * @param method
+	 * @return The FIRE number (in pixels)
+	 */
+	public double calculateFireNumber(ImageProcessor ip1, ImageProcessor ip2, ThresholdMethod method)
+	{
+		double[][] frcCurve = calculateFrcCurve(ip1, ip2);
+		return calculateFireNumber(frcCurve, method);
+	}
+
+	/**
+	 * Utility function that calculates the Fourier Image Resolution (FIRE) number using the provided FRC curve data.
+	 * 
+	 * @param frcCurve
+	 * @param method
+	 * @return The FIRE number (in pixels)
+	 */
+	public double calculateFireNumber(double[][] frcCurve, ThresholdMethod method)
+	{
+		double[] thresholdCurve = calculateThresholdCurve(frcCurve, method);
+		double[] intersections = getIntersections(frcCurve, thresholdCurve);
+
+		double fire = Double.NaN;
+		if (intersections == null || intersections.length != 0)
+		{
+			double spatialFrequency = getCorrectIntersection(intersections, method);
+			// Since the Fourier calculation only uses half of the image (from centre to the edge) 
+			// we must double the curve length to get the original maximum image width. In addition
+			// the computation was up to the edge-1 pixels so add back a pixel to the curve length.
+			fire = 2 * (frcCurve.length + 1) / spatialFrequency;
+		}
+		return fire;
+	}
+	
+
+
+	/**
+	 * Convenience function to batch calculate the FRC by going though file on two separate directories. 
+	 * When one file name is identical in both directories, it opens them and computes the FRC.
+	 * If is_save_plot is true, then a plot of the FRC and FIRE calculations is saved in the parent of 
+	 * directory1 in a folder named "Graphs"
+	 * NOTE: This only works for TIFF files at the moment.
+	 * 
+	 * @param directory1 is the first directory containing data
+	 * @param directory2 is the first directory containing data
+	 * @param method is the {@link ThresholdMethod} used
+	 * @param rt is the {@link ij.measure.ResultsTable} where you want the batch to write the FIRE number to
+	 * @param is_save_plot defines whether the plot will be saved, as described above.
+	 */
+	public void batchCalculateFireNumber(File directory1, File directory2, ThresholdMethod method, ResultsTable rt, boolean is_save_plot)
+	{
+		// Navigate folder for tiffs
+		String[] the_files = directory1.list(new FilenameFilter() {
+			
+			@Override
+			public boolean accept(File dir, String name) {
+				// check extensions
+				return true;				
+			}
+		});
+		
+		// Prepare saving if necessary
+		File save_dir = new File(directory1.getParentFile(), "Graphs");
+		save_dir.mkdir();
+		
+		// For each file, open one in each directory
+		for(String the_file : the_files) {
+			File f1 = new File(directory1.getAbsolutePath()+File.separator+the_file);
+			File f2 = new File(directory2.getAbsolutePath()+File.separator+the_file);
+			
+			if(f2.exists()) {
+				ImagePlus i1 = IJ.openImage(f1.getAbsolutePath());
+				ImagePlus i2 = IJ.openImage(f2.getAbsolutePath());
+				
+				// Finally calculate FRC
+				double[][] frc_curve = calculateFrcCurve(i1.getProcessor(), i2.getProcessor());
+				
+				double[][] smooth_frc = getSmoothedCurve(frc_curve);
+
+				// Fourier Image REsolution number ("FIRE")
+				double fire = calculateFireNumber(smooth_frc, method);
+				
+				rt.incrementCounter();
+				rt.addLabel(i1.getTitle());
+				rt.addValue("FRC ["+method+"]", fire);
+				rt.addValue("FRC ["+method+"] Calibrated ["+i1.getCalibration().getUnit()+"]", fire*i1.getCalibration().pixelHeight);
+				rt.show("FRC Results");
+				
+				if(is_save_plot) {
+					Plot p = doPlot(frc_curve, smooth_frc, method, fire, i1.getTitle());
+					
+					ImagePlus plot_image = p.makeHighResolution("FRC", 3, true, false);
+					String plot_name = save_dir.getAbsolutePath()+File.separator+f1.getName().substring(0,f1.getName().lastIndexOf("."));
+					plot_name += "_"+method.toString().replaceAll("/", " over ")+".tif";
+					IJ.save(plot_image, plot_name);
+				}
+			}
+			
+		}
+		
+	}
+	
+	/**
+	 * Simple method to plot the data using {@link ij.gui.Plot #Plot}
+	 * @param frc_curve is the raw FRC curve
+	 * @param smooth_frc is the smoothed curve computed using {@link #getSmoothedCurve(double[][])} or {@link #getSmoothedCurve(double[][], double, int)}
+	 * @param tm is the {@link ThresholdMethod} used
+	 * @param fire is the Fourier Image Resolution given by {@link #calculateFireNumber(double[][], ThresholdMethod)}
+	 * @param name is the title of the image that was used for the FRC calculations
+	 * @return a Plot object to use at your convenience
+	 */
+	public Plot doPlot(double[][] frc_curve, double[][] smooth_frc, ThresholdMethod tm, double fire, String name) {
+		
+		FRCbiop frc = new FRCbiop();
+		
+		
+		// Prepare arrays for Plot class
+		double[] x  = new double[frc_curve.length];
+		double[] y  = new double[frc_curve.length];
+		double[] sy = new double[frc_curve.length];
+		
+		
+		for(int i=0; i<frc_curve.length; i++) {
+			// Since the Fourier calculation only uses half of the image (from centre to the edge) 
+			// we must double the curve length to get the original maximum image width. In addition
+			// the computation was up to the edge-1 pixels so add back a pixel to the curve length.
+			// If we divide the value of the x axes by the highest spatial frequency (representing 1 pixel^-1)
+			// we can get a calibrated frequency axis.
+			x[i]  =  frc_curve[i][0]/(2*(frc_curve.length+1));
+			
+			// Original FRC curve
+			y[i]  =  frc_curve[i][1];
+			
+			// Smoothed FRC Curve
+			sy[i] = smooth_frc[i][1];
+		}
+		// Get Curve of Threshold
+		
+		// Curve representing the Threshold method calculation that should intercept with the FRC Curve
+		double[] thr_curve = frc.calculateThresholdCurve(smooth_frc, tm);
+		
+		// Plot the data
+		Plot p = new Plot("FRC Of "+name, "Spatial Frequency", "Correlation");
+
+		p.setLineWidth(1);
+
+		// Set Limits
+		p.setLimits(0, x[x.length-1], 0, 1);
+		
+		// Add Original Data in black
+		p.setColor(new Color(0, 0, 0));
+		p.addPoints(x, y, PlotWindow.LINE);
+		
+		// Add Smoothed Curve in clear red
+		p.setColor(new Color(255, 120, 120));
+		p.addPoints(x, sy, PlotWindow.LINE);
+
+		// Add Threshold Curve in clear blue
+		p.setColor(new Color(120, 120, 255));
+		p.addPoints(x, thr_curve, PlotWindow.LINE);
+
+		// Add FIRE number line in gray
+		p.setColor(new Color(69, 69, 69));
+		p.drawLine(1/fire, 0, 1/fire, 1);
+		
+		// Add the FIRE number on the graph
+		p.addLabel(0.02, 0.3, "FIRE = "+String.format("%.3f", fire));
+		
+		// Add legend to the plot
+		p.addLegend("FRC\nSmoothed FRC\nThreshold");
+		
+		return p;
+	}
+        
+        
+        
+        
+        
+    
+    
+    
+    public void plot(double [] y,String title){
+        
+        double yMin=y[0]; double yMax=y[0];
+        double [] x = new double [y.length];
+        for (int i=0;i<y.length;i++){
+            x[i]=i;
+            if (y[i]<yMin){
+                yMin=y[i];
+            }
+            if (y[i]>yMax){
+                yMax=y[i];
+            }
+            
+        }
+        
+        
+        Plot p = new Plot(""+title,"","");
+        p.setLimits(0, x.length, yMin-.01*(yMax-yMin), yMax+.01*(yMax-yMin));
+        p.setLineWidth(2);
+        p.setFont(0, 18);
+        
+        p.addPoints(x, y,Plot.LINE);
+        p.show();
+        
+    }
+    
+	
+}
diff --git a/src/main/java/org/pasteur/imagej/utils/FSC.java b/src/main/java/org/pasteur/imagej/utils/FSC.java
old mode 100644
new mode 100755
index 312e07d..1cf3f4e
--- a/src/main/java/org/pasteur/imagej/utils/FSC.java
+++ b/src/main/java/org/pasteur/imagej/utils/FSC.java
@@ -145,7 +145,7 @@ public void run(){
                     else{
                         ztmp=1;
                     }
-                    weight[iii][i][ii]=1;//xtmp*ytmp*ztmp;
+                    weight[iii][i][ii]=xtmp*ytmp*ztmp;
                 }
             }
         }
@@ -163,6 +163,7 @@ public void run(){
         double [][][] F = new double [size][size][size];
         double [][][] G = new double [size][size][size];
         double [][][] H = new double [size][size][size];
+        double [][][] numeratorI = new double [size][size][size];
         double [][][] tmpointer;
         FastFourierTransform fft1 = new FastFourierTransform(size,size,size);
         FastFourierTransform fft2 = new FastFourierTransform(size,size,size);
@@ -235,6 +236,7 @@ public void run(){
                     F[iii][i][ii]=(C[iii][i][ii]*C[iii][i][ii]+D[iii][i][ii]*D[iii][i][ii]);
                     G[iii][i][ii]=(A[iii][i][ii]*C[iii][i][ii]+B[iii][i][ii]*D[iii][i][ii]);
                     H[iii][i][ii]=(B[iii][i][ii]*C[iii][i][ii]-A[iii][i][ii]*D[iii][i][ii]);
+                    numeratorI[iii][i][ii]=(G[iii][i][ii]*G[iii][i][ii]+H[iii][i][ii]*H[iii][i][ii]);
                 }
             }
         }
@@ -259,12 +261,13 @@ public void run(){
         
         
         
-        double step=1;//pixel step to convert in vector
+        double step=2;//pixel step to convert in vector
         int sizeVect=(int)Math.ceil(size*Math.sqrt(3)/(2*step));//sqrt(3) because 3D
         double [] vectE = new double[sizeVect];
         double [] vectF = new double[sizeVect];
         double [] vectG = new double[sizeVect];
         double [] vectH = new double[sizeVect];
+        double [] numerator = new double[sizeVect];
         double [] vectN = new double[sizeVect];//number of pixels
         
         
@@ -278,7 +281,7 @@ public void run(){
             vectG[bin]+=G[(int)v[i][1]][(int)v[i][2]][(int)v[i][3]];
             vectH[bin]+=H[(int)v[i][1]][(int)v[i][2]][(int)v[i][3]];
             vectN[bin]++;
-            
+            numerator[bin]+=G[(int)v[i][1]][(int)v[i][2]][(int)v[i][3]]*G[(int)v[i][1]][(int)v[i][2]][(int)v[i][3]]+H[(int)v[i][1]][(int)v[i][2]][(int)v[i][3]]*H[(int)v[i][1]][(int)v[i][2]][(int)v[i][3]];
         }
         
         
@@ -289,9 +292,9 @@ public void run(){
         double [] xaxisXY = new double[sizeVect];
         double [] xaxisZ = new double[sizeVect];
         for (int i=0;i<sizeVect;i++){
-            vectFinal[i]=(vectG[i]/Math.sqrt(vectE[i]*vectF[i]));// - 2/Math.sqrt(vectN[i]/2);
-            xaxisXY[i]=i*step/(size*this.pixSizeXY);
-            xaxisZ[i]=i*step/(size*this.pixSizeZ);
+            vectFinal[i]=(Math.sqrt(numerator[i])/Math.sqrt(vectE[i]*vectF[i]));// - 2/Math.sqrt(vectN[i]/2);
+            xaxisXY[i]=i*step/((size)*this.pixSizeXY);//(size/2) because we use half of image ?
+            xaxisZ[i]=i*step/((size)*this.pixSizeZ);
             
         }
         
diff --git a/src/main/java/org/pasteur/imagej/utils/FastFourierTransform.java b/src/main/java/org/pasteur/imagej/utils/FastFourierTransform.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/utils/FileVectorLoader.java b/src/main/java/org/pasteur/imagej/utils/FileVectorLoader.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/utils/FourierTransform.java b/src/main/java/org/pasteur/imagej/utils/FourierTransform.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/utils/GaussianElimination.java b/src/main/java/org/pasteur/imagej/utils/GaussianElimination.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/utils/ImageShow.java b/src/main/java/org/pasteur/imagej/utils/ImageShow.java
old mode 100644
new mode 100755
index 8aadf37..f1ec682
--- a/src/main/java/org/pasteur/imagej/utils/ImageShow.java
+++ b/src/main/java/org/pasteur/imagej/utils/ImageShow.java
@@ -329,4 +329,70 @@ public static void imshow(float [][] im,String title){
         imp.show();
     }
     
+    
+    
+    
+    
+    
+    
+    
+    
+    
+    public static void imshow(double [] im,int width,int height){
+        imshow(im, width, height, "");
+    }
+    public static void imshow(double [] im,int width,int height,String title){
+        
+            ImageProcessor ip = new FloatProcessor(width,height);
+            for (int i=0;i<width;i++){
+                for (int ii=0;ii<height;ii++){
+                    ip.putPixelValue(i, ii, (float)im[ii*width+i]);
+                }
+            }
+        
+        ImagePlus imp=new ImagePlus(""+title,ip);
+        imp.show();
+    }
+    
+    
+    
+    public static void imshow(float [] im,int width,int height){
+        imshow(im, width, height, "");
+    }
+    public static void imshow(float [] im,int width,int height,String title){
+        
+            ImageProcessor ip = new FloatProcessor(width,height);
+            for (int i=0;i<width;i++){
+                for (int ii=0;ii<height;ii++){
+                    ip.putPixelValue(i, ii, (float)im[ii*width+i]);
+                }
+            }
+        
+        ImagePlus imp=new ImagePlus(""+title,ip);
+        imp.show();
+    }
+    
+    public static void imshow(int [] im,int width,int height){
+        imshow(im, width, height, "");
+    }
+    public static void imshow(int [] im,int width,int height,String title){
+        
+            ImageProcessor ip = new FloatProcessor(width,height);
+            for (int i=0;i<width;i++){
+                for (int ii=0;ii<height;ii++){
+                    ip.putPixelValue(i, ii, (float)im[ii*width+i]);
+                }
+            }
+        
+        ImagePlus imp=new ImagePlus(""+title,ip);
+        imp.show();
+    }
+    
+    
+    
+    
+    
+    
+    
+    
 }
diff --git a/src/main/java/org/pasteur/imagej/utils/JSONformat.java b/src/main/java/org/pasteur/imagej/utils/JSONformat.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/utils/Matrixe.java b/src/main/java/org/pasteur/imagej/utils/Matrixe.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/utils/PolynomialFit.java b/src/main/java/org/pasteur/imagej/utils/PolynomialFit.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/utils/Protocol_localization.java b/src/main/java/org/pasteur/imagej/utils/Protocol_localization.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/utils/Regression.java b/src/main/java/org/pasteur/imagej/utils/Regression.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/utils/SCMOScamera.java b/src/main/java/org/pasteur/imagej/utils/SCMOScamera.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/utils/SplineFit.java b/src/main/java/org/pasteur/imagej/utils/SplineFit.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/utils/Statistic.java b/src/main/java/org/pasteur/imagej/utils/Statistic.java
old mode 100644
new mode 100755
diff --git a/src/main/java/org/pasteur/imagej/utils/Zernike.java b/src/main/java/org/pasteur/imagej/utils/Zernike.java
old mode 100644
new mode 100755
diff --git a/src/main/resources/kernels/ZOLA_CudaCode_ReadMe_file.txt b/src/main/resources/kernels/ZOLA_CudaCode_ReadMe_file.txt
old mode 100644
new mode 100755
diff --git a/src/main/resources/kernels/ZOLA_kernels_30.ptx b/src/main/resources/kernels/ZOLA_kernels_30.ptx
old mode 100644
new mode 100755
diff --git a/src/main/resources/kernels/ZOLA_kernels_60.ptx b/src/main/resources/kernels/ZOLA_kernels_60.ptx
old mode 100644
new mode 100755
diff --git a/src/main/resources/kernels/ZOLA_kernels_61.ptx b/src/main/resources/kernels/ZOLA_kernels_61.ptx
old mode 100644
new mode 100755
diff --git a/src/main/resources/plugins.config b/src/main/resources/plugins.config
old mode 100644
new mode 100755