diff --git a/README.md b/README.md
index 8620829..eb7d1fc 100644
--- a/README.md
+++ b/README.md
@@ -21,7 +21,7 @@ Big-FISH requires Python 3.6 or newer. Additionally, it has the following depend
 
 - numpy (>= 1.16.0)
 - scipy (>= 1.4.1)
-- scikit-learn (>= 0.21.0)
+- scikit-learn (>= 0.24.0)
 - scikit-image (>= 0.14.2)
 - matplotlib (>= 3.0.2)
 - pandas (>= 0.24.0)
diff --git a/bigfish/__init__.py b/bigfish/__init__.py
index 39b1340..f946f79 100644
--- a/bigfish/__init__.py
+++ b/bigfish/__init__.py
@@ -12,4 +12,4 @@
 # MINOR: new features
 # PATCH: backwards compatible bug fixes
 # MAJOR.MINOR.PATCHdev means a version under development
-__version__ = "0.6.1"
+__version__ = "0.6.2"
diff --git a/bigfish/classification/features.py b/bigfish/classification/features.py
index 53ed099..ab38a49 100644
--- a/bigfish/classification/features.py
+++ b/bigfish/classification/features.py
@@ -11,24 +11,37 @@
 from scipy import ndimage as ndi
 
 from skimage.morphology import binary_opening
-from skimage.morphology.selem import disk
+import skimage
+from sklearn.utils.fixes import parse_version
+if parse_version(skimage.__version__) < parse_version("0.19.0"):
+    from skimage.morphology.selem import disk
+else:
+    from skimage.morphology.footprints import disk
 
 import bigfish.stack as stack
 from .input_preparation import prepare_extracted_data
 
 
-# TODO allow RNA coordinates in float64 and int64
-
 # ### Main functions ###
 
-def compute_features(cell_mask, nuc_mask, ndim, rna_coord, smfish=None,
-                     voxel_size_yx=None, foci_coord=None,
-                     centrosome_coord=None,
-                     compute_distance=False, compute_intranuclear=False,
-                     compute_protrusion=False, compute_dispersion=False,
-                     compute_topography=False, compute_foci=False,
-                     compute_area=False, compute_centrosome=False,
-                     return_names=False):
+def compute_features(
+        cell_mask,
+        nuc_mask,
+        ndim,
+        rna_coord,
+        smfish=None,
+        voxel_size_yx=None,
+        foci_coord=None,
+        centrosome_coord=None,
+        compute_distance=False,
+        compute_intranuclear=False,
+        compute_protrusion=False,
+        compute_dispersion=False,
+        compute_topography=False,
+        compute_foci=False,
+        compute_area=False,
+        compute_centrosome=False,
+        return_names=False):
     """Compute requested features.
 
     Parameters
@@ -39,7 +52,7 @@ def compute_features(cell_mask, nuc_mask, ndim, rna_coord, smfish=None,
         Surface of the nucleus with shape (y, x).
     ndim : int
         Number of spatial dimensions to consider (2 or 3).
-    rna_coord : np.ndarray, np.int64
+    rna_coord : np.ndarray, np.int
         Coordinates of the detected spots with shape (nb_spots, 4) or
         (nb_spots, 3). One coordinate per dimension (zyx or yx dimensions)
         plus the index of the cluster assigned to the spot. If no cluster was
@@ -49,11 +62,11 @@ def compute_features(cell_mask, nuc_mask, ndim, rna_coord, smfish=None,
         Image of RNAs, with shape (y, x).
     voxel_size_yx : int, float or None
         Size of a voxel on the yx plan, in nanometer.
-    foci_coord : np.ndarray, np.int64
+    foci_coord : np.ndarray, np.int
         Array with shape (nb_foci, 5) or (nb_foci, 4). One coordinate per
         dimension for the foci centroid (zyx or yx coordinates), the number of
         spots detected in the foci and its index.
-    centrosome_coord : np.ndarray, np.int64
+    centrosome_coord : np.ndarray, np.int
         Coordinates of the detected centrosome with shape (nb_elements, 3) or
         (nb_elements, 2). One coordinate per dimension (zyx or yx dimensions).
         These coordinates are mandatory to compute centrosome related features.
@@ -83,22 +96,23 @@ def compute_features(cell_mask, nuc_mask, ndim, rna_coord, smfish=None,
 
     """
     # check parameters
-    stack.check_parameter(voxel_size_yx=(int, float, type(None)),
-                          compute_distance=bool,
-                          compute_intranuclear=bool,
-                          compute_protrusion=bool,
-                          compute_dispersion=bool,
-                          compute_topography=bool,
-                          compute_foci=bool,
-                          compute_area=bool,
-                          compute_centrosome=bool,
-                          return_names=bool)
+    stack.check_parameter(
+        voxel_size_yx=(int, float, type(None)),
+        compute_distance=bool,
+        compute_intranuclear=bool,
+        compute_protrusion=bool,
+        compute_dispersion=bool,
+        compute_topography=bool,
+        compute_foci=bool,
+        compute_area=bool,
+        compute_centrosome=bool,
+        return_names=bool)
     if smfish is not None:
         stack.check_array(smfish, ndim=[2, 3], dtype=[np.uint8, np.uint16])
         if smfish.ndim == 3:
             smfish = stack.maximum_projection(smfish)
     if foci_coord is not None:
