1 #include <stdlib.h>
2 #include <grass/gis.h>
3 #include <grass/raster.h>
4 #include <grass/imagery.h>
5 #include <grass/segment.h> /* segmentation library */
6 #include <grass/glocale.h>
7 #include "iseg.h"
8
write_ids(struct globals * globals)9 int write_ids(struct globals *globals)
10 {
11 int out_fd, row, col, maxid;
12 CELL *outbuf, rid;
13 struct Colors colors;
14 struct History hist;
15
16 outbuf = Rast_allocate_c_buf();
17
18 G_debug(1, "preparing output raster");
19 /* open output raster map */
20 out_fd = Rast_open_new(globals->out_name, CELL_TYPE);
21
22 G_debug(1, "start data transfer from segmentation file to raster");
23
24 G_message(_("Writing out segment IDs..."));
25 maxid = 0;
26 for (row = 0; row < globals->nrows; row++) {
27
28 G_percent(row, globals->nrows, 9);
29
30 Rast_set_c_null_value(outbuf, globals->ncols);
31 for (col = 0; col < globals->ncols; col++) {
32
33 if (!(FLAG_GET(globals->null_flag, row, col))) {
34 Segment_get(&globals->rid_seg, (void *) &rid, row, col);
35
36 if (rid > 0) {
37 if (globals->method == ORM_RG)
38 rid = globals->new_id[rid];
39 outbuf[col] = rid;
40 if (maxid < rid)
41 maxid = rid;
42 }
43 }
44 }
45 Rast_put_row(out_fd, outbuf, CELL_TYPE);
46 }
47 G_percent(1, 1, 1);
48
49 /* close and save segment id file */
50 Rast_close(out_fd);
51 G_free(outbuf);
52
53 /* set colors */
54 Rast_init_colors(&colors);
55 Rast_make_random_colors(&colors, 1, maxid);
56 Rast_write_colors(globals->out_name, G_mapset(), &colors);
57
58 Rast_short_history(globals->out_name, "raster", &hist);
59 Rast_command_history(&hist);
60 Rast_write_history(globals->out_name, &hist);
61
62 /* free memory */
63 Rast_free_colors(&colors);
64
65 return TRUE;
66 }
67
68 /* write goodness of fit */
write_gof_rg(struct globals * globals)69 int write_gof_rg(struct globals *globals)
70 {
71 int row, col;
72 int mean_fd;
73 CELL rid;
74 FCELL *meanbuf;
75 double thresh, maxdev, sim, mingood;
76 struct ngbr_stats Ri, Rk;
77 DCELL **inbuf; /* buffers to store lines from each of the imagery group rasters */
78 int n, *in_fd;
79 struct FPRange *fp_range; /* min/max values of each input raster */
80 struct Colors colors;
81 struct History hist;
82 DCELL *min, *max;
83
84 mean_fd = Rast_open_new(globals->gof, FCELL_TYPE);
85 meanbuf = Rast_allocate_f_buf();
86
87 /* goodness of fit for each cell: 1 = good fit, 0 = bad fit */
88 /* similarity of each cell to region mean
89 * max possible difference: globals->threshold
90 * if similarity < globals->alpha * globals->alpha * globals->threshold
91 * 1
92 * else
93 * (similarity - globals->alpha * globals->alpha * globals->threshold) /
94 * (globals->threshold * (1 - globals->alpha * globals->alpha) */
95
96 thresh = globals->alpha * globals->alpha * globals->max_diff;
97 maxdev = globals->max_diff * (1 - globals->alpha * globals->alpha);
98 mingood = 1;
99
100 /* open input bands */
101 in_fd = G_malloc(globals->Ref.nfiles * sizeof(int));
102 inbuf = (DCELL **) G_malloc(globals->Ref.nfiles * sizeof(DCELL *));
103 fp_range = G_malloc(globals->Ref.nfiles * sizeof(struct FPRange));
104 min = G_malloc(globals->Ref.nfiles * sizeof(DCELL));
105 max = G_malloc(globals->Ref.nfiles * sizeof(DCELL));
106
107 G_debug(1, "Opening input rasters...");
108 for (n = 0; n < globals->Ref.nfiles; n++) {
109 inbuf[n] = Rast_allocate_d_buf();
110 in_fd[n] = Rast_open_old(globals->Ref.file[n].name, globals->Ref.file[n].mapset);
111
112 /* returns -1 on error, 2 on empty range, quitting either way. */
113 if (Rast_read_fp_range(globals->Ref.file[n].name, globals->Ref.file[n].mapset, &fp_range[n]) != 1)
114 G_fatal_error(_("No min/max found in raster map <%s>"),
115 globals->Ref.file[n].name);
116 Rast_get_fp_range_min_max(&(fp_range[n]), &min[n], &max[n]);
117
118 G_debug(1, "Range for layer %d: min = %f, max = %f",
119 n, min[n], max[n]);
120 }
121
122 G_message(_("Writing out goodness of fit"));
123 for (row = 0; row < globals->nrows; row++) {
124
125 G_percent(row, globals->nrows, 9);
126
127 Rast_set_f_null_value(meanbuf, globals->ncols);
128
129 for (n = 0; n < globals->Ref.nfiles; n++) {
130 Rast_get_d_row(in_fd[n], inbuf[n], row);
131 }
132
133 for (col = 0; col < globals->ncols; col++) {
134
135 if (!(FLAG_GET(globals->null_flag, row, col))) {
136
137 Segment_get(&globals->rid_seg, (void *) &rid, row, col);
138
139 if (rid > 0) {
140
141 Ri.row = Rk.row = row;
142 Ri.col = Rk.col = col;
143
144 /* get values for Ri = this region */
145 globals->rs.id = rid;
146 fetch_reg_stats(row, col, &globals->rs, globals);
147 Ri.mean = globals->rs.mean;
148 Ri.count = globals->rs.count;
149
150 sim = 0.;
151 /* region consists of more than one cell */
152 if (Ri.count > 1) {
153
154 /* get values for Rk = this cell */
155 for (n = 0; n < globals->Ref.nfiles; n++) {
156 if (globals->weighted == FALSE)
157 /* scaled version */
158 globals->second_val[n] = (inbuf[n][col] - min[n]) / (max[n] - min[n]);
159 else
160 globals->second_val[n] = inbuf[n][col];
161 }
162
163 Rk.mean = globals->second_val;
164
165 /* calculate similarity */
166 sim = (*globals->calculate_similarity) (&Ri, &Rk, globals);
167 }
168
169 if (0) {
170 if (sim < thresh)
171 meanbuf[col] = 1;
172 else {
173 sim = 1. - (sim - thresh) / maxdev;
174 meanbuf[col] = sim;
175 if (mingood > sim)
176 mingood = sim;
177 }
178 }
179 else {
180 sim = 1 - sim;
181 meanbuf[col] = sim;
182 if (mingood > sim)
183 mingood = sim;
184 }
185 }
186 }
187 }
188 Rast_put_row(mean_fd, meanbuf, FCELL_TYPE);
189 }
190
191 Rast_close(mean_fd);
192
193 Rast_init_colors(&colors);
194 Rast_make_grey_scale_fp_colors(&colors, mingood, 1);
195 Rast_write_colors(globals->gof, G_mapset(), &colors);
196
197 Rast_short_history(globals->gof, "raster", &hist);
198 Rast_command_history(&hist);
199 Rast_write_history(globals->gof, &hist);
200
201 G_free(meanbuf);
202
203 G_debug(1, "Closing input rasters...");
204 for (n = 0; n < globals->Ref.nfiles; n++) {
205 Rast_close(in_fd[n]);
206 G_free(inbuf[n]);
207 }
208
209 G_free(inbuf);
210 G_free(in_fd);
211 G_free(fp_range);
212 G_free(min);
213 G_free(max);
214
215 return TRUE;
216 }
217
write_bands_ms(struct globals * globals)218 int write_bands_ms(struct globals *globals)
219 {
220 int *out_fd, row, col, n;
221 DCELL **outbuf;
222 char **name;
223 struct Colors colors;
224 struct History hist;
225 struct ngbr_stats Rout;
226
227 out_fd = G_malloc(sizeof(int) * globals->nbands);
228 name = G_malloc(sizeof(char *) * globals->nbands);
229 outbuf = G_malloc(sizeof(DCELL) * globals->nbands);
230 for (n = 0; n < globals->nbands; n++) {
231 outbuf[n] = Rast_allocate_d_buf();
232 G_asprintf(&name[n], "%s%s", globals->Ref.file[n].name, globals->bsuf);
233 out_fd[n] = Rast_open_new(name[n], DCELL_TYPE);
234 }
235
236 Rout.mean = G_malloc(globals->datasize);
237
238 G_message(_("Writing out shifted band values..."));
239
240 for (row = 0; row < globals->nrows; row++) {
241
242 G_percent(row, globals->nrows, 9);
243
244 for (n = 0; n < globals->nbands; n++)
245 Rast_set_d_null_value(outbuf[n], globals->ncols);
246 for (col = 0; col < globals->ncols; col++) {
247
248 if (!(FLAG_GET(globals->null_flag, row, col))) {
249 Segment_get(globals->bands_out, (void *) Rout.mean, row, col);
250
251 for (n = 0; n < globals->nbands; n++) {
252 outbuf[n][col] = Rout.mean[n];
253 if (globals->weighted == FALSE)
254 outbuf[n][col] = Rout.mean[n] * (globals->max[n] - globals->min[n]) + globals->min[n];
255 }
256 }
257 }
258 for (n = 0; n < globals->nbands; n++)
259 Rast_put_row(out_fd[n], outbuf[n], DCELL_TYPE);
260 }
261
262 for (n = 0; n < globals->nbands; n++) {
263 Rast_close(out_fd[n]);
264
265 Rast_read_colors(globals->Ref.file[n].name, globals->Ref.file[n].mapset, &colors);
266 Rast_write_colors(name[n], G_mapset(), &colors);
267
268 Rast_short_history(name[n], "raster", &hist);
269 Rast_command_history(&hist);
270 Rast_write_history(name[n], &hist);
271 }
272
273 /* free */
274
275 return TRUE;
276 }
277
close_files(struct globals * globals)278 int close_files(struct globals *globals)
279 {
280 G_debug(1, "Closing input rasters...");
281
282 /* close segmentation files and output raster */
283 G_debug(1, "closing files");
284 Segment_close(&globals->bands_seg);
285 if (globals->method == ORM_MS)
286 Segment_close(&globals->bands_seg2);
287 if (globals->bounds_map)
288 Segment_close(&globals->bounds_seg);
289
290 G_free(globals->bands_val);
291 G_free(globals->second_val);
292
293 Segment_close(&globals->rid_seg);
294
295 flag_destroy(globals->null_flag);
296 flag_destroy(globals->candidate_flag);
297
298 rgtree_destroy(globals->reg_tree);
299
300 /* anything else left to clean up? */
301
302 return TRUE;
303 }
304