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 * Copyright (C) 2013 Nathaniel Hunter Clay <clay.nathaniel@gmail.com>
14 *
15 * This program is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU General Public License
17 * as published by the Free Software Foundation; either version 2
18 * of the License, or (at your option) any later version.
19 *
20 * This program is distributed in the hope that it will be useful,
21 * but WITHOUT ANY WARRANTY; without even the implied warranty of
22 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
23 * GNU General Public License for more details.
24 *
25 * You should have received a copy of the GNU General Public License
26 * along with this program; if not, write to the Free Software Foundation,
27 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
28 *
29 */
30
31 #include "librtcore.h"
32 #include "librtcore_internal.h"
33
34 /******************************************************************************
35 * rt_pixeltype
36 ******************************************************************************/
37
38 int
rt_pixtype_size(rt_pixtype pixtype)39 rt_pixtype_size(rt_pixtype pixtype) {
40 int pixbytes = -1;
41
42 switch (pixtype) {
43 case PT_1BB:
44 case PT_2BUI:
45 case PT_4BUI:
46 case PT_8BSI:
47 case PT_8BUI:
48 pixbytes = 1;
49 break;
50 case PT_16BSI:
51 case PT_16BUI:
52 pixbytes = 2;
53 break;
54 case PT_32BSI:
55 case PT_32BUI:
56 case PT_32BF:
57 pixbytes = 4;
58 break;
59 case PT_64BF:
60 pixbytes = 8;
61 break;
62 default:
63 rterror("rt_pixtype_size: Unknown pixeltype %d", pixtype);
64 pixbytes = -1;
65 break;
66 }
67
68 RASTER_DEBUGF(3, "Pixel type = %s and size = %d bytes",
69 rt_pixtype_name(pixtype), pixbytes);
70
71 return pixbytes;
72 }
73
74 int
rt_pixtype_alignment(rt_pixtype pixtype)75 rt_pixtype_alignment(rt_pixtype pixtype) {
76 return rt_pixtype_size(pixtype);
77 }
78
79 rt_pixtype
rt_pixtype_index_from_name(const char * pixname)80 rt_pixtype_index_from_name(const char* pixname) {
81 assert(pixname);
82
83 if (strcmp(pixname, "1BB") == 0)
84 return PT_1BB;
85 else if (strcmp(pixname, "2BUI") == 0)
86 return PT_2BUI;
87 else if (strcmp(pixname, "4BUI") == 0)
88 return PT_4BUI;
89 else if (strcmp(pixname, "8BSI") == 0)
90 return PT_8BSI;
91 else if (strcmp(pixname, "8BUI") == 0)
92 return PT_8BUI;
93 else if (strcmp(pixname, "16BSI") == 0)
94 return PT_16BSI;
95 else if (strcmp(pixname, "16BUI") == 0)
96 return PT_16BUI;
97 else if (strcmp(pixname, "32BSI") == 0)
98 return PT_32BSI;
99 else if (strcmp(pixname, "32BUI") == 0)
100 return PT_32BUI;
101 else if (strcmp(pixname, "32BF") == 0)
102 return PT_32BF;
103 else if (strcmp(pixname, "64BF") == 0)
104 return PT_64BF;
105
106 return PT_END;
107 }
108
109 const char*
rt_pixtype_name(rt_pixtype pixtype)110 rt_pixtype_name(rt_pixtype pixtype) {
111 switch (pixtype) {
112 case PT_1BB:
113 return "1BB";
114 case PT_2BUI:
115 return "2BUI";
116 case PT_4BUI:
117 return "4BUI";
118 case PT_8BSI:
119 return "8BSI";
120 case PT_8BUI:
121 return "8BUI";
122 case PT_16BSI:
123 return "16BSI";
124 case PT_16BUI:
125 return "16BUI";
126 case PT_32BSI:
127 return "32BSI";
128 case PT_32BUI:
129 return "32BUI";
130 case PT_32BF:
131 return "32BF";
132 case PT_64BF:
133 return "64BF";
134 default:
135 rterror("rt_pixtype_name: Unknown pixeltype %d", pixtype);
136 return "Unknown";
137 }
138 }
139
140 /**
141 * Return minimum value possible for pixel type
142 *
143 * @param pixtype : the pixel type to get minimum possible value for
144 *
145 * @return the minimum possible value for the pixel type.
146 */
147 double
rt_pixtype_get_min_value(rt_pixtype pixtype)148 rt_pixtype_get_min_value(rt_pixtype pixtype) {
149 switch (pixtype) {
150 case PT_1BB: {
151 return (double) rt_util_clamp_to_1BB((double) CHAR_MIN);
152 }
153 case PT_2BUI: {
154 return 0;
155 }
156 case PT_4BUI: {
157 return 0;
158 }
159 case PT_8BUI: {
160 return 0;
161 }
162 case PT_8BSI: {
163 return (double) rt_util_clamp_to_8BSI((double) SCHAR_MIN);
164 }
165 case PT_16BSI: {
166 return (double) rt_util_clamp_to_16BSI((double) SHRT_MIN);
167 }
168 case PT_16BUI: {
169 return 0;
170 }
171 case PT_32BSI: {
172 return (double) rt_util_clamp_to_32BSI((double) INT_MIN);
173 }
174 case PT_32BUI: {
175 return 0;
176 }
177 case PT_32BF: {
178 return (double) -FLT_MAX;
179 }
180 case PT_64BF: {
181 return (double) -DBL_MAX;
182 }
183 default: {
184 rterror("rt_pixtype_get_min_value: Unknown pixeltype %d", pixtype);
185 return (double) rt_util_clamp_to_8BUI((double) CHAR_MIN);
186 }
187 }
188 }
189
190 /**
191 * Returns 1 if clamped values are equal, 0 if not equal, -1 if error
192 *
193 * @param pixtype : the pixel type to clamp the provided values
194 * @param val : value to compare to reference value
195 * @param refval : reference value to be compared with
196 * @param isequal : non-zero if clamped values are equal, 0 otherwise
197 *
198 * @return ES_NONE on success, ES_ERROR on error
199 */
rt_pixtype_compare_clamped_values(rt_pixtype pixtype,double val,double refval,int * isequal)200 rt_errorstate rt_pixtype_compare_clamped_values(
201 rt_pixtype pixtype,
202 double val, double refval,
203 int *isequal
204 ) {
205 assert(isequal != NULL);
206 *isequal = 0;
207
208 switch (pixtype) {
209 case PT_1BB:
210 if (rt_util_clamp_to_1BB(val) == rt_util_clamp_to_1BB(refval))
211 *isequal = 1;
212 break;
213 case PT_2BUI:
214 if (rt_util_clamp_to_2BUI(val) == rt_util_clamp_to_2BUI(refval))
215 *isequal = 1;
216 break;
217 case PT_4BUI:
218 if (rt_util_clamp_to_4BUI(val) == rt_util_clamp_to_4BUI(refval))
219 *isequal = 1;
220 break;
221 case PT_8BSI:
222 if (rt_util_clamp_to_8BSI(val) == rt_util_clamp_to_8BSI(refval))
223 *isequal = 1;
224 break;
225 case PT_8BUI:
226 if (rt_util_clamp_to_8BUI(val) == rt_util_clamp_to_8BUI(refval))
227 *isequal = 1;
228 break;
229 case PT_16BSI:
230 if (rt_util_clamp_to_16BSI(val) == rt_util_clamp_to_16BSI(refval))
231 *isequal = 1;
232 break;
233 case PT_16BUI:
234 if (rt_util_clamp_to_16BUI(val) == rt_util_clamp_to_16BUI(refval))
235 *isequal = 1;
236 break;
237 case PT_32BSI:
238 if (rt_util_clamp_to_32BSI(val) == rt_util_clamp_to_32BSI(refval))
239 *isequal = 1;
240 break;
241 case PT_32BUI:
242 if (rt_util_clamp_to_32BUI(val) == rt_util_clamp_to_32BUI(refval))
243 *isequal = 1;
244 break;
245 case PT_32BF:
246 if (FLT_EQ(rt_util_clamp_to_32F(val), rt_util_clamp_to_32F(refval)))
247 *isequal = 1;
248 break;
249 case PT_64BF:
250 if (FLT_EQ(val, refval))
251 *isequal = 1;
252 break;
253 default:
254 rterror("rt_pixtype_compare_clamped_values: Unknown pixeltype %d", pixtype);
255 return ES_ERROR;
256 }
257
258 return ES_NONE;
259 }
260
261 /******************************************************************************
262 * rt_pixel
263 ******************************************************************************/
264
265 /*
266 * Convert an array of rt_pixel objects to two 2D arrays of value and NODATA.
267 * The dimensions of the returned 2D array are [Y][X], going by row Y and
268 * then column X.
269 *
270 * @param npixel : array of rt_pixel objects
271 * @param count : number of elements in npixel
272 * @param mask : mask to be respected when retruning array
273 * @param x : the column of the center pixel (0-based)
274 * @param y : the line of the center pixel (0-based)
275 * @param distancex : the number of pixels around the specified pixel
276 * along the X axis
277 * @param distancey : the number of pixels around the specified pixel
278 * along the Y axis
279 * @param value : pointer to pointer for 2D value array
280 * @param nodata : pointer to pointer for 2D NODATA array
281 * @param dimx : size of value and nodata along the X axis
282 * @param dimy : size of value and nodata along the Y axis
283 *
284 * @return ES_NONE on success, ES_ERROR on error
285 */
rt_pixel_set_to_array(rt_pixel npixel,uint32_t count,rt_mask mask,int x,int y,uint16_t distancex,uint16_t distancey,double *** value,int *** nodata,int * dimx,int * dimy)286 rt_errorstate rt_pixel_set_to_array(
287 rt_pixel npixel, uint32_t count, rt_mask mask,
288 int x, int y,
289 uint16_t distancex, uint16_t distancey,
290 double ***value,
291 int ***nodata,
292 int *dimx, int *dimy
293 ) {
294 uint32_t i;
295 uint32_t j;
296 uint32_t dim[2] = {0};
297 double **values = NULL;
298 int **nodatas = NULL;
299 int zero[2] = {0};
300 int _x;
301 int _y;
302
303 assert(npixel != NULL && count > 0);
304 assert(value != NULL);
305 assert(nodata != NULL);
306
307 /* dimensions */
308 dim[0] = distancex * 2 + 1;
309 dim[1] = distancey * 2 + 1;
310 RASTER_DEBUGF(4, "dimensions = %d x %d", dim[0], dim[1]);
311
312 /* make sure that the dimx and dimy match mask */
313 if (mask != NULL) {
314 if (dim[0] != mask->dimx || dim[1] != mask->dimy) {
315 rterror("rt_pixel_set_array: mask dimensions %d x %d do not match given dims %d x %d", mask->dimx, mask->dimy, dim[0], dim[1]);
316 return ES_ERROR;
317 }
318
319 if (mask->values == NULL || mask->nodata == NULL) {
320 rterror("rt_pixel_set_array: Invalid mask");
321 return ES_ERROR;
322 }
323
324 }
325
326 /* establish 2D arrays (Y axis) */
327 values = rtalloc(sizeof(double *) * dim[1]);
328 nodatas = rtalloc(sizeof(int *) * dim[1]);
329
330 if (values == NULL || nodatas == NULL) {
331 rterror("rt_pixel_set_to_array: Could not allocate memory for 2D array");
332 return ES_ERROR;
333 }
334
335 /* initialize X axis */
336 for (i = 0; i < dim[1]; i++) {
337 values[i] = rtalloc(sizeof(double) * dim[0]);
338 nodatas[i] = rtalloc(sizeof(int) * dim[0]);
339
340 if (values[i] == NULL || nodatas[i] == NULL) {
341 rterror("rt_pixel_set_to_array: Could not allocate memory for dimension of 2D array");
342
343 if (values[i] == NULL) {
344 for (j = 0; j < i; j++) {
345 rtdealloc(values[j]);
346 rtdealloc(nodatas[j]);
347 }
348 }
349 else {
350 for (j = 0; j <= i; j++) {
351 rtdealloc(values[j]);
352 if (j < i)
353 rtdealloc(nodatas[j]);
354 }
355 }
356
357 rtdealloc(values);
358 rtdealloc(nodatas);
359
360 return ES_ERROR;
361 }
362
363 /* set values to 0 */
364 memset(values[i], 0, sizeof(double) * dim[0]);
365
366 /* set nodatas to 1 */
367 for (j = 0; j < dim[0]; j++)
368 nodatas[i][j] = 1;
369 }
370
371 /* get 0,0 of grid */
372 zero[0] = x - distancex;
373 zero[1] = y - distancey;
374
375 /* populate 2D arrays */
376 for (i = 0; i < count; i++) {
377 if (npixel[i].nodata)
378 continue;
379
380 _x = npixel[i].x - zero[0];
381 _y = npixel[i].y - zero[1];
382
383 RASTER_DEBUGF(4, "absolute x,y: %d x %d", npixel[i].x, npixel[i].y);
384 RASTER_DEBUGF(4, "relative x,y: %d x %d", _x, _y);
385
386 /* no mask */
387 if (mask == NULL) {
388 values[_y][_x] = npixel[i].value;
389 nodatas[_y][_x] = 0;
390 }
391 /* mask */
392 else {
393 /* unweighted (boolean) mask */
394 if (mask->weighted == 0) {
395 /* pixel is set to zero or nodata */
396 if (FLT_EQ(mask->values[_y][_x], 0.0) || mask->nodata[_y][_x] == 1)
397 {
398 values[_y][_x] = 0;
399 nodatas[_y][_x] = 1;
400 }
401 /* use pixel */
402 else {
403 values[_y][_x] = npixel[i].value;
404 nodatas[_y][_x] = 0;
405 }
406 }
407 /* weighted mask */
408 else {
409 /* nodata */
410 if(mask->nodata[_y][_x] == 1) {
411 values[_y][_x] = 0;
412 nodatas[_y][_x] = 1;
413 }
414 /* apply weight to pixel value */
415 else {
416 values[_y][_x] = npixel[i].value * mask->values[_y][_x];
417 nodatas[_y][_x] = 0;
418 }
419 }
420 }
421
422 RASTER_DEBUGF(4, "(x, y, nodata, value) = (%d, %d, %d, %f)", _x, _y, nodatas[_y][_x], values[_y][_x]);
423 }
424
425 *value = &(*values);
426 *nodata = &(*nodatas);
427 if (dimx != NULL)
428 *dimx = dim[0];
429 if (dimy != NULL)
430 *dimy = dim[1];
431
432 return ES_NONE;
433 }
434