1 /*
2 *
3 * WKTRaster - Raster Types for PostGIS
4 * http://trac.osgeo.org/postgis/wiki/WKTRaster
5 *
6 * Copyright (C) 2011-2013 Regents of the University of California
7 * <bkpark@ucdavis.edu>
8 * Copyright (C) 2010-2011 Jorge Arevalo <jorge.arevalo@deimos-space.com>
9 * Copyright (C) 2010-2011 David Zwarg <dzwarg@azavea.com>
10 * Copyright (C) 2009-2011 Pierre Racine <pierre.racine@sbf.ulaval.ca>
11 * Copyright (C) 2009-2011 Mateusz Loskot <mateusz@loskot.net>
12 * Copyright (C) 2008-2009 Sandro Santilli <strk@kbt.io>
13 *
14 * This program is free software; you can redistribute it and/or
15 * modify it under the terms of the GNU General Public License
16 * as published by the Free Software Foundation; either version 2
17 * of the License, or (at your option) any later version.
18 *
19 * This program is distributed in the hope that it will be useful,
20 * but WITHOUT ANY WARRANTY; without even the implied warranty of
21 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22 * GNU General Public License for more details.
23 *
24 * You should have received a copy of the GNU General Public License
25 * along with this program; if not, write to the Free Software Foundation,
26 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
27 *
28 */
29
30 #include "librtcore.h"
31 #include "librtcore_internal.h"
32
33 /******************************************************************************
34 * rt_band_reclass()
35 ******************************************************************************/
36
37 /**
38 * Returns new band with values reclassified
39 *
40 * @param srcband : the band who's values will be reclassified
41 * @param pixtype : pixel type of the new band
42 * @param hasnodata : indicates if the band has a nodata value
43 * @param nodataval : nodata value for the new band
44 * @param exprset : array of rt_reclassexpr structs
45 * @param exprcount : number of elements in expr
46 *
47 * @return a new rt_band or NULL on error
48 */
49 rt_band
rt_band_reclass(rt_band srcband,rt_pixtype pixtype,uint32_t hasnodata,double nodataval,rt_reclassexpr * exprset,int exprcount)50 rt_band_reclass(
51 rt_band srcband, rt_pixtype pixtype,
52 uint32_t hasnodata, double nodataval,
53 rt_reclassexpr *exprset, int exprcount
54 ) {
55 rt_band band = NULL;
56 uint32_t width = 0;
57 uint32_t height = 0;
58 int numval = 0;
59 int memsize = 0;
60 void *mem = NULL;
61 uint32_t src_hasnodata = 0;
62 double src_nodataval = 0.0;
63 int isnodata = 0;
64
65 int rtn;
66 uint32_t x;
67 uint32_t y;
68 int i;
69 double or = 0;
70 double ov = 0;
71 double nr = 0;
72 double nv = 0;
73 int do_nv = 0;
74 rt_reclassexpr expr = NULL;
75
76 assert(NULL != srcband);
77 assert(NULL != exprset && exprcount > 0);
78 RASTER_DEBUGF(4, "exprcount = %d", exprcount);
79 RASTER_DEBUGF(4, "exprset @ %p", exprset);
80
81 /* source nodata */
82 src_hasnodata = rt_band_get_hasnodata_flag(srcband);
83 if (src_hasnodata)
84 rt_band_get_nodata(srcband, &src_nodataval);
85
86 /* size of memory block to allocate */
87 width = rt_band_get_width(srcband);
88 height = rt_band_get_height(srcband);
89 numval = width * height;
90 memsize = rt_pixtype_size(pixtype) * numval;
91 mem = (int *) rtalloc(memsize);
92 if (!mem) {
93 rterror("rt_band_reclass: Could not allocate memory for band");
94 return 0;
95 }
96
97 /* initialize to zero */
98 if (!hasnodata) {
99 memset(mem, 0, memsize);
100 }
101 /* initialize to nodataval */
102 else {
103 int32_t checkvalint = 0;
104 uint32_t checkvaluint = 0;
105 double checkvaldouble = 0;
106 float checkvalfloat = 0;
107
108 switch (pixtype) {
109 case PT_1BB:
110 {
111 uint8_t *ptr = mem;
112 uint8_t clamped_initval = rt_util_clamp_to_1BB(nodataval);
113 for (i = 0; i < numval; i++)
114 ptr[i] = clamped_initval;
115 checkvalint = ptr[0];
116 break;
117 }
118 case PT_2BUI:
119 {
120 uint8_t *ptr = mem;
121 uint8_t clamped_initval = rt_util_clamp_to_2BUI(nodataval);
122 for (i = 0; i < numval; i++)
123 ptr[i] = clamped_initval;
124 checkvalint = ptr[0];
125 break;
126 }
127 case PT_4BUI:
128 {
129 uint8_t *ptr = mem;
130 uint8_t clamped_initval = rt_util_clamp_to_4BUI(nodataval);
131 for (i = 0; i < numval; i++)
132 ptr[i] = clamped_initval;
133 checkvalint = ptr[0];
134 break;
135 }
136 case PT_8BSI:
137 {
138 int8_t *ptr = mem;
139 int8_t clamped_initval = rt_util_clamp_to_8BSI(nodataval);
140 for (i = 0; i < numval; i++)
141 ptr[i] = clamped_initval;
142 checkvalint = ptr[0];
143 break;
144 }
145 case PT_8BUI:
146 {
147 uint8_t *ptr = mem;
148 uint8_t clamped_initval = rt_util_clamp_to_8BUI(nodataval);
149 for (i = 0; i < numval; i++)
150 ptr[i] = clamped_initval;
151 checkvalint = ptr[0];
152 break;
153 }
154 case PT_16BSI:
155 {
156 int16_t *ptr = mem;
157 int16_t clamped_initval = rt_util_clamp_to_16BSI(nodataval);
158 for (i = 0; i < numval; i++)
159 ptr[i] = clamped_initval;
160 checkvalint = ptr[0];
161 break;
162 }
163 case PT_16BUI:
164 {
165 uint16_t *ptr = mem;
166 uint16_t clamped_initval = rt_util_clamp_to_16BUI(nodataval);
167 for (i = 0; i < numval; i++)
168 ptr[i] = clamped_initval;
169 checkvalint = ptr[0];
170 break;
171 }
172 case PT_32BSI:
173 {
174 int32_t *ptr = mem;
175 int32_t clamped_initval = rt_util_clamp_to_32BSI(nodataval);
176 for (i = 0; i < numval; i++)
177 ptr[i] = clamped_initval;
178 checkvalint = ptr[0];
179 break;
180 }
181 case PT_32BUI:
182 {
183 uint32_t *ptr = mem;
184 uint32_t clamped_initval = rt_util_clamp_to_32BUI(nodataval);
185 for (i = 0; i < numval; i++)
186 ptr[i] = clamped_initval;
187 checkvaluint = ptr[0];
188 break;
189 }
190 case PT_32BF:
191 {
192 float *ptr = mem;
193 float clamped_initval = rt_util_clamp_to_32F(nodataval);
194 for (i = 0; i < numval; i++)
195 ptr[i] = clamped_initval;
196 checkvalfloat = ptr[0];
197 break;
198 }
199 case PT_64BF:
200 {
201 double *ptr = mem;
202 for (i = 0; i < numval; i++)
203 ptr[i] = nodataval;
204 checkvaldouble = ptr[0];
205 break;
206 }
207 default:
208 {
209 rterror("rt_band_reclass: Unknown pixeltype %d", pixtype);
210 rtdealloc(mem);
211 return 0;
212 }
213 }
214
215 /* Overflow checking */
216 rt_util_dbl_trunc_warning(
217 nodataval,
218 checkvalint, checkvaluint,
219 checkvalfloat, checkvaldouble,
220 pixtype
221 );
222 }
223 RASTER_DEBUGF(3, "rt_band_reclass: width = %d height = %d", width, height);
224
225 band = rt_band_new_inline(width, height, pixtype, hasnodata, nodataval, mem);
226 if (!band) {
227 rterror("rt_band_reclass: Could not create new band");
228 rtdealloc(mem);
229 return 0;
230 }
231 rt_band_set_ownsdata_flag(band, 1); /* we DO own this data!!! */
232 RASTER_DEBUGF(3, "rt_band_reclass: new band @ %p", band);
233
234 for (x = 0; x < width; x++) {
235 for (y = 0; y < height; y++) {
236 rtn = rt_band_get_pixel(srcband, x, y, &ov, &isnodata);
237
238 /* error getting value, skip */
239 if (rtn != ES_NONE) {
240 RASTER_DEBUGF(3, "Cannot get value at %d, %d", x, y);
241 continue;
242 }
243 RASTER_DEBUGF(4, "(x, y, ov, isnodata) = (%d, %d, %f, %d)", x, y, ov, isnodata);
244
245 do {
246 do_nv = 0;
247
248 /* no data*/
249 if (hasnodata && isnodata) {
250 do_nv = 1;
251 break;
252 }
253
254 for (i = 0; i < exprcount; i++) {
255 expr = exprset[i];
256
257 /* ov matches min and max*/
258 if (
259 FLT_EQ(expr->src.min, ov) &&
260 FLT_EQ(expr->src.max, ov)
261 ) {
262 do_nv = 1;
263 break;
264 }
265
266 /* process min */
267 if ((
268 expr->src.exc_min && (
269 expr->src.min > ov ||
270 FLT_EQ(expr->src.min, ov)
271 )) || (
272 expr->src.inc_min && (
273 expr->src.min < ov ||
274 FLT_EQ(expr->src.min, ov)
275 )) || (
276 expr->src.min < ov
277 )) {
278 /* process max */
279 if ((
280 expr->src.exc_max && (
281 ov > expr->src.max ||
282 FLT_EQ(expr->src.max, ov)
283 )) || (
284 expr->src.inc_max && (
285 ov < expr->src.max ||
286 FLT_EQ(expr->src.max, ov)
287 )) || (
288 ov < expr->src.max
289 )) {
290 do_nv = 1;
291 break;
292 }
293 }
294 }
295 }
296 while (0);
297
298 /* no expression matched, do not continue */
299 if (!do_nv) continue;
300 RASTER_DEBUGF(3, "Using exprset[%d] unless NODATA", i);
301
302 /* converting a value from one range to another range
303 OldRange = (OldMax - OldMin)
304 NewRange = (NewMax - NewMin)
305 NewValue = (((OldValue - OldMin) * NewRange) / OldRange) + NewMin
306 */
307
308 /* NODATA */
309 if (hasnodata && isnodata) {
310 nv = nodataval;
311 }
312 /*
313 "src" min and max is the same, prevent division by zero
314 set nv to "dst" min, which should be the same as "dst" max
315 */
316 else if (FLT_EQ(expr->src.max, expr->src.min)) {
317 nv = expr->dst.min;
318 }
319 else {
320 or = expr->src.max - expr->src.min;
321 nr = expr->dst.max - expr->dst.min;
322 nv = (((ov - expr->src.min) * nr) / or) + expr->dst.min;
323
324 /* if dst range is from high to low */
325 if (expr->dst.min > expr->dst.max) {
326 if (nv > expr->dst.min)
327 nv = expr->dst.min;
328 else if (nv < expr->dst.max)
329 nv = expr->dst.max;
330 }
331 /* if dst range is from low to high */
332 else {
333 if (nv < expr->dst.min)
334 nv = expr->dst.min;
335 else if (nv > expr->dst.max)
336 nv = expr->dst.max;
337 }
338 }
339
340 /* round the value for integers */
341 switch (pixtype) {
342 case PT_1BB:
343 case PT_2BUI:
344 case PT_4BUI:
345 case PT_8BSI:
346 case PT_8BUI:
347 case PT_16BSI:
348 case PT_16BUI:
349 case PT_32BSI:
350 case PT_32BUI:
351 nv = round(nv);
352 break;
353 default:
354 break;
355 }
356
357 RASTER_DEBUGF(4, "(%d, %d) ov: %f or: %f - %f nr: %f - %f nv: %f"
358 , x
359 , y
360 , ov
361 , (NULL != expr) ? expr->src.min : 0
362 , (NULL != expr) ? expr->src.max : 0
363 , (NULL != expr) ? expr->dst.min : 0
364 , (NULL != expr) ? expr->dst.max : 0
365 , nv
366 );
367 if (rt_band_set_pixel(band, x, y, nv, NULL) != ES_NONE) {
368 rterror("rt_band_reclass: Could not assign value to new band");
369 rt_band_destroy(band);
370 rtdealloc(mem);
371 return 0;
372 }
373
374 expr = NULL;
375 }
376 }
377
378 return band;
379 }
380
381 /******************************************************************************
382 * rt_raster_iterator()
383 ******************************************************************************/
384
385 typedef struct _rti_iterator_arg_t* _rti_iterator_arg;
386 struct _rti_iterator_arg_t {
387 uint32_t count;
388
389 rt_raster *raster;
390 int *isempty;
391 double **offset;
392 int *width;
393 int *height;
394
395 struct {
396 rt_band *rtband;
397 int *hasnodata;
398 int *isnodata;
399 double *nodataval;
400 double *minval;
401 } band;
402
403 struct {
404 uint16_t x;
405 uint16_t y;
406 } distance;
407
408 struct {
409 uint32_t rows;
410 uint32_t columns;
411 } dimension;
412
413 struct {
414 double **values;
415 int **nodata;
416 } empty;
417
418 rt_iterator_arg arg;
419 };
420
421 static _rti_iterator_arg
_rti_iterator_arg_init()422 _rti_iterator_arg_init() {
423 _rti_iterator_arg _param;
424
425 _param = rtalloc(sizeof(struct _rti_iterator_arg_t));
426 if (_param == NULL) {
427 rterror("_rti_iterator_arg_init: Could not allocate memory for _rti_iterator_arg");
428 return NULL;
429 }
430
431 _param->count = 0;
432
433 _param->raster = NULL;
434 _param->isempty = NULL;
435 _param->offset = NULL;
436 _param->width = NULL;
437 _param->height = NULL;
438
439 _param->band.rtband = NULL;
440 _param->band.hasnodata = NULL;
441 _param->band.isnodata = NULL;
442 _param->band.nodataval = NULL;
443 _param->band.minval = NULL;
444
445 _param->distance.x = 0;
446 _param->distance.y = 0;
447
448 _param->dimension.rows = 0;
449 _param->dimension.columns = 0;
450
451 _param->empty.values = NULL;
452 _param->empty.nodata = NULL;
453
454 _param->arg = NULL;
455
456 return _param;
457 }
458
459 static void
_rti_iterator_arg_destroy(_rti_iterator_arg _param)460 _rti_iterator_arg_destroy(_rti_iterator_arg _param) {
461 uint32_t i = 0;
462
463 if (_param->raster != NULL)
464 rtdealloc(_param->raster);
465 if (_param->isempty != NULL)
466 rtdealloc(_param->isempty);
467 if (_param->width != NULL)
468 rtdealloc(_param->width);
469 if (_param->height != NULL)
470 rtdealloc(_param->height);
471
472 if (_param->band.rtband != NULL)
473 rtdealloc(_param->band.rtband);
474 if (_param->band.hasnodata != NULL)
475 rtdealloc(_param->band.hasnodata);
476 if (_param->band.isnodata != NULL)
477 rtdealloc(_param->band.isnodata);
478 if (_param->band.nodataval != NULL)
479 rtdealloc(_param->band.nodataval);
480 if (_param->band.minval != NULL)
481 rtdealloc(_param->band.minval);
482
483 if (_param->offset != NULL) {
484 for (i = 0; i < _param->count; i++) {
485 if (_param->offset[i] == NULL)
486 continue;
487 rtdealloc(_param->offset[i]);
488 }
489 rtdealloc(_param->offset);
490 }
491
492 if (_param->empty.values != NULL) {
493 for (i = 0; i < _param->dimension.rows; i++) {
494 if (_param->empty.values[i] == NULL)
495 continue;
496 rtdealloc(_param->empty.values[i]);
497 }
498 rtdealloc(_param->empty.values);
499 }
500 if (_param->empty.nodata != NULL) {
501 for (i = 0; i < _param->dimension.rows; i++) {
502 if (_param->empty.nodata[i] == NULL)
503 continue;
504 rtdealloc(_param->empty.nodata[i]);
505 }
506 rtdealloc(_param->empty.nodata);
507 }
508
509 if (_param->arg != NULL) {
510 if (_param->arg->values != NULL)
511 rtdealloc(_param->arg->values);
512 if (_param->arg->nodata != NULL)
513 rtdealloc(_param->arg->nodata);
514 if (_param->arg->src_pixel != NULL) {
515 for (i = 0; i < _param->count; i++) {
516 if (_param->arg->src_pixel[i] == NULL)
517 continue;
518 rtdealloc(_param->arg->src_pixel[i]);
519 }
520
521 rtdealloc(_param->arg->src_pixel);
522 }
523
524 rtdealloc(_param->arg);
525 }
526
527 rtdealloc(_param);
528 }
529
530 static int
_rti_iterator_arg_populate(_rti_iterator_arg _param,rt_iterator itrset,uint16_t itrcount,uint16_t distancex,uint16_t distancey,int * allnull,int * allempty)531 _rti_iterator_arg_populate(
532 _rti_iterator_arg _param,
533 rt_iterator itrset, uint16_t itrcount,
534 uint16_t distancex, uint16_t distancey,
535 int *allnull, int *allempty
536 ) {
537 int i = 0;
538 int hasband = 0;
539
540 _param->count = itrcount;
541 _param->distance.x = distancex;
542 _param->distance.y = distancey;
543 _param->dimension.columns = distancex * 2 + 1;
544 _param->dimension.rows = distancey * 2 + 1;
545
546 /* allocate memory for children */
547 _param->raster = rtalloc(sizeof(rt_raster) * itrcount);
548 _param->isempty = rtalloc(sizeof(int) * itrcount);
549 _param->width = rtalloc(sizeof(int) * itrcount);
550 _param->height = rtalloc(sizeof(int) * itrcount);
551
552 _param->offset = rtalloc(sizeof(double *) * itrcount);
553
554 _param->band.rtband = rtalloc(sizeof(rt_band) * itrcount);
555 _param->band.hasnodata = rtalloc(sizeof(int) * itrcount);
556 _param->band.isnodata = rtalloc(sizeof(int) * itrcount);
557 _param->band.nodataval = rtalloc(sizeof(double) * itrcount);
558 _param->band.minval = rtalloc(sizeof(double) * itrcount);
559
560 if (
561 _param->raster == NULL ||
562 _param->isempty == NULL ||
563 _param->width == NULL ||
564 _param->height == NULL ||
565 _param->offset == NULL ||
566 _param->band.rtband == NULL ||
567 _param->band.hasnodata == NULL ||
568 _param->band.isnodata == NULL ||
569 _param->band.nodataval == NULL ||
570 _param->band.minval == NULL
571 ) {
572 rterror("_rti_iterator_arg_populate: Could not allocate memory for children of _rti_iterator_arg");
573 return 0;
574 }
575
576 *allnull = 0;
577 *allempty = 0;
578
579 /*
580 check input rasters
581 not empty, band # is valid
582 copy raster pointers and set flags
583 */
584 for (i = 0; i < itrcount; i++) {
585 /* initialize elements */
586 _param->raster[i] = NULL;
587 _param->isempty[i] = 0;
588 _param->width[i] = 0;
589 _param->height[i] = 0;
590
591 _param->offset[i] = NULL;
592
593 _param->band.rtband[i] = NULL;
594 _param->band.hasnodata[i] = 0;
595 _param->band.isnodata[i] = 0;
596 _param->band.nodataval[i] = 0;
597 _param->band.minval[i] = 0;
598
599 /* set isempty */
600 if (itrset[i].raster == NULL) {
601 _param->isempty[i] = 1;
602
603 (*allnull)++;
604 (*allempty)++;
605
606 continue;
607 }
608 else if (rt_raster_is_empty(itrset[i].raster)) {
609 _param->isempty[i] = 1;
610
611 (*allempty)++;
612
613 continue;
614 }
615
616 /* check band number */
617 hasband = rt_raster_has_band(itrset[i].raster, itrset[i].nband);
618 if (!hasband) {
619 if (!itrset[i].nbnodata) {
620 rterror("_rti_iterator_arg_populate: Band %d not found for raster %d", itrset[i].nband, i);
621 return 0;
622 }
623 else {
624 RASTER_DEBUGF(4, "Band %d not found for raster %d. Using NODATA", itrset[i].nband, i);
625 }
626 }
627
628 _param->raster[i] = itrset[i].raster;
629 if (hasband) {
630 _param->band.rtband[i] = rt_raster_get_band(itrset[i].raster, itrset[i].nband);
631 if (_param->band.rtband[i] == NULL) {
632 rterror("_rti_iterator_arg_populate: Could not get band %d for raster %d", itrset[i].nband, i);
633 return 0;
634 }
635
636 /* hasnodata */
637 _param->band.hasnodata[i] = rt_band_get_hasnodata_flag(_param->band.rtband[i]);
638
639 /* hasnodata = TRUE */
640 if (_param->band.hasnodata[i]) {
641 /* nodataval */
642 rt_band_get_nodata(_param->band.rtband[i], &(_param->band.nodataval[i]));
643
644 /* isnodata */
645 _param->band.isnodata[i] = rt_band_get_isnodata_flag(_param->band.rtband[i]);
646 }
647 /* hasnodata = FALSE */
648 else {
649 /* minval */
650 _param->band.minval[i] = rt_band_get_min_value(_param->band.rtband[i]);
651 }
652 }
653
654 /* width, height */
655 _param->width[i] = rt_raster_get_width(_param->raster[i]);
656 _param->height[i] = rt_raster_get_height(_param->raster[i]);
657
658 /* init offset */
659 _param->offset[i] = rtalloc(sizeof(double) * 2);
660 if (_param->offset[i] == NULL) {
661 rterror("_rti_iterator_arg_populate: Could not allocate memory for offsets");
662 return 0;
663 }
664 }
665
666 return 1;
667 }
668
669 static int
_rti_iterator_arg_empty_init(_rti_iterator_arg _param)670 _rti_iterator_arg_empty_init(_rti_iterator_arg _param) {
671 uint32_t x = 0;
672 uint32_t y = 0;
673
674 _param->empty.values = rtalloc(sizeof(double *) * _param->dimension.rows);
675 _param->empty.nodata = rtalloc(sizeof(int *) * _param->dimension.rows);
676 if (_param->empty.values == NULL || _param->empty.nodata == NULL) {
677 rterror("_rti_iterator_arg_empty_init: Could not allocate memory for empty values and NODATA");
678 return 0;
679 }
680
681 for (y = 0; y < _param->dimension.rows; y++) {
682 _param->empty.values[y] = rtalloc(sizeof(double) * _param->dimension.columns);
683 _param->empty.nodata[y] = rtalloc(sizeof(int) * _param->dimension.columns);
684
685 if (_param->empty.values[y] == NULL || _param->empty.nodata[y] == NULL) {
686 rterror("_rti_iterator_arg_empty_init: Could not allocate memory for elements of empty values and NODATA");
687 return 0;
688 }
689
690 for (x = 0; x < _param->dimension.columns; x++) {
691 _param->empty.values[y][x] = 0;
692 _param->empty.nodata[y][x] = 1;
693 }
694 }
695
696 return 1;
697 }
698
699 static int
_rti_iterator_arg_callback_init(_rti_iterator_arg _param)700 _rti_iterator_arg_callback_init(_rti_iterator_arg _param) {
701 uint32_t i = 0;
702
703 _param->arg = rtalloc(sizeof(struct rt_iterator_arg_t));
704 if (_param->arg == NULL) {
705 rterror("_rti_iterator_arg_callback_init: Could not allocate memory for rt_iterator_arg");
706 return 0;
707 }
708
709 _param->arg->values = NULL;
710 _param->arg->nodata = NULL;
711 _param->arg->src_pixel = NULL;
712
713 /* initialize argument components */
714 _param->arg->values = rtalloc(sizeof(double **) * _param->count);
715 _param->arg->nodata = rtalloc(sizeof(int **) * _param->count);
716 _param->arg->src_pixel = rtalloc(sizeof(int *) * _param->count);
717 if (_param->arg->values == NULL || _param->arg->nodata == NULL || _param->arg->src_pixel == NULL) {
718 rterror("_rti_iterator_arg_callback_init: Could not allocate memory for element of rt_iterator_arg");
719 return 0;
720 }
721 memset(_param->arg->values, 0, sizeof(double **) * _param->count);
722 memset(_param->arg->nodata, 0, sizeof(int **) * _param->count);
723
724 /* initialize pos */
725 for (i = 0; i < _param->count; i++) {
726
727 _param->arg->src_pixel[i] = rtalloc(sizeof(int) * 2);
728 if (_param->arg->src_pixel[i] == NULL) {
729 rterror("_rti_iterator_arg_callback_init: Could not allocate memory for position elements of rt_iterator_arg");
730 return 0;
731 }
732 memset(_param->arg->src_pixel[i], 0, sizeof(int) * 2);
733 }
734
735 _param->arg->rasters = _param->count;
736 _param->arg->rows = _param->dimension.rows;
737 _param->arg->columns = _param->dimension.columns;
738
739 _param->arg->dst_pixel[0] = 0;
740 _param->arg->dst_pixel[1] = 0;
741
742 return 1;
743 }
744
745 static void
_rti_iterator_arg_callback_clean(_rti_iterator_arg _param)746 _rti_iterator_arg_callback_clean(_rti_iterator_arg _param) {
747 uint32_t i = 0;
748 uint32_t y = 0;
749
750 for (i = 0; i < _param->count; i++) {
751 RASTER_DEBUGF(5, "empty at @ %p", _param->empty.values);
752 RASTER_DEBUGF(5, "values at @ %p", _param->arg->values[i]);
753
754 if (_param->arg->values[i] == _param->empty.values) {
755 _param->arg->values[i] = NULL;
756 _param->arg->nodata[i] = NULL;
757
758 continue;
759 }
760
761 for (y = 0; y < _param->dimension.rows; y++) {
762 rtdealloc(_param->arg->values[i][y]);
763 rtdealloc(_param->arg->nodata[i][y]);
764 }
765
766 rtdealloc(_param->arg->values[i]);
767 rtdealloc(_param->arg->nodata[i]);
768
769 _param->arg->values[i] = NULL;
770 _param->arg->nodata[i] = NULL;
771 }
772 }
773
774 /**
775 * n-raster iterator.
776 * The raster returned should be freed by the caller
777 *
778 * @param itrset : set of rt_iterator objects.
779 * @param itrcount : number of objects in itrset.
780 * @param extenttype : type of extent for the output raster.
781 * @param customextent : raster specifying custom extent.
782 * is only used if extenttype is ET_CUSTOM.
783 * @param pixtype : the desired pixel type of the output raster's band.
784 * @param hasnodata : indicates if the band has nodata value
785 * @param nodataval : the nodata value, will be appropriately
786 * truncated to fit the pixtype size.
787 * @param distancex : the number of pixels around the specified pixel
788 * along the X axis
789 * @param distancey : the number of pixels around the specified pixel
790 * along the Y axis
791 * @param mask : the object of mask
792 * @param userarg : pointer to any argument that is passed as-is to callback.
793 * @param callback : callback function for actual processing of pixel values.
794 * @param *rtnraster : return one band raster from iterator process
795 *
796 * The callback function _must_ have the following signature.
797 *
798 * int FNAME(rt_iterator_arg arg, void *userarg, double *value, int *nodata)
799 *
800 * The callback function _must_ return zero (error) or non-zero (success)
801 * indicating whether the function ran successfully.
802 * The parameters passed to the callback function are as follows.
803 *
804 * - rt_iterator_arg arg: struct containing pixel values, NODATA flags and metadata
805 * - void *userarg: NULL or calling function provides to rt_raster_iterator() for use by callback function
806 * - double *value: value of pixel to be burned by rt_raster_iterator()
807 * - int *nodata: flag (0 or 1) indicating that pixel to be burned is NODATA
808 *
809 * @return ES_NONE on success, ES_ERROR on error
810 */
811 rt_errorstate
rt_raster_iterator(rt_iterator itrset,uint16_t itrcount,rt_extenttype extenttype,rt_raster customextent,rt_pixtype pixtype,uint8_t hasnodata,double nodataval,uint16_t distancex,uint16_t distancey,rt_mask mask,void * userarg,int (* callback)(rt_iterator_arg arg,void * userarg,double * value,int * nodata),rt_raster * rtnraster)812 rt_raster_iterator(
813 rt_iterator itrset, uint16_t itrcount,
814 rt_extenttype extenttype, rt_raster customextent,
815 rt_pixtype pixtype,
816 uint8_t hasnodata, double nodataval,
817 uint16_t distancex, uint16_t distancey,
818 rt_mask mask,
819 void *userarg,
820 int (*callback)(
821 rt_iterator_arg arg,
822 void *userarg,
823 double *value,
824 int *nodata
825 ),
826 rt_raster *rtnraster
827 ) {
828 /* output raster */
829 rt_raster rtnrast = NULL;
830 /* output raster's band */
831 rt_band rtnband = NULL;
832
833 /* working raster */
834 rt_raster rast = NULL;
835
836 _rti_iterator_arg _param = NULL;
837 int allnull = 0;
838 int allempty = 0;
839 int aligned = 0;
840 double offset[4] = {0.};
841 rt_pixel npixels;
842
843 int i = 0;
844 int status = 0;
845 int inextent = 0;
846 int x = 0;
847 int y = 0;
848 int _x = 0;
849 int _y = 0;
850
851 int _width = 0;
852 int _height = 0;
853
854 double minval;
855 double value;
856 int isnodata;
857 int nodata;
858
859 RASTER_DEBUG(3, "Starting...");
860
861 assert(itrset != NULL && itrcount > 0);
862 assert(rtnraster != NULL);
863
864 /* init rtnraster to NULL */
865 *rtnraster = NULL;
866
867 /* check that callback function is not NULL */
868 if (callback == NULL) {
869 rterror("rt_raster_iterator: Callback function not provided");
870 return ES_ERROR;
871 }
872
873 /* check that custom extent is provided if extenttype = ET_CUSTOM */
874 if (extenttype == ET_CUSTOM && rt_raster_is_empty(customextent)) {
875 rterror("rt_raster_iterator: Custom extent cannot be empty if extent type is ET_CUSTOM");
876 return ES_ERROR;
877 }
878
879 /* check that pixtype != PT_END */
880 if (pixtype == PT_END) {
881 rterror("rt_raster_iterator: Pixel type cannot be PT_END");
882 return ES_ERROR;
883 }
884
885 /* initialize _param */
886 if ((_param = _rti_iterator_arg_init()) == NULL) {
887 rterror("rt_raster_iterator: Could not initialize internal variables");
888 return ES_ERROR;
889 }
890
891 /* fill _param */
892 if (!_rti_iterator_arg_populate(_param, itrset, itrcount, distancex, distancey, &allnull, &allempty)) {
893 rterror("rt_raster_iterator: Could not populate for internal variables");
894 _rti_iterator_arg_destroy(_param);
895 return ES_ERROR;
896 }
897
898 /* shortcut if all null, return NULL */
899 if (allnull == itrcount) {
900 RASTER_DEBUG(3, "all rasters are NULL, returning NULL");
901
902 _rti_iterator_arg_destroy(_param);
903
904 return ES_NONE;
905 }
906 /* shortcut if all empty, return empty raster */
907 else if (allempty == itrcount) {
908 RASTER_DEBUG(3, "all rasters are empty, returning empty raster");
909
910 _rti_iterator_arg_destroy(_param);
911
912 rtnrast = rt_raster_new(0, 0);
913 if (rtnrast == NULL) {
914 rterror("rt_raster_iterator: Could not create empty raster");
915 return ES_ERROR;
916 }
917 rt_raster_set_scale(rtnrast, 0, 0);
918
919 *rtnraster = rtnrast;
920 return ES_NONE;
921 }
922
923 /* check that all rasters are aligned */
924 RASTER_DEBUG(3, "checking alignment of all rasters");
925 rast = NULL;
926
927 /* find raster to use as reference */
928 /* use custom if provided */
929 if (extenttype == ET_CUSTOM) {
930 RASTER_DEBUG(4, "using custom extent as reference raster");
931 rast = customextent;
932 }
933 /* use first valid one in _param->raster */
934 else {
935 for (i = 0; i < itrcount; i++) {
936 if (!_param->isempty[i]) {
937 RASTER_DEBUGF(4, "using raster at index %d as reference raster", i);
938 rast = _param->raster[i];
939 break;
940 }
941 }
942 }
943
944 /* no rasters found, SHOULD NEVER BE HERE! */
945 if (rast == NULL) {
946 rterror("rt_raster_iterator: Could not find reference raster to use for alignment tests");
947
948 _rti_iterator_arg_destroy(_param);
949
950 return ES_ERROR;
951 }
952
953 do {
954 aligned = 1;
955
956 /* check custom first if set. also skip if rasters are the same */
957 if (extenttype == ET_CUSTOM && rast != customextent) {
958 if (rt_raster_same_alignment(rast, customextent, &aligned, NULL) != ES_NONE) {
959 rterror("rt_raster_iterator: Could not test for alignment between reference raster and custom extent");
960
961 _rti_iterator_arg_destroy(_param);
962
963 return ES_ERROR;
964 }
965
966 RASTER_DEBUGF(5, "custom extent alignment: %d", aligned);
967 if (!aligned)
968 break;
969 }
970
971 for (i = 0; i < itrcount; i++) {
972 /* skip NULL rasters and if rasters are the same */
973 if (_param->isempty[i] || rast == _param->raster[i])
974 continue;
975
976 if (rt_raster_same_alignment(rast, _param->raster[i], &aligned, NULL) != ES_NONE) {
977 rterror("rt_raster_iterator: Could not test for alignment between reference raster and raster %d", i);
978
979 _rti_iterator_arg_destroy(_param);
980
981 return ES_ERROR;
982 }
983 RASTER_DEBUGF(5, "raster at index %d alignment: %d", i, aligned);
984
985 /* abort checking since a raster isn't aligned */
986 if (!aligned)
987 break;
988 }
989 }
990 while (0);
991
992 /* not aligned, error */
993 if (!aligned) {
994 rterror("rt_raster_iterator: The set of rasters provided (custom extent included, if appropriate) do not have the same alignment");
995
996 _rti_iterator_arg_destroy(_param);
997
998 return ES_ERROR;
999 }
1000
1001 /* use extenttype to build output raster (no bands though) */
1002 i = -1;
1003 switch (extenttype) {
1004 case ET_INTERSECTION:
1005 case ET_UNION:
1006 /* make copy of first "real" raster */
1007 rtnrast = rtalloc(sizeof(struct rt_raster_t));
1008 if (rtnrast == NULL) {
1009 rterror("rt_raster_iterator: Could not allocate memory for output raster");
1010
1011 _rti_iterator_arg_destroy(_param);
1012
1013 return ES_ERROR;
1014 }
1015
1016 for (i = 0; i < itrcount; i++) {
1017 if (!_param->isempty[i]) {
1018 memcpy(rtnrast, _param->raster[i], sizeof(struct rt_raster_t));
1019 break;
1020 }
1021 }
1022 rtnrast->numBands = 0;
1023 rtnrast->bands = NULL;
1024
1025 /* get extent of output raster */
1026 rast = NULL;
1027 for (i = i + 1; i < itrcount; i++) {
1028 if (_param->isempty[i])
1029 continue;
1030
1031 status = rt_raster_from_two_rasters(rtnrast, _param->raster[i], extenttype, &rast, NULL);
1032 rtdealloc(rtnrast);
1033
1034 if (rast == NULL || status != ES_NONE) {
1035 rterror("rt_raster_iterator: Could not compute %s extent of rasters",
1036 extenttype == ET_UNION ? "union" : "intersection"
1037 );
1038
1039 _rti_iterator_arg_destroy(_param);
1040
1041 return ES_ERROR;
1042 }
1043 else if (rt_raster_is_empty(rast)) {
1044 rtinfo("rt_raster_iterator: Computed raster for %s extent is empty",
1045 extenttype == ET_UNION ? "union" : "intersection"
1046 );
1047
1048 _rti_iterator_arg_destroy(_param);
1049
1050 *rtnraster = rast;
1051 return ES_NONE;
1052 }
1053
1054 rtnrast = rast;
1055 rast = NULL;
1056 }
1057
1058 break;
1059 /*
1060 first, second and last have similar checks
1061 and continue into custom
1062 */
1063 case ET_FIRST:
1064 i = 0;
1065 /* FALLTHROUGH */
1066 case ET_SECOND:
1067 if (i < 0) {
1068 if (itrcount < 2)
1069 i = 0;
1070 else
1071 i = 1;
1072 }
1073 /* FALLTHROUGH */
1074 case ET_LAST:
1075 if (i < 0) i = itrcount - 1;
1076
1077 /* input raster is null, return NULL */
1078 if (_param->raster[i] == NULL) {
1079 RASTER_DEBUGF(3, "returning NULL as %s raster is NULL and extent type is ET_%s",
1080 (i == 0 ? "first" : (i == 1 ? "second" : "last")),
1081 (i == 0 ? "FIRST" : (i == 1 ? "SECOND" : "LAST"))
1082 );
1083
1084 _rti_iterator_arg_destroy(_param);
1085
1086 return ES_NONE;
1087 }
1088 /* input raster is empty, return empty raster */
1089 else if (_param->isempty[i]) {
1090 RASTER_DEBUGF(3, "returning empty raster as %s raster is empty and extent type is ET_%s",
1091 (i == 0 ? "first" : (i == 1 ? "second" : "last")),
1092 (i == 0 ? "FIRST" : (i == 1 ? "SECOND" : "LAST"))
1093 );
1094
1095 _rti_iterator_arg_destroy(_param);
1096
1097 rtnrast = rt_raster_new(0, 0);
1098 if (rtnrast == NULL) {
1099 rterror("rt_raster_iterator: Could not create empty raster");
1100 return ES_ERROR;
1101 }
1102 rt_raster_set_scale(rtnrast, 0, 0);
1103
1104 *rtnraster = rtnrast;
1105 return ES_NONE;
1106 }
1107 /* FALLTHROUGH */
1108 /* copy the custom extent raster */
1109 case ET_CUSTOM:
1110 rtnrast = rtalloc(sizeof(struct rt_raster_t));
1111 if (rtnrast == NULL) {
1112 rterror("rt_raster_iterator: Could not allocate memory for output raster");
1113
1114 _rti_iterator_arg_destroy(_param);
1115
1116 return ES_ERROR;
1117 }
1118
1119 switch (extenttype) {
1120 case ET_CUSTOM:
1121 memcpy(rtnrast, customextent, sizeof(struct rt_raster_t));
1122 break;
1123 /* first, second, last */
1124 default:
1125 memcpy(rtnrast, _param->raster[i], sizeof(struct rt_raster_t));
1126 break;
1127 }
1128 rtnrast->numBands = 0;
1129 rtnrast->bands = NULL;
1130 break;
1131 }
1132
1133 _width = rt_raster_get_width(rtnrast);
1134 _height = rt_raster_get_height(rtnrast);
1135
1136 RASTER_DEBUGF(4, "rtnrast (width, height, ulx, uly, scalex, scaley, skewx, skewy, srid) = (%d, %d, %f, %f, %f, %f, %f, %f, %d)",
1137 _width,
1138 _height,
1139 rt_raster_get_x_offset(rtnrast),
1140 rt_raster_get_y_offset(rtnrast),
1141 rt_raster_get_x_scale(rtnrast),
1142 rt_raster_get_y_scale(rtnrast),
1143 rt_raster_get_x_skew(rtnrast),
1144 rt_raster_get_y_skew(rtnrast),
1145 rt_raster_get_srid(rtnrast)
1146 );
1147
1148 /* init values and NODATA for use with empty rasters */
1149 if (!_rti_iterator_arg_empty_init(_param)) {
1150 rterror("rt_raster_iterator: Could not initialize empty values and NODATA");
1151
1152 _rti_iterator_arg_destroy(_param);
1153 rt_raster_destroy(rtnrast);
1154
1155 return ES_ERROR;
1156 }
1157
1158 /* create output band */
1159 if (rt_raster_generate_new_band(
1160 rtnrast,
1161 pixtype,
1162 nodataval,
1163 hasnodata, nodataval,
1164 0
1165 ) < 0) {
1166 rterror("rt_raster_iterator: Could not add new band to output raster");
1167
1168 _rti_iterator_arg_destroy(_param);
1169 rt_raster_destroy(rtnrast);
1170
1171 return ES_ERROR;
1172 }
1173
1174 /* get output band */
1175 rtnband = rt_raster_get_band(rtnrast, 0);
1176 if (rtnband == NULL) {
1177 rterror("rt_raster_iterator: Could not get new band from output raster");
1178
1179 _rti_iterator_arg_destroy(_param);
1180 rt_raster_destroy(rtnrast);
1181
1182 return ES_ERROR;
1183 }
1184
1185 /* output band's minimum value */
1186 minval = rt_band_get_min_value(rtnband);
1187
1188 /* initialize argument for callback function */
1189 if (!_rti_iterator_arg_callback_init(_param)) {
1190 rterror("rt_raster_iterator: Could not initialize callback function argument");
1191
1192 _rti_iterator_arg_destroy(_param);
1193 rt_band_destroy(rtnband);
1194 rt_raster_destroy(rtnrast);
1195
1196 return ES_ERROR;
1197 }
1198
1199 /* fill _param->offset */
1200 for (i = 0; i < itrcount; i++) {
1201 if (_param->isempty[i])
1202 continue;
1203
1204 status = rt_raster_from_two_rasters(rtnrast, _param->raster[i], ET_FIRST, &rast, offset);
1205 rtdealloc(rast);
1206 if (status != ES_NONE) {
1207 rterror("rt_raster_iterator: Could not compute raster offsets");
1208
1209 _rti_iterator_arg_destroy(_param);
1210 rt_band_destroy(rtnband);
1211 rt_raster_destroy(rtnrast);
1212
1213 return ES_ERROR;
1214 }
1215
1216 _param->offset[i][0] = offset[2];
1217 _param->offset[i][1] = offset[3];
1218 RASTER_DEBUGF(4, "rast %d offset: %f %f", i, offset[2], offset[3]);
1219 }
1220
1221 /* loop over each pixel (POI) of output raster */
1222 /* _x,_y are for output raster */
1223 /* x,y are for input raster */
1224 for (_y = 0; _y < _height; _y++) {
1225 for (_x = 0; _x < _width; _x++) {
1226 RASTER_DEBUGF(4, "iterating output pixel (x, y) = (%d, %d)", _x, _y);
1227 _param->arg->dst_pixel[0] = _x;
1228 _param->arg->dst_pixel[1] = _y;
1229
1230 /* loop through each input raster */
1231 for (i = 0; i < itrcount; i++) {
1232 RASTER_DEBUGF(4, "raster %d", i);
1233
1234 /*
1235 empty raster
1236 OR band does not exist and flag set to use NODATA
1237 OR band is NODATA
1238 */
1239 if (
1240 _param->isempty[i] ||
1241 (_param->band.rtband[i] == NULL && itrset[i].nbnodata) ||
1242 _param->band.isnodata[i]
1243 ) {
1244 RASTER_DEBUG(4, "empty raster, band does not exist or band is NODATA. using empty values and NODATA");
1245
1246 x = _x;
1247 y = _y;
1248
1249 _param->arg->values[i] = _param->empty.values;
1250 _param->arg->nodata[i] = _param->empty.nodata;
1251
1252 continue;
1253 }
1254
1255 /* input raster's X,Y */
1256 x = _x - (int) _param->offset[i][0];
1257 y = _y - (int) _param->offset[i][1];
1258 RASTER_DEBUGF(4, "source pixel (x, y) = (%d, %d)", x, y);
1259
1260 _param->arg->src_pixel[i][0] = x;
1261 _param->arg->src_pixel[i][1] = y;
1262
1263 /* neighborhood */
1264 npixels = NULL;
1265 status = 0;
1266 if (distancex > 0 && distancey > 0) {
1267 RASTER_DEBUG(4, "getting neighborhood");
1268
1269 status = rt_band_get_nearest_pixel(
1270 _param->band.rtband[i],
1271 x, y,
1272 distancex, distancey,
1273 1,
1274 &npixels
1275 );
1276 if (status < 0) {
1277 rterror("rt_raster_iterator: Could not get pixel neighborhood");
1278
1279 _rti_iterator_arg_destroy(_param);
1280 rt_band_destroy(rtnband);
1281 rt_raster_destroy(rtnrast);
1282
1283 return ES_ERROR;
1284 }
1285 }
1286
1287 /* get value of POI */
1288 /* get pixel's value */
1289 if (
1290 (x >= 0 && x < _param->width[i]) &&
1291 (y >= 0 && y < _param->height[i])
1292 ) {
1293 RASTER_DEBUG(4, "getting value of POI");
1294 if (rt_band_get_pixel(
1295 _param->band.rtband[i],
1296 x, y,
1297 &value,
1298 &isnodata
1299 ) != ES_NONE) {
1300 rterror("rt_raster_iterator: Could not get the pixel value of band");
1301
1302 _rti_iterator_arg_destroy(_param);
1303 rt_band_destroy(rtnband);
1304 rt_raster_destroy(rtnrast);
1305
1306 return ES_ERROR;
1307 }
1308 inextent = 1;
1309 }
1310 /* outside band extent, set to NODATA */
1311 else {
1312 RASTER_DEBUG(4, "Outside band extent, setting value to NODATA");
1313 /* has NODATA, use NODATA */
1314 if (_param->band.hasnodata[i])
1315 value = _param->band.nodataval[i];
1316 /* no NODATA, use min possible value */
1317 else
1318 value = _param->band.minval[i];
1319
1320 inextent = 0;
1321 isnodata = 1;
1322 }
1323
1324 /* add pixel to neighborhood */
1325 status++;
1326 if (status > 1)
1327 npixels = (rt_pixel) rtrealloc(npixels, sizeof(struct rt_pixel_t) * status);
1328 else
1329 npixels = (rt_pixel) rtalloc(sizeof(struct rt_pixel_t));
1330
1331 if (npixels == NULL) {
1332 rterror("rt_raster_iterator: Could not reallocate memory for neighborhood");
1333
1334 _rti_iterator_arg_destroy(_param);
1335 rt_band_destroy(rtnband);
1336 rt_raster_destroy(rtnrast);
1337
1338 return ES_ERROR;
1339 }
1340
1341 npixels[status - 1].x = x;
1342 npixels[status - 1].y = y;
1343 npixels[status - 1].nodata = 1;
1344 npixels[status - 1].value = value;
1345
1346 /* set nodata flag */
1347 if ((!_param->band.hasnodata[i] && inextent) || !isnodata) {
1348 npixels[status - 1].nodata = 0;
1349 }
1350 RASTER_DEBUGF(4, "value, nodata: %f, %d", value, npixels[status - 1].nodata);
1351
1352 /* convert set of rt_pixel to 2D array */
1353 status = rt_pixel_set_to_array(
1354 npixels, status, mask,
1355 x, y,
1356 distancex, distancey,
1357 &(_param->arg->values[i]),
1358 &(_param->arg->nodata[i]),
1359 NULL, NULL
1360 );
1361 rtdealloc(npixels);
1362 if (status != ES_NONE) {
1363 rterror("rt_raster_iterator: Could not create 2D array of neighborhood");
1364
1365 _rti_iterator_arg_destroy(_param);
1366 rt_band_destroy(rtnband);
1367 rt_raster_destroy(rtnrast);
1368
1369 return ES_ERROR;
1370 }
1371 }
1372
1373 /* callback */
1374 RASTER_DEBUG(4, "calling callback function");
1375 value = 0;
1376 nodata = 0;
1377 status = callback(_param->arg, userarg, &value, &nodata);
1378
1379 /* free memory from callback */
1380 _rti_iterator_arg_callback_clean(_param);
1381
1382 /* handle callback status */
1383 if (status == 0) {
1384 rterror("rt_raster_iterator: Callback function returned an error");
1385
1386 _rti_iterator_arg_destroy(_param);
1387 rt_band_destroy(rtnband);
1388 rt_raster_destroy(rtnrast);
1389
1390 return ES_ERROR;
1391 }
1392
1393 /* burn value to pixel */
1394 status = 0;
1395 if (!nodata) {
1396 status = rt_band_set_pixel(rtnband, _x, _y, value, NULL);
1397 RASTER_DEBUGF(4, "burning pixel (%d, %d) with value: %f", _x, _y, value);
1398 }
1399 else if (!hasnodata) {
1400 status = rt_band_set_pixel(rtnband, _x, _y, minval, NULL);
1401 RASTER_DEBUGF(4, "burning pixel (%d, %d) with minval: %f", _x, _y, minval);
1402 }
1403 else {
1404 RASTER_DEBUGF(4, "NOT burning pixel (%d, %d)", _x, _y);
1405 }
1406 if (status != ES_NONE) {
1407 rterror("rt_raster_iterator: Could not set pixel value");
1408
1409 _rti_iterator_arg_destroy(_param);
1410 rt_band_destroy(rtnband);
1411 rt_raster_destroy(rtnrast);
1412
1413 return ES_ERROR;
1414 }
1415 }
1416 }
1417
1418 /* lots of cleanup */
1419 _rti_iterator_arg_destroy(_param);
1420
1421 *rtnraster = rtnrast;
1422 return ES_NONE;
1423 }
1424
1425 /******************************************************************************
1426 * rt_raster_colormap()
1427 ******************************************************************************/
1428
1429 typedef struct _rti_colormap_arg_t* _rti_colormap_arg;
1430 struct _rti_colormap_arg_t {
1431 rt_raster raster;
1432 rt_band band;
1433
1434 rt_colormap_entry nodataentry;
1435 int hasnodata;
1436 double nodataval;
1437
1438 int nexpr;
1439 rt_reclassexpr *expr;
1440
1441 int npos;
1442 uint16_t *pos;
1443
1444 };
1445
1446 static _rti_colormap_arg
_rti_colormap_arg_init(rt_raster raster)1447 _rti_colormap_arg_init(rt_raster raster) {
1448 _rti_colormap_arg arg = NULL;
1449
1450 arg = rtalloc(sizeof(struct _rti_colormap_arg_t));
1451 if (arg == NULL) {
1452 rterror("_rti_colormap_arg_init: Could not allocate memory for _rti_color_arg");
1453 return NULL;
1454 }
1455
1456 arg->band = NULL;
1457 arg->nodataentry = NULL;
1458 arg->hasnodata = 0;
1459 arg->nodataval = 0;
1460
1461 if (raster == NULL)
1462 arg->raster = NULL;
1463 /* raster provided */
1464 else {
1465 arg->raster = rt_raster_clone(raster, 0);
1466 if (arg->raster == NULL) {
1467 rterror("_rti_colormap_arg_init: Could not create output raster");
1468 return NULL;
1469 }
1470 }
1471
1472 arg->nexpr = 0;
1473 arg->expr = NULL;
1474
1475 arg->npos = 0;
1476 arg->pos = NULL;
1477
1478 return arg;
1479 }
1480
1481 static void
_rti_colormap_arg_destroy(_rti_colormap_arg arg)1482 _rti_colormap_arg_destroy(_rti_colormap_arg arg) {
1483 int i = 0;
1484
1485 if (arg->raster != NULL) {
1486 rt_band band = NULL;
1487
1488 for (i = rt_raster_get_num_bands(arg->raster) - 1; i >= 0; i--) {
1489 band = rt_raster_get_band(arg->raster, i);
1490 if (band != NULL)
1491 rt_band_destroy(band);
1492 }
1493
1494 rt_raster_destroy(arg->raster);
1495 }
1496
1497 if (arg->nexpr) {
1498 for (i = 0; i < arg->nexpr; i++) {
1499 if (arg->expr[i] != NULL)
1500 rtdealloc(arg->expr[i]);
1501 }
1502 rtdealloc(arg->expr);
1503 }
1504
1505 if (arg->npos)
1506 rtdealloc(arg->pos);
1507
1508 rtdealloc(arg);
1509 arg = NULL;
1510 }
1511
1512 /**
1513 * Returns a new raster with up to four 8BUI bands (RGBA) from
1514 * applying a colormap to the user-specified band of the
1515 * input raster.
1516 *
1517 * @param raster: input raster
1518 * @param nband: 0-based index of the band to process with colormap
1519 * @param colormap: rt_colormap object of colormap to apply to band
1520 *
1521 * @return new raster or NULL on error
1522 */
rt_raster_colormap(rt_raster raster,int nband,rt_colormap colormap)1523 rt_raster rt_raster_colormap(
1524 rt_raster raster, int nband,
1525 rt_colormap colormap
1526 ) {
1527 _rti_colormap_arg arg = NULL;
1528 rt_raster rtnraster = NULL;
1529 rt_band band = NULL;
1530 int i = 0;
1531 int j = 0;
1532 int k = 0;
1533
1534 assert(colormap != NULL);
1535
1536 /* empty raster */
1537 if (rt_raster_is_empty(raster))
1538 return NULL;
1539
1540 /* no colormap entries */
1541 if (colormap->nentry < 1) {
1542 rterror("rt_raster_colormap: colormap must have at least one entry");
1543 return NULL;
1544 }
1545
1546 /* nband is valid */
1547 if (!rt_raster_has_band(raster, nband)) {
1548 rterror("rt_raster_colormap: raster has no band at index %d", nband);
1549 return NULL;
1550 }
1551
1552 band = rt_raster_get_band(raster, nband);
1553 if (band == NULL) {
1554 rterror("rt_raster_colormap: Could not get band at index %d", nband);
1555 return NULL;
1556 }
1557
1558 /* init internal variables */
1559 arg = _rti_colormap_arg_init(raster);
1560 if (arg == NULL) {
1561 rterror("rt_raster_colormap: Could not initialize internal variables");
1562 return NULL;
1563 }
1564
1565 /* handle NODATA */
1566 if (rt_band_get_hasnodata_flag(band)) {
1567 arg->hasnodata = 1;
1568 rt_band_get_nodata(band, &(arg->nodataval));
1569 }
1570
1571 /* # of colors */
1572 if (colormap->ncolor < 1) {
1573 rterror("rt_raster_colormap: At least one color must be provided");
1574 _rti_colormap_arg_destroy(arg);
1575 return NULL;
1576 }
1577 else if (colormap->ncolor > 4) {
1578 rtinfo("More than four colors indicated. Using only the first four colors");
1579 colormap->ncolor = 4;
1580 }
1581
1582 /* find non-NODATA entries */
1583 arg->npos = 0;
1584 arg->pos = rtalloc(sizeof(uint16_t) * colormap->nentry);
1585 if (arg->pos == NULL) {
1586 rterror("rt_raster_colormap: Could not allocate memory for valid entries");
1587 _rti_colormap_arg_destroy(arg);
1588 return NULL;
1589 }
1590 for (i = 0, j = 0; i < colormap->nentry; i++) {
1591 /* special handling for NODATA entries */
1592 if (colormap->entry[i].isnodata) {
1593 /* first NODATA entry found, use it */
1594 if (arg->nodataentry == NULL)
1595 arg->nodataentry = &(colormap->entry[i]);
1596 else
1597 rtwarn("More than one colormap entry found for NODATA value. Only using first NOTDATA entry");
1598
1599 continue;
1600 }
1601
1602 (arg->npos)++;
1603 arg->pos[j++] = i;
1604 }
1605
1606 /* INTERPOLATE and only one non-NODATA entry */
1607 if (colormap->method == CM_INTERPOLATE && arg->npos < 2) {
1608 rtwarn("Method INTERPOLATE requires at least two non-NODATA colormap entries. Using NEAREST instead");
1609 colormap->method = CM_NEAREST;
1610 }
1611
1612 /* NODATA entry but band has no NODATA value */
1613 if (!arg->hasnodata && arg->nodataentry != NULL) {
1614 rtinfo("Band at index %d has no NODATA value. Ignoring NODATA entry", nband);
1615 arg->nodataentry = NULL;
1616 }
1617
1618 /* allocate expr */
1619 arg->nexpr = arg->npos;
1620
1621 /* INTERPOLATE needs one less than the number of entries */
1622 if (colormap->method == CM_INTERPOLATE)
1623 arg->nexpr -= 1;
1624 /* EXACT requires a no matching expression */
1625 else if (colormap->method == CM_EXACT)
1626 arg->nexpr += 1;
1627
1628 /* NODATA entry exists, add expression */
1629 if (arg->nodataentry != NULL)
1630 arg->nexpr += 1;
1631 arg->expr = rtalloc(sizeof(rt_reclassexpr) * arg->nexpr);
1632 if (arg->expr == NULL) {
1633 rterror("rt_raster_colormap: Could not allocate memory for reclass expressions");
1634 _rti_colormap_arg_destroy(arg);
1635 return NULL;
1636 }
1637 RASTER_DEBUGF(4, "nexpr = %d", arg->nexpr);
1638 RASTER_DEBUGF(4, "expr @ %p", arg->expr);
1639
1640 for (i = 0; i < arg->nexpr; i++) {
1641 arg->expr[i] = rtalloc(sizeof(struct rt_reclassexpr_t));
1642 if (arg->expr[i] == NULL) {
1643 rterror("rt_raster_colormap: Could not allocate memory for reclass expression");
1644 _rti_colormap_arg_destroy(arg);
1645 return NULL;
1646 }
1647 }
1648
1649 /* reclassify bands */
1650 /* by # of colors */
1651 for (i = 0; i < colormap->ncolor; i++) {
1652 k = 0;
1653
1654 /* handle NODATA entry first */
1655 if (arg->nodataentry != NULL) {
1656 arg->expr[k]->src.min = arg->nodataentry->value;
1657 arg->expr[k]->src.max = arg->nodataentry->value;
1658 arg->expr[k]->src.inc_min = 1;
1659 arg->expr[k]->src.inc_max = 1;
1660 arg->expr[k]->src.exc_min = 0;
1661 arg->expr[k]->src.exc_max = 0;
1662
1663 arg->expr[k]->dst.min = arg->nodataentry->color[i];
1664 arg->expr[k]->dst.max = arg->nodataentry->color[i];
1665
1666 arg->expr[k]->dst.inc_min = 1;
1667 arg->expr[k]->dst.inc_max = 1;
1668 arg->expr[k]->dst.exc_min = 0;
1669 arg->expr[k]->dst.exc_max = 0;
1670
1671 RASTER_DEBUGF(4, "NODATA expr[%d]->src (min, max, in, ix, en, ex) = (%f, %f, %d, %d, %d, %d)",
1672 k,
1673 arg->expr[k]->src.min,
1674 arg->expr[k]->src.max,
1675 arg->expr[k]->src.inc_min,
1676 arg->expr[k]->src.inc_max,
1677 arg->expr[k]->src.exc_min,
1678 arg->expr[k]->src.exc_max
1679 );
1680 RASTER_DEBUGF(4, "NODATA expr[%d]->dst (min, max, in, ix, en, ex) = (%f, %f, %d, %d, %d, %d)",
1681 k,
1682 arg->expr[k]->dst.min,
1683 arg->expr[k]->dst.max,
1684 arg->expr[k]->dst.inc_min,
1685 arg->expr[k]->dst.inc_max,
1686 arg->expr[k]->dst.exc_min,
1687 arg->expr[k]->dst.exc_max
1688 );
1689
1690 k++;
1691 }
1692
1693 /* by non-NODATA entry */
1694 for (j = 0; j < arg->npos; j++) {
1695 if (colormap->method == CM_INTERPOLATE) {
1696 if (j == arg->npos - 1)
1697 continue;
1698
1699 arg->expr[k]->src.min = colormap->entry[arg->pos[j + 1]].value;
1700 arg->expr[k]->src.inc_min = 1;
1701 arg->expr[k]->src.exc_min = 0;
1702
1703 arg->expr[k]->src.max = colormap->entry[arg->pos[j]].value;
1704 arg->expr[k]->src.inc_max = 1;
1705 arg->expr[k]->src.exc_max = 0;
1706
1707 arg->expr[k]->dst.min = colormap->entry[arg->pos[j + 1]].color[i];
1708 arg->expr[k]->dst.max = colormap->entry[arg->pos[j]].color[i];
1709
1710 arg->expr[k]->dst.inc_min = 1;
1711 arg->expr[k]->dst.exc_min = 0;
1712
1713 arg->expr[k]->dst.inc_max = 1;
1714 arg->expr[k]->dst.exc_max = 0;
1715 }
1716 else if (colormap->method == CM_NEAREST) {
1717
1718 /* NOT last entry */
1719 if (j != arg->npos - 1) {
1720 arg->expr[k]->src.min = ((colormap->entry[arg->pos[j]].value - colormap->entry[arg->pos[j + 1]].value) / 2.) + colormap->entry[arg->pos[j + 1]].value;
1721 arg->expr[k]->src.inc_min = 0;
1722 arg->expr[k]->src.exc_min = 0;
1723 }
1724 /* last entry */
1725 else {
1726 arg->expr[k]->src.min = colormap->entry[arg->pos[j]].value;
1727 arg->expr[k]->src.inc_min = 1;
1728 arg->expr[k]->src.exc_min = 1;
1729 }
1730
1731 /* NOT first entry */
1732 if (j > 0) {
1733 arg->expr[k]->src.max = arg->expr[k - 1]->src.min;
1734 arg->expr[k]->src.inc_max = 1;
1735 arg->expr[k]->src.exc_max = 0;
1736 }
1737 /* first entry */
1738 else {
1739 arg->expr[k]->src.max = colormap->entry[arg->pos[j]].value;
1740 arg->expr[k]->src.inc_max = 1;
1741 arg->expr[k]->src.exc_max = 1;
1742 }
1743
1744 arg->expr[k]->dst.min = colormap->entry[arg->pos[j]].color[i];
1745 arg->expr[k]->dst.inc_min = 1;
1746 arg->expr[k]->dst.exc_min = 0;
1747
1748 arg->expr[k]->dst.max = colormap->entry[arg->pos[j]].color[i];
1749 arg->expr[k]->dst.inc_max = 1;
1750 arg->expr[k]->dst.exc_max = 0;
1751 }
1752 else if (colormap->method == CM_EXACT) {
1753 arg->expr[k]->src.min = colormap->entry[arg->pos[j]].value;
1754 arg->expr[k]->src.inc_min = 1;
1755 arg->expr[k]->src.exc_min = 0;
1756
1757 arg->expr[k]->src.max = colormap->entry[arg->pos[j]].value;
1758 arg->expr[k]->src.inc_max = 1;
1759 arg->expr[k]->src.exc_max = 0;
1760
1761 arg->expr[k]->dst.min = colormap->entry[arg->pos[j]].color[i];
1762 arg->expr[k]->dst.inc_min = 1;
1763 arg->expr[k]->dst.exc_min = 0;
1764
1765 arg->expr[k]->dst.max = colormap->entry[arg->pos[j]].color[i];
1766 arg->expr[k]->dst.inc_max = 1;
1767 arg->expr[k]->dst.exc_max = 0;
1768 }
1769
1770 RASTER_DEBUGF(4, "expr[%d]->src (min, max, in, ix, en, ex) = (%f, %f, %d, %d, %d, %d)",
1771 k,
1772 arg->expr[k]->src.min,
1773 arg->expr[k]->src.max,
1774 arg->expr[k]->src.inc_min,
1775 arg->expr[k]->src.inc_max,
1776 arg->expr[k]->src.exc_min,
1777 arg->expr[k]->src.exc_max
1778 );
1779
1780 RASTER_DEBUGF(4, "expr[%d]->dst (min, max, in, ix, en, ex) = (%f, %f, %d, %d, %d, %d)",
1781 k,
1782 arg->expr[k]->dst.min,
1783 arg->expr[k]->dst.max,
1784 arg->expr[k]->dst.inc_min,
1785 arg->expr[k]->dst.inc_max,
1786 arg->expr[k]->dst.exc_min,
1787 arg->expr[k]->dst.exc_max
1788 );
1789
1790 k++;
1791 }
1792
1793 /* EXACT has one last expression for catching all uncaught values */
1794 if (colormap->method == CM_EXACT) {
1795 arg->expr[k]->src.min = 0;
1796 arg->expr[k]->src.inc_min = 1;
1797 arg->expr[k]->src.exc_min = 1;
1798
1799 arg->expr[k]->src.max = 0;
1800 arg->expr[k]->src.inc_max = 1;
1801 arg->expr[k]->src.exc_max = 1;
1802
1803 arg->expr[k]->dst.min = 0;
1804 arg->expr[k]->dst.inc_min = 1;
1805 arg->expr[k]->dst.exc_min = 0;
1806
1807 arg->expr[k]->dst.max = 0;
1808 arg->expr[k]->dst.inc_max = 1;
1809 arg->expr[k]->dst.exc_max = 0;
1810
1811 RASTER_DEBUGF(4, "expr[%d]->src (min, max, in, ix, en, ex) = (%f, %f, %d, %d, %d, %d)",
1812 k,
1813 arg->expr[k]->src.min,
1814 arg->expr[k]->src.max,
1815 arg->expr[k]->src.inc_min,
1816 arg->expr[k]->src.inc_max,
1817 arg->expr[k]->src.exc_min,
1818 arg->expr[k]->src.exc_max
1819 );
1820
1821 RASTER_DEBUGF(4, "expr[%d]->dst (min, max, in, ix, en, ex) = (%f, %f, %d, %d, %d, %d)",
1822 k,
1823 arg->expr[k]->dst.min,
1824 arg->expr[k]->dst.max,
1825 arg->expr[k]->dst.inc_min,
1826 arg->expr[k]->dst.inc_max,
1827 arg->expr[k]->dst.exc_min,
1828 arg->expr[k]->dst.exc_max
1829 );
1830
1831 k++;
1832 }
1833
1834 /* call rt_band_reclass */
1835 arg->band = rt_band_reclass(band, PT_8BUI, 0, 0, arg->expr, arg->nexpr);
1836 if (arg->band == NULL) {
1837 rterror("rt_raster_colormap: Could not reclassify band");
1838 _rti_colormap_arg_destroy(arg);
1839 return NULL;
1840 }
1841
1842 /* add reclassified band to raster */
1843 if (rt_raster_add_band(arg->raster, arg->band, rt_raster_get_num_bands(arg->raster)) < 0) {
1844 rterror("rt_raster_colormap: Could not add reclassified band to output raster");
1845 _rti_colormap_arg_destroy(arg);
1846 return NULL;
1847 }
1848 }
1849
1850 rtnraster = arg->raster;
1851 arg->raster = NULL;
1852 _rti_colormap_arg_destroy(arg);
1853
1854 return rtnraster;
1855 }
1856