-        stack.check_array(foci_coord, ndim=2, dtype=np.int64)
+        stack.check_array(foci_coord, ndim=2, dtype=[np.int32, np.int64])
 
     # prepare input data
     (cell_mask,
@@ -213,14 +227,15 @@ def compute_features(cell_mask, nuc_mask, ndim, rna_coord, smfish=None,
     return features
 
 
-def get_features_name(names_features_distance=False,
-                      names_features_intranuclear=False,
-                      names_features_protrusion=False,
-                      names_features_dispersion=False,
-                      names_features_topography=False,
-                      names_features_foci=False,
-                      names_features_area=False,
-                      names_features_centrosome=False):
+def get_features_name(
+        names_features_distance=False,
+        names_features_intranuclear=False,
+        names_features_protrusion=False,
+        names_features_dispersion=False,
+        names_features_topography=False,
+        names_features_foci=False,
+        names_features_area=False,
+        names_features_centrosome=False):
     """Return the current list of features names.
 
     Parameters
@@ -251,84 +266,99 @@ def get_features_name(names_features_distance=False,
 
     """
     # check parameters
-    stack.check_parameter(names_features_distance=bool,
-                          names_features_intranuclear=bool,
-                          names_features_protrusion=bool,
-                          names_features_dispersion=bool,
-                          names_features_topography=bool,
-                          names_features_foci=bool,
-                          names_features_area=bool,
-                          names_features_centrosome=bool)
+    stack.check_parameter(
+        names_features_distance=bool,
+        names_features_intranuclear=bool,
+        names_features_protrusion=bool,
+        names_features_dispersion=bool,
+        names_features_topography=bool,
+        names_features_foci=bool,
+        names_features_area=bool,
+        names_features_centrosome=bool)
 
     # initialization
     features_name = []
 
     # get feature names
     if names_features_distance:
-        features_name += ["index_mean_distance_cell",
-                          "index_median_distance_cell",
-                          "index_mean_distance_nuc",
-                          "index_median_distance_nuc"]
+        features_name += [
+            "index_mean_distance_cell",
+            "index_median_distance_cell",
+            "index_mean_distance_nuc",
+            "index_median_distance_nuc"]
 
     if names_features_intranuclear:
-        features_name += ["proportion_rna_in_nuc",
-                          "nb_rna_out_nuc",
-                          "nb_rna_in_nuc"]
+        features_name += [
+            "proportion_rna_in_nuc",
+            "nb_rna_out_nuc",
+            "nb_rna_in_nuc"]
 
     if names_features_protrusion:
-        features_name += ["index_rna_protrusion",
-                          "proportion_rna_protrusion",
-                          "protrusion_area"]
+        features_name += [
+            "index_rna_protrusion",
+            "proportion_rna_protrusion",
+            "protrusion_area"]
 
     if names_features_dispersion:
-        features_name += ["index_polarization",
-                          "index_dispersion",
-                          "index_peripheral_distribution"]
+        features_name += [
+            "index_polarization",
+            "index_dispersion",
+            "index_peripheral_distribution"]
 
     if names_features_topography:
-        features_name += ["index_rna_nuc_edge",
-                          "proportion_rna_nuc_edge"]
+        features_name += [
+            "index_rna_nuc_edge",
+            "proportion_rna_nuc_edge"]
 
         a = 500
         for b in range(1000, 3001, 500):
-            features_name += ["index_rna_nuc_radius_{}_{}".format(a, b),
-                              "proportion_rna_nuc_radius_{}_{}".format(a, b)]
+            features_name += [
+                "index_rna_nuc_radius_{}_{}".format(a, b),
+                "proportion_rna_nuc_radius_{}_{}".format(a, b)]
             a = b
 
         a = 0
         for b in range(500, 3001, 500):
-            features_name += ["index_rna_cell_radius_{}_{}".format(a, b),
-                              "proportion_rna_cell_radius_{}_{}".format(a, b)]
+            features_name += [
+                "index_rna_cell_radius_{}_{}".format(a, b),
+                "proportion_rna_cell_radius_{}_{}".format(a, b)]
             a = b
 
     if names_features_foci:
         features_name += ["proportion_rna_in_foci"]
 
     if names_features_area:
-        features_name += ["proportion_nuc_area",
-                          "cell_area",
-                          "nuc_area",
-                          "cell_area_out_nuc"]
+        features_name += [
+            "proportion_nuc_area",
+            "cell_area",
+            "nuc_area",
+            "cell_area_out_nuc"]
 
     if names_features_centrosome:
-        features_name += ["index_mean_distance_centrosome",
-                          "index_median_distance_centrosome",
-                          "index_rna_centrosome",
-                          "proportion_rna_centrosome",
-                          "index_centrosome_dispersion"]
+        features_name += [
+            "index_mean_distance_centrosome",
+            "index_median_distance_centrosome",
+            "index_rna_centrosome",
+            "proportion_rna_centrosome",
+            "index_centrosome_dispersion"]
 
     return features_name
 
 
 # ### Features functions ###
 
-def features_distance(rna_coord, distance_cell, distance_nuc, cell_mask, ndim,
-                      check_input=True):
+def features_distance(
+        rna_coord,
+        distance_cell,
+        distance_nuc,
+        cell_mask,
+        ndim,
+        check_input=True):
     """Compute distance related features.
 
     Parameters
     ----------
-    rna_coord : np.ndarray, np.int64
+    rna_coord : np.ndarray, np.int
         Coordinates of the detected RNAs with zyx or yx coordinates in the
         first 3 or 2 columns.
     distance_cell : np.ndarray, np.float32
@@ -360,7 +390,10 @@ def features_distance(rna_coord, distance_cell, distance_nuc, cell_mask, ndim,
         stack.check_parameter(ndim=int)
         if ndim not in [2, 3]:
             raise ValueError("'ndim' should be 2 or 3, not {0}.".format(ndim))
-        stack.check_array(rna_coord, ndim=2, dtype=np.int64)
+        stack.check_array(
+            rna_coord,
+            ndim=2,
+            dtype=[np.int32, np.int64])
         stack.check_array(
             distance_cell,
             ndim=2,
@@ -406,10 +439,10 @@ def features_in_out_nucleus(rna_coord, rna_coord_out_nuc, check_input=True):
 
     Parameters
     ----------
-    rna_coord : np.ndarray, np.int64
+    rna_coord : np.ndarray, np.int
         Coordinates of the detected RNAs with zyx or yx coordinates in the
         first 3 or 2 columns.
-    rna_coord_out_nuc : np.ndarray, np.int64
+    rna_coord_out_nuc : np.ndarray, np.int
         Coordinates of the detected RNAs with zyx or yx coordinates in the
         first 3 or 2 columns. Spots detected inside the nucleus are removed.
     check_input : bool
@@ -428,8 +461,9 @@ def features_in_out_nucleus(rna_coord, rna_coord_out_nuc, check_input=True):
     # check parameters
     stack.check_parameter(check_input=bool)
     if check_input:
-        stack.check_array(rna_coord, ndim=2, dtype=np.int64)
-        stack.check_array(rna_coord_out_nuc, ndim=2, dtype=np.int64)
+        stack.check_array(rna_coord, ndim=2, dtype=[np.int32, np.int64])
+        stack.check_array(
+            rna_coord_out_nuc, ndim=2, dtype=[np.int32, np.int64])
 
     # count total number of rna
     nb_rna = float(len(rna_coord))
@@ -451,13 +485,18 @@ def features_in_out_nucleus(rna_coord, rna_coord_out_nuc, check_input=True):
     return features
 
 
-def features_protrusion(rna_coord, cell_mask, nuc_mask, ndim, voxel_size_yx,
-                        check_input=True):
+def features_protrusion(
+        rna_coord,
+        cell_mask,
+        nuc_mask,
+        ndim,
+        voxel_size_yx,
+        check_input=True):
     """Compute protrusion related features.
 
     Parameters
     ----------
-    rna_coord : np.ndarray, np.int64
+    rna_coord : np.ndarray, np.int
         Coordinates of the detected RNAs with zyx or yx coordinates in the
         first 3 or 2 columns.
     cell_mask : np.ndarray, bool
@@ -485,7 +524,6 @@ def features_protrusion(rna_coord, cell_mask, nuc_mask, ndim, voxel_size_yx,
     # TODO fin a better feature for the protrusion (idea: dilate region from
     #  centroid and stop when a majority of new pixels do not belong to the
     #  cell).
-
     # check parameters
     stack.check_parameter(check_input=bool)
     if check_input:
@@ -494,7 +532,7 @@ def features_protrusion(rna_coord, cell_mask, nuc_mask, ndim, voxel_size_yx,
             voxel_size_yx=(int, float))
         if ndim not in [2, 3]:
             raise ValueError("'ndim' should be 2 or 3, not {0}.".format(ndim))
-        stack.check_array(rna_coord, ndim=2, dtype=np.int64)
+        stack.check_array(rna_coord, ndim=2, dtype=[np.int32, np.int64])
         stack.check_array(cell_mask, ndim=2, dtype=bool)
         stack.check_array(nuc_mask, ndim=2, dtype=bool)
 
@@ -534,8 +572,15 @@ def features_protrusion(rna_coord, cell_mask, nuc_mask, ndim, voxel_size_yx,
     return features
 
 
-def features_dispersion(smfish, rna_coord, centroid_rna, cell_mask,
-                        centroid_cell, centroid_nuc, ndim, check_input=True):
+def features_dispersion(
+        smfish,
+        rna_coord,
+        centroid_rna,
+        cell_mask,
+        centroid_cell,
+        centroid_nuc,
+        ndim,
+        check_input=True):
     """Compute RNA Distribution Index features (RDI) described in:
 
     RDI Calculator: An analysis Tool to assess RNA distributions in cells,
@@ -545,16 +590,16 @@ def features_dispersion(smfish, rna_coord, centroid_rna, cell_mask,
     ----------
     smfish : np.ndarray, np.uint
         Image of RNAs, with shape (y, x).
-    rna_coord : np.ndarray, np.int64
+    rna_coord : np.ndarray, np.int
         Coordinates of the detected RNAs with zyx or yx coordinates in the
         first 3 or 2 columns.
-    centroid_rna : np.ndarray, np.int64
+    centroid_rna : np.ndarray, np.int
         Coordinates of the rna centroid with shape (2,) or (3,).
     cell_mask : np.ndarray, bool
         Surface of the cell with shape (y, x).
-    centroid_cell : np.ndarray, np.int64
+    centroid_cell : np.ndarray, np.int
         Coordinates of the cell centroid with shape (2,).
-    centroid_nuc : np.ndarray, np.int64
+    centroid_nuc : np.ndarray, np.int
         Coordinates of the nucleus centroid with shape (2,).
     ndim : int
         Number of spatial dimensions to consider.
@@ -578,11 +623,11 @@ def features_dispersion(smfish, rna_coord, centroid_rna, cell_mask,
         if ndim not in [2, 3]:
             raise ValueError("'ndim' should be 2 or 3, not {0}.".format(ndim))
         stack.check_array(smfish, ndim=2, dtype=[np.uint8, np.uint16])
-        stack.check_array(rna_coord, ndim=2, dtype=np.int64)
-        stack.check_array(centroid_rna, ndim=1, dtype=np.int64)
+        stack.check_array(rna_coord, ndim=2, dtype=[np.int32, np.int64])
+        stack.check_array(centroid_rna, ndim=1, dtype=[np.int32, np.int64])
         stack.check_array(cell_mask, ndim=2, dtype=bool)
-        stack.check_array(centroid_cell, ndim=1, dtype=np.int64)
-        stack.check_array(centroid_nuc, ndim=1, dtype=np.int64)
+        stack.check_array(centroid_cell, ndim=1, dtype=[np.int32, np.int64])
+        stack.check_array(centroid_nuc, ndim=1, dtype=[np.int32, np.int64])
 
     # case where no mRNAs are detected
     if len(rna_coord) == 0:
@@ -645,7 +690,7 @@ def _rmsd(coord, reference_coord):
 
     Parameters
     ----------
-    coord : np.ndarray, np.int64
+    coord : np.ndarray, np.int
         Coordinates with shape (nb_points, 2).
     reference_coord : np.ndarray, np.int64
         Reference coordinate to compute the distance from, with shape (2,).
@@ -664,14 +709,19 @@ def _rmsd(coord, reference_coord):
     return rmsd
 
 
-def features_topography(rna_coord, cell_mask, nuc_mask, cell_mask_out_nuc,
-                        ndim, voxel_size_yx,
-                        check_input=True):
+def features_topography(
+        rna_coord,
+        cell_mask,
+        nuc_mask,
+        cell_mask_out_nuc,
+        ndim,
+        voxel_size_yx,
+        check_input=True):
     """Compute topographic features.
 
     Parameters
     ----------
-    rna_coord : np.ndarray, np.int64
+    rna_coord : np.ndarray, np.int
         Coordinates of the detected RNAs with zyx or yx coordinates in the
         first 3 or 2 columns.
     cell_mask : np.ndarray, bool
@@ -717,7 +767,7 @@ def features_topography(rna_coord, cell_mask, nuc_mask, cell_mask_out_nuc,
             voxel_size_yx=(int, float))
         if ndim not in [2, 3]:
             raise ValueError("'ndim' should be 2 or 3, not {0}.".format(ndim))
-        stack.check_array(rna_coord, ndim=2, dtype=np.int64)
+        stack.check_array(rna_coord, ndim=2, dtype=[np.int32, np.int64])
         stack.check_array(cell_mask, ndim=2, dtype=bool)
         stack.check_array(nuc_mask, ndim=2, dtype=bool)
         stack.check_array(cell_mask_out_nuc, ndim=2, dtype=bool)
@@ -807,10 +857,10 @@ def features_foci(rna_coord, foci_coord, ndim, check_input=True):
 
     Parameters
     ----------
-    rna_coord : np.ndarray, np.int64
+    rna_coord : np.ndarray, np.int
         Coordinates of the detected RNAs with zyx or yx coordinates in the
         first 3 or 2 columns.
-    foci_coord : np.ndarray, np.int64
+    foci_coord : np.ndarray, np.int
         Array with shape (nb_foci, 5) or (nb_foci, 4). One coordinate per
         dimension for the foci centroid (zyx or yx coordinates), the number of
         spots detected in the foci and its index.
@@ -831,8 +881,8 @@ def features_foci(rna_coord, foci_coord, ndim, check_input=True):
         stack.check_parameter(ndim=int)
         if ndim not in [2, 3]:
             raise ValueError("'ndim' should be 2 or 3, not {0}.".format(ndim))
-        stack.check_array(rna_coord, ndim=2, dtype=np.int64)
-        stack.check_array(foci_coord, ndim=2, dtype=np.int64)
+        stack.check_array(rna_coord, ndim=2, dtype=[np.int32, np.int64])
+        stack.check_array(foci_coord, ndim=2, dtype=[np.int32, np.int64])
 
     if len(rna_coord) == 0 or len(foci_coord) == 0:
         features = (0.,)
@@ -897,15 +947,21 @@ def features_area(cell_mask, nuc_mask, cell_mask_out_nuc, check_input=True):
     return features
 
 
-def features_centrosome(smfish, rna_coord, distance_centrosome, cell_mask,
-                        ndim, voxel_size_yx, check_input=True):
+def features_centrosome(
+        smfish,
+        rna_coord,
+        distance_centrosome,
+        cell_mask,
+        ndim,
+        voxel_size_yx,
+        check_input=True):
     """Compute centrosome related features (in 2 dimensions).
 
     Parameters
     ----------
     smfish : np.ndarray, np.uint
         Image of RNAs, with shape (y, x).
-    rna_coord : np.ndarray, np.int64
+    rna_coord : np.ndarray, np.int
         Coordinates of the detected RNAs with zyx or yx coordinates in the
         first 3 or 2 columns.
     distance_centrosome : np.ndarray, np.float32
@@ -944,7 +1000,7 @@ def features_centrosome(smfish, rna_coord, distance_centrosome, cell_mask,
         if ndim not in [2, 3]:
             raise ValueError("'ndim' should be 2 or 3, not {0}.".format(ndim))
         stack.check_array(smfish, ndim=2, dtype=[np.uint8, np.uint16])
-        stack.check_array(rna_coord, ndim=2, dtype=np.int64)
+        stack.check_array(rna_coord, ndim=2, dtype=[np.int32, np.int64])
         stack.check_array(
             distance_centrosome,
             ndim=2,
diff --git a/bigfish/classification/input_preparation.py b/bigfish/classification/input_preparation.py
index 514b5b9..e1725e3 100644
--- a/bigfish/classification/input_preparation.py
+++ b/bigfish/classification/input_preparation.py
@@ -16,8 +16,12 @@
 
 # ### Input data ###
 
-def prepare_extracted_data(cell_mask, nuc_mask=None, ndim=None, rna_coord=None,
-                           centrosome_coord=None):
+def prepare_extracted_data(
+        cell_mask,
+        nuc_mask=None,
+        ndim=None,
+        rna_coord=None,
+        centrosome_coord=None):
     """Prepare data extracted from images.
 
     Parameters
@@ -29,12 +33,12 @@ def prepare_extracted_data(cell_mask, nuc_mask=None, ndim=None, rna_coord=None,
     ndim : int
         Number of spatial dimensions to consider (2 or 3). Mandatory if
         `rna_coord` is provided.
-    rna_coord : np.ndarray, np.int64
+    rna_coord : np.ndarray, np.int
         Coordinates of the detected spots with shape (nb_spots, 4) or
         (nb_spots, 3). One coordinate per dimension (zyx or yx dimensions)
         plus the index of the cluster assigned to the spot. If no cluster was
         assigned, value is -1.
-    centrosome_coord : np.ndarray, np.int64
+    centrosome_coord : np.ndarray, np.int
         Coordinates of the detected centrosome with shape (nb_elements, 3) or
         (nb_elements, 2). One coordinate per dimension (zyx or yx dimensions).
 
@@ -46,7 +50,7 @@ def prepare_extracted_data(cell_mask, nuc_mask=None, ndim=None, rna_coord=None,
         Distance map from the cell with shape (y, x), in pixels.
     distance_cell_normalized : np.ndarray, np.float32
         Normalized distance map from the cell with shape (y, x).
-    centroid_cell : np.ndarray, np.int64
+    centroid_cell : np.ndarray, np.int
         Coordinates of the cell centroid with shape (2,).
     distance_centroid_cell : np.ndarray, np.float32
         Distance map from the cell centroid with shape (y, x), in pixels.
@@ -58,20 +62,20 @@ def prepare_extracted_data(cell_mask, nuc_mask=None, ndim=None, rna_coord=None,
         Distance map from the nucleus with shape (y, x), in pixels.
     distance_nuc_normalized : np.ndarray, np.float32
         Normalized distance map from the nucleus with shape (y, x).
-    centroid_nuc : np.ndarray, np.int64
+    centroid_nuc : np.ndarray, np.int
         Coordinates of the nucleus centroid with shape (2,).
     distance_centroid_nuc : np.ndarray, np.float32
         Distance map from the nucleus centroid with shape (y, x), in pixels.
-    rna_coord_out_nuc : np.ndarray, np.int64
+    rna_coord_out_nuc : np.ndarray, np.int
         Coordinates of the detected spots with shape (nb_spots, 4) or
         (nb_spots, 3). One coordinate per dimension (zyx or yx dimensions)
         plus the index of the cluster assigned to the spot. If no cluster was
         assigned, value is -1. Spots detected inside the nucleus are removed.
-    centroid_rna : np.ndarray, np.int64
+    centroid_rna : np.ndarray, np.int
         Coordinates of the rna centroid with shape (2,) or (3,).
     distance_centroid_rna : np.ndarray, np.float32
         Distance map from the rna centroid with shape (y, x), in pixels.
-    centroid_rna_out_nuc : np.ndarray, np.int64
+    centroid_rna_out_nuc : np.ndarray, np.int
         Coordinates of the rna centroid (outside the nucleus) with shape (2,)
         or (3,).
     distance_centroid_rna_out_nuc : np.ndarray, np.float32
@@ -81,7 +85,6 @@ def prepare_extracted_data(cell_mask, nuc_mask=None, ndim=None, rna_coord=None,
         Distance map from the centrosome with shape (y, x), in pixels.
 
     """
-    # TODO allow RNA coordinates in float64 and int64
     # check parameters
     stack.check_parameter(ndim=(int, type(None)))
     if rna_coord is not None and ndim is None:
@@ -91,18 +94,24 @@ def prepare_extracted_data(cell_mask, nuc_mask=None, ndim=None, rna_coord=None,
     stack.check_array(
         cell_mask,
         ndim=2,
-        dtype=[np.uint8, np.uint16, np.int64, bool])
+        dtype=[np.uint8, np.uint16, np.int32, np.int64, bool])
     cell_mask = cell_mask.astype(bool)
     if nuc_mask is not None:
         stack.check_array(
             nuc_mask,
             ndim=2,
-            dtype=[np.uint8, np.uint16, np.int64, bool])
+            dtype=[np.uint8, np.uint16, np.int32, np.int64, bool])
         nuc_mask = nuc_mask.astype(bool)
     if rna_coord is not None:
-        stack.check_array(rna_coord, ndim=2, dtype=np.int64)
+        stack.check_array(
+            rna_coord,
+            ndim=2,
+            dtype=[np.int32, np.int64])
     if centrosome_coord is not None:
-        stack.check_array(centrosome_coord, ndim=2, dtype=np.int64)
+        stack.check_array(
+            centrosome_coord,
+            ndim=2,
+            dtype=[np.int32, np.int64])
 
     # build distance map from the cell boundaries
     distance_cell = ndi.distance_transform_edt(cell_mask)
@@ -110,7 +119,7 @@ def prepare_extracted_data(cell_mask, nuc_mask=None, ndim=None, rna_coord=None,
     distance_cell_normalized = distance_cell / distance_cell.max()
 
     # get cell centroid and a distance map from its localisation
-    centroid_cell = _get_centroid_surface(cell_mask)
+    centroid_cell = _get_centroid_surface(cell_mask).astype(rna_coord.dtype)
     distance_centroid_cell = _get_centroid_distance_map(
         centroid_cell,
         cell_mask)
@@ -129,7 +138,7 @@ def prepare_extracted_data(cell_mask, nuc_mask=None, ndim=None, rna_coord=None,
         distance_nuc_normalized = distance_nuc / distance_nuc.max()
 
         # get nucleus centroid and a distance map from its localisation
-        centroid_nuc = _get_centroid_surface(nuc_mask)
+        centroid_nuc = _get_centroid_surface(nuc_mask).astype(rna_coord.dtype)
         distance_centroid_nuc = _get_centroid_distance_map(
             centroid_nuc,
             cell_mask)
@@ -146,7 +155,7 @@ def prepare_extracted_data(cell_mask, nuc_mask=None, ndim=None, rna_coord=None,
 
         # get rna centroid
         if len(rna_coord) == 0:
-            centroid_rna = np.array([0] * ndim, dtype=np.int64)
+            centroid_rna = np.array([0] * ndim, dtype=rna_coord.dtype)
         else:
             centroid_rna = _get_centroid_rna(rna_coord, ndim)
 
@@ -164,7 +173,8 @@ def prepare_extracted_data(cell_mask, nuc_mask=None, ndim=None, rna_coord=None,
 
             # get rna centroid (outside nucleus)
             if len(rna_coord_out_nuc) == 0:
-                centroid_rna_out_nuc = np.array([0] * ndim, dtype=np.int64)
+                centroid_rna_out_nuc = np.array(
+                    [0] * ndim, dtype=rna_coord.dtype)
             else:
                 centroid_rna_out_nuc = _get_centroid_rna(
                     rna_coord_out_nuc,
@@ -234,7 +244,7 @@ def _get_centroid_surface(mask):
 
     Returns
     -------
-    centroid : np.ndarray, np.int64
+    centroid : np.ndarray, np.int
         Coordinates of the centroid with shape (2,).
 
     """
@@ -250,7 +260,7 @@ def _get_centroid_rna(rna_coord, ndim):
 
     Parameters
     ----------
-    rna_coord : np.ndarray, np.int64
+    rna_coord : np.ndarray, np.int
         Coordinates of the detected spots with shape (nb_spots, 4) or
         (nb_spots, 3). One coordinate per dimension (zyx or yx dimensions)
         plus the index of the cluster assigned to the spot. If no cluster was
@@ -260,12 +270,12 @@ def _get_centroid_rna(rna_coord, ndim):
 
     Returns
     -------
-    centroid_rna : np.ndarray, np.int64
+    centroid_rna : np.ndarray, np.int
         Coordinates of the rna centroid with shape (2,) or (3,).
 
     """
     # get rna centroids
-    centroid_rna = np.mean(rna_coord[:, :ndim], axis=0, dtype=np.int64)
+    centroid_rna = np.mean(rna_coord[:, :ndim], axis=0, dtype=rna_coord.dtype)
 
     return centroid_rna
 
@@ -275,7 +285,7 @@ def _get_centroid_distance_map(centroid, cell_mask):
 
     Parameters
     ----------
-    centroid : np.ndarray, np.int64
+    centroid : np.ndarray, np.int
         Coordinates of the centroid with shape (2,) or (3,).
     cell_mask : np.ndarray, bool
         Binary surface of the cell with shape (y, x).
@@ -308,7 +318,7 @@ def _get_centrosome_distance_map(centrosome_coord, cell_mask):
 
     Parameters
     ----------
-    centrosome_coord : np.ndarray, np.int64
+    centrosome_coord : np.ndarray, np.int
         Coordinates of the detected centrosome with shape (nb_elements, 3) or
         (nb_elements, 2). One coordinate per dimension (zyx or yx dimensions).
     cell_mask : np.ndarray, bool
diff --git a/bigfish/detection/__init__.py b/bigfish/detection/__init__.py
index d48d7d8..7b90428 100644
--- a/bigfish/detection/__init__.py
+++ b/bigfish/detection/__init__.py
@@ -30,6 +30,8 @@
 from .utils import get_object_radius_pixel
 from .utils import get_object_radius_nm
 from .utils import build_reference_spot
+from .utils import get_spot_volume
+from .utils import get_spot_surface
 from .utils import compute_snr_spots
 from .utils import get_breaking_point
 
@@ -62,6 +64,8 @@
     "get_object_radius_pixel",
     "get_object_radius_nm",
     "build_reference_spot",
+    "get_spot_volume",
+    "get_spot_surface",
     "compute_snr_spots",
     "get_breaking_point"]
 
diff --git a/bigfish/detection/cluster_detection.py b/bigfish/detection/cluster_detection.py
index 92e9487..73aad26 100644
--- a/bigfish/detection/cluster_detection.py
+++ b/bigfish/detection/cluster_detection.py
@@ -58,7 +58,10 @@ def detect_clusters(spots, voxel_size, radius=350, nb_min_spots=4):
     # TODO check that the behavior is the same with float64 and int64
     #  coordinates
     # check parameters
-    stack.check_array(spots, ndim=2, dtype=[np.float64, np.int64])
+    stack.check_array(
+        spots,
+        ndim=2,
+        dtype=[np.float32, np.float64, np.int32, np.int64])
     stack.check_parameter(
         voxel_size=(int, float, tuple, list),
         radius=int,
@@ -158,6 +161,9 @@ def _extract_information(clustered_spots):
         number of spots detected in the cluster and its index.
 
     """
+    # store dtype
+    dtype = clustered_spots.dtype
+
     # extract information for 3-d cluster...
     if clustered_spots.shape[1] == 4:
 
@@ -165,7 +171,7 @@ def _extract_information(clustered_spots):
         labels_clusters = np.unique(
             clustered_spots[clustered_spots[:, 3] != -1, 3])
         if labels_clusters.size == 0:
-            clusters = np.array([], dtype=np.int64).reshape((0, 5))
+            clusters = np.array([], dtype=dtype).reshape((0, 5))
             return clusters
 
         # shape information
@@ -177,7 +183,7 @@ def _extract_information(clustered_spots):
             nb_spots_cluster = len(spots_in_cluster)
             clusters.append([z_cluster, y_cluster, x_cluster,
                              nb_spots_cluster, label])
-        clusters = np.array(clusters, dtype=np.int64)
+        clusters = np.array(clusters, dtype=dtype)
 
     # ... and 2-d cluster
     else:
@@ -186,7 +192,7 @@ def _extract_information(clustered_spots):
         labels_clusters = np.unique(
             clustered_spots[clustered_spots[:, 2] != -1, 2])
         if labels_clusters.size == 0:
-            clusters = np.array([], dtype=np.int64).reshape((0, 4))
+            clusters = np.array([], dtype=dtype).reshape((0, 4))
             return clusters
 
         # shape information
@@ -197,6 +203,6 @@ def _extract_information(clustered_spots):
             y_cluster, x_cluster = spots_in_cluster.mean(axis=0)
             nb_spots_cluster = len(spots_in_cluster)
             clusters.append([y_cluster, x_cluster, nb_spots_cluster, label])
-        clusters = np.array(clusters, dtype=np.int64)
+        clusters = np.array(clusters, dtype=dtype)
 
     return clusters
diff --git a/bigfish/detection/dense_decomposition.py b/bigfish/detection/dense_decomposition.py
index 88054fa..1a15292 100644
--- a/bigfish/detection/dense_decomposition.py
+++ b/bigfish/detection/dense_decomposition.py
@@ -28,8 +28,15 @@
 
 # ### Main function ###
 
-def decompose_dense(image, spots, voxel_size, spot_radius, kernel_size=None,
-                    alpha=0.5, beta=1, gamma=5):
+def decompose_dense(
+        image,
+        spots,
+        voxel_size,
+        spot_radius,
+        kernel_size=None,
+        alpha=0.5,
+        beta=1,
+        gamma=5):
     """Detect dense and bright regions with potential clustered spots and
     simulate a more realistic number of spots in these regions.
 
@@ -44,7 +51,7 @@ def decompose_dense(image, spots, voxel_size, spot_radius, kernel_size=None,
     ----------
     image : np.ndarray
         Image with shape (z, y, x) or (y, x).
-    spots : np.ndarray, np.int64
+    spots : np.ndarray
         Coordinate of the spots with shape (nb_spots, 3) or (nb_spots, 2)
         for 3-d or 2-d images respectively.
     voxel_size : int, float, Tuple(int, float) or List(int, float)
@@ -93,7 +100,7 @@ def decompose_dense(image, spots, voxel_size, spot_radius, kernel_size=None,
 
     Returns
     -------
-    spots : np.ndarray, np.int64
+    spots : np.ndarray
         Coordinate of the spots detected, with shape (nb_spots, 3) or
         (nb_spots, 2). One coordinate per dimension (zyx or yx coordinates).
     dense_regions : np.ndarray, np.int64
@@ -105,13 +112,15 @@ def decompose_dense(image, spots, voxel_size, spot_radius, kernel_size=None,
         Reference spot in 3-d or 2-d.
 
     """
-    # TODO allow/return float64 spots
     # check parameters
     stack.check_array(
         image,
         ndim=[2, 3],
         dtype=[np.uint8, np.uint16, np.float32, np.float64])
-    stack.check_array(spots, ndim=2, dtype=np.int64)
+    stack.check_array(
+        spots,
+        ndim=2,
+        dtype=[np.float32, np.float64, np.int32, np.int64])
     stack.check_parameter(
         voxel_size=(int, float, tuple, list),
         spot_radius=(int, float, tuple, list),
@@ -129,6 +138,9 @@ def decompose_dense(image, spots, voxel_size, spot_radius, kernel_size=None,
         raise ValueError("'gamma' should be a positive value, not {0}"
                          .format(gamma))
 
+    # store spots dtype
+    dtype = spots.dtype
+
     # check consistency between parameters
     ndim = image.ndim
     if ndim != spots.shape[1]:
@@ -158,7 +170,7 @@ def decompose_dense(image, spots, voxel_size, spot_radius, kernel_size=None,
 
     # case where no spot were detected
     if spots.size == 0:
-        dense_regions = np.array([], dtype=np.int64).reshape((0, ndim + 4))
+        dense_regions = np.array([], dtype=dtype).reshape((0, ndim + 4))
         reference_spot = np.zeros((5,) * ndim, dtype=image.dtype)
         return spots, dense_regions, reference_spot
 
@@ -189,7 +201,7 @@ def decompose_dense(image, spots, voxel_size, spot_radius, kernel_size=None,
 
     # case with an empty frame as reference spot
     if reference_spot.sum() == 0:
-        dense_regions = np.array([], dtype=np.int64).reshape((0, ndim + 4))
+        dense_regions = np.array([], dtype=dtype).reshape((0, ndim + 4))
         return spots, dense_regions, reference_spot
 
     # fit a gaussian function on the reference spot to be able to simulate it
@@ -214,7 +226,7 @@ def decompose_dense(image, spots, voxel_size, spot_radius, kernel_size=None,
 
     # case where no region where detected
     if regions_to_decompose.size == 0:
-        dense_regions = np.array([], dtype=np.int64).reshape((0, ndim + 4))
+        dense_regions = np.array([], dtype=dtype).reshape((0, ndim + 4))
         return spots, dense_regions, reference_spot
 
     # precompute gaussian function values
@@ -234,6 +246,7 @@ def decompose_dense(image, spots, voxel_size, spot_radius, kernel_size=None,
         amplitude=amplitude,
         background=background,
         precomputed_gaussian=precomputed_gaussian)
+    spots_in_regions = spots_in_regions.astype(dtype)
 
     # normally the number of detected spots should increase
     if len(spots_out_regions) + len(spots_in_regions) < len(spots):
@@ -261,7 +274,7 @@ def get_dense_region(image, spots, voxel_size, spot_radius, beta=1):
     ----------
     image : np.ndarray
         Image with shape (z, y, x) or (y, x).
-    spots : np.ndarray, np.int64
+    spots : np.ndarray
         Coordinate of the spots with shape (nb_spots, 3) or (nb_spots, 2).
     voxel_size : int, float, Tuple(int, float) or List(int, float)
         Size of a voxel, in nanometer. One value per spatial dimension (zyx or
@@ -286,7 +299,7 @@ def get_dense_region(image, spots, voxel_size, spot_radius, beta=1):
     dense_regions : np.ndarray
         Array with selected ``skimage.measure._regionprops._RegionProperties``
         objects.
-    spots_out_region : np.ndarray, np.int64
+    spots_out_region : np.ndarray
         Coordinate of the spots detected out of dense regions, with shape
         (nb_spots, 3) or (nb_spots, 2). One coordinate per dimension (zyx or
         yx coordinates).
@@ -294,13 +307,15 @@ def get_dense_region(image, spots, voxel_size, spot_radius, beta=1):
         Maximum size of the regions.
 
     """
-    # TODO allow/return float64 spots
     # check parameters
     stack.check_array(
         image,
         ndim=[2, 3],
         dtype=[np.uint8, np.uint16, np.float32, np.float64])
-    stack.check_array(spots, ndim=2, dtype=np.int64)
+    stack.check_array(
+        spots,
+        ndim=2,
+        dtype=[np.float32, np.float64, np.int32, np.int64])
     stack.check_parameter(
         voxel_size=(int, float, tuple, list),
         spot_radius=(int, float, tuple, list),
@@ -387,14 +402,14 @@ def _filter_connected_region(image, connected_component, spots):
         Image with shape (z, y, x) or (y, x).
     connected_component : np.ndarray, np.int64
         Image labelled with shape (z, y, x) or (y, x).
-    spots : np.ndarray, np.int64
+    spots : np.ndarray
         Coordinate of the spots with shape (nb_spots, 3) or (nb_spots, 2).
 
     Returns
     -------
     regions_filtered : np.ndarray
         Array with filtered skimage.measure._regionprops._RegionProperties.
-    spots_out_region : np.ndarray, np.int64
+    spots_out_region : np.ndarray
         Coordinate of the spots outside the regions with shape (nb_spots, 3)
         or (nb_spots, 2).
     max_region_size : int
@@ -437,14 +452,14 @@ def _filter_spot_out_candidate_regions(candidate_bbox, spots, nb_dim):
     ----------
     candidate_bbox : List[Tuple]
         List of Tuples with the bounding box coordinates.
-    spots : np.ndarray, np.int64
+    spots : np.ndarray
         Coordinate of the spots with shape (nb_spots, 3) or (nb_spots, 2).
     nb_dim : int
         Number of dimensions to consider (2 or 3).
 
     Returns
     -------
-    spots_out_region : np.ndarray, np.int64
+    spots_out_region : np.ndarray
         Coordinate of the spots outside the regions with shape (nb_spots, 3)
         or (nb_spots, 2).
     max_region_size : int
@@ -501,9 +516,14 @@ def _filter_spot_out_candidate_regions(candidate_bbox, spots, nb_dim):
 
 # ### Gaussian simulation ###
 
-def simulate_gaussian_mixture(image, candidate_regions, voxel_size, sigma,
-                              amplitude=100, background=0,
-                              precomputed_gaussian=None):
+def simulate_gaussian_mixture(
+        image,
+        candidate_regions,
+        voxel_size,
+        sigma,
+        amplitude=100,
+        background=0,
+        precomputed_gaussian=None):
     """Simulate as many gaussians as possible in the candidate dense regions in
     order to get a more realistic number of spots.
 
@@ -542,7 +562,6 @@ def simulate_gaussian_mixture(image, candidate_regions, voxel_size, sigma,
         average intensity value and its index.
 
     """
-    # TODO allow/return float64 spots
     # check parameters
     stack.check_array(
         image,
@@ -643,14 +662,23 @@ def simulate_gaussian_mixture(image, candidate_regions, voxel_size, sigma,
                             region_area, region_intensity, i_region])
 
     spots_in_regions = np.concatenate(spots_in_regions, axis=0)
+    spots_in_regions = spots_in_regions.astype(np.int64)
     regions = np.array(regions, dtype=np.int64)
 
     return spots_in_regions, regions
 
 
-def _gaussian_mixture_3d(image, region, voxel_size_z, voxel_size_yx, sigma_z,
-                         sigma_yx, amplitude, background, precomputed_gaussian,
-                         limit_gaussian=1000):
+def _gaussian_mixture_3d(
+        image,
+        region,
+        voxel_size_z,
+        voxel_size_yx,
+        sigma_z,
+        sigma_yx,
+        amplitude,
+        background,
+        precomputed_gaussian,
+        limit_gaussian=1000):
     """Fit as many 3-d gaussians as possible in a candidate region.
 
     Parameters
@@ -749,9 +777,15 @@ def _gaussian_mixture_3d(image, region, voxel_size_z, voxel_size_yx, sigma_z,
     return image_region, best_simulation, positions_gaussian
 
 
-def _gaussian_mixture_2d(image, region, voxel_size_yx, sigma_yx, amplitude,
-                         background, precomputed_gaussian,
-                         limit_gaussian=1000):
+def _gaussian_mixture_2d(
+        image,
+        region,
+        voxel_size_yx,
+        sigma_yx,
+        amplitude,
+        background,
+        precomputed_gaussian,
+        limit_gaussian=1000):
     """Fit as many 2-d gaussians as possible in a candidate region.
 
     Parameters
diff --git a/bigfish/detection/spot_detection.py b/bigfish/detection/spot_detection.py
index 30c6ccb..1d1ed90 100644
--- a/bigfish/detection/spot_detection.py
+++ b/bigfish/detection/spot_detection.py
@@ -22,9 +22,15 @@
 
 # ### Main function ###
 
-def detect_spots(images, threshold=None, remove_duplicate=True,
-                 return_threshold=False, voxel_size=None, spot_radius=None,
-                 log_kernel_size=None, minimum_distance=None):
+def detect_spots(
+        images,
+        threshold=None,
+        remove_duplicate=True,
+        return_threshold=False,
+        voxel_size=None,
+        spot_radius=None,
+        log_kernel_size=None,
+        minimum_distance=None):
     """Apply LoG filter followed by a Local Maximum algorithm to detect spots
     in a 2-d or 3-d image.
 
@@ -209,9 +215,13 @@ def detect_spots(images, threshold=None, remove_duplicate=True,
         return spots
 
 
-def _detect_spots_from_images(images, threshold=None, remove_duplicate=True,
-                              return_threshold=False, log_kernel_size=None,
-                              min_distance=None):
+def _detect_spots_from_images(
+        images,
+        threshold=None,
+        remove_duplicate=True,
+        return_threshold=False,
+        log_kernel_size=None,
+        min_distance=None):
     """Apply LoG filter followed by a Local Maximum algorithm to detect spots
     in a 2-d or 3-d image.
 
@@ -311,8 +321,8 @@ def _detect_spots_from_images(images, threshold=None, remove_duplicate=True,
         mask_local_max = masks[i]
 
         # detection
-        spots, _ = spots_thresholding(image_filtered, mask_local_max,
-                                      threshold, remove_duplicate)
+        spots, _ = spots_thresholding(
+            image_filtered, mask_local_max, threshold, remove_duplicate)
         all_spots.append(spots)
 
     # return threshold or not
@@ -382,8 +392,11 @@ def local_maximum_detection(image, min_distance):
     return mask
 
 
-def spots_thresholding(image, mask_local_max, threshold,
-                       remove_duplicate=True):
+def spots_thresholding(
+        image,
+        mask_local_max,
+        threshold,
+        remove_duplicate=True):
     """Filter detected spots and get coordinates of the remaining spots.
 
     In order to make the thresholding robust, it should be applied to a
@@ -422,7 +435,7 @@ def spots_thresholding(image, mask_local_max, threshold,
     stack.check_array(
         mask_local_max,
         ndim=[2, 3],
-        dtype=[bool])
+        dtype=bool)
     stack.check_parameter(
         threshold=(float, int, type(None)),
         remove_duplicate=bool)
@@ -499,7 +512,10 @@ def automated_threshold_setting(image, mask_local_max):
         image,
         ndim=[2, 3],
         dtype=[np.uint8, np.uint16, np.float32, np.float64])
-    stack.check_array(mask_local_max, ndim=[2, 3], dtype=[bool])
+    stack.check_array(
+        mask_local_max,
+        ndim=[2, 3],
+        dtype=bool)
 
     # get threshold values we want to test
     thresholds = _get_candidate_thresholds(image.ravel())
@@ -578,8 +594,12 @@ def _get_spot_counts(thresholds, value_spots):
     return thresholds, count_spots
 
 
-def get_elbow_values(images, voxel_size=None, spot_radius=None,
-                     log_kernel_size=None, minimum_distance=None):
+def get_elbow_values(
+        images,
+        voxel_size=None,
+        spot_radius=None,
+        log_kernel_size=None,
+        minimum_distance=None):
     """Get values to plot the elbow curve used to automatically set the
     threshold to detect spots.
 
diff --git a/bigfish/detection/spot_modeling.py b/bigfish/detection/spot_modeling.py
index 43ef9c4..d851974 100644
--- a/bigfish/detection/spot_modeling.py
+++ b/bigfish/detection/spot_modeling.py
@@ -10,8 +10,8 @@
 
 import bigfish.stack as stack
 
-from .utils import _get_spot_volume
-from .utils import _get_spot_surface
+from .utils import get_spot_volume
+from .utils import get_spot_surface
 from .utils import get_object_radius_pixel
 
 from scipy.special import erf
@@ -215,8 +215,11 @@ def initialize_grid(image_spot, voxel_size, return_centroid=False):
             return grid
 
 
-def _initialize_grid_3d(image_spot, voxel_size_z, voxel_size_yx,
-                        return_centroid=False):
+def _initialize_grid_3d(
+        image_spot,
+        voxel_size_z,
+        voxel_size_yx,
+        return_centroid=False):
     """Build a grid in nanometer to compute gaussian function values over a
     full volume.
 
@@ -351,8 +354,7 @@ def _initialize_background_amplitude(image_spot):
 
 # ### Pixel fitting ###
 
-def _objective_function(ndim, voxel_size, sigma_z, sigma_yx,
-                        amplitude):
+def _objective_function(ndim, voxel_size, sigma_z, sigma_yx, amplitude):
     """Design the objective function used to fit the gaussian function.
 
     Parameters
@@ -393,8 +395,12 @@ def _objective_function(ndim, voxel_size, sigma_z, sigma_yx,
     return f
 
 
-def _objective_function_3d(voxel_size_z, voxel_size_yx, sigma_z, sigma_yx,
-                           amplitude):
+def _objective_function_3d(
+        voxel_size_z,
+        voxel_size_yx,
+        sigma_z,
+        sigma_yx,
+        amplitude):
     """Design the objective function used to fit the gaussian function.
 
     Parameters
@@ -497,9 +503,18 @@ def f(grid, mu_z, mu_y, mu_x, sigma_z, sigma_yx, amplitude,
 
 
 # TODO add equations in the docstring
-def gaussian_3d(grid, mu_z, mu_y, mu_x, sigma_z, sigma_yx, voxel_size_z,
-                voxel_size_yx, amplitude, background,
-                precomputed=None):
+def gaussian_3d(
+        grid,
+        mu_z,
+        mu_y,
+        mu_x,
+        sigma_z,
+        sigma_yx,
+        voxel_size_z,
+        voxel_size_yx,
+        amplitude,
+        background,
+        precomputed=None):
     """Compute the gaussian function over the grid representing a volume V
     with shape (V_z, V_y, V_x).
 
@@ -668,8 +683,15 @@ def f(grid, mu_y, mu_x, sigma_yx, amplitude, background):
 
 
 # TODO add equations in the docstring
-def gaussian_2d(grid, mu_y, mu_x, sigma_yx, voxel_size_yx, amplitude,
-                background, precomputed=None):
+def gaussian_2d(
+        grid,
+        mu_y,
+        mu_x,
+        sigma_yx,
+        voxel_size_yx,
+        amplitude,
+        background,
+        precomputed=None):
     """Compute the gaussian function over the grid representing a surface S
     with shape (S_y, S_x).
 
@@ -747,11 +769,10 @@ def gaussian_2d(grid, mu_y, mu_x, sigma_yx, voxel_size_yx, amplitude,
     return values
 
 
+# TODO add equations in the docstring
 def _rescaled_erf(low, high, mu, sigma):
     """Rescaled the Error function along a specific axis.
 
-    # TODO add equations
-
     Parameters
     ----------
     low : np.ndarray, np.float
@@ -778,11 +799,10 @@ def _rescaled_erf(low, high, mu, sigma):
     return rescaled_erf
 
 
+# TODO add equations and optimization algorithm in the docstring
 def _fit_gaussian(f, grid, image_spot, p0, lower_bound=None, upper_bound=None):
     """Fit a gaussian function to a 3-d or 2-d image.
 
-    # TODO add equations and algorithm
-
     Parameters
     ----------
     f : func
@@ -810,7 +830,6 @@ def _fit_gaussian(f, grid, image_spot, p0, lower_bound=None, upper_bound=None):
     """
     # TODO check that we do not fit a 2-d gaussian function to a 3-d image or
     #  the opposite
-
     # compute lower bound and upper bound
     if lower_bound is None:
         lower_bound = [-np.inf for _ in p0]
@@ -939,7 +958,7 @@ def fit_subpixel(image, spots, voxel_size, spot_radius):
 
     Returns
     -------
-    spots_subpixel : np.ndarray, np.float64
+    spots_subpixel : np.ndarray
         Coordinate of the spots detected, with shape (nb_spots, 3) or
         (nb_spots, 2). One coordinate per dimension (zyx or yx coordinates).
 
@@ -949,7 +968,10 @@ def fit_subpixel(image, spots, voxel_size, spot_radius):
         image,
         ndim=[2, 3],
         dtype=[np.uint8, np.uint16, np.float32, np.float64])
-    stack.check_array(spots, ndim=2, dtype=[np.float64, np.int64])
+    stack.check_array(
+        spots,
+        ndim=2,
+        dtype=[np.float32, np.float64, np.int32, np.int64])
     stack.check_parameter(
         voxel_size=(int, float, tuple, list),
         spot_radius=(int, float, tuple, list))
@@ -1006,19 +1028,29 @@ def fit_subpixel(image, spots, voxel_size, spot_radius):
 
     # format results
     spots_subpixel = np.stack(spots_subpixel)
+    if spots.dtype == np.float32:
+        spots_subpixel = spots_subpixel.astype(np.float32)
+    else:
+        spots_subpixel = spots_subpixel.astype(np.float64)
 
     return spots_subpixel
 
 
-def _fit_subpixel_3d(image, coord, radius_to_crop, voxel_size_z, voxel_size_yx,
-                     spot_radius_z, spot_radius_yx):
+def _fit_subpixel_3d(
+        image,
+        coord,
+        radius_to_crop,
+        voxel_size_z,
+        voxel_size_yx,
+        spot_radius_z,
+        spot_radius_yx):
     """Fit a 3-d gaussian on a detected spot.
 
     Parameters
     ----------
     image : np.ndarray
         Image with shape (z, y, x).
-    coord : np.ndarray, np.int64
+    coord : np.ndarray
         Coordinate of the spot detected, with shape (3,). One coordinate per
         dimension (zyx coordinates).
     radius_to_crop : Tuple[float]
@@ -1041,7 +1073,7 @@ def _fit_subpixel_3d(image, coord, radius_to_crop, voxel_size_z, voxel_size_yx,
 
     """
     # extract spot image
-    image_spot, bbox_low = _get_spot_volume(
+    image_spot, bbox_low = get_spot_volume(
         image=image,
         spot_z=coord[0],
         spot_y=coord[1],
@@ -1083,15 +1115,19 @@ def _fit_subpixel_3d(image, coord, radius_to_crop, voxel_size_z, voxel_size_yx,
     return new_coord
 
 
-def _fit_subpixel_2d(image, coord, radius_to_crop, voxel_size_yx,
-                     spot_radius_yx):
+def _fit_subpixel_2d(
+        image,
+        coord,
+        radius_to_crop,
+        voxel_size_yx,
+        spot_radius_yx):
     """Fit a 2-d gaussian on a detected spot.
 
     Parameters
     ----------
     image : np.ndarray
         Image with shape (y, x).
-    coord : np.ndarray, np.int64
+    coord : np.ndarray
         Coordinate of the spot detected, with shape (2,). One coordinate per
         dimension (yx coordinates).
     radius_to_crop : Tuple[float]
@@ -1110,7 +1146,7 @@ def _fit_subpixel_2d(image, coord, radius_to_crop, voxel_size_yx,
 
     """
     # extract spot image
-    image_spot, bbox_low = _get_spot_surface(
+    image_spot, bbox_low = get_spot_surface(
         image=image,
         spot_y=coord[0],
         spot_x=coord[1],
diff --git a/bigfish/detection/tests/test_utils.py b/bigfish/detection/tests/test_utils.py
index e115147..f9a028a 100644
--- a/bigfish/detection/tests/test_utils.py
+++ b/bigfish/detection/tests/test_utils.py
@@ -10,5 +10,7 @@
 # TODO add test for bigfish.detection.get_object_radius_pixel
 # TODO add test for bigfish.detection.get_object_radius_nm
 # TODO add test for bigfish.detection.build_reference_spot
+# TODO add test for bigfish.detection.get_spot_volume
+# TODO add test for bigfish.detection.get_spot_surface
 # TODO add test for bigfish.detection.compute_snr_spots
 # TODO add test for bigfish.detection.get_breaking_point
diff --git a/bigfish/detection/utils.py b/bigfish/detection/utils.py
index 919bee1..c921a86 100644
--- a/bigfish/detection/utils.py
+++ b/bigfish/detection/utils.py
@@ -37,7 +37,10 @@ def convert_spot_coordinates(spots, voxel_size):
     """
     # check parameters
     stack.check_parameter(voxel_size=(int, float, tuple, list))
-    stack.check_array(spots, ndim=2, dtype=[np.float64, np.int64])
+    stack.check_array(
+        spots,
+        ndim=2,
+        dtype=[np.float32, np.float64, np.int32, np.int64])
     dtype = spots.dtype
 
     # check consistency between parameters
@@ -211,7 +214,10 @@ def build_reference_spot(image, spots, voxel_size, spot_radius, alpha=0.5):
         image,
         ndim=[2, 3],
         dtype=[np.uint8, np.uint16, np.float32, np.float64])
-    stack.check_array(spots, ndim=2, dtype=[np.float64, np.int64])
+    stack.check_array(
+        spots,
+        ndim=2,
+        dtype=[np.float32, np.float64, np.int32, np.int64])
     stack.check_parameter(
         voxel_size=(int, float, tuple, list),
         spot_radius=(int, float, tuple, list),
@@ -303,8 +309,8 @@ def _build_reference_spot_3d(image, spots, radius, alpha):
         spot_z, spot_y, spot_x = candidate_spots[i_spot, :]
 
         # get the volume of the spot
-        image_spot, _, = _get_spot_volume(image, spot_z, spot_y, spot_x,
-                                          radius_z, radius_yx)
+        image_spot, _, = get_spot_volume(
+            image, spot_z, spot_y, spot_x, radius_z, radius_yx)
 
         # keep images that are not cropped by the borders
         if image_spot.shape == (z_shape, yx_shape, yx_shape):
@@ -329,22 +335,22 @@ def _build_reference_spot_3d(image, spots, radius, alpha):
     return reference_spot
 
 
-def _get_spot_volume(image, spot_z, spot_y, spot_x, radius_z, radius_yx):
+def get_spot_volume(image, spot_z, spot_y, spot_x, radius_z, radius_yx):
     """Get a subimage of a detected spot in 3 dimensions.
 
     Parameters
     ----------
     image : np.ndarray
         Image with shape (z, y, x).
-    spot_z : np.int64
+    spot_z : int or float
         Coordinate of the detected spot along the z axis.
-    spot_y : np.int64
+    spot_y : int or float
         Coordinate of the detected spot along the y axis.
-    spot_x : np.int64
+    spot_x : int or float
         Coordinate of the detected spot along the x axis.
-    radius_z : int
+    radius_z : int or float
         Radius in pixel of the detected spot, along the z axis.
-    radius_yx : int
+    radius_yx : int or float
         Radius in pixel of the detected spot, on the yx plan.
 
     Returns
@@ -414,7 +420,7 @@ def _build_reference_spot_2d(image, spots, radius, alpha):
         spot_y, spot_x = candidate_spots[i_spot, :]
 
         # get the volume of the spot
-        image_spot, _ = _get_spot_surface(image, spot_y, spot_x, radius_yx)
+        image_spot, _ = get_spot_surface(image, spot_y, spot_x, radius_yx)
 
         # keep images that are not cropped by the borders
         if image_spot.shape == (yx_shape, yx_shape):
@@ -438,18 +444,18 @@ def _build_reference_spot_2d(image, spots, radius, alpha):
     return reference_spot
 
 
-def _get_spot_surface(image, spot_y, spot_x, radius_yx):
+def get_spot_surface(image, spot_y, spot_x, radius_yx):
     """Get a subimage of a detected spot in 2 dimensions.
 
     Parameters
     ----------
     image : np.ndarray
         Image with shape (y, x).
-    spot_y : np.int64
+    spot_y : int or float
         Coordinate of the detected spot along the y axis.
-    spot_x : np.int64
+    spot_x : int or float
         Coordinate of the detected spot along the x axis.
-    radius_yx : int
+    radius_yx : int or float
         Radius in pixel of the detected spot, on the yx plan.
 
     Returns
@@ -516,7 +522,10 @@ def compute_snr_spots(image, spots, voxel_size, spot_radius):
         ndim=[2, 3],
         dtype=[np.uint8, np.uint16, np.float32, np.float64])
     stack.check_range_value(image, min_=0)
-    stack.check_array(spots, ndim=2, dtype=[np.float64, np.int64])
+    stack.check_array(
+        spots,
+        ndim=2,
+        dtype=[np.float32, np.float64, np.int32, np.int64])
     stack.check_parameter(
         voxel_size=(int, float, tuple, list),
         spot_radius=(int, float, tuple, list))
@@ -587,7 +596,7 @@ def compute_snr_spots(image, spots, voxel_size, spot_radius):
             spot_z = spot[0]
             radius_background_z = radius_background[0]
             max_signal = image_to_process[spot_z, spot_y, spot_x]
-            spot_background_, _ = _get_spot_volume(
+            spot_background_, _ = get_spot_volume(
                 image_to_process, spot_z, spot_y, spot_x,
                 radius_background_z, radius_background_yx)
             spot_background = spot_background_.copy()
@@ -606,7 +615,7 @@ def compute_snr_spots(image, spots, voxel_size, spot_radius):
 
         else:
             max_signal = image_to_process[spot_y, spot_x]
-            spot_background_, _ = _get_spot_surface(
+            spot_background_, _ = get_spot_surface(
                 image_to_process, spot_y, spot_x, radius_background_yx)
             spot_background = spot_background_.copy()
 
@@ -665,8 +674,14 @@ def get_breaking_point(x, y):
 
     """
     # check parameters
-    stack.check_array(x, ndim=1, dtype=[np.float64, np.int64])
-    stack.check_array(y, ndim=1, dtype=[np.float64, np.int64])
+    stack.check_array(
+        x,
+        ndim=1,
+        dtype=[np.float32, np.float64, np.int32, np.int64])
+    stack.check_array(
+        y,
+        ndim=1,
+        dtype=[np.float32, np.float64, np.int32, np.int64])
 
     # select threshold where curve break
     slope = (y[-1] - y[0]) / len(y)
diff --git a/bigfish/multistack/colocalization.py b/bigfish/multistack/colocalization.py
index b2ae091..2ec570b 100644
--- a/bigfish/multistack/colocalization.py
+++ b/bigfish/multistack/colocalization.py
@@ -20,8 +20,13 @@
 
 # ### Main function ###
 
-def detect_spots_colocalization(spots_1, spots_2, voxel_size, threshold=None,
-                                return_indices=False, return_threshold=False):
+def detect_spots_colocalization(
+        spots_1,
+        spots_2,
+        voxel_size,
+        threshold=None,
+        return_indices=False,
+        return_threshold=False):
     """Detect colocalized spots between two arrays of spot coordinates
     'spots_1' and 'spots_2'. Pairs of spots below a specific threshold are
     defined as colocalized.
@@ -76,8 +81,14 @@ def detect_spots_colocalization(spots_1, spots_2, voxel_size, threshold=None,
         return_threshold=bool)
 
     # check spots coordinates
-    stack.check_array(spots_1, ndim=2, dtype=[np.float64, np.int64])
-    stack.check_array(spots_2, ndim=2, dtype=[np.float64, np.int64])
+    stack.check_array(
+        spots_1,
+        ndim=2,
+        dtype=[np.float32, np.float64, np.int32, np.int64])
+    stack.check_array(
+        spots_2,
+        ndim=2,
+        dtype=[np.float32, np.float64, np.int32, np.int64])
 
     # convert spots coordinates in nanometer
     spots_1_nanometer = detection.convert_spot_coordinates(
@@ -194,8 +205,14 @@ def get_elbow_value_colocalized(spots_1, spots_2, voxel_size):
     """
     # check parameters
     stack.check_parameter(voxel_size=(int, float, tuple, list))
-    stack.check_array(spots_1, ndim=2, dtype=[np.float64, np.int64])
-    stack.check_array(spots_2, ndim=2, dtype=[np.float64, np.int64])
+    stack.check_array(
+        spots_1,
+        ndim=2,
+        dtype=[np.float32, np.float64, np.int32, np.int64])
+    stack.check_array(
+        spots_2,
+        ndim=2,
+        dtype=[np.float32, np.float64, np.int32, np.int64])
 
     # check consistency between parameters
     ndim = spots_1.shape[1]
diff --git a/bigfish/multistack/postprocess.py b/bigfish/multistack/postprocess.py
index 27786f0..b550c40 100644
--- a/bigfish/multistack/postprocess.py
+++ b/bigfish/multistack/postprocess.py
@@ -48,8 +48,11 @@ def identify_objects_in_region(mask, coord, ndim):
     stack.check_array(
         mask,
         ndim=2,
-        dtype=[np.uint8, np.uint16, np.int64, bool])
-    stack.check_array(coord, ndim=2, dtype=[np.int64, np.float64])
+        dtype=[np.uint8, np.uint16, np.int32, np.int64, bool])
+    stack.check_array(
+        coord,
+        ndim=2,
+        dtype=[np.float32, np.float64, np.int32, np.int64])
 
     # check number of dimensions
     if ndim not in [2, 3]:
@@ -66,7 +69,7 @@ def identify_objects_in_region(mask, coord, ndim):
         mask = mask.astype(bool)
 
     # cast coordinates dtype if necessary
-    if coord.dtype == np.int64:
+    if coord.dtype in [np.int32, np.int64]:
         coord_int = coord
     else:
         coord_int = np.round(coord).astype(np.int64)
@@ -119,7 +122,10 @@ def remove_transcription_site(rna, clusters, nuc_mask, ndim):
 
     """
     # check parameters
-    stack.check_array(rna, ndim=2, dtype=[np.int64, np.float64])
+    stack.check_array(
+        rna,
+        ndim=2,
+        dtype=[np.float32, np.float64, np.int32, np.int64])
 
     # discriminate foci from transcription sites
     ts, foci = identify_objects_in_region(
@@ -161,11 +167,11 @@ def match_nuc_cell(nuc_label, cell_label, single_nuc, cell_alone):
     stack.check_array(
         nuc_label,
         ndim=[2, 3],
-        dtype=[np.uint8, np.uint16, np.int64])
+        dtype=[np.uint8, np.uint16, np.int32, np.int64])
     stack.check_array(
         cell_label,
         ndim=[2, 3],
-        dtype=[np.uint8, np.uint16, np.int64])
+        dtype=[np.uint8, np.uint16, np.int32, np.int64])
 
     # initialize new labelled images
     new_nuc_label = np.zeros_like(nuc_label)
@@ -252,9 +258,16 @@ def _get_most_frequent_value(array):
 
 # ### Cell extraction ###
 
-def extract_cell(cell_label, ndim, nuc_label=None, rna_coord=None,
-                 others_coord=None, image=None, others_image=None,
-                 remove_cropped_cell=True, check_nuc_in_cell=True):
+def extract_cell(
+        cell_label,
+        ndim,
+        nuc_label=None,
+        rna_coord=None,
+        others_coord=None,
+        image=None,
+        others_image=None,
+        remove_cropped_cell=True,
+        check_nuc_in_cell=True):
     """Extract cell-level results for an image.
 
     The function gathers different segmentation and detection results obtained
@@ -315,14 +328,17 @@ def extract_cell(cell_label, ndim, nuc_label=None, rna_coord=None,
     stack.check_array(
         cell_label,
         ndim=2,
-        dtype=[np.uint8, np.uint16, np.int64])
+        dtype=[np.uint8, np.uint16, np.int32, np.int64])
     if nuc_label is not None:
         stack.check_array(
             nuc_label,
             ndim=2,
-            dtype=[np.uint8, np.uint16, np.int64])
+            dtype=[np.uint8, np.uint16, np.int32, np.int64])
     if rna_coord is not None:
-        stack.check_array(rna_coord, ndim=2, dtype=[np.int64, np.float64])
+        stack.check_array(
+            rna_coord,
+            ndim=2,
+            dtype=[np.float32, np.float64, np.int32, np.int64])
     if image is not None:
         stack.check_array(image, ndim=2, dtype=[np.uint8, np.uint16])
     actual_keys = ["cell_id", "bbox", "cell_coord", "cell_mask", "nuc_coord",
@@ -335,16 +351,17 @@ def extract_cell(cell_label, ndim, nuc_label=None, rna_coord=None,
             else:
                 actual_keys.append(key)
             array = others_coord[key]
-            stack.check_array(array, ndim=2, dtype=[np.int64, np.float64])
+            stack.check_array(
+                array,
+                ndim=2,
+                dtype=[np.float32, np.float64, np.int32, np.int64])
             if array.shape[1] < ndim:
                 warnings.warn("Array in 'others_coord' have less coordinates "
                               "({0}) than the minimum number of spatial "
                               "dimension we consider ({1})."
                               .format(array.shape[1], ndim),
                               UserWarning)
-    # TODO allow boolean for 'others_image'
-    # TODO bug if 'image' is None but not 'others_image'
-    if others_image is not None:
+    if others_image is not None and image is not None:
         for key in others_image:
             if key in actual_keys:
                 raise KeyError("Key {0} in 'others_image' is already taken. "
@@ -352,7 +369,10 @@ def extract_cell(cell_label, ndim, nuc_label=None, rna_coord=None,
             else:
                 actual_keys.append(key)
             image_ = others_image[key]
-            stack.check_array(image_, ndim=2, dtype=[np.uint8, np.uint16])
+            stack.check_array(
+                image_,
+                ndim=2,
+                dtype=[np.uint8, np.uint16, np.int32, np.int64, bool])
             if image_.shape != image.shape:
                 warnings.warn("Image in 'others_image' does not have the same "
                               "shape ({0}) than original image ({1})."
@@ -518,16 +538,21 @@ def _check_nucleus_in_cell(cell_mask, nuc_mask):
         return True
 
 
-# TODO fix docstring
-def extract_spots_from_frame(spots, z_lim=None, y_lim=None, x_lim=None):
+def extract_spots_from_frame(
+        spots,
+        ndim,
+        z_lim=None,
+        y_lim=None,
+        x_lim=None):
     """Get spots coordinates within a given frame.
 
     Parameters
     ----------
     spots : np.ndarray
-        Coordinate of the spots detected inside foci, with shape (nb_spots, 3)
-        or (nb_spots, 4). One coordinate per dimension (zyx coordinates) plus
-        the index of the foci if necessary.
+        Coordinate of the spots. One coordinate per dimension first (zyx
+        coordinates or yx coordinates) plus additional dimensions if necessary.
+    ndim : {2, 3}
+        Number of spatial dimension to consider.
     z_lim : tuple[int, int]
         Minimum and maximum coordinate of the frame along the z axis.
     y_lim : tuple[int, int]
@@ -538,17 +563,23 @@ def extract_spots_from_frame(spots, z_lim=None, y_lim=None, x_lim=None):
     Returns
     -------
     extracted_spots : np.ndarray
-        Coordinate of the spots detected inside foci, with shape (nb_spots, 3)
-        or (nb_spots, 4). One coordinate per dimension (zyx coordinates) plus
-        the index of the foci if necessary.
+        Coordinate of the spots. One coordinate per dimension first (zyx
+        coordinates or yx coordinates) plus additional dimensions if necessary.
 
     """
     # check parameters
-    stack.check_array(spots, ndim=2, dtype=[np.int64, np.float64])
+    stack.check_array(
+        spots,
+        ndim=2,
+        dtype=[np.float32, np.float64, np.int32, np.int64])
     stack.check_parameter(
         z_lim=(tuple, type(None)),
         y_lim=(tuple, type(None)),
         x_lim=(tuple, type(None)))
+    if ndim not in [2, 3]:
+        raise ValueError("'ndim' should be 2 or 3.")
+    if ndim == 2:
+        z_lim = None
 
     # extract spots
     extracted_spots = spots.copy()
@@ -557,19 +588,26 @@ def extract_spots_from_frame(spots, z_lim=None, y_lim=None, x_lim=None):
         extracted_spots = extracted_spots[z_lim[0] < extracted_spots[:, 0]]
         extracted_spots[:, 0] -= z_lim[0]
     if y_lim is not None:
-        extracted_spots = extracted_spots[extracted_spots[:, 1] < y_lim[1]]
-        extracted_spots = extracted_spots[y_lim[0] < extracted_spots[:, 1]]
-        extracted_spots[:, 1] -= y_lim[0]
+        extracted_spots = extracted_spots[
+            extracted_spots[:, ndim - 2] < y_lim[1]]
+        extracted_spots = extracted_spots[
+            y_lim[0] < extracted_spots[:, ndim - 2]]
+        extracted_spots[:, ndim - 2] -= y_lim[0]
     if x_lim is not None:
-        extracted_spots = extracted_spots[extracted_spots[:, 2] < x_lim[1]]
-        extracted_spots = extracted_spots[x_lim[0] < extracted_spots[:, 2]]
-        extracted_spots[:, 2] -= x_lim[0]
+        extracted_spots = extracted_spots[
+            extracted_spots[:, ndim - 1] < x_lim[1]]
+        extracted_spots = extracted_spots[
+            x_lim[0] < extracted_spots[:, ndim - 1]]
+        extracted_spots[:, ndim - 1] -= x_lim[0]
 
     return extracted_spots
 
 
-def summarize_extraction_results(fov_results, ndim, path_output=None,
-                                 delimiter=";"):
+def summarize_extraction_results(
+        fov_results,
+        ndim,
+        path_output=None,
+        delimiter=";"):
     """Summarize results extracted from an image and store them in a dataframe.
 
     Parameters
@@ -596,9 +634,11 @@ def summarize_extraction_results(fov_results, ndim, path_output=None,
     -------
     df : pd.DataFrame
         Dataframe with summarized results from the field of view, at the cell
-        level. At least `cell_id` (Unique id of the cell) is returned. Other
-        indicators are summarized if available:
+        level. At least `cell_id` (Unique id of the cell) and 'cell_area' (2-d
+        area of the cell, in pixel) are returned. Other indicators are
+        summarized if available:
 
+        * `nuc_area`: 2-d area of the nucleus, in pixel.
         * `nb_rna`: Number of detected rna in the cell.
         * `nb_rna_in_nuc`: Number of detected rna inside the nucleus.
         * `nb_rna_out_nuc`: Number of detected rna outside the nucleus.
@@ -630,12 +670,15 @@ def summarize_extraction_results(fov_results, ndim, path_output=None,
             continue
         others_coord = cell_results[key]
         if (not isinstance(others_coord, np.ndarray)
-                or others_coord.dtype not in [np.int64, np.float64]):
+                or others_coord.dtype not in [
+                    np.float32, np.float64, np.int32, np.int64]):
             continue
         _extra_coord[key] = []
 
     # summarize results at the cell level
     _cell_id = []
+    _cell_area = []
+    _nuc_area = []
     _nb_rna = []
     _nb_rna_in_nuc = []
     _nb_rna_out_nuc = []
@@ -643,20 +686,27 @@ def summarize_extraction_results(fov_results, ndim, path_output=None,
         # get cell id
         _cell_id.append(cell_results["cell_id"])
 
+        # get cell area
+        cell_mask = cell_results["cell_mask"]
+        cell_area = int(cell_mask.sum())
+        _cell_area.append(cell_area)
+
         # get rna coordinates and relative results
         if "rna_coord" in cell_results:
             rna_coord = cell_results["rna_coord"]
             _nb_rna.append(len(rna_coord))
 
-            # get rna in nucleus
+            # get rna in nucleus and nucleus area
             if "nuc_mask" in cell_results:
                 nuc_mask = cell_results["nuc_mask"]
                 rna_in_nuc, rna_out_nuc = identify_objects_in_region(
                     nuc_mask,
                     rna_coord,
                     ndim)
+                nuc_area = int(nuc_mask.sum())
                 _nb_rna_in_nuc.append(len(rna_in_nuc))
                 _nb_rna_out_nuc.append(len(rna_out_nuc))
+                _nuc_area.append(nuc_area)
 
         # get others coordinates
         for key in _extra_coord:
@@ -671,21 +721,17 @@ def summarize_extraction_results(fov_results, ndim, path_output=None,
         _nb_rna_in_nuc = [np.nan] * n
     if len(_nb_rna_out_nuc) == 0:
         _nb_rna_out_nuc = [np.nan] * n
+    if len(_nuc_area) == 0:
+        _nuc_area = [np.nan] * n
 
     # store minimum results in a dataframe
     result_summary = {"cell_id": _cell_id,
+                      "cell_area": _cell_area,
+                      "nuc_area": _nuc_area,
                       "nb_rna": _nb_rna,
                       "nb_rna_in_nuc": _nb_rna_in_nuc,
                       "nb_rna_out_nuc": _nb_rna_out_nuc}
 
-    # store available results on nucleus and rna
-    if len(_nb_rna) > 0:
-        result_summary["nb_rna"] = _nb_rna
-    if len(_nb_rna_in_nuc) > 0:
-        result_summary["nb_rna_in_nuc"] = _nb_rna_in_nuc
-    if len(_nb_rna_out_nuc) > 0:
-        result_summary["nb_rna_out_nuc"] = _nb_rna_out_nuc
-
     # store results from others elements detected in the cell
     for key in _extra_coord:
         result_summary["nb_{0}".format(key)] = _extra_coord[key]
@@ -734,7 +780,7 @@ def center_mask_coord(main, others=None):
     stack.check_array(
         main,
         ndim=2,
-        dtype=[np.uint8, np.uint16, np.int64, bool])
+        dtype=[np.uint8, np.uint16, np.int32, np.int64, bool])
     stack.check_parameter(others=(list, type(None)))
     if others is not None and len(others) != 0:
         for x in others:
@@ -743,7 +789,7 @@ def center_mask_coord(main, others=None):
             stack.check_array(
                 x,
                 ndim=2,
-                dtype=[np.uint8, np.uint16, np.int64, bool])
+                dtype=[np.uint8, np.uint16, np.int32, np.int64, bool])
 
     # initialize parameter
     marge = stack.get_margin_value()
@@ -825,7 +871,7 @@ def from_boundaries_to_surface(binary_boundaries):
     stack.check_array(
         binary_boundaries,
         ndim=2,
-        dtype=[np.uint8, np.uint16, np.int64, bool])
+        dtype=[np.uint8, np.uint16, np.int32, np.int64, bool])
 
     # from binary boundaries to binary surface
     binary_surface = ndi.binary_fill_holes(binary_boundaries)
@@ -851,7 +897,7 @@ def from_surface_to_boundaries(binary_surface):
     stack.check_array(
         binary_surface,
         ndim=2,
-        dtype=[np.uint8, np.uint16, np.int64, bool])
+        dtype=[np.uint8, np.uint16, np.int32, np.int64, bool])
     original_dtype = binary_surface.dtype
 
     # pad the binary surface in case object if on the edge
@@ -883,7 +929,7 @@ def from_binary_to_coord(binary):
 
     Returns
     -------
-    coord : np.ndarray, np.int64
+    coord : np.ndarray, np.int
         Array of boundaries coordinates with shape (nb_points, 2).
 
     """
@@ -891,14 +937,20 @@ def from_binary_to_coord(binary):
     stack.check_array(
         binary,
         ndim=2,
-        dtype=[np.uint8, np.uint16, np.int64, bool])
+        dtype=[np.uint8, np.uint16, np.int32, np.int64, bool])
+
+    # store dtype
+    if binary.dtype == np.int32:
+        dtype = np.int32
+    else:
+        dtype = np.int64
 
     # we enlarge the binary mask with one pixel to be sure the external
     # boundaries of the object still fit within the frame
     binary_ = np.pad(binary, [(1, 1)], mode="constant")
 
     # get external boundaries coordinates
-    coord = find_contours(binary_, level=0)[0].astype(np.int64)
+    coord = find_contours(binary_, level=0)[0].astype(dtype)
 
     # remove the pad
     coord -= 1
@@ -912,17 +964,17 @@ def complete_coord_boundaries(coord):
 
     Parameters
     ----------
-    coord : np.ndarray, np.int64
+    coord : np.ndarray, np.int
         Array of coordinates to complete, with shape (nb_points, 2).
 
     Returns
     -------
-    coord_completed : np.ndarray, np.int64
+    coord_completed : np.ndarray, np.int
         Completed coordinates arrays, with shape (nb_points, 2).
 
     """
     # check parameters
-    stack.check_array(coord, ndim=2, dtype=[np.int64])
+    stack.check_array(coord, ndim=2, dtype=[np.int32, np.int64])
 
     # for each array in the list, complete its coordinates using the scikit
     # image method 'polygon_perimeter'
@@ -930,6 +982,7 @@ def complete_coord_boundaries(coord):
     coord_y = coord_y[:, np.newaxis]
     coord_x = coord_x[:, np.newaxis]
     coord_completed = np.concatenate((coord_y, coord_x), axis=-1)
+    coord_completed = coord_completed.astype(coord.dtype)
 
     return coord_completed
 
@@ -943,7 +996,7 @@ def from_coord_to_frame(coord, external_coord=True):
 
     Parameters
     ----------
-    coord : np.ndarray, np.int64
+    coord : np.ndarray, np.int
         Array of cell boundaries coordinates with shape (nb_points, 2) or
         (nb_points, 3).
     external_coord : bool
@@ -981,8 +1034,11 @@ def from_coord_to_frame(coord, external_coord=True):
     return frame_shape, min_y, min_x, marge
 
 
-def from_coord_to_surface(cell_coord, nuc_coord=None, rna_coord=None,
-                          external_coord=True):
+def from_coord_to_surface(
+        cell_coord,
+        nuc_coord=None,
+        rna_coord=None,
+        external_coord=True):
     """Convert 2-d coordinates to a binary matrix with the surface of the
     object.
 
@@ -995,11 +1051,11 @@ def from_coord_to_surface(cell_coord, nuc_coord=None, rna_coord=None,
 
     Parameters
     ----------
-    cell_coord : np.ndarray, np.int64
+    cell_coord : np.ndarray, np.int
         Array of cell boundaries coordinates with shape (nb_points, 2).
-    nuc_coord : np.ndarray, np.int64
+    nuc_coord : np.ndarray, np.int
         Array of nucleus boundaries coordinates with shape (nb_points, 2).
-    rna_coord : np.ndarray, np.int64
+    rna_coord : np.ndarray, np.int
         Array of mRNAs coordinates with shape (nb_points, 2) or
         (nb_points, 3).
     external_coord : bool
@@ -1013,16 +1069,16 @@ def from_coord_to_surface(cell_coord, nuc_coord=None, rna_coord=None,
         Binary image of nucleus surface with shape (y, x).
     rna_binary : np.ndarray, bool
         Binary image of mRNAs localizations with shape (y, x).
-    new_rna_coord : np.ndarray, np.int64
+    new_rna_coord : np.ndarray, np.int
         Array of mRNAs coordinates with shape (nb_points, 2) or (nb_points, 3).
 
     """
     # check parameters
-    stack.check_array(cell_coord, ndim=2, dtype=[np.int64])
+    stack.check_array(cell_coord, ndim=2, dtype=[np.int32, np.int64])
     if nuc_coord is not None:
-        stack.check_array(nuc_coord, ndim=2, dtype=[np.int64])
+        stack.check_array(nuc_coord, ndim=2, dtype=[np.int32, np.int64])
     if rna_coord is not None:
-        stack.check_array(rna_coord, ndim=2, dtype=[np.int64])
+        stack.check_array(rna_coord, ndim=2, dtype=[np.int32, np.int64])
     stack.check_parameter(external_coord=bool)
 
     # center coordinates
@@ -1065,7 +1121,7 @@ def from_coord_to_surface(cell_coord, nuc_coord=None, rna_coord=None,
             rna_binary[rna_coord_[:, 0], rna_coord_[:, 1]] = True
         else:
             rna_binary[rna_coord_[:, 1], rna_coord_[:, 2]] = True
-        new_rna_coord = rna_coord_.copy()
+        new_rna_coord = rna_coord_.copy().astype(rna_coord.dtype)
 
     else:
         rna_binary = None
diff --git a/bigfish/multistack/preprocess.py b/bigfish/multistack/preprocess.py
index 5e18f97..2597e58 100644
--- a/bigfish/multistack/preprocess.py
+++ b/bigfish/multistack/preprocess.py
@@ -23,8 +23,11 @@
 
 # ### Building stack ###
 
-def build_stacks(data_map, input_dimension=None, sanity_check=False,
-                 return_origin=False):
+def build_stacks(
+        data_map,
+        input_dimension=None,
+        sanity_check=False,
+        return_origin=False):
     """Generator to build several stacks from recipe-folder pairs.
 
     To build a stack, a recipe should be linked to a directory including all
@@ -145,8 +148,12 @@ def build_stacks(data_map, input_dimension=None, sanity_check=False,
                 yield tensor
 
 
-def build_stack(recipe, input_folder, input_dimension=None, sanity_check=False,
-                i_fov=0):
+def build_stack(
+        recipe,
+        input_folder,
+        input_dimension=None,
+        sanity_check=False,
+        i_fov=0):
     """Build a 5-d stack from the same field of view (fov).
 
     The recipe dictionary for one field of view takes the form:
@@ -577,8 +584,8 @@ def build_stack_no_recipe(paths, input_dimension=None, sanity_check=False):
     if sanity_check:
         stack.check_array(
             tensor,
-            dtype=[np.uint8, np.uint16, np.uint32,
-                   np.int8, np.int16, np.int32,
+            dtype=[np.uint8, np.uint16, np.uint32, np.uint64,
+                   np.int8, np.int16, np.int32, np.int64,
                    np.float16, np.float32, np.float64,
                    bool],
             ndim=5,
diff --git a/bigfish/multistack/tests/test_postprocess.py b/bigfish/multistack/tests/test_postprocess.py
index cb9e21e..2574aa7 100644
--- a/bigfish/multistack/tests/test_postprocess.py
+++ b/bigfish/multistack/tests/test_postprocess.py
@@ -22,9 +22,9 @@
 
 @pytest.mark.parametrize("ndim", [2, 3])
 @pytest.mark.parametrize("mask_dtype", [
-    np.uint8, np.uint16, np.int64, bool])
+    np.uint8, np.uint16, np.int32, np.int64, bool])
 @pytest.mark.parametrize("spot_dtype", [
-    np.int64, np.float64])
+    np.float32, np.float64, np.int32, np.int64])
 def test_identify_objects_in_region(ndim, mask_dtype, spot_dtype):
     # simulate mask and coordinates
     mask = np.zeros((10, 10), dtype=mask_dtype)
@@ -49,9 +49,9 @@ def test_identify_objects_in_region(ndim, mask_dtype, spot_dtype):
 
 @pytest.mark.parametrize("ndim", [2, 3])
 @pytest.mark.parametrize("mask_dtype", [
-    np.uint8, np.uint16, np.int64, bool])
+    np.uint8, np.uint16, np.int32, np.int64, bool])
 @pytest.mark.parametrize("spot_dtype", [
-    np.int64, np.float64])
+    np.float32, np.float64, np.int32, np.int64])
 def test_remove_transcription_site(ndim, mask_dtype, spot_dtype):
     # simulate mask and coordinates
     nuc_mask = np.zeros((10, 10), dtype=mask_dtype)
diff --git a/bigfish/plot/__init__.py b/bigfish/plot/__init__.py
index 65b5f15..29c1039 100644
--- a/bigfish/plot/__init__.py
+++ b/bigfish/plot/__init__.py
@@ -11,10 +11,10 @@
 from .plot_images import plot_segmentation
 from .plot_images import plot_segmentation_boundary
 from .plot_images import plot_segmentation_diff
-
 from .plot_images import plot_detection
 from .plot_images import plot_reference_spot
 from .plot_images import plot_cell
+from .plot_images import plot_cell_coordinates
 
 from .plot_quality import plot_sharpness
 from .plot_quality import plot_elbow
@@ -33,7 +33,8 @@
     "plot_segmentation_diff",
     "plot_detection",
     "plot_reference_spot",
-    "plot_cell"]
+    "plot_cell",
+    "plot_cell_coordinates"]
 
 _quality = [
     "plot_sharpness",
diff --git a/bigfish/plot/plot_images.py b/bigfish/plot/plot_images.py
index 1088035..aaf6672 100644
--- a/bigfish/plot/plot_images.py
+++ b/bigfish/plot/plot_images.py
@@ -23,9 +23,19 @@
 
 # ### General plot ###
 
-def plot_yx(image, r=0, c=0, z=0, rescale=False, contrast=False,
-            title=None, framesize=(10, 10), remove_frame=True,
-            path_output=None, ext="png", show=True):
+def plot_yx(
+        image,
+        r=0,
+        c=0,
+        z=0,
+        rescale=False,
+        contrast=False,
+        title=None,
+        framesize=(10, 10),
+        remove_frame=True,
+        path_output=None,
+        ext="png",
+        show=True):
     """Plot the selected yx plan of the selected dimensions of an image.
 
     Parameters
@@ -62,7 +72,9 @@ def plot_yx(image, r=0, c=0, z=0, rescale=False, contrast=False,
     stack.check_array(
         image,
         ndim=[2, 3, 4, 5],
-        dtype=[np.uint8, np.uint16, np.int64, np.float32, np.float64, bool])
+        dtype=[np.uint8, np.uint16, np.int32, np.int64,
+               np.float32, np.float64,
+               bool])
     stack.check_parameter(
         r=int, c=int, z=int,
         rescale=bool,
@@ -112,9 +124,16 @@ def plot_yx(image, r=0, c=0, z=0, rescale=False, contrast=False,
         plt.close()
 
 
-def plot_images(images, rescale=False, contrast=False, titles=None,
-                framesize=(15, 10), remove_frame=True, path_output=None,
-                ext="png", show=True):
+def plot_images(
+        images,
+        rescale=False,
+        contrast=False,
+        titles=None,
+        framesize=(15, 10),
+        remove_frame=True,
+        path_output=None,
+        ext="png",
+        show=True):
     """Plot or subplot of 2-d images.
 
     Parameters
@@ -139,7 +158,6 @@ def plot_images(images, rescale=False, contrast=False, titles=None,
     show : bool, default=True
         Show the figure or not.
 
-
     """
     # enlist image if necessary
     if isinstance(images, np.ndarray):
@@ -160,10 +178,14 @@ def plot_images(images, rescale=False, contrast=False, titles=None,
         stack.check_array(
             image,
             ndim=2,
-            dtype=[np.uint8, np.uint16, np.int64,
+            dtype=[np.uint8, np.uint16, np.int32, np.int64,
                    np.float32, np.float64,
                    bool])
 
+    # enlist 'titles' if needed
+    if titles is not None and isinstance(titles, str):
+        titles = [titles]
+
     # we plot 3 images by row maximum
     nrow = int(np.ceil(len(images)/3))
     ncol = min(len(images), 3)
@@ -247,9 +269,17 @@ def plot_images(images, rescale=False, contrast=False, titles=None,
 
 # ### Segmentation plot ###
 
-def plot_segmentation(image, mask, rescale=False, contrast=False, title=None,
-                      framesize=(15, 10), remove_frame=True,
-                      path_output=None, ext="png", show=True):
+def plot_segmentation(
+        image,
+        mask,
+        rescale=False,
+        contrast=False,
+        title=None,
+        framesize=(15, 10),
+        remove_frame=True,
+        path_output=None,
+        ext="png",
+        show=True):
     """Plot result of a 2-d segmentation, with labelled instances if available.
 
     Parameters
@@ -281,11 +311,13 @@ def plot_segmentation(image, mask, rescale=False, contrast=False, title=None,
     stack.check_array(
         image,
         ndim=2,
-        dtype=[np.uint8, np.uint16, np.int64, np.float32, np.float64, bool])
+        dtype=[np.uint8, np.uint16, np.int32, np.int64,
+               np.float32, np.float64,
+               bool])
     stack.check_array(
         mask,
         ndim=2,
-        dtype=[np.uint8, np.uint16, np.int64, bool])
+        dtype=[np.uint8, np.uint16, np.int32, np.int64, bool])
     stack.check_parameter(
         rescale=bool,
         contrast=bool,
@@ -346,11 +378,19 @@ def plot_segmentation(image, mask, rescale=False, contrast=False, title=None,
         plt.close()
 
 
-def plot_segmentation_boundary(image, cell_label=None, nuc_label=None,
-                               boundary_size=1, rescale=False, contrast=False,
-                               title=None, framesize=(10, 10),
-                               remove_frame=True, path_output=None, ext="png",
-                               show=True):
+def plot_segmentation_boundary(
+        image,
+        cell_label=None,
+        nuc_label=None,
+        boundary_size=1,
+        rescale=False,
+        contrast=False,
+        title=None,
+        framesize=(10, 10),
+        remove_frame=True,
+        path_output=None,
+        ext="png",
+        show=True):
     """Plot the boundary of the segmented objects.
 
     Parameters
@@ -386,17 +426,19 @@ def plot_segmentation_boundary(image, cell_label=None, nuc_label=None,
     stack.check_array(
         image,
         ndim=2,
-        dtype=[np.uint8, np.uint16, np.int64, np.float32, np.float64, bool])
+        dtype=[np.uint8, np.uint16, np.int32, np.int64,
+               np.float32, np.float64,
+               bool])
     if cell_label is not None:
         stack.check_array(
             cell_label,
             ndim=2,
-            dtype=[np.uint8, np.uint16, np.int64, bool])
+            dtype=[np.uint8, np.uint16, np.int32, np.int64, bool])
     if nuc_label is not None:
         stack.check_array(
             nuc_label,
             ndim=2,
-            dtype=[np.uint8, np.uint16, np.int64, bool])
+            dtype=[np.uint8, np.uint16, np.int32, np.int64, bool])
     stack.check_parameter(
         rescale=bool,
         contrast=bool,
@@ -461,10 +503,18 @@ def plot_segmentation_boundary(image, cell_label=None, nuc_label=None,
         plt.close()
 
 
-def plot_segmentation_diff(image, mask_pred, mask_gt, rescale=False,
-                           contrast=False, title=None, framesize=(15, 10),
-                           remove_frame=True, path_output=None, ext="png",
-                           show=True):
+def plot_segmentation_diff(
+        image,
+        mask_pred,
+        mask_gt,
+        rescale=False,
+        contrast=False,
+        title=None,
+        framesize=(15, 10),
+        remove_frame=True,
+        path_output=None,
+        ext="png",
+        show=True):
     """Plot segmentation results along with ground truth to compare.
 
     Parameters
@@ -507,7 +557,9 @@ def plot_segmentation_diff(image, mask_pred, mask_gt, rescale=False,
     stack.check_array(
         image,
         ndim=2,
-        dtype=[np.uint8, np.uint16, np.int64, np.float32, np.float64, bool])
+        dtype=[np.uint8, np.uint16, np.int32, np.int64,
+               np.float32, np.float64,
+               bool])
     stack.check_array(
         mask_pred,
         ndim=2,
@@ -569,10 +621,22 @@ def plot_segmentation_diff(image, mask_pred, mask_gt, rescale=False,
 
 # ### Detection plot ###
 
-def plot_detection(image, spots, shape="circle", radius=3, color="red",
-                   linewidth=1, fill=False, rescale=False, contrast=False,
-                   title=None, framesize=(15, 10), remove_frame=True,
-                   path_output=None, ext="png", show=True):
+def plot_detection(
+        image,
+        spots,
+        shape="circle",
+        radius=3,
+        color="red",
+        linewidth=1,
+        fill=False,
+        rescale=False,
+        contrast=False,
+        title=None,
+        framesize=(15, 10),
+        remove_frame=True,
+        path_output=None,
+        ext="png",
+        show=True):
     """Plot detected spots and foci on a 2-d image.
 
     Parameters
@@ -626,7 +690,8 @@ def plot_detection(image, spots, shape="circle", radius=3, color="red",
     stack.check_array(
         image,
         ndim=2,
-        dtype=[np.uint8, np.uint16, np.int64, np.float32, np.float64])
+        dtype=[np.uint8, np.uint16, np.int32, np.int64,
+               np.float32, np.float64])
     stack.check_parameter(
         spots=(list, np.ndarray),
         shape=(list, str),
@@ -644,9 +709,15 @@ def plot_detection(image, spots, shape="circle", radius=3, color="red",
         show=bool)
     if isinstance(spots, list):
         for spots_ in spots:
-            stack.check_array(spots_, ndim=2, dtype=[np.int64, np.float64])
+            stack.check_array(
+                spots_,
+                ndim=2,
+                dtype=[np.float32, np.float64, np.int32, np.int64])
     else:
-        stack.check_array(spots, ndim=2, dtype=[np.int64, np.float64])
+        stack.check_array(
+            spots,
+            ndim=2,
+            dtype=[np.float32, np.float64, np.int32, np.int64])
 
     # enlist and format parameters
     if not isinstance(spots, list):
@@ -713,8 +784,8 @@ def plot_detection(image, spots, shape="circle", radius=3, color="red",
 
         # plot symbols
         for y, x in coordinates_2d:
-            x = _define_patch(x, y, shape[i], radius[i], color[i],
-                              linewidth[i], fill[i])
+            x = _define_patch(
+                x, y, shape[i], radius[i], color[i], linewidth[i], fill[i])
             ax[1].add_patch(x)
 
     # titles and frames
@@ -794,9 +865,16 @@ def _define_patch(x, y, shape, radius, color, linewidth, fill):
     return x
 
 
-def plot_reference_spot(reference_spot, rescale=False, contrast=False,
-                        title=None, framesize=(5, 5), remove_frame=True,
-                        path_output=None, ext="png", show=True):
+def plot_reference_spot(
+        reference_spot,
+        rescale=False,
+        contrast=False,
+        title=None,
+        framesize=(5, 5),
+        remove_frame=True,
+        path_output=None,
+        ext="png",
+        show=True):
     """Plot the selected yx plan of the selected dimensions of an image.
 
     Parameters
@@ -826,7 +904,8 @@ def plot_reference_spot(reference_spot, rescale=False, contrast=False,
     stack.check_array(
         reference_spot,
         ndim=[2,  3],
-        dtype=[np.uint8, np.uint16, np.int64, np.float32, np.float64])
+        dtype=[np.uint8, np.uint16, np.int32, np.int64,
+               np.float32, np.float64])
     stack.check_parameter(
         rescale=bool,
         contrast=bool,
@@ -873,11 +952,25 @@ def plot_reference_spot(reference_spot, rescale=False, contrast=False,
 
 # ### Individual cell plot ###
 
-def plot_cell(ndim, cell_coord=None, nuc_coord=None, rna_coord=None,
-              foci_coord=None, other_coord=None, image=None, cell_mask=None,
-              nuc_mask=None, boundary_size=1, title=None, remove_frame=True,
-              rescale=False, contrast=False, framesize=(15, 10),
-              path_output=None, ext="png", show=True):
+def plot_cell(
+        ndim,
+        cell_coord=None,
+        nuc_coord=None,
+        rna_coord=None,
+        foci_coord=None,
+        other_coord=None,
+        image=None,
+        cell_mask=None,
+        nuc_mask=None,
+        boundary_size=1,
+        title=None,
+        remove_frame=True,
+        rescale=False,
+        contrast=False,
+        framesize=(15, 10),
+        path_output=None,
+        ext="png",
+        show=True):
     """
     Plot image and coordinates extracted for a specific cell.
 
@@ -885,25 +978,25 @@ def plot_cell(ndim, cell_coord=None, nuc_coord=None, rna_coord=None,
     ----------
     ndim : {2, 3}
         Number of spatial dimensions to consider in the coordinates.
-    cell_coord : np.ndarray, np.int64, optional
+    cell_coord : np.ndarray, optional
         Coordinates of the cell border with shape (nb_points, 2). If None,
         coordinate representation of the cell is not shown.
-    nuc_coord : np.ndarray, np.int64, optional
+    nuc_coord : np.ndarray, optional
         Coordinates of the nucleus border with shape (nb_points, 2).
-    rna_coord : np.ndarray, np.int64, optional
+    rna_coord : np.ndarray, optional
         Coordinates of the detected spots with shape (nb_spots, 4) or
         (nb_spots, 3). One coordinate per dimension (zyx or yx dimensions)
         plus the index of the cluster assigned to the spot. If no cluster was
         assigned, value is -1. If only coordinates of spatial dimensions are
         available, only centroid of foci can be shown.
-    foci_coord : np.ndarray, np.int64, optional
+    foci_coord : np.ndarray, optional
         Array with shape (nb_foci, 5) or (nb_foci, 4). One coordinate per
         dimension for the foci centroid (zyx or yx dimensions), the number of
         spots detected in the foci and its index.
-    other_coord : np.ndarray, np.int64, optional
+    other_coord : np.ndarray, optional
         Coordinates of the detected elements with shape (nb_elements, 3) or
         (nb_elements, 2). One coordinate per dimension (zyx or yx dimensions).
-    image : np.ndarray, np.uint, optional
+    image : np.ndarray, optional
         Original image of the cell with shape (y, x). If None, original image
         of the cell is not shown.
     cell_mask : np.ndarray, optional
@@ -936,30 +1029,46 @@ def plot_cell(ndim, cell_coord=None, nuc_coord=None, rna_coord=None,
 
     # check parameters
     if cell_coord is not None:
-        stack.check_array(cell_coord, ndim=2, dtype=np.int64)
+        stack.check_array(
+            cell_coord,
+            ndim=2,
+            dtype=[np.float32, np.float64, np.int32, np.int64])
     if nuc_coord is not None:
-        stack.check_array(nuc_coord, ndim=2, dtype=np.int64)
+        stack.check_array(
+            nuc_coord,
+            ndim=2,
+            dtype=[np.float32, np.float64, np.int32, np.int64])
     if rna_coord is not None:
-        stack.check_array(rna_coord, ndim=2, dtype=np.int64)
+        stack.check_array(
+            rna_coord,
+            ndim=2,
+            dtype=[np.float32, np.float64, np.int32, np.int64])
     if foci_coord is not None:
-        stack.check_array(foci_coord, ndim=2, dtype=np.int64)
+        stack.check_array(
+            foci_coord,
+            ndim=2,
+            dtype=[np.float32, np.float64, np.int32, np.int64])
     if other_coord is not None:
-        stack.check_array(other_coord, ndim=2, dtype=np.int64)
+        stack.check_array(
+            other_coord,
+            ndim=2,
+            dtype=[np.float32, np.float64, np.int32, np.int64])
     if image is not None:
         stack.check_array(
             image,
             ndim=2,
-            dtype=[np.uint8, np.uint16, np.int64, np.float32, np.float64])
+            dtype=[np.uint8, np.uint16, np.int32, np.int64,
+                   np.float32, np.float64])
     if cell_mask is not None:
         stack.check_array(
             cell_mask,
             ndim=2,
-            dtype=[np.uint8, np.uint16, np.int64, bool])
+            dtype=[np.uint8, np.uint16, np.int32, np.int64, bool])
     if nuc_mask is not None:
         stack.check_array(
             nuc_mask,
             ndim=2,
-            dtype=[np.uint8, np.uint16, np.int64, bool])
+            dtype=[np.uint8, np.uint16, np.int32, np.int64, bool])
     stack.check_parameter(
         ndim=int,
         boundary_size=int,
@@ -1180,3 +1289,240 @@ def plot_cell(ndim, cell_coord=None, nuc_coord=None, rna_coord=None,
             path_output=path_output,
             ext=ext,
             show=show)
+
+
+def plot_cell_coordinates(
+        ndim,
+        cell_coord,
+        nuc_coord,
+        rna_coord,
+        titles=None,
+        remove_frame=True,
+        framesize=(10, 5),
+        path_output=None,
+        ext="png",
+        show=True):
+    """
+    Plot cell coordinates for one or several cells.
+
+    Parameters
+    ----------
+    ndim : {2, 3}
+        Number of spatial dimensions to consider in the coordinates.
+    cell_coord : np.ndarray or list
+        Coordinates or list of coordinates of the cell border with shape
+        (nb_points, 2).
+    nuc_coord : np.ndarray or list
+        Coordinates or list of coordinates of the nucleus border with shape
+        (nb_points, 2).
+    rna_coord : np.ndarray or list
+        Coordinates or list of coordinates of the detected spots with shape
+        (nb_spots, 3) or (nb_spots, 2). One coordinate per dimension (zyx or
+        yx dimensions).
+    titles : str or list, optional
+        Title or list of titles.
+    remove_frame : bool, default=True
+        Remove axes and frame.
+    framesize : tuple, default=(10, 5)
+        Size of the frame.
+    path_output : str, optional
+        Path to save the image (without extension).
+    ext : str or list, default='png'
+        Extension used to save the plot. If it is a list of strings, the plot
+        will be saved several times.
+    show : bool, default=True
+        Show the figure or not.
+
+    """
+    # enlist coordinates if necessary
+    if isinstance(cell_coord, np.ndarray):
+        cell_coord = [cell_coord]
+    if isinstance(nuc_coord, np.ndarray):
+        nuc_coord = [nuc_coord]
+    if isinstance(rna_coord, np.ndarray):
+        rna_coord = [rna_coord]
+
+    # check parameters
+    stack.check_parameter(
+        ndim=int,
+        titles=(str, list, type(None)),
+        remove_frame=bool,
+        framesize=tuple,
+        path_output=(str, type(None)),
+        ext=(str, list))
+
+    # check coordinates
+    for i in range(len(cell_coord)):
+        stack.check_array(
+            cell_coord[i],
+            ndim=2,
+            dtype=[np.float32, np.float64, np.int32, np.int64])
+        stack.check_array(
+            nuc_coord[i],
+            ndim=2,
+            dtype=[np.float32, np.float64, np.int32, np.int64])
+        stack.check_array(
+            rna_coord[i],
+            ndim=2,
+            dtype=[np.float32, np.float64, np.int32, np.int64])
+
+    # enlist 'titles' if needed
+    if titles is not None and isinstance(titles, str):
+        titles = [titles]
+
+    # we plot 3 images by row maximum
+    nrow = int(np.ceil(len(cell_coord) / 3))
+    ncol = min(len(cell_coord), 3)
+
+    # plot one image
+    if len(cell_coord) == 1:
+
+        # frame
+        if remove_frame:
+            fig = plt.figure(figsize=framesize, frameon=False)
+            ax = fig.add_axes([0, 0, 1, 1])
+            ax.axis('off')
+        else:
+            plt.figure(figsize=framesize)
+
+        # coordinate image
+        plt.plot(
+            cell_coord[0][:, 1],
+            cell_coord[0][:, 0],
+            c="black",
+            linewidth=2)
+        plt.plot(
+            nuc_coord[0][:, 1],
+            nuc_coord[0][:, 0],
+            c="steelblue",
+            linewidth=2)
+        plt.scatter(
+            rna_coord[0][:, ndim - 1],
+            rna_coord[0][:, ndim - 2],
+            s=25,
+            c="firebrick",
+            marker=".")
+
+        # titles and frames
+        _, _, min_y, max_y = plt.axis()
+        plt.ylim(max_y, min_y)
+        plt.use_sticky_edges = True
+        plt.margins(0.01, 0.01)
+        plt.axis('scaled')
+        if titles is not None:
+            plt.title(
+                titles[0],
+                fontweight="bold",
+                fontsize=10)
+        if not remove_frame:
+            plt.tight_layout()
+
+        # output
+        if path_output is not None:
+            save_plot(path_output, ext)
+        if show:
+            plt.show()
+        else:
+            plt.close()
+
+        return
+
+    # plot multiple images
+    fig, ax = plt.subplots(nrow, ncol, figsize=framesize)
+
+    # one row
+    if len(cell_coord) in [2, 3]:
+
+        # loop over instance coordinates
+        for i in range(len(cell_coord)):
+
+            # coordinate image
+            ax[i].plot(
+                cell_coord[i][:, 1],
+                cell_coord[i][:, 0],
+                c="black",
+                linewidth=2)
+            ax[i].plot(
+                nuc_coord[i][:, 1],
+                nuc_coord[i][:, 0],
+                c="steelblue",
+                linewidth=2)
+            ax[i].scatter(
+                rna_coord[i][:, ndim - 1],
+                rna_coord[i][:, ndim - 2],
+                s=25,
+                c="firebrick",
+                marker=".")
+
+            # titles and frames
+            _, _, min_y, max_y = ax[i].axis()
+            ax[i].set_ylim(max_y, min_y)
+            ax[i].use_sticky_edges = True
+            ax[i].margins(0.01, 0.01)
+            ax[i].axis('scaled')
+            if remove_frame:
+                ax[i].axis("off")
+            if titles is not None:
+                ax[i].set_title(
+                    titles[i],
+                    fontweight="bold",
+                    fontsize=10)
+
+    # several rows
+    else:
+
+        # we complete the row with empty frames
+        r = nrow * 3 - len(cell_coord)
+        cell_coord_completed = [cell_coord_ for cell_coord_ in cell_coord]
+        cell_coord_completed += [None] * r
+
+        # loop over instance coordinates
+        for i in range(len(cell_coord_completed)):
+            row = i // 3
+            col = i % 3
+
+            # empty subplot
+            if cell_coord_completed[i] is None:
+                ax[row, col].set_visible(False)
+                continue
+
+            # coordinate image
+            ax[row, col].plot(
+                cell_coord_completed[i][:, 1],
+                cell_coord_completed[i][:, 0],
+                c="black",
+                linewidth=2)
+            ax[row, col].plot(
+                nuc_coord[i][:, 1],
+                nuc_coord[i][:, 0],
+                c="steelblue",
+                linewidth=2)
+            ax[row, col].scatter(
+                rna_coord[i][:, ndim - 1],
+                rna_coord[i][:, ndim - 2],
+                s=25,
+                c="firebrick",
+                marker=".")
+
+            # titles and frames
+            _, _, min_y, max_y = ax[row, col].axis()
+            ax[row, col].set_ylim(max_y, min_y)
+            ax[row, col].use_sticky_edges = True
+            ax[row, col].margins(0.01, 0.01)
+            ax[row, col].axis('scaled')
+            if remove_frame:
+                ax[row, col].axis("off")
+            if titles is not None:
+                ax[row, col].set_title(
+                    titles[i],
+                    fontweight="bold",
+                    fontsize=10)
+
+    # output
+    plt.tight_layout()
+    if path_output is not None:
+        save_plot(path_output, ext)
+    if show:
+        plt.show()
+    else:
+        plt.close()
diff --git a/bigfish/plot/plot_quality.py b/bigfish/plot/plot_quality.py
index bfeb31e..8b45061 100644
--- a/bigfish/plot/plot_quality.py
+++ b/bigfish/plot/plot_quality.py
@@ -18,9 +18,17 @@
 
 # ### Focus - sharpness ###
 
-def plot_sharpness(focus_measures, labels=None, title=None, framesize=(5, 5),
-                   size_title=20, size_axes=15, size_legend=15,
-                   path_output=None, ext="png", show=True):
+def plot_sharpness(
+        focus_measures,
+        labels=None,
+        title=None,
+        framesize=(5, 5),
+        size_title=20,
+        size_axes=15,
+        size_legend=15,
+        path_output=None,
+        ext="png",
+        show=True):
     """Plot focus measures of a 3-d image, at the z-slice level.
 
     A measure of focus for each z-slice can be computed by averaging the
@@ -105,10 +113,20 @@ def plot_sharpness(focus_measures, labels=None, title=None, framesize=(5, 5),
 
 # ### Elbow plots ###
 
-def plot_elbow(images, voxel_size=None, spot_radius=None, log_kernel_size=None,
-               minimum_distance=None, title=None, framesize=(5, 5),
-               size_title=20, size_axes=15, size_legend=15, path_output=None,
-               ext="png", show=True):
+def plot_elbow(
+        images,
+        voxel_size=None,
+        spot_radius=None,
+        log_kernel_size=None,
+        minimum_distance=None,
+        title=None,
+        framesize=(5, 5),
+        size_title=20,
+        size_axes=15,
+        size_legend=15,
+        path_output=None,
+        ext="png",
+        show=True):
     """Plot the elbow curve that allows an automated spot detection.
 
     Parameters
@@ -208,10 +226,19 @@ def plot_elbow(images, voxel_size=None, spot_radius=None, log_kernel_size=None,
         plt.close()
 
 
-def plot_elbow_colocalized(spots_1, spots_2, voxel_size, threshold_max=None,
-                           title=None, framesize=(5, 5), size_title=20,
-                           size_axes=15, size_legend=15, path_output=None,
-                           ext="png", show=True):
+def plot_elbow_colocalized(
+        spots_1,
+        spots_2,
+        voxel_size,
+        threshold_max=None,
+        title=None,
+        framesize=(5, 5),
+        size_title=20,
+        size_axes=15,
+        size_legend=15,
+        path_output=None,
+        ext="png",
+        show=True):
     """Plot the elbow curve that allows an automated colocalized spot
     detection.
 
diff --git a/bigfish/segmentation/cell_segmentation.py b/bigfish/segmentation/cell_segmentation.py
index a28a441..35ab48d 100644
--- a/bigfish/segmentation/cell_segmentation.py
+++ b/bigfish/segmentation/cell_segmentation.py
@@ -44,8 +44,13 @@ def unet_distance_edge_double():
     return model
 
 
-def apply_unet_distance_double(model, nuc, cell, nuc_label, target_size=None,
-                               test_time_augmentation=False):
+def apply_unet_distance_double(
+        model,
+        nuc,
+        cell,
+        nuc_label,
+        target_size=None,
+        test_time_augmentation=False):
     """Segment cell with a pretrained model to predict distance map and  use
     it with a watershed algorithm.
 
@@ -203,8 +208,12 @@ def apply_unet_distance_double(model, nuc, cell, nuc_label, target_size=None,
     return cell_label_pred
 
 
-def from_distance_to_instances(label_x_nuc, label_2_cell, label_distance,
-                               nuc_3_classes=False, compute_nuc_label=False):
+def from_distance_to_instances(
+        label_x_nuc,
+        label_2_cell,
+        label_distance,
+        nuc_3_classes=False,
+        compute_nuc_label=False):
     """Extract instance labels from a distance map and a binary surface
     prediction with a watershed algorithm.
 
@@ -237,8 +246,8 @@ def from_distance_to_instances(label_x_nuc, label_2_cell, label_distance,
         nuc_3_classes=bool,
         compute_nuc_label=bool)
     stack.check_array(label_x_nuc, ndim=2, dtype=[np.float32, np.int64])
-    stack.check_array(label_2_cell, ndim=2, dtype=[np.float32])
-    stack.check_array(label_distance, ndim=2, dtype=[np.uint16])
+    stack.check_array(label_2_cell, ndim=2, dtype=np.float32)
+    stack.check_array(label_distance, ndim=2, dtype=np.uint16)
 
     # get nuclei labels
     if nuc_3_classes and compute_nuc_label:
@@ -451,7 +460,7 @@ def apply_watershed(watershed_relief, nuc_label, cell_mask):
     stack.check_array(
         watershed_relief,
         ndim=2,
-        dtype=[np.uint8, np.uint16, np.int64])
+        dtype=[np.uint8, np.uint16, np.int32, np.int64])
     stack.check_array(nuc_label, ndim=2, dtype=np.int64)
     stack.check_array(cell_mask, ndim=2, dtype=bool)
 
diff --git a/bigfish/segmentation/nuc_segmentation.py b/bigfish/segmentation/nuc_segmentation.py
index 9d459d9..0c87e98 100644
--- a/bigfish/segmentation/nuc_segmentation.py
+++ b/bigfish/segmentation/nuc_segmentation.py
@@ -38,8 +38,11 @@ def unet_3_classes_nuc():
     return model
 
 
-def apply_unet_3_classes(model, image, target_size=None,
-                         test_time_augmentation=False):
+def apply_unet_3_classes(
+        model,
+        image,
+        target_size=None,
+        test_time_augmentation=False):
     """Segment image with a 3-classes trained model.
 
     Parameters
@@ -168,7 +171,7 @@ def from_3_classes_to_instances(label_3_classes):
 
     """
     # check parameters
-    stack.check_array(label_3_classes, ndim=3, dtype=[np.float32])
+    stack.check_array(label_3_classes, ndim=3, dtype=np.float32)
 
     # get classes indices
     label_3_classes = np.argmax(label_3_classes, axis=-1)
diff --git a/bigfish/segmentation/postprocess.py b/bigfish/segmentation/postprocess.py
index 6050257..c74a7ad 100644
--- a/bigfish/segmentation/postprocess.py
+++ b/bigfish/segmentation/postprocess.py
@@ -15,8 +15,6 @@
 from skimage.morphology import remove_small_objects
 
 
-# TODO make functions compatible with different type of integers
-
 # ### Labelled images ###
 
 def label_instances(image_binary):
@@ -38,7 +36,7 @@ def label_instances(image_binary):
     stack.check_array(image_binary, ndim=[2, 3], dtype=bool)
 
     # label instances
-    image_label = label(image_binary)
+    image_label = label(image_binary).astype(np.int64)
 
     return image_label
 
@@ -84,8 +82,12 @@ def merge_labels(image_label_1, image_label_2):
 
 # ### Clean segmentation ###
 # TODO make it available for 3D images
-def clean_segmentation(image, small_object_size=None, fill_holes=False,
-                       smoothness=None, delimit_instance=False):
+def clean_segmentation(
+        image,
+        small_object_size=None,
+        fill_holes=False,
+        smoothness=None,
+        delimit_instance=False):
     """Clean segmentation results (binary masks or integer labels).
 
     Parameters
@@ -279,7 +281,7 @@ def remove_disjoint(image):
     stack.check_array(
         image,
         ndim=[2, 3],
-        dtype=[np.uint8, np.uint16, np.int64, bool])
+        dtype=[np.uint8, np.uint16, np.int32, np.int64, bool])
 
     # handle boolean array
     cast_to_bool = False
diff --git a/bigfish/segmentation/utils.py b/bigfish/segmentation/utils.py
index 4e768a3..1a3b2ec 100644
--- a/bigfish/segmentation/utils.py
+++ b/bigfish/segmentation/utils.py
@@ -13,8 +13,6 @@
 from skimage.measure import regionprops
 
 
-# TODO make functions compatible with different type of integers
-
 # ### Thresholding method ###
 
 def thresholding(image, threshold):
@@ -23,7 +21,7 @@ def thresholding(image, threshold):
 
     Parameters
     ----------
-    image : np.ndarray, np.uint
+    image : np.ndarray
         A 2-d image to segment with shape (y, x).
     threshold : int or float
         Pixel intensity threshold used to discriminate foreground from
@@ -36,7 +34,11 @@ def thresholding(image, threshold):
 
     """
     # check parameters
-    stack.check_array(image, ndim=2, dtype=[np.uint8, np.uint16])
+    stack.check_array(
+        image,
+        ndim=2,
+        dtype=[np.uint8, np.uint16, np.int32, np.int64,
+               np.float32, np.float64])
     stack.check_parameter(threshold=(float, int))
 
     # discriminate nuclei from background, applying a threshold.
@@ -68,7 +70,7 @@ def compute_mean_diameter(image_label):
     stack.check_array(
         image_label,
         ndim=2,
-        dtype=[np.uint8, np.uint16, np.int64])
+        dtype=[np.uint8, np.uint16, np.int32, np.int64])
 
     # compute properties of the segmented instances
     props = regionprops(image_label)
@@ -107,7 +109,7 @@ def compute_mean_convexity_ratio(image_label):
     stack.check_array(
         image_label,
         ndim=2,
-        dtype=[np.uint8, np.uint16, np.int64])
+        dtype=[np.uint8, np.uint16, np.int32, np.int64])
 
     # compute properties of the segmented instances
     props = regionprops(image_label)
@@ -145,7 +147,7 @@ def compute_surface_ratio(image_label):
     stack.check_array(
         image_label,
         ndim=2,
-        dtype=[np.uint8, np.uint16, np.int64])
+        dtype=[np.uint8, np.uint16, np.int32, np.int64])
 
     # compute surface ratio
     surface_instances = image_label > 0
@@ -173,7 +175,7 @@ def count_instances(image_label):
     stack.check_array(
         image_label,
         ndim=2,
-        dtype=[np.uint8, np.uint16, np.int64])
+        dtype=[np.uint8, np.uint16, np.int32, np.int64])
 
     indices = set(image_label.ravel())
     if 0 in indices:
diff --git a/bigfish/stack/filter.py b/bigfish/stack/filter.py
index c3d2c09..3b0d133 100644
--- a/bigfish/stack/filter.py
+++ b/bigfish/stack/filter.py
@@ -14,10 +14,18 @@
 from .preprocess import cast_img_uint8
 from .preprocess import cast_img_uint16
 
-from skimage.morphology.selem import square
-from skimage.morphology.selem import diamond
-from skimage.morphology.selem import rectangle
-from skimage.morphology.selem import disk
+import skimage
+from sklearn.utils.fixes import parse_version
+if parse_version(skimage.__version__) < parse_version("0.19.0"):
+    from skimage.morphology.selem import square
+    from skimage.morphology.selem import diamond
+    from skimage.morphology.selem import rectangle
+    from skimage.morphology.selem import disk
+else:
+    from skimage.morphology.footprints import square
+    from skimage.morphology.footprints import diamond
+    from skimage.morphology.footprints import rectangle
+    from skimage.morphology.footprints import disk
 from skimage.morphology import binary_dilation
 from skimage.morphology import dilation
 from skimage.morphology import binary_erosion
diff --git a/bigfish/stack/io.py b/bigfish/stack/io.py
index d932516..a1f5065 100644
--- a/bigfish/stack/io.py
+++ b/bigfish/stack/io.py
@@ -120,8 +120,12 @@ def read_array(path):
     return array
 
 
-def read_array_from_csv(path, dtype=None, delimiter=";", encoding="utf-8",
-                        skiprows=0):
+def read_array_from_csv(
+        path,
+        dtype=None,
+        delimiter=";",
+        encoding="utf-8",
+        skiprows=0):
     """Read a numpy array saved in a ``csv`` file.
 
     Parameters
@@ -381,7 +385,7 @@ def save_array(array, path):
                np.int8, np.int16, np.int32, np.int64,
                np.float16, np.float32, np.float64,
                bool],
-        ndim=[2, 3, 4, 5])
+        ndim=[1, 2, 3, 4, 5])
 
     # add extension if necessary
     if ".npy" not in path:
diff --git a/bigfish/stack/preprocess.py b/bigfish/stack/preprocess.py
index 83a30a3..0577e4e 100644
--- a/bigfish/stack/preprocess.py
+++ b/bigfish/stack/preprocess.py
@@ -43,7 +43,10 @@ def compute_image_standardization(image):
 
     """
     # check parameters
-    check_array(image, ndim=2, dtype=[np.uint8, np.uint16, np.float32])
+    check_array(
+        image,
+        ndim=2,
+        dtype=[np.uint8, np.uint16, np.float32, np.float64])
 
     # check image is in 2D
     if len(image.shape) != 2:
@@ -454,7 +457,10 @@ def resize_image(image, output_shape, method="bilinear"):
     """
     # check parameters
     check_parameter(output_shape=tuple, method=str)
-    check_array(image, ndim=[2, 3], dtype=[np.uint8, np.uint16, np.float32])
+    check_array(
+        image,
+        ndim=[2, 3],
+        dtype=[np.uint8, np.uint16, np.float32, np.float64])
 
     # resize image
     if method == "bilinear":
diff --git a/bigfish/stack/projection.py b/bigfish/stack/projection.py
index e2aefbd..947dc1f 100644
--- a/bigfish/stack/projection.py
+++ b/bigfish/stack/projection.py
@@ -32,7 +32,8 @@ def maximum_projection(image):
     check_array(
         image,
         ndim=3,
-        dtype=[np.uint8, np.uint16, np.float32, np.float64])
+        dtype=[np.uint8, np.uint16, np.int32, np.int64,
+               np.float32, np.float64])
 
     # project image along the z axis
     projected_image = image.max(axis=0)
@@ -61,7 +62,8 @@ def mean_projection(image, return_float=False):
     check_array(
         image,
         ndim=3,
-        dtype=[np.uint8, np.uint16, np.float32, np.float64])
+        dtype=[np.uint8, np.uint16, np.int32, np.int64,
+               np.float32, np.float64])
 
     # project image along the z axis
     if return_float:
@@ -91,7 +93,8 @@ def median_projection(image):
     check_array(
         image,
         ndim=3,
-        dtype=[np.uint8, np.uint16, np.float32, np.float64])
+        dtype=[np.uint8, np.uint16, np.int32, np.int64,
+               np.float32, np.float64])
 
     # project image along the z axis
     projected_image = np.median(image, axis=0)
@@ -100,8 +103,11 @@ def median_projection(image):
     return projected_image
 
 
-def focus_projection(image, proportion=0.75, neighborhood_size=7,
-                     method="median"):
+def focus_projection(
+        image,
+        proportion=0.75,
+        neighborhood_size=7,
+        method="median"):
     """Project the z-dimension of an image.
 
     Inspired from Samacoits Aubin's thesis (part 5.3, strategy 5). Compare to
@@ -137,7 +143,8 @@ def focus_projection(image, proportion=0.75, neighborhood_size=7,
     check_array(
         image,
         ndim=3,
-        dtype=[np.uint8, np.uint16, np.float32, np.float64])
+        dtype=[np.uint8, np.uint16, np.int32, np.int64,
+               np.float32, np.float64])
 
     # compute focus measure for each pixel
     focus = compute_focus(image, neighborhood_size)
@@ -247,7 +254,8 @@ def in_focus_selection(image, focus, proportion):
     check_array(
         image,
         ndim=3,
-        dtype=[np.uint8, np.uint16, np.float32, np.float64])
+        dtype=[np.uint8, np.uint16, np.int32, np.int64,
+               np.float32, np.float64])
 
     # select and keep best z-slices
     indices_to_keep = get_in_focus_indices(focus, proportion)
@@ -263,7 +271,7 @@ def get_in_focus_indices(focus, proportion):
 
     Parameters
     ----------
-    focus : np.ndarray, np.float64
+    focus : np.ndarray, np.float
         A 3-d tensor with a focus metric computed for each pixel of the
         original image. See :func:`bigfish.stack.compute_focus`.
     proportion : float or int
@@ -278,7 +286,7 @@ def get_in_focus_indices(focus, proportion):
     """
     # check parameters
     check_parameter(proportion=(float, int))
-    check_array(focus, ndim=3, dtype=np.float64)
+    check_array(focus, ndim=3, dtype=[np.float32, np.float64])
     if isinstance(proportion, float) and 0 <= proportion <= 1:
         n = int(focus.shape[0] * proportion)
     elif isinstance(proportion, int) and 0 <= proportion:
diff --git a/bigfish/stack/quality.py b/bigfish/stack/quality.py
index b47d4d7..53b82f6 100644
--- a/bigfish/stack/quality.py
+++ b/bigfish/stack/quality.py
@@ -50,7 +50,8 @@ def compute_focus(image, neighborhood_size=31):
     check_array(
         image,
         ndim=[2, 3],
-        dtype=[np.uint8, np.uint16, np.float32, np.float64])
+        dtype=[np.uint8, np.uint16, np.int32, np.int64,
+               np.float32, np.float64])
     check_parameter(neighborhood_size=(int, tuple, list))
     if (isinstance(neighborhood_size, (tuple, list))
             and len(neighborhood_size) != 2):
@@ -61,7 +62,7 @@ def compute_focus(image, neighborhood_size=31):
                          .format(len(neighborhood_size)))
 
     # cast image in float if necessary
-    if image.dtype in [np.uint8, np.uint16]:
+    if image.dtype in [np.uint8, np.uint16, np.int32, np.int64]:
         image_float = image.astype(np.float64)
     else:
         image_float = image
@@ -89,7 +90,7 @@ def _compute_focus_3d(image_3d, kernel_size):
 
     Returns
     -------
-    focus : np.ndarray, np.float64
+    focus : np.ndarray, np.float
         A 3-d tensor with the R(z, y, x) computed for each pixel of the
         original image.
 
@@ -101,7 +102,7 @@ def _compute_focus_3d(image_3d, kernel_size):
         focus.append(focus_2d)
 
     #  stack focus metrics
-    focus = np.stack(focus)
+    focus = np.stack(focus).astype(image_3d.dtype)
 
     return focus
 
@@ -120,7 +121,7 @@ def _compute_focus_2d(image_2d, kernel_size):
 
     Returns
     -------
-    focus : np.ndarray, np.float64
+    focus : np.ndarray, np.float
         A 2-d tensor with the R(y, x) computed for each pixel of the original
         image.
 
@@ -132,8 +133,8 @@ def _compute_focus_2d(image_2d, kernel_size):
         image_filtered_mean = mean_filter(image_2d, "rectangle", kernel_size)
 
     # compute focus metric
-    ratio_default_1 = np.ones_like(image_2d, dtype=np.float64)
-    ratio_default_2 = np.ones_like(image_filtered_mean, dtype=np.float64)
+    ratio_default_1 = np.ones_like(image_2d)
+    ratio_default_2 = np.ones_like(image_filtered_mean)
     ratio_1 = np.divide(
         image_2d, image_filtered_mean,
         out=ratio_default_1,
@@ -144,4 +145,7 @@ def _compute_focus_2d(image_2d, kernel_size):
         where=image_2d > 0)
     focus = np.where(image_2d >= image_filtered_mean, ratio_1, ratio_2)
 
+    # cast focus dtype (np.float32 or np.float64)
+    focus = focus.astype(image_2d.dtype)
+
     return focus
diff --git a/bigfish/stack/tests/test_io.py b/bigfish/stack/tests/test_io.py
index f4392c2..1230813 100644
--- a/bigfish/stack/tests/test_io.py
+++ b/bigfish/stack/tests/test_io.py
@@ -120,7 +120,7 @@ def test_dv(dtype):
 
 
 @pytest.mark.parametrize("shape", [
-    (8, 8), (8, 8, 8), (8, 8, 8, 8), (8, 8, 8, 8, 8)])
+    (8,), (8, 8), (8, 8, 8), (8, 8, 8, 8), (8, 8, 8, 8, 8)])
 @pytest.mark.parametrize("dtype", [
     np.uint8, np.uint16, np.uint32, np.uint64,
     np.int8, np.int16, np.int32, np.int64,
diff --git a/bigfish/stack/tests/test_projection.py b/bigfish/stack/tests/test_projection.py
index c4bebee..f9d7c34 100644
--- a/bigfish/stack/tests/test_projection.py
+++ b/bigfish/stack/tests/test_projection.py
@@ -80,7 +80,7 @@
 # ### test 2-d projection ###
 
 @pytest.mark.parametrize("dtype", [
-    np.uint8, np.uint16, np.float32, np.float64])
+    np.uint8, np.uint16, np.int32, np.int64, np.float32, np.float64])
 def test_maximum_projection(dtype):
     x = x_3d.astype(dtype)
     expected_y = np.array([[3, 0, 0, 0, 0],
@@ -94,7 +94,7 @@ def test_maximum_projection(dtype):
 
 
 @pytest.mark.parametrize("dtype", [
-    np.uint8, np.uint16, np.float32, np.float64])
+    np.uint8, np.uint16, np.int32, np.int64, np.float32, np.float64])
 def test_mean_projection(dtype):
     # 'return_float' == False
     x = x_3d.astype(dtype)
@@ -119,7 +119,7 @@ def test_mean_projection(dtype):
 
 
 @pytest.mark.parametrize("dtype", [
-    np.uint8, np.uint16, np.float32, np.float64])
+    np.uint8, np.uint16, np.int32, np.int64, np.float32, np.float64])
 def test_median_projection(dtype):
     x = x_3d.astype(dtype)
     expected_y = np.array([[2, 0, 0, 0, 0],
@@ -215,7 +215,7 @@ def test_get_in_focus_indices():
 
 
 @pytest.mark.parametrize("dtype", [
-    np.uint8, np.uint16, np.float32, np.float64])
+    np.uint8, np.uint16, np.int32, np.int64, np.float32, np.float64])
 def test_in_focus_selection(dtype):
     x = x_3d_out_focus.astype(dtype)
     expected_y = x_3d.astype(dtype)
@@ -226,7 +226,7 @@ def test_in_focus_selection(dtype):
 
 
 @pytest.mark.parametrize("dtype", [
-    np.uint8, np.uint16, np.float32, np.float64])
+    np.uint8, np.uint16, np.int32, np.int64, np.float32, np.float64])
 def test_focus_projection(dtype):
     x = x_3d_out_focus.astype(dtype)
 
diff --git a/bigfish/stack/tests/test_quality.py b/bigfish/stack/tests/test_quality.py
index fe57efd..068a67d 100644
--- a/bigfish/stack/tests/test_quality.py
+++ b/bigfish/stack/tests/test_quality.py
@@ -52,7 +52,7 @@
 
 
 @pytest.mark.parametrize("dtype", [
-    np.uint8, np.uint16, np.float32, np.float64])
+    np.uint8, np.uint16, np.int32, np.int64, np.float32, np.float64])
 @pytest.mark.parametrize("neighborhood_size", [
     31, (5, 11), [11, 5]])
 def test_compute_focus_format(dtype, neighborhood_size):
@@ -60,7 +60,10 @@ def test_compute_focus_format(dtype, neighborhood_size):
     image = np.random.normal(loc=0.0, scale=1.0, size=10000)
     image = np.reshape(image, (100, 100)).astype(dtype)
     focus = stack.compute_focus(image, neighborhood_size=neighborhood_size)
-    assert focus.dtype == np.float64
+    if dtype == np.float32:
+        assert focus.dtype == np.float32
+    else:
+        assert focus.dtype == np.float64
     assert focus.shape == image.shape
     assert focus.min() >= 1
 
@@ -68,7 +71,10 @@ def test_compute_focus_format(dtype, neighborhood_size):
     image = np.random.normal(loc=0.0, scale=1.0, size=100000)
     image = np.reshape(image, (10, 100, 100)).astype(dtype)
     focus = stack.compute_focus(image, neighborhood_size=neighborhood_size)
-    assert focus.dtype == np.float64
+    if dtype == np.float32:
+        assert focus.dtype == np.float32
+    else:
+        assert focus.dtype == np.float64
     assert focus.shape == image.shape
     assert focus.min() >= 1
 
diff --git a/bigfish/stack/utils.py b/bigfish/stack/utils.py
index dcf2851..4a6b8e9 100644
--- a/bigfish/stack/utils.py
+++ b/bigfish/stack/utils.py
@@ -160,7 +160,6 @@ def _check_dtype_array(array, dtype):
     if isinstance(dtype, type):
         dtype = [dtype]
 
-    # TODO simplify
     # check the dtype of the array
     error = True
     for dtype_expected in dtype:
diff --git a/docs/README.md b/docs/README.md
new file mode 100644
index 0000000..288923e
--- /dev/null
+++ b/docs/README.md
@@ -0,0 +1,5 @@
+# Documentation for Big-FISH
+
+[![Documentation Status](https://readthedocs.org/projects/big-fish/badge/?version=stable)](https://big-fish.readthedocs.io/en/latest/?badge=stable)
+
+This folder gathers the materials for the online [documentation](https://big-fish.readthedocs.io/en/stable/).
diff --git a/docs/source/plot/plot_coordinate.rst b/docs/source/plot/plot_coordinate.rst
index 86de745..867e041 100644
--- a/docs/source/plot/plot_coordinate.rst
+++ b/docs/source/plot/plot_coordinate.rst
@@ -7,4 +7,8 @@ Single-cell plot
 
 Function to visualize coordinate representation of individual cells:
 
+* :func:`bigfish.plot.plot_cell`
+* :func:`bigfish.plot.plot_cell_coordinates`
+
 .. autofunction:: plot_cell
+.. autofunction:: plot_cell_coordinates
diff --git a/docs/source/utils/utils.rst b/docs/source/utils/utils.rst
index e9a4d57..41f38e4 100644
--- a/docs/source/utils/utils.rst
+++ b/docs/source/utils/utils.rst
@@ -72,6 +72,17 @@ Convert pixels and nanometers
 
 ------------
 
+Extract a spot image
+====================
+
+* :func:`bigfish.detection.get_spot_volume`
+* :func:`bigfish.detection.get_spot_surface`
+
+.. autofunction:: get_spot_volume
+.. autofunction:: get_spot_surface
+
+------------
+
 Format and save plots
 =====================
 
diff --git a/requirements.txt b/requirements.txt
index 74c7fa7..3e97e58 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,5 +1,5 @@
 numpy >= 1.16.0
-scikit-learn >= 0.21.0
+scikit-learn >= 0.24.0
 scikit-image >= 0.14.2
 scipy >= 1.4.1
 matplotlib >= 3.0.2