new distances
netzel

netzel commited on 2024-03-25 09:42:52
Showing 3 changed files, with 62 additions and 12 deletions.

... ...
@@ -32,15 +32,40 @@ void init_cumul_sum(double* cumul_sum, long index, int nodata,
32 32
 }
33 33
 
34 34
 double calc_distance(double** layers, int n_layers, double* cumul_sum,
35
-                     long index, long superpixel_size) {
36
-	double diff, dist = 0.0;
35
+                     long index, long superpixel_size, int distance, double minkowski_p) {
36
+    double diff, dist = 0.0, A = 1.0, B = 1.0, AB = 1.0;
37 37
     int i;
38 38
     
39
+    switch(distance) {
40
+        case 0: // minkowski
39 41
             for(i=0; i<n_layers; i++) {
40
-		diff = cumul_sum[i] - layers[i][index]*superpixel_size;
41
-		dist += diff*diff;
42
+                diff = fabs(cumul_sum[i] - layers[i][index]*superpixel_size);
43
+                dist += pow(diff, minkowski_p);
42 44
             }
43
-	return sqrt(dist)/superpixel_size;
45
+            dist = pow(dist, 1.0/minkowski_p)/superpixel_size;
46
+            break;
47
+        case 1: // cosine
48
+            for(i=0; i<n_layers; i++) {
49
+                diff = cumul_sum[i]/superpixel_size;
50
+                A += pow(diff, 2.0);
51
+                B += pow(layers[i][index],2.0);
52
+                AB += layers[i][index]*diff;
53
+            }
54
+            dist = AB/sqrt(A*B);
55
+            break;
56
+        case 2: // angular
57
+            for(i=0; i<n_layers; i++) {
58
+                diff = cumul_sum[i]/superpixel_size;
59
+                A += pow(diff, 2.0);
60
+                B += pow(layers[i][index],2.0);
61
+                AB += layers[i][index]*diff;
62
+            }
63
+            dist = acos(AB/sqrt(A*B))/M_PI;
64
+            break;
65
+        default:
66
+            dist = 0.0;
67
+    }
68
+    return dist;
44 69
 }
45 70
 
