diff --git a/scripts/functions_for_rasterio.py b/scripts/functions_for_rasterio.py
index 5ed85cb8f5cf2ada01f1ea7d49484253582e94cb..37199f1c077894ad14d10f92d368f267c32e3010 100644
--- a/scripts/functions_for_rasterio.py
+++ b/scripts/functions_for_rasterio.py
@@ -21,8 +21,10 @@ def rcs_pansharpen(xs, pan):
     print(0)
     xs_zoom = ndimage.zoom(xs, zoom=(1, 4, 4), order=0)
     print(1)
+    print(xs_zoom.shape)
     pan_low_pass = ndimage.gaussian_filter(pan, sigma=1)
     print(2)
+    print(pan_low_pass.shape)
 
     pansharpened = xs_zoom * (pan / pan_low_pass)
     print(3)
@@ -86,14 +88,15 @@ def pansharpen(m, pan, method='browley', W=0.1, all_data=False):
     image = None
 
     if method == 'simple_browley':
-        all_in = R+G+B
-        prod = np.multiply(all_in, pan)
+        all_in = R+G+B+I
+        prod = np.multiply(all_in, np.squeeze(pan))
 
-        r = np.multiply(R, pan/all_in)[:, :, np.newaxis]
-        g = np.multiply(G, pan/all_in)[:, :, np.newaxis]
-        b = np.multiply(B, pan/all_in)[:, :, np.newaxis]
+        r = np.multiply(R, np.squeeze(pan)/all_in)[:, :, np.newaxis]
+        g = np.multiply(G, np.squeeze(pan)/all_in)[:, :, np.newaxis]
+        b = np.multiply(B, np.squeeze(pan)/all_in)[:, :, np.newaxis]
+        i = np.multiply(I, np.squeeze(pan)/all_in)[:, :, np.newaxis]
 
-        image = np.concatenate([r,g,b], axis=2)
+        image = np.concatenate([r,g,b,i], axis=2)
 
     if method == 'sample_mean':
         r = 0.5 * (R + pan)[:, :, np.newaxis]
@@ -103,7 +106,7 @@ def pansharpen(m, pan, method='browley', W=0.1, all_data=False):
         image = np.concatenate([r,g,b], axis=2)
 
     if method == 'esri':
-        ADJ = pan-rgbn_scaled.mean(axis=2)
+        ADJ = np.squeeze(pan)-rgbn_scaled.mean(axis=2)
         r = (R + ADJ)[:, :, np.newaxis]
         g = (G + ADJ)[:, :, np.newaxis]
         b = (B + ADJ)[:, :, np.newaxis]
diff --git a/scripts/pyotb_usecase.py b/scripts/pyotb_usecase.py
new file mode 100644
index 0000000000000000000000000000000000000000..db171fad796f83e859d850472affc300985a2900
--- /dev/null
+++ b/scripts/pyotb_usecase.py
@@ -0,0 +1,30 @@
+import pyotb
+import rasterio
+from scipy import ndimage
+
+# Pansharpening
+pxs = pyotb.BundleToPerfectSensor(inp='pan.vrt', inxs='xs.tif') 
+
+# Computing NDVI
+ndvi = (pxs[:, :, -1] - pxs[:, :, 0]) / (pxs[:, :, -1] + pxs[:, :, 0])
+
+# Computing a boolean raster, describing the bare soils i.e. pixels without vegetation
+bare_soils = (ndvi < 0.3)
+
+# Creating a boolean cloud mask from the GML vector
+cloud_mask = pyotb.Rasterization('cloud_mask.GML', im=bare_soils)
+
+# Masking clouds (i.e. assigning to 0) on the result
+bare_soils_masked = pyotb.where(cloud_mask == 1, 0, bare_soils)
+
+# Exporting to rasterio
+bare_soils_arr, profile = bare_soils_masked.to_rasterio()  # warning, this loads the whole array in memory!
+
+# Gathering the contiguous groups of pixels with an integer index
+labeled_bare_soils, _ = ndimage.label(bare_soils_arr)
+
+# Writing the result to disk
+with rasterio.open('labeled_bare_soil.tif', 'w', **profile) as f:
+    f.write(labeled_bare_soils)
+
+
diff --git a/scripts/with_pyotb.py b/scripts/with_pyotb.py
index 8d99260c83b292310f5b1a6fbcc1c63c8f71dcd8..6f2582f07fdd9bee6b20df3fa892331b016d5c68 100644
--- a/scripts/with_pyotb.py
+++ b/scripts/with_pyotb.py
@@ -1,27 +1,8 @@
 import pyotb
