diff --git a/src/main/java/StartingOSStats.java b/src/main/java/StartingOSStats.java
index 97dcd05ad18109a56ba87db87ed532b66bd292f5..9c0c4c3fe8c0765ec3df25a9884e4b1fffe8879b 100644
--- a/src/main/java/StartingOSStats.java
+++ b/src/main/java/StartingOSStats.java
@@ -100,7 +100,7 @@ public class StartingOSStats
             ImageJFunctions.show( TestMIP.findMIP( stackImage, stackImage.factory() ), "MIP" );
 
 
-            SurfacesExtraction.extract(stackImage,amplitude, otsu, sigmas, false, delta, false,
+            SurfacesExtraction.extract(stackImage, stackImage.factory(),amplitude, otsu, sigmas, false, delta, false,
                     false, true, "MIP", k1, percent1, k2, percent2);
         }
         else
diff --git a/src/main/java/SurfacesExtractionAnalyse.java b/src/main/java/SurfacesExtractionAnalyse.java
index bdf8e2f20360acde36786b057d0385737ec189dc..ed2a32fd47b36d1083a6653e2e3c833d44c645b9 100644
--- a/src/main/java/SurfacesExtractionAnalyse.java
+++ b/src/main/java/SurfacesExtractionAnalyse.java
@@ -27,375 +27,349 @@ import java.util.ArrayList;
 
 public class SurfacesExtractionAnalyse
 {
-    public static < T extends RealType < T > & NativeType < T > > void main( String[] args )
-    {
-
-        /* User Parameters AKA arguments to run the program*/
-        double amplitude = Double.parseDouble( args[ 0 ] );
-        double otsu = Double.parseDouble( args[ 1 ] );
-        int sigmas = Integer.parseInt( args[ 2 ] );
-        int k1 = Integer.parseInt( args[ 3 ] );
-        double percent1 = Double.parseDouble( args[ 4 ] );
-        int k2 = Integer.parseInt( args[ 5 ] );
-        double percent2 = Double.parseDouble( args[ 6 ] );
-        int delta = Integer.parseInt( args[ 7 ] );
-
-        final String imagePath = args[8].trim();
+//    public static < T extends RealType < T > & NativeType < T > > void main( String[] args )
+//    {
 //
-        System.out.println(imagePath);
-        // Creation of the image : version with unsigned type. */
-        /* JY version for opening files. */
-        final SCIFIOImgPlus< ? > imgPlus = IO.openImgs( imagePath ).get( 0 );
-
-        final Img < T > stackImage = ( Img < T > ) imgPlus.getImg();
-        /* End of parameters. */
-
-
-
-    }
-
-    public static < T extends RealType< T > & NativeType< T > > void extract(
-            Img< T > input, double amplitude, double threshold, int sigmas,
-            boolean zMapDisplay, int delta,
-            boolean extractedStackDisplay,
-            boolean heightMapDisplay,
-            boolean projectionDisplay,
-            String projectionType ,
-            int k1, double percent1, int k2, double percent2)
-    {
-        /* First step : Pixel selection */
-        Pixels[][] maximums = SurfacePixelSelection.run( input, amplitude, threshold, sigmas );
-        int width = ( int ) input.dimension( 0 );
-        int height = ( int ) input.dimension( 1 );
-        try
-        {
-            /*  First round construction*/
-            ArrayList< Surface > surfaces = firstRoundConstruction( maximums, width, height , percent1, k1);
-            System.out.println( "first round surfaces = " + surfaces.size() );
-//
-            for (Surface surface : surfaces) {
-                displaySurface(surface);
-            }
-            // Second round construction
-            ArrayList< Surface > finalSurfaces = secondRoundConstruction( surfaces, width, height, percent2, k2 );
-//            if ( !finalSurfaces.isEmpty() )
-//            {
-            int index = 0;
-            for ( Surface surface : finalSurfaces )
-            {
-                System.out.println( "Surface number " + index );
+//        /* User Parameters AKA arguments to run the program*/
+//        double amplitude = Double.parseDouble( args[ 0 ] );
+//        double otsu = Double.parseDouble( args[ 1 ] );
+//        int sigmas = Integer.parseInt( args[ 2 ] );
+//        int k1 = Integer.parseInt( args[ 3 ] );
+//        double percent1 = Double.parseDouble( args[ 4 ] );
+//        int k2 = Integer.parseInt( args[ 5 ] );
+//        double percent2 = Double.parseDouble( args[ 6 ] );
+//        int delta = Integer.parseInt( args[ 7 ] );
+//
+//        final String imagePath = args[8].trim();
+////
+//        System.out.println(imagePath);
+//        // Creation of the image : version with unsigned type. */
+//        /* JY version for opening files. */
+//        final SCIFIOImgPlus< ? > imgPlus = IO.openImgs( imagePath ).get( 0 );
+//
+//        final Img < T > stackImage = ( Img < T > ) imgPlus.getImg();
+//        /* End of parameters. */
+//
+//
+//
+//    }
+//
+//    public static < T extends RealType< T > & NativeType< T > > void extract(
+//            Img< T > input, double amplitude, double threshold, int sigmas,
+//            boolean zMapDisplay, int delta,
+//            boolean extractedStackDisplay,
+//            boolean heightMapDisplay,
+//            boolean projectionDisplay,
+//            String projectionType ,
+//            int k1, double percent1, int k2, double percent2)
+//    {
+//        /* First step : Pixel selection */
+//        Pixels[][] maximums = SurfacePixelSelection.run( input, amplitude, threshold, sigmas );
+//        int width = ( int ) input.dimension( 0 );
+//        int height = ( int ) input.dimension( 1 );
+//        try
+//        {
+//            /*  First round construction*/
+//            ArrayList< Surface > surfaces = firstRoundConstruction( maximums, width, height , percent1, k1);
+//            System.out.println( "first round surfaces = " + surfaces.size() );
+////
+//            for (Surface surface : surfaces) {
 //                displaySurface(surface);
-                Img< UnsignedShortType > zMap = surface.getZMap();
-                projection( input, delta, zMap, zMapDisplay,
-                        extractedStackDisplay, heightMapDisplay,
-                        projectionDisplay, projectionType, index );
-                index++;
-            }
-
-//            // Display the projection
-            System.out.println( "The end" );
-        }
-        catch ( Exception e )
-        {
-            e.printStackTrace();
-        }
-    }
-
-
-    //TODO Bad documentation
-
-    /**
-     * This method is used to binarize a grey-level image using Otsu Thresholding method in order to classify its pixels
-     * into foreground/background and then apply a anisotropic gaussian blur.
-     */
-
-    private static < T extends RealType< T > & NativeType< T > > Img< T > maximumSearch( Img< T > input, double percent, int sigma )
-    {
-        /* Local maximum detection using partial derivative */
-        Img< T > maximums = LocalMaximumDetection.findMaximums( input, input.factory() );
-        ImageJFunctions.show( maximums, "maximums" );
-
-        /* Grid of local maximums*/
-        Img< BitType > threshold = LocalOtsuClassification.find( input, percent );
-        //ImageJFunctions.show( threshold, "Local threshold " );
-
-        /* Maxima selection according to threshold*/
-        maximums = Threshold.classification( maximums.copy(), maximums.factory(), threshold );
-
-        /* Dilatation of the resulting image*/
-        Utils.gaussConvolution( maximums.copy(), maximums, new double[]{ sigma, sigma } );
-//        ImageJFunctions.show( maximums, "blurred maximums" );
-
-        /* New local maximum detection due to previous smoothing*/
-        maximums = LocalMaximumDetection.findMaximums( maximums.copy(), maximums.factory() );
-        return maximums;
-    }
-
-
-    /**
-     * Setting of the unfiltered local maximums : first the original stack is filtered with an anisotropic
-     * gaussian convolution,  the pixel values are normalized between 0 and 255 then the local maximums
-     * are detected.
-     * During these steps, the temporary local, sectional and global thresholds are computed.
-     */
-
-
-    /**
-     * Creates the OS for each orthogonal section from the maximum coordinates found according to the dimension.
-     *
-     * @param maximums     - the local maximums found with the {@link LocalMaximumDetection} class.
-     * @param osListsArray - the OS corresponding to a 1D surface.
-     */
-    private static void set( Pixels[][] maximums, OSList[] osListsArray, int dimension, OSEStartingStatus startingStatus )
-    {
-
-
-        for ( int i = 0; i <= maximums.length - 1; i++ )
-        {
-            osListsArray[ i ] = OSConstruction.findOS( maximums[ i ], dimension, startingStatus );
-        }
-        startingStatus.setStartingStatus();
-    }
-
-    /* ----- First and second round construction methods. ----- */
-    private static ArrayList< Surface > firstRoundConstruction( Pixels[][] maximums, int width, int height, double percent, int k ) throws NoSurfaceFoundException
-    {
-        int dimension = 0;
-        // Setting and storage of the OS build in the X dimension according to the user thresholds settings.
-        OSEStartingStatus startingStatus = new OSEStartingStatus( dimension );
-        OSList[] osListsArrayX = new OSList[ height ];
-        set( maximums, osListsArrayX, dimension, startingStatus );
-//        displayOS(osListsArrayX);
-        System.out.println( "OS.count = " + OSE.getCount() );
-        /* Construction of the tempSurfaces*/
-        ArrayList< Surface > surfaces =
-                SurfacesReconstruction.buildSurfaces( 0, osListsArrayX, width, height, percent, k );
-        if ( ! surfaces.isEmpty() )
-        {
-            mergeReferenceSurface( surfaces );
-        }
-        else
-        {
-            System.out.println( "No Surface Exception" );
-            throw new FirstRoundConstructionException();
-        }
-        return surfaces;
-    }
-
-    /**
-     * Rebuilds the double pixel array from the {@link SurfaceLine} array
-     * of a {@link Surface} for a construction in the height dimension.
-     *
-     * @param surface - the {@link SurfaceLine}  array of a {@link Surface}
-     * @return a {@link Surface} specific double pixel array.
-     */
-    private static Pixels[][] rebuildPixelsArray( Surface surface, int width, int height )
-    {
-        Pixels[][] tempCoordinates = new Pixels[ width ][ height ]; // Transposed array
-        for ( int i = 0; i <= height - 1; i++ )
-        {
-            for ( int j = 0; j <= width - 1; j++ )
-            {
-                if ( surface.get( i ) != null )
-                {
-                    tempCoordinates[ j ][ i ] = surface.get( i ).get( j );
-                }
-            }
-        }
-        displayMaximums( tempCoordinates );
-        return tempCoordinates;
-    }
-
-    /**
-     * Converts and transposes the specified Surface object into the specified OSList array.
-     *
-     * @param surface       - the Surface object to process
-     * @param osListsArrayY -the output OSList array
-     */
-    private static void transposeSurfaceLines( Surface surface, OSList[] osListsArrayY, int width, int height )
-    {
-        int dimension = 1;
-        Pixels[][] pixels = rebuildPixelsArray( surface, width, height );
-        /* OS construction in dimension height.*/
-        OSEStartingStatus startingStatus = new OSEStartingStatus( dimension );
-        set( pixels, osListsArrayY, dimension, startingStatus );
-    }
-
-
-    /* -----  Other displaying methods -----*/
-
-    /**
-     * Reconstructs the TempSurface objects  of the specified list in dimension width.
-     *
-     * @param surfaces - the list of TempSurface objects.
-     * @return a list of reconstructed Tempsurface objects
-     */
-    private static ArrayList< Surface > secondRoundConstruction( ArrayList< Surface > surfaces, int width, int height, double percent, int k ) throws NoSurfaceFoundException
-    {
-        System.out.println( "========================================" );
-        ArrayList< Surface > finalSurfaces = new ArrayList<>();
-
-        for ( Surface surface : surfaces )
-        {
-            if ( surface.hasDuplicate() )// CoordinateList instead of Coordinate
-            {
-                OSList[] osListsArrayY = new OSList[ width ];
-
-                transposeSurfaceLines( surface, osListsArrayY, width, height );
-                OSE.setStartingStatus( 1 );
-                System.out.println( "OS.count = " + OSE.getCount() );
-
-//                displayOS( osListsArrayY );
-                ArrayList< Surface > temps = SurfacesReconstruction.buildSurfaces( 1, osListsArrayY, width, height, percent, k );
-                if ( temps.isEmpty() )
-                {
-                    System.out.println( "small surface in second round" );
-//                    throw new SecondRoundConstructionException();
-                }
-                else
-                {
-                    finalSurfaces.addAll( temps );
-                }
-            }
-            else
-            {
-                finalSurfaces.add( surface.transpose() );
-            }
-        }
-        mergeReferenceSurface( finalSurfaces );
-        return finalSurfaces;
-    }
-    //TODO Handle the NOSurfaceFoundException dynamically : this particular surface has to be constructed with a
-    // smaller StartingOSMinimumSize but not necessary the other ones.
-    // May be add a static variable to store the maximum OS size per surfaces
-    //
-
-    /**
-     * Merges the TempSurface objects of the specified list.
-     *
-     * @param surfaces - the list to merge
-     */
-    private static void mergeReferenceSurface( ArrayList< Surface > surfaces )
-    {
-        // Merging step.
-        ArrayList< Surface > toRemoved = new ArrayList<>();
-        for ( int i = surfaces.size() - 1; i > 0; i-- )
-        {
-            Surface one = surfaces.get( i );
-            for ( int j = i - 1; j >= 0; j-- )
-            {
-                Surface two = surfaces.get( j );
-                if ( one.sameSurface( two ) )
-                {
-                    two.merge( one );
-                    toRemoved.add( one );
-                }
-            }
-        }
-        surfaces.removeAll( toRemoved );
-    }
-
-
-    private static < R extends RealType< R > & NativeType< R > > void projection(
-            Img< R > original, int delta, Img< UnsignedShortType > zMap,
-            boolean zMapDisplay,
-            boolean extractedStackDisplay,
-            boolean heightMapDisplay,
-            boolean projectionDisplay,
-            String projectionType, int index )
-    {
-
-        // First step : fill in the holes. */
-        Img< UnsignedShortType > interpolated = Interpolation.execute( zMap );
-        Img< UnsignedShortType > smoothed = interpolated.copy();
-        Utils.gaussConvolution( interpolated, smoothed, new double[]{ 1.0, 1.0 } );
-//         Second step smoothed the elevation map
-        if ( zMapDisplay )
-        {
-            ImageJFunctions.show( smoothed, "Elevation map" + "RS n°" + index );
-        }
-
-        Img< R > extractedStack = StackProjection.getExtractedStack( original, smoothed, delta );
-        if ( extractedStackDisplay )
-        {
-            ImageJFunctions.show( extractedStack,
-                    "Extracted stack (delta = " + delta + ")" + "RS n°" + index );
-        }
-        if ( projectionDisplay )
-        {
-            // A heightmap can be displayed
-            if ( projectionType.equals( "MIP" ) ||
-                    projectionType.equals( "Minimum Intensity" ) )
-            {
-                Img< UnsignedShortType > heightMap = StackProjection.getHeightMap
-                        ( original, smoothed, projectionType, delta );
-                if ( heightMapDisplay )
-                {
-                    ImageJFunctions.show( heightMap, "Height Map with " +
-                            projectionType + " RS n°" + index );
-                }
-
-                Img< R > projection = StackProjection.projection1( original, heightMap, projectionType );
-                ImageJFunctions.show( projection, "Projection with " + projectionType + " RS n°" + index );
-            }
-            else
-            {
-
-                Img< R > projection = StackProjection.projection2( StackProjection.getHeightMapStack( original, zMap, delta ), projectionType );
-                ImageJFunctions.show( projection, "Projection with " + projectionType + " RS n°" + index );
-
-            }
-//            ImageJFunctions.show(StackProjection.getHeightMapStack(original, zMap, delta), "HMS");
-        }
-
-    }
-
-    /**
-     * Displays the local maximums found using jzy3D package.
-     */
-    public static void displayMaximums( Pixels[][] maximums )
-    {
-        try
-        {
-            LocalMaximumsDisplay localMaximumsDisplay = new LocalMaximumsDisplay( maximums );
-            AnalysisLauncher.open( localMaximumsDisplay );
-        }
-        catch ( Exception e )
-        {
-            e.printStackTrace();
-        }
-    }
-
-    /**
-     * Displays the local maximums found using jzy3D package.
-     */
-    public static void displayOS( OSList[] osLists )
-    {
-        try
-        {
-            OSDisplay display = new OSDisplay( osLists );
-            AnalysisLauncher.open( display );
-        }
-        catch ( Exception e )
-        {
-            e.printStackTrace();
-        }
-    }
-
-
-    public static void displaySurface( Surface surface )
-    {
-        try
-        {
-            SurfaceDisplay s = new SurfaceDisplay( surface );
-            AnalysisLauncher.open( s );
-        }
-        catch ( Exception e )
-        {
-            e.printStackTrace();
-        }
-
-    }
+//            }
+//            // Second round construction
+//            ArrayList< Surface > finalSurfaces = secondRoundConstruction( surfaces, width, height, percent2, k2 );
+////            if ( !finalSurfaces.isEmpty() )
+////            {
+//            int index = 0;
+//            for ( Surface surface : finalSurfaces )
+//            {
+//                System.out.println( "Surface number " + index );
+////                displaySurface(surface);
+//                Img< UnsignedShortType > zMap = surface.getZMap();
+//                projection( input, delta, zMap, zMapDisplay,
+//                        extractedStackDisplay, heightMapDisplay,
+//                        projectionDisplay, projectionType, index );
+//                index++;
+//            }
+//
+////            // Display the projection
+//            System.out.println( "The end" );
+//        }
+//        catch ( Exception e )
+//        {
+//            e.printStackTrace();
+//        }
+//    }
+//
+//
+//    //TODO Bad documentation
+//
+//
+//
+//    /**
+//     * Setting of the unfiltered local maximums : first the original stack is filtered with an anisotropic
+//     * gaussian convolution,  the pixel values are normalized between 0 and 255 then the local maximums
+//     * are detected.
+//     * During these steps, the temporary local, sectional and global thresholds are computed.
+//     */
+//
+//
+//    /**
+//     * Creates the OS for each orthogonal section from the maximum coordinates found according to the dimension.
+//     *
+//     * @param maximums     - the local maximums found with the {@link LocalMaximumDetection} class.
+//     * @param osListsArray - the OS corresponding to a 1D surface.
+//     */
+//    private static void set( Pixels[][] maximums, OSList[] osListsArray, int dimension, OSEStartingStatus startingStatus )
+//    {
+//
+//
+//        for ( int i = 0; i <= maximums.length - 1; i++ )
+//        {
+//            osListsArray[ i ] = OSConstruction.findOS( maximums[ i ], dimension, startingStatus );
+//        }
+//        startingStatus.setStartingStatus();
+//    }
+//
+//    /* ----- First and second round construction methods. ----- */
+//    private static ArrayList< Surface > firstRoundConstruction( Pixels[][] maximums, int width, int height, double percent, int k ) throws NoSurfaceFoundException
+//    {
+//        int dimension = 0;
+//        // Setting and storage of the OS build in the X dimension according to the user thresholds settings.
+//        OSEStartingStatus startingStatus = new OSEStartingStatus( dimension );
+//        OSList[] osListsArrayX = new OSList[ height ];
+//        set( maximums, osListsArrayX, dimension, startingStatus );
+////        displayOS(osListsArrayX);
+//        System.out.println( "OS.count = " + OSE.getCount() );
+//        /* Construction of the tempSurfaces*/
+//        ArrayList< Surface > surfaces =
+//                SurfacesReconstruction.buildSurfaces( 0, osListsArrayX, width, height, percent, k );
+//        if ( ! surfaces.isEmpty() )
+//        {
+//            mergeReferenceSurface( surfaces );
+//        }
+//        else
+//        {
+//            System.out.println( "No Surface Exception" );
+//            throw new FirstRoundConstructionException();
+//        }
+//        return surfaces;
+//    }
+//
+//    /**
+//     * Rebuilds the double pixel array from the {@link SurfaceLine} array
+//     * of a {@link Surface} for a construction in the height dimension.
+//     *
+//     * @param surface - the {@link SurfaceLine}  array of a {@link Surface}
+//     * @return a {@link Surface} specific double pixel array.
+//     */
+//    private static Pixels[][] rebuildPixelsArray( Surface surface, int width, int height )
+//    {
+//        Pixels[][] tempCoordinates = new Pixels[ width ][ height ]; // Transposed array
+//        for ( int i = 0; i <= height - 1; i++ )
+//        {
+//            for ( int j = 0; j <= width - 1; j++ )
+//            {
+//                if ( surface.get( i ) != null )
+//                {
+//                    tempCoordinates[ j ][ i ] = surface.get( i ).get( j );
+//                }
+//            }
+//        }
+//        displayMaximums( tempCoordinates );
+//        return tempCoordinates;
+//    }
+//
+//    /**
+//     * Converts and transposes the specified Surface object into the specified OSList array.
+//     *
+//     * @param surface       - the Surface object to process
+//     * @param osListsArrayY -the output OSList array
+//     */
+//    private static void transposeSurfaceLines( Surface surface, OSList[] osListsArrayY, int width, int height )
+//    {
+//        int dimension = 1;
+//        Pixels[][] pixels = rebuildPixelsArray( surface, width, height );
+//        /* OS construction in dimension height.*/
+//        OSEStartingStatus startingStatus = new OSEStartingStatus( dimension );
+//        set( pixels, osListsArrayY, dimension, startingStatus );
+//    }
+//
+//
+//    /* -----  Other displaying methods -----*/
+//
+//    /**
+//     * Reconstructs the TempSurface objects  of the specified list in dimension width.
+//     *
+//     * @param surfaces - the list of TempSurface objects.
+//     * @return a list of reconstructed Tempsurface objects
+//     */
+//    private static ArrayList< Surface > secondRoundConstruction( ArrayList< Surface > surfaces, int width, int height, double percent, int k ) throws NoSurfaceFoundException
+//    {
+//        System.out.println( "========================================" );
+//        ArrayList< Surface > finalSurfaces = new ArrayList<>();
+//
+//        for ( Surface surface : surfaces )
+//        {
+//            if ( surface.hasDuplicate() )// CoordinateList instead of Coordinate
+//            {
+//                OSList[] osListsArrayY = new OSList[ width ];
+//
+//                transposeSurfaceLines( surface, osListsArrayY, width, height );
+//                OSE.setStartingStatus( 1 );
+//                System.out.println( "OS.count = " + OSE.getCount() );
+//
+////                displayOS( osListsArrayY );
+//                ArrayList< Surface > temps = SurfacesReconstruction.buildSurfaces( 1, osListsArrayY, width, height, percent, k );
+//                if ( temps.isEmpty() )
+//                {
+//                    System.out.println( "small surface in second round" );
+////                    throw new SecondRoundConstructionException();
+//                }
+//                else
+//                {
+//                    finalSurfaces.addAll( temps );
+//                }
+//            }
+//            else
+//            {
+//                finalSurfaces.add( surface.transpose() );
+//            }
+//        }
+//        mergeReferenceSurface( finalSurfaces );
+//        return finalSurfaces;
+//    }
+//    //TODO Handle the NOSurfaceFoundException dynamically : this particular surface has to be constructed with a
+//    // smaller StartingOSMinimumSize but not necessary the other ones.
+//    // May be add a static variable to store the maximum OS size per surfaces
+//    //
+//
+//    /**
+//     * Merges the TempSurface objects of the specified list.
+//     *
+//     * @param surfaces - the list to merge
+//     */
+//    private static void mergeReferenceSurface( ArrayList< Surface > surfaces )
+//    {
+//        // Merging step.
+//        ArrayList< Surface > toRemoved = new ArrayList<>();
+//        for ( int i = surfaces.size() - 1; i > 0; i-- )
+//        {
+//            Surface one = surfaces.get( i );
+//            for ( int j = i - 1; j >= 0; j-- )
+//            {
+//                Surface two = surfaces.get( j );
+//                if ( one.sameSurface( two ) )
+//                {
+//                    two.merge( one );
+//                    toRemoved.add( one );
+//                }
+//            }
+//        }
+//        surfaces.removeAll( toRemoved );
+//    }
+//
+//
+//    private static < R extends RealType< R > & NativeType< R > > void projection(
+//            Img< R > original, int delta, Img< UnsignedShortType > zMap,
+//            boolean zMapDisplay,
+//            boolean extractedStackDisplay,
+//            boolean heightMapDisplay,
+//            boolean projectionDisplay,
+//            String projectionType, int index )
+//    {
+//
+//        // First step : fill in the holes. */
+//        Img< UnsignedShortType > interpolated = Interpolation.execute( zMap );
+//        Img< UnsignedShortType > smoothed = interpolated.copy();
+//        Utils.gaussConvolution( interpolated, smoothed, new double[]{ 1.0, 1.0 } );
+////         Second step smoothed the elevation map
+//        if ( zMapDisplay )
+//        {
+//            ImageJFunctions.show( smoothed, "Elevation map" + "RS n°" + index );
+//        }
+//
+//        Img< R > extractedStack = StackProjection.getExtractedStack( original, smoothed, delta );
+//        if ( extractedStackDisplay )
+//        {
+//            ImageJFunctions.show( extractedStack,
+//                    "Extracted stack (delta = " + delta + ")" + "RS n°" + index );
+//        }
+//        if ( projectionDisplay )
+//        {
+//            // A heightmap can be displayed
+//            if ( projectionType.equals( "MIP" ) ||
+//                    projectionType.equals( "Minimum Intensity" ) )
+//            {
+//                Img< UnsignedShortType > heightMap = StackProjection.getHeightMap
+//                        ( original, smoothed, projectionType, delta );
+//                if ( heightMapDisplay )
+//                {
+//                    ImageJFunctions.show( heightMap, "Height Map with " +
+//                            projectionType + " RS n°" + index );
+//                }
+//
+//                Img< R > projection = StackProjection.projection1( original, heightMap, projectionType );
+//                ImageJFunctions.show( projection, "Projection with " + projectionType + " RS n°" + index );
+//            }
+//            else
+//            {
+//
+//                Img< R > projection = StackProjection.projection2( StackProjection.getHeightMapStack( original, zMap, delta ), projectionType );
+//                ImageJFunctions.show( projection, "Projection with " + projectionType + " RS n°" + index );
+//
+//            }
+////            ImageJFunctions.show(StackProjection.getHeightMapStack(original, zMap, delta), "HMS");
+//        }
+//
+//    }
+//
+//    /**
+//     * Displays the local maximums found using jzy3D package.
+//     */
+//    public static void displayMaximums( Pixels[][] maximums )
+//    {
+//        try
+//        {
+//            LocalMaximumsDisplay localMaximumsDisplay = new LocalMaximumsDisplay( maximums );
+//            AnalysisLauncher.open( localMaximumsDisplay );
+//        }
+//        catch ( Exception e )
+//        {
+//            e.printStackTrace();
+//        }
+//    }
+//
+//    /**
+//     * Displays the local maximums found using jzy3D package.
+//     */
+//    public static void displayOS( OSList[] osLists )
+//    {
+//        try
+//        {
+//            OSDisplay display = new OSDisplay( osLists );
+//            AnalysisLauncher.open( display );
+//        }
+//        catch ( Exception e )
+//        {
+//            e.printStackTrace();
+//        }
+//    }
+//
+//
+//    public static void displaySurface( Surface surface )
+//    {
+//        try
+//        {
+//            SurfaceDisplay s = new SurfaceDisplay( surface );
+//            AnalysisLauncher.open( s );
+//        }
+//        catch ( Exception e )
+//        {
+//            e.printStackTrace();
+//        }
+//
+//    }
 
 }
 
