diff --git a/Common/otb_numpy_proc.py b/Common/otb_numpy_proc.py
index 09122ac60440455608cd317a5b8eef921e8d87d2..b784ce42a0a0b891495e33b2ccc824af43e18114 100644
--- a/Common/otb_numpy_proc.py
+++ b/Common/otb_numpy_proc.py
@@ -194,7 +194,7 @@ def check_function_padding(function, test_tile_size=(100,100)):
def to_otb_pipeline(obj):
if isinstance(obj, str):
- ppl = otb.Registry.CreateApplication('ExtractROI')
+ ppl = otb.Registry.CreateApplicationWithoutLogger('ExtractROI')
ppl.SetParameterString('in',obj)
ppl.Execute()
elif isinstance(obj,otb.Application):
diff --git a/Learning/ObjectBased.py b/Learning/ObjectBased.py
index 17c70e3b6e401d9702dcc035532f2c7629b9982f..46c459e8a4dd697f0e7ad827eeff77f79b957382 100644
--- a/Learning/ObjectBased.py
+++ b/Learning/ObjectBased.py
@@ -94,13 +94,41 @@ class ObjectBasedClassifier:
}
return models, summary, results
+ def classify(self, model):
+ if isinstance(model, list):
+ for t, L, X in self.obia_base.tiled_data(normalize=[self.training_base['perc2'],
+ self.training_base['perc98']]):
+ prob = []
+ for m in model:
+ prob.append(m.predict_proba(X))
+ prob = np.prod(prob, axis=0)
+ c = model[0].classes_[np.argmax(prob, axis=1)]
+ self.obia_base.populate_map(t, L, c)
+ else:
+ for t,L,X in self.obia_base.tiled_data(normalize=[self.training_base['perc2'],
+ self.training_base['perc98']]):
+ c = model.predict(X)
+ self.obia_base.populate_map(t, L, c)
+ return
+
#TEST CODE
def run_test():
+ '''
obc = ObjectBasedClassifier('/DATA/Moringa_Sample/Parakou/output/segmentation/segmentation.tif',
'/DATA/Moringa_Sample/Parakou/input/REF/ref.shp',
['/DATA/Moringa_Sample/Parakou/output/S2_processed/T31PDL/*/*FEAT.tif'],
['/DATA/Moringa_Sample/Parakou/input/THR/THR_SPOT6.tif'])
+ '''
+ print('Initializing the OBIA training base...')
+ obc = ObjectBasedClassifier('/DATA/Benin/OBSYDYA_data/MORINGA/SEGMENTATION/segmentation.tif',
+ '/DATA/Benin/OBSYDYA_data/MORINGA/reference/BD_OBSYDYA_2022_ParakouNdali_v0.2.shp',
+ ['/DATA/Benin/OBSYDYA_data/MORINGA/basefolder/FEAT/S2_THEIA_FEAT/S2_THEIA_MOSAIC_*.tif'],
+ glob.glob('/DATA/Benin/OBSYDYA_data/MORINGA/ext_features'))
obc.gen_k_folds(5)
- m,s,r = obc.train_RF(100)
+ print('Performing Training and Cross-Validation...')
+ m,s,r = obc.train_RF(400)
print(s)
- return
+ print('Performing Classification...')
+ obc.classify(m)
+ obc.obia_base.export_raster_map('/DATA/Benin/OBSYDYA_data/MORINGAv2/firstmap.tif')
+ return obc
diff --git a/OBIA/OBIABase.py b/OBIA/OBIABase.py
index e900d4aacf14cba9fbd7627ccf8330fcb1f09b21..03e5c6dacdefa2fc49db6f8ae583ea685930f3e1 100644
--- a/OBIA/OBIABase.py
+++ b/OBIA/OBIABase.py
@@ -48,7 +48,11 @@ class OBIABase:
self.raster_groups = []
self.n_vars = 0
+ # Output raster
+ self.output_map = None
+
def init_ref_db(self, vector_file, id_field, class_field):
+ print("Enter init")
ras_id = otb.Registry.CreateApplication('Rasterization')
ras_id.SetParameterString('in', vector_file)
ras_id.SetParameterString('im', self.object_layer)
@@ -56,6 +60,7 @@ class OBIABase:
ras_id.SetParameterString('mode.attribute.field', id_field)
ras_id.Execute()
ids = ras_id.GetImageAsNumpyArray('out')
+ print("First raster")
ras_cl = otb.Registry.CreateApplication('Rasterization')
ras_cl.SetParameterString('in', vector_file)
@@ -64,15 +69,18 @@ class OBIABase:
ras_cl.SetParameterString('mode.attribute.field', class_field)
ras_cl.Execute()
cls = ras_cl.GetImageAsNumpyArray('out')
+ print("Second raster")
in_seg = to_otb_pipeline(self.object_layer)
obj = in_seg.GetImageAsNumpyArray('out')
+ print("Load segmentation")
mask = ids>0
# This was to relabel object being disconnected by intersecting the GT
# Re-enable if needed...
# self.ref_obj_layer = label(obj * mask, connectivity=1).astype(np.uint32)
self.ref_obj_layer = (obj * mask).astype(np.uint32)
+ print("Before Regionprops")
rp = regionprops(self.ref_obj_layer, intensity_image=np.dstack((obj,ids,cls)).astype(np.int))
self.ref_db = pd.DataFrame(data=[np.insert(o.intensity_min,0,o.area) for o in rp],
columns=['area','orig_label','polygon_id','class'],
@@ -80,11 +88,35 @@ class OBIABase:
return
def generate_tile_info(self, nominal_tile_size, background_label):
+ print('Enter generate_tile')
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')
+ r = otb.itkRegion()
+ obj_bbox = {}
+ for i in range(0, W, nominal_tile_size[0]):
+ for j in range(0, H, nominal_tile_size[1]):
+ r['index'][0], r['index'][1] = i, j
+ r['size'][0], r['size'][1] = min(nominal_tile_size[0], W-i), min(nominal_tile_size[1], H-j)
+ in_seg.PropagateRequestedRegion('out', r)
+ tile = in_seg.GetImageAsNumpyArray('out')
+ rp = regionprops(tile.astype(np.uint32))
+ for o in rp:
+ bbox = np.array([o.bbox[0]+j, o.bbox[1]+i, o.bbox[2]+j, o.bbox[3]+i])
+ if o.label in obj_bbox.keys():
+ obj_bbox[o.label].append(bbox)
+ else:
+ obj_bbox[o.label] = [bbox]
+
+ for o in obj_bbox.keys():
+ obj_bbox[o] = np.array(obj_bbox[o])
+ mn, mx = np.min(obj_bbox[o], axis=0), np.max(obj_bbox[o], axis=0)
+ obj_bbox[o] = [mn[0], mn[1], mx[2], mx[3]]
+
+ print("Finished BBOX.")
+ print("BBOX:")
+ print(obj_bbox)
+ input("Pippo...")
+
tx, ty = int(W / nominal_tile_size[0]) + 1, int(H / nominal_tile_size[1]) + 1
obj_to_tile = {}
@@ -104,6 +136,8 @@ class OBIABase:
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])}
+ self.W, self.H = W, H
+
return tiles, obj_to_tile
def compute_stats_on_tile(self, tile_num, input_image, stats=[OStats.MEAN], on_ref=False):
@@ -246,15 +280,41 @@ class OBIABase:
vars = [item for sublist in self.raster_var_names for item in sublist]
for tilenum in self.tiles.keys():
tile_data = self.get_full_stats_on_tile(tilenum)
- L = tile_data.index.to_numpy(dtype=np.int)
+ L = tile_data.index.to_numpy(dtype=np.int32)
X = tile_data[vars].to_numpy()
if normalize is not None:
X = (X - normalize[0]) / (normalize[1] - normalize[0])
- yield L,X
+ yield tilenum,L,X
- def populate_map(self, tilenum, labels):
+ def populate_map(self, tilenum, obj_id, classes):
+ if self.output_map is None:
+ self.output_map = np.zeros((self.H,self.W), dtype=np.int)
+ r = otb.itkRegion()
+ r['index'][0], r['index'][1], r['size'][0], r['size'][1] = self.tiles[tilenum]
+ 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)
+ tmp = np.zeros(np.max(tile_obj) + 1, dtype=int)
+ tmp[obj_id] = classes
+ self.output_map[self.tiles[tilenum][1]:self.tiles[tilenum][1]+self.tiles[tilenum][3],
+ self.tiles[tilenum][0]:self.tiles[tilenum][0]+self.tiles[tilenum][2]] += tmp[tile_obj]
return
+ def export_raster_map(self, output_file, compress=False):
+ assert (self.output_map is not None)
+ if not os.path.exists(os.path.dirname(output_file)):
+ os.makedirs(os.path.dirname(output_file))
+ obj = to_otb_pipeline(self.object_layer)
+ arr = obj.ExportImage('out')
+ arr['array'] = self.output_map.astype(np.float32)
+ obj.ImportImage('in', arr)
+ if compress:
+ output_file = output_file + "?gdal:co:compress=deflate"
+ obj.SetParameterString('out', output_file)
+ obj.SetParameterOutputImagePixelType('out', otb.ImagePixelType_uint16)
+ obj.ExecuteAndWriteOutput()
+
def true_pred_bypixel(self, labels, predicted_classes):
pred_c = np.zeros(np.max(self.ref_obj_layer)+1)
pred_c[labels] = predicted_classes
diff --git a/moringa.py b/moringa.py
index de7a3d60256e3490877f9192b61bc0b4104a0bfb..840250ff03f0cdfb52a91cb5a1aba1f7b774ab0f 100755
--- a/moringa.py
+++ b/moringa.py
@@ -91,7 +91,7 @@ def main(args):
vhrprep.add_argument("--align_to_band", type=int, default=3, help="Band of reference image used for alignment.")
vhrprep.add_argument("--align_using_band", type=int, default=1, help="Band of current image used for alignment.")
vhrprep.add_argument("--skip_ps", help="Skip pansharpening step.", action='store_true')
- vhrprep.add_argument("--compress", help="Use lossless compresion on outputs.", action='store_true')
+ vhrprep.add_argument("--compress", help="Use lossless compression on outputs.", action='store_true')