-import rasterio
-from scipy import ndimage
-
-xs_path = 'xs.tif'
-pan_path = 'pan.vrt'
-cloud_path = 'cloud_mask.GML'
-
-pxs = pyotb.BundleToPerfectSensor(inp='pan.vrt', inxs='xs.tif') 
-
-ndvi = (pxs[:, :, -1] - pxs[:, :, 0]) / (pxs[:, :, -1] + pxs[:, :, 0])
-
-bare_soils = (ndvi < 0.3)
-
-cloud_mask = pyotb.Rasterization('cloud_mask.GML', im=bare_soils)
-
-bare_soils_masked = pyotb.where(cloud_mask == 1, 0, bare_soils)
-
-bare_soils_arr, profile = bare_soils_masked.to_rasterio()  # warning, this loads the whole array in memory!
-
-labeled_bare_soils, _ = ndimage.label(bare_soils_arr)
-
-# Writing the result to disk
-with rasterio.open('labeled_bare_soil.tif', 'w', **profile) as f:
-    f.write(labeled_bare_soils)
+import time
 
+start = time.time()
+pxs = pyotb.BundleToPerfectSensor(inp='pan.vrt', inxs='xs.tif', ram=190000) 
+pxs.write('pxs_pyotb.tif', pixel_type='uint16')
 
+print('ELAPSED TIME:', time.time() - start)
diff --git a/scripts/with_rasterio.py b/scripts/with_rasterio.py
index 30bf471b1ead83068ec647bde57434a4522211e7..320d949a10e2b5ca3b025ef4f5a3d1265fdde318 100644
--- a/scripts/with_rasterio.py
+++ b/scripts/with_rasterio.py
@@ -1,5 +1,11 @@
 import rasterio
 from functions_for_rasterio import rcs_pansharpen, pansharpen, dummy_pansharpen
+from scipy import ndimage
+
+
+import time
+
+start = time.time()
 
 xs_path = 'xs.tif'
 pan_path = 'pan.vrt'
@@ -17,10 +23,24 @@ with rasterio.open(pan_path) as pan_src:
 print(xs.shape)
 print(pan.shape)
 
+""" Just to try the gaussian filter
+pan_low_pass = ndimage.gaussian_filter(pan, sigma=1)
+
+profile.update(driver='GTiff')
+with rasterio.open('pan_gaussian.tif', 'w', **profile) as out_ds:
+    out_ds.write(pan_low_pass)
+"""
+
+print('Read time', time.time() - start)
+
 # Pansharpening
-pxs = pansharpen(xs[:, :1000,:1000], pan[:, :4000,:4000])
+pxs = rcs_pansharpen(xs, pan)
+
+print('After pan sharpen', time.time() - start)
 
 # Writing pansharpened output
 profile.update(count=4, height=pxs.shape[1], width=pxs.shape[2], driver='GTiff')
-with rasterio.open('/tmp/pansharpened_rasterio.tif', 'w', **profile) as out_ds:
+with rasterio.open('pansharpened_rasterio_rcs.tif', 'w', **profile) as out_ds:
     out_ds.write(pxs)
+
+print('ELAPSED TIME', time.time() - start)