diff --git a/src/main/java/fr/pasteur/ida/zellige/main/Main.java b/src/main/java/fr/pasteur/ida/zellige/main/Main.java
index 8f01c10bf21870099e8539099e878f5f7a1bc159..f014ae61912712e9f1ce429887bdc6982da94899 100644
--- a/src/main/java/fr/pasteur/ida/zellige/main/Main.java
+++ b/src/main/java/fr/pasteur/ida/zellige/main/Main.java
@@ -102,7 +102,7 @@ public class Main
             ImageJFunctions.show( stackImage, "original" );
             ImageJFunctions.show( TestMIP.findMIP( stackImage, stackImage.factory() ), "MIP" );
 
-            SurfacesExtraction.extract(stackImage,amplitude,otsu, sigmas, false, delta, false,
+            SurfacesExtraction.extract(stackImage, stackImage.factory(),amplitude,otsu, sigmas, false, delta, false,
                     false, true, "MIP", k1, percent1, k2, percent2);
         }
         else
diff --git a/src/main/java/fr/pasteur/ida/zellige/surfaceConstruction/construction/SurfacesExtraction.java b/src/main/java/fr/pasteur/ida/zellige/surfaceConstruction/construction/SurfacesExtraction.java
index 688c16de66d026e3537328702abcc8aada6e64c1..e62720161e0951a3b7ee6fe2f08fd8bd731dbf0b 100644
--- a/src/main/java/fr/pasteur/ida/zellige/surfaceConstruction/construction/SurfacesExtraction.java
+++ b/src/main/java/fr/pasteur/ida/zellige/surfaceConstruction/construction/SurfacesExtraction.java
@@ -1,12 +1,5 @@
 package fr.pasteur.ida.zellige.surfaceConstruction.construction;
 
