1 
2 ///////////////////////////////////////////////////////////
3 //                                                       //
4 //                         SAGA                          //
5 //                                                       //
6 //      System for Automated Geoscientific Analyses      //
7 //                                                       //
8 //                     Tool Library                      //
9 //                  imagery_isocluster                   //
10 //                                                       //
11 //-------------------------------------------------------//
12 //                                                       //
13 //                classify_isocluster.cpp                //
14 //                                                       //
15 //                 Copyright (C) 2016 by                 //
16 //                      Olaf Conrad                      //
17 //                                                       //
18 //-------------------------------------------------------//
19 //                                                       //
20 // This file is part of 'SAGA - System for Automated     //
21 // Geoscientific Analyses'. SAGA is free software; you   //
22 // can redistribute it and/or modify it under the terms  //
23 // of the GNU General Public License as published by the //
24 // Free Software Foundation, either version 3 of the     //
25 // License, or (at your option) any later version.       //
26 //                                                       //
27 // SAGA is distributed in the hope that it will be       //
28 // useful, but WITHOUT ANY WARRANTY; without even the    //
29 // implied warranty of MERCHANTABILITY or FITNESS FOR A  //
30 // PARTICULAR PURPOSE. See the GNU General Public        //
31 // License for more details.                             //
32 //                                                       //
33 // You should have received a copy of the GNU General    //
34 // Public License along with this program; if not,       //
35 // see <http://www.gnu.org/licenses/>.                   //
36 //                                                       //
37 //-------------------------------------------------------//
38 //                                                       //
39 //    e-mail:     oconrad@saga-gis.org                   //
40 //                                                       //
41 //    contact:    Olaf Conrad                            //
42 //                Institute of Geography                 //
43 //                University of Hamburg                  //
44 //                Germany                                //
45 //                                                       //
46 ///////////////////////////////////////////////////////////
47 
48 ///////////////////////////////////////////////////////////
49 //														 //
50 //														 //
51 //														 //
52 ///////////////////////////////////////////////////////////
53 
54 //---------------------------------------------------------
55 #include "classify_isocluster.h"
56 
57 //---------------------------------------------------------
58 #include "cluster_isodata.h"
59 
60 
61 ///////////////////////////////////////////////////////////
62 //														 //
63 //														 //
64 //														 //
65 ///////////////////////////////////////////////////////////
66 
67 //---------------------------------------------------------
CGrid_Cluster_ISODATA(void)68 CGrid_Cluster_ISODATA::CGrid_Cluster_ISODATA(void)
69 {
70 	//-----------------------------------------------------
71 	Set_Name		(_TL("ISODATA Clustering for Grids"));
72 
73 	Set_Author		("O.Conrad (c) 2016");
74 
75 	Set_Description	(_TW(
76 		"This tool executes the Isodata unsupervised "
77 		"classification - clustering algorithm. Isodata "
78 		"stands for Iterative Self-Organizing Data Analysis "
79 		"Techniques. This is a more sophisticated algorithm "
80 		"which allows the number of clusters to be "
81 		"automatically adjusted during the iteration by "
82 		"merging similar clusters and splitting clusters "
83 		"with large standard deviations. "
84 		"The tool is based on Christos Iosifidis' Isodata implementation. "
85 	));
86 
87 	Add_Reference("http://users.ntua.gr/chiossif/Free_As_Freedom_Software/isodata.c",
88 		SG_T("isodata.c (Christos Iosifidis)")
89 	);
90 
91 	Add_Reference("https://www.cs.umd.edu/~mount/Projects/ISODATA",
92 		SG_T("A Fast Implementation of the ISODATA Clustering Algorithm")
93 	);
94 
95 	Add_Reference("Memarsadeghi, N., Mount, D. M., Netanyahu, N. S., Le Moigne, J.", "2007",
96 		"A Fast Implementation of the ISODATA Clustering Algorithm",
97 		"International Journal of Computational Geometry and Applications, 17, 71-103.",
98 		SG_T("https://www.cs.umd.edu/~mount/Projects/ISODATA/ijcga07-isodata.pdf"), SG_T("online")
99 	);
100 
101 	//-----------------------------------------------------
102 	Parameters.Add_Grid_List("",
103 		"FEATURES"		, _TL("Features"),
104 		_TL(""),
105 		PARAMETER_INPUT
106 	);
107 
108 	Parameters.Add_Grid("",
109 		"CLUSTER"		, _TL("Clusters"),
110 		_TL(""),
111 		PARAMETER_OUTPUT, true, SG_DATATYPE_Byte
112 	);
113 
114 	Parameters.Add_Table("",
115 		"STATISTICS"	, _TL("Statistics"),
116 		_TL(""),
117 		PARAMETER_OUTPUT
118 	);
119 
120 	//-----------------------------------------------------
121 	Parameters.Add_Bool("",
122 		"NORMALIZE"		, _TL("Normalize"),
123 		_TL(""),
124 		false
125 	);
126 
127 	Parameters.Add_Int("",
128 		"ITERATIONS"	, _TL("Maximum Number of Iterations"),
129 		_TL(""),
130 		20, 3, true
131 	);
132 
133 	Parameters.Add_Int("",
134 		"CLUSTER_INI"	, _TL("Initial Number of Clusters"),
135 		_TL(""),
136 		5, 0, true
137 	);
138 
139 	Parameters.Add_Int("",
140 		"CLUSTER_MAX"	, _TL("Maximum Number of Clusters"),
141 		_TL(""),
142 		16, 3, true
143 	);
144 
145 	Parameters.Add_Int("",
146 		"SAMPLES_MIN"	, _TL("Minimum Number of Samples in a Cluster"),
147 		_TL(""),
148 		5, 2, true
149 	);
150 
151 	//Parameters.Add_Double("",
152 	//	"DIST_MAX"	, _TL("Distance Threshold"),
153 	//	_TL("Clusters, which are closer than this distance to each other, are merged."),
154 	//	0.001, 0.0, true
155 	//);
156 
157 	//Parameters.Add_Double("",
158 	//	"STDV_MAX"	, _TL("Maximum Standard Deviation within a Cluster"),
159 	//	_TL(""),
160 	//	10.0, 0.0, true
161 	//);
162 
163 	Parameters.Add_Bool("",
164 		"RGB_COLORS"	, _TL("Update Colors from Features"),
165 		_TL("Use the first three features in list to obtain blue, green, red components for class colour in look-up table."),
166 		false
167 	)->Set_UseInCMD(false);
168 
169 	Parameters.Add_Choice("",
170 		"INITIALIZE"	, _TL("Start Partition"),
171 		_TL(""),
172 		CSG_String::Format("%s|%s|%s",
173 			_TL("random"),
174 			_TL("periodical"),
175 			_TL("keep values")
176 		), 0
177 	);
178 }
179 
180 
181 ///////////////////////////////////////////////////////////
182 //														 //
183 ///////////////////////////////////////////////////////////
184 
185 //---------------------------------------------------------
On_Parameter_Changed(CSG_Parameters * pParameters,CSG_Parameter * pParameter)186 int CGrid_Cluster_ISODATA::On_Parameter_Changed(CSG_Parameters *pParameters, CSG_Parameter *pParameter)
187 {
188 	return(CSG_Tool_Grid::On_Parameter_Changed(pParameters, pParameter));
189 }
190 
191 //---------------------------------------------------------
On_Parameters_Enable(CSG_Parameters * pParameters,CSG_Parameter * pParameter)192 int CGrid_Cluster_ISODATA::On_Parameters_Enable(CSG_Parameters *pParameters, CSG_Parameter *pParameter)
193 {
194 	if( pParameter->Cmp_Identifier("FEATURES") )
195 	{
196 		pParameters->Set_Enabled("RGB_COLORS", pParameter->asGridList()->Get_Grid_Count() >= 3);
197 	}
198 
199 	return( CSG_Tool_Grid::On_Parameters_Enable(pParameters, pParameter) );
200 }
201 
202 
203 ///////////////////////////////////////////////////////////
204 //														 //
205 ///////////////////////////////////////////////////////////
206 
207 //---------------------------------------------------------
On_Execute(void)208 bool CGrid_Cluster_ISODATA::On_Execute(void)
209 {
210 	int		iFeature;
211 	sLong	iCell;
212 	size_t	iSample, iCluster;
213 
214 	//-----------------------------------------------------
215 	CSG_Parameter_Grid_List	*pFeatures	= Parameters("FEATURES")->asGridList();
216 
217 	CSG_Grid	*pCluster	= Parameters("CLUSTER")->asGrid();
218 
219 	pCluster->Set_NoData_Value(0.0);
220 
221 	bool	bNormalize	= Parameters("NORMALIZE")->asBool();
222 
223 	//-----------------------------------------------------
224 	TSG_Data_Type	Data_Type;
225 
226 	if( bNormalize )
227 	{
228 		Data_Type	= SG_DATATYPE_Float;
229 	}
230 	else
231 	{
232 		Data_Type	= SG_DATATYPE_Char;
233 
234 		for(iFeature=0; iFeature<pFeatures->Get_Grid_Count(); iFeature++)
235 		{
236 			if( Data_Type < pFeatures->Get_Grid(iFeature)->Get_Type() )
237 			{
238 				Data_Type	= pFeatures->Get_Grid(iFeature)->Get_Type();
239 			}
240 		}
241 
242 		Message_Fmt("\n%s: %s", _TL("internal data type"), SG_Data_Type_Get_Name(Data_Type).c_str());
243 	}
244 
245 	//-----------------------------------------------------
246 	CCluster_ISODATA	Cluster(pFeatures->Get_Grid_Count(), Data_Type);
247 
248 	Cluster.Set_Max_Iterations(Parameters("ITERATIONS" )->asInt   ());
249 	Cluster.Set_Ini_Clusters  (Parameters("CLUSTER_INI")->asInt   ());
250 	Cluster.Set_Max_Clusters  (Parameters("CLUSTER_MAX")->asInt   ());
251 	Cluster.Set_Min_Samples   (Parameters("SAMPLES_MIN")->asInt   ());
252 //	Cluster.Set_Max_Distance  (Parameters("DIST_MAX"   )->asDouble());
253 //	Cluster.Set_Max_StdDev    (Parameters("STDV_MAX"   )->asDouble());
254 
255 	//-----------------------------------------------------
256 	for(iCell=0; iCell<Get_NCells() && Set_Progress_NCells(iCell); iCell++)
257 	{
258 		CSG_Vector	Features(pFeatures->Get_Grid_Count());
259 
260 		for(iFeature=0; Features.Get_Size() && iFeature<pFeatures->Get_Grid_Count(); iFeature++)
261 		{
262 			if( pFeatures->Get_Grid(iFeature)->is_NoData(iCell) )
263 			{
264 				Features.Destroy();
265 			}
266 			else
267 			{
268 				Features[iFeature]	= pFeatures->Get_Grid(iFeature)->asDouble(iCell);
269 
270 				if( bNormalize )
271 				{
272 					Features[iFeature]	= (Features[iFeature] - pFeatures->Get_Grid(iFeature)->Get_Mean()) / pFeatures->Get_Grid(iFeature)->Get_StdDev();
273 				}
274 			}
275 		}
276 
277 		if( Features.Get_Size() )
278 		{
279 			Cluster.Add_Sample(Features);
280 
281 			pCluster->Set_Value(iCell, 1.0);
282 		}
283 		else
284 		{
285 			pCluster->Set_Value(iCell, 0.0);
286 		}
287 	}
288 
289 	//-----------------------------------------------------
290 	if( !Cluster.Run(Parameters("INITIALIZE")->asInt()) )
291 	{
292 		return( false );
293 	}
294 
295 	//-----------------------------------------------------
296 	for(iCell=0, iSample=0; iCell<Get_NCells() && Set_Progress_NCells(iCell); iCell++)
297 	{
298 		if( pCluster->asInt(iCell) )
299 		{
300 			pCluster->Set_Value(iCell, 1 + Cluster.Get_Cluster(iSample++));
301 		}
302 	}
303 
304 	//-----------------------------------------------------
305 	CSG_Table	&Statistics	= *Parameters("STATISTICS")->asTable();
306 
307 	Statistics.Destroy();
308 	Statistics.Set_Name(_TL("ISODATA Cluster Statistics"));
309 
310 	Statistics.Add_Field("CLUSTER" , SG_DATATYPE_Int);
311 	Statistics.Add_Field("ELEMENTS", SG_DATATYPE_Int);
312 	Statistics.Add_Field("MEANDIST", SG_DATATYPE_Double);
313 
314 	for(iFeature=0; iFeature<pFeatures->Get_Grid_Count(); iFeature++)
315 	{
316 		Statistics.Add_Field(CSG_String::Format("MEAN.%s", pFeatures->Get_Grid(iFeature)->Get_Name()), SG_DATATYPE_Double);
317 		Statistics.Add_Field(CSG_String::Format("STDV.%s", pFeatures->Get_Grid(iFeature)->Get_Name()), SG_DATATYPE_Double);
318 	}
319 
320 	for(iCluster=0; iCluster<Cluster.Get_Cluster_Count(); iCluster++)
321 	{
322 		CSG_Table_Record	&Record	= *Statistics.Add_Record();
323 
324 		Record.Set_Value(0, iCluster + 1);
325 		Record.Set_Value(1, Cluster.Get_Cluster_Count (iCluster));
326 		Record.Set_Value(2, Cluster.Get_Cluster_StdDev(iCluster));
327 
328 		for(iFeature=0; iFeature<pFeatures->Get_Grid_Count(); iFeature++)
329 		{
330 			double	Mean	= Cluster.Get_Cluster_Mean  (iCluster, iFeature);
331 			double	Stdv	= Cluster.Get_Cluster_StdDev(iCluster, iFeature);
332 
333 			if( bNormalize )
334 			{
335 				Mean	= Mean * pFeatures->Get_Grid(iFeature)->Get_StdDev() + pFeatures->Get_Grid(iFeature)->Get_Mean();
336 				Stdv	= Stdv * pFeatures->Get_Grid(iFeature)->Get_StdDev();
337 			}
338 
339 			Record.Set_Value(3 + 2 * iFeature + 0, Mean);
340 			Record.Set_Value(3 + 2 * iFeature + 1, Stdv);
341 		}
342 	}
343 
344 	//-----------------------------------------------------
345 	CSG_Parameter	*pLUT	= DataObject_Get_Parameter(pCluster, "LUT");
346 
347 	if( pLUT && pLUT->asTable() )
348 	{
349 		bool	bRGB	= pFeatures->Get_Grid_Count() >= 3 && Parameters("RGB_COLORS")->asBool();
350 
351 		for(iCluster=0; iCluster<Statistics.Get_Count(); iCluster++)
352 		{
353 			CSG_Table_Record	*pClass	= pLUT->asTable()->Get_Record(iCluster);
354 
355 			if( !pClass )
356 			{
357 				(pClass	= pLUT->asTable()->Add_Record())->Set_Value(0, SG_Color_Get_Random());
358 			}
359 
360 			pClass->Set_Value(1, CSG_String::Format("%s %d", _TL("Cluster"), iCluster + 1));
361 			pClass->Set_Value(2, "");
362 			pClass->Set_Value(3, iCluster + 1);
363 			pClass->Set_Value(4, iCluster + 1);
364 
365 			if( bRGB )
366 			{
367 				#define SET_COLOR_COMPONENT(c, i)	c = (int)(127 + (Statistics[iCluster].asDouble(3 + 2 * i) - pFeatures->Get_Grid(i)->Get_Mean()) * 127 / pFeatures->Get_Grid(i)->Get_StdDev()); if( c < 0 ) c = 0; else if( c > 255 ) c = 255;
368 
369 				int	r; SET_COLOR_COMPONENT(r, 2);
370 				int	g; SET_COLOR_COMPONENT(g, 1);
371 				int	b; SET_COLOR_COMPONENT(b, 0);
372 
373 				pClass->Set_Value(0, SG_GET_RGB(r, g, b));
374 			}
375 		}
376 
377 		pLUT->asTable()->Set_Record_Count(Statistics.Get_Count());
378 
379 		DataObject_Set_Parameter(pCluster, pLUT);
380 		DataObject_Set_Parameter(pCluster, "COLORS_TYPE", 1);	// Color Classification Type: Lookup Table
381 	}
382 
383 	//-----------------------------------------------------
384 	return( true );
385 }
386 
387 
388 ///////////////////////////////////////////////////////////
389 //														 //
390 //														 //
391 //														 //
392 ///////////////////////////////////////////////////////////
393 
394 //---------------------------------------------------------
395