Repository: fish-quant/big-fish (language: python)

Issue:
v.0.6.1 thresholding type error
Using big-fish 0.6.1 and running the following code from your example notebook

```
spots, threshold = detection.detect_spots(
    images=rna, 
    return_threshold=True, 
    voxel_size=(300, 103, 103),  # in nanometer (one value per dimension zyx)
    spot_radius=(300, 150, 150))  # in nanometer (one value per dimension zyx)
```
Generates the following error:

```
TypeError                                 Traceback (most recent call last)

File ~\anaconda3\envs\fish\lib\site-packages\bigfish\detection\spot_detection.py:516, in automated_threshold_setting(image, mask_local_max)
    514 # select threshold where the break of the distribution is located
    515 if count_spots.size > 0:
--> 516     optimal_threshold, _, _ = get_breaking_point(thresholds, count_spots)
    518 # case where no spots were detected
    519 else:
    520     optimal_threshold = None

File ~\anaconda3\envs\fish\lib\site-packages\bigfish\detection\utils.py:668, in get_breaking_point(x, y)
    645 """Select the x-axis value where a L-curve has a kink.
    646 
    647 Assuming a L-curve from A to B, the 'breaking_point' is the more distant
   (...)
    665 
    666 """
    667 # check parameters
--> 668 stack.check_array(x, ndim=1, dtype=[np.float64, np.int64])
    669 stack.check_array(y, ndim=1, dtype=[np.float64, np.int64])
    671 # select threshold where curve break

File ~\anaconda3\envs\fish\lib\site-packages\bigfish\stack\utils.py:135, in check_array(array, ndim, dtype, allow_nan)
    133 # check the dtype
    134 if dtype is not None:
--> 135     _check_dtype_array(array, dtype)
    137 # check the number of dimension
    138 if ndim is not None:

File ~\anaconda3\envs\fish\lib\site-packages\bigfish\stack\utils.py:172, in _check_dtype_array(array, dtype)
    169         break
    171 if error:
--> 172     raise TypeError("{0} is not supported yet. Use one of those dtypes "
    173                     "instead: {1}.".format(array.dtype, dtype))

TypeError: int32 is not supported yet. Use one of those dtypes instead: [<class 'numpy.float64'>, <class 'numpy.int64'>].

```
The error was suppressed by changing line 546 of spot_detection.py

from: `thresholds = np.array(thresholds)`
to: `thresholds = np.array(thresholds, dtype=np.int64)`




Discussion / hints:
Hi @owigg, do you use windows ? If yes, the error should be related to this [stackoverflow post](https://stackoverflow.com/questions/36278590/numpy-array-dtype-is-coming-as-int32-by-default-in-a-windows-10-64-bit-machine). By default your numpy arrays are in 32bit instead of 64bit. 

I will add a fix to solve this error whatever platform you use, thanks!
See 8fe67e0d16d6929a6edf4d27e9d1924ead481b14.
Hi @Henley13 ,

I was using the code on windows and I think I installed the last update but I still get the following error.
![image](https://user-images.githubusercontent.com/13469801/161114315-16ce51b6-ea66-4219-b5ca-d062d9a0440d.png)
 
Thanks for the fix. Yes, I am using Windows 64 bit.
Hi @luisub,

Since you work on windows, when you use numpy, arrays are defined in 32bit by default, instead of 64bit... 

See this [stackoverflow post](https://stackoverflow.com/questions/36278590/numpy-array-dtype-is-coming-as-int32-by-default-in-a-windows-10-64-bit-machine).

With some functions and intermediate functions, I am a bit too strict concerning the dtypes allowed, I will need to loose it for the next release (v0.6.2 I guess).

I will let this issue opened until I merge a PR solving it.
Hi, @Henley13 and @owigg. Thanks a lot for the feedback. I usually use Linux, and it has been working great! 

But we have undergrad students in the group who want to run codes on their personal Windows computers, and they are getting the error.
Hello Arthur :-),  same problem here on a windows machine with version 0.6.1.

happy to test the new version once released. 

Fix patch:
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/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..8b60cf2 100644
--- a/bigfish/detection/spot_modeling.py
+++ b/bigfish/detection/spot_modeling.py
@@ -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]
@@ -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]
diff --git a/bigfish/detection/utils.py b/bigfish/detection/utils.py
index 2a77b10..1671e0d 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),
@@ -336,11 +342,11 @@ def _get_spot_volume(image, spot_z, spot_y, spot_x, radius_z, radius_yx):
     ----------
     image : np.ndarray
         Image with shape (z, y, x).
-    spot_z : np.int64
+    spot_z : scalar
         Coordinate of the detected spot along the z axis.
-    spot_y : np.int64
+    spot_y : scalar
         Coordinate of the detected spot along the y axis.
-    spot_x : np.int64
+    spot_x : scalar
         Coordinate of the detected spot along the x axis.
     radius_z : int
         Radius in pixel of the detected spot, along the z axis.
@@ -445,9 +451,9 @@ def _get_spot_surface(image, spot_y, spot_x, radius_yx):
     ----------
     image : np.ndarray
         Image with shape (y, x).
-    spot_y : np.int64
+    spot_y : scalar
         Coordinate of the detected spot along the y axis.
-    spot_x : np.int64
+    spot_x : scalar
         Coordinate of the detected spot along the x axis.
     radius_yx : int
         Radius in pixel of the detected spot, on the yx plan.
@@ -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))
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..e78db5b 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,14 +351,16 @@ 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:
         for key in others_image:
@@ -352,7 +370,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 +539,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 +564,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 +589,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
@@ -734,7 +773,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 +782,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 +864,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 +890,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 +922,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 +930,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 +957,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 +975,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 +989,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 +1027,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 +1044,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 +1062,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 +1114,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/plot/plot_images.py b/bigfish/plot/plot_images.py
index 452b7e2..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,7 +178,7 @@ 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])
 
@@ -251,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
@@ -285,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,
@@ -350,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
@@ -390,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,
@@ -465,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
@@ -511,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,
@@ -573,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
@@ -630,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),
@@ -648,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):
@@ -717,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
@@ -798,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
@@ -830,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,
@@ -877,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.
 
@@ -889,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
@@ -939,32 +1028,47 @@ def plot_cell(ndim, cell_coord=None, nuc_coord=None, rna_coord=None,
         return
 
     # check parameters
-    # TODO allow np.float
     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,
@@ -1187,9 +1291,17 @@ def plot_cell(ndim, cell_coord=None, nuc_coord=None, rna_coord=None,
             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):
+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.
 
@@ -1241,9 +1353,18 @@ def plot_cell_coordinates(ndim, cell_coord, nuc_coord, rna_coord, titles=None,
 
     # check coordinates
     for i in range(len(cell_coord)):
-        stack.check_array(cell_coord[i], ndim=2, dtype=[np.int64, np.float64])
-        stack.check_array(nuc_coord[i], ndim=2, dtype=[np.int64, np.float64])
-        stack.check_array(rna_coord[i], ndim=2, dtype=[np.int64, np.float64])
+        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):
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 eaa7651..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
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/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:


Test patch:
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/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