46 71
 void get_cell_neighbour(CELL* c, int n, int* x, int* y) {
... ...
@@ -55,7 +80,7 @@ void grow_adaptel(int i, LABELS* labels, MINLIST* minlist, SEEDS* seeds,
55 80
                   double* distances, double** layers, int n_layers,
56 81
                   unsigned char* mask, int cols, int rows,
57 82
                   long size, double threshold, int connectivity,
58
-		          int* dIdx) {
83
+                  int* dIdx, int distance, double minkowski_p) {
59 84
 
60 85
     long idx1, idx2, index, superpixel_size;
61 86
     double* cumul_sum = (double*)malloc(sizeof(double)*n_layers);
... ...
@@ -91,7 +116,8 @@ void grow_adaptel(int i, LABELS* labels, MINLIST* minlist, SEEDS* seeds,
91 116
                             else
92 117
                                 dist = distances[idx1] + 
93 118
                                           calc_distance(layers, n_layers, cumul_sum,
94
-                                                        idx2, superpixel_size);
119
+                                                        idx2, superpixel_size,
120
+                                                        distance, minkowski_p);
95 121
 
96 122
                             if(dist < distances[idx2] || labels->labels[idx2] < 0) {
97 123
                                 distances[idx2] = dist;
... ...
@@ -150,7 +176,9 @@ LABELS* create_adaptels(
150 176
              int            cols,
151 177
              int            rows,
152 178
              const double   threshold,
153
-             const int      connectivity) {
179
+             const int      connectivity,
180
+             const int      distance,
181
+             const double   minkowski_p) {
154 182
     
155 183
     const long size = (long)cols*rows;
156 184
     long i;
... ...
@@ -171,7 +199,8 @@ LABELS* create_adaptels(
171 199
     for(i = 0; i < seeds->N; i++)
172 200
         grow_adaptel(i, labels, minlist, seeds, distances, 
173 201
                      layers, n_layers, mask, cols, rows,
174
-		             size, threshold, connectivity, dIdx);
202
+                     size, threshold, connectivity, dIdx,
203
+                     distance, minkowski_p);
175 204
 
176 205
     free_minlist(minlist);
177 206
     free_seeds(seeds);
... ...
@@ -10,6 +10,8 @@ LABELS* create_adaptels(
10 10
              int            cols, 
11 11
              int            rows,
12 12
              const double   threshold,
13
-             const int      connectivity);
13
+             const int      connectivity,
14
+             const int      distance,
15
+             const double   minkowski_p);
14 16
 
15 17
 #endif
... ...
@@ -43,18 +43,23 @@ void usage(char *progname, void *argtable) {
43 43
 int main(int argc, char **argv) {
44 44
 
45 45
     long i;
46
+    char* dist_names[] = {"minkowski","cosine","angular"};
46 47
 
47 48
     double thr_val = 60.0;
48 49
     int con_val = 4;
50
+    int distance_val = 0;
51
+    double minkowski_p_val = 2.0;
49 52
 
50 53
     struct arg_str  *inp  = arg_strn("i","input","<file_name>",1,99,"name of input file(s) (GeoTIFF)");
51 54
     struct arg_str  *out  = arg_str1("o","output","<file_name>","name of output file (GeoTIFF) with adaptels");
52 55
     struct arg_dbl  *thr  = arg_dbl0("t","threshold","<val>","threshold value for adaptels algorithm (default: 60.0)");
56
+    struct arg_str  *dist = arg_str0("d","distance","<distance name>","name of the distance measure: minkowki, cosine, angular (default: minkowski)");
57
+    struct arg_dbl  *mink = arg_dbl0("p","minkowski_p","<val>","value for minkowski distance (default: 2.0)");
53 58
     struct arg_lit  *con  = arg_lit0("8","Queen topology","cells connectivity (default: 4, Rook)");
54 59
     struct arg_lit  *norm = arg_lit0("n","normalize","normalize cells values to the interval [0,1] (default: no)");
55 60
     struct arg_lit  *help = arg_lit0("h","help","print this help and exit");
56 61
     struct arg_end  *end  = arg_end(20);
57
-    void* argtable[] = {inp,out,thr,con,norm,help,end};
62
+    void* argtable[] = {inp,out,thr,dist,mink,con,norm,help,end};
58 63
 
59 64
     int nerrors = arg_parse(argc,argv,argtable);
60 65
 
... ...
@@ -74,6 +79,16 @@ int main(int argc, char **argv) {
74 79
       }
75 80
     }
76 81
     
82
+    if(mink->count>0)
83
+      minkowski_p_val = mink->dval[0];
84
+
85
+    if(dist->count>0) {
86
+      if(strcmp(dist_names[1],dist->sval[0])==0)
87
+        distance_val = 1;
88
+      else if(strcmp(dist_names[2],dist->sval[0])==0)
89
+        distance_val = 2;
90
+    }
91
+
77 92
     if(con->count>0)
78 93
       con_val = 8;
79 94
 
... ...
@@ -163,6 +178,8 @@ int main(int argc, char **argv) {
163 178
       ezgdal_show_message(stderr,"OK");
164 179
     }
165 180
     
181
+    ezgdal_show_message_nolf(stderr,"Distance:             ");
182
+    ezgdal_show_message(stderr,dist_names[distance_val]);
166 183
     ezgdal_show_message_nolf(stderr,"Adaptels creating:    ");
167 184
 
168 185
     clock_t time_start, time_end;
... ...
@@ -172,7 +189,9 @@ int main(int argc, char **argv) {
172 189
     LABELS* labels = create_adaptels(input_buffers, inp->count, mask,
173 190
                                      input_layers[0]->cols, 
174 191
                                      input_layers[0]->rows, 
175
-                                     thr_val, con_val);
192
+                                     thr_val, con_val,
193
+                                     distance_val,
194
+                                     minkowski_p_val);
176 195
     time_end = clock();
177 196
     time_elapsed = ((double) (time_end - time_start)) / CLOCKS_PER_SEC;
178 197
     fprintf(stderr,"%d adaptels created in %.4lf seconds\n",labels->max_label, time_elapsed);
179 198