SAGA API Version 9.13
Loading...
Searching...
No Matches
mat_tools.cpp
Go to the documentation of this file.
1
3// //
4// SAGA //
5// //
6// System for Automated Geoscientific Analyses //
7// //
8// Application Programming Interface //
9// //
10// Library: SAGA_API //
11// //
12//-------------------------------------------------------//
13// //
14// mat_tools.cpp //
15// //
16// Copyright (C) 2005 by Olaf Conrad //
17// //
18//-------------------------------------------------------//
19// //
20// This file is part of 'SAGA - System for Automated //
21// Geoscientific Analyses'. //
22// //
23// This library is free software; you can redistribute //
24// it and/or modify it under the terms of the GNU Lesser //
25// General Public License as published by the Free //
26// Software Foundation, either version 2.1 of the //
27// License, or (at your option) any later version. //
28// //
29// This library is distributed in the hope that it will //
30// be useful, but WITHOUT ANY WARRANTY; without even the //
31// implied warranty of MERCHANTABILITY or FITNESS FOR A //
32// PARTICULAR PURPOSE. See the GNU Lesser General Public //
33// License for more details. //
34// //
35// You should have received a copy of the GNU Lesser //
36// General Public License along with this program; if //
37// not, see <http://www.gnu.org/licenses/>. //
38// //
39//-------------------------------------------------------//
40// //
41// contact: Olaf Conrad //
42// Institute of Geography //
43// University of Goettingen //
44// Goldschmidtstr. 5 //
45// 37077 Goettingen //
46// Germany //
47// //
48// e-mail: oconrad@saga-gis.org //
49// //
51
52//---------------------------------------------------------
53#include <time.h>
54#include <cfloat>
55
56#include "mat_tools.h"
57
58#include "table.h"
59#include "grid.h"
60#include "grids.h"
61
62
64// //
65// //
66// //
68
69//---------------------------------------------------------
70double SG_Get_Square(double Value)
71{
72 return( Value * Value );
73}
74
75
77// //
78// //
79// //
81
82//---------------------------------------------------------
83double SG_Get_Rounded(double Value, int Decimals)
84{
85 if( Decimals < 0 )
86 {
87 return( Value );
88 }
89
90 if( Decimals == 0 )
91 {
92 return( floor(0.5 + Value) );
93 }
94
95 double d = pow(10., Decimals);
96 double v = Value * d;
97
98 if( fabs(v - floor(v)) > 0. )
99 {
100 return( floor(0.5 + v) / d );
101 }
102
103 return( Value );
104}
105
106//---------------------------------------------------------
107double SG_Get_Rounded_To_SignificantFigures(double Value, int Decimals)
108{
109 if( Decimals <= 0 || Value == 0. )
110 {
111 return( (int)(0.5 + Value) );
112 }
113
114 Decimals = (int)(-(ceil(log10(fabs(Value))) - Decimals));
115
116 if( Decimals > 0 )
117 {
118 double d = pow(10., Decimals);
119
120 return( Value < 0.
121 ? -((int)(0.5 - Value * d)) / d
122 : ((int)(0.5 + Value * d)) / d
123 );
124 }
125 else
126 {
127 double d = pow(10., -Decimals);
128
129 return( Value < 0.
130 ? -((int)(0.5 - Value / d)) * d
131 : ((int)(0.5 + Value / d)) * d
132 );
133 }
134}
135
136
138// //
139// //
140// //
142
143//---------------------------------------------------------
144int SG_Get_Digit_Count(int Number)
145{
146 Number = abs(Number);
147
148 return( Number < 10 ? 1 : 1 + (int)log10((double)Number) );
149}
150
151
153// //
154// //
155// //
157
158//---------------------------------------------------------
159CSG_String SG_Get_Double_asString(double Number, int Width, int Precision, bool bScientific)
160{
161 if( bScientific )
162 {
163 if( Width > 0 && Precision >= 0 ) return( CSG_String::Format("%*.*e", Width, Precision, Number) );
164 if( Width > 0 ) return( CSG_String::Format("%*e" , Width , Number) );
165 if( Precision >= 0 ) return( CSG_String::Format("%.*e" , Precision, Number) );
166
167 return( CSG_String::Format("%e", Number) );
168 }
169 else
170 {
171 if( Width > 0 && Precision >= 0 ) return( CSG_String::Format("%*.*f", Width, Precision, Number) );
172 if( Width > 0 ) return( CSG_String::Format("%*f" , Width , Number) );
173 if( Precision >= 0 ) return( CSG_String::Format("%.*f" , Precision, Number) );
174
175 return( CSG_String::Format("%f", Number) );
176 }
177}
178
179
181// //
182// //
183// //
185
186//---------------------------------------------------------
187int SG_Compare_Int(const void *a, const void *b)
188{
189 if( *((int *)a) < *((int *)b) )
190 return( -1 );
191
192 if( *((int *)a) > *((int *)b) )
193 return( 1 );
194
195 return( 0 );
196}
197
198//---------------------------------------------------------
199int SG_Compare_Double(const void *a, const void *b)
200{
201 if( *((double *)a) < *((double *)b) )
202 return( -1 );
203
204 if( *((double *)a) > *((double *)b) )
205 return( 1 );
206
207 return( 0 );
208}
209
210//---------------------------------------------------------
211int SG_Compare_Char_Ptr(const void *a, const void *b)
212{
213 return( strcmp((const char *)a, (const char *)b) );
214}
215
216
218// //
219// //
220// //
222
223//---------------------------------------------------------
224double SG_Degree_To_Decimal(double Deg, double Min, double Sec)
225{
226 if( Deg < 0. ) { Min = -fabs(Min); Sec = -fabs(Sec); }
227
228 return( Deg + Min / 60. + Sec / 3600. );
229}
230
231//---------------------------------------------------------
232void SG_Decimal_To_Degree(double Value, double &Deg, double &Min, double &Sec)
233{
234 Sec = fmod(Value < 0. ? -Value : Value, 360.);
235
236 Deg = (int)Sec; Sec = 60. * (Sec - Deg);
237 Min = (int)Sec; Sec = 60. * (Sec - Min);
238
239 if( Value < 0. )
240 {
241 Deg = -Deg;
242 }
243}
244
245
247// //
248// //
249// //
251
252//---------------------------------------------------------
254{
255 Initialize();
256}
257
258//---------------------------------------------------------
260{
261 Initialize((unsigned)time(NULL));
262}
263
264//---------------------------------------------------------
265void CSG_Random::Initialize(unsigned int Value)
266{
267 srand(Value);
268}
269
270//---------------------------------------------------------
271// Uniform distributed pseudo-random numbers in the range from 0 to 1.
272//
274{
275 return( 1. * rand() / (double)RAND_MAX );
276}
277
278//---------------------------------------------------------
279// Uniform distributed pseudo-random numbers in the range from min to max.
280//
281double CSG_Random::Get_Uniform(double min, double max)
282{
283 return( min + (max - min) * rand() / (double)RAND_MAX );
284}
285
286//---------------------------------------------------------
287// Generating Gaussian pseudo-random numbers using
288// the polar form of the Box-Muller transformation.
289//
290// Box, G.E.P, Muller, M.E. (1958):
291// 'A note on the generation of random normal deviates',
292// Annals Math. Stat, V. 29, pp. 610-611
293//
294// Link: http://www.taygeta.com/random/gaussian.html
295//
296//---------------------------------------------------------
297double CSG_Random::Get_Gaussian(double mean, double stddev)
298{
299 double x1, x2, w;
300
301 do
302 {
303 x1 = 2. * Get_Uniform() - 1.;
304 x2 = 2. * Get_Uniform() - 1.;
305
306 w = x1 * x1 + x2 * x2;
307 }
308 while( w >= 1. );
309
310 w = sqrt((-2. * log(w)) / w);
311
312 return( mean + stddev * x1 * w );
313}
314
315
317// //
318// //
319// //
321
322//---------------------------------------------------------
327
329{
330 Create(bHoldValues);
331}
332
337
338CSG_Simple_Statistics::CSG_Simple_Statistics(double Mean, double StdDev, sLong Count)
339{
340 Create(Mean, StdDev, Count);
341}
342
344{
345 Create(Values, bHoldValues);
346}
347
348//---------------------------------------------------------
349bool CSG_Simple_Statistics::Create(bool bHoldValues)
350{
351 Invalidate();
352
353 m_Values.Create(bHoldValues ? sizeof(double) : 0, 0, TSG_Array_Growth::SG_ARRAY_GROWTH_1);
354
355 return( true );
356}
357
359{
360 m_bEvaluated = Statistics.m_bEvaluated;
361
362 m_nValues = Statistics.m_nValues;
363 m_Weights = Statistics.m_Weights;
364 m_Sum = Statistics.m_Sum;
365 m_Sum2 = Statistics.m_Sum2;
366
367 m_Minimum = Statistics.m_Minimum;
368 m_Maximum = Statistics.m_Maximum;
369 m_Range = Statistics.m_Range;
370 m_Mean = Statistics.m_Mean;
371 m_Variance = Statistics.m_Variance;
372 m_StdDev = Statistics.m_StdDev;
373
374 m_Kurtosis = Statistics.m_Kurtosis;
375 m_Skewness = Statistics.m_Skewness;
376
377 m_Gini = Statistics.m_Gini;
378
379 m_bSorted = Statistics.m_bSorted;
380 m_Values .Create(Statistics.m_Values);
381
382 return( true );
383}
384
385bool CSG_Simple_Statistics::Create(double Mean, double StdDev, sLong Count)
386{
387 Invalidate();
388
389 m_bEvaluated = 1;
390
391 m_Mean = Mean;
392 m_StdDev = StdDev;
393 m_Variance = StdDev*StdDev;
394 m_nValues = Count;
395 m_Weights = (double)Count;
396
399
400 m_Minimum = m_Mean - 1.5 * m_StdDev;
401 m_Maximum = m_Mean + 1.5 * m_StdDev;
403
404 return( true );
405}
406
407bool CSG_Simple_Statistics::Create(const CSG_Vector &Values, bool bHoldValues)
408{
409 if( Create(bHoldValues) )
410 {
411 for(sLong i=0; i<Values.Get_Size(); i++)
412 {
413 Add_Value(Values[i]);
414 }
415
416 return( true );
417 }
418
419 return( false );
420}
421
422//---------------------------------------------------------
424{
425 if( m_nValues < 1 || nValues < 1 || m_nValues == nValues )
426 {
427 return( false );
428 }
429
430 double Scale = nValues / (double)m_nValues;
431
432 m_Weights *= Scale;
433 m_Sum *= Scale;
434 m_Sum2 *= Scale;
435
436 m_nValues = nValues;
437
438 m_bEvaluated = 0;
439
440 m_Values.Destroy();
441
442 return( true );
443}
444
445//---------------------------------------------------------
447{
448 m_bEvaluated = 0;
449
450 m_nValues = 0;
451 m_Weights = 0.;
452 m_Sum = 0.;
453 m_Sum2 = 0.;
454
455 m_Minimum = 0.;
456 m_Maximum = 0.;
457 m_Range = 0.;
458 m_Mean = 0.;
459 m_Variance = 0.;
460 m_StdDev = 0.;
461
462 m_Kurtosis = 0.;
463 m_Skewness = 0.;
464
465 m_Gini = -1.;
466
467 m_bSorted = false;
468 m_Values .Destroy();
469}
470
471//---------------------------------------------------------
473{
474 _Evaluate();
475
476 return( is_Evaluated() > 0 );
477}
478
479//---------------------------------------------------------
481{
482 if( Statistics.m_nValues <= 0 )
483 {
484 return;
485 }
486
487 if( m_nValues == 0 )
488 {
489 Create(Statistics);
490
491 return;
492 }
493
494 //--------------------------------------------------------
495 if( (sLong)m_Values.Get_Size() == m_nValues && (sLong)Statistics.m_Values.Get_Size() == Statistics.m_nValues && m_Values.Set_Array((size_t)(m_nValues + Statistics.m_nValues)) )
496 {
497 for(sLong i=0, j=m_nValues; i<Statistics.m_nValues; i++, j++)
498 {
499 ((double *)m_Values.Get_Array())[j] = Statistics.Get_Value(i);
500 }
501 }
502 else
503 {
504 m_Values.Destroy();
505 }
506
507 m_nValues += Statistics.m_nValues;
508 m_Weights += Statistics.m_Weights;
509 m_Sum += Statistics.m_Sum;
510 m_Sum2 += Statistics.m_Sum2;
511
512 if( m_Minimum > Statistics.m_Minimum )
513 m_Minimum = Statistics.m_Minimum;
514
515 if( m_Maximum < Statistics.m_Maximum )
516 m_Maximum = Statistics.m_Maximum;
517
518 m_Kurtosis = 0.;
519 m_Skewness = 0.;
520
521 m_bEvaluated = 0;
522 m_bSorted = false;
523}
524
525//---------------------------------------------------------
526void CSG_Simple_Statistics::Add_Value(double Value, double Weight)
527{
528 if( m_nValues == 0 )
529 {
530 m_Minimum = m_Maximum = Value;
531 }
532 else if( m_Minimum > Value )
533 {
534 m_Minimum = Value;
535 }
536 else if( m_Maximum < Value )
537 {
538 m_Maximum = Value;
539 }
540
541 if( Weight )
542 {
543 m_Weights += fabs(Weight);
544 m_Sum += Weight * Value;
545 m_Sum2 += Weight * Value*Value;
546
547 m_bEvaluated = 0;
548
549 if( m_Values.Get_Value_Size() > 0 && m_Values.Inc_Array() )
550 {
551 m_bSorted = false;
552
553 ((double *)m_Values.Get_Array())[m_nValues] = Value;
554 }
555
556 m_nValues++;
557 }
558}
559
560//---------------------------------------------------------
562{
563 if( m_bEvaluated == 0 && m_Weights > 0. )
564 {
565 m_bEvaluated = 1;
566
570 m_StdDev = m_Variance > 0. ? sqrt(m_Variance) : 0.;
571 }
572
573 //-----------------------------------------------------
574 if( m_bEvaluated == 1 && Level > 1 )
575 {
576 m_bEvaluated = 2;
577
578 m_Kurtosis = 0.;
579 m_Skewness = 0.;
580
581 if( Get_StdDev() > 0. && m_Values.Get_Size() > 0 )
582 {
583 for(sLong i=0; i<Get_Count(); i++)
584 {
585 double d = (Get_Value(i) - Get_Mean()) / Get_StdDev();
586
587 m_Kurtosis += d*d*d*d;
588 m_Skewness += d*d*d;
589 }
590
593 // m_Skewness *= Get_Count() / ((Get_Count() - 1) * (Get_Count() - 2));
594 }
595 }
596}
597
598//---------------------------------------------------------
606{
607 return( Get_StdDev() != 0. ? (Get_Mean() - Get_Median()) / Get_StdDev() : 0. );
608}
609
610//---------------------------------------------------------
618{
619 if( m_Values.Get_Size() < 1 )
620 {
621 return( m_Mean );
622 }
623
624 //-----------------------------------------------------
625 if( !m_bSorted )
626 {
627 qsort(m_Values.Get_Array(), m_Values.Get_Size(), sizeof(double), SG_Compare_Double);
628
629 m_bSorted = true;
630 }
631
632 if( Quantile <= 0. || m_Values.Get_Size() == 1 )
633 {
634 return( Get_Values()[0] );
635 }
636
637 if( Quantile >= 1. )
638 {
639 return( Get_Values()[m_Values.Get_Size() - 1] );
640 }
641
642 //-----------------------------------------------------
643 double r = Quantile * (m_Values.Get_Size() - 1);
644
645 sLong i = (sLong)r; r -= i;
646
647 return( r == 0. ? Get_Values()[i] : ((1. - r) * Get_Values()[i] + r * Get_Values()[i + 1]) );
648}
649
650//---------------------------------------------------------
658{
659 return( Get_Quantile(Percentile / 100.) );
660}
661
662//---------------------------------------------------------
669{
670 return( Get_Quantile(0.5) );
671}
672
673//---------------------------------------------------------
681{
682 if( m_Gini < 0. && m_Values.Get_Size() > 1 )
683 {
684 if( !m_bSorted )
685 {
686 qsort(m_Values.Get_Array(), m_Values.Get_Size(), sizeof(double), SG_Compare_Double);
687
688 m_bSorted = true;
689 }
690
691 m_Gini = 0.;
692
693 for(sLong i=0; i<Get_Count(); i++)
694 {
695 m_Gini += (i + 1.) * Get_Value(i);
696 }
697
698 m_Gini = 2. * m_Gini / (Get_Count() * Get_Sum()) - (Get_Count() + 1.) / Get_Count();
699 }
700
701 return( m_Gini );
702}
703
704//---------------------------------------------------------
710{
711 if( m_Values.Get_Size() == 0 ) { return( -1 ); }
712
713 size_t Index = 0;
714 double Value = Get_Values()[Index];
715
716 for(size_t i=1; i<(size_t)m_Values.Get_Size(); i++)
717 {
718 if( Value > Get_Values()[i] )
719 {
720 Index = i;
721 Value = Get_Values()[i];
722 }
723 }
724
725 return( (sLong)Index );
726}
727
728//---------------------------------------------------------
734{
735 if( m_Values.Get_Size() == 0 ) { return( -1 ); }
736
737 size_t Index = 0;
738 double Value = Get_Values()[Index];
739
740 for(size_t i=1; i<(size_t)m_Values.Get_Size(); i++)
741 {
742 if( Value < Get_Values()[i] )
743 {
744 Index = i;
745 Value = Get_Values()[i];
746 }
747 }
748
749 return( (sLong)Index );
750}
751
752//---------------------------------------------------------
757sLong CSG_Simple_Statistics::Get_nValues_Above(double Threshold, bool bEquals)
758{
759 if( m_Values.Get_Size() == 0 ) { return( -1 ); }
760
761 sLong n = 0;
762
763 for(sLong i=0; i<Get_Count(); i++)
764 {
765 if( (bEquals && Get_Value(i) >= Threshold) || Get_Value(i) > Threshold )
766 {
767 n++;
768 }
769 }
770
771 return( n );
772}
773
774//---------------------------------------------------------
779sLong CSG_Simple_Statistics::Get_nValues_Below(double Threshold, bool bEquals)
780{
781 if( m_Values.Get_Size() == 0 ) { return( -1 ); }
782
783 sLong n = 0;
784
785 for(sLong i=0; i<Get_Count(); i++)
786 {
787 if( (bEquals && Get_Value(i) <= Threshold) || Get_Value(i) < Threshold )
788 {
789 n++;
790 }
791 }
792
793 return( n );
794}
795
796
798// //
799// //
800// //
802
803//---------------------------------------------------------
805{
806 int Index = 0;
807
808 bWeighted = bWeighted && m_bWeights;
809
810 for(int i=1; i<Get_Count(); i++)
811 {
812 if( bWeighted )
813 {
814 if( m_Weight[i] > m_Weight[Index] )
815 {
816 Index = i;
817 }
818 }
819 else
820 {
821 if( m_Count[i] > m_Count[Index] )
822 {
823 Index = i;
824 }
825 }
826 }
827
828 return( Index );
829}
830
831//---------------------------------------------------------
833{
834 int Index = 0;
835
836 bWeighted = bWeighted && m_bWeights;
837
838 for(int i=1; i<Get_Count(); i++)
839 {
840 if( bWeighted )
841 {
842 if( m_Weight[i] < m_Weight[Index] )
843 {
844 Index = i;
845 }
846 }
847 else
848 {
849 if( m_Count[i] < m_Count[Index] )
850 {
851 Index = i;
852 }
853 }
854 }
855
856 return( Index );
857}
858
859
861// //
863
864//---------------------------------------------------------
866{
867 m_bWeights = bWeights;
868
869 m_Count.Destroy();
870 m_Value.Destroy();
871}
872
873//---------------------------------------------------------
874void CSG_Unique_Number_Statistics::Add_Value(double Value, double Weight)
875{
876 for(int i=0; i<Get_Count(); i++)
877 {
878 if( Value == m_Value[i] )
879 {
880 m_Count[i]++;
881
882 if( m_bWeights && Weight > 0. )
883 {
884 m_Weight[i] += Weight;
885 }
886
887 return;
888 }
889 }
890
891 m_Count.Add(1);
892 m_Value.Add_Row(Value);
893
894 if( m_bWeights && Weight > 0. )
895 {
896 m_Weight.Add_Row(Weight);
897 }
898}
899
900//---------------------------------------------------------
902{
903 for(int i=0; i<Get_Count(); i++)
904 {
905 if( Value == m_Value[i] )
906 {
907 return( i );
908 }
909 }
910
911 return( -1 );
912}
913
914
916// //
918
919//---------------------------------------------------------
921{
922 m_bWeights = bWeights;
923
924 m_Count.Destroy();
925 m_Value.Clear();
926}
927
928//---------------------------------------------------------
930{
931 for(int i=0; i<Get_Count(); i++)
932 {
933 if( Value.Cmp(m_Value[i]) == 0 )
934 {
935 m_Count[i]++;
936
937 if( m_bWeights && Weight > 0. )
938 {
939 m_Weight[i] += Weight;
940 }
941
942 return;
943 }
944 }
945
946 m_Count.Add(1);
947 m_Value.Add(Value);
948
949 if( m_bWeights && Weight > 0. )
950 {
951 m_Weight.Add_Row(Weight);
952 }
953}
954
955//---------------------------------------------------------
957{
958 for(int i=0; i<Get_Count(); i++)
959 {
960 if( Value.Cmp(m_Value[i]) == 0 )
961 {
962 return( i );
963 }
964 }
965
966 return( -1 );
967}
968
969
971// //
972// //
973// //
975
976//---------------------------------------------------------
978{
979 m_pTable = new CSG_Table;
980
981 Create(Type);
982}
983
984//---------------------------------------------------------
986{
987 delete(m_pTable);
988}
989
990//---------------------------------------------------------
992{
993 m_pTable->Destroy();
994
995 m_pTable->Add_Field("VALUE", Type);
996 m_pTable->Add_Field("COUNT", SG_DATATYPE_ULong);
997}
998
999//---------------------------------------------------------
1001{
1002 m_pTable->Del_Records();
1003}
1004
1005//---------------------------------------------------------
1007{
1008 return( m_pTable->Get_Field_Type(0) );
1009}
1010
1011
1013// //
1015
1016//---------------------------------------------------------
1018{
1019 CSG_Table_Record *pRecord = m_pTable->Find_Record(0, Value);
1020
1021 return( pRecord ? (int)pRecord->Get_Index() : -1);
1022}
1023
1024//---------------------------------------------------------
1026{
1027 CSG_Table_Record *pRecord = m_pTable->Find_Record(0, Value);
1028
1029 return( pRecord ? (int)pRecord->Get_Index() : -1);
1030}
1031
1032//---------------------------------------------------------
1034{
1035 CSG_Table_Record *pRecord = m_pTable->Find_Record(0, Value);
1036
1037 return( pRecord ? (int)pRecord->Get_Index() : -1);
1038}
1039
1040
1042// //
1044
1045//---------------------------------------------------------
1047{
1048 CSG_Table_Record *pRecord = m_pTable->Find_Record(0, Value);
1049
1050 if( !pRecord )
1051 {
1052 (pRecord = m_pTable->Add_Record())->Set_Value(0, Value);
1053
1054 if( m_pTable->Get_Count() > 16 || m_pTable->is_Indexed() )
1055 {
1056 Sort();
1057 }
1058 }
1059
1060 pRecord->Add_Value(1, 1);
1061
1062 return( (int)(m_pTable->Get_Count() - 1) );
1063}
1064
1065//---------------------------------------------------------
1067{
1068 CSG_Table_Record *pRecord = m_pTable->Find_Record(0, Value);
1069
1070 if( !pRecord )
1071 {
1072 (pRecord = m_pTable->Add_Record())->Set_Value(0, Value);
1073
1074 if( m_pTable->Get_Count() > 16 || m_pTable->is_Indexed() )
1075 {
1076 Sort();
1077 }
1078 }
1079
1080 pRecord->Add_Value(1, 1);
1081
1082 return( (int)(m_pTable->Get_Count() - 1) );
1083}
1084
1085//---------------------------------------------------------
1087{
1088 CSG_Table_Record *pRecord = m_pTable->Find_Record(0, Value);
1089
1090 if( !pRecord )
1091 {
1092 (pRecord = m_pTable->Add_Record())->Set_Value(0, Value);
1093
1094 if( m_pTable->Get_Count() > 16 || m_pTable->is_Indexed() )
1095 {
1096 Sort();
1097 }
1098 }
1099
1100 pRecord->Add_Value(1, 1);
1101
1102 return( (int)(m_pTable->Get_Count() - 1) );
1103}
1104
1105//---------------------------------------------------------
1106// sort categories ascending
1108{
1109 CSG_Table Table(*m_pTable);
1110
1111 if( Table.Set_Index(0, TABLE_INDEX_Ascending) )
1112 {
1113 for(sLong i=0; i<m_pTable->Get_Count(); i++)
1114 {
1115 m_pTable->Set_Record(i, Table.Get_Record_byIndex(i));
1116 }
1117
1118 return( m_pTable->Set_Index(0, TABLE_INDEX_Ascending) );
1119 }
1120
1121 return( false );
1122}
1123
1124
1126// //
1128
1129//---------------------------------------------------------
1130// returns the number of categories.
1132{
1133 return( (int)m_pTable->Get_Count() );
1134}
1135
1136//---------------------------------------------------------
1137// returns the number of observations for the i'th category.
1139{
1140 CSG_Table_Record *pRecord = m_pTable->Get_Record(i);
1141
1142 return( pRecord ? pRecord->asInt(1) : 0 );
1143}
1144
1145//---------------------------------------------------------
1147{
1148 CSG_Table_Record *pRecord = m_pTable->Get_Record(i);
1149
1150 return( pRecord ? pRecord->asInt(0) : 0 );
1151}
1152
1153//---------------------------------------------------------
1155{
1156 CSG_Table_Record *pRecord = m_pTable->Get_Record(i);
1157
1158 return( pRecord ? pRecord->asDouble(0) : 0 );
1159}
1160
1161//---------------------------------------------------------
1163{
1164 CSG_Table_Record *pRecord = m_pTable->Get_Record(i);
1165
1166 return( pRecord ? pRecord->asString(0) : SG_T("") );
1167}
1168
1169
1171// //
1173
1174//---------------------------------------------------------
1176{
1177 if( m_pTable->Get_Count() > 0 )
1178 {
1179 int Index = 0, Count = m_pTable->Get_Record_byIndex(0)->asInt(1);
1180
1181 for(int i=1; i<m_pTable->Get_Count(); i++)
1182 {
1183 if( Count < m_pTable->Get_Record_byIndex(i)->asInt(1) )
1184 {
1185 Index = i; Count = m_pTable->Get_Record_byIndex(i)->asInt(1);
1186 }
1187 }
1188
1189 return( Index );
1190 }
1191
1192 return( -1 );
1193}
1194
1195//---------------------------------------------------------
1197{
1198 if( m_pTable->Get_Count() > 0 )
1199 {
1200 int Index = 0, Count = m_pTable->Get_Record_byIndex(0)->asInt(1);
1201
1202 for(int i=1; i<m_pTable->Get_Count(); i++)
1203 {
1204 if( Count > m_pTable->Get_Record_byIndex(i)->asInt(1) )
1205 {
1206 Index = i; Count = m_pTable->Get_Record_byIndex(i)->asInt(1);
1207 }
1208 }
1209
1210 return( Index );
1211 }
1212
1213 return( -1 );
1214}
1215
1216
1218// //
1219// //
1220// //
1222
1223//---------------------------------------------------------
1225{
1226 _On_Construction();
1227}
1228
1229//---------------------------------------------------------
1231{
1232 _On_Construction();
1233
1234 Create(Histogram);
1235}
1236
1237//---------------------------------------------------------
1238CSG_Histogram::CSG_Histogram(size_t nClasses, double Minimum, double Maximum)
1239{
1240 _On_Construction();
1241
1242 Create(nClasses, Minimum, Maximum);
1243}
1244
1245//---------------------------------------------------------
1246CSG_Histogram::CSG_Histogram(size_t nClasses, const CSG_Vector &Values, double Minimum, double Maximum, size_t maxSamples)
1247{
1248 _On_Construction();
1249
1250 Create(nClasses, Values, Minimum, Maximum, maxSamples);
1251}
1252
1253//---------------------------------------------------------
1254CSG_Histogram::CSG_Histogram(size_t nClasses, CSG_Table *pTable, int Field, double Minimum, double Maximum, size_t maxSamples, int Normalize, double Scale)
1255{
1256 _On_Construction();
1257
1258 Create(nClasses, pTable, Field, Minimum, Maximum, maxSamples, Normalize, Scale);
1259}
1260
1261//---------------------------------------------------------
1262CSG_Histogram::CSG_Histogram(size_t nClasses, class CSG_Grid *pGrid, double Minimum, double Maximum, size_t maxSamples)
1263{
1264 _On_Construction();
1265
1266 Create(nClasses, pGrid, Minimum, Maximum, maxSamples);
1267}
1268
1269//---------------------------------------------------------
1270CSG_Histogram::CSG_Histogram(size_t nClasses, CSG_Grids *pGrids, double Minimum, double Maximum, size_t maxSamples)
1271{
1272 _On_Construction();
1273
1274 Create(nClasses, pGrids, Minimum, Maximum, maxSamples);
1275}
1276
1277//---------------------------------------------------------
1279{
1280 Destroy();
1281}
1282
1283//---------------------------------------------------------
1285{
1286 m_Statistics.Create();
1287
1288 SG_FREE_SAFE(m_Elements );
1289 SG_FREE_SAFE(m_Cumulative);
1290
1291 _On_Construction();
1292
1293 return( true );
1294}
1295
1296
1298// //
1300
1301//---------------------------------------------------------
1302void CSG_Histogram::_On_Construction(void)
1303{
1304 m_nClasses = 0;
1305 m_Elements = NULL;
1306 m_Cumulative = NULL;
1307 m_Minimum = 0.;
1308 m_Maximum = 0.;
1309 m_ClassWidth = 1.;
1310}
1311
1312//---------------------------------------------------------
1313bool CSG_Histogram::_Create(size_t nClasses, double Minimum, double Maximum)
1314{
1315 if( nClasses > 0 && Minimum < Maximum )
1316 {
1317 Destroy();
1318
1319 m_Elements = (size_t *)SG_Calloc(nClasses, sizeof(size_t));
1320 m_Cumulative = (size_t *)SG_Calloc(nClasses, sizeof(size_t));
1321
1322 if( m_Elements && m_Cumulative )
1323 {
1324 m_nClasses = nClasses;
1325 m_Minimum = Minimum;
1326 m_Maximum = Maximum;
1327 m_ClassWidth = (Maximum - Minimum) / (double)m_nClasses;
1328
1329 return( true );
1330 }
1331 }
1332
1333 Destroy();
1334
1335 return( false );
1336}
1337
1338//---------------------------------------------------------
1340{
1341 m_Statistics += Value;
1342
1343 if( m_nClasses > 0 && m_Minimum <= Value && Value <= m_Maximum )
1344 {
1345 size_t Class = (size_t)((Value - m_Minimum) / m_ClassWidth);
1346
1347 if( Class >= m_nClasses )
1348 {
1349 Class = m_nClasses - 1;
1350 }
1351
1352 m_Elements[Class]++;
1353 }
1354}
1355
1356//---------------------------------------------------------
1358{
1359 if( m_nClasses > 0 && Scale > 0. )
1360 {
1361 m_Statistics.Set_Count((sLong)(Scale * Get_Element_Count()));
1362
1363 for(size_t i=0; i<m_nClasses; i++)
1364 {
1365 m_Elements[i] = (size_t)(Scale * m_Elements[i]);
1366 }
1367
1368 return( Update() );
1369 }
1370
1371 return( false );
1372}
1373
1374//---------------------------------------------------------
1376{
1377 if( m_nClasses > 0 )
1378 {
1379 m_Statistics.Get_Mean(); // _Evaluate()
1380
1381 m_nMaximum = m_Cumulative[0] = m_Elements[0];
1382
1383 for(size_t i=1; i<m_nClasses; i++)
1384 {
1385 m_Cumulative[i] = m_Cumulative[i - 1] + m_Elements[i];
1386
1387 if( m_nMaximum < m_Elements[i] )
1388 {
1389 m_nMaximum = m_Elements[i];
1390 }
1391 }
1392
1393 return( Get_Element_Count() > 0 );
1394 }
1395
1396 return( false );
1397}
1398
1399//---------------------------------------------------------
1400bool CSG_Histogram::_Update(sLong nElements)
1401{
1402 if( nElements > 0 && m_Statistics.Get_Count() > 0 )
1403 {
1404 double Scale = (double)nElements / (double)m_Statistics.Get_Count();
1405
1406 m_Statistics.Create(m_Statistics.Get_Mean(), m_Statistics.Get_StdDev(), nElements);
1407
1408 for(size_t i=1; i<m_nClasses; i++)
1409 {
1410 m_Elements[i] = (size_t)(0.5 + Scale * m_Elements[i]);
1411 }
1412 }
1413
1414 return( Update() );
1415}
1416
1418{
1419 if( m_nClasses == 0
1420 || m_nClasses != Histogram.m_nClasses
1421 || m_Minimum != Histogram.m_Minimum
1422 || m_Maximum != Histogram.m_Maximum )
1423 {
1424 return( false );
1425 }
1426
1427 m_Statistics += Histogram.m_Statistics;
1428
1429 for( size_t i=0; i<m_nClasses; i++ )
1430 {
1431 m_Elements[i] += Histogram.m_Elements[i];
1432 }
1433
1434 return( true );
1435}
1436
1437//---------------------------------------------------------
1441double CSG_Histogram::Get_Quantile(double Quantile) const
1442{
1443 if( m_nClasses < 2 ) { return( 0. ); }
1444
1445 if( Quantile <= 0. ) { return( m_Minimum ); }
1446 if( Quantile >= 1. ) { return( m_Maximum ); }
1447
1448 size_t n = (size_t)(Quantile * Get_Element_Count()); // number of elements
1449
1450 for(size_t i=0, n0=0; i<m_nClasses; n0=m_Cumulative[i++])
1451 {
1452 if( n < m_Cumulative[i] )
1453 {
1454 if( m_Cumulative[i] >= n0 )
1455 {
1456 return( Get_Center(i) );
1457 }
1458
1459 double d = (n - n0) / (double)(m_Cumulative[i] - n0);
1460
1461 return( Get_Break(i) + d * m_ClassWidth );
1462 }
1463 else if( n == m_Cumulative[i] )
1464 {
1465 return( Get_Break(i + 1) );
1466 }
1467 }
1468
1469 return( m_Maximum );
1470}
1471
1472//---------------------------------------------------------
1476double CSG_Histogram::Get_Percentile(double Percentile) const
1477{
1478 return( Get_Quantile(Percentile / 100.) );
1479}
1480
1481//---------------------------------------------------------
1485double CSG_Histogram::Get_Quantile_Value(double Value) const
1486{
1487 if( m_nClasses < 2 ) { return( 0. ); }
1488
1489 if( Value <= m_Minimum ) { return( 0. ); }
1490 if( Value >= m_Maximum ) { return( 1. ); }
1491
1492 size_t Class = (size_t)(m_nClasses * (Value - m_Minimum) / (m_Maximum - m_Minimum));
1493
1494 if( Class >= m_nClasses )
1495 {
1496 return( 1. );
1497 }
1498
1499 if( Class < 1 )
1500 {
1501 double dq = m_Cumulative[Class] / (double)Get_Element_Count();
1502
1503 return( dq * (Value - m_Minimum) / m_ClassWidth );
1504 }
1505
1506 double q0 = m_Cumulative[Class - 1] / (double)Get_Element_Count();
1507 double dq = (m_Cumulative[Class ] / (double)Get_Element_Count()) - q0;
1508
1509 return( q0 + dq * (Value - Get_Break(Class)) / m_ClassWidth );
1510}
1511
1512//---------------------------------------------------------
1516double CSG_Histogram::Get_Percentile_Value(double Value) const
1517{
1518 return( Get_Quantile_Value(Value) * 100. );
1519}
1520
1521
1523// //
1525
1526//---------------------------------------------------------
1528{
1529 if( !_Create(Histogram.m_nClasses, Histogram.m_Minimum, Histogram.m_Maximum) )
1530 {
1531 return( false );
1532 }
1533
1534 m_Statistics = Histogram.m_Statistics;
1535 m_ClassWidth = Histogram.m_ClassWidth;
1536 m_nMaximum = Histogram.m_nMaximum ;
1537
1538 for(size_t i=0; i<m_nClasses; i++)
1539 {
1540 m_Cumulative[i] = Histogram.m_Cumulative[i];
1541 m_Elements [i] = Histogram.m_Elements [i];
1542 }
1543
1544 return( true );
1545}
1546
1547//---------------------------------------------------------
1548bool CSG_Histogram::Create(size_t nClasses, double Minimum, double Maximum)
1549{
1550 return( _Create(nClasses, Minimum, Maximum) );
1551}
1552
1553//---------------------------------------------------------
1554bool CSG_Histogram::Create(size_t nClasses, const CSG_Vector &Values, double Minimum, double Maximum, size_t maxSamples)
1555{
1556 if( Minimum >= Maximum )
1557 {
1558 CSG_Simple_Statistics Statistics(Values);
1559
1560 Minimum = Statistics.Get_Minimum();
1561 Maximum = Statistics.Get_Maximum();
1562 }
1563
1564 if( !_Create(nClasses, Minimum, Maximum) )
1565 {
1566 return( false );
1567 }
1568
1569 //-----------------------------------------------------
1570 if( maxSamples > 0 && maxSamples < (size_t)Values.Get_N() )
1571 {
1572 double d = (double)Values.Get_N() / (double)maxSamples;
1573
1574 for(double i=0; i<(double)Values.Get_N(); i+=d)
1575 {
1576 Add_Value(Values[(sLong)i]);
1577 }
1578
1579 d = (double)m_Statistics.Get_Count() / (double)maxSamples;
1580
1581 return( _Update(d < 1. ? (int)(d * (double)Values.Get_N()) : Values.Get_N()) );
1582 }
1583
1584 //-----------------------------------------------------
1585 for(int i=0; i<Values.Get_N(); i++)
1586 {
1587 Add_Value(Values[i]);
1588 }
1589
1590 return( Update() );
1591}
1592
1593//---------------------------------------------------------
1594bool CSG_Histogram::Create(size_t nClasses, CSG_Table *pTable, int Field, double Minimum, double Maximum, size_t maxSamples, int Normalize, double Scale)
1595{
1596 if( !pTable || Field < 0 || Field >= pTable->Get_Field_Count() || !_Create(nClasses,
1598 Minimum < Maximum ? Maximum : pTable->Get_Maximum(Field)) )
1599 {
1600 return( false );
1601 }
1602
1603 //-----------------------------------------------------
1604 if( maxSamples > 0 && maxSamples < (size_t)pTable->Get_Count() )
1605 {
1606 double Value, Dividend, d = (double)pTable->Get_Count() / (double)maxSamples;
1607
1608 for(double i=0; i<(double)pTable->Get_Count(); i+=d)
1609 {
1610 if( pTable->Get_Value((sLong)i, Field, Value) )
1611 {
1612 if( Normalize < 0 )
1613 {
1614 Add_Value(Value);
1615 }
1616 else if( pTable->Get_Value((sLong)i, Normalize, Dividend) && Dividend != 0. )
1617 {
1618 Add_Value(Scale * Value / Dividend);
1619 }
1620 }
1621 }
1622
1623 d = (double)m_Statistics.Get_Count() / (double)maxSamples;
1624
1625 return( _Update(d < 1. ? (int)(d * pTable->Get_Count()) : pTable->Get_Count()) );
1626 }
1627
1628 //-----------------------------------------------------
1629 for(sLong i=0; i<pTable->Get_Count(); i++)
1630 {
1631 double Value, Dividend;
1632
1633 if( pTable->Get_Value(i, Field, Value) )
1634 {
1635 if( Normalize < 0 )
1636 {
1637 Add_Value(Value);
1638 }
1639 else if( pTable->Get_Value((sLong)i, Normalize, Dividend) && Dividend != 0. )
1640 {
1641 Add_Value(Scale * Value / Dividend);
1642 }
1643 }
1644 }
1645
1646 return( Update() );
1647}
1648
1649//---------------------------------------------------------
1650bool CSG_Histogram::Create(size_t nClasses, CSG_Grid *pGrid, double Minimum, double Maximum, size_t maxSamples)
1651{
1652 if( !pGrid || !_Create(nClasses,
1655 {
1656 return( false );
1657 }
1658
1659 //-----------------------------------------------------
1660 if( maxSamples > 0 && (sLong)maxSamples < pGrid->Get_NCells() )
1661 {
1662 double d = (double)pGrid->Get_NCells() / (double)maxSamples;
1663
1664 for(double i=0; i<(double)pGrid->Get_NCells(); i+=d)
1665 {
1666 if( !pGrid->is_NoData((sLong)i) )
1667 {
1668 Add_Value(pGrid->asDouble((sLong)i));
1669 }
1670 }
1671
1672 d = (double)m_Statistics.Get_Count() / (double)maxSamples;
1673
1674 return( _Update(d < 1. ? (sLong)(d * (double)pGrid->Get_NCells()) : pGrid->Get_NCells()) );
1675 }
1676
1677 //-----------------------------------------------------
1678 for(sLong i=0; i<pGrid->Get_NCells(); i++)
1679 {
1680 if( !pGrid->is_NoData(i) )
1681 {
1682 Add_Value(pGrid->asDouble(i));
1683 }
1684 }
1685
1686 return( Update() );
1687}
1688
1689//---------------------------------------------------------
1690bool CSG_Histogram::Create(size_t nClasses, CSG_Grids *pGrids, double Minimum, double Maximum, size_t maxSamples)
1691{
1692 if( !pGrids || !_Create(nClasses,
1695 {
1696 return( false );
1697 }
1698
1699 //-----------------------------------------------------
1700 if( maxSamples > 0 && (sLong)maxSamples < pGrids->Get_NCells() )
1701 {
1702 double d = (double)pGrids->Get_NCells() / (double)maxSamples;
1703
1704 for(double i=0; i<(double)pGrids->Get_NCells(); i+=d)
1705 {
1706 if( !pGrids->is_NoData((sLong)i) )
1707 {
1708 Add_Value(pGrids->asDouble((sLong)i));
1709 }
1710 }
1711
1712 d = (double)m_Statistics.Get_Count() / (double)maxSamples;
1713
1714 return( _Update(d < 1. ? (sLong)(d * (double)pGrids->Get_NCells()) : pGrids->Get_NCells()) );
1715 }
1716
1717 //-----------------------------------------------------
1718 for(sLong i=0; i<pGrids->Get_NCells(); i++)
1719 {
1720 if( !pGrids->is_NoData(i) )
1721 {
1722 Add_Value(pGrids->asDouble(i));
1723 }
1724 }
1725
1726 return( Update() );
1727}
1728
1729
1731// //
1733
1734//---------------------------------------------------------
1736{
1737 Create(Histogram);
1738
1739 return( *this );
1740}
1741
1742
1744// //
1745// //
1746// //
1748
1749//---------------------------------------------------------
1752
1753//---------------------------------------------------------
1756
1757//---------------------------------------------------------
1758CSG_Natural_Breaks::CSG_Natural_Breaks(CSG_Table *pTable, int Field, int nClasses, int Histogram)
1759{
1760 Create(pTable, Field, nClasses, Histogram);
1761}
1762
1763//---------------------------------------------------------
1764CSG_Natural_Breaks::CSG_Natural_Breaks(CSG_Grid *pGrid, int nClasses, int Histogram)
1765{
1766 Create(pGrid, nClasses, Histogram);
1767}
1768
1769//---------------------------------------------------------
1770CSG_Natural_Breaks::CSG_Natural_Breaks(CSG_Grids *pGrids, int nClasses, int Histogram)
1771{
1772 Create(pGrids, nClasses, Histogram);
1773}
1774
1775//---------------------------------------------------------
1776CSG_Natural_Breaks::CSG_Natural_Breaks(const CSG_Vector &Values, int nClasses, int Histogram)
1777{
1778 Create(Values, nClasses, Histogram);
1779}
1780
1781
1783// //
1785
1786//---------------------------------------------------------
1787bool CSG_Natural_Breaks::Create(CSG_Table *pTable, int Field, int nClasses, int Histogram)
1788{
1789 bool bResult = false;
1790
1791 if( Histogram > 0 )
1792 {
1793 bResult = m_Histogram.Create(Histogram, pTable, Field) && _Histogram(nClasses);
1794 }
1795 else if( Field >= 0 && Field < pTable->Get_Field_Count() )
1796 {
1797 for(sLong i=0; i<pTable->Get_Count(); i++)
1798 {
1799 double Value;
1800
1801 if( pTable->Get_Value(i, Field, Value) )
1802 {
1803 m_Values.Add_Row(Value);
1804 }
1805 }
1806
1807 bResult = m_Values.Sort() && _Calculate(nClasses);
1808
1809 m_Values.Destroy();
1810 }
1811
1812 return( bResult );
1813}
1814
1815//---------------------------------------------------------
1816bool CSG_Natural_Breaks::Create(CSG_Grid *pGrid, int nClasses, int Histogram)
1817{
1818 bool bResult = false;
1819
1820 if( Histogram > 0 )
1821 {
1822 bResult = m_Histogram.Create(Histogram, pGrid) && _Histogram(nClasses);
1823 }
1824 else
1825 {
1826 for(sLong i=0; i<pGrid->Get_NCells(); i++)
1827 {
1828 if( !pGrid->is_NoData(i) )
1829 {
1830 m_Values.Add_Row(pGrid->asDouble(i));
1831 }
1832 }
1833
1834 bResult = m_Values.Sort() && _Calculate(nClasses);
1835
1836 m_Values.Destroy();
1837 }
1838
1839 return( bResult );
1840}
1841
1842//---------------------------------------------------------
1843bool CSG_Natural_Breaks::Create(CSG_Grids *pGrids, int nClasses, int Histogram)
1844{
1845 bool bResult = false;
1846
1847 if( Histogram > 0 )
1848 {
1849 bResult = m_Histogram.Create(Histogram, pGrids) && _Histogram(nClasses);
1850 }
1851 else
1852 {
1853 for(sLong i=0; i<pGrids->Get_NCells(); i++)
1854 {
1855 if( !pGrids->is_NoData(i) )
1856 {
1857 m_Values.Add_Row(pGrids->asDouble(i));
1858 }
1859 }
1860
1861 bResult = m_Values.Sort() && _Calculate(nClasses);
1862
1863 m_Values.Destroy();
1864 }
1865
1866 return( bResult );
1867}
1868
1869//---------------------------------------------------------
1870bool CSG_Natural_Breaks::Create(const CSG_Vector &Values, int nClasses, int Histogram)
1871{
1872 bool bResult = false;
1873
1874 if( Histogram > 0 )
1875 {
1876 bResult = m_Histogram.Create(Histogram, Values) && _Histogram(nClasses);
1877 }
1878 else
1879 {
1880 bResult = m_Values.Create(Values) && m_Values.Sort() && _Calculate(nClasses);
1881
1882 m_Values.Destroy();
1883 }
1884
1885 return( bResult );
1886}
1887
1888
1890// //
1892
1893//---------------------------------------------------------
1894bool CSG_Natural_Breaks::_Histogram(int nClasses)
1895{
1896 if( _Calculate(nClasses) )
1897 {
1898 double d = (double)m_Histogram.Get_Class_Count() / m_Histogram.Get_Cumulative((int)(m_Histogram.Get_Class_Count() - 1));
1899
1900 m_Breaks[0] = m_Histogram.Get_Break(0);
1901
1902 for(int i=1; i<Get_Count(); i++)
1903 {
1904 m_Breaks[i] = m_Histogram.Get_Value(m_Breaks[i] * d);
1905 }
1906
1907 m_Breaks[nClasses] = m_Histogram.Get_Break((int)m_Histogram.Get_Class_Count());
1908
1909 m_Histogram.Destroy();
1910
1911 return( true );
1912 }
1913
1914 m_Histogram.Destroy();
1915
1916 return( false );
1917}
1918
1919//---------------------------------------------------------
1920inline double CSG_Natural_Breaks::_Get_Value(int i)
1921{
1922 if( m_Histogram.Get_Class_Count() > 0 )
1923 {
1924 return( (double)m_Histogram.Get_Cumulative(i) );
1925 }
1926
1927 return( m_Values[i] );
1928}
1929
1930//---------------------------------------------------------
1931bool CSG_Natural_Breaks::_Calculate(int nClasses)
1932{
1933 if( m_Histogram.Get_Class_Count() == 0 && m_Values.Get_Size() == 0 )
1934 {
1935 return( false );
1936 }
1937
1938 int nValues = m_Histogram.Get_Class_Count() > 0 ? (int)m_Histogram.Get_Class_Count() : m_Values.Get_N();
1939
1940 CSG_Matrix mv(nClasses, nValues); mv.Assign(FLT_MAX);
1941
1942 int **mc = (int **)SG_Malloc(nValues * sizeof(int *));
1943
1944 mc[0] = (int *)SG_Calloc((size_t)nClasses * nValues, sizeof(int));
1945
1946 for(int i=0; i<nValues; i++)
1947 {
1948 mc[i] = mc[0] + i * (size_t)nClasses;
1949 }
1950
1951 //-----------------------------------------------------
1952 for(int i=1; i<nValues; i++)
1953 {
1954 double v = 0., s1 = 0., s2 = 0., w = 0.;
1955
1956 for(int m=0, n=i+1; m<=i; m++, n--)
1957 {
1958 v = _Get_Value(n);
1959 s2 += v*v;
1960 s1 += v;
1961 w ++;
1962 v = s2 - (s1 * s1) / w;
1963
1964 if( n > 0 )
1965 {
1966 for(int j=1; j<nClasses; j++)
1967 {
1968 if( mv[i][j] >= (v + mv[n - 1][j - 1]) )
1969 {
1970 mc[i][j] = n;
1971 mv[i][j] = v + mv[n - 1][j - 1];
1972 }
1973 }
1974 }
1975 }
1976
1977 mc[i][0] = 0;
1978 mv[i][0] = v;
1979 }
1980
1981 //-----------------------------------------------------
1982 CSG_Array_Int Class(nClasses);
1983
1984 for(int i=0; i<nClasses; i++)
1985 {
1986 Class[i] = i;
1987 }
1988
1989 int j = Class[(size_t)nClasses - 1] = nValues - 1;
1990
1991 for(int i=nClasses-1; i>0; i--)
1992 {
1993 Class[(size_t)i - 1] = j = mc[j - 1][i];
1994 }
1995
1996 //-----------------------------------------------------
1997 m_Breaks.Create((size_t)nClasses + 1);
1998
1999 m_Breaks[0] = _Get_Value(0);
2000
2001 for(int i=1; i<nClasses; i++)
2002 {
2003 m_Breaks[i] = _Get_Value(Class[i - 1]);
2004 }
2005
2006 m_Breaks[nClasses] = _Get_Value(nValues - 1);
2007
2008 SG_Free(mc[0]); SG_Free(mc);
2009
2010 return( true );
2011}
2012
2013
2015// //
2016// //
2017// //
2019
2020//---------------------------------------------------------
2022{
2023 m_nFeatures = 0;
2024 m_Iteration = 0;
2025}
2026
2027//---------------------------------------------------------
2032
2033//---------------------------------------------------------
2035{
2036 m_Centroid.Destroy();
2037 m_Variance.Destroy();
2038 m_nMembers.Destroy();
2039 m_Clusters.Destroy();
2040 m_Features.Destroy();
2041 m_nFeatures = 0;
2042 m_Iteration = 0;
2043
2044 return( true );
2045}
2046
2047//---------------------------------------------------------
2049{
2050 Destroy();
2051
2052 if( nFeatures > 0 )
2053 {
2054 m_nFeatures = nFeatures;
2055
2056 m_Features.Create(m_nFeatures * sizeof(double), 0, TSG_Array_Growth::SG_ARRAY_GROWTH_3);
2057
2058 return( true );
2059 }
2060
2061 return( false );
2062}
2063
2064//---------------------------------------------------------
2066{
2067 return( m_nFeatures > 0 && m_Features.Inc_Array() );
2068}
2069
2070//---------------------------------------------------------
2071bool CSG_Cluster_Analysis::Set_Feature(sLong iElement, int iFeature, double Value)
2072{
2073 if( iElement >= 0 && iElement < Get_nElements() && iFeature >= 0 && iFeature < m_nFeatures )
2074 {
2075 ((double *)m_Features.Get_Entry(iElement))[iFeature] = Value;
2076
2077 return( true );
2078 }
2079
2080 return( false );
2081}
2082
2083//---------------------------------------------------------
2093//---------------------------------------------------------
2094bool CSG_Cluster_Analysis::Execute(int Method, int nClusters, int nMaxIterations, int Initialization)
2095{
2096 if( Get_nElements() < 2 || nClusters < 1 )
2097 {
2098 return( false );
2099 }
2100
2101 //-----------------------------------------------------
2102 m_nMembers.Create(nClusters);
2103 m_Variance.Create(nClusters);
2104 m_Centroid.Create(m_nFeatures, nClusters);
2105 m_Clusters.Create(Get_nElements());
2106
2107 //-----------------------------------------------------
2108 bool bResult = true;
2109
2110 if( nClusters == 1 )
2111 {
2112 for(int iElement=0; iElement<Get_nElements(); iElement++)
2113 {
2114 m_Clusters[iElement] = 0;
2115 }
2116
2117 _Update_Clusters();
2118 }
2119 else
2120 {
2121 for(int iElement=0; iElement<Get_nElements(); iElement++)
2122 {
2123 switch( Initialization )
2124 {
2125 default: // random
2126 if( (m_Clusters[iElement] = (int)CSG_Random::Get_Uniform(0, nClusters)) >= nClusters )
2127 {
2128 m_Clusters[iElement] = nClusters - 1;
2129 }
2130 break;
2131
2132 case 1: // periodic
2133 {
2134 m_Clusters[iElement] = iElement % nClusters;
2135 }
2136 break;
2137
2138 case 2: // keep as is, but check for valid cluster ids
2139 if( 0 > m_Clusters[iElement] || m_Clusters[iElement] >= nClusters )
2140 {
2141 m_Clusters[iElement] = iElement % nClusters;
2142 }
2143 break;
2144 }
2145 }
2146
2147 //-----------------------------------------------------
2148 m_Iteration = 0;
2149
2150 switch( Method )
2151 {
2152 default: bResult = _Minimum_Distance(true , nMaxIterations); break;
2153 case 1: bResult = _Hill_Climbing (true , nMaxIterations); break;
2154 case 2: bResult = _Minimum_Distance(true , nMaxIterations)
2155 && _Hill_Climbing (false, nMaxIterations); break;
2156 }
2157 }
2158
2159 //-----------------------------------------------------
2160 if( bResult )
2161 {
2162 m_Inertia = 0.;
2163
2164 for(int iElement=0; iElement<Get_nElements(); iElement++)
2165 {
2166 int iCluster = m_Clusters[iElement]; double d = 0., *Feature = (double *)m_Features.Get_Entry(iElement);
2167
2168 for(int iFeature=0; iFeature<m_nFeatures; iFeature++)
2169 {
2170 d += SG_Get_Square(Feature[iFeature] - m_Centroid[iCluster][iFeature]);
2171 }
2172
2173 m_Inertia += sqrt(d);
2174 }
2175
2176 m_Distortion = m_Inertia / Get_nElements();
2177
2178 //-------------------------------------------------
2179 m_WCSS = 0.;
2180
2181 for(int iCluster=0; iCluster<nClusters; iCluster++)
2182 {
2183 m_WCSS += m_Variance[iCluster];
2184
2185 m_Variance[iCluster] = m_nMembers[iCluster] <= 0 ? 0. : m_Variance[iCluster] / m_nMembers[iCluster];
2186 }
2187 }
2188
2189 return( bResult );
2190}
2191
2192//---------------------------------------------------------
2193bool CSG_Cluster_Analysis::_Update_Clusters(void)
2194{
2195 m_nMembers = 0; m_Centroid = 0.; m_Variance = 0.;
2196
2197 for(int iElement=0; iElement<Get_nElements(); iElement++)
2198 {
2199 int iCluster = m_Clusters[iElement]; m_nMembers[iCluster]++;
2200
2201 double *Feature = (double *)m_Features.Get_Entry(iElement);
2202
2203 for(int iFeature=0; iFeature<m_nFeatures; iFeature++)
2204 {
2205 double d = Feature[iFeature];
2206
2207 m_Centroid[iCluster][iFeature] += d; m_Variance[iCluster] += d*d;
2208 }
2209 }
2210
2211 //-------------------------------------------------
2212 for(int iCluster=0; iCluster<Get_nClusters(); iCluster++)
2213 {
2214 double v = 0., d = m_nMembers[iCluster] < 1 ? 0. : 1. / (double)m_nMembers[iCluster];
2215
2216 for(int iFeature=0; iFeature<m_nFeatures; iFeature++)
2217 {
2218 m_Centroid[iCluster][iFeature] *= d; v += SG_Get_Square(m_Centroid[iCluster][iFeature]);
2219 }
2220
2221 m_Variance[iCluster] -= v * m_nMembers[iCluster];
2222 }
2223
2224 return( true );
2225}
2226
2227//---------------------------------------------------------
2228bool CSG_Cluster_Analysis::_Minimum_Distance(bool bInitialize, int nMaxIterations)
2229{
2230 double SP_Last = -1.;
2231
2232 for(m_Iteration=1; SG_UI_Process_Get_Okay(); m_Iteration++)
2233 {
2234 _Update_Clusters();
2235
2236 int nShifts = 0; m_SP = 0.;
2237
2238 for(int iElement=0; iElement<Get_nElements(); iElement++)
2239 {
2240 double *Feature = (double *)m_Features.Get_Entry(iElement);
2241
2242 double minVariance = -1.;
2243 int minCluster = -1;
2244
2245 for(int iCluster=0; iCluster<Get_nClusters(); iCluster++)
2246 {
2247 double Variance = 0.;
2248
2249 for(int iFeature=0; iFeature<m_nFeatures; iFeature++)
2250 {
2251 Variance += SG_Get_Square(m_Centroid[iCluster][iFeature] - Feature[iFeature]);
2252 }
2253
2254 if( minVariance < 0. || Variance < minVariance )
2255 {
2256 minVariance = Variance;
2257 minCluster = iCluster;
2258 }
2259 }
2260
2261 if( m_Clusters[iElement] != minCluster )
2262 {
2263 m_Clusters[iElement] = minCluster;
2264
2265 nShifts++;
2266 }
2267
2268 m_SP += minVariance;
2269 m_Variance[minCluster] += minVariance;
2270 }
2271
2272 //-------------------------------------------------
2273 m_SP /= Get_nElements();
2274
2275 SG_UI_StatusBar_Set_Text(CSG_String::Format("[k = %d] %s: %d >> %s %f", Get_nClusters(),
2276 _TL("pass" ), m_Iteration,
2277 _TL("change"), m_Iteration < 2 ? m_SP : SP_Last - m_SP
2278 ));
2279
2280 SP_Last = m_SP;
2281
2282 if( nShifts == 0 || (nMaxIterations > 0 && nMaxIterations <= m_Iteration) )
2283 {
2284 if( nShifts > 0 )
2285 {
2286 _Update_Clusters();
2287 }
2288
2289 return( true );
2290 }
2291 }
2292
2293 return( true );
2294}
2295
2296//---------------------------------------------------------
2297bool CSG_Cluster_Analysis::_Hill_Climbing(bool bInitialize, int nMaxIterations)
2298{
2299 _Update_Clusters();
2300
2301 //-----------------------------------------------------
2302 double SP_Last = -1.;
2303
2304 for(m_Iteration=1; SG_UI_Process_Get_Okay(false); m_Iteration++)
2305 {
2306 int nShifts = 0;
2307
2308 for(int iElement=0; iElement<Get_nElements(); iElement++)
2309 {
2310 int iCluster = m_Clusters[iElement];
2311
2312 if( m_nMembers[iCluster] > 1 )
2313 {
2314 double Variance = 0., *Feature = (double *)m_Features.Get_Entry(iElement);
2315
2316 for(int iFeature=0; iFeature<m_nFeatures; iFeature++)
2317 {
2318 Variance += SG_Get_Square(m_Centroid[iCluster][iFeature] - Feature[iFeature]);
2319 }
2320
2321 double V1 = Variance * m_nMembers[iCluster] / (m_nMembers[iCluster] - 1.);
2322
2323 //-----------------------------------------
2324 int kCluster = 0; double VMin = -1.;
2325
2326 for(int jCluster=0; jCluster<Get_nClusters(); jCluster++)
2327 {
2328 if( jCluster != iCluster )
2329 {
2330 Variance = 0.;
2331
2332 for(int iFeature=0; iFeature<m_nFeatures; iFeature++)
2333 {
2334 Variance += SG_Get_Square(m_Centroid[jCluster][iFeature] - Feature[iFeature]);
2335 }
2336
2337 double V2 = Variance * m_nMembers[jCluster] / (m_nMembers[jCluster] + 1.);
2338
2339 if( VMin < 0. || V2 < VMin )
2340 {
2341 VMin = V2;
2342 kCluster = jCluster;
2343 }
2344 }
2345 }
2346
2347 //-----------------------------------------
2348 if( VMin >= 0 && VMin < V1 )
2349 {
2350 m_Variance[iCluster] -= V1;
2351 m_Variance[kCluster] += VMin;
2352 V1 = 1. / (m_nMembers[iCluster] - 1.);
2353 double V2 = 1. / (m_nMembers[kCluster] + 1.);
2354
2355 for(int iFeature=0; iFeature<m_nFeatures; iFeature++)
2356 {
2357 double d = Feature[iFeature];
2358
2359 m_Centroid[iCluster][iFeature] = (m_nMembers[iCluster] * m_Centroid[iCluster][iFeature] - d) * V1;
2360 m_Centroid[kCluster][iFeature] = (m_nMembers[kCluster] * m_Centroid[kCluster][iFeature] + d) * V2;
2361 }
2362
2363 m_Clusters[iElement] = kCluster;
2364
2365 m_nMembers[iCluster]--;
2366 m_nMembers[kCluster]++;
2367
2368 nShifts++;
2369 }
2370 }
2371 }
2372
2373 //-------------------------------------------------
2374 m_SP = 0.;
2375
2376 for(int iCluster=0; iCluster<Get_nClusters(); iCluster++)
2377 {
2378 m_SP += m_Variance[iCluster];
2379 }
2380
2381 m_SP /= Get_nElements();
2382
2383 SG_UI_StatusBar_Set_Text(CSG_String::Format("[k = %d] %s: %d >> %s %f", Get_nClusters(),
2384 _TL("pass" ), m_Iteration,
2385 _TL("change"), m_Iteration <= 1 ? m_SP : SP_Last - m_SP
2386 ));
2387
2388 SP_Last = m_SP;
2389
2390 if( nShifts == 0 || (nMaxIterations > 0 && nMaxIterations <= m_Iteration) )
2391 {
2392 return( true );
2393 }
2394 }
2395
2396 return( true );
2397}
2398
2399
2401// //
2402// //
2403// //
2405
2406//---------------------------------------------------------
2408{
2409 m_nFeatures = 0; m_pClasses = NULL; m_nClasses = 0;
2410
2411 m_Threshold_Distance = 0.;
2412 m_Threshold_Angle = 0.;
2413 m_Threshold_Probability = 0.;
2414 m_Probability_Relative = false;
2415
2416 for(int i=0; i<SG_CLASSIFY_SUPERVISED_WTA; i++)
2417 {
2419 // || i == SG_CLASSIFY_SUPERVISED_Mahalonobis
2422 }
2423}
2424
2425//---------------------------------------------------------
2430
2431//---------------------------------------------------------
2433{
2434 Destroy();
2435
2436 if( nFeatures > 0 )
2437 {
2438 m_nFeatures = nFeatures;
2439 }
2440}
2441
2442//---------------------------------------------------------
2444{
2445 if( m_nClasses > 0 )
2446 {
2447 for(int i=0; i<m_nClasses; i++)
2448 {
2449 delete(m_pClasses[i]);
2450 }
2451
2452 SG_FREE_SAFE(m_pClasses);
2453 }
2454
2455 m_nFeatures = 0;
2456
2457 m_Info.Clear();
2458}
2459
2460
2462// //
2464
2465//---------------------------------------------------------
2466void CSG_Classifier_Supervised::Set_Threshold_Distance (double Value) { m_Threshold_Distance = Value; }
2467double CSG_Classifier_Supervised::Get_Threshold_Distance (void) { return( m_Threshold_Distance ); }
2468
2469//---------------------------------------------------------
2470void CSG_Classifier_Supervised::Set_Threshold_Angle (double Value) { m_Threshold_Angle = Value; }
2471double CSG_Classifier_Supervised::Get_Threshold_Angle (void) { return( m_Threshold_Angle ); }
2472
2473//---------------------------------------------------------
2474void CSG_Classifier_Supervised::Set_Threshold_Probability(double Value) { m_Threshold_Probability = Value; }
2475double CSG_Classifier_Supervised::Get_Threshold_Probability(void) { return( m_Threshold_Probability ); }
2476
2477//---------------------------------------------------------
2478void CSG_Classifier_Supervised::Set_Probability_Relative (bool Value) { m_Probability_Relative = Value; }
2479bool CSG_Classifier_Supervised::Get_Probability_Relative (void) { return( m_Probability_Relative ); }
2480
2481//---------------------------------------------------------
2482void CSG_Classifier_Supervised::Set_WTA(int Method, bool bOn)
2483{
2484 if( Method >= 0 && Method < SG_CLASSIFY_SUPERVISED_WTA )
2485 {
2486 m_bWTA[Method] = bOn;
2487 }
2488}
2489
2491{
2492 return( Method >= 0 && Method < SG_CLASSIFY_SUPERVISED_WTA ? m_bWTA[Method] : false );
2493}
2494
2495
2497// //
2499
2500//---------------------------------------------------------
2501#include "saga_api.h"
2502
2503//---------------------------------------------------------
2505{
2506 int nFeatures = m_nFeatures; Destroy(); m_nFeatures = nFeatures;
2507
2508 //-----------------------------------------------------
2509 CSG_MetaData Data;
2510
2511 if( !Data.Load(File) || !Data.Cmp_Name("supervised_classifier") || SG_Compare_Version(Data.Get_Property("saga-version"), "2.1.4") < 0 )
2512 {
2513 return( false );
2514 }
2515
2516 if( !Data("classes") || !Data("features") || !Data["features"]("count") || Data["features"]["count"].Get_Content().asInt() != m_nFeatures || m_nFeatures == 0 )
2517 {
2518 return( false );
2519 }
2520
2521 if( Data["features"]("info") )
2522 {
2523 m_Info = Data["features"]["info"].Get_Content();
2524 }
2525
2526 //-----------------------------------------------------
2527 CSG_MetaData &Classes = *Data.Get_Child("CLASSES");
2528
2529 for(int i=0; i<Classes.Get_Children_Count(); i++)
2530 {
2531 if( Classes[i].Cmp_Name("class") && Classes[i].Get_Child("id") )
2532 {
2533 bool bAdd = true;
2534
2535 CClass *pClass = new CClass(Classes[i]["id"].Get_Content());
2536
2537 if( !pClass->m_Cov .from_String(Classes[i]["cov" ].Get_Content()) || pClass->m_Cov .Get_NX() != m_nFeatures || !pClass->m_Cov.is_Square() ) { bAdd = false; }
2538 if( !pClass->m_Mean.from_String(Classes[i]["mean"].Get_Content()) || pClass->m_Mean.Get_N () != m_nFeatures ) { bAdd = false; }
2539 if( !pClass->m_Min .from_String(Classes[i]["min" ].Get_Content()) || pClass->m_Min .Get_N () != m_nFeatures ) { bAdd = false; }
2540 if( !pClass->m_Max .from_String(Classes[i]["max" ].Get_Content()) || pClass->m_Max .Get_N () != m_nFeatures ) { bAdd = false; }
2541
2542 //---------------------------------------------
2543 if( !bAdd )
2544 {
2545 delete(pClass);
2546 }
2547 else
2548 {
2549 m_pClasses = (CClass **)SG_Realloc(m_pClasses, ((size_t)m_nClasses + 1) * sizeof(CClass *));
2550 m_pClasses[m_nClasses++] = pClass;
2551
2552 pClass->m_Cov_Det = pClass->m_Cov.Get_Determinant();
2553 pClass->m_Cov_Inv = pClass->m_Cov.Get_Inverse();
2554
2555 pClass->m_Mean_Spectral = CSG_Simple_Statistics(pClass->m_Mean).Get_Mean();
2556 }
2557 }
2558 }
2559
2560 return( m_nClasses > 0 );
2561}
2562
2563//---------------------------------------------------------
2564bool CSG_Classifier_Supervised::Save(const CSG_String &File, const SG_Char *Feature_Info)
2565{
2566 if( m_nFeatures < 1 || m_nClasses < 1 || File.is_Empty() )
2567 {
2568 return( false );
2569 }
2570
2571 CSG_MetaData Data;
2572
2573 Data.Set_Name ("supervised_classifier");
2574 Data.Add_Property("saga-version", SAGA_VERSION);
2575
2576 CSG_MetaData &Features = *Data.Add_Child("features");
2577
2578 Features.Add_Child("count", m_nFeatures);
2579
2580 if( Feature_Info && *Feature_Info )
2581 {
2582 Features.Add_Child("info", Feature_Info);
2583 }
2584
2585 CSG_MetaData &Classes = *Data.Add_Child("classes");
2586
2587 Classes.Add_Property("count", m_nClasses);
2588
2589 for(int i=0; i<m_nClasses; i++)
2590 {
2591 CSG_MetaData &Class = *Classes.Add_Child("class");
2592
2593 CClass *pClass = m_pClasses[i];
2594
2595 Class.Add_Child("id" , pClass->m_ID );
2596 Class.Add_Child("mean", pClass->m_Mean.to_String());
2597 Class.Add_Child("min" , pClass->m_Min .to_String());
2598 Class.Add_Child("max" , pClass->m_Max .to_String());
2599 Class.Add_Child("cov" , pClass->m_Cov .to_String());
2600 }
2601
2602 return( Data.Save(File) );
2603}
2604
2605
2607// //
2609
2610//---------------------------------------------------------
2612{
2613 CSG_String s;
2614
2615 if( m_nFeatures > 0 && m_nClasses > 0 )
2616 {
2617 s += "\n";
2618
2619 for(int iClass=0; iClass<m_nClasses; iClass++)
2620 {
2621 CClass *pClass = m_pClasses[iClass];
2622
2623 s += "\n____\n" + pClass->m_ID + "\nFeature\tMean\tMin\tMax\tStdDev";
2624
2625 for(int i=0; i<m_nFeatures; i++)
2626 {
2627 s += CSG_String::Format("\n%3d.", i + 1);
2628 s += "\t" + SG_Get_String(pClass->m_Mean[i]);
2629 s += "\t" + SG_Get_String(pClass->m_Min [i]);
2630 s += "\t" + SG_Get_String(pClass->m_Max [i]);
2631 s += "\t" + SG_Get_String(sqrt(pClass->m_Cov[i][i]));
2632 }
2633
2634 s += "\n";
2635 }
2636 }
2637
2638 return( s );
2639}
2640
2641
2643// //
2645
2646//---------------------------------------------------------
2647bool CSG_Classifier_Supervised::Add_Class(const CSG_String &Class_ID, const CSG_Vector &Mean, const CSG_Vector &Min, const CSG_Vector &Max, const CSG_Matrix &Cov)
2648{
2649 if( m_nFeatures < 1 || Mean.Get_N() != m_nFeatures || Min.Get_N() != m_nFeatures || Max.Get_N() != m_nFeatures || Cov.Get_NCols() != m_nFeatures || Cov.Get_NRows() != m_nFeatures )
2650 {
2651 return( false );
2652 }
2653
2654 CClass *pClass, **pClasses = (CClass **)SG_Realloc(m_pClasses, ((size_t)m_nClasses + 1) * sizeof(CClass *));
2655
2656 if( pClasses )
2657 {
2658 m_pClasses = pClasses;
2659
2660 m_pClasses[m_nClasses++] = pClass = new CClass(Class_ID);
2661
2662 pClass->m_ID = Class_ID;
2663 pClass->m_Mean = Mean;
2664 pClass->m_Min = Min;
2665 pClass->m_Max = Max;
2666 pClass->m_Cov = Cov;
2667 pClass->m_Cov_Inv = Cov.Get_Inverse();
2668 pClass->m_Cov_Det = Cov.Get_Determinant();
2669
2670 pClass->m_Mean_Spectral = CSG_Simple_Statistics(Mean).Get_Mean();
2671
2672 return( true );
2673 }
2674
2675 return( false );
2676}
2677
2678
2680// //
2682
2683//---------------------------------------------------------
2685{
2686 for(int i=0; i<m_nClasses; i++)
2687 {
2688 m_pClasses[i]->m_Samples.Destroy();
2689 }
2690
2691 return( true );
2692}
2693
2694//---------------------------------------------------------
2696{
2697 if( m_nFeatures > 0 && m_nFeatures == Features.Get_N() )
2698 {
2699 int iClass = Get_Class(Class_ID);
2700
2701 if( iClass < 0 )
2702 {
2703 CClass **pClasses = (CClass **)SG_Realloc(m_pClasses, ((size_t)m_nClasses + 1) * sizeof(CClass *));
2704
2705 if( pClasses )
2706 {
2707 m_pClasses = pClasses;
2708
2709 m_pClasses[iClass = m_nClasses++] = new CClass(Class_ID);
2710 }
2711 }
2712
2713 if( iClass >= 0 )
2714 {
2715 return( m_pClasses[iClass]->m_Samples.Add_Row(Features) );
2716 }
2717 }
2718
2719 return( false );
2720}
2721
2722//---------------------------------------------------------
2723bool CSG_Classifier_Supervised::Train(bool bClear_Samples)
2724{
2725 if( m_nFeatures < 1 || m_nClasses < 1 )
2726 {
2727 return( false );
2728 }
2729
2730 for(int iClass=0; iClass<m_nClasses; iClass++)
2731 {
2732 if( !m_pClasses[iClass]->Train() )
2733 {
2734 return( false );
2735 }
2736 }
2737
2738 if( bClear_Samples )
2739 {
2741 }
2742
2743 return( true );
2744}
2745
2746
2748// //
2750
2751//---------------------------------------------------------
2752bool CSG_Classifier_Supervised::CClass::Train(void)
2753{
2754 if( m_Samples.Get_NCols() < 1 || m_Samples.Get_NRows() < 1 )
2755 {
2756 return( false );
2757 }
2758
2759 //-----------------------------------------------------
2760 m_Mean.Create(m_Samples.Get_NCols());
2761 m_Min .Create(m_Samples.Get_NCols());
2762 m_Max .Create(m_Samples.Get_NCols());
2763
2764 for(int iFeature=0; iFeature<m_Samples.Get_NCols(); iFeature++)
2765 {
2766 CSG_Simple_Statistics s;
2767
2768 for(int iSample=0; iSample<m_Samples.Get_NRows(); iSample++)
2769 {
2770 s += m_Samples[iSample][iFeature];
2771 }
2772
2773 m_Mean[iFeature] = s.Get_Mean ();
2774 m_Min [iFeature] = s.Get_Minimum();
2775 m_Max [iFeature] = s.Get_Maximum();
2776 }
2777
2778 //-----------------------------------------------------
2779 m_Cov.Create(m_Samples.Get_NCols(), m_Samples.Get_NCols());
2780
2781 for(int iFeature=0; iFeature<m_Samples.Get_NCols(); iFeature++)
2782 {
2783 for(int jFeature=iFeature; jFeature<m_Samples.Get_NCols(); jFeature++)
2784 {
2785 double cov = 0.;
2786
2787 for(int iSample=0; iSample<m_Samples.Get_NRows(); iSample++)
2788 {
2789 cov += (m_Samples[iSample][iFeature] - m_Mean[iFeature]) * (m_Samples[iSample][jFeature] - m_Mean[jFeature]);
2790 }
2791
2792 if( m_Samples.Get_NRows() > 1 )
2793 {
2794 cov /= m_Samples.Get_NRows() - 1;
2795 }
2796
2797 m_Cov[iFeature][jFeature] = m_Cov[jFeature][iFeature] = cov;
2798 }
2799 }
2800
2801 m_Cov_Inv = m_Cov.Get_Inverse ();
2802 m_Cov_Det = m_Cov.Get_Determinant();
2803
2804 m_Mean_Spectral = CSG_Simple_Statistics(m_Mean).Get_Mean();
2805
2806 //-----------------------------------------------------
2807 return( true );
2808}
2809
2810
2812// //
2814
2815//---------------------------------------------------------
2817{
2818 if( m_nFeatures > 0 )
2819 {
2820 for(int iClass=0; iClass<Get_Class_Count(); iClass++)
2821 {
2822 if( !Get_Class_ID(iClass).Cmp(Class_ID) )
2823 {
2824 return( iClass );
2825 }
2826 }
2827 }
2828
2829 return( -1 );
2830}
2831
2832//---------------------------------------------------------
2833bool CSG_Classifier_Supervised::Get_Class(const CSG_Vector &Features, int &Class, double &Quality, int Method)
2834{
2835 Class = -1; Quality = 0.;
2836
2837 if( Get_Feature_Count() == Features.Get_N() )
2838 {
2839 switch( Method )
2840 {
2841 case SG_CLASSIFY_SUPERVISED_BinaryEncoding : _Get_Binary_Encoding (Features, Class, Quality); break;
2842 case SG_CLASSIFY_SUPERVISED_ParallelEpiped : _Get_Parallel_Epiped (Features, Class, Quality); break;
2843 case SG_CLASSIFY_SUPERVISED_MinimumDistance : _Get_Minimum_Distance (Features, Class, Quality); break;
2844 case SG_CLASSIFY_SUPERVISED_Mahalonobis : _Get_Mahalanobis_Distance (Features, Class, Quality); break;
2845 case SG_CLASSIFY_SUPERVISED_MaximumLikelihood: _Get_Maximum_Likelihood (Features, Class, Quality); break;
2846 case SG_CLASSIFY_SUPERVISED_SAM : _Get_Spectral_Angle_Mapping(Features, Class, Quality); break;
2847 case SG_CLASSIFY_SUPERVISED_SID : _Get_Spectral_Divergence (Features, Class, Quality); break;
2848 case SG_CLASSIFY_SUPERVISED_WTA : _Get_Winner_Takes_All (Features, Class, Quality); break;
2849 }
2850
2851 return( Class >= 0 );
2852 }
2853
2854 return( false );
2855}
2856
2857
2859// //
2861
2862//---------------------------------------------------------
2864{
2865 switch( Method )
2866 {
2867 case SG_CLASSIFY_SUPERVISED_BinaryEncoding : return( _TL("Binary Encoding") );
2868 case SG_CLASSIFY_SUPERVISED_ParallelEpiped : return( _TL("Parallelepiped") );
2869 case SG_CLASSIFY_SUPERVISED_MinimumDistance : return( _TL("Minimum Distance") );
2870 case SG_CLASSIFY_SUPERVISED_Mahalonobis : return( _TL("Mahalanobis Distance") );
2871 case SG_CLASSIFY_SUPERVISED_MaximumLikelihood: return( _TL("Maximum Likelihood") );
2872 case SG_CLASSIFY_SUPERVISED_SAM : return( _TL("Spectral Angle Mapping") );
2873 case SG_CLASSIFY_SUPERVISED_SID : return( _TL("Spectral Information Divergence") );
2874 case SG_CLASSIFY_SUPERVISED_SVM : return( _TL("Support Vector Machine") );
2875 case SG_CLASSIFY_SUPERVISED_WTA : return( _TL("Winner Takes All") );
2876 }
2877
2878 return( SG_T("") );
2879}
2880
2881//---------------------------------------------------------
2883{
2884 switch( Method )
2885 {
2886 case SG_CLASSIFY_SUPERVISED_BinaryEncoding : return( _TL("Difference") );
2887 case SG_CLASSIFY_SUPERVISED_ParallelEpiped : return( _TL("Memberships") );
2888 case SG_CLASSIFY_SUPERVISED_MinimumDistance : return( _TL("Distance") );
2889 case SG_CLASSIFY_SUPERVISED_Mahalonobis : return( _TL("Distance") );
2890 case SG_CLASSIFY_SUPERVISED_MaximumLikelihood: return( _TL("Proximity") );
2891 case SG_CLASSIFY_SUPERVISED_SAM : return( _TL("Angle") );
2892 case SG_CLASSIFY_SUPERVISED_SID : return( _TL("Divergence") );
2893 case SG_CLASSIFY_SUPERVISED_SVM : return( _TL("") );
2894 case SG_CLASSIFY_SUPERVISED_WTA : return( _TL("Votes") );
2895 }
2896
2897 return( SG_T("") );
2898}
2899
2900
2902// //
2904
2905//---------------------------------------------------------
2906// Mazer, A. S., Martin, M., Lee, M., and Solomon, J. E. (1988):
2907// Image Processing Software for Imaging Spectrometry Analysis.
2908// Remote Sensing of Environment, v. 24, no. 1, p. 201-210.
2909//
2910void CSG_Classifier_Supervised::_Get_Binary_Encoding(const CSG_Vector &Features, int &Class, double &Quality)
2911{
2912 for(int iClass=0; iClass<Get_Class_Count(); iClass++)
2913 {
2914 CClass *pClass = m_pClasses[iClass];
2915
2916 double Mean_Spectral = CSG_Simple_Statistics(Features).Get_Mean();
2917
2918 int d = 0;
2919
2920 for(int iFeature=0; iFeature<Get_Feature_Count(); iFeature++)
2921 {
2922 d += (Features(iFeature) < Mean_Spectral) == (pClass->m_Mean[iFeature] < pClass->m_Mean_Spectral) ? 0 : 1;
2923
2924 if( iFeature == 0 ) // spectral slopes
2925 {
2926 d += (Features[iFeature ] < Features[iFeature + 1]) == (pClass->m_Mean[iFeature ] < pClass->m_Mean[iFeature + 1]) ? 0 : 1;
2927 }
2928 else if( iFeature == Get_Feature_Count() - 1 )
2929 {
2930 d += (Features[iFeature - 1] < Features[iFeature ]) == (pClass->m_Mean[iFeature - 1] < pClass->m_Mean[iFeature ]) ? 0 : 1;
2931 }
2932 else
2933 {
2934 d += (Features[iFeature - 1] < Features[iFeature + 1]) == (pClass->m_Mean[iFeature - 1] < pClass->m_Mean[iFeature + 1]) ? 0 : 1;
2935 }
2936 }
2937
2938 if( Class < 0 || Quality > d ) // find the minimum 'Hamming' distance
2939 {
2940 Class = iClass; Quality = d;
2941 }
2942 }
2943}
2944
2945//---------------------------------------------------------
2946void CSG_Classifier_Supervised::_Get_Parallel_Epiped(const CSG_Vector &Features, int &Class, double &Quality)
2947{
2948 for(int iClass=0; iClass<Get_Class_Count(); iClass++)
2949 {
2950 CClass *pClass = m_pClasses[iClass];
2951
2952 bool bMember = true;
2953
2954 for(int iFeature=0; bMember && iFeature<Get_Feature_Count(); iFeature++)
2955 {
2956 bMember = pClass->m_Min[iFeature] <= Features[iFeature] && Features[iFeature] <= pClass->m_Max[iFeature];
2957 }
2958
2959 if( bMember )
2960 {
2961 Class = iClass; Quality++;
2962 }
2963 }
2964}
2965
2966//---------------------------------------------------------
2967void CSG_Classifier_Supervised::_Get_Minimum_Distance(const CSG_Vector &Features, int &Class, double &Quality)
2968{
2969 for(int iClass=0; iClass<Get_Class_Count(); iClass++)
2970 {
2971 CClass *pClass = m_pClasses[iClass];
2972
2973 double Distance = (Features - pClass->m_Mean).Get_Length();
2974
2975 if( Class < 0 || Quality > Distance )
2976 {
2977 Class = iClass; Quality = Distance;
2978 }
2979 }
2980
2981 if( m_Threshold_Distance > 0. && Quality > m_Threshold_Distance )
2982 {
2983 Class = -1;
2984 }
2985}
2986
2987//---------------------------------------------------------
2988void CSG_Classifier_Supervised::_Get_Mahalanobis_Distance(const CSG_Vector &Features, int &Class, double &Quality)
2989{
2990 for(int iClass=0; iClass<Get_Class_Count(); iClass++)
2991 {
2992 CClass *pClass = m_pClasses[iClass];
2993
2994 CSG_Vector D = Features - pClass->m_Mean;
2995
2996 double Distance = D * (pClass->m_Cov_Inv * D);
2997
2998 if( Class < 0 || Quality > Distance )
2999 {
3000 Class = iClass; Quality = Distance;
3001 }
3002 }
3003
3004 if( m_Threshold_Distance > 0. && Quality > m_Threshold_Distance )
3005 {
3006 Class = -1;
3007 }
3008}
3009
3010//---------------------------------------------------------
3011void CSG_Classifier_Supervised::_Get_Maximum_Likelihood(const CSG_Vector &Features, int &Class, double &Quality)
3012{
3013 double dSum = 0.;
3014
3015 for(int iClass=0; iClass<Get_Class_Count(); iClass++)
3016 {
3017 CClass *pClass = m_pClasses[iClass];
3018
3019 CSG_Vector D = Features - pClass->m_Mean;
3020
3021 double Distance = D * (pClass->m_Cov_Inv * D);
3022
3023 double Probability = pow(2. * M_PI, -0.5 * m_nFeatures) * pow(pClass->m_Cov_Det, -0.5) * exp(-0.5 * Distance);
3024 // double Probability = -log(pClass->m_Cov_Det) - Distance;
3025
3026 dSum += Probability;
3027
3028 if( Class < 0 || Quality < Probability )
3029 {
3030 Class = iClass; Quality = Probability;
3031 }
3032 }
3033
3034 if( Class >= 0 )
3035 {
3036 if( m_Probability_Relative )
3037 {
3038 Quality = 100. * Quality / dSum;
3039 }
3040
3041 if( m_Threshold_Probability > 0. && Quality < m_Threshold_Probability )
3042 {
3043 Class = -1;
3044 }
3045 }
3046}
3047
3048//---------------------------------------------------------
3049void CSG_Classifier_Supervised::_Get_Spectral_Angle_Mapping(const CSG_Vector &Features, int &Class, double &Quality)
3050{
3051 for(int iClass=0; iClass<Get_Class_Count(); iClass++)
3052 {
3053 CClass *pClass = m_pClasses[iClass];
3054
3055 double Angle = Features.Get_Angle(pClass->m_Mean);
3056
3057 if( Class < 0 || Quality > Angle )
3058 {
3059 Class = iClass; Quality = Angle;
3060 }
3061 }
3062
3063 Quality *= M_RAD_TO_DEG;
3064
3065 if( m_Threshold_Angle > 0. && Quality > m_Threshold_Angle )
3066 {
3067 Class = -1;
3068 }
3069}
3070
3071//---------------------------------------------------------
3072void CSG_Classifier_Supervised::_Get_Spectral_Divergence(const CSG_Vector &Features, int &Class, double &Quality)
3073{
3074}
3075
3076//---------------------------------------------------------
3077void CSG_Classifier_Supervised::_Get_Winner_Takes_All(const CSG_Vector &Features, int &Class, double &Quality)
3078{
3079 int *Votes = (int *)SG_Calloc(Get_Class_Count(), sizeof(int));
3080
3081 for(int iMethod=0; iMethod<SG_CLASSIFY_SUPERVISED_WTA; iMethod++)
3082 {
3083 int iClass; double iQuality;
3084
3085 if( m_bWTA[iMethod] && Get_Class(Features, iClass, iQuality, iMethod) && ++Votes[iClass] > Quality )
3086 {
3087 Class = iClass; Quality = Votes[iClass];
3088 }
3089 }
3090
3091 SG_Free(Votes);
3092}
3093
3094
3096// //
3097// //
3098// //
3100
3101//---------------------------------------------------------
3102// source: http://psydok.sulb.uni-saarland.de/volltexte/2004/268/html/
3103
3104//---------------------------------------------------------
3106{ // Hill's approx. to cumulative t-dist, Commun.A.C.M. 13,617-619.
3107 // See: J.H.Maindonald, Computational Statistics, p.295.
3108 // Calculates p given t and tail type.
3109
3110 if( !T || !df || df < 1. )
3111 {
3112 return( -1. );
3113 }
3114
3115 return( _Change_Tail_Type(Get_T_P(T, df), TESTDIST_TYPE_TwoTail, Type, T < 0.) );
3116}
3117
3118//---------------------------------------------------------
3120{ // Keith Dear & Robert Brennan.
3121 // Returns an accurate t to tol sig. fig.'s given p & df.
3122
3123 if( p <= 0. || p >= 1. || df < 1 )
3124 {
3125 return( -1. );
3126 }
3127
3128 bool bNegative = (Type == TESTDIST_TYPE_Left && p < 0.5) || (Type == TESTDIST_TYPE_Right && p > 0.5);
3129 double t = 0, p0, p1, diff = 1.;
3130
3131 p0 = p1 = _Change_Tail_Type(p, Type, TESTDIST_TYPE_TwoTail, bNegative);
3132
3133 while( fabs(diff) > 0.0001 )
3134 {
3135 t = Get_T_Inv(p1, df); // initial rough value
3136 diff = Get_T_P(t, df) - p0; // compare result with forward fn
3137 p1 = p1 - diff; // small adjustment to p1
3138 }
3139
3140 return( bNegative ? -t : t );
3141}
3142
3143//---------------------------------------------------------
3144double CSG_Test_Distribution::_Change_Tail_Type(double p, TSG_Test_Distribution_Type from, TSG_Test_Distribution_Type to, bool bNegative)
3145{
3146 if( from != to )
3147 {
3148 switch( from ) // convert any tail type to 'left'
3149 {
3150 case TESTDIST_TYPE_Left : break;
3151 case TESTDIST_TYPE_Right : p = 1. - p; break;
3152 case TESTDIST_TYPE_Middle : p = p / 2. + 0.5; if( bNegative ) p = 1. - p; break;
3153 case TESTDIST_TYPE_TwoTail: p = 1. - p / 2. ; if( bNegative ) p = 1. - p; break;
3154 // case TESTDIST_TYPE_Half : p = p + 0.5 ; if( bNegative ) p = 1. - p; break;
3155 }
3156
3157 switch( to ) // convert p from tail type 'left' to any other
3158 {
3159 case TESTDIST_TYPE_Left : break;
3160 case TESTDIST_TYPE_Right : p = 1. - p; break;
3161 case TESTDIST_TYPE_Middle : if( bNegative ) p = 1. - p; p = 2. * (1. - p); break;
3162 case TESTDIST_TYPE_TwoTail: if( bNegative ) p = 1. - p; p = 2. * p - 1. ; break;
3163 // case TESTDIST_TYPE_Half : if( bNegative ) p = 1. - p; p = p - 0.5 ; break;
3164 }
3165 }
3166
3167 return( p );
3168}
3169
3170//---------------------------------------------------------
3172{ // Returns the two-tailed standard normal probability of z
3173 const double a1 = 0.0000053830, a2 = 0.0000488906, a3 = 0.0000380036;
3174 const double a4 = 0.0032776263, a5 = 0.0211410061, a6 = 0.0498673470;
3175
3176 z = fabs(z);
3177
3178 double p = (((((a1 * z + a2) * z + a3) * z + a4) * z + a5) * z + a6) * z + 1.;
3179
3180 return( pow(p, -16) );
3181}
3182
3183//---------------------------------------------------------
3185{ // Returns z given a half-middle tail type p.
3186 const double a0 = 2.5066282, a1 = -18.6150006, a2 = 41.3911977, a3 = -25.4410605;
3187 const double b1 = -8.4735109, b2 = 23.0833674, b3 = -21.0622410, b4 = 3.1308291;
3188 const double c0 = -2.7871893, c1 = -2.2979648, c2 = 4.8501413, c3 = 2.3212128;
3189 const double d1 = 3.5438892, d2 = 1.6370678;
3190
3191 if( p > 0.42 )
3192 {
3193 double r = sqrt(-log(0.5 - p));
3194
3195 return( (((c3 * r + c2) * r + c1) * r + c0) / ((d2 * r + d1) * r + 1.) );
3196 }
3197 else
3198 {
3199 double r = p * p;
3200
3201 return( p * (((a3 * r + a2) * r + a1) * r + a0) / ((((b4 * r + b3) * r + b2) * r + b1) * r + 1.) );
3202 }
3203}
3204
3205//---------------------------------------------------------
3206double CSG_Test_Distribution::Get_T_P(double T, int df)
3207{ // Returns two-tail probability level given t and df.
3208 return( df == 1 ? 1. - 2. * atan(fabs(T)) / M_PI
3209 : df == 2 ? 1. - fabs(T) / sqrt(T*T + 2.)
3210 : df == 3 ? 1. - 2. * (atan(fabs(T) / sqrt(3.)) + fabs(T) * sqrt(3.) / (T*T + 3.)) / M_PI
3211 : df == 4 ? 1. - fabs(T) * (1. + 2. / (T*T + 4.)) / sqrt(T*T + 4.)
3212 : Get_Norm_P(Get_T_Z(fabs(T), df))
3213 );
3214}
3215
3216//---------------------------------------------------------
3217double CSG_Test_Distribution::Get_T_Z(double T, int df)
3218{ // Converts a t value to an approximate z value w.r.t the given df
3219 // s.t. std.norm.(z) = t(z, df) at the two-tail probability level.
3220
3221 double A9, B9, T9, Z8, P7, B7, z;
3222
3223 A9 = df - 0.5;
3224 B9 = 48. * A9*A9,
3225 T9 = T*T / df;
3226 Z8 = T9 >= 0.04
3227 ? A9 * log(1. + T9)
3228 : A9 * (((1. - T9 * 0.75) * T9 / 3. - 0.5) * T9 + 1.) * T9;
3229 P7 = ((0.4 * Z8 + 3.3) * Z8 + 24.) * Z8 + 85.5;
3230 B7 = 0.8 * pow(Z8, 2.) + 100. + B9;
3231 z = (1. + (-P7 / B7 + Z8 + 3.) / B9) * sqrt(Z8);
3232
3233 return( z );
3234}
3235
3236//---------------------------------------------------------
3237double CSG_Test_Distribution::Get_T_Inv(double p, int df)
3238{ // Hill's approx. inverse t-dist.: Comm. of A.C.M Vol.13 No.10 1970 pg 620.
3239 // Calculates t given df and two-tail probability.
3240
3241 if( df == 1 )
3242 {
3243 return( cos(p * M_PI / 2.) / sin(p * M_PI / 2.) );
3244 }
3245
3246 if( df == 2 )
3247 {
3248 return( sqrt(2. / (p * (2. - p)) - 2.) );
3249 }
3250
3251 double a = 1. / (df - 0.5);
3252 double b = 48. / (a*a);
3253 double c = ((20700. * a / b - 98.) * a - 16.) * a + 96.36;
3254 double d = ((94.5 / (b + c) - 3.) / b + 1.) * sqrt(a * M_PI / 2.) * df;
3255 double x = d * p;
3256 double y = pow(x, 2. / df);
3257
3258 if( y > 0.05 + a )
3259 {
3260 x = Get_Norm_Z(0.5 * (1. - p));
3261 y = x*x;
3262
3263 if( df < 5 )
3264 {
3265 c = c + 0.3 * (df - 4.5) * (x + 0.6);
3266 }
3267
3268 c = (((0.05 * d * x - 5) * x - 7.) * x - 2.) * x + b + c;
3269 y = (((((0.4 * y + 6.3) * y + 36.) * y + 94.5) / c - y - 3.) / b + 1.) * x;
3270 y = a * y*y;
3271
3272 if( y > 0.002 )
3273 {
3274 y = exp(y) - 1.;
3275 }
3276 else
3277 {
3278 y = 0.5 * y*y + y;
3279 }
3280 }
3281 else
3282 {
3283 y = ((1. / (((df + 6.) / (df * y) - 0.089 * d - 0.822) * (df + 2.) * 3.)
3284 + 0.5 / (df + 4.)) * y - 1.) * (df + 1.) / (df + 2.) + 1. / y;
3285 }
3286
3287 return( sqrt(df * y) );
3288}
3289
3290
3292// //
3294
3295//---------------------------------------------------------
3296double CSG_Test_Distribution::Get_F_Tail_from_R2(double R2, int nPredictors, int nSamples, TSG_Test_Distribution_Type Type)
3297{
3298 double F = ((sLong)nSamples - (sLong)nPredictors - 1) * (R2 / nPredictors) / (1. - R2);
3299
3300 return( CSG_Test_Distribution::Get_F_Tail(F, nPredictors, nSamples - nPredictors - 1, Type) );
3301}
3302
3303//---------------------------------------------------------
3305{
3306 // calculates for F, dfn(ominator) and dfd(enominator) the "tail" of the F-distribution
3307
3308 double p = 1.;
3309
3310 if( F >= 0.00001 && dfn > 0 && dfd > 0 )
3311 {
3312 if( F * dfn >= dfd || F > 1. + 20. / dfn + 10. / sqrt((double)dfn) )
3313 {
3314 p = Get_Gamma(F, dfn, dfd);
3315 }
3316 else
3317 {
3318 p = 1. - Get_Gamma(1. / F, dfd, dfn);
3319 }
3320 }
3321
3322 if( p <= 0. || p >= 1. )
3323 {
3324 p = F > 1. ? 0. : F < 1. ? 1. : 0.5;
3325 }
3326
3327 return( Type == TESTDIST_TYPE_Right ? p : 1. - p );
3328}
3329
3330//---------------------------------------------------------
3331double CSG_Test_Distribution::Get_F_Inverse(double alpha, int dfn, int dfd, TSG_Test_Distribution_Type Type)
3332{
3333 if( alpha < 0. || alpha > 1. || dfd < 0 || dfn < 0 )
3334 {
3335 return( -1 );
3336 }
3337
3338 if( Type != TESTDIST_TYPE_Right )
3339 {
3340 alpha = 1. - alpha;
3341 }
3342
3343 const int ITERMAX = 100;
3344 const double EPSILON = 0.0001;
3345
3346 int i;
3347 double lo, hi, mid, p;
3348
3349 if( alpha <= 0.5 )
3350 {
3351 lo = 0.5;
3352 hi = lo;
3353
3354 for(i=0; i<ITERMAX; i++)
3355 {
3356 hi *= 2.;
3357 p = Get_F_Tail(hi, dfn, dfd);
3358
3359 if( p > alpha )
3360 {
3361 lo = hi;
3362 }
3363 else
3364 {
3365 break;
3366 }
3367 }
3368
3369 if( p > alpha )
3370 {
3371 return( hi );
3372 }
3373 }
3374 else
3375 {
3376 hi = 2;
3377 lo = hi;
3378
3379 for(i=0; i<ITERMAX; i++)
3380 {
3381 lo /= 2.;
3382 p = Get_F_Tail(lo, dfn, dfd);
3383
3384 if( p < alpha )
3385 {
3386 hi = lo;
3387 }
3388 else
3389 {
3390 break;
3391 }
3392 }
3393
3394 if( p < alpha )
3395 {
3396 return( lo );
3397 }
3398 }
3399
3400 mid = (hi + lo) / 2.;
3401
3402 for(i=0; i<ITERMAX && (hi-lo)>EPSILON*mid; i++)
3403 {
3404 mid = (hi + lo) / 2.;
3405 p = Get_F_Tail(mid, dfn, dfd);
3406
3407 if( p < alpha )
3408 hi = mid;
3409 else if( p > alpha )
3410 lo = mid;
3411 else
3412 break;
3413 }
3414
3415 return( mid );
3416}
3417
3418//---------------------------------------------------------
3419double CSG_Test_Distribution::Get_Gamma(double F, double dfn, double dfd)
3420{
3421 // calculates for F, dfn(ominator) and dfd(enominator) the incomplete Gamma-function
3422
3423 const double EXPMIN = -30.;
3424 const double SMALL = 0.00000000001;
3425
3426 double x, c, er, s, n, t1, t;
3427
3428 dfn /= 2.;
3429 dfd /= 2.;
3430
3431 x = dfd / (dfd + dfn * F);
3432 c = Get_Log_Gamma(dfn + dfd) - Get_Log_Gamma(dfn) - Get_Log_Gamma(dfd + 1.) + dfd * log(x) + dfn * log(1. - x);
3433
3434 if( c < EXPMIN )
3435 {
3436 return( -1. );
3437 }
3438
3439 dfn += dfd;
3440 dfd += 1.;
3441 c = exp(c);
3442 er = SMALL / c;
3443 t = dfn * x / dfd;
3444 t1 = 0.;
3445 s = t + 1.;
3446 n = 0;
3447
3448 while( t > er || t > t1 )
3449 {
3450 n += 1;
3451 t1 = t;
3452 t *= ((dfn + n) * x / (dfd + n));
3453 s += t;
3454 }
3455
3456 return( s * c );
3457}
3458
3459//---------------------------------------------------------
3460double CSG_Test_Distribution::Get_Log_Gamma(double a)
3461{
3462 // calculates the logarithm of the Gamma-function
3463
3464 const int ARGMIN = 6;
3465
3466 const double HL2PI = 0.91893853320467275; // = log(2. * M_PI) / 2.
3467
3468 int n = (int)floor(ARGMIN - a + 0.0001);
3469
3470 if( n > 0 )
3471 {
3472 a += n;
3473 }
3474
3475 double g;
3476
3477 g = 1. / (a*a);
3478 g = (1. - g * (1. / 30. - g * (1. / 105. - g * (1. / 140. - g / 99.)))) / (12. * a);
3479 g = g + ((a - 0.5) * log(a) - a + HL2PI);
3480
3481 for(int i=0; i<n; i++)
3482 {
3483 a = a - 1.;
3484 g = g - log(a);
3485 }
3486
3487 return( g );
3488}
3489
3490
3492// //
3493// //
3494// //
3496
3497//---------------------------------------------------------
3498CSG_Matrix SG_Get_Correlation_Matrix (const CSG_Matrix &Values, bool bCovariances)
3499{
3500 int nVariables = Values.Get_NX();
3501 int nSamples = Values.Get_NY();
3502
3503 //-----------------------------------------------------
3504 int i, j, k;
3506 CSG_Matrix C;
3507
3508 C.Create(nVariables, nVariables);
3509
3510 //-----------------------------------------------------
3511 S = new CSG_Simple_Statistics[nVariables];
3512
3513 for(j=0; j<nVariables; j++)
3514 {
3515 for(i=0; i<nSamples; i++)
3516 {
3517 S[j] += Values[i][j];
3518 }
3519 }
3520
3521 //-----------------------------------------------------
3522 for(k=0; k<nVariables; k++)
3523 {
3524 for(j=k; j<nVariables; j++)
3525 {
3526 double cov = 0.;
3527
3528 for(i=0; i<nSamples; i++)
3529 {
3530 cov += (Values[i][j] - S[j].Get_Mean()) * (Values[i][k] - S[k].Get_Mean());
3531 }
3532
3533 cov /= nSamples;
3534
3535 if( !bCovariances )
3536 {
3537 cov /= (S[j].Get_StdDev() * S[k].Get_StdDev());
3538 }
3539
3540 C[j][k] = C[k][j] = cov;
3541 }
3542 }
3543
3544 //-----------------------------------------------------
3545 delete[](S);
3546
3547 return( C );
3548}
3549
3550
3552// //
3553// //
3554// //
3556
3557//---------------------------------------------------------
bool SG_UI_Process_Get_Okay(bool bBlink)
void SG_UI_StatusBar_Set_Text(const CSG_String &Text)
SAGA_API_DLL_EXPORT void * SG_Malloc(size_t size)
signed long long sLong
Definition api_core.h:158
#define SG_T(s)
Definition api_core.h:537
SAGA_API_DLL_EXPORT void SG_Free(void *memblock)
SAGA_API_DLL_EXPORT void * SG_Realloc(void *memblock, size_t size)
SAGA_API_DLL_EXPORT void * SG_Calloc(size_t num, size_t size)
#define SG_FREE_SAFE(PTR)
Definition api_core.h:205
TSG_Data_Type
Definition api_core.h:1043
@ SG_DATATYPE_ULong
Definition api_core.h:1051
#define SG_Char
Definition api_core.h:536
#define _TL(s)
Definition api_core.h:1618
SAGA_API_DLL_EXPORT CSG_String SG_Get_String(double Value, int Precision=-99)
sLong Get_Size(void) const
Definition api_core.h:327
void * Get_Entry(sLong Index) const
Returns a pointer to the memory address of the requested variable. You have to type cast and derefere...
Definition api_core.h:331
void Create(TSG_Data_Type Type=SG_DATATYPE_String)
virtual ~CSG_Category_Statistics(void)
int Add_Value(int Value)
int Get_Category(int Value) const
int asInt(int iCategory) const
double asDouble(int iCategory) const
int Get_Count(void) const
CSG_Category_Statistics(TSG_Data_Type Type=SG_DATATYPE_String)
CSG_String asString(int iCategory) const
TSG_Data_Type Get_Category_Type(void) const
void Set_Probability_Relative(bool Value)
virtual ~CSG_Classifier_Supervised(void)
void Set_Threshold_Distance(double Value)
double Get_Threshold_Probability(void)
bool Get_WTA(int Method)
void Create(int nFeatures)
bool Add_Class(const CSG_String &Class_ID, const CSG_Vector &Mean, const CSG_Vector &Min, const CSG_Vector &Max, const CSG_Matrix &Cov)
void Set_Threshold_Probability(double Value)
bool Get_Probability_Relative(void)
bool Train(bool bClr_Samples=false)
const CSG_String & Get_Class_ID(int iClass)
Definition mat_tools.h:1267
void Set_Threshold_Angle(double Value)
bool Train_Add_Sample(const CSG_String &Class_ID, const CSG_Vector &Features)
double Get_Threshold_Distance(void)
static CSG_String Get_Name_of_Quality(int Method)
int Get_Class(const CSG_String &Class_ID)
void Set_WTA(int Method, bool bOn)
bool Save(const CSG_String &File, const SG_Char *Feature_Info=NULL)
static CSG_String Get_Name_of_Method(int Method)
bool Load(const CSG_String &File)
double Get_Threshold_Angle(void)
bool Create(int nFeatures)
int Get_nClusters(void) const
Definition mat_tools.h:1183
bool Execute(int Method, int nClusters, int nMaxIterations=0, int Initialization=0)
bool Set_Feature(sLong iElement, int iFeature, double Value)
sLong Get_nElements(void) const
Definition mat_tools.h:1181
virtual ~CSG_Cluster_Analysis(void)
virtual bool is_NoData(int x, int y) const
Definition grid.h:727
virtual double asDouble(sLong i, bool bScaled=true) const
Definition grid.h:789
sLong Get_NCells(void) const
Definition grid.h:565
virtual bool is_NoData(int x, int y, int z) const
Definition grids.h:375
sLong Get_NCells(void) const
Definition grids.h:191
virtual double asDouble(sLong i, bool bScaled=true) const
Definition grids.h:404
CSG_Histogram & operator=(const CSG_Histogram &Histogram)
bool Scale_Element_Count(double Scale)
virtual ~CSG_Histogram(void)
double Get_Center(int i) const
Definition mat_tools.h:1065
void Add_Value(double Value)
double Get_Value(double i) const
Definition mat_tools.h:1060
size_t Get_Class_Count(void) const
Definition mat_tools.h:1049
size_t Get_Cumulative(int i) const
Definition mat_tools.h:1057
bool Create(const CSG_Histogram &Histogram)
bool Destroy(void)
double Get_Quantile(double Quantile) const
double Get_Percentile_Value(double Value) const
double Get_Percentile(double Percentile) const
bool Add_Histogram(const CSG_Histogram &Histogram)
double Get_Break(int i) const
Definition mat_tools.h:1062
bool Update(void)
size_t Get_Element_Count(void) const
Definition mat_tools.h:1051
double Get_Quantile_Value(double Value) const
sLong Get_NRows(void) const
Definition mat_tools.h:525
bool from_String(const CSG_String &String)
int Get_NX(void) const
Definition mat_tools.h:522
sLong Get_NCols(void) const
Definition mat_tools.h:524
double Get_Determinant(void) const
int Get_NY(void) const
Definition mat_tools.h:523
bool is_Square(void) const
Definition mat_tools.h:538
CSG_String to_String(int Width=-1, int Precision=-1, bool bScientific=false, const SG_Char *Separator=NULL) const
CSG_Matrix Get_Inverse(bool bSilent=true, int nSubSquare=0) const
bool Cmp_Name(const CSG_String &String, bool bNoCase=true) const
Definition metadata.cpp:484
bool Save(const CSG_String &File, const SG_Char *Extension=NULL) const
Definition metadata.cpp:904
int Get_Children_Count(void) const
Definition metadata.h:154
CSG_MetaData * Get_Child(int Index) const
Definition metadata.h:155
const CSG_String & Get_Content(void) const
Definition metadata.h:139
const SG_Char * Get_Property(int Index) const
Definition metadata.h:188
void Set_Name(const CSG_String &Name)
Definition metadata.h:136
CSG_MetaData * Add_Child(void)
Definition metadata.cpp:166
bool Load(const CSG_String &File, const SG_Char *Extension=NULL)
Definition metadata.cpp:809
bool Add_Property(const CSG_String &Name, const CSG_String &Value)
Definition metadata.cpp:582
virtual ~CSG_Natural_Breaks(void)
bool Create(class CSG_Table *pTable, int Field, int nClasses, int Histogram=0)
int Get_Count(void) const
Definition mat_tools.h:1128
CSG_Random(void)
static double Get_Gaussian(double mean, double stddev)
static void Initialize(void)
static double Get_Uniform(void)
sLong Get_IndexOfMinimum(void)
void Add(const CSG_Simple_Statistics &Statistics)
double Get_Percentile(double Percentile)
double Get_Median(void)
double Get_Value(sLong i) const
Definition mat_tools.h:777
double Get_Gini(void)
sLong Get_nValues_Above(double Threshold, bool bEquals=false)
sLong Get_nValues_Below(double Threshold, bool bEquals=false)
int is_Evaluated(void) const
Definition mat_tools.h:741
double Get_Maximum(void)
Definition mat_tools.h:749
double * Get_Values(void) const
Definition mat_tools.h:776
double Get_Mean(void)
Definition mat_tools.h:753
double Get_Sum(void)
Definition mat_tools.h:751
double Get_SkewnessPearson(void)
sLong Get_IndexOfMaximum(void)
bool Set_Count(sLong Count)
void _Evaluate(int Level=1)
bool Create(bool bHoldValues=false)
double Get_StdDev(void)
Definition mat_tools.h:755
sLong Get_Count(void) const
Definition mat_tools.h:745
void Add_Value(double Value, double Weight=1.)
double Get_Minimum(void)
Definition mat_tools.h:748
double Get_Quantile(double Quantile)
int Cmp(const CSG_String &String) const
static CSG_String Format(const char *Format,...)
bool is_Empty(void) const
bool Create(const class wxString *pString)
sLong Get_Index(void) const
Definition table.h:136
double asDouble(int Field) const
int asInt(int Field) const
Definition table.h:219
bool Add_Value(int Field, const double &Value)
const SG_Char * asString(int Field, int Decimals=-99) const
sLong Get_Count(void) const
Definition table.h:410
virtual bool Get_Value(sLong Index, int Field, CSG_String &Value) const
Definition table.cpp:1159
bool Set_Index(CSG_Index &Index, int Field, bool bAscending=true) const
Definition table.cpp:1485
int Get_Field_Count(void) const
Definition table.h:371
CSG_Table_Record * Get_Record_byIndex(sLong Index) const
Definition table.h:417
static double Get_Norm_P(double Z)
static double Get_F_Inverse(double alpha, int dfn, int dfd, TSG_Test_Distribution_Type Type=TESTDIST_TYPE_Right)
static double Get_F_Tail_from_R2(double R2, int nPredictors, int nSamples, TSG_Test_Distribution_Type Type=TESTDIST_TYPE_Right)
static double Get_T_Tail(double T, int df, TSG_Test_Distribution_Type Type=TESTDIST_TYPE_Right)
static double Get_F_Tail(double F, int dfn, int dfd, TSG_Test_Distribution_Type Type=TESTDIST_TYPE_Right)
static double Get_Norm_Z(double P)
static double Get_T_Inverse(double alpha, int df, TSG_Test_Distribution_Type Type=TESTDIST_TYPE_Right)
void Add_Value(double Value, double Weight=1.)
virtual void Create(bool bWeights=false)
int Get_Class_Index(double Value) const
void Add_Value(const CSG_String &Value, double Weight=1.)
virtual void Create(bool bWeights=false)
int Get_Class_Index(const CSG_String &Value) const
virtual int Get_Minority(bool bWeighted=false) const
int Get_Count(void) const
Definition mat_tools.h:819
virtual int Get_Majority(bool bWeighted=false) const
CSG_String to_String(int Width=-1, int Precision=-1, bool bScientific=false, const SG_Char *Separator=NULL) const
sLong Get_Size(void) const
Definition mat_tools.h:382
int Get_N(void) const
Definition mat_tools.h:384
double Get_Angle(const CSG_Vector &Vector) const
bool from_String(const CSG_String &String)
#define C
#define EPSILON
int SG_Compare_Double(const void *a, const void *b)
double SG_Degree_To_Decimal(double Deg, double Min, double Sec)
int SG_Compare_Char_Ptr(const void *a, const void *b)
double SG_Get_Rounded(double Value, int Decimals)
Definition mat_tools.cpp:83
double SG_Get_Rounded_To_SignificantFigures(double Value, int Decimals)
int SG_Get_Digit_Count(int Number)
CSG_String SG_Get_Double_asString(double Number, int Width, int Precision, bool bScientific)
int SG_Compare_Int(const void *a, const void *b)
void SG_Decimal_To_Degree(double Value, double &Deg, double &Min, double &Sec)
double SG_Get_Square(double Value)
Definition mat_tools.cpp:70
TSG_Test_Distribution_Type
Definition mat_tools.h:1530
@ TESTDIST_TYPE_Right
Definition mat_tools.h:1532
@ TESTDIST_TYPE_Left
Definition mat_tools.h:1531
@ TESTDIST_TYPE_Middle
Definition mat_tools.h:1533
@ TESTDIST_TYPE_TwoTail
Definition mat_tools.h:1534
SAGA_API_DLL_EXPORT CSG_Matrix SG_Get_Correlation_Matrix(const CSG_Matrix &Values, bool bCovariances=false)
#define M_RAD_TO_DEG
Definition mat_tools.h:108
@ SG_CLASSIFY_SUPERVISED_ParallelEpiped
Definition mat_tools.h:1231
@ SG_CLASSIFY_SUPERVISED_SVM
Definition mat_tools.h:1238
@ SG_CLASSIFY_SUPERVISED_MinimumDistance
Definition mat_tools.h:1232
@ SG_CLASSIFY_SUPERVISED_SAM
Definition mat_tools.h:1235
@ SG_CLASSIFY_SUPERVISED_SID
Definition mat_tools.h:1237
@ SG_CLASSIFY_SUPERVISED_Mahalonobis
Definition mat_tools.h:1233
@ SG_CLASSIFY_SUPERVISED_WTA
Definition mat_tools.h:1236
@ SG_CLASSIFY_SUPERVISED_MaximumLikelihood
Definition mat_tools.h:1234
@ SG_CLASSIFY_SUPERVISED_BinaryEncoding
Definition mat_tools.h:1230
SAGA_API_DLL_EXPORT int SG_Compare_Double(const void *a, const void *b)
SAGA_API_DLL_EXPORT double SG_Get_Square(double Value)
Definition mat_tools.cpp:70
#define M_PI
Definition mat_tools.h:96
int SG_Compare_Version(const CSG_String &Version, int Major, int Minor, int Release)
Definition saga_api.cpp:82
#define SAGA_VERSION
@ TABLE_INDEX_Ascending
Definition table.h:105