-import java.util.ArrayList;
-
-import fr.pasteur.ida.zellige.surfaceConstruction.element.os.OSE;
-import fr.pasteur.ida.zellige.surfaceConstruction.element.os.OSEStartingStatus;
-import fr.pasteur.ida.zellige.utils.*;
-import org.jzy3d.analysis.AnalysisLauncher;
-
 import fr.pasteur.ida.zellige.exception.FirstRoundConstructionException;
 import fr.pasteur.ida.zellige.exception.NoSurfaceFoundException;
 import fr.pasteur.ida.zellige.jzy3D.LocalMaximumsDisplay;
@@ -14,14 +7,21 @@ import fr.pasteur.ida.zellige.jzy3D.OSDisplay;
 import fr.pasteur.ida.zellige.jzy3D.SurfaceDisplay;
 import fr.pasteur.ida.zellige.surfaceConstruction.element.Pixels;
 import fr.pasteur.ida.zellige.surfaceConstruction.element.Surface;
+import fr.pasteur.ida.zellige.surfaceConstruction.element.os.OSE;
+import fr.pasteur.ida.zellige.surfaceConstruction.element.os.OSEStartingStatus;
 import fr.pasteur.ida.zellige.surfaceConstruction.element.os.OSList;
 import fr.pasteur.ida.zellige.surfaceConstruction.element.surfaceLine.SurfaceLine;
