1 /**********************************************************
2 * Version $Id: classify_cluster_analysis.cpp 1921 2014-01-09 10:24:11Z oconrad $
3 *********************************************************/
4
5 ///////////////////////////////////////////////////////////
6 // //
7 // SAGA //
8 // //
9 // System for Automated Geoscientific Analyses //
10 // //
11 // Tool Library //
12 // imagery_classification //
13 // //
14 //-------------------------------------------------------//
15 // //
16 // Grid_Cluster_Analysis.cpp //
17 // //
18 // Copyright (C) 2003 by //
19 // Olaf Conrad //
20 // //
21 //-------------------------------------------------------//
22 // //
23 // This file is part of 'SAGA - System for Automated //
24 // Geoscientific Analyses'. SAGA is free software; you //
25 // can redistribute it and/or modify it under the terms //
26 // of the GNU General Public License as published by the //
27 // Free Software Foundation, either version 2 of the //
28 // License, or (at your option) any later version. //
29 // //
30 // SAGA is distributed in the hope that it will be //
31 // useful, but WITHOUT ANY WARRANTY; without even the //
32 // implied warranty of MERCHANTABILITY or FITNESS FOR A //
33 // PARTICULAR PURPOSE. See the GNU General Public //
34 // License for more details. //
35 // //
36 // You should have received a copy of the GNU General //
37 // Public License along with this program; if not, see //
38 // <http://www.gnu.org/licenses/>. //
39 // //
40 //-------------------------------------------------------//
41 // //
42 // e-mail: oconrad@saga-gis.org //
43 // //
44 // contact: Olaf Conrad //
45 // Institute of Geography //
46 // University of Goettingen //
47 // Goldschmidtstr. 5 //
48 // 37077 Goettingen //
49 // Germany //
50 // //
51 ///////////////////////////////////////////////////////////
52
53 //---------------------------------------------------------
54
55
56 ///////////////////////////////////////////////////////////
57 // //
58 // //
59 // //
60 ///////////////////////////////////////////////////////////
61
62 //---------------------------------------------------------
63 #include "classify_cluster_analysis.h"
64
65
66 ///////////////////////////////////////////////////////////
67 // //
68 // //
69 // //
70 ///////////////////////////////////////////////////////////
71
72 //---------------------------------------------------------
CGrid_Cluster_Analysis(void)73 CGrid_Cluster_Analysis::CGrid_Cluster_Analysis(void)
74 {
75 //-----------------------------------------------------
76 Set_Name (_TL("K-Means Clustering for Grids"));
77
78 Set_Author ("O.Conrad (c) 2001");
79
80 Set_Description (_TW(
81 "This tool implements the K-Means cluster analysis for grids "
82 "in two variants, iterative minimum distance (Forgy 1965) "
83 "and hill climbing (Rubin 1967). "
84 ));
85
86 Add_Reference("Forgy, E.", "1965",
87 "Cluster analysis of multivariate data: efficiency vs. interpretability of classifications",
88 "Biometrics 21:768."
89 );
90
91 Add_Reference("Rubin, J.", "1967",
92 "Optimal classification into groups: an approach for solving the taxonomy problem",
93 "J. Theoretical Biology, 15:103-144."
94 );
95
96 //-----------------------------------------------------
97 Parameters.Add_Grid_List("",
98 "GRIDS" , _TL("Grids"),
99 _TL(""),
100 PARAMETER_INPUT
101 );
102
103 Parameters.Add_Grid("",
104 "CLUSTER" , _TL("Clusters"),
105 _TL(""),
106 PARAMETER_OUTPUT, true, SG_DATATYPE_Byte
107 );
108
109 Parameters.Add_Table("",
110 "STATISTICS" , _TL("Statistics"),
111 _TL(""),
112 PARAMETER_OUTPUT
113 );
114
115 Parameters.Add_Choice("",
116 "METHOD" , _TL("Method"),
117 _TL(""),
118 CSG_String::Format("%s|%s|%s",
119 _TL("Iterative Minimum Distance (Forgy 1965)"),
120 _TL("Hill-Climbing (Rubin 1967)"),
121 _TL("Combined Minimum Distance / Hillclimbing")
122 ), 1
123 );
124
125 Parameters.Add_Int("",
126 "NCLUSTER" , _TL("Clusters"),
127 _TL("Number of clusters"),
128 10, 2, true
129 );
130
131 Parameters.Add_Int("",
132 "MAXITER" , _TL("Maximum Iterations"),
133 _TL("maximum number of iterations, ignored if set to zero (default)"),
134 10, 0, true
135 );
136
137 Parameters.Add_Bool("",
138 "NORMALISE" , _TL("Normalise"),
139 _TL("Automatically normalise grids by standard deviation before clustering."),
140 false
141 );
142
143 Parameters.Add_Bool("",
144 "RGB_COLORS" , _TL("Update Colors from Features"),
145 _TL("Use the first three features in list to obtain blue, green, red components for class colour in look-up table."),
146 false
147 )->Set_UseInCMD(false);
148
149 Parameters.Add_Choice("",
150 "INITIALIZE" , _TL("Start Partition"),
151 _TL(""),
152 CSG_String::Format("%s|%s|%s",
153 _TL("random"),
154 _TL("periodical"),
155 _TL("keep values")
156 ), 0
157 );
158
159 //-----------------------------------------------------
160 Parameters.Add_Bool("" , "OLDVERSION", _TL("Old Version"), _TL("slower but memory saving"), false);
161 Parameters.Add_Bool("OLDVERSION", "UPDATEVIEW", _TL("Update View"), _TL(""), true);
162 }
163
164
165 ///////////////////////////////////////////////////////////
166 // //
167 ///////////////////////////////////////////////////////////
168
169 //---------------------------------------------------------
On_Parameters_Enable(CSG_Parameters * pParameters,CSG_Parameter * pParameter)170 int CGrid_Cluster_Analysis::On_Parameters_Enable(CSG_Parameters *pParameters, CSG_Parameter *pParameter)
171 {
172 if( pParameter->Cmp_Identifier("OLDVERSION") )
173 {
174 pParameters->Set_Enabled("INITIALIZE", pParameter->asBool() == false);
175 pParameters->Set_Enabled("UPDATEVIEW", pParameter->asBool() == true );
176 }
177
178 if( pParameter->Cmp_Identifier("GRIDS") )
179 {
180 pParameters->Set_Enabled("RGB_COLORS", pParameter->asGridList()->Get_Grid_Count() >= 3);
181 }
182
183 return( CSG_Tool_Grid::On_Parameters_Enable(pParameters, pParameter) );
184 }
185
186
187 ///////////////////////////////////////////////////////////
188 // //
189 ///////////////////////////////////////////////////////////
190
191 //---------------------------------------------------------
On_Execute(void)192 bool CGrid_Cluster_Analysis::On_Execute(void)
193 {
194 if( Parameters("OLDVERSION")->asBool() ) { return( _On_Execute() ); }
195
196 //-----------------------------------------------------
197 CSG_Parameter_Grid_List *pGrids = Parameters("GRIDS")->asGridList();
198
199 CSG_Grid *pCluster = Parameters("CLUSTER")->asGrid();
200
201 bool bNormalize = Parameters("NORMALISE")->asBool();
202
203 //-----------------------------------------------------
204 CSG_Cluster_Analysis Analysis;
205
206 if( !Analysis.Create(pGrids->Get_Grid_Count()) )
207 {
208 return( false );
209 }
210
211 //-----------------------------------------------------
212 int iFeature;
213 sLong iElement, nElements;
214
215 pCluster->Set_NoData_Value(0);
216
217 //-----------------------------------------------------
218 for(iElement=0, nElements=0; iElement<Get_NCells() && Set_Progress_NCells(iElement); iElement++)
219 {
220 bool bNoData = false;
221
222 for(iFeature=0; iFeature<pGrids->Get_Grid_Count() && !bNoData; iFeature++)
223 {
224 if( pGrids->Get_Grid(iFeature)->is_NoData(iElement) )
225 {
226 bNoData = true;
227 }
228 }
229
230 if( bNoData || !Analysis.Add_Element() )
231 {
232 pCluster->Set_Value(iElement, 0);
233 }
234 else
235 {
236 pCluster->Set_Value(iElement, 1);
237
238 for(iFeature=0; iFeature<pGrids->Get_Grid_Count(); iFeature++)
239 {
240 double d = pGrids->Get_Grid(iFeature)->asDouble(iElement);
241
242 if( bNormalize )
243 {
244 d = (d - pGrids->Get_Grid(iFeature)->Get_Mean()) / pGrids->Get_Grid(iFeature)->Get_StdDev();
245 }
246
247 Analysis.Set_Feature(nElements, iFeature, d);
248 }
249
250 nElements++;
251 }
252 }
253
254 if( Analysis.Get_nElements() < 2 )
255 {
256 return( false );
257 }
258
259 //-----------------------------------------------------
260 bool bResult = Analysis.Execute(
261 Parameters("METHOD" )->asInt(),
262 Parameters("NCLUSTER" )->asInt(),
263 Parameters("MAXITER" )->asInt(),
264 Parameters("INITIALIZE")->asInt()
265 );
266
267 for(iElement=0, nElements=0; iElement<Get_NCells(); iElement++)
268 {
269 Set_Progress_NCells(iElement);
270
271 if( !pCluster->is_NoData(iElement) )
272 {
273 pCluster->Set_Value(iElement, 1 + Analysis.Get_Cluster(nElements++));
274 }
275 }
276
277 Save_Statistics(pGrids, bNormalize, Analysis);
278
279 Save_LUT(pCluster);
280
281 return( bResult );
282 }
283
284
285 ///////////////////////////////////////////////////////////
286 // //
287 ///////////////////////////////////////////////////////////
288
289 //---------------------------------------------------------
Save_Statistics(CSG_Parameter_Grid_List * pGrids,bool bNormalize,const CSG_Cluster_Analysis & Analysis)290 void CGrid_Cluster_Analysis::Save_Statistics(CSG_Parameter_Grid_List *pGrids, bool bNormalize, const CSG_Cluster_Analysis &Analysis)
291 {
292 int iCluster, iFeature;
293 CSG_String s;
294 CSG_Table *pTable = Parameters("STATISTICS")->asTable();
295
296 pTable->Destroy();
297 pTable->Set_Name(_TL("Cluster Analysis"));
298
299 pTable->Add_Field(_TL("ClusterID"), SG_DATATYPE_Int );
300 pTable->Add_Field(_TL("Elements" ), SG_DATATYPE_Int );
301 pTable->Add_Field(_TL("Std.Dev." ), SG_DATATYPE_Double);
302
303 s.Printf("\n%s:\t%d \n%s:\t%ld \n%s:\t%d \n%s:\t%d \n%s:\t%f\n\n%s\t%s\t%s",
304 _TL("Number of Iterations"), Analysis.Get_Iteration(),
305 _TL("Number of Elements" ), Analysis.Get_nElements(),
306 _TL("Number of Variables" ), Analysis.Get_nFeatures(),
307 _TL("Number of Clusters" ), Analysis.Get_nClusters(),
308 _TL("Standard Deviation" ), sqrt(Analysis.Get_SP()),
309 _TL("Cluster"), _TL("Elements"), _TL("Std.Dev.")
310 );
311
312 for(iFeature=0; iFeature<Analysis.Get_nFeatures(); iFeature++)
313 {
314 s += CSG_String::Format("\t%s", pGrids->Get_Grid(iFeature)->Get_Name());
315
316 pTable->Add_Field(pGrids->Get_Grid(iFeature)->Get_Name(), SG_DATATYPE_Double);
317 }
318
319 Message_Add(s);
320
321 for(iCluster=0; iCluster<Analysis.Get_nClusters(); iCluster++)
322 {
323 s.Printf("\n%d\t%d\t%f", iCluster, Analysis.Get_nMembers(iCluster), sqrt(Analysis.Get_Variance(iCluster)));
324
325 CSG_Table_Record *pRecord = pTable->Add_Record();
326
327 pRecord->Set_Value(0, iCluster);
328 pRecord->Set_Value(1, Analysis.Get_nMembers(iCluster));
329 pRecord->Set_Value(2, sqrt(Analysis.Get_Variance(iCluster)));
330
331 for(iFeature=0; iFeature<Analysis.Get_nFeatures(); iFeature++)
332 {
333 double Centroid = Analysis.Get_Centroid(iCluster, iFeature);
334
335 if( bNormalize )
336 {
337 Centroid = pGrids->Get_Grid(iFeature)->Get_Mean() + Centroid * pGrids->Get_Grid(iFeature)->Get_StdDev();
338 }
339
340 s += CSG_String::Format("\t%f", Centroid);
341
342 pRecord->Set_Value(iFeature + 3, Centroid);
343 }
344
345 Message_Add(s, false);
346 }
347 }
348
349 //---------------------------------------------------------
Save_LUT(CSG_Grid * pCluster)350 void CGrid_Cluster_Analysis::Save_LUT(CSG_Grid *pCluster)
351 {
352 CSG_Parameter *pLUT = DataObject_Get_Parameter(pCluster, "LUT");
353
354 if( pLUT && pLUT->asTable() )
355 {
356 CSG_Parameter_Grid_List *pGrids = Parameters("GRIDS")->asGridList();
357
358 CSG_Table &Statistics = *Parameters("STATISTICS")->asTable();
359
360 bool bRGB = pGrids->Get_Grid_Count() >= 3 && Parameters("RGB_COLORS")->asBool();
361
362 for(int iCluster=0; iCluster<Statistics.Get_Count(); iCluster++)
363 {
364 CSG_Table_Record *pClass = pLUT->asTable()->Get_Record(iCluster);
365
366 if( !pClass )
367 {
368 (pClass = pLUT->asTable()->Add_Record())->Set_Value(0, SG_Color_Get_Random());
369 }
370
371 pClass->Set_Value(1, CSG_String::Format("%s %d", _TL("Cluster"), iCluster + 1));
372 pClass->Set_Value(2, "");
373 pClass->Set_Value(3, iCluster + 1);
374 pClass->Set_Value(4, iCluster + 1);
375
376 if( bRGB )
377 {
378 #define SET_COLOR_COMPONENT(c, i) c = (int)(127 + (Statistics[iCluster].asDouble(3 + i) - pGrids->Get_Grid(i)->Get_Mean()) * 127 / pGrids->Get_Grid(i)->Get_StdDev()); if( c < 0 ) c = 0; else if( c > 255 ) c = 255;
379
380 int r; SET_COLOR_COMPONENT(r, 2);
381 int g; SET_COLOR_COMPONENT(g, 1);
382 int b; SET_COLOR_COMPONENT(b, 0);
383
384 pClass->Set_Value(0, SG_GET_RGB(r, g, b));
385 }
386 }
387
388 pLUT->asTable()->Set_Record_Count(Statistics.Get_Count());
389
390 DataObject_Set_Parameter(pCluster, pLUT);
391 DataObject_Set_Parameter(pCluster, "COLORS_TYPE", 1); // Color Classification Type: Lookup Table
392
393 DataObject_Update(pCluster);
394 }
395 }
396
397
398 ///////////////////////////////////////////////////////////
399 // //
400 // //
401 // //
402 // //
403 // //
404 // //
405 ///////////////////////////////////////////////////////////
406 // //
407 // Deprecated Old Version //
408 // //
409 ///////////////////////////////////////////////////////////
410 // //
411 // slow, but saves memory ! //
412 // //
413 // //
414 // //
415 // //
416 ///////////////////////////////////////////////////////////
417
418 //---------------------------------------------------------
_On_Execute(void)419 bool CGrid_Cluster_Analysis::_On_Execute(void)
420 {
421 int i, j, *nMembers, nCluster, nElements;
422 double *Variances, **Centroids, SP;
423 CSG_Grid **Grids, *pCluster;
424 CSG_Parameter_Grid_List *pGrids;
425
426 //-----------------------------------------------------
427 pGrids = Parameters("GRIDS") ->asGridList();
428 pCluster = Parameters("CLUSTER") ->asGrid();
429 nCluster = Parameters("NCLUSTER")->asInt();
430
431 if( pGrids->Get_Grid_Count() < 1 )
432 {
433 return( false );
434 }
435
436 //-----------------------------------------------------
437 Grids = (CSG_Grid **)SG_Malloc(pGrids->Get_Grid_Count() * sizeof(CSG_Grid *));
438
439 if( Parameters("NORMALISE")->asBool() )
440 {
441 for(i=0; i<pGrids->Get_Grid_Count(); i++)
442 {
443 Grids[i] = SG_Create_Grid(pGrids->Get_Grid(i), SG_DATATYPE_Float);
444 Grids[i] ->Assign(pGrids->Get_Grid(i));
445 Grids[i] ->Standardise();
446 }
447 }
448 else
449 {
450 for(i=0; i<pGrids->Get_Grid_Count(); i++)
451 {
452 Grids[i] = pGrids->Get_Grid(i);
453 }
454 }
455
456 pCluster->Set_NoData_Value(-1.0);
457 pCluster->Assign_NoData();
458
459 nMembers = (int *)SG_Malloc(nCluster * sizeof(int));
460 Variances = (double *)SG_Malloc(nCluster * sizeof(double));
461 Centroids = (double **)SG_Malloc(nCluster * sizeof(double *));
462
463 for(i=0; i<nCluster; i++)
464 {
465 Centroids[i] = (double *)SG_Malloc(pGrids->Get_Grid_Count() * sizeof(double));
466 }
467
468 //-------------------------------------------------
469 switch( Parameters("METHOD")->asInt() )
470 {
471 case 0: SP = _MinimumDistance (Grids, pGrids->Get_Grid_Count(), pCluster, nCluster, nMembers, Variances, Centroids, nElements = Get_NCells()); break;
472 case 1: SP = _HillClimbing (Grids, pGrids->Get_Grid_Count(), pCluster, nCluster, nMembers, Variances, Centroids, nElements = Get_NCells()); break;
473 case 2: SP = _MinimumDistance (Grids, pGrids->Get_Grid_Count(), pCluster, nCluster, nMembers, Variances, Centroids, nElements = Get_NCells());
474 SP = _HillClimbing (Grids, pGrids->Get_Grid_Count(), pCluster, nCluster, nMembers, Variances, Centroids, nElements = Get_NCells()); break;
475 }
476
477 //-------------------------------------------------
478 if( Parameters("NORMALISE")->asBool() )
479 {
480 for(i=0; i<pGrids->Get_Grid_Count(); i++)
481 {
482 delete(Grids[i]);
483
484 for(j=0; j<nCluster; j++)
485 {
486 Centroids[j][i] = pGrids->Get_Grid(i)->Get_StdDev() * Centroids[j][i] + pGrids->Get_Grid(i)->Get_Mean();
487 }
488 }
489 }
490
491 //-------------------------------------------------
492 int iCluster, iFeature;
493 CSG_String s;
494 CSG_Table_Record *pRecord;
495 CSG_Table *pTable;
496
497 pTable = Parameters("STATISTICS")->asTable();
498
499 pTable->Destroy();
500 pTable->Set_Name(_TL("Cluster Analysis"));
501
502 pTable->Add_Field(_TL("ClusterID"), SG_DATATYPE_Int );
503 pTable->Add_Field(_TL("Elements" ), SG_DATATYPE_Int );
504 pTable->Add_Field(_TL("Std.Dev." ), SG_DATATYPE_Double);
505
506 s.Printf("\n%s:\t%ld \n%s:\t%d \n%s:\t%d \n%s:\t%f\n\n%s\t%s\t%s",
507 _TL("Number of Elements" ), nElements,
508 _TL("Number of Variables"), pGrids->Get_Grid_Count(),
509 _TL("Number of Clusters" ), nCluster,
510 _TL("Standard Deviation" ), sqrt(SP),
511 _TL("Cluster"), _TL("Elements"), _TL("Std.Dev.")
512 );
513
514 for(iFeature=0; iFeature<pGrids->Get_Grid_Count(); iFeature++)
515 {
516 s += CSG_String::Format("\t%s", pGrids->Get_Grid(iFeature)->Get_Name());
517
518 pTable->Add_Field(pGrids->Get_Grid(iFeature)->Get_Name(), SG_DATATYPE_Double);
519 }
520
521 Message_Add(s);
522
523 for(iCluster=0; iCluster<nCluster; iCluster++)
524 {
525 Variances[iCluster] = nMembers[iCluster] ? Variances[iCluster] / nMembers[iCluster] : 0.0;
526
527 s.Printf("\n%d\t%d\t%f", iCluster, nMembers[iCluster], sqrt(Variances[iCluster]));
528
529 pRecord = pTable->Add_Record();
530 pRecord->Set_Value(0, iCluster);
531 pRecord->Set_Value(1, nMembers[iCluster]);
532 pRecord->Set_Value(2, sqrt(Variances[iCluster]));
533
534 for(iFeature=0; iFeature<pGrids->Get_Grid_Count(); iFeature++)
535 {
536 double Centroid = Centroids[iCluster][iFeature];
537
538 if( Parameters("NORMALISE")->asBool() )
539 {
540 Centroid = pGrids->Get_Grid(iFeature)->Get_Mean() + Centroid * pGrids->Get_Grid(iFeature)->Get_StdDev();
541 }
542
543 s += CSG_String::Format("\t%f", Centroid);
544
545 pRecord->Set_Value(iFeature + 3, Centroid);
546 }
547
548 Message_Add(s, false);
549 }
550
551 //-------------------------------------------------
552 Save_LUT(pCluster);
553
554 //-------------------------------------------------
555 for(i=0; i<nCluster; i++)
556 {
557 SG_Free(Centroids[i]);
558 }
559
560 SG_Free(Centroids);
561 SG_Free(Variances);
562 SG_Free(nMembers);
563 SG_Free(Grids);
564
565 return( true );
566 }
567
568 //---------------------------------------------------------
_MinimumDistance(CSG_Grid ** Grids,int nGrids,CSG_Grid * pCluster,int nCluster,int * nMembers,double * Variances,double ** Centroids,int & nElements)569 double CGrid_Cluster_Analysis::_MinimumDistance(CSG_Grid **Grids, int nGrids, CSG_Grid *pCluster, int nCluster, int *nMembers, double *Variances, double **Centroids, int &nElements)
570 {
571 bool bContinue;
572 int iElement, iGrid, iCluster, nClusterElements, nShifts, minCluster, nPasses;
573 double d, Variance, minVariance, SP, SP_Last = -1;
574
575 //-----------------------------------------------------
576 for(iElement=0, nClusterElements=0; iElement<nElements; iElement++)
577 {
578 for(iGrid=0, bContinue=true; iGrid<nGrids && bContinue; iGrid++)
579 {
580 if( Grids[iGrid]->is_NoData(iElement) )
581 {
582 bContinue = false;
583 }
584 }
585
586 if( bContinue )
587 {
588 if( pCluster->asInt(iElement) < 0 || pCluster->asInt(iElement) >= nCluster )
589 {
590 pCluster->Set_Value(iElement, iElement % nCluster);
591 }
592
593 nClusterElements++;
594 }
595 else
596 {
597 pCluster->Set_Value(iElement, -1);
598 }
599 }
600
601 if( Parameters("UPDATEVIEW")->asBool() )
602 {
603 DataObject_Update(pCluster, 0, nCluster, 1);
604 }
605
606 //-----------------------------------------------------
607 int maxIter = Parameters("MAXITER")->asInt();
608
609 for(nPasses=1, bContinue=true; bContinue && (maxIter==0 || nPasses<=maxIter) && Process_Get_Okay(false); nPasses++)
610 {
611 for(iCluster=0; iCluster<nCluster; iCluster++)
612 {
613 Variances[iCluster] = 0;
614 nMembers [iCluster] = 0;
615
616 for(iGrid=0; iGrid<nGrids; iGrid++)
617 {
618 Centroids[iCluster][iGrid] = 0;
619 }
620 }
621
622 //-------------------------------------------------
623 for(iElement=0; iElement<nElements; iElement++)
624 {
625 if( pCluster->asInt(iElement) >= 0 )
626 {
627 iCluster = pCluster->asInt(iElement);
628 nMembers[iCluster]++;
629
630 for(iGrid=0; iGrid<nGrids; iGrid++)
631 {
632 Centroids[iCluster][iGrid] += Grids[iGrid]->asDouble(iElement);
633 }
634 }
635 }
636
637 //-------------------------------------------------
638 for(iCluster=0; iCluster<nCluster; iCluster++)
639 {
640 d = nMembers[iCluster] > 0 ? 1.0 / (double)nMembers[iCluster] : 0;
641
642 for(iGrid=0; iGrid<nGrids; iGrid++)
643 {
644 Centroids[iCluster][iGrid] *= d;
645 }
646 }
647
648 //-------------------------------------------------
649 SP = 0;
650 nShifts = 0;
651
652 for(iElement=0; iElement<nElements && bContinue; iElement++)
653 {
654 if( pCluster->asInt(iElement) >= 0 )
655 {
656 minVariance = -1;
657
658 for(iCluster=0; iCluster<nCluster; iCluster++)
659 {
660 Variance = 0;
661
662 for(iGrid=0; iGrid<nGrids; iGrid++)
663 {
664 d = Centroids[iCluster][iGrid] - Grids[iGrid]->asDouble(iElement);
665 Variance += d * d;
666 }
667
668 if( minVariance<0 || Variance<minVariance )
669 {
670 minVariance = Variance;
671 minCluster = iCluster;
672 }
673 }
674
675 if( pCluster->asInt(iElement) != minCluster )
676 {
677 pCluster->Set_Value(iElement, minCluster);
678 nShifts++;
679 }
680
681 SP += minVariance;
682 Variances[minCluster] += minVariance;
683 }
684 }
685
686 //-------------------------------------------------
687 if( nShifts == 0 )//|| (SP_Last >= 0 && SP >= SP_Last) )
688 {
689 bContinue = false;
690 }
691
692 SP /= nElements;
693
694 Process_Set_Text("%s: %d >> %s %f", _TL("pass"), nPasses, _TL("change"), SP_Last < 0.0 ? SP : SP_Last - SP);
695
696 SP_Last = SP;
697
698 if( Parameters("UPDATEVIEW")->asBool() )
699 {
700 DataObject_Update(pCluster);
701 }
702 }
703
704 nElements = nClusterElements;
705
706 return( SP );
707 }
708
709 //---------------------------------------------------------
_HillClimbing(CSG_Grid ** Grids,int nGrids,CSG_Grid * pCluster,int nCluster,int * nMembers,double * Variances,double ** Centroids,int & nElements)710 double CGrid_Cluster_Analysis::_HillClimbing(CSG_Grid **Grids, int nGrids, CSG_Grid *pCluster, int nCluster, int *nMembers, double *Variances, double **Centroids, int &nElements)
711 {
712 bool bContinue;
713 int iElement, iGrid, iCluster, jCluster, kCluster, nClusterElements, noShift, nPasses;
714 double d, e, n_iK, n_jK, V, VMin, V1, V2, SP, SP_Last = -1;
715
716 //-----------------------------------------------------
717 for(iCluster=0; iCluster<nCluster; iCluster++)
718 {
719 Variances[iCluster] = 0;
720 nMembers [iCluster] = 0;
721
722 for(iGrid=0; iGrid<nGrids; iGrid++)
723 {
724 Centroids[iCluster][iGrid] = 0;
725 }
726 }
727
728 //-----------------------------------------------------
729 for(iElement=0, nClusterElements=0; iElement<nElements; iElement++)
730 {
731 for(iGrid=0, bContinue=true; iGrid<nGrids && bContinue; iGrid++)
732 {
733 if( Grids[iGrid]->is_NoData(iElement) )
734 {
735 bContinue = false;
736 }
737 }
738
739 if( bContinue )
740 {
741 if( pCluster->asInt(iElement) < 0 || pCluster->asInt(iElement) >= nCluster )
742 {
743 pCluster->Set_Value(iElement, iElement % nCluster);
744 }
745
746 nClusterElements++;
747
748 iCluster = pCluster->asInt(iElement);
749
750 nMembers[iCluster]++;
751
752 V = 0.0;
753
754 for(iGrid=0; iGrid<nGrids; iGrid++)
755 {
756 d = Grids[iGrid]->asDouble(iElement);
757 Centroids[iCluster][iGrid] += d;
758 V += d * d;
759 }
760
761 Variances[iCluster] += V;
762 }
763 else
764 {
765 pCluster->Set_Value(iElement, -1);
766 }
767 }
768
769 //-----------------------------------------------------
770 for(iCluster=0; iCluster<nCluster; iCluster++)
771 {
772 d = nMembers[iCluster] != 0 ? 1.0 / (double)nMembers[iCluster] : 0;
773 V = 0.0;
774
775 for(iGrid=0; iGrid<nGrids; iGrid++)
776 {
777 Centroids[iCluster][iGrid] *= d;
778 e = Centroids[iCluster][iGrid];
779 V += e * e;
780 }
781
782 Variances[iCluster] -= nMembers [iCluster] * V;
783 }
784
785 if( Parameters("UPDATEVIEW")->asBool() )
786 {
787 DataObject_Update(pCluster, 0, nCluster, 1);
788 }
789
790 //-----------------------------------------------------
791 noShift = 0;
792
793 int maxIter = Parameters("MAXITER")->asInt();
794
795 for(nPasses=1, bContinue=true; bContinue && (maxIter==0 || nPasses<=maxIter) && Process_Get_Okay(false); nPasses++)
796 {
797 for(iElement=0; iElement<nElements && bContinue; iElement++)
798 {
799 if( pCluster->asInt(iElement) >= 0 )
800 {
801 if( noShift++ >= nElements )
802 {
803 bContinue = false;
804 }
805 else
806 {
807
808 //-------------------------------------
809 iCluster = pCluster->asInt(iElement);
810
811 if( nMembers[iCluster] > 1 )
812 {
813 V = 0.0;
814
815 for(iGrid=0; iGrid<nGrids; iGrid++)
816 {
817 d = Centroids[iCluster][iGrid] - Grids[iGrid]->asDouble(iElement);
818 V += d * d;
819 }
820
821 n_iK = nMembers[iCluster];
822 V1 = V * n_iK / (n_iK - 1.0);
823 VMin = -1.0;
824
825 //---------------------------------
826 for(jCluster=0; jCluster<nCluster; jCluster++)
827 {
828 if( jCluster != iCluster )
829 {
830 V = 0.0;
831
832 for(iGrid=0; iGrid<nGrids; iGrid++)
833 {
834 d = Centroids[jCluster][iGrid] - Grids[iGrid]->asDouble(iElement);
835 V += d * d;
836 }
837
838 n_jK = nMembers[jCluster];
839 V2 = V * n_jK / (n_jK + 1.0);
840
841 if( VMin < 0 || V2 < VMin )
842 {
843 VMin = V2;
844 kCluster = jCluster;
845 }
846 }
847 }
848
849 //---------------------------------
850 if( VMin >= 0 && VMin < V1 )
851 {
852 noShift = 0;
853 Variances[iCluster] -= V1;
854 Variances[kCluster] += VMin;
855 V1 = 1.0 / (n_iK - 1.0);
856 n_jK = nMembers[kCluster];
857 V2 = 1.0 / (n_jK + 1.0);
858
859 for(iGrid=0; iGrid<nGrids; iGrid++)
860 {
861 d = Grids[iGrid]->asDouble(iElement);
862 Centroids[iCluster][iGrid] = (n_iK * Centroids[iCluster][iGrid] - d) * V1;
863 Centroids[kCluster][iGrid] = (n_jK * Centroids[kCluster][iGrid] + d) * V2;
864 }
865
866 pCluster->Set_Value(iElement, kCluster);
867
868 nMembers[iCluster]--;
869 nMembers[kCluster]++;
870 }
871 }
872 }
873 }
874 }
875
876 //-------------------------------------------------
877 for(iCluster=0, SP=0.0; iCluster<nCluster; iCluster++)
878 {
879 SP += Variances[iCluster];
880 }
881
882 SP /= nElements;
883
884 Process_Set_Text("%s: %d >> %s %f", _TL("pass"), nPasses, _TL("change"), SP_Last < 0.0 ? SP : SP_Last - SP);
885
886 SP_Last = SP;
887
888 if( Parameters("UPDATEVIEW")->asBool() )
889 {
890 DataObject_Update(pCluster);
891 }
892 }
893
894 nElements = nClusterElements;
895
896 return( SP );
897 }
898
899
900 ///////////////////////////////////////////////////////////
901 // //
902 // //
903 // //
904 ///////////////////////////////////////////////////////////
905
906 //---------------------------------------------------------
907