1 
2 ///////////////////////////////////////////////////////////
3 //                                                       //
4 //                         SAGA                          //
5 //                                                       //
6 //      System for Automated Geoscientific Analyses      //
7 //                                                       //
8 //                     Tool Library:                     //
9 //                 imagery_segmentation                  //
10 //                                                       //
11 //-------------------------------------------------------//
12 //                                                       //
13 //                       slic.cpp                        //
14 //                                                       //
15 //                  Olaf Conrad (C) 2019                 //
16 //                                                       //
17 //-------------------------------------------------------//
18 //                                                       //
19 // This file is part of 'SAGA - System for Automated     //
20 // Geoscientific Analyses'. SAGA is free software; you   //
21 // can redistribute it and/or modify it under the terms  //
22 // of the GNU General Public License as published by the //
23 // Free Software Foundation, either version 2 of the     //
24 // License, or (at your option) any later version.       //
25 //                                                       //
26 // SAGA is distributed in the hope that it will be       //
27 // useful, but WITHOUT ANY WARRANTY; without even the    //
28 // implied warranty of MERCHANTABILITY or FITNESS FOR A  //
29 // PARTICULAR PURPOSE. See the GNU General Public        //
30 // License for more details.                             //
31 //                                                       //
32 // You should have received a copy of the GNU General    //
33 // Public License along with this program; if not, see   //
34 // <http://www.gnu.org/licenses/>.                       //
35 //                                                       //
36 //-------------------------------------------------------//
37 //                                                       //
38 //    e-mail:     oconrad@saga-gis.org                   //
39 //                                                       //
40 //    contact:    Olaf Conrad                            //
41 //                Institute of Geography                 //
42 //                University of Hamburg                  //
43 //                Germany                                //
44 //                                                       //
45 ///////////////////////////////////////////////////////////
46 
47 //---------------------------------------------------------
48 #include "slic.h"
49 
50 
51 ///////////////////////////////////////////////////////////
52 //														 //
53 //														 //
54 //														 //
55 ///////////////////////////////////////////////////////////
56 
57 //---------------------------------------------------------
CSLIC(void)58 CSLIC::CSLIC(void)
59 {
60 	Set_Name		(_TL("Superpixel Segmentation"));
61 
62 	Set_Author		("O.Conrad (c) 2019");
63 
64 	Set_Description	(_TW(
65 		"The Superpixel Segmentation tool implements the "
66 		"'Simple Linear Iterative Clustering' (SLIC) algorithm, "
67 		"an image segmentation method described in Achanta et al. (2010). "
68 		"\n\n"
69 		"SLIC is a simple and efficient method to decompose an image in "
70 		"visually homogeneous regions. It is based on a spatially "
71 		"localized version of k-means clustering. Similar to mean shift or "
72 		"quick shift, each pixel is associated to a feature vector. "
73 		"\n\n"
74 		"This tool is follows the SLIC implementation created by "
75 		"Vedaldi and Fulkerson as part of the VLFeat library. "
76 	));
77 
78 	Add_Reference("Achanta, R., Shaji, A., Smith, K., Lucchi, A., Fua, P., & S�sstrunk, S.", "2010",
79 		"Slic Superpixels",
80 		"EPFL Technical Report no. 149300, June 2010.",
81 		SG_T("https://infoscience.epfl.ch/record/149300/files/SLIC_Superpixels_TR_2.pdf"), SG_T("epfl.ch")
82 	);
83 
84 	Add_Reference("Achanta, R., Shaji, A., Smith, K., Lucchi, A., Fua, P., & S�sstrunk, S.", "2012",
85 		"SLIC Superpixels compared to state-of-the-art superpixel methods",
86 		"IEEE transactions on pattern analysis and machine intelligence, 34(11), 2274-2282.",
87 		SG_T("https://ieeexplore.ieee.org/iel5/34/4359286/06205760.pdf"), SG_T("ieee.org")
88 	);
89 
90 	Add_Reference(
91 		SG_T("http://www.vlfeat.org/overview/slic.html"), SG_T("SLIC at VLFeat.org")
92 	);
93 
94 	//-----------------------------------------------------
95 	Parameters.Add_Grid_List("",
96 		"FEATURES"		, _TL("Features"),
97 		_TL(""),
98 		PARAMETER_INPUT
99 	);
100 
101 	Parameters.Add_Bool("FEATURES",
102 		"NORMALIZE"		, _TL("Normalize"),
103 		_TL(""),
104 		false
105 	);
106 
107 	Parameters.Add_Shapes("",
108 		"POLYGONS"		, _TL("Segments"),
109 		_TL(""),
110 		PARAMETER_OUTPUT, SHAPE_TYPE_Polygon
111 	);
112 
113 	//-----------------------------------------------------
114 	Parameters.Add_Int("",
115 		"MAX_ITERATIONS", _TL("Maximum Iterations"),
116 		_TL(""),
117 		100, 1, true
118 	);
119 
120 	Parameters.Add_Double("",
121 		"REGULARIZATION", _TL("Regularization"),
122 		_TL(""),
123 		1., 0., true
124 	);
125 
126 	Parameters.Add_Int("",
127 		"SIZE"			, _TL("Region Size"),
128 		_TL("Starting 'cell size' of the superpixels given as number of cells."),
129 		10, 1, true
130 	);
131 
132 	Parameters.Add_Int("",
133 		"SIZE_MIN"		, _TL("Minimum Region Size"),
134 		_TL("In postprocessing join segments, which cover less cells than specified here, to a larger neighbour segment."),
135 		1, 1, true
136 	);
137 
138 	//-----------------------------------------------------
139 	Parameters.Add_Bool("",
140 		"SUPERPIXELS_DO", _TL("Create Superpixel Grids"),
141 		_TL(""),
142 		false
143 	);
144 
145 	Parameters.Add_Grid_List("",
146 		"SUPERPIXELS"	, _TL("Superpixels"),
147 		_TL(""),
148 		PARAMETER_OUTPUT_OPTIONAL
149 	);
150 
151 	//-----------------------------------------------------
152 	Parameters.Add_Choice("",
153 		"POSTPROCESSING", _TL("Post-Processing"),
154 		_TL(""),
155 		CSG_String::Format("%s|%s",
156 			_TL("none"),
157 			_TL("unsupervised classification")
158 		), 0
159 	);
160 
161 	Parameters.Add_Int("POSTPROCESSING",
162 		"NCLUSTER"		, _TL("Number of Clusters"),
163 		_TL(""),
164 		12, 2, true
165 	);
166 
167 	Parameters.Add_Bool("POSTPROCESSING",
168 		"SPLIT_CLUSTERS", _TL("Split Clusters"),
169 		_TL(""),
170 		true
171 	);
172 
173 	//-----------------------------------------------------
174 	m_Centroid	= NULL;
175 }
176 
177 
178 ///////////////////////////////////////////////////////////
179 //														 //
180 ///////////////////////////////////////////////////////////
181 
182 //---------------------------------------------------------
On_Parameter_Changed(CSG_Parameters * pParameters,CSG_Parameter * pParameter)183 int CSLIC::On_Parameter_Changed(CSG_Parameters *pParameters, CSG_Parameter *pParameter)
184 {
185 	return( CSG_Tool_Grid::On_Parameter_Changed(pParameters, pParameter) );
186 }
187 
188 //---------------------------------------------------------
On_Parameters_Enable(CSG_Parameters * pParameters,CSG_Parameter * pParameter)189 int CSLIC::On_Parameters_Enable(CSG_Parameters *pParameters, CSG_Parameter *pParameter)
190 {
191 	if( pParameter->Cmp_Identifier("POSTPROCESSING") )
192 	{
193 		pParameter->Set_Children_Enabled(pParameter->asInt() == 1);
194 	}
195 
196 	return( CSG_Tool_Grid::On_Parameters_Enable(pParameters, pParameter) );
197 }
198 
199 
200 ///////////////////////////////////////////////////////////
201 //														 //
202 ///////////////////////////////////////////////////////////
203 
204 //---------------------------------------------------------
On_Execute(void)205 bool CSLIC::On_Execute(void)
206 {
207 	m_pGrids	= Parameters("FEATURES")->asGridList();
208 
209 	m_bNormalize	= Parameters("NORMALIZE")->asBool();
210 
211 	//-----------------------------------------------------
212 	CSG_Grid	Segments;
213 
214 	if( !Get_Segments(Segments) )
215 	{
216 		Del_Centroids();
217 
218 		return( false );
219 	}
220 
221 	Get_Generalized(Segments);
222 
223 	//-----------------------------------------------------
224 	Get_Grids(Segments);
225 
226 	bool	bResult	= Get_Polygons(Segments);
227 
228 	Parameters("POLYGONS")->asShapes()->Set_Name(_TL("Segments"));
229 
230 	Del_Centroids();
231 
232 	return( bResult );
233 }
234 
235 
236 ///////////////////////////////////////////////////////////
237 //														 //
238 ///////////////////////////////////////////////////////////
239 
240 //---------------------------------------------------------
Get_Feature_Count(void)241 inline int CSLIC::Get_Feature_Count(void)
242 {
243 	return( m_pGrids->Get_Grid_Count() );
244 }
245 
246 //---------------------------------------------------------
Get_Feature(int k,int x,int y)247 inline double CSLIC::Get_Feature(int k, int x, int y)
248 {
249 	CSG_Grid	*pGrid	= m_pGrids->Get_Grid(k);
250 
251 	double	Value	= pGrid->asDouble(x, y);
252 
253 	if( m_bNormalize && pGrid->Get_StdDev() > 0. )
254 	{
255 		Value	= (Value - pGrid->Get_Min()) / pGrid->Get_StdDev();
256 	}
257 
258 	return( Value );
259 }
260 
261 //---------------------------------------------------------
Fit_To_Grid_System(double Value,int Coordinate)262 inline int CSLIC::Fit_To_Grid_System(double Value, int Coordinate)
263 {
264 	int	i	= (int)floor(Value + 0.5);
265 
266 	switch( Coordinate )
267 	{
268 	default: return( M_GET_MAX(M_GET_MIN(i, Get_NX() - 1), 0) );
269 	case  1: return( M_GET_MAX(M_GET_MIN(i, Get_NY() - 1), 0) );
270 	}
271 }
272 
273 
274 ///////////////////////////////////////////////////////////
275 //														 //
276 ///////////////////////////////////////////////////////////
277 
278 //---------------------------------------------------------
Get_Polygons(CSG_Grid & Segments)279 bool CSLIC::Get_Polygons(CSG_Grid &Segments)
280 {
281 	CSG_Shapes	*pPolygons	= Parameters("POLYGONS")->asShapes();
282 
283 	CSG_Tool	*pTool	= SG_Get_Tool_Library_Manager().Create_Tool("shapes_grid", 6);	// Vectorising Grid Classes
284 
285 	if( (pPolygons && pTool && pTool->Set_Manager(NULL)
286 		&&  pTool->Set_Parameter("CLASS_ALL"  , 1        )
287 		&&  pTool->Set_Parameter("SPLIT"      , 0        )
288 		&&  pTool->Set_Parameter("ALLVERTICES", false    )
289 		&&  pTool->Set_Parameter("GRID"       , &Segments)
290 		&&  pTool->Set_Parameter("POLYGONS"   , pPolygons)
291 		&&  pTool->Execute()) == false )
292 	{
293 		SG_Get_Tool_Library_Manager().Delete_Tool(pTool);
294 
295 		return( false );
296 	}
297 
298 	SG_Get_Tool_Library_Manager().Delete_Tool(pTool);
299 
300 	//-----------------------------------------------------
301 	pPolygons->Del_Field(pPolygons->Get_Field("NAME"));
302 	pPolygons->Del_Field(pPolygons->Get_Field("ID"  ));
303 
304 	for(int k=0; k<Get_Feature_Count(); k++)
305 	{
306 		pPolygons->Add_Field(m_pGrids->Get_Grid(k)->Get_Name(), SG_DATATYPE_Double);
307 	}
308 
309 	for(int i=0; i<pPolygons->Get_Count(); i++)
310 	{
311 		CSG_Shape	*pPolygon	= pPolygons->Get_Shape(i);
312 
313 		sLong	ID	= pPolygon->asInt(0);
314 
315 		for(int k=0; k<Get_Feature_Count(); k++)
316 		{
317 			if( ID >= 0 && ID < m_Centroid->Get_NCells() )
318 			{
319 				pPolygon->Set_Value(1 + k, m_Centroid[2 + k].asDouble(ID));
320 			}
321 			else
322 			{
323 				pPolygon->Set_NoData(1 + k);
324 			}
325 		}
326 	}
327 
328 	//-----------------------------------------------------
329 	if( Parameters("POSTPROCESSING")->asInt() == 0 )
330 	{
331 		return( true );
332 	}
333 
334 	CSG_Table	Statistics;
335 
336 	if( ((pTool = SG_Get_Tool_Library_Manager().Create_Tool("table_calculus", 14)) && pTool->Set_Manager(NULL) // Cluster Analysis (Shapes)
337 		&&  pTool->Set_Parameter("NCLUSTER"  , Parameters("NCLUSTER" ))
338 		&&  pTool->Set_Parameter("NORMALISE" , Parameters("NORMALIZE"))
339 		&&  pTool->Set_Parameter("STATISTICS", &Statistics)
340 		&&  pTool->Set_Parameter("INPUT"     , pPolygons)
341 		&&  pTool->Set_Parameter("FIELDS"    , "1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32")
342 		&&  pTool->Execute()) == false )
343 	{
344 		SG_Get_Tool_Library_Manager().Delete_Tool(pTool);
345 
346 		return( false );
347 	}
348 
349 	SG_Get_Tool_Library_Manager().Delete_Tool(pTool);
350 
351 	//-----------------------------------------------------
352 	CSG_Shapes	Polygons(SHAPE_TYPE_Polygon);
353 
354 	if( ((pTool = SG_Get_Tool_Library_Manager().Create_Tool("shapes_polygons", 5)) && pTool->Set_Manager(NULL) // Polygon Dissolve
355 		&&  pTool->Set_Parameter("POLYGONS"  , pPolygons)
356 		&&  pTool->Set_Parameter("DISSOLVED" , &Polygons)
357 		&&  pTool->Set_Parameter("FIELDS"    , "CLUSTER")
358 		&&  pTool->Execute()) == false )
359 	{
360 		SG_Get_Tool_Library_Manager().Delete_Tool(pTool);
361 
362 		return( false );
363 	}
364 
365 	SG_Get_Tool_Library_Manager().Delete_Tool(pTool);
366 
367 	//-----------------------------------------------------
368 	if( Parameters("SPLIT_CLUSTERS")->asBool() == false )
369 	{
370 		return( pPolygons->Create(Polygons) );
371 	}
372 
373 	if( ((pTool = SG_Get_Tool_Library_Manager().Create_Tool("shapes_polygons", 10)) && pTool->Set_Manager(NULL) // Polygon Parts to Separate Polygons
374 		&&  pTool->Set_Parameter("POLYGONS"  , &Polygons)
375 		&&  pTool->Set_Parameter("PARTS"     , pPolygons)
376 		&&  pTool->Execute()) == false )
377 	{
378 		SG_Get_Tool_Library_Manager().Delete_Tool(pTool);
379 
380 		return( false );
381 	}
382 
383 	SG_Get_Tool_Library_Manager().Delete_Tool(pTool);
384 
385 	return( true );
386 }
387 
388 //---------------------------------------------------------
Get_Grids(CSG_Grid & Segments)389 bool CSLIC::Get_Grids(CSG_Grid &Segments)
390 {
391 	if( Parameters("SUPERPIXELS_DO")->asBool() == false )
392 	{
393 		return( true ); // nothing to do !
394 	}
395 
396 	//-----------------------------------------------------
397 	CSG_Parameter_Grid_List	*pGrids	= Parameters("SUPERPIXELS")->asGridList();
398 
399 	pGrids->Del_Items();
400 
401 	for(int i=0; i<m_pGrids->Get_Item_Count(); i++)
402 	{
403 		CSG_Data_Object	*pItem = m_pGrids->Get_Item(i), *pNew = NULL;
404 
405 		switch( pItem->Get_ObjectType() )
406 		{
407 		case SG_DATAOBJECT_TYPE_Grid : pNew = SG_Create_Grid (*pItem->asGrid ()); break;
408 		case SG_DATAOBJECT_TYPE_Grids: pNew = SG_Create_Grids(*pItem->asGrids()); break;
409 		default                      : pNew = NULL;
410 		}
411 
412 		if( !pNew )
413 		{
414 			return( false );
415 		}
416 
417 		pNew->Fmt_Name("%s [%s]", pItem->Get_Name(), _TL("SLIC"));
418 
419 		pGrids->Add_Item(pNew);
420 	}
421 
422 	//-----------------------------------------------------
423 	for(int k=0; k<pGrids->Get_Grid_Count(); k++)
424 	{
425 		CSG_Grid	*pGrid = pGrids->Get_Grid(k), *pCentroid = m_Centroid + 2 + k;
426 
427 		for(sLong Cell=0; Cell<Get_NCells(); Cell++)
428 		{
429 			sLong	ID	= Segments.asLong(Cell);
430 
431 			for(int k=0; k<Get_Feature_Count(); k++)
432 			{
433 				if( ID >= 0 && ID < m_Centroid->Get_NCells() )
434 				{
435 					pGrid->Set_Value(Cell, pCentroid->asDouble(ID));
436 				}
437 				else
438 				{
439 					pGrid->Set_NoData(Cell);
440 				}
441 			}
442 		}
443 	}
444 
445 	//-----------------------------------------------------
446 	for(int i=0; i<m_pGrids->Get_Item_Count(); i++)
447 	{
448 		DataObject_Add           (pGrids->Get_Item(i));
449 		DataObject_Set_Parameters(pGrids->Get_Item(i), m_pGrids->Get_Item(i));
450 	}
451 
452 	return( true );
453 }
454 
455 
456 ///////////////////////////////////////////////////////////
457 //														 //
458 ///////////////////////////////////////////////////////////
459 
460 //---------------------------------------------------------
Get_Segments(CSG_Grid & Segments)461 bool CSLIC::Get_Segments(CSG_Grid &Segments)
462 {
463 	int	Size	= Parameters("SIZE")->asInt();
464 
465 	if( Size < 1 || !Get_Centroids(Size) )
466 	{
467 		Error_Set(_TL("failed to initialize centroids"));
468 
469 		return( false );
470 	}
471 
472 	//-----------------------------------------------------
473 	Process_Set_Text(_TL("running k-means iterations"));
474 
475 	Segments.Create(Get_System(), SG_DATATYPE_Word);
476 
477 	CSG_Grid	Masses(Get_System(), SG_DATATYPE_Word);
478 
479 	double	Energy_0, Energy_Last = -1.;
480 
481 	double	factor	= Parameters("REGULARIZATION")->asDouble() / (Size*Size);
482 
483 	int	max_Iterations	= Parameters("MAX_ITERATIONS")->asInt();
484 
485 	//-----------------------------------------------------
486 	for(int Iteration=0; Iteration<max_Iterations && Process_Get_Okay(); Iteration++)
487 	{
488 		double	Energy	= 0.;
489 
490 		for(int y=0; y<Get_NY(); y++)	// assign pixels to centers
491 		{
492 			for(int x=0; x<Get_NX(); x++)
493 			{
494 				int	cx	= (int)floor((double)x / Size - 0.5);
495 				int	cy	= (int)floor((double)y / Size - 0.5);
496 
497 				double	min_Distance	= -1.;
498 
499 				for(int iy=M_GET_MAX(0, cy); iy<=M_GET_MIN(m_Centroid->Get_NY()-1, cy+1); iy++)
500 				{
501 					for(int ix=M_GET_MAX(0, cx); ix<=M_GET_MIN(m_Centroid->Get_NX()-1, cx+1); ix++)
502 					{
503 						int	min_x	= m_Centroid[0].asInt(ix, iy);
504 						int	min_y	= m_Centroid[1].asInt(ix, iy);	//	double min_y	= centers[(2 + Get_Feature_Count()) * region + 1];
505 
506 						double	appearance	= 0., spatial = (x - min_x)*(x - min_x) + (y - min_y)*(y - min_y);
507 
508 						for(int k=0; k<Get_Feature_Count(); k++)
509 						{
510 							double	cz	= m_Centroid[2 + k].asDouble(ix, iy);
511 							double	z	= Get_Feature(k, x, y);
512 							appearance	+= (z - cz)*(z - cz);
513 						}
514 
515 						double	Distance	= appearance + factor * spatial;
516 
517 						if( min_Distance > Distance || min_Distance < 0. )
518 						{
519 							min_Distance	= Distance;
520 
521 							Segments.Set_Value(x, y, ix + iy * m_Centroid->Get_NX());
522 						}
523 					}
524 				}
525 
526 				Energy	+= min_Distance;
527 			}
528 		}
529 
530 		//-------------------------------------------------
531 		// check energy termination conditions
532 
533 		Process_Set_Text(CSG_String::Format("%s %d, %s: %f", _TL("iteration"), 1 + Iteration, _TL("energy"), Energy));
534 
535 		if( Iteration < 1 )
536 		{
537 			Energy_0	= Energy;
538 		}
539 		else if( (Energy_Last - Energy) < 1e-5 * (Energy_0 - Energy) )
540 		{
541 			break;
542 		}
543 
544 		Energy_Last	= Energy;
545 
546 		//-------------------------------------------------
547 		// recompute centers
548 
549 		Masses.Assign(0.);
550 
551 		for(int i=0; i<Get_Feature_Count()+2; i++)
552 		{
553 			m_Centroid[i].Assign(0.);
554 		}
555 
556 		for(int y=0; y<Get_NY(); y++)
557 		{
558 			for(int x=0; x<Get_NX(); x++)
559 			{
560 				sLong	i	= Segments.asInt(x, y);
561 
562 				Masses.Add_Value(i, 1.);
563 
564 				m_Centroid[0].Add_Value(i, x);
565 				m_Centroid[1].Add_Value(i, y);
566 
567 				for(int k=0; k<Get_Feature_Count(); k++)
568 				{
569 					m_Centroid[2 + k].Add_Value(i, Get_Feature(k, x, y));
570 				}
571 			}
572 		}
573 
574 		#pragma omp parallel for
575 		for(sLong i=0; i<m_Centroid->Get_NCells(); i++)
576 		{
577 			double	Mass	= 1. / M_GET_MAX(Masses.asDouble(i), 1e-8);
578 
579 			for(int k=0; k<Get_Feature_Count()+2; k++)
580 			{
581 				m_Centroid[k].Mul_Value(i, Mass);
582 			}
583 		}
584 	}
585 
586 	return( true );
587 }
588 
589 
590 ///////////////////////////////////////////////////////////
591 //														 //
592 ///////////////////////////////////////////////////////////
593 
594 //---------------------------------------------------------
Get_Edge(CSG_Grid & Edge)595 bool CSLIC::Get_Edge(CSG_Grid &Edge)
596 {
597 	if( !Edge.Create(Get_System(), SG_DATATYPE_Float) )
598 	{
599 		Error_Set(_TL("failed to create edge map"));
600 
601 		return( false );
602 	}
603 
604 	//-----------------------------------------------------
605 	Process_Set_Text(_TL("computing edge map"));	// compute edge map (gradient strength)
606 
607 	#pragma omp parallel for
608 	for(int y=1; y<Get_NY()-1; y++)
609 	{
610 		for(int x=1; x<Get_NX()-1; x++)
611 		{
612 			for(int k=0; k<Get_Feature_Count(); k++)
613 			{
614 				double	a	= Get_Feature(k, x - 1, y    );
615 				double	b	= Get_Feature(k, x + 1, y    );
616 				double	c	= Get_Feature(k, x    , y + 1);
617 				double	d	= Get_Feature(k, x    , y - 1);
618 
619 				Edge.Add_Value(x, y, (a - b)*(a - b) + (c - d)*(c - d));
620 			}
621 		}
622 	}
623 
624 	return( true );
625 }
626 
627 //---------------------------------------------------------
Get_Centroids(int Size)628 bool CSLIC::Get_Centroids(int Size)
629 {
630 	CSG_Grid	Edge;
631 
632 	if( !Get_Edge(Edge) )
633 	{
634 		return( false );
635 	}
636 
637 	//-----------------------------------------------------
638 	Process_Set_Text(_TL("initializing k-means centroids"));	// initialize k-means centroids
639 
640 	m_Centroid	= new CSG_Grid[2 + Get_Feature_Count()];
641 
642 	if( !m_Centroid )
643 	{
644 		return( false );
645 	}
646 
647 	CSG_Grid_System	System(Get_Cellsize() / Size, Get_XMin(), Get_YMin(),
648 		(int)ceil((double)Get_NX() / Size),
649 		(int)ceil((double)Get_NY() / Size)
650 	);
651 
652 	if( !m_Centroid[0].Create(System, SG_DATATYPE_Word)
653 	||  !m_Centroid[1].Create(System, SG_DATATYPE_Word) )
654 	{
655 		return( false );
656 	}
657 
658 	for(int k=0; k<Get_Feature_Count(); k++)
659 	{
660 		if( !m_Centroid[2 + k].Create(System, SG_DATATYPE_Float) )
661 		{
662 			return( false );
663 		}
664 	}
665 
666 	//-----------------------------------------------------
667 	#pragma omp parallel for
668 	for(int cy=0; cy<System.Get_NY(); cy++)
669 	{
670 		for(int cx=0; cx<System.Get_NX(); cx++)
671 		{
672 			int	x	= Fit_To_Grid_System(Size * (cx + 0.5), 0);
673 			int	y	= Fit_To_Grid_System(Size * (cy + 0.5), 1);
674 
675 			double	min_e	= -1.;
676 			int		min_x	= 0;
677 			int		min_y	= 0;
678 
679 			for(int iy=M_GET_MAX(0, y-1); iy<=M_GET_MIN(Get_NY()-1, y+1); iy++)	// search in a 3x3 neighbourhood the smallest edge response
680 			{
681 				for(int ix=M_GET_MAX(0, x-1); ix<=M_GET_MIN(Get_NX()-1, x+1); ix++)
682 				{
683 					double	ie	= Edge.asDouble(ix, iy);
684 
685 					if( min_e > ie || min_e < 0. )
686 					{
687 						min_e	= ie;
688 						min_x	= ix;
689 						min_y	= iy;
690 					}
691 				}
692 			}
693 
694 			m_Centroid[0].Set_Value(cx, cy, min_x);	// initialize the new center at this location
695 			m_Centroid[1].Set_Value(cx, cy, min_y);
696 
697 			for(int k=0; k<Get_Feature_Count(); k++)
698 			{
699 				m_Centroid[2 + k].Set_Value(cx, cy, Get_Feature(k, min_x, min_y));
700 			}
701 		}
702 	}
703 
704 	//-----------------------------------------------------
705 	return( true );
706 }
707 
708 //---------------------------------------------------------
Del_Centroids(void)709 bool CSLIC::Del_Centroids(void)
710 {
711 	if( m_Centroid )
712 	{
713 		for(int i=0; i<Get_Feature_Count()+2; i++)
714 		{
715 			m_Centroid[i].Destroy();
716 		}
717 
718 		delete[](m_Centroid);
719 
720 		m_Centroid	= NULL;
721 
722 		return( true );
723 	}
724 
725 	return( false );
726 }
727 
728 
729 ///////////////////////////////////////////////////////////
730 //														 //
731 ///////////////////////////////////////////////////////////
732 
733 //---------------------------------------------------------
Get_Generalized(CSG_Grid & Segments)734 bool CSLIC::Get_Generalized(CSG_Grid &Segments)
735 {
736 	int	min_Size	= Parameters("SIZE_MIN")->asInt();
737 
738 	if( min_Size <= 1 )
739 	{
740 		return( true );	// nothing to do !
741 	}
742 
743 	//-----------------------------------------------------
744 	Process_Set_Text(_TL("eliminating small regions"));
745 
746 	CSG_Grid	Cleaned(Get_System(), SG_DATATYPE_Int);
747 	CSG_Grid	Segment(Get_System(), SG_DATATYPE_Int);
748 
749 	//-----------------------------------------------------
750 	for(sLong Cell=0; Cell<Get_NCells(); Cell++)
751 	{
752 		if( Cleaned.asInt(Cell) )
753 		{
754 			continue;
755 		}
756 
757 		//-------------------------------------------------
758 		int	ID	= Segments.asInt(Cell), Size = 0;
759 
760 		int	ID_Cleaned	= ID + 1;
761 
762 		Segment.Set_Value(Size++, (int)Cell);
763 
764 		Cleaned.Set_Value(Cell, ID + 1);
765 
766 		for(int i=0; i<8; i+=2)	// find ID_Cleaned as the ID of an already cleaned region neihbour of this pixel
767 		{
768 			int	ix	= Get_xTo(i, (int)(Cell % Get_NX()));
769 			int	iy	= Get_yTo(i, (int)(Cell / Get_NX()));
770 
771 			if( is_InGrid(ix, iy) && Cleaned.asInt(ix, iy) != 0 )
772 			{
773 				ID_Cleaned	= Cleaned.asInt(ix, iy);
774 			}
775 		}
776 
777 		//-------------------------------------------------
778 		for(int	Expanded=0; Expanded<Size; Expanded++)	// expand the segment
779 		{
780 			sLong	open	= Segment.asLong(Expanded);
781 
782 			for(int i=0; i<8; i+=2)
783 			{
784 				int	ix	= Get_xTo(i, (int)(open % Get_NX()));
785 				int	iy	= Get_yTo(i, (int)(open / Get_NX()));
786 
787 				if( is_InGrid(ix, iy) && Cleaned.asInt(ix, iy) == 0 && Segments.asInt(ix, iy) == ID )
788 				{
789 					Cleaned.Set_Value(ix, iy, ID + 1);
790 
791 					Segment.Set_Value(Size++, ix + iy * Get_NX());
792 				}
793 			}
794 		}
795 
796 		//-------------------------------------------------
797 		if( Size < min_Size )	// change label to ID_Cleaned if the segment is too small
798 		{
799 			while( Size > 0 )
800 			{
801 				Cleaned.Set_Value(Segment.asInt(--Size), ID_Cleaned);
802 			}
803 		}
804 	}
805 
806 	//-----------------------------------------------------
807 	for(sLong i=0; i<Get_NCells(); i++)	// restore base 0 indexing of the regions
808 	{
809 		Segments.Set_Value(i, Cleaned.asDouble(i) - 1);
810 	}
811 
812 	return( true );
813 }
814 
815 
816 ///////////////////////////////////////////////////////////
817 //														 //
818 //														 //
819 //														 //
820 ///////////////////////////////////////////////////////////
821 
822 //---------------------------------------------------------
823