diff --git a/OBIA/OBIABase.py b/OBIA/OBIABase.py
new file mode 100644
index 0000000000000000000000000000000000000000..eead7fb425d7155778d0f9933a05e32f438574d7
--- /dev/null
+++ b/OBIA/OBIABase.py
@@ -0,0 +1,122 @@
+import os.path
+import sys
+import otbApplication as otb
+import gdal
+import numpy as np
+import pandas as pd
+from enum import Enum
+from Common.otb_numpy_proc import to_otb_pipeline
+from skimage.measure import regionprops
+
+class OStats(Enum):
+ MEAN = 'mean'
+ STD = 'std'
+ MIN = 'min'
+ MAX = 'max'
+
+def intensity_std(region, intensities):
+ # note the ddof arg to get the sample var if you so desire!
+ return np.std(intensities[region], ddof=1)
+
+class OBIABase:
+
+ def __init__(self, object_layer, nominal_tile_size=5000, background_label=0):
+ self.object_layer = object_layer
+ self.background_label = background_label
+ if type(nominal_tile_size) is tuple:
+ self.nominal_tile_size = nominal_tile_size
+ elif type(nominal_tile_size) is int:
+ self.nominal_tile_size = (nominal_tile_size, nominal_tile_size)
+ else:
+ sys.exit('Invalid nominal tile size')
+ self.tiles, self.obj_to_tile = self.generate_tile_info(self.nominal_tile_size, self.background_label)
+ # Init database
+ self.reset_db()
+
+ def reset_db(self):
+ self.database = pd.DataFrame(index=self.obj_to_tile.keys())
+
+ def generate_tile_info(self, nominal_tile_size, background_label):
+ in_seg = to_otb_pipeline(self.object_layer)
+ full = in_seg.GetImageAsNumpyArray('out')
+ rp = regionprops(full.astype(np.uint32))
+
+ W, H = in_seg.GetImageSize('out')
+ tx, ty = int(W / nominal_tile_size[0]) + 1, int(H / nominal_tile_size[1]) + 1
+
+ obj_to_tile = {}
+ tiles = dict.fromkeys(range(tx * ty))
+ for i in range(tx * ty):
+ tiles[i] = [np.inf, np.inf, 0, 0]
+
+ for o in rp:
+ if o.label != background_label:
+ ix, iy = int(o.bbox[1] / nominal_tile_size[0]), int(o.bbox[0] / nominal_tile_size[1])
+ idx = ix * ty + iy
+ obj_to_tile[o.label] = idx
+ tiles[idx][0] = min(o.bbox[1], tiles[idx][0])
+ tiles[idx][1] = min(o.bbox[0], tiles[idx][1])
+ tiles[idx][2] = max(o.bbox[3], tiles[idx][2])
+ tiles[idx][3] = max(o.bbox[2], tiles[idx][3])
+
+ tiles = {k: [v[0],v[1],v[2]-v[0],v[3]-v[1]] for k,v in tiles.items() if not np.isinf(v[0])}
+
+ return tiles, obj_to_tile
+
+ def compute_stats_on_tile(self, tile_num, input_image, stats=[OStats.MEAN]):
+ r = otb.itkRegion()
+ r['index'][0], r['index'][1] = self.tiles[tile_num][0], self.tiles[tile_num][1]
+ r['size'][0], r['size'][1] = self.tiles[tile_num][2], self.tiles[tile_num][3]
+
+ obj = to_otb_pipeline(self.object_layer)
+ obj.PropagateRequestedRegion('out', r)
+ clip_obj = obj.ExportImage('out')
+ tile_obj = np.squeeze(clip_obj['array']).astype(np.uint32)
+
+ si = otb.Registry.CreateApplication('Superimpose')
+ si.SetParameterString('inm', input_image)
+ si.SetParameterString('inr', self.object_layer)
+ si.SetParameterString('interpolator', 'nn')
+ si.Execute()
+ si.PropagateRequestedRegion('out', r)
+ clip_img = si.ExportImage('out')
+ tile_img = np.squeeze(clip_img['array'])
+
+ if OStats.STD in stats:
+ op = regionprops(tile_obj, intensity_image=tile_img, extra_properties=[intensity_std])
+ else:
+ op = regionprops(tile_obj, intensity_image=tile_img)
+
+ out_stats = {}
+ for o in op:
+ if self.obj_to_tile[o.label] == tile_num:
+ out_stats[o.label] = []
+ if OStats.MEAN in stats:
+ out_stats[o.label].append(o.intensity_mean)
+ if OStats.STD in stats:
+ out_stats[o.label].append(o.intensity_std)
+ if OStats.MIN in stats:
+ out_stats[o.label].append(o.intensity_min)
+ if OStats.MAX in stats:
+ out_stats[o.label].append(o.intensity_max)
+ out_stats[o.label] = np.concatenate(out_stats[o.label])
+
+ return out_stats.keys(), np.vstack(list(out_stats.values()))
+
+ def add_raster_stats(self, raster_file, raster_name=None, stats=[OStats.MEAN]):
+ if raster_name is None:
+ raster_name = os.path.splitext(os.path.basename(raster_file))[0]
+
+ ds = gdal.Open(raster_file)
+ n_bands = ds.RasterCount
+ ds = None
+
+ var_names = [raster_name + '_{}_{}'.format(i+1,j.value) for i in range(n_bands) for j in stats]
+ for t in self.tiles.keys():
+ k,v = self.compute_stats_on_tile(t, raster_file, stats)
+ self.database.loc[k,var_names] = v
+
+ return
+
+
+
diff --git a/VHR/segmentation.py b/OBIA/segmentation.py
similarity index 98%
rename from VHR/segmentation.py
rename to OBIA/segmentation.py
index af815013435b52ce359a1f686fa96b718037566a..5be2e7e19f89e60f0e25d633dc0a829b83d14e92 100644
--- a/VHR/segmentation.py
+++ b/OBIA/segmentation.py
@@ -20,7 +20,7 @@ class LSGRMParams:
color_weight: float
spatial_weight: float
n_first_iter: int = 8
- margin : int = 100
+ margin: int = 100
def lsgrm_process_tile(input_image, params : LSGRMParams, tile_width, tile_height, tile_idx, out_graph, roi=None):
if roi is None:
@@ -104,6 +104,7 @@ def lsgrm(input_image, params : LSGRMParams, out_seg, n_proc=None, memory=None,
grm.SetParameterFloat('cw', params.color_weight)
grm.SetParameterFloat('sw', params.spatial_weight)
grm.SetParameterString('out', out_seg)
+ grm.SetParameterOutputImagePixelType('out', otb.ImagePixelType_uint32)
grm.ExecuteAndWriteOutput()
else:
@@ -134,6 +135,7 @@ def lsgrm(input_image, params : LSGRMParams, out_seg, n_proc=None, memory=None,
write_qgis_seg_style(out_seg.replace(ext, '.qml'))
else:
agg_app.SetParameterString('out', out_seg)
+ agg_app.SetParameterOutputImagePixelType('out', otb.ImagePixelType_uint32)
agg_app.ExecuteAndWriteOutput()
if remove_graph:
@@ -265,7 +267,6 @@ def get_bounding_boxes(input_segm):
bboxes[lbl[i]][2] = max(bboxes[lbl[i]][2], rpos[i])
bboxes[lbl[i]][3] = y
-
return bboxes
diff --git a/moringa.py b/moringa.py
index 25fecbe8340c27764476c3168d1fbd999096b9c9..ca764aade1256040f2eb25efa0e833b4456c9598 100755
--- a/moringa.py
+++ b/moringa.py
@@ -1,13 +1,13 @@
import sys
import argparse
-import VHR.segmentation
+import OBIA.segmentation
import VHR.vhrbase
def run_segmentation(img, threshold, cw, sw , out_seg,
n_first_iter, margin, n_proc, memory,
remove_graph, force_parallel):
- params = VHR.segmentation.LSGRMParams(threshold, cw, sw, n_first_iter, margin)
- VHR.segmentation.lsgrm(img,params,out_seg,n_proc,memory,None,remove_graph,force_parallel)
+ params = OBIA.segmentation.LSGRMParams(threshold, cw, sw, n_first_iter, margin)
+ OBIA.segmentation.lsgrm(img, params, out_seg, n_proc, memory, None, remove_graph, force_parallel)
return
def preprocess_spot67(in_fld, out_fld, dem_fld, geoid_file, skip_ps, compress):