+import fr.pasteur.ida.zellige.utils.*;
+import net.imglib2.RandomAccessibleInterval;
 import net.imglib2.img.Img;
+import net.imglib2.img.ImgFactory;
 import net.imglib2.img.display.imagej.ImageJFunctions;
 import net.imglib2.type.NativeType;
-import net.imglib2.type.logic.BitType;
 import net.imglib2.type.numeric.RealType;
 import net.imglib2.type.numeric.integer.UnsignedShortType;
+import org.jzy3d.analysis.AnalysisLauncher;
+
+import java.util.ArrayList;
 
 
 
@@ -29,7 +29,7 @@ public class SurfacesExtraction
 {
 
     public static < T extends RealType< T > & NativeType< T > > void extract(
-            Img< T > input, double amplitude, double threshold, int sigmas,
+            RandomAccessibleInterval< T > input, ImgFactory<T> factory,double amplitude, double threshold, int sigmas,
             boolean zMapDisplay, int delta,
             boolean extractedStackDisplay,
             boolean heightMapDisplay,
@@ -60,7 +60,7 @@ public class SurfacesExtraction
                 System.out.println( "Surface number " + index );
 //                displaySurface(surface);
                 Img< UnsignedShortType > zMap = surface.getZMap();
-                projection( input, delta, zMap, zMapDisplay,
+                projection( input, factory,delta, zMap, zMapDisplay,
                         extractedStackDisplay, heightMapDisplay,
                         projectionDisplay, projectionType, index );
                 index++;
@@ -76,44 +76,6 @@ public class SurfacesExtraction
     }
 
 
-    //TODO Bad documentation
-
-    /**
-     * This method is used to binarize a grey-level image using Otsu Thresholding method in order to classify its pixels
-     * into foreground/background and then apply a anisotropic gaussian blur.
-     */
-
-    private static < T extends RealType< T > & NativeType< T > > Img< T > maximumSearch( Img< T > input, double percent, int sigma )
-    {
-        /* Local maximum detection using partial derivative */
-        Img< T > maximums = LocalMaximumDetection.findMaximums( input, input.factory() );
-        ImageJFunctions.show( maximums, "maximums" );
-
-        /* Grid of local maximums*/
-        Img< BitType > threshold = LocalOtsuClassification.find( input, percent );
-        //ImageJFunctions.show( threshold, "Local threshold " );
-
-        /* Maxima selection according to threshold*/
-        maximums = Threshold.classification( maximums.copy(), maximums.factory(), threshold );
-
-        /* Dilatation of the resulting image*/
-        Utils.gaussConvolution( maximums.copy(), maximums, new double[]{ sigma, sigma } );
-//        ImageJFunctions.show( maximums, "blurred maximums" );
-
-        /* New local maximum detection due to previous smoothing*/
-        maximums = LocalMaximumDetection.findMaximums( maximums.copy(), maximums.factory() );
-        return maximums;
-    }
-
-
-    /**
-     * Setting of the unfiltered local maximums : first the original stack is filtered with an anisotropic
-     * gaussian convolution,  the pixel values are normalized between 0 and 255 then the local maximums
-     * are detected.
-     * During these steps, the temporary local, sectional and global thresholds are computed.
-     */
-
-
     /**
      * Creates the OS for each orthogonal section from the maximum coordinates found according to the dimension.
      *
@@ -271,7 +233,7 @@ public class SurfacesExtraction
 
 
     private static < R extends RealType< R > & NativeType< R > > void projection(
-            Img< R > original, int delta, Img< UnsignedShortType > zMap,
+            RandomAccessibleInterval< R > original, ImgFactory<R> factory, int delta, Img< UnsignedShortType > zMap,
             boolean zMapDisplay,
             boolean extractedStackDisplay,
             boolean heightMapDisplay,
@@ -289,7 +251,7 @@ public class SurfacesExtraction
             ImageJFunctions.show( smoothed, "Elevation map" + "RS n°" + index );
         }
 
-        Img< R > extractedStack = StackProjection.getExtractedStack( original, smoothed, delta );
+        Img< R > extractedStack = StackProjection.getExtractedStack( original, factory, smoothed, delta );
         if ( extractedStackDisplay )
         {
             ImageJFunctions.show( extractedStack,
@@ -309,13 +271,13 @@ public class SurfacesExtraction
                             projectionType + " RS n°" + index );
                 }
 
-                Img< R > projection = StackProjection.projection1( original, heightMap, projectionType );
+                Img< R > projection = StackProjection.projection1( original,factory, heightMap, projectionType );
                 ImageJFunctions.show( projection, "Projection with " + projectionType + " RS n°" + index );
             }
             else
             {
 
-                Img< R > projection = StackProjection.projection2( StackProjection.getHeightMapStack( original, zMap, delta ), projectionType );
+                Img< R > projection = StackProjection.projection2( StackProjection.getHeightMapStack( original, factory,zMap, delta ),factory, projectionType );
                 ImageJFunctions.show( projection, "Projection with " + projectionType + " RS n°" + index );
 
             }
diff --git a/src/main/java/fr/pasteur/ida/zellige/utils/HistogramZ.java b/src/main/java/fr/pasteur/ida/zellige/utils/HistogramZ.java
new file mode 100644
index 0000000000000000000000000000000000000000..a094d3e335acb47b2937acebe213b51165f6a863
--- /dev/null
+++ b/src/main/java/fr/pasteur/ida/zellige/utils/HistogramZ.java
@@ -0,0 +1,206 @@
+package fr.pasteur.ida.zellige.utils;
+
+
+import java.util.ArrayList;
+import java.util.Arrays;
+
+import net.imglib2.IterableInterval;
+import net.imglib2.RealCursor;
+import net.imglib2.algorithm.Algorithm;
+import net.imglib2.algorithm.Benchmark;
+import net.imglib2.algorithm.stats.HistogramBinMapper;
+import net.imglib2.type.Type;
+
+/**
+ * Implements a Histogram over an Image.
+ *
+ * @author 2011 Larry Lindsey
+ * @author Larry Lindsey
+ */
+public class HistogramZ< T > implements Algorithm, Benchmark
+{
+    /**
+     * Hold the histogram itself.
+     */
+    private final int[] histogram;
+    /**
+     * The Cursor from which the histogram is to be calculated.
+     */
+    private final RealCursor< T > cursor;
+    /**
+     * The HistogramBinMapper, used to map Type values to histogram bin indices.
+     */
+    private final HistogramBinMapper< T > binMapper;
+    /**
+     * Processing time, milliseconds.
+     */
+    private long pTime = 0;
+
+    /**
+     * Create a Histogram using the given mapper, calculating from the given
+     * Cursor.
+     *
+     * @param mapper the HistogramBinMapper used to map Type values to histogram
+     *               bin indices.
+     * @param c      a Cursor corresponding to the Image from which the Histogram
+     *               will be calculated
+     */
+    public HistogramZ( final HistogramBinMapper< T > mapper,
+                       final RealCursor< T > c )
+    {
+        cursor = c;
+        binMapper = mapper;
+        histogram = new int[ binMapper.getNumBins() ];
+    }
+
+    /**
+     * Create a Histogram using the given mapper, calculating from the given
+     * Image.
+     *
+     * @param mapper the HistogramBinMapper used to map Type values to histogram
+     *               bin indices.
+     * @param image  an Image from which the Histogram will be calculated
+     */
+    public HistogramZ( final HistogramBinMapper< T > mapper,
+                       final IterableInterval< T > image )
+    {
+        this( mapper, image.cursor() );
+    }
+
+    /**
+     * Resets the histogram array and the Cursor.
+     */
+    public void reset()
+    {
+        Arrays.fill( histogram, 0 );
+        cursor.reset();
+    }
+
+    /**
+     * Returns the bin count corresponding to a given {@link Type}.
+     *
+     * @param t the Type corresponding to the requested
+     * @return The requested bin count.
+     */
+    public int getBin( final T t )
+    {
+        return getHistogram()[ binMapper.map( t ) ];
+    }
+
+    /**
+     * Returns the bin count given by the indicated bin index.
+     *
+     * @param i the index of the requested bin
+     * @return the bin count at the given index
+     */
+    public int getBin( final int i )
+    {
+        return getHistogram()[ i ];
+    }
+
+    /**
+     * Returns this Histogram's HistogramBinMapper.
+     *
+     * @return the HistogramBinMapper associated with this Histogram.
+     */
+    public HistogramBinMapper< T > getBinMapper()
+    {
+        return binMapper;
+    }
+
+    /**
+     * Returns the histogram array.
+     *
+     * @return the histogram array.
+     */
+    public int[] getHistogram()
+    {
+        return histogram;
+    }
+
+    /**
+     * Creates and returns the a Type whose value corresponds to the center of
+     * the bin indexed by i.
+     *
+     * @param i the requested bin index.
+     * @return a Type whose value corresponds to the requested bin center.
+     */
+    public T getBinCenter( final int i )
+    {
+        return getBinMapper().invMap( i );
+    }
+
+    /**
+     * Creates and returns a List containing Types that correspond to the
+     * centers of the histogram bins.
+     *
+     * @return a List containing Types that correspond to the centers of the
+     * histogram bins.
+     */
+    public ArrayList< T > getBinCenters()
+    {
+        final ArrayList< T > binCenters = new ArrayList< T >( histogram.length );
+        for ( int i = 0; i < histogram.length; ++ i )
+        {
+            binCenters.add( i, getBinMapper().invMap( i ) );
+        }
+
+        return binCenters;
+    }
+
+    /**
+     * Returns the number of bins in this Histogram.
+     *
+     * @return the number of bins in this Histogram
+     */
+    public int getNumBins()
+    {
+        return getBinMapper().getNumBins();
+    }
+
+    @Override
+    public boolean checkInput()
+    {
+        return true;
+    }
+
+    @Override
+    public String getErrorMessage()
+    {
+        return null;
+    }
+
+    @Override
+    public boolean process()
+    {
+        final long startTime = System.currentTimeMillis();
+        int index;
+
+        while ( cursor.hasNext() )
+        {
+            cursor.fwd();
+            index = binMapper.map( cursor.get() );
+            /*
+             * The following check makes this run for IntegerTypes at 3 to 4
+             * longer than the manual case on my machine. This is a necessary
+             * check, but if this takes too long, it might be worthwhile to
+             * separate out an UncheckedHistogram, which would instead throw an
+             * ArrayOutOfBoundsException.
+             */
+            if ( index >= 0 && index < histogram.length )
+            {
+                ++ histogram[ index ];
+            }
+        }
+
+        pTime = System.currentTimeMillis() - startTime;
+        return true;
+    }
+
+    @Override
+    public long getProcessingTime()
+    {
+        return pTime;
+    }
+
+}
\ No newline at end of file
diff --git a/src/main/java/fr/pasteur/ida/zellige/utils/LocalOtsuClassification.java b/src/main/java/fr/pasteur/ida/zellige/utils/LocalOtsuClassification.java
index 8adad8ef99f8f7fc944bc1e98f35a7a30a1b723d..c7a4ecb329d9af6847222aaf45bad4fac5673e09 100644
--- a/src/main/java/fr/pasteur/ida/zellige/utils/LocalOtsuClassification.java
+++ b/src/main/java/fr/pasteur/ida/zellige/utils/LocalOtsuClassification.java
@@ -1,8 +1,6 @@
 package fr.pasteur.ida.zellige.utils;
 
-import net.imglib2.Cursor;
-import net.imglib2.Interval;
-import net.imglib2.RandomAccess;
+import net.imglib2.*;
 import net.imglib2.algorithm.util.Grids;
 import net.imglib2.img.Img;
 import net.imglib2.img.ImgFactory;
@@ -23,45 +21,43 @@ public class LocalOtsuClassification
 {
 
     /**
-     *
-     * @param input the input {@link Img}
+     * @param input   the input {@link Img}
      * @param percent the percentage of global otsu value
-     * @param <T> the type on the input
+     * @param <T>     the type on the input
      * @return a 3D binary {@link Img<BitType>}
      */
-    public static  < T extends RealType< T > & NativeType< T > > Img< BitType > find( final Img< T > input , double percent)
+    public static < T extends RealType< T > & NativeType< T > > Img< BitType > find( final RandomAccessibleInterval< T > input,ImgFactory<T> factory, double percent )
     {
         // Prepare output.
-        final ImgFactory< BitType > factory =  Util.getArrayOrCellImgFactory( input, new BitType() );
-        Img< BitType > binary = factory.create( input );
-        new OtsuThreshold<>( input, binary, percent);
+        final ImgFactory< BitType > bitTypeImgFactory = Util.getArrayOrCellImgFactory( input, new BitType() );
+        Img< BitType > binary = bitTypeImgFactory.create( input );
+        new OtsuThreshold<>( input,factory, binary, percent );
         return binary;
     }
 
 
     /**
-     *
      * @param <T>
      */
     private static final class OtsuThreshold< T extends RealType< T > & NativeType< T > >
     {
-        private final Img< T > source;
-        private final Img< T > grid;
+        private final RandomAccessibleInterval< T > source;
+        private final ImgFactory<T> factory;
+        private final RandomAccessibleInterval< T > grid;
         private final Img< BitType > binary;
         private final double percent;
 
         /**
-         *
-         * @param source the input {@link Img}
-         * @param binary the resulting binary image as a {@link Img<BitType> }
+         * @param source  the input {@link Img}
+         * @param binary  the resulting binary image as a {@link Img<BitType> }
          * @param percent the percentage of global otsu value
          */
-        private OtsuThreshold( final Img< T > source, final Img< BitType > binary, double percent )
+        private OtsuThreshold( final RandomAccessibleInterval< T > source, final ImgFactory<T> factory, final Img< BitType > binary, double percent )
         {
             this.source = source;
+            this.factory = factory;
             this.binary = binary;
             this.percent = percent;
-            final ImgFactory< T > factory = Util.getArrayOrCellImgFactory( source, Util.getTypeFromInterval( source ) );
             this.grid = factory.create( source );
             run();
         }
@@ -71,27 +67,18 @@ public class LocalOtsuClassification
          */
         public void run()
         {
-            computeLocalThreshold( source, grid);
-            applyLocalThreshold( source, grid, binary, percent );
+            computeLocalThreshold( source, factory, grid );
+            applyLocalThreshold( Views.iterable( source) ,Views.iterable(  grid), binary, percent );
         }
 
-        /**
-         *
-         * @return the image containing the grid of local thresholds
-         */
-        public Img< BitType > getGrid()
-        {
-            return binary;
-        }
 
         /**
-         *
          * @param input the input {@link Img}
-         * @param grid the image containing the grid of local thresholds
-         * @param <T> the type on the input
+         * @param grid  the image containing the grid of local thresholds
+         * @param <T>   the type on the input
          */
         public static < T extends RealType< T > & NativeType< T > > void
-        computeLocalThreshold( Img< T > input, Img< T > grid )
+        computeLocalThreshold( RandomAccessibleInterval< T > input, ImgFactory<T> factory,RandomAccessibleInterval< T > grid )
         {
             long width = input.dimension( 0 );
             long height = input.dimension( 1 );
@@ -101,13 +88,13 @@ public class LocalOtsuClassification
 
             for ( Interval interval : intervals )
             {
-                computeLocalThreshold( input,  grid , interval );
+                computeLocalThreshold( input, grid, interval );
             }
+            grid = Utils.gaussConvolution(  grid,factory,  new double[]{ 5, 5, 1 } );
         }
 
-
         public static < T extends RealType< T > & NativeType< T > > void
-        computeLocalThreshold( Img< T > input, Img< T > grid , Interval interval )
+        computeLocalThreshold( RandomAccessibleInterval< T > input, RandomAccessibleInterval< T > grid, Interval interval )
         {
             IntervalView< T > viewSource = Views.offsetInterval( input, interval );
             IntervalView< T > viewGrid = Views.offsetInterval( grid, interval );
@@ -116,30 +103,27 @@ public class LocalOtsuClassification
         }
 
         /**
-         *
-         * @param input the input {@link Img}
-         * @param grid the image containing the grid of local thresholds
-         * @param binary the resulting binary image as a {@link Img<BitType> }
+         * @param input   the input image as an {@link IterableInterval}
+         * @param grid    the image containing the grid of local thresholds as an {@link IterableInterval}
+         * @param binary  the resulting binary image as a {@link RandomAccessibleInterval<BitType> }
          * @param percent the percentage of global otsu value
-         * @param <T> the input type
+         * @param <T>     the input type
          */
         public static < T extends RealType< T > & NativeType< T > > void
-        applyLocalThreshold( Img< T > input, Img< T > grid, Img<BitType> binary , double percent)
+        applyLocalThreshold( IterableInterval< T > input, IterableInterval< T > grid, RandomAccessibleInterval< BitType > binary, double percent )
         {
             double threshold = Threshold.getThreshold( input, percent ).getRealDouble();
-            Utils.gaussConvolution( grid.copy(), grid, new double[]{5, 5, 1} );
-            ImageJFunctions.show( grid.copy(), "grid" );
-            Cursor<T> sourceCursor = input.localizingCursor();
-            Cursor<T> outputCursor = grid.cursor();
+            Cursor< T > sourceCursor = input.localizingCursor();
+            Cursor< T > outputCursor = grid.cursor();
             RandomAccess< BitType > randomAccess = binary.randomAccess();
-            while(sourceCursor.hasNext())
+            while ( sourceCursor.hasNext() )
             {
                 sourceCursor.fwd();
                 outputCursor.fwd();
                 final double s = sourceCursor.get().getRealDouble();
                 final double o = outputCursor.get().getRealDouble();
                 {
-                if (s >= Math.max( o, threshold))// background subtraction
+                    if ( s >= Math.max( o, threshold ) )// background subtraction
                     {
                         randomAccess.setPosition( sourceCursor );
                         randomAccess.get().set( true );
@@ -150,9 +134,17 @@ public class LocalOtsuClassification
         }
 
 
+        /**
+         * @return the image containing the grid of local thresholds
+         */
+        public Img< BitType > getGrid()
+        {
+            return binary;
+        }
     }
 
 
+}
+
 
 
-}
diff --git a/src/main/java/fr/pasteur/ida/zellige/utils/MaximumAmplitudeClassification.java b/src/main/java/fr/pasteur/ida/zellige/utils/MaximumAmplitudeClassification.java
index 74e9bc4941cf457a3435b4043d4b202b3686e0da..bb15d9423c986d0e0fd71cf42bd76a12adf4e5a5 100644
--- a/src/main/java/fr/pasteur/ida/zellige/utils/MaximumAmplitudeClassification.java
+++ b/src/main/java/fr/pasteur/ida/zellige/utils/MaximumAmplitudeClassification.java
@@ -30,54 +30,57 @@ public class MaximumAmplitudeClassification
 {
 
     /**
-     *
-     * @param input - the input image as a {@link Img}
+     * @param input              - the input image as a {@link RandomAccessibleInterval}
      * @param amplitudeThreshold - the parameter user value to apply to the first mode of the image histogram.
-     * @param sizePercent the minimum percentage of pixels for isolated pixel elimination
-     * @param <T> - the type of the input
+     * @param sizePercent        the minimum percentage of pixels for isolated pixel elimination
+     * @param <T>                - the type of the input
      * @return a binary image of background foreground classification
      */
-    public static < T extends RealType< T > & NativeType< T > > Img< FloatType > find( final Img< T > input, double amplitudeThreshold, double sizePercent )
+    public static < T extends RealType< T > & NativeType< T > > Img< FloatType > find( final RandomAccessibleInterval< T > input,
+                                                                                       final ImgFactory< T > factory, double amplitudeThreshold, double sizePercent )
     {
         // Prepare output.
         Img< BitType > amplitude = new ArrayImgFactory<>( new BitType() ).create( input );
-        MaximumAmplitude< T > maximumAmplitude = new MaximumAmplitude<>( input, amplitude, amplitudeThreshold );
-        final ImgFactory< FloatType > factory = Util.getArrayOrCellImgFactory( input, new FloatType() );
-        Img< FloatType > output = factory.create( input );
+        MaximumAmplitude< T > maximumAmplitude = new MaximumAmplitude<>( input, factory, amplitude, amplitudeThreshold );
+        final ImgFactory< FloatType > floatTypeImgFactory = Util.getArrayOrCellImgFactory( input, new FloatType() );
+        Img< FloatType > output = floatTypeImgFactory.create( input );
         new BackgroundForegroundGrid( maximumAmplitude.getAmplitude(), output, sizePercent );
+        Utils.gaussConvolution( output.copy(), output, new double[]{ 0, 0, 0.5 } );
         ImageJFunctions.show( output, "new class amplitude" );
         return output;
     }
 
     /**
      * This class generate a binary image according to a percentage of the smallest mode value  of the input image
+     *
      * @param <T> the type of the input {@link Img}
      */
     private static final class MaximumAmplitude< T extends RealType< T > & NativeType< T > >
     {
-        private final Img< T > input;
+        private final RandomAccessibleInterval< T > input;
+        private final ImgFactory< T > factory;
         private final double amplitudeThreshold;
-        private Img< BitType > amplitude;
+        private RandomAccessibleInterval< BitType > amplitude;
 
-        public MaximumAmplitude( Img< T > input, Img< BitType > amplitude, double amplitudeThreshold )
+        public MaximumAmplitude( RandomAccessibleInterval< T > input, ImgFactory< T > factory, RandomAccessibleInterval< BitType > amplitude, double amplitudeThreshold )
         {
             this.input = input;
+            this.factory = factory;
             this.amplitude = amplitude;
             this.amplitudeThreshold = amplitudeThreshold;
             run();
         }
 
         /**
-         *
-         * @param max  the image containing only the values of local maximums
+         * @param max the image containing only the values of local maximums
          * @param min the image containing only the values of local minimums
          * @param <T> the type of both images {@link Img}
          * @return an image containing the value of the maximums amplitude
          */
         private static < T extends RealType< T > & NativeType< T > > Img< T > getAmplitude(
-                Img< T > max, RandomAccessibleInterval< T > min )
+                RandomAccessibleInterval< T > max, ImgFactory< T > factory, RandomAccessibleInterval< T > min )
         {
-            Img< T > amp = max.factory().create( max );
+            Img< T > amp = factory.create( max );
             RandomAccess< T > maxAccess = max.randomAccess();
             RandomAccess< T > minAccess = min.randomAccess();
             RandomAccess< T > ampAccess = amp.randomAccess();
@@ -95,7 +98,7 @@ public class MaximumAmplitudeClassification
                         if ( maxValue != 0 )
                         {
                             minAccess.setPosition( maxAccess );
-                            double amplitude = getAmplitude( maxValue, minAccess, (int) max.dimension( 2 ) );
+                            double amplitude = getAmplitude( maxValue, minAccess, ( int ) max.dimension( 2 ) );
                             ampAccess.setPosition( maxAccess );
                             ampAccess.get().setReal( amplitude );
                         }
@@ -106,16 +109,15 @@ public class MaximumAmplitudeClassification
         }
 
         /**
-         *
-         * @param maxValue the intensity value at local maximum positioned at minAccess location
+         * @param maxValue  the intensity value at local maximum positioned at minAccess location
          * @param minAccess - the random access positioned at  the same specific local maximum location
-         * @param <T> the type of the {@link RandomAccess}
+         * @param <T>       the type of the {@link RandomAccess}
          * @return the chosen amplitude value of the local maximum positioned at minAccess location
          */
-        private static < T extends RealType< T > & NativeType< T > > double getAmplitude( double maxValue, RandomAccess< T > minAccess , int depth)
+        private static < T extends RealType< T > & NativeType< T > > double getAmplitude( double maxValue, RandomAccess< T > minAccess, int depth )
         {
             double up = findValueUp( minAccess, maxValue );
-            double down = findValueDown( minAccess, maxValue , depth);
+            double down = findValueDown( minAccess, maxValue, depth );
             double result = Math.max( Math.abs( maxValue - up ), Math.abs( maxValue - down ) );
             if ( result == 0 )
             {
@@ -126,10 +128,9 @@ public class MaximumAmplitudeClassification
         }
 
         /**
-         *
          * @param minAccess the random access of the local minimum image
-         * @param maxValue the intensity value of the local maximum located at minAccess position
-         * @param <T> the type of the {@link RandomAccess}
+         * @param maxValue  the intensity value of the local maximum located at minAccess position
+         * @param <T>       the type of the {@link RandomAccess}
          * @return the amplitude value above the local maximum location in the Z dimension
          */
         private static < T extends RealType< T > & NativeType< T > > double findValueUp(
@@ -149,10 +150,9 @@ public class MaximumAmplitudeClassification
         }
 
         /**
-         *
          * @param minAccess the random access of the local minimum image
-         * @param maxValue the intensity value of the local maximum located at minAccess position
-         * @param <T> the type of the {@link RandomAccess}
+         * @param maxValue  the intensity value of the local maximum located at minAccess position
+         * @param <T>       the type of the {@link RandomAccess}
          * @return the amplitude value below the local maximum location in the Z dimension
          */
         private static < T extends RealType< T > & NativeType< T > > double findValueDown(
@@ -173,16 +173,16 @@ public class MaximumAmplitudeClassification
 
         public void run()
         {
-            Img< T > maximums = LocalMaximumDetection.findMaximums( input, input.factory() );
-            Img< T > minimums = LocalMinimumDetection.findMinimums( input, input.factory() );
-            Img< T > amp = getAmplitude( maximums, minimums );
+            Img< T > maximums = LocalMaximumDetection.findMaximums( input, factory );
+            Img< T > minimums = LocalMinimumDetection.findMinimums( input, factory );
+            Img< T > amp = getAmplitude( maximums, factory, minimums );
             T TMax = Threshold.getFirstMaxValue( maximums, false );
             TMax.mul( amplitudeThreshold );
             this.amplitude = Thresholder.threshold( amp.copy(), TMax, true, 2 );
 
         }
 
-        public Img< BitType > getAmplitude()
+        public RandomAccessibleInterval< BitType > getAmplitude()
         {
             return amplitude;
         }
@@ -194,17 +194,16 @@ public class MaximumAmplitudeClassification
      */
     private static final class BackgroundForegroundGrid
     {
-        private final Img< BitType > source;
-        private final Img< FloatType > output;
+        private final RandomAccessibleInterval< BitType > source;
+        private final RandomAccessibleInterval< FloatType > output;
         private final double sizePercent;
 
         /**
-         *
-         * @param source - the input  as a 2D {@link Img<BitType> }
-         * @param output - the output as a 2D {@link Img<FloatType> }
+         * @param source      - the input  as a 2D {@link Img<BitType> }
+         * @param output      - the output as a 2D {@link Img<FloatType> }
          * @param sizePercent - the minimum percentage of positive pixels
          */
-        private BackgroundForegroundGrid( final Img< BitType > source, Img< FloatType > output, double sizePercent )
+        private BackgroundForegroundGrid( final RandomAccessibleInterval< BitType > source, Img< FloatType > output, double sizePercent )
         {
             this.source = source;
             this.output = output;
@@ -213,10 +212,9 @@ public class MaximumAmplitudeClassification
         }
 
         /**
-         *
-         * @param intervalView  a sub 2D image
-         * @param sizePercent the minimum percentage of pixels for isolated pixel elimination
-         * @param <T> the type of the sub image
+         * @param intervalView a sub 2D image
+         * @param sizePercent  the minimum percentage of pixels for isolated pixel elimination
+         * @param <T>          the type of the sub image
          * @return true if the sub image does contain enough foreground pixels false otherwise
          */
         public static < T extends RealType< T > & NativeType< T > > boolean isForeground( IntervalView< T > intervalView, double sizePercent )
@@ -253,14 +251,13 @@ public class MaximumAmplitudeClassification
                 viewOutput.forEach( pixel -> pixel.setReal( foreground ) );
             }
             // Slight dilatation in z dimension
-            Utils.gaussConvolution( output.copy(), output, new double[]{ 0, 0, 0.5 } );
+
         }
 
         /**
-         *
          * @return the smoothed amplitude classified image
          */
-        public Img< FloatType > getOutput()
+        public RandomAccessibleInterval< FloatType > getOutput()
         {
             return output;
         }
diff --git a/src/main/java/fr/pasteur/ida/zellige/utils/StackProjection.java b/src/main/java/fr/pasteur/ida/zellige/utils/StackProjection.java
index cd78b4543e1b0ca920aff92268a9a0c3ac891d04..de605624d3050b6c4a99ae92d9c36c1182b9a0ce 100644
--- a/src/main/java/fr/pasteur/ida/zellige/utils/StackProjection.java
+++ b/src/main/java/fr/pasteur/ida/zellige/utils/StackProjection.java
@@ -79,7 +79,7 @@ public class StackProjection
      * @return
      */
     public static < T extends RealType< T > & NativeType< T > > Img< UnsignedShortType > getHeightMap
-    ( Img< T > stack, RandomAccessibleInterval< UnsignedShortType > zMap, String method, int delta )
+    ( RandomAccessibleInterval< T > stack, RandomAccessibleInterval< UnsignedShortType > zMap, String method, int delta )
     {
         ImgFactory< UnsignedShortType > imgFactory = new ArrayImgFactory<>( new UnsignedShortType() );
         Img< UnsignedShortType > heightMap =
@@ -120,11 +120,11 @@ public class StackProjection
      * @return
      */
     public static < T extends RealType< T > & NativeType< T > > Img< T > getHeightMapStack
-    ( Img< T > stack, RandomAccessibleInterval< UnsignedShortType > zMap,
+    ( RandomAccessibleInterval< T > stack, ImgFactory<T> factory, RandomAccessibleInterval< UnsignedShortType > zMap,
       int delta )
     {
         //the stack slice number depends of the user "delta" parameter.
-        Img< T > heightMapStack = stack.factory().create( stack.dimension( 0 ), stack.dimension( 1 ), ( delta * 2 ) + 1 );
+        Img< T > heightMapStack = factory.create( stack.dimension( 0 ), stack.dimension( 1 ), ( delta * 2 ) + 1 );
         RandomAccess< T > stackAccess = stack.randomAccess();
         RandomAccess< T > heightMapAccess = heightMapStack.randomAccess();
         RandomAccess< UnsignedShortType > zMapAccess = zMap.randomAccess();
@@ -176,9 +176,9 @@ public class StackProjection
 
 
     public static < T extends RealType< T > & NativeType< T > > Img< T > getExtractedStack
-            ( Img< T > stack, RandomAccessibleInterval< UnsignedShortType > zMap, int delta )
+            ( RandomAccessibleInterval< T > stack, ImgFactory<T> factory,RandomAccessibleInterval< UnsignedShortType > zMap, int delta )
     {
-        Img< T > heightMapStack = stack.factory().create( stack.dimension( 0 ), stack.dimension( 1 ), stack.dimension( 2 ) );
+        Img< T > heightMapStack = factory.create( stack.dimension( 0 ), stack.dimension( 1 ), stack.dimension( 2 ) );
         RandomAccess< T > stackAccess = stack.randomAccess();
         RandomAccess< T > heightmapAccess = heightMapStack.randomAccess();
         RandomAccess< UnsignedShortType > zMapAccess = zMap.randomAccess();
@@ -222,12 +222,12 @@ public class StackProjection
      * @return the projection of the original stack image according to yhe height map.
      */
     public static < T extends RealType< T > & NativeType< T > > Img< T > projection1(
-            Img< T > stack, Img< UnsignedShortType > heightMap, String projectionType )
+            RandomAccessibleInterval< T > stack,ImgFactory<T> factory, RandomAccessibleInterval< UnsignedShortType > heightMap, String projectionType )
     {
         //TODO replace stack by extracted stack
         if ( projectionType.equals( "MIP" ) )
         {
-            Img< T > projection = stack.factory().create( stack.dimension( 0 ), stack.dimension( 1 ) );
+            Img< T > projection = factory.create( stack.dimension( 0 ), stack.dimension( 1 ) );
             RandomAccess< T > stackAccess = stack.randomAccess();
             RandomAccess< T > projectionAccess = projection.randomAccess();
             RandomAccess< UnsignedShortType > heightmapAccess = heightMap.randomAccess();
@@ -259,9 +259,9 @@ public class StackProjection
     }
 
     public static < T extends RealType< T > & NativeType< T > > Img< T > projection2(
-            Img< T > extractedStack, String projectionType )
+            RandomAccessibleInterval< T > extractedStack, ImgFactory<T> factory,String projectionType )
     {
-        Img< T > projection = extractedStack.factory().create( extractedStack.dimension( 0 ), extractedStack.dimension( 1 ) );
+        Img< T > projection = factory.create( extractedStack.dimension( 0 ), extractedStack.dimension( 1 ) );
         if ( projectionType.equals( "median" ) )// sumSlices, median, standardDeviation
         {
             medianProjection( extractedStack, projection );
@@ -306,7 +306,7 @@ public class StackProjection
     }
 
     public static < T extends RealType< T > & NativeType< T > > void medianProjection
-            ( Img< T > input, Img< T > output )
+            ( RandomAccessibleInterval< T > input,RandomAccessibleInterval< T > output )
     {
         RandomAccess< T > inputAccess = input.randomAccess();
         RandomAccess< T > outputAccess = output.randomAccess();
diff --git a/src/main/java/fr/pasteur/ida/zellige/utils/SurfacePixelSelection.java b/src/main/java/fr/pasteur/ida/zellige/utils/SurfacePixelSelection.java
index c28fccf223904999e4a667db425bd5eb94e816d9..7691c3e86b9ba9d70b83c2bf3a844e38533da7aa 100644
--- a/src/main/java/fr/pasteur/ida/zellige/utils/SurfacePixelSelection.java
+++ b/src/main/java/fr/pasteur/ida/zellige/utils/SurfacePixelSelection.java
@@ -5,12 +5,18 @@ import fr.pasteur.ida.zellige.surfaceConstruction.element.Coordinate;
 import fr.pasteur.ida.zellige.surfaceConstruction.element.Pixels;
 import net.imglib2.RandomAccess;
 import net.imglib2.RandomAccessibleInterval;
+import net.imglib2.converter.Converters;
+import net.imglib2.converter.RealFloatConverter;
+import net.imglib2.converter.readwrite.RealFloatSamplerConverter;
 import net.imglib2.img.Img;
+import net.imglib2.img.ImgFactory;
+import net.imglib2.img.array.ArrayImgFactory;
 import net.imglib2.img.display.imagej.ImageJFunctions;
 import net.imglib2.type.NativeType;
 import net.imglib2.type.logic.BitType;
 import net.imglib2.type.numeric.RealType;
 import net.imglib2.type.numeric.real.FloatType;
+import net.imglib2.util.ImgUtil;
 import org.jzy3d.analysis.AnalysisLauncher;
 
 public class SurfacePixelSelection
@@ -26,14 +32,14 @@ public class SurfacePixelSelection
      * @param <T>
      * @return
      */
-    public static < T extends RealType< T > & NativeType< T > > Pixels[][] run( Img< T > source, double amplitude, double threshold, int sigmas )
+    public static < T extends RealType< T > & NativeType< T > > Pixels[][] run( RandomAccessibleInterval< T > source, double amplitude, double threshold, int sigmas )
     {
         /* Pretreatment of the image.*/
         Img< FloatType > normalized = pretreatment( source );
         ImageJFunctions.show( normalized.copy(), "converted, blurred and normalized " );
-
+        ImgFactory<FloatType> factory = normalized.factory();
         /* Local maximum search and selection */
-        Pixels[][] maximums  = maximumSearch( normalized, amplitude,threshold, sigmas );
+        Pixels[][] maximums  = maximumSearch( normalized, factory, amplitude,threshold, sigmas );
 
         /* Output  */
         displayMaximums( maximums );
@@ -46,14 +52,17 @@ public class SurfacePixelSelection
      * @param <T>
      * @return
      */
-    private static < T extends RealType< T > & NativeType< T > >Img<FloatType> pretreatment(Img< T > source)
+    private static < T extends RealType< T > & NativeType< T > >Img<FloatType> pretreatment( RandomAccessibleInterval< T > source )
     {
         // Conversion into FloatType for the derivative computation (negative values)
-        Img< FloatType > converted = Utils.convertIntoFloatType( source.copy() );
+        RandomAccessibleInterval< FloatType > converted = Converters.convert( source, new RealFloatSamplerConverter< T >() );
+
         //Image denoising with an anisotropic convolution
-        Utils.gaussConvolution( converted.copy(), converted, new double[]{1, 1, 0 } );
+        ImgFactory<FloatType> factory = new ArrayImgFactory<>(new FloatType());
+        Img <FloatType> c = Utils.gaussConvolution(  converted,factory, new double[]{1, 1, 0 } );
+        ImageJFunctions.show( c );
         // Normalization
-       return Utils.normalizeImage( converted );
+       return Utils.normalizeImage( c, c.factory());
     }
 
 
@@ -66,17 +75,17 @@ public class SurfacePixelSelection
      * @param <T>
      * @return
      */
-    private static < T extends RealType< T > & NativeType< T > > Pixels[][] maximumSearch( Img< T > input, double amplitude, double threshold, int sigma  )
+    private static < T extends RealType< T > & NativeType< T > > Pixels[][] maximumSearch( RandomAccessibleInterval< T > input,ImgFactory<T> factory, double amplitude, double threshold, int sigma  )
     {
         /* Grid of amplitude thresholds */
-        Img< FloatType > amplitudeImg = MaximumAmplitudeClassification.find(input,amplitude, 0.10);
+        Img< FloatType > amplitudeImg = MaximumAmplitudeClassification.find(input,factory, amplitude, 0.10);
         ImageJFunctions.show( amplitudeImg," amplitude" );
 
         /* Grid of local thresholds*/
-        Img< BitType > thresholds = LocalOtsuClassification.find( input, threshold );
+        Img< BitType > thresholds = LocalOtsuClassification.find( input, factory, threshold );
 
         /* Pixel selection */
-        Img<FloatType> finalList = Threshold.classification( amplitudeImg,  amplitudeImg.factory(), thresholds );
+        Img<FloatType> finalList = Threshold.classification( amplitudeImg, amplitudeImg.factory(), thresholds );
 
         /* Dilatation of the resulting image*/
         Utils.gaussConvolution( finalList.copy(), finalList, new double[]{ sigma, sigma, 1 } );
diff --git a/src/main/java/fr/pasteur/ida/zellige/utils/Threshold.java b/src/main/java/fr/pasteur/ida/zellige/utils/Threshold.java
index 0a93f3086d70448c697471594eb18bd512f4fd49..81afcb9a2eed4c52b89aa0808503bdf528d2e2c4 100644
--- a/src/main/java/fr/pasteur/ida/zellige/utils/Threshold.java
+++ b/src/main/java/fr/pasteur/ida/zellige/utils/Threshold.java
@@ -4,13 +4,11 @@ package fr.pasteur.ida.zellige.utils;
 import net.imglib2.Cursor;
 import net.imglib2.IterableInterval;
 import net.imglib2.RandomAccess;
-import net.imglib2.algorithm.stats.Histogram;
 import net.imglib2.algorithm.stats.Max;
 import net.imglib2.algorithm.stats.RealBinMapper;
 import net.imglib2.img.Img;
 import net.imglib2.img.ImgFactory;
 import net.imglib2.img.display.imagej.ImageJFunctions;
-import net.imglib2.script.slice.Slice2D;
 import net.imglib2.type.NativeType;
 import net.imglib2.type.logic.BitType;
 import net.imglib2.type.numeric.RealType;
@@ -102,9 +100,9 @@ public class Threshold
 // TODO implements Otsu-2D ???
 
 
-    public static < T extends RealType< T > & NativeType< T > > int[] getHistogram( Img< T > image, T min, T max, int numBins )
+    public static < T extends RealType< T > & NativeType< T > > int[] getHistogram( IterableInterval< T > image, T min, T max, int numBins )
     {
-        Histogram< T > H = new Histogram<>( new RealBinMapper<>( min, max, numBins ), image );
+        HistogramZ< T > H = new HistogramZ<>( new RealBinMapper<>( min, max, numBins ), image );
         H.process();
 
         return H.getHistogram();
@@ -176,25 +174,33 @@ public class Threshold
     }
 
 
-
-    public static < T extends RealType< T > & NativeType< T > > T getThreshold( Img< T > input, double percent )
+    public static < T extends RealType< T > & NativeType< T > > T getThreshold( IterableInterval< T > input, double percent )
     {
-
         T max = Threshold.findMax( input );
         T min = Threshold.findMin( input );
         int[] histogram = getHistogram( input, min, max, 256 );
+        T threshold = getThreshold( min, max, histogram );
+        threshold.mul( percent );
+        return threshold;
+    }
+
+    public static < T extends RealType< T > & NativeType< T > > T getThreshold( T min, T max, int[] histogram )
+    {
         histogram[ 0 ] = 0;
         int thresholdIndex = OtsuCelebiIndex( histogram );
+       return  getBinValueFromIndex( min, max, histogram, thresholdIndex );
+    }
 
+
+    public static < T extends RealType< T > & NativeType< T > > T getBinValueFromIndex( T min, T max, int[] histogram, int thresholdIndex )
+    {
         double binWidth = ( max.getRealDouble() - min.getRealDouble() ) / ( double ) histogram.length;
         double result = ( ( double ) ( thresholdIndex + 1 ) * binWidth + min.getRealDouble() );
-        T threshold = Utils.getTypeValue( result, input.firstElement() );
-        threshold.mul( percent );
+        T threshold = min.createVariable();
+        threshold.setReal( result );
         return threshold;
     }
 
-
-
     public static < T extends RealType< T > & NativeType< T > > Img< T > classification(
             IterableInterval< T > amplitude, ImgFactory< T > factory, Img< BitType > thresholds )
     {
@@ -202,13 +208,13 @@ public class Threshold
         Cursor< T > amplitudeCursor = amplitude.localizingCursor();
         RandomAccess< BitType > thresholdAccess = thresholds.randomAccess();
         RandomAccess< T > outputAccess = output.randomAccess();
-        while( amplitudeCursor.hasNext())
+        while ( amplitudeCursor.hasNext() )
         {
             amplitudeCursor.fwd();
-            if(amplitudeCursor.get().getRealDouble()!= 0)
+            if ( amplitudeCursor.get().getRealDouble() != 0 )
             {
                 thresholdAccess.setPosition( amplitudeCursor );
-                if(thresholdAccess.get().get())
+                if ( thresholdAccess.get().get() )
                 {
                     outputAccess.setPosition( amplitudeCursor );
                     outputAccess.get().setOne();
@@ -216,18 +222,15 @@ public class Threshold
             }
         }
         ImageJFunctions.show( output.copy(), "amp+thresh" );
-      return output;
+        return output;
     }
 
 
-
-
-
     public static < T extends RealType< T > & NativeType< T > > T getFirstMaxValue( Img< T > input, boolean zero )
     {
         T max = Threshold.findMax( input );
         T min = Threshold.findMin( input );
-        Histogram< T > H = new Histogram<>( new RealBinMapper<>( min, max, 256 ), input );
+        HistogramZ< T > H = new HistogramZ<>( new RealBinMapper<>( min, max, 256 ), input );
         H.process();
         int[] histogram = H.getHistogram();
 
diff --git a/src/main/java/fr/pasteur/ida/zellige/utils/Utils.java b/src/main/java/fr/pasteur/ida/zellige/utils/Utils.java
index 9463814acbcc1aa72e9005d947df0fa7d8d7cf19..ec80bf79177cabde79c4a4363030b0eede17590e 100644
--- a/src/main/java/fr/pasteur/ida/zellige/utils/Utils.java
+++ b/src/main/java/fr/pasteur/ida/zellige/utils/Utils.java
@@ -19,6 +19,7 @@ import net.imglib2.type.numeric.integer.IntType;
 import net.imglib2.type.numeric.integer.ShortType;
 import net.imglib2.type.numeric.real.DoubleType;
 import net.imglib2.type.numeric.real.FloatType;
+import net.imglib2.util.ImgUtil;
 import net.imglib2.view.Views;
 
 
@@ -48,6 +49,15 @@ public class Utils
         Gauss3.gauss( sigma, infiniteInput, output );
     }
 
+    public static < T extends RealType< T > & NativeType< T > > Img< T >
+    gaussConvolution( RandomAccessibleInterval< T > input, ImgFactory<T> factory, double[] sigma )
+    {
+        Img<T> output = factory.create( input );
+        RandomAccessible< T > infiniteInput = Views.extendValue( input, input.randomAccess().get() );
+        Gauss3.gauss( sigma, infiniteInput, output );
+        return output;
+    }
+
     /**
      * Method to apply a centered partial derivative function to a {@link RandomAccessibleInterval}
      * in the second dimension of a {@link RandomAccessibleInterval}.
@@ -75,12 +85,13 @@ public class Utils
      * @return a new Img with normalize pixel intensity values.
      */
     public static < T extends RealType< T > & NativeType< T > > Img< T > normalizeImage(
-            Img< T > image )
+            RandomAccessibleInterval< T > image, ImgFactory<T> factory )
     {
-        Img< T > normalizedImage = image.copy();
-        T min = image.firstElement().copy().createVariable();
+        Img< T > normalizedImage = factory.create( image );
+        ImgUtil.copy( image, normalizedImage );
+        T min = normalizedImage.firstElement().createVariable();
         min.setReal( 0 );
-        T max = image.firstElement().copy().createVariable();
+        T max = normalizedImage.firstElement().copy().createVariable();
         max.setReal( 255 );
         Normalize.normalize( normalizedImage, min, max );
         return normalizedImage;