WIP: solving holes problem with light lsgrm

OBIA/segmentation.py

@@ -182,7 +182,9 @@ def lsgrm_light(input_image, params : LSGRMParams, out_seg, n_proc=None, memory=
 
     print('[INFO] Using a layout of {} x {} tiles.'.format(n_tiles_x,n_tiles_y))
     print('[INFO] Nominal tile size w/margin: {} x {} pixels'.format(nominal_tw+2*params.margin,
-                                                                     nominal_th+2*params.margin))
+                                                                         nominal_th+2*params.margin))
+    
+    cumul_label, nholes = None, None
 
     if (n_tiles_x == 1 and n_tiles_y == 1):
         # Fallback to classical GRM
@@ -212,29 +214,50 @@ def lsgrm_light(input_image, params : LSGRMParams, out_seg, n_proc=None, memory=
 
         print('[INFO] Output mode {}.'.format(mode))
         if mode == 'vrt':
+            # Attention: cannot fill holes here!
             vrtopt = gdal.BuildVRTOptions(separate=False, srcNodata=0, VRTNodata=0)
             gdal.BuildVRT(out_seg, [x[0] for x in out_tiles], options=vrtopt)
         elif mode in ['raster', 'vector']:
             mos = otb.Registry.CreateApplication('Mosaic')
             mos.SetParameterStringList('il', [x[0] for x in out_tiles])
             mos.SetParameterInt('nodata', 0)
+            mos.Execute()
+
+            bm = otb.Registry.CreateApplication('BandMath')
+            bm.AddImageToParameterInputImageList('il', mos.GetParameterOutputImage('out'))
+            bm.AddImageToParameterInputImageList('il', in_img.GetParameterOutputImage('out'))
+            bm.SetParameterString('exp', 'im1b1==0 && im2b1!=0')
+            bm.Execute()
+
+            holes = bm.ExportImage('out')
+            holes['array'], nholes = label(holes['array'], background=0, connectivity=1, return_num=True)
+            holes['array'][holes['array']!=0] += cumul_label
+            holes['array'] = np.ascontiguousarray(holes['array']).astype(float)
+
+            # assuming nodata for input image always 0 here
+            holes_img = otb.Registry.CreateApplication('ExtractROI')
+            holes_img.ImportImage('in', holes)
+            holes_img.Execute()
+
+            out_img = otb.Registry.CreateApplication('BandMath')
+            out_img.AddImageToParameterInputImageList('il', mos.GetParameterOutputImage('out'))
+            out_img.AddImageToParameterInputImageList('il', holes_img.GetParameterOutputImage('out'))
+            out_img.SetParameterString('exp', 'im1b1+im2b1')
+
             if mode == 'raster':
-                mos.SetParameterOutputImagePixelType('out', otb.ImagePixelType_uint32)
-                mos.SetParameterString('out', out_seg)
-                mos.ExecuteAndWriteOutput()
+                out_img.SetParameterOutputImagePixelType('out', otb.ImagePixelType_uint32)
+                out_img.SetParameterString('out', out_seg)
+                out_img.ExecuteAndWriteOutput()
             elif mode == 'vector':
-                mos.Execute()
+                out_img.Execute()
                 vec = otb.Registry.CreateApplication('SimpleVectorization')
-                vec.SetParameterInputImage('in', mos.GetParameterOutputImage('out'))
+                vec.SetParameterInputImage('in', out_img.GetParameterOutputImage('out'))
                 vec.SetParameterString('out', out_seg)
                 vec.ExecuteAndWriteOutput()
             if remove_tiles_if_useless:
                 [os.remove(x[0]) for x in out_tiles]
 
-    return out_seg
-
-
-
+    return out_seg, cumul_label + nholes
 
 def lsgrm(input_image, params : LSGRMParams, out_seg, n_proc=None, memory=None, roi=None, remove_graph=True, force_parallel=False):
     # Check output file type