1 /******************************************************************************
2 * $Id$
3 *
4 * Project: MapServer
5 * Purpose: Code for drawing GDAL raster layers. Called from
6 * msDrawRasterLayerLow() in mapraster.c.
7 * Author: Frank Warmerdam, warmerdam@pobox.com
8 *
9 ******************************************************************************
10 * Copyright (c) 1996-2005 Regents of the University of Minnesota.
11 *
12 * Permission is hereby granted, free of charge, to any person obtaining a
13 * copy of this software and associated documentation files (the "Software"),
14 * to deal in the Software without restriction, including without limitation
15 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
16 * and/or sell copies of the Software, and to permit persons to whom the
17 * Software is furnished to do so, subject to the following conditions:
18 *
19 * The above copyright notice and this permission notice shall be included in
20 * all copies of this Software or works derived from this Software.
21 *
22 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
23 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
24 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
25 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
26 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
27 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
28 * DEALINGS IN THE SOFTWARE.
29 *****************************************************************************/
30
31 #include <assert.h>
32 #include "mapserver.h"
33 #include "mapresample.h"
34 #include "mapthread.h"
35 #include "maptime.h"
36
37
38 extern int InvGeoTransform( double *gt_in, double *gt_out );
39
40 #define MAXCOLORS 256
41 #define GEO_TRANS(tr,x,y) ((tr)[0]+(tr)[1]*(x)+(tr)[2]*(y))
42 #define SKIP_MASK(x,y) (mask_rb && !*(mask_rb->data.rgba.a+(y)*mask_rb->data.rgba.row_step+(x)*mask_rb->data.rgba.pixel_step))
43
44 #include "gdal.h"
45 #include "cpl_string.h"
46
47 #include "gdal_alg.h"
48
49 static int
50 LoadGDALImages( GDALDatasetH hDS, int band_numbers[4], int band_count,
51 layerObj *layer,
52 int src_xoff, int src_yoff, int src_xsize, int src_ysize,
53 GByte *pabyBuffer,
54 int dst_xsize, int dst_ysize,
55 int *pbHaveRGBNoData,
56 int *pnNoData1, int *pnNoData2, int *pnNoData3 );
57 static int
58 msDrawRasterLayerGDAL_RawMode(
59 mapObj *map, layerObj *layer, imageObj *image, GDALDatasetH hDS,
60 int src_xoff, int src_yoff, int src_xsize, int src_ysize,
61 int dst_xoff, int dst_yoff, int dst_xsize, int dst_ysize );
62
63 static int
64 msDrawRasterLayerGDAL_16BitClassification(
65 mapObj *map, layerObj *layer, rasterBufferObj *rb,
66 GDALDatasetH hDS, GDALRasterBandH hBand,
67 int src_xoff, int src_yoff, int src_xsize, int src_ysize,
68 int dst_xoff, int dst_yoff, int dst_xsize, int dst_ysize );
69
70
71 /*
72 * rasterBufferObj setting macros.
73 */
74
75
76
77 /************************************************************************/
78 /* msDrawRasterLayerGDAL() */
79 /************************************************************************/
80
msDrawRasterLayerGDAL(mapObj * map,layerObj * layer,imageObj * image,rasterBufferObj * rb,void * hDSVoid)81 int msDrawRasterLayerGDAL(mapObj *map, layerObj *layer, imageObj *image,
82 rasterBufferObj *rb, void *hDSVoid )
83
84 {
85 int i,j, k; /* loop counters */
86 int cmap[MAXCOLORS];
87 #ifndef NDEBUG
88 int cmap_set = FALSE;
89 #endif
90 unsigned char rb_cmap[4][MAXCOLORS];
91 double adfGeoTransform[6], adfInvGeoTransform[6];
92 int dst_xoff, dst_yoff, dst_xsize, dst_ysize;
93 int src_xoff, src_yoff, src_xsize, src_ysize;
94 double llx, lly, urx, ury;
95 rectObj copyRect, mapRect;
96 unsigned char *pabyRaw1=NULL, *pabyRaw2=NULL, *pabyRaw3=NULL,
97 *pabyRawAlpha = NULL;
98 int classified = FALSE;
99 int red_band=0, green_band=0, blue_band=0, alpha_band=0;
100 int band_count, band_numbers[4];
101 GDALDatasetH hDS = hDSVoid;
102 GDALColorTableH hColorMap=NULL;
103 GDALRasterBandH hBand1=NULL, hBand2=NULL, hBand3=NULL, hBandAlpha=NULL;
104 int bHaveRGBNoData = FALSE;
105 int nNoData1=-1,nNoData2=-1,nNoData3=-1;
106 rasterBufferObj *mask_rb = NULL;
107 rasterBufferObj s_mask_rb;
108 if(layer->mask) {
109 int ret;
110 layerObj *maskLayer = GET_LAYER(map, msGetLayerIndex(map,layer->mask));
111 mask_rb = &s_mask_rb;
112 memset(mask_rb, 0, sizeof(s_mask_rb));
113 ret = MS_IMAGE_RENDERER(maskLayer->maskimage)->getRasterBufferHandle(maskLayer->maskimage,mask_rb);
114 if(ret != MS_SUCCESS)
115 return -1;
116 }
117
118 /*only support rawdata and pluggable renderers*/
119 assert(MS_RENDERER_RAWDATA(image->format) || (MS_RENDERER_PLUGIN(image->format) && rb));
120
121
122 memset( cmap, 0xff, MAXCOLORS * sizeof(int) );
123 memset( rb_cmap, 0, sizeof(rb_cmap) );
124
125 /* -------------------------------------------------------------------- */
126 /* Test the image format instead of the map format. */
127 /* Normally the map format and the image format should be the */
128 /* same but In somes cases like swf and pdf support, a temporary */
129 /* GD image object is created and used to render raster layers */
130 /* and then dumped into the SWF or the PDF file. */
131 /* -------------------------------------------------------------------- */
132
133
134 src_xsize = GDALGetRasterXSize( hDS );
135 src_ysize = GDALGetRasterYSize( hDS );
136
137 /*
138 * If the RAW_WINDOW attribute is set, use that to establish what to
139 * load. This is normally just set by the mapresample.c module to avoid
140 * problems with rotated maps.
141 */
142
143 if( CSLFetchNameValue( layer->processing, "RAW_WINDOW" ) != NULL ) {
144 char **papszTokens =
145 CSLTokenizeString(
146 CSLFetchNameValue( layer->processing, "RAW_WINDOW" ) );
147
148 if( layer->debug )
149 msDebug( "msDrawGDAL(%s): using RAW_WINDOW=%s, dst=0,0,%d,%d\n",
150 layer->name,
151 CSLFetchNameValue( layer->processing, "RAW_WINDOW" ),
152 image->width, image->height );
153
154 if( CSLCount(papszTokens) != 4 ) {
155 CSLDestroy( papszTokens );
156 msSetError( MS_IMGERR, "RAW_WINDOW PROCESSING directive corrupt.",
157 "msDrawGDAL()" );
158 return -1;
159 }
160
161 src_xoff = atoi(papszTokens[0]);
162 src_yoff = atoi(papszTokens[1]);
163 src_xsize = atoi(papszTokens[2]);
164 src_ysize = atoi(papszTokens[3]);
165
166 dst_xoff = 0;
167 dst_yoff = 0;
168 dst_xsize = image->width;
169 dst_ysize = image->height;
170
171 CSLDestroy( papszTokens );
172 }
173
174 /*
175 * Compute the georeferenced window of overlap, and do nothing if there
176 * is no overlap between the map extents, and the file we are operating on.
177 */
178 else if( layer->transform ) {
179 int dst_lrx, dst_lry;
180
181 if( layer->debug )
182 msDebug( "msDrawRasterLayerGDAL(): Entering transform.\n" );
183
184 msGetGDALGeoTransform( hDS, map, layer, adfGeoTransform );
185 InvGeoTransform( adfGeoTransform, adfInvGeoTransform );
186
187 mapRect = map->extent;
188
189 mapRect.minx -= map->cellsize*0.5;
190 mapRect.maxx += map->cellsize*0.5;
191 mapRect.miny -= map->cellsize*0.5;
192 mapRect.maxy += map->cellsize*0.5;
193
194 copyRect = mapRect;
195
196 if( copyRect.minx < GEO_TRANS(adfGeoTransform,0,src_ysize) )
197 copyRect.minx = GEO_TRANS(adfGeoTransform,0,src_ysize);
198 if( copyRect.maxx > GEO_TRANS(adfGeoTransform,src_xsize,0) )
199 copyRect.maxx = GEO_TRANS(adfGeoTransform,src_xsize,0);
200
201 if( copyRect.miny < GEO_TRANS(adfGeoTransform+3,0,src_ysize) )
202 copyRect.miny = GEO_TRANS(adfGeoTransform+3,0,src_ysize);
203 if( copyRect.maxy > GEO_TRANS(adfGeoTransform+3,src_xsize,0) )
204 copyRect.maxy = GEO_TRANS(adfGeoTransform+3,src_xsize,0);
205
206 if( copyRect.minx >= copyRect.maxx || copyRect.miny >= copyRect.maxy ) {
207 if( layer->debug )
208 msDebug( "msDrawRasterLayerGDAL(): Error in overlap calculation.\n" );
209 return 0;
210 }
211
212 /*
213 * Copy the source and destination raster coordinates.
214 */
215 llx = GEO_TRANS(adfInvGeoTransform+0,copyRect.minx,copyRect.miny);
216 lly = GEO_TRANS(adfInvGeoTransform+3,copyRect.minx,copyRect.miny);
217 urx = GEO_TRANS(adfInvGeoTransform+0,copyRect.maxx,copyRect.maxy);
218 ury = GEO_TRANS(adfInvGeoTransform+3,copyRect.maxx,copyRect.maxy);
219
220 src_xoff = MS_MAX(0,(int) floor(llx+0.5));
221 src_yoff = MS_MAX(0,(int) floor(ury+0.5));
222 src_xsize = MS_MIN(MS_MAX(0,(int) (urx - llx + 0.5)),
223 GDALGetRasterXSize(hDS) - src_xoff);
224 src_ysize = MS_MIN(MS_MAX(0,(int) (lly - ury + 0.5)),
225 GDALGetRasterYSize(hDS) - src_yoff);
226
227 /* We want very small windows to use at least one source pixel (#4172) */
228 if( src_xsize == 0 && (urx - llx) > 0.0 ) {
229 src_xsize = 1;
230 src_xoff = MS_MIN(src_xoff,GDALGetRasterXSize(hDS)-1);
231 }
232 if( src_ysize == 0 && (lly - ury) > 0.0 ) {
233 src_ysize = 1;
234 src_yoff = MS_MIN(src_yoff,GDALGetRasterYSize(hDS)-1);
235 }
236
237 if( src_xsize == 0 || src_ysize == 0 ) {
238 if( layer->debug )
239 msDebug( "msDrawRasterLayerGDAL(): no apparent overlap between map view and this window(1).\n" );
240 return 0;
241 }
242
243 if (map->cellsize == 0) {
244 if( layer->debug )
245 msDebug( "msDrawRasterLayerGDAL(): Cellsize can't be 0.\n" );
246 return 0;
247 }
248
249 dst_xoff = (int) ((copyRect.minx - mapRect.minx) / map->cellsize);
250 dst_yoff = (int) ((mapRect.maxy - copyRect.maxy) / map->cellsize);
251
252 dst_lrx = (int) ((copyRect.maxx - mapRect.minx) / map->cellsize + 0.5);
253 dst_lry = (int) ((mapRect.maxy - copyRect.miny) / map->cellsize + 0.5);
254 dst_lrx = MS_MAX(0,MS_MIN(image->width,dst_lrx));
255 dst_lry = MS_MAX(0,MS_MIN(image->height,dst_lry));
256
257 dst_xsize = MS_MAX(0,MS_MIN(image->width,dst_lrx - dst_xoff));
258 dst_ysize = MS_MAX(0,MS_MIN(image->height,dst_lry - dst_yoff));
259
260 if( dst_xsize == 0 || dst_ysize == 0 ) {
261 if( layer->debug )
262 msDebug( "msDrawRasterLayerGDAL(): no apparent overlap between map view and this window(2).\n" );
263 return 0;
264 }
265
266 if( layer->debug )
267 msDebug( "msDrawRasterLayerGDAL(): src=%d,%d,%d,%d, dst=%d,%d,%d,%d\n",
268 src_xoff, src_yoff, src_xsize, src_ysize,
269 dst_xoff, dst_yoff, dst_xsize, dst_ysize );
270 #ifndef notdef
271 if( layer->debug ) {
272 double d_src_xoff, d_src_yoff, geo_x, geo_y;
273
274 geo_x = mapRect.minx + dst_xoff * map->cellsize;
275 geo_y = mapRect.maxy - dst_yoff * map->cellsize;
276
277 d_src_xoff = (geo_x - adfGeoTransform[0]) / adfGeoTransform[1];
278 d_src_yoff = (geo_y - adfGeoTransform[3]) / adfGeoTransform[5];
279
280 msDebug( "msDrawRasterLayerGDAL(): source raster PL (%.3f,%.3f) for dst PL (%d,%d).\n",
281 d_src_xoff, d_src_yoff,
282 dst_xoff, dst_yoff );
283 }
284 #endif
285 }
286
287 /*
288 * If layer transforms are turned off, just map 1:1.
289 */
290 else {
291 dst_xoff = src_xoff = 0;
292 dst_yoff = src_yoff = 0;
293 dst_xsize = src_xsize = MS_MIN(image->width,src_xsize);
294 dst_ysize = src_ysize = MS_MIN(image->height,src_ysize);
295 }
296
297 /*
298 * In RAWDATA mode we don't fool with colors. Do the raw processing,
299 * and return from the function early.
300 */
301 if( MS_RENDERER_RAWDATA( image->format ) ) {
302 return msDrawRasterLayerGDAL_RawMode(
303 map, layer, image, hDS,
304 src_xoff, src_yoff, src_xsize, src_ysize,
305 dst_xoff, dst_yoff, dst_xsize, dst_ysize );
306 }
307
308 /*
309 * Is this image classified? We consider it classified if there are
310 * classes with an expression string *or* a color range. We don't want
311 * to treat the raster as classified if there is just a bogus class here
312 * to force inclusion in the legend.
313 */
314 for( i = 0; i < layer->numclasses; i++ ) {
315 int s;
316
317 /* change colour based on colour range? */
318 for(s=0; s<layer->class[i]->numstyles; s++) {
319 if( MS_VALID_COLOR(layer->class[i]->styles[s]->mincolor)
320 && MS_VALID_COLOR(layer->class[i]->styles[s]->maxcolor) ) {
321 classified = TRUE;
322 break;
323 }
324 }
325
326 if( layer->class[i]->expression.string != NULL ) {
327 classified = TRUE;
328 break;
329 }
330 }
331
332 /*
333 * Set up the band selection. We look for a BANDS directive in the
334 * the PROCESSING options. If not found we default to grey=1, grey=1,alpha=2,
335 * red=1,green=2,blue=3 or red=1,green=2,blue=3,alpha>=4.
336 */
337
338 if( CSLFetchNameValue( layer->processing, "BANDS" ) == NULL ) {
339 const int gdal_band_count = GDALGetRasterCount( hDS );
340 red_band = 1;
341
342 if( gdal_band_count >= 4 )
343 {
344 /* The alpha band is not necessarily the 4th one */
345 for( i = 4; i <= gdal_band_count; i ++ ) {
346 if( GDALGetRasterColorInterpretation(
347 GDALGetRasterBand( hDS, i ) ) == GCI_AlphaBand ) {
348 alpha_band = i;
349 break;
350 }
351 }
352 }
353
354 if( gdal_band_count >= 3 ) {
355 green_band = 2;
356 blue_band = 3;
357 }
358
359 if( gdal_band_count == 2
360 && GDALGetRasterColorInterpretation(
361 GDALGetRasterBand( hDS, 2 ) ) == GCI_AlphaBand )
362 alpha_band = 2;
363
364 hBand1 = GDALGetRasterBand( hDS, red_band );
365 if( classified
366 || GDALGetRasterColorTable( hBand1 ) != NULL ) {
367 alpha_band = 0;
368 green_band = 0;
369 blue_band = 0;
370 }
371 } else {
372 int *band_list;
373
374 band_list = msGetGDALBandList( layer, hDS, 4, &band_count );
375 if( band_list == NULL )
376 return -1;
377
378 if( band_count > 0 )
379 red_band = band_list[0];
380 else
381 red_band = 0;
382 if( band_count > 2 ) {
383 green_band = band_list[1];
384 blue_band = band_list[2];
385 } else {
386 green_band = 0;
387 blue_band = 0;
388 }
389
390 if( band_count > 3 )
391 alpha_band = band_list[3];
392 else
393 alpha_band = 0;
394
395 free( band_list );
396 }
397
398 band_numbers[0] = red_band;
399 band_numbers[1] = green_band;
400 band_numbers[2] = blue_band;
401 band_numbers[3] = 0;
402
403 if( blue_band != 0 && alpha_band != 0 ) {
404 band_numbers[3] = alpha_band;
405 band_count = 4;
406 } else if( blue_band != 0 && alpha_band == 0 )
407 band_count = 3;
408 else if( alpha_band != 0 ) {
409 band_numbers[1] = alpha_band;
410 band_count = 2;
411 } else
412 band_count = 1;
413
414 if( layer->debug > 1 || (layer->debug > 0 && green_band != 0) ) {
415 msDebug( "msDrawRasterLayerGDAL(): red,green,blue,alpha bands = %d,%d,%d,%d\n",
416 red_band, green_band, blue_band, alpha_band );
417 }
418
419 /*
420 * Get band handles for PC256, RGB or RGBA cases.
421 */
422 hBand1 = GDALGetRasterBand( hDS, red_band );
423 if( hBand1 == NULL )
424 return -1;
425
426 hBand2 = hBand3 = hBandAlpha = NULL;
427
428 if( green_band != 0 ) {
429 hBand1 = GDALGetRasterBand( hDS, red_band );
430 hBand2 = GDALGetRasterBand( hDS, green_band );
431 hBand3 = GDALGetRasterBand( hDS, blue_band );
432 if( hBand1 == NULL || hBand2 == NULL || hBand3 == NULL )
433 return -1;
434 }
435
436 if( alpha_band != 0 )
437 hBandAlpha = GDALGetRasterBand( hDS, alpha_band );
438
439 /*
440 * The logic for a classification rendering of non-8bit raster bands
441 * is sufficiently different than the normal mechanism of loading
442 * into an 8bit buffer, that we isolate it into it's own subfunction.
443 */
444 if( classified
445 && hBand1 != NULL && GDALGetRasterDataType( hBand1 ) != GDT_Byte ) {
446 return msDrawRasterLayerGDAL_16BitClassification(
447 map, layer, rb, hDS, hBand1,
448 src_xoff, src_yoff, src_xsize, src_ysize,
449 dst_xoff, dst_yoff, dst_xsize, dst_ysize );
450 }
451
452 /*
453 * Get colormap for this image. If there isn't one, and we have only
454 * one band create a greyscale colormap.
455 */
456 if( hBand2 != NULL )
457 hColorMap = NULL;
458 else {
459 hColorMap = GDALGetRasterColorTable( hBand1 );
460 if( hColorMap != NULL )
461 hColorMap = GDALCloneColorTable( hColorMap );
462 else if( hBand2 == NULL ) {
463 hColorMap = GDALCreateColorTable( GPI_RGB );
464
465 for( i = 0; i < 256; i++ ) {
466 colorObj pixel;
467 GDALColorEntry sEntry;
468
469 pixel.red = i;
470 pixel.green = i;
471 pixel.blue = i;
472
473 if(MS_COMPARE_COLORS(pixel, layer->offsite)) {
474 sEntry.c1 = 0;
475 sEntry.c2 = 0;
476 sEntry.c3 = 0;
477 sEntry.c4 = 0; /* alpha set to zero */
478 } else {
479 sEntry.c1 = i;
480 sEntry.c2 = i;
481 sEntry.c3 = i;
482 sEntry.c4 = 255;
483 }
484
485 GDALSetColorEntry( hColorMap, i, &sEntry );
486 }
487 }
488
489 /*
490 ** If we have a known NODATA value, mark it now as transparent.
491 */
492 {
493 int bGotNoData;
494 double dfNoDataValue = msGetGDALNoDataValue( layer, hBand1,
495 &bGotNoData);
496
497 if( bGotNoData && dfNoDataValue >= 0
498 && dfNoDataValue < GDALGetColorEntryCount( hColorMap ) ) {
499 GDALColorEntry sEntry;
500
501 memcpy( &sEntry,
502 GDALGetColorEntry( hColorMap, (int) dfNoDataValue ),
503 sizeof(GDALColorEntry) );
504
505 sEntry.c4 = 0;
506 GDALSetColorEntry( hColorMap, (int) dfNoDataValue, &sEntry );
507 }
508 }
509 }
510
511 /*
512 * Setup the mapping between source eight bit pixel values, and the
513 * output images color table. There are two general cases, where the
514 * class colors are provided by the MAP file, or where we use the native
515 * color table.
516 */
517 if( classified ) {
518 int c, color_count;
519 const char* pszRangeColorspace = msLayerGetProcessingKey( layer, "RANGE_COLORSPACE" );
520 colorspace iRangeColorspace;
521
522 #ifndef NDEBUG
523 cmap_set = TRUE;
524 #endif
525
526 if( hColorMap == NULL ) {
527 msSetError(MS_IOERR,
528 "Attempt to classify 24bit image, this is unsupported.",
529 "drawGDAL()");
530 return -1;
531 }
532
533 if(!pszRangeColorspace || !strcasecmp(pszRangeColorspace, "RGB")) {
534 iRangeColorspace = MS_COLORSPACE_RGB;
535 } else if(!strcasecmp(pszRangeColorspace, "HSL")) {
536 iRangeColorspace = MS_COLORSPACE_HSL;
537 } else {
538 msSetError(MS_MISCERR,
539 "Unknown RANGE_COLORSPACE \"%s\", expecting RGB or HSL",
540 "drawGDAL()", pszRangeColorspace);
541 GDALDestroyColorTable( hColorMap );
542 return -1;
543 }
544
545 color_count = MS_MIN(256,GDALGetColorEntryCount(hColorMap));
546 for(i=0; i < color_count; i++) {
547 colorObj pixel;
548 int colormap_index;
549 GDALColorEntry sEntry;
550
551 GDALGetColorEntryAsRGB( hColorMap, i, &sEntry );
552
553 pixel.red = sEntry.c1;
554 pixel.green = sEntry.c2;
555 pixel.blue = sEntry.c3;
556 colormap_index = i;
557
558 if(!MS_COMPARE_COLORS(pixel, layer->offsite)) {
559 c = msGetClass(layer, &pixel, colormap_index);
560
561 if(c != -1) { /* belongs to any class */
562 int s;
563
564 /* change colour based on colour range? Currently we
565 only address the greyscale case properly. */
566
567 if( MS_VALID_COLOR(layer->class[c]->styles[0]->mincolor) ) {
568 for(s=0; s<layer->class[c]->numstyles; s++) {
569 if( MS_VALID_COLOR(layer->class[c]->styles[s]->mincolor)
570 && MS_VALID_COLOR(layer->class[c]->styles[s]->maxcolor)) {
571 if( layer->class[c]->numstyles == 1 || (sEntry.c1 >= layer->class[c]->styles[s]->minvalue
572 && sEntry.c1 <= layer->class[c]->styles[s]->maxvalue )) {
573 msValueToRange(layer->class[c]->styles[s], sEntry.c1, iRangeColorspace);
574 if(MS_VALID_COLOR(layer->class[c]->styles[s]->color)) {
575 rb_cmap[0][i] = layer->class[c]->styles[s]->color.red;
576 rb_cmap[1][i] = layer->class[c]->styles[s]->color.green;
577 rb_cmap[2][i] = layer->class[c]->styles[s]->color.blue;
578 rb_cmap[3][i] = (layer->class[c]->styles[s]->color.alpha != 255)?(layer->class[c]->styles[s]->color.alpha):(255*layer->class[c]->styles[0]->opacity / 100);
579 break;
580 }
581 }
582 }
583 }
584 }
585 else if( MS_TRANSPARENT_COLOR(layer->class[c]->styles[0]->color))
586 /* leave it transparent */;
587
588 else if( MS_VALID_COLOR(layer->class[c]->styles[0]->color)) {
589 rb_cmap[0][i] = layer->class[c]->styles[0]->color.red;
590 rb_cmap[1][i] = layer->class[c]->styles[0]->color.green;
591 rb_cmap[2][i] = layer->class[c]->styles[0]->color.blue;
592 rb_cmap[3][i] = (255*layer->class[c]->styles[0]->opacity / 100);
593 }
594
595 else { /* Use raster color */
596 rb_cmap[0][i] = pixel.red;
597 rb_cmap[1][i] = pixel.green;
598 rb_cmap[2][i] = pixel.blue;
599 rb_cmap[3][i] = 255;
600 }
601 }
602 }
603 }
604 } else if( hBand2 == NULL && hColorMap != NULL && rb->type == MS_BUFFER_BYTE_RGBA ) {
605 int color_count;
606 #ifndef NDEBUG
607 cmap_set = TRUE;
608 #endif
609
610 color_count = MS_MIN(256,GDALGetColorEntryCount(hColorMap));
611
612 for(i=0; i < color_count; i++) {
613 GDALColorEntry sEntry;
614
615 GDALGetColorEntryAsRGB( hColorMap, i, &sEntry );
616
617 if( sEntry.c4 != 0
618 && (!MS_VALID_COLOR( layer->offsite )
619 || layer->offsite.red != sEntry.c1
620 || layer->offsite.green != sEntry.c2
621 || layer->offsite.blue != sEntry.c3 ) ) {
622 rb_cmap[0][i] = sEntry.c1;
623 rb_cmap[1][i] = sEntry.c2;
624 rb_cmap[2][i] = sEntry.c3;
625 rb_cmap[3][i] = sEntry.c4;
626 }
627 }
628 }
629 /*
630 * Allocate imagery buffers.
631 */
632 pabyRaw1 = (unsigned char *) malloc(dst_xsize * dst_ysize * band_count);
633 if( pabyRaw1 == NULL ) {
634 msSetError(MS_MEMERR, "Allocating work image of size %dx%dx%d failed.",
635 "msDrawRasterLayerGDAL()", dst_xsize, dst_ysize, band_count );
636
637 if( hColorMap != NULL )
638 GDALDestroyColorTable( hColorMap );
639 return -1;
640 }
641
642 if( hBand2 != NULL && hBand3 != NULL ) {
643 pabyRaw2 = pabyRaw1 + dst_xsize * dst_ysize * 1;
644 pabyRaw3 = pabyRaw1 + dst_xsize * dst_ysize * 2;
645 }
646
647 if( hBandAlpha != NULL ) {
648 if( hBand2 != NULL )
649 pabyRawAlpha = pabyRaw1 + dst_xsize * dst_ysize * 3;
650 else
651 pabyRawAlpha = pabyRaw1 + dst_xsize * dst_ysize * 1;
652 }
653
654 /*
655 * Load image data into buffers with scaling, etc.
656 */
657 if( LoadGDALImages( hDS, band_numbers, band_count, layer,
658 src_xoff, src_yoff, src_xsize, src_ysize,
659 pabyRaw1, dst_xsize, dst_ysize,
660 &bHaveRGBNoData,
661 &nNoData1, &nNoData2, &nNoData3 ) == -1 ) {
662 free( pabyRaw1 );
663 if( hColorMap != NULL )
664 GDALDestroyColorTable( hColorMap );
665 return -1;
666 }
667
668 if( bHaveRGBNoData && layer->debug )
669 msDebug( "msDrawGDAL(): using RGB nodata values from GDAL dataset.\n" );
670
671 /* -------------------------------------------------------------------- */
672 /* If there was no alpha band, but we have a dataset level mask */
673 /* load it as massage it so it will function as our alpha for */
674 /* transparency purposes. */
675 /* -------------------------------------------------------------------- */
676 if( hBandAlpha == NULL ) {
677 int nMaskFlags = GDALGetMaskFlags(hBand1);
678
679 if( (CSLFetchNameValue( layer->processing, "BANDS" ) == NULL ) &&
680 (nMaskFlags & GMF_PER_DATASET) != 0 &&
681 (nMaskFlags & (GMF_NODATA|GMF_ALL_VALID)) == 0 ) {
682 CPLErr eErr;
683 unsigned char * pabyOrig = pabyRaw1;
684
685
686 if( layer->debug )
687 msDebug( "msDrawGDAL(): using GDAL mask band for alpha.\n" );
688
689 band_count++;
690
691 pabyRaw1 = (unsigned char *)
692 realloc(pabyOrig,dst_xsize * dst_ysize * band_count);
693
694 if( pabyRaw1 == NULL ) {
695 msSetError(MS_MEMERR,
696 "Allocating work image of size %dx%dx%d failed.",
697 "msDrawRasterLayerGDAL()",
698 dst_xsize, dst_ysize, band_count );
699 free(pabyOrig);
700 if( hColorMap != NULL )
701 GDALDestroyColorTable( hColorMap );
702 return -1;
703 }
704
705 if( hBand2 != NULL ) {
706 pabyRaw2 = pabyRaw1 + dst_xsize * dst_ysize * 1;
707 pabyRaw3 = pabyRaw1 + dst_xsize * dst_ysize * 2;
708 pabyRawAlpha = pabyRaw1 + dst_xsize * dst_ysize * 3;
709 } else {
710 pabyRawAlpha = pabyRaw1 + dst_xsize * dst_ysize * 1;
711 }
712
713 hBandAlpha = GDALGetMaskBand(hBand1);
714
715 eErr = GDALRasterIO( hBandAlpha, GF_Read,
716 src_xoff, src_yoff, src_xsize, src_ysize,
717 pabyRawAlpha, dst_xsize, dst_ysize, GDT_Byte, 0,0);
718
719 if( eErr != CE_None ) {
720 msSetError( MS_IOERR, "GDALRasterIO() failed: %s",
721 "drawGDAL()", CPLGetLastErrorMsg() );
722 free( pabyRaw1 );
723 if( hColorMap != NULL )
724 GDALDestroyColorTable( hColorMap );
725 return -1;
726 }
727
728 /* In case the mask is not an alpha channel, expand values of 1 to 255, */
729 /* so we can deal as it was an alpha band afterwards */
730 if ((nMaskFlags & GMF_ALPHA) == 0) {
731 for(i=0; i<dst_xsize * dst_ysize; i++)
732 if (pabyRawAlpha[i])
733 pabyRawAlpha[i] = 255;
734 }
735 }
736 }
737
738 /* -------------------------------------------------------------------- */
739 /* Single band plus colormap and alpha to truecolor. (RB) */
740 /* -------------------------------------------------------------------- */
741 if( hBand2 == NULL && rb->type == MS_BUFFER_BYTE_RGBA && hBandAlpha != NULL ) {
742 assert( cmap_set );
743
744 k = 0;
745 for( i = dst_yoff; i < dst_yoff + dst_ysize; i++ ) {
746 for( j = dst_xoff; j < dst_xoff + dst_xsize; j++ ) {
747 int src_pixel, src_alpha, cmap_alpha, merged_alpha;
748 if(SKIP_MASK(j,i)) {
749 k++;
750 continue;
751 }
752
753 src_pixel = pabyRaw1[k];
754 src_alpha = pabyRawAlpha[k];
755 cmap_alpha = rb_cmap[3][src_pixel];
756
757 merged_alpha = (src_alpha * cmap_alpha) / 255;
758
759 if( merged_alpha < 2 )
760 /* do nothing - transparent */;
761 else if( merged_alpha > 253 ) {
762 RB_SET_PIXEL( rb, j, i,
763 rb_cmap[0][src_pixel],
764 rb_cmap[1][src_pixel],
765 rb_cmap[2][src_pixel],
766 cmap_alpha );
767 } else {
768 RB_MIX_PIXEL( rb, j, i,
769 rb_cmap[0][src_pixel],
770 rb_cmap[1][src_pixel],
771 rb_cmap[2][src_pixel],
772 merged_alpha );
773 }
774 k++;
775 }
776 }
777 }
778
779 /* -------------------------------------------------------------------- */
780 /* Single band plus colormap (no alpha) to truecolor (RB) */
781 /* -------------------------------------------------------------------- */
782 else if( hBand2 == NULL && rb->type == MS_BUFFER_BYTE_RGBA ) {
783 assert( cmap_set );
784
785 k = 0;
786 for( i = dst_yoff; i < dst_yoff + dst_ysize; i++ ) {
787 for( j = dst_xoff; j < dst_xoff + dst_xsize; j++ ) {
788 int src_pixel = pabyRaw1[k++];
789 if(SKIP_MASK(j,i)) {
790 continue;
791 }
792
793 if( rb_cmap[3][src_pixel] > 253 ) {
794 RB_SET_PIXEL( rb, j, i,
795 rb_cmap[0][src_pixel],
796 rb_cmap[1][src_pixel],
797 rb_cmap[2][src_pixel],
798 rb_cmap[3][src_pixel] );
799 } else if( rb_cmap[3][src_pixel] > 1 ) {
800 RB_MIX_PIXEL( rb, j, i,
801 rb_cmap[0][src_pixel],
802 rb_cmap[1][src_pixel],
803 rb_cmap[2][src_pixel],
804 rb_cmap[3][src_pixel] );
805 }
806 }
807 }
808 }
809
810 /* -------------------------------------------------------------------- */
811 /* Input is 3 band RGB. Alpha blending is mixed into the loop */
812 /* since this case is less commonly used and has lots of other */
813 /* overhead. (RB) */
814 /* -------------------------------------------------------------------- */
815 else if( hBand3 != NULL && rb->type == MS_BUFFER_BYTE_RGBA ) {
816 k = 0;
817 for( i = dst_yoff; i < dst_yoff + dst_ysize; i++ ) {
818 for( j = dst_xoff; j < dst_xoff + dst_xsize; j++, k++ ) {
819 if(SKIP_MASK(j,i)) {
820 continue;
821 }
822 if( MS_VALID_COLOR( layer->offsite )
823 && pabyRaw1[k] == layer->offsite.red
824 && pabyRaw2[k] == layer->offsite.green
825 && pabyRaw3[k] == layer->offsite.blue )
826 continue;
827
828 if( bHaveRGBNoData
829 && pabyRaw1[k] == nNoData1
830 && pabyRaw2[k] == nNoData2
831 && pabyRaw3[k] == nNoData3 )
832 continue;
833
834 if( pabyRawAlpha == NULL || pabyRawAlpha[k] == 255 ) {
835 RB_SET_PIXEL( rb, j, i,
836 pabyRaw1[k],
837 pabyRaw2[k],
838 pabyRaw3[k],
839 255 );
840 } else if( pabyRawAlpha[k] != 0 ) {
841 RB_MIX_PIXEL( rb, j, i,
842 pabyRaw1[k],
843 pabyRaw2[k],
844 pabyRaw3[k],
845 pabyRawAlpha[k] );
846 }
847 }
848 }
849 }
850
851
852 /*
853 ** Cleanup
854 */
855 free( pabyRaw1 );
856
857 if( hColorMap != NULL )
858 GDALDestroyColorTable( hColorMap );
859
860 return 0;
861 }
862
863 /************************************************************************/
864 /* ParseDefaultLUT() */
865 /************************************************************************/
866
ParseDefaultLUT(const char * lut_def,GByte * lut,int nMaxValIn)867 static int ParseDefaultLUT( const char *lut_def, GByte *lut, int nMaxValIn )
868
869 {
870 const char *lut_read;
871 int last_in=0, last_out=0, all_done=FALSE;
872
873 /* -------------------------------------------------------------------- */
874 /* Parse definition, applying to lut on the fly. */
875 /* -------------------------------------------------------------------- */
876 lut_read = lut_def;
877 while( !all_done ) {
878 int this_in=0, this_out=0;
879 int lut_i;
880
881 while( isspace(*lut_read) )
882 lut_read++;
883
884 /* if we are at end, assume nMaxValIn:255 */
885 if( *lut_read == '\0' ) {
886 all_done = TRUE;
887 if ( last_in != nMaxValIn ) {
888 this_in = nMaxValIn;
889 this_out = 255;
890 }
891 }
892
893 /* otherwise read "in:out", and skip past */
894 else {
895 this_in = atoi(lut_read);
896 while( *lut_read != ':' && *lut_read )
897 lut_read++;
898 if( *lut_read == ':' )
899 lut_read++;
900 while( isspace(*lut_read) )
901 lut_read++;
902 this_out = atoi(lut_read);
903 while( *lut_read && *lut_read != ',' && *lut_read != '\n' )
904 lut_read++;
905 if( *lut_read == ',' || *lut_read == '\n' )
906 lut_read++;
907 while( isspace(*lut_read) )
908 lut_read++;
909 }
910
911 this_in = MS_MAX(0,MS_MIN(nMaxValIn,this_in));
912 this_out = MS_MAX(0,MS_MIN(255,this_out));
913
914 /* apply linear values from last in:out to this in:out */
915 for( lut_i = last_in; lut_i <= this_in; lut_i++ ) {
916 double ratio;
917
918 if( last_in == this_in )
919 ratio = 1.0;
920 else
921 ratio = (lut_i - last_in) / (double) (this_in - last_in);
922
923 lut[lut_i] = (int)
924 floor(((1.0 - ratio)*last_out + ratio*this_out) + 0.5);
925 }
926
927 last_in = this_in;
928 last_out = this_out;
929 }
930
931 return 0;
932 }
933
934 /************************************************************************/
935 /* LutFromGimpLine() */
936 /************************************************************************/
937
LutFromGimpLine(char * lut_line,GByte * lut)938 static int LutFromGimpLine( char *lut_line, GByte *lut )
939
940 {
941 char wrkLUTDef[1000];
942 int i, count = 0;
943 char **tokens;
944
945 /* cleanup white space at end of line (DOS LF, etc) */
946 i = strlen(lut_line) - 1;
947 while( i > 0 && isspace(lut_line[i]) )
948 lut_line[i--] = '\0';
949
950 while( *lut_line == 10 || *lut_line == 13 )
951 lut_line++;
952
953 /* tokenize line */
954 tokens = CSLTokenizeString( lut_line );
955 if( CSLCount(tokens) != 17 * 2 ) {
956 CSLDestroy( tokens );
957 msSetError(MS_MISCERR,
958 "GIMP curve file appears corrupt.",
959 "LutFromGimpLine()" );
960 return -1;
961 }
962
963 /* Convert to our own format */
964 wrkLUTDef[0] = '\0';
965 for( i = 0; i < 17; i++ ) {
966 if( atoi(tokens[i*2]) >= 0 ) {
967 if( count++ > 0 )
968 strlcat( wrkLUTDef, ",", sizeof(wrkLUTDef));
969
970 snprintf( wrkLUTDef + strlen(wrkLUTDef), sizeof(wrkLUTDef)-strlen(wrkLUTDef),
971 "%s:%s", tokens[i*2], tokens[i*2+1] );
972 }
973 }
974
975 CSLDestroy( tokens );
976
977 return ParseDefaultLUT( wrkLUTDef, lut, 255 );
978 }
979
980 /************************************************************************/
981 /* ParseGimpLUT() */
982 /* */
983 /* Parse a Gimp style LUT. */
984 /************************************************************************/
985
ParseGimpLUT(const char * lut_def,GByte * lut,int iColorIndex)986 static int ParseGimpLUT( const char *lut_def, GByte *lut, int iColorIndex )
987
988 {
989 int i;
990 GByte lutValue[256];
991 GByte lutColorBand[256];
992 char **lines =
993 CSLTokenizeStringComplex( lut_def, "\n", FALSE, FALSE );
994
995 if( !EQUALN(lines[0],"# GIMP Curves File",18)
996 || CSLCount(lines) < 6 ) {
997 msSetError(MS_MISCERR,
998 "GIMP curve file appears corrupt.",
999 "ParseGimpLUT()" );
1000 CSLDestroy( lines );
1001 return -1;
1002 }
1003
1004 /*
1005 * Convert the overall curve, and the color band specific curve into LUTs.
1006 */
1007 if( LutFromGimpLine( lines[1], lutValue ) != 0
1008 || LutFromGimpLine( lines[iColorIndex + 1], lutColorBand ) != 0 ) {
1009 CSLDestroy( lines );
1010 return -1;
1011 }
1012 CSLDestroy( lines );
1013
1014 /*
1015 * Compose the two luts as if the raw value passed through the color band
1016 * specific lut, and then the general value lut. Usually one or the
1017 * other will be the identity mapping, but not always.
1018 */
1019
1020 for( i = 0; i < 256; i++ ) {
1021 lut[i] = lutValue[lutColorBand[i]];
1022 }
1023
1024 return 0;
1025 }
1026
1027 /************************************************************************/
1028 /* LoadLUT() */
1029 /* */
1030 /* Load a LUT according to RFC 21. */
1031 /************************************************************************/
1032
LoadLUT(layerObj * layer,int iColorIndex,char ** ppszLutDef)1033 static int LoadLUT( layerObj *layer, int iColorIndex, char** ppszLutDef )
1034
1035 {
1036 const char *lut_def;
1037 char key[20], lut_def_fromfile[2500];
1038
1039 *ppszLutDef = NULL;
1040
1041 /* -------------------------------------------------------------------- */
1042 /* Get lut specifier from processing directives. Do nothing if */
1043 /* none are found. */
1044 /* -------------------------------------------------------------------- */
1045 sprintf( key, "LUT_%d", iColorIndex );
1046 lut_def = msLayerGetProcessingKey( layer, key );
1047 if( lut_def == NULL )
1048 lut_def = msLayerGetProcessingKey( layer, "LUT" );
1049 if( lut_def == NULL )
1050 return 0;
1051
1052 /* -------------------------------------------------------------------- */
1053 /* Does this look like a file? If so, read it into memory. */
1054 /* -------------------------------------------------------------------- */
1055 if( strspn(lut_def,"0123456789:, ") != strlen(lut_def) ) {
1056 FILE *fp;
1057 char path[MS_MAXPATHLEN];
1058 int len;
1059 msBuildPath(path, layer->map->mappath, lut_def);
1060 fp = fopen( path, "rb" );
1061 if( fp == NULL ) {
1062 msSetError(MS_IOERR,
1063 "Failed to open LUT file '%s'.",
1064 "drawGDAL()", path );
1065 return -1;
1066 }
1067
1068 len = fread( lut_def_fromfile, 1, sizeof(lut_def_fromfile), fp );
1069 fclose( fp );
1070
1071 if( len == sizeof(lut_def_fromfile) ) {
1072 msSetError(MS_IOERR,
1073 "LUT definition from file %s longer than maximum buffer size (%d bytes).",
1074 "drawGDAL()",
1075 path, (int)sizeof(lut_def_fromfile) );
1076 return -1;
1077 }
1078
1079 lut_def_fromfile[len] = '\0';
1080
1081 lut_def = lut_def_fromfile;
1082 }
1083
1084 *ppszLutDef = msStrdup(lut_def);
1085
1086 return 0;
1087
1088 }
1089
1090 /************************************************************************/
1091 /* FreeLUTs() */
1092 /************************************************************************/
1093
FreeLUTs(char ** apszLUTs)1094 static void FreeLUTs(char** apszLUTs)
1095 {
1096 int i;
1097 for( i = 0; i < 4; i++ )
1098 msFree(apszLUTs[i]);
1099 msFree(apszLUTs);
1100 }
1101
1102 /************************************************************************/
1103 /* LoadLUTs() */
1104 /* */
1105 /* Return an array of 4 strings (some possibly NULL) with loaded LUTs */
1106 /* or NULL in case of failure. */
1107 /************************************************************************/
1108
LoadLUTs(layerObj * layer,int band_count)1109 static char** LoadLUTs(layerObj *layer, int band_count)
1110 {
1111 int i;
1112 char** apszLUTs;
1113
1114 assert(band_count <= 4);
1115
1116 apszLUTs = (char**) msSmallCalloc( 4, sizeof(char*) );
1117 for( i = 0; i < band_count; i++ )
1118 {
1119 if( LoadLUT( layer, i+1, &apszLUTs[i] ) != 0 )
1120 {
1121 FreeLUTs(apszLUTs);
1122 return NULL;
1123 }
1124 }
1125
1126 return apszLUTs;
1127 }
1128
1129 /************************************************************************/
1130 /* GetDataTypeAppropriateForLUTS() */
1131 /* */
1132 /* This does a quick examination of the LUT strings to determine */
1133 /* if they have input values > 255, in which case the raster data */
1134 /* must be queries on 16-bits. */
1135 /************************************************************************/
1136
GetDataTypeAppropriateForLUTS(char ** apszLUTs)1137 static GDALDataType GetDataTypeAppropriateForLUTS(char** apszLUTs)
1138 {
1139 GDALDataType eDT = GDT_Byte;
1140 int i;
1141 for( i = 0; i < 4; i++ )
1142 {
1143 const char* pszLastTuple;
1144 int nLastInValue;
1145 if( apszLUTs[i] == NULL )
1146 continue;
1147 if( EQUALN(apszLUTs[i],"# GIMP",6) )
1148 continue;
1149 /* Find last in:out tuple in string */
1150 pszLastTuple = strrchr( apszLUTs[i], ',' );
1151 if( pszLastTuple == NULL )
1152 pszLastTuple = apszLUTs[i];
1153 else
1154 pszLastTuple ++;
1155 while( *pszLastTuple == ' ' )
1156 pszLastTuple ++;
1157 nLastInValue = atoi(pszLastTuple);
1158 if( nLastInValue > 255 )
1159 {
1160 eDT = GDT_UInt16;
1161 break;
1162 }
1163 }
1164 return eDT;
1165 }
1166
1167 /************************************************************************/
1168 /* ApplyLUT() */
1169 /************************************************************************/
1170
ApplyLUT(int iColorIndex,const char * lut_def,const void * pInBuffer,GDALDataType eDT,GByte * pabyOutBuffer,int nPixelCount)1171 static int ApplyLUT( int iColorIndex, const char* lut_def,
1172 const void* pInBuffer, GDALDataType eDT,
1173 GByte* pabyOutBuffer,
1174 int nPixelCount )
1175 {
1176 int i, err;
1177 GByte byteLut[256];
1178 GByte* uint16Lut = NULL;
1179
1180 assert( eDT == GDT_Byte || eDT == GDT_UInt16 );
1181
1182 if( lut_def == NULL )
1183 {
1184 if( pInBuffer != pabyOutBuffer )
1185 {
1186 GDALCopyWords( (void*)pInBuffer, eDT, GDALGetDataTypeSize(eDT) / 8,
1187 pabyOutBuffer, GDT_Byte, 1,
1188 nPixelCount );
1189 }
1190 return 0;
1191 }
1192
1193 /* -------------------------------------------------------------------- */
1194 /* Parse the LUT description. */
1195 /* -------------------------------------------------------------------- */
1196 if( EQUALN(lut_def,"# GIMP",6) ) {
1197 if( eDT != GDT_Byte ) {
1198 msSetError(MS_MISCERR,
1199 "Cannot apply a GIMP LUT on a 16-bit buffer",
1200 "ApplyLUT()");
1201 return -1;
1202 }
1203 err = ParseGimpLUT( lut_def, byteLut, iColorIndex );
1204 } else {
1205 if( eDT == GDT_Byte )
1206 err = ParseDefaultLUT( lut_def, byteLut, 255 );
1207 else
1208 {
1209 uint16Lut = (GByte*) malloc( 65536 );
1210 if( uint16Lut == NULL )
1211 {
1212 msSetError(MS_MEMERR,
1213 "Cannot allocate 16-bit LUT",
1214 "ApplyLUT()" );
1215 return -1;
1216 }
1217 err = ParseDefaultLUT( lut_def, uint16Lut, 65535 );
1218 }
1219 }
1220
1221 /* -------------------------------------------------------------------- */
1222 /* Apply LUT. */
1223 /* -------------------------------------------------------------------- */
1224 if( eDT == GDT_Byte )
1225 {
1226 for( i = 0; i < nPixelCount; i++ )
1227 pabyOutBuffer[i] = byteLut[((GByte*)pInBuffer)[i]];
1228 }
1229 else
1230 {
1231 for( i = 0; i < nPixelCount; i++ )
1232 pabyOutBuffer[i] = uint16Lut[((GUInt16*)pInBuffer)[i]];
1233 free(uint16Lut);
1234 }
1235
1236 return err;
1237 }
1238
1239 /************************************************************************/
1240 /* LoadGDALImages() */
1241 /* */
1242 /* This call will load and process 1-4 bands of input for the */
1243 /* selected rectangle, loading the result into the passed 8bit */
1244 /* buffer. The processing options include scaling. */
1245 /************************************************************************/
1246
1247 static int
LoadGDALImages(GDALDatasetH hDS,int band_numbers[4],int band_count,layerObj * layer,int src_xoff,int src_yoff,int src_xsize,int src_ysize,GByte * pabyWholeBuffer,int dst_xsize,int dst_ysize,int * pbHaveRGBNoData,int * pnNoData1,int * pnNoData2,int * pnNoData3)1248 LoadGDALImages( GDALDatasetH hDS, int band_numbers[4], int band_count,
1249 layerObj *layer,
1250 int src_xoff, int src_yoff, int src_xsize, int src_ysize,
1251 GByte *pabyWholeBuffer,
1252 int dst_xsize, int dst_ysize,
1253 int *pbHaveRGBNoData,
1254 int *pnNoData1, int *pnNoData2, int *pnNoData3 )
1255
1256 {
1257 int iColorIndex, result_code=0;
1258 CPLErr eErr;
1259 float *pafWholeRawData;
1260 char** papszLUTs;
1261
1262 /* -------------------------------------------------------------------- */
1263 /* If we have no alpha band, but we do have three input */
1264 /* bands, then check for nodata values. If we only have one */
1265 /* input band, then nodata will already have been adderssed as */
1266 /* part of the real or manufactured color table. */
1267 /* -------------------------------------------------------------------- */
1268 if( band_count == 3 ) {
1269 *pnNoData1 = (int)
1270 msGetGDALNoDataValue( layer,
1271 GDALGetRasterBand(hDS,band_numbers[0]),
1272 pbHaveRGBNoData);
1273
1274 if( *pbHaveRGBNoData )
1275 *pnNoData2 = (int)
1276 msGetGDALNoDataValue( layer,
1277 GDALGetRasterBand(hDS,band_numbers[1]),
1278 pbHaveRGBNoData);
1279 if( *pbHaveRGBNoData )
1280 *pnNoData3 = (int)
1281 msGetGDALNoDataValue( layer,
1282 GDALGetRasterBand(hDS,band_numbers[2]),
1283 pbHaveRGBNoData);
1284 }
1285
1286 /* -------------------------------------------------------------------- */
1287 /* Are we doing a simple, non-scaling case? If so, read directly */
1288 /* and return. */
1289 /* -------------------------------------------------------------------- */
1290 if( CSLFetchNameValue( layer->processing, "SCALE" ) == NULL
1291 && CSLFetchNameValue( layer->processing, "SCALE_1" ) == NULL
1292 && CSLFetchNameValue( layer->processing, "SCALE_2" ) == NULL
1293 && CSLFetchNameValue( layer->processing, "SCALE_3" ) == NULL
1294 && CSLFetchNameValue( layer->processing, "SCALE_4" ) == NULL ) {
1295
1296 GDALDataType eDT;
1297 GUInt16* panBuffer = NULL;
1298 void* pBuffer;
1299
1300 papszLUTs = LoadLUTs(layer, band_count);
1301 if( papszLUTs == NULL ) {
1302 return -1;
1303 }
1304
1305 eDT = GetDataTypeAppropriateForLUTS(papszLUTs);
1306 if( eDT == GDT_UInt16 ) {
1307 panBuffer =
1308 (GUInt16 *) malloc(sizeof(GUInt16) * dst_xsize * dst_ysize * band_count );
1309
1310 if( panBuffer == NULL ) {
1311 msSetError(MS_MEMERR,
1312 "Allocating work uint16 image of size %dx%dx%d failed.",
1313 "msDrawRasterLayerGDAL()",
1314 dst_xsize, dst_ysize, band_count );
1315 FreeLUTs(papszLUTs);
1316 return -1;
1317 }
1318 pBuffer = panBuffer;
1319 }
1320 else
1321 pBuffer = pabyWholeBuffer;
1322
1323 eErr = GDALDatasetRasterIO( hDS, GF_Read,
1324 src_xoff, src_yoff, src_xsize, src_ysize,
1325 pBuffer,
1326 dst_xsize, dst_ysize, eDT,
1327 band_count, band_numbers,
1328 0,0,0);
1329
1330 if( eErr != CE_None ) {
1331 msSetError( MS_IOERR,
1332 "GDALDatasetRasterIO() failed: %s",
1333 "drawGDAL()",
1334 CPLGetLastErrorMsg() );
1335 FreeLUTs(papszLUTs);
1336 msFree(panBuffer);
1337 return -1;
1338 }
1339
1340 for( iColorIndex = 0;
1341 iColorIndex < band_count && result_code == 0; iColorIndex++ ) {
1342 result_code = ApplyLUT( iColorIndex+1,
1343 papszLUTs[iColorIndex],
1344 (GByte*)pBuffer
1345 + dst_xsize*dst_ysize*iColorIndex*(GDALGetDataTypeSize(eDT)/8),
1346 eDT,
1347 pabyWholeBuffer + dst_xsize*dst_ysize*iColorIndex,
1348 dst_xsize * dst_ysize );
1349 }
1350
1351 FreeLUTs(papszLUTs);
1352 msFree(panBuffer);
1353 return result_code;
1354 }
1355
1356 /* -------------------------------------------------------------------- */
1357 /* Disable use of nodata if we are doing scaling. */
1358 /* -------------------------------------------------------------------- */
1359 *pbHaveRGBNoData = FALSE;
1360
1361 /* -------------------------------------------------------------------- */
1362 /* We need to do some scaling. Will load into either a 16bit */
1363 /* unsigned or a floating point buffer depending on the source */
1364 /* data. We offer a special case for 16U data because it is */
1365 /* common and it is a substantial win to avoid alot of floating */
1366 /* point operations on it. */
1367 /* -------------------------------------------------------------------- */
1368 /* TODO */
1369
1370 /* -------------------------------------------------------------------- */
1371 /* Allocate the raw imagery buffer, and load into it (band */
1372 /* interleaved). */
1373 /* -------------------------------------------------------------------- */
1374 pafWholeRawData =
1375 (float *) malloc(sizeof(float) * dst_xsize * dst_ysize * band_count );
1376
1377 if( pafWholeRawData == NULL ) {
1378 msSetError(MS_MEMERR,
1379 "Allocating work float image of size %dx%dx%d failed.",
1380 "msDrawRasterLayerGDAL()",
1381 dst_xsize, dst_ysize, band_count );
1382 return -1;
1383 }
1384
1385 eErr = GDALDatasetRasterIO(
1386 hDS, GF_Read,
1387 src_xoff, src_yoff, src_xsize, src_ysize,
1388 pafWholeRawData, dst_xsize, dst_ysize, GDT_Float32,
1389 band_count, band_numbers,
1390 0, 0, 0 );
1391
1392 if( eErr != CE_None ) {
1393 msSetError( MS_IOERR, "GDALDatasetRasterIO() failed: %s",
1394 "drawGDAL()",
1395 CPLGetLastErrorMsg() );
1396
1397 free( pafWholeRawData );
1398 return -1;
1399 }
1400
1401 papszLUTs = LoadLUTs(layer, band_count);
1402 if( papszLUTs == NULL ) {
1403 free( pafWholeRawData );
1404 return -1;
1405 }
1406
1407 if( GetDataTypeAppropriateForLUTS(papszLUTs) != GDT_Byte ) {
1408 msDebug( "LoadGDALImage(%s): One of the LUT contains a input value > 255.\n"
1409 "This is not properly supported in combination with SCALE\n",
1410 layer->name );
1411 }
1412
1413 /* -------------------------------------------------------------------- */
1414 /* Fetch the scale processing option. */
1415 /* -------------------------------------------------------------------- */
1416 for( iColorIndex = 0; iColorIndex < band_count; iColorIndex++ ) {
1417 unsigned char *pabyBuffer;
1418 double dfScaleMin=0.0, dfScaleMax=255.0, dfScaleRatio, dfNoDataValue;
1419 const char *pszScaleInfo;
1420 float *pafRawData;
1421 int nPixelCount = dst_xsize * dst_ysize, i, bGotNoData = FALSE;
1422 GDALRasterBandH hBand =GDALGetRasterBand(hDS,band_numbers[iColorIndex]);
1423 pszScaleInfo = CSLFetchNameValue( layer->processing, "SCALE" );
1424 if( pszScaleInfo == NULL ) {
1425 char szBandScalingName[20];
1426
1427 sprintf( szBandScalingName, "SCALE_%d", iColorIndex+1 );
1428 pszScaleInfo = CSLFetchNameValue( layer->processing,
1429 szBandScalingName);
1430 }
1431
1432 if( pszScaleInfo != NULL ) {
1433 char **papszTokens;
1434
1435 papszTokens = CSLTokenizeStringComplex( pszScaleInfo, " ,",
1436 FALSE, FALSE );
1437 if( CSLCount(papszTokens) == 1
1438 && EQUAL(papszTokens[0],"AUTO") ) {
1439 dfScaleMin = dfScaleMax = 0.0;
1440 } else if( CSLCount(papszTokens) != 2 ) {
1441 free( pafWholeRawData );
1442 CSLDestroy( papszTokens );
1443 FreeLUTs( papszLUTs );
1444 msSetError( MS_MISCERR,
1445 "SCALE PROCESSING option unparsable for layer %s.",
1446 "msDrawGDAL()",
1447 layer->name );
1448 return -1;
1449 } else {
1450 dfScaleMin = atof(papszTokens[0]);
1451 dfScaleMax = atof(papszTokens[1]);
1452 }
1453 CSLDestroy( papszTokens );
1454 }
1455
1456 /* -------------------------------------------------------------------- */
1457 /* If we are using autoscaling, then compute the max and min */
1458 /* now. Perhaps we should eventually honour the offsite value */
1459 /* as a nodata value, or get it from GDAL. */
1460 /* -------------------------------------------------------------------- */
1461 pafRawData = pafWholeRawData + iColorIndex * dst_xsize * dst_ysize;
1462
1463 dfNoDataValue = msGetGDALNoDataValue( layer, hBand, &bGotNoData );
1464
1465 if( dfScaleMin == dfScaleMax ) {
1466 int bMinMaxSet = 0;
1467
1468 /* we force assignment to a float rather than letting pafRawData[i]
1469 get promoted to double later to avoid float precision issues. */
1470 float fNoDataValue = (float) dfNoDataValue;
1471
1472 for( i = 0; i < nPixelCount; i++ ) {
1473 if( bGotNoData && pafRawData[i] == fNoDataValue )
1474 continue;
1475
1476 if( CPLIsNan(pafRawData[i]) )
1477 continue;
1478
1479 if( !bMinMaxSet ) {
1480 dfScaleMin = dfScaleMax = pafRawData[i];
1481 bMinMaxSet = TRUE;
1482 }
1483
1484 dfScaleMin = MS_MIN(dfScaleMin,pafRawData[i]);
1485 dfScaleMax = MS_MAX(dfScaleMax,pafRawData[i]);
1486 }
1487
1488 if( dfScaleMin == dfScaleMax )
1489 dfScaleMax = dfScaleMin + 1.0;
1490 }
1491
1492 if( layer->debug > 0 )
1493 msDebug( "msDrawGDAL(%s): scaling to 8bit, src range=%g,%g\n",
1494 layer->name, dfScaleMin, dfScaleMax );
1495
1496 /* -------------------------------------------------------------------- */
1497 /* Now process the data. */
1498 /* -------------------------------------------------------------------- */
1499 dfScaleRatio = 256.0 / (dfScaleMax - dfScaleMin);
1500 pabyBuffer = pabyWholeBuffer + iColorIndex * nPixelCount;
1501
1502 for( i = 0; i < nPixelCount; i++ ) {
1503 float fScaledValue = (float) ((pafRawData[i]-dfScaleMin)*dfScaleRatio);
1504
1505 if( fScaledValue < 0.0 )
1506 pabyBuffer[i] = 0;
1507 else if( fScaledValue > 255.0 )
1508 pabyBuffer[i] = 255;
1509 else
1510 pabyBuffer[i] = (int) fScaledValue;
1511 }
1512
1513 /* -------------------------------------------------------------------- */
1514 /* Report a warning if NODATA keyword was applied. We are */
1515 /* unable to utilize it since we can't return any pixels marked */
1516 /* as nodata from this function. Need to fix someday. */
1517 /* -------------------------------------------------------------------- */
1518 if( bGotNoData )
1519 msDebug( "LoadGDALImage(%s): NODATA value %g in GDAL\n"
1520 "file or PROCESSING directive largely ignored. Not yet fully supported for\n"
1521 "unclassified scaled data. The NODATA value is excluded from auto-scaling\n"
1522 "min/max computation, but will not be transparent.\n",
1523 layer->name, dfNoDataValue );
1524
1525 /* -------------------------------------------------------------------- */
1526 /* Apply LUT if there is one. */
1527 /* -------------------------------------------------------------------- */
1528 result_code = ApplyLUT( iColorIndex+1,
1529 papszLUTs[iColorIndex],
1530 pabyBuffer,
1531 GDT_Byte,
1532 pabyBuffer,
1533 dst_xsize * dst_ysize );
1534 if( result_code == -1 ) {
1535 free( pafWholeRawData );
1536 FreeLUTs( papszLUTs );
1537 return result_code;
1538 }
1539 }
1540
1541 free( pafWholeRawData );
1542 FreeLUTs( papszLUTs );
1543
1544 return result_code;
1545 }
1546
1547
1548 /************************************************************************/
1549 /* msGetGDALGeoTransform() */
1550 /* */
1551 /* Cover function that tries GDALGetGeoTransform(), a world */
1552 /* file or OWS extents. It is assumed that TLOCK_GDAL is held */
1553 /* before this function is called. */
1554 /************************************************************************/
1555
msGetGDALGeoTransform(GDALDatasetH hDS,mapObj * map,layerObj * layer,double * padfGeoTransform)1556 int msGetGDALGeoTransform( GDALDatasetH hDS, mapObj *map, layerObj *layer,
1557 double *padfGeoTransform )
1558
1559 {
1560 const char *extent_priority = NULL;
1561 const char *value;
1562 const char *gdalDesc;
1563 const char *fullPath;
1564 char *fileExtension = NULL;
1565 char szPath[MS_MAXPATHLEN];
1566 int fullPathLen;
1567 int success = 0;
1568 rectObj rect;
1569
1570 /* -------------------------------------------------------------------- */
1571 /* some GDAL drivers (ie. GIF) don't set geotransform on failure. */
1572 /* -------------------------------------------------------------------- */
1573 padfGeoTransform[0] = 0.0;
1574 padfGeoTransform[1] = 1.0;
1575 padfGeoTransform[2] = 0.0;
1576 padfGeoTransform[3] = GDALGetRasterYSize(hDS);
1577 padfGeoTransform[4] = 0.0;
1578 padfGeoTransform[5] = -1.0;
1579
1580 /* -------------------------------------------------------------------- */
1581 /* Do we want to override GDAL with a worldfile if one is present? */
1582 /* -------------------------------------------------------------------- */
1583 extent_priority = CSLFetchNameValue( layer->processing,
1584 "EXTENT_PRIORITY" );
1585
1586 if( extent_priority != NULL
1587 && EQUALN(extent_priority,"WORLD",5) ) {
1588 /* is there a user configured worldfile? */
1589 value = CSLFetchNameValue( layer->processing, "WORLDFILE" );
1590
1591 if( value != NULL ) {
1592 /* at this point, fullPath may be a filePath or path only */
1593 fullPath = msBuildPath(szPath, map->mappath, value);
1594
1595 /* fullPath is a path only, so let's append the dataset filename */
1596 if( strrchr(szPath,'.') == NULL ) {
1597 fullPathLen = strlen(fullPath);
1598 gdalDesc = msStripPath((char*)GDALGetDescription(hDS));
1599
1600 if ( (fullPathLen + strlen(gdalDesc)) < MS_MAXPATHLEN ) {
1601 strcpy((char*)(fullPath + fullPathLen), gdalDesc);
1602 } else {
1603 msDebug("Path length to alternative worldfile exceeds MS_MAXPATHLEN.\n");
1604 fullPath = NULL;
1605 }
1606 }
1607 /* fullPath has a filename included, so get the extension */
1608 else {
1609 fileExtension = msStrdup(strrchr(szPath,'.') + 1);
1610 }
1611 }
1612 /* common behaviour with worldfile generated from basename + .wld */
1613 else {
1614 fullPath = GDALGetDescription(hDS);
1615 fileExtension = msStrdup("wld");
1616 }
1617
1618 if( fullPath )
1619 success = GDALReadWorldFile(fullPath, fileExtension, padfGeoTransform);
1620
1621 if( success && layer->debug >= MS_DEBUGLEVEL_VVV ) {
1622 msDebug("Worldfile location: %s.\n", fullPath);
1623 } else if( layer->debug >= MS_DEBUGLEVEL_VVV ) {
1624 msDebug("Failed using worldfile location: %s.\n", fullPath);
1625 }
1626
1627 msFree(fileExtension);
1628
1629 if( success )
1630 return MS_SUCCESS;
1631 }
1632
1633 /* -------------------------------------------------------------------- */
1634 /* Try GDAL. */
1635 /* */
1636 /* Make sure that ymax is always at the top, and ymin at the */
1637 /* bottom ... that is flip any files without the usual */
1638 /* orientation. This is intended to enable display of "raw" */
1639 /* files with no coordinate system otherwise they break down in */
1640 /* many ways. */
1641 /* -------------------------------------------------------------------- */
1642 if (GDALGetGeoTransform( hDS, padfGeoTransform ) == CE_None ) {
1643 if( padfGeoTransform[5] == 1.0 && padfGeoTransform[3] == 0.0 ) {
1644 padfGeoTransform[5] = -1.0;
1645 padfGeoTransform[3] = GDALGetRasterYSize(hDS);
1646 }
1647
1648 return MS_SUCCESS;
1649 }
1650
1651 /* -------------------------------------------------------------------- */
1652 /* Try worldfile. */
1653 /* -------------------------------------------------------------------- */
1654 if( GDALGetDescription(hDS) != NULL
1655 && GDALReadWorldFile(GDALGetDescription(hDS), "wld",
1656 padfGeoTransform) ) {
1657 return MS_SUCCESS;
1658 }
1659
1660 /* -------------------------------------------------------------------- */
1661 /* Do we have the extent keyword on the layer? We only use */
1662 /* this if we have a single raster associated with the layer as */
1663 /* opposed to a tile index. */
1664 /* */
1665 /* Arguably this ought to take priority over all the other */
1666 /* stuff. */
1667 /* -------------------------------------------------------------------- */
1668 if (MS_VALID_EXTENT(layer->extent) && layer->data != NULL) {
1669 rect = layer->extent;
1670
1671 padfGeoTransform[0] = rect.minx;
1672 padfGeoTransform[1] = (rect.maxx - rect.minx) /
1673 (double) GDALGetRasterXSize( hDS );
1674 padfGeoTransform[2] = 0;
1675 padfGeoTransform[3] = rect.maxy;
1676 padfGeoTransform[4] = 0;
1677 padfGeoTransform[5] = (rect.miny - rect.maxy) /
1678 (double) GDALGetRasterYSize( hDS );
1679
1680 return MS_SUCCESS;
1681 }
1682
1683 /* -------------------------------------------------------------------- */
1684 /* We didn't find any info ... use the default. */
1685 /* Reset our default geotransform. GDALGetGeoTransform() may */
1686 /* have altered it even if GDALGetGeoTransform() failed. */
1687 /* -------------------------------------------------------------------- */
1688 padfGeoTransform[0] = 0.0;
1689 padfGeoTransform[1] = 1.0;
1690 padfGeoTransform[2] = 0.0;
1691 padfGeoTransform[3] = GDALGetRasterYSize(hDS);
1692 padfGeoTransform[4] = 0.0;
1693 padfGeoTransform[5] = -1.0;
1694
1695 return MS_FAILURE;
1696 }
1697
1698 /************************************************************************/
1699 /* msDrawRasterLayerGDAL_RawMode() */
1700 /* */
1701 /* Handle the loading and application of data in rawmode. */
1702 /************************************************************************/
1703
1704 static int
msDrawRasterLayerGDAL_RawMode(mapObj * map,layerObj * layer,imageObj * image,GDALDatasetH hDS,int src_xoff,int src_yoff,int src_xsize,int src_ysize,int dst_xoff,int dst_yoff,int dst_xsize,int dst_ysize)1705 msDrawRasterLayerGDAL_RawMode(
1706 mapObj *map, layerObj *layer, imageObj *image, GDALDatasetH hDS,
1707 int src_xoff, int src_yoff, int src_xsize, int src_ysize,
1708 int dst_xoff, int dst_yoff, int dst_xsize, int dst_ysize )
1709
1710 {
1711 void *pBuffer;
1712 GDALDataType eDataType;
1713 int *band_list, band_count;
1714 int i, j, k, band;
1715 CPLErr eErr;
1716 float *f_nodatas = NULL;
1717 unsigned char *b_nodatas = NULL;
1718 GInt16 *i_nodatas = NULL;
1719 int got_nodata=FALSE;
1720 rasterBufferObj *mask_rb = NULL;
1721 rasterBufferObj s_mask_rb;
1722 if(layer->mask) {
1723 int ret;
1724 layerObj *maskLayer = GET_LAYER(map, msGetLayerIndex(map,layer->mask));
1725 mask_rb = &s_mask_rb;
1726 memset(mask_rb, 0, sizeof(s_mask_rb));
1727 ret = MS_IMAGE_RENDERER(maskLayer->maskimage)->getRasterBufferHandle(maskLayer->maskimage,mask_rb);
1728 if(ret != MS_SUCCESS)
1729 return -1;
1730 }
1731
1732 if( image->format->bands > 256 ) {
1733 msSetError( MS_IMGERR, "Too many bands (more than 256).",
1734 "msDrawRasterLayerGDAL_RawMode()" );
1735 return -1;
1736 }
1737
1738 /* -------------------------------------------------------------------- */
1739 /* What data type do we need? */
1740 /* -------------------------------------------------------------------- */
1741 if( image->format->imagemode == MS_IMAGEMODE_INT16 )
1742 eDataType = GDT_Int16;
1743 else if( image->format->imagemode == MS_IMAGEMODE_FLOAT32 )
1744 eDataType = GDT_Float32;
1745 else if( image->format->imagemode == MS_IMAGEMODE_BYTE )
1746 eDataType = GDT_Byte;
1747 else
1748 return -1;
1749
1750 /* -------------------------------------------------------------------- */
1751 /* Work out the band list. */
1752 /* -------------------------------------------------------------------- */
1753 band_list = msGetGDALBandList( layer, hDS, image->format->bands,
1754 &band_count );
1755 if( band_list == NULL )
1756 return -1;
1757
1758 if( band_count != image->format->bands ) {
1759 free( band_list );
1760 msSetError( MS_IMGERR, "BANDS PROCESSING directive has wrong number of bands, expected %d, got %d.",
1761 "msDrawRasterLayerGDAL_RawMode()",
1762 image->format->bands, band_count );
1763 return -1;
1764 }
1765
1766 /* -------------------------------------------------------------------- */
1767 /* Do we have nodata values? */
1768 /* -------------------------------------------------------------------- */
1769 f_nodatas = (float *) calloc(sizeof(float),band_count);
1770 if (f_nodatas == NULL) {
1771 msSetError(MS_MEMERR, "%s: %d: Out of memory allocating %u bytes.\n", "msDrawRasterLayerGDAL_RawMode()",
1772 __FILE__, __LINE__, (unsigned int)(sizeof(float)*band_count));
1773 free( band_list );
1774 return -1;
1775 }
1776
1777 if( band_count <= 3
1778 && (layer->offsite.red != -1
1779 || layer->offsite.green != -1
1780 || layer->offsite.blue != -1) ) {
1781 if( band_count > 0 )
1782 f_nodatas[0] = layer->offsite.red;
1783 if( band_count > 1 )
1784 f_nodatas[1] = layer->offsite.green;
1785 if( band_count > 2 )
1786 f_nodatas[2] = layer->offsite.blue;
1787 got_nodata = TRUE;
1788 }
1789
1790 if( !got_nodata ) {
1791 got_nodata = TRUE;
1792 for( band = 0; band < band_count && got_nodata; band++ ) {
1793 f_nodatas[band] = msGetGDALNoDataValue(
1794 layer, GDALGetRasterBand(hDS,band_list[band]), &got_nodata );
1795 }
1796 }
1797
1798 if( !got_nodata ) {
1799 free( f_nodatas );
1800 f_nodatas = NULL;
1801 } else if( eDataType == GDT_Byte ) {
1802 b_nodatas = (unsigned char *) f_nodatas;
1803 for( band = 0; band < band_count; band++ )
1804 b_nodatas[band] = (unsigned char) f_nodatas[band];
1805 } else if( eDataType == GDT_Int16 ) {
1806 i_nodatas = (GInt16 *) f_nodatas;
1807 for( band = 0; band < band_count; band++ )
1808 i_nodatas[band] = (GInt16) f_nodatas[band];
1809 }
1810
1811 /* -------------------------------------------------------------------- */
1812 /* Allocate buffer, and read data into it. */
1813 /* -------------------------------------------------------------------- */
1814 pBuffer = malloc(dst_xsize * dst_ysize * image->format->bands
1815 * (GDALGetDataTypeSize(eDataType)/8) );
1816 if( pBuffer == NULL ) {
1817 msSetError(MS_MEMERR,
1818 "Allocating work image of size %dx%d failed.",
1819 "msDrawRasterLayerGDAL()", dst_xsize, dst_ysize );
1820 free( band_list );
1821 free( f_nodatas );
1822 return -1;
1823 }
1824
1825 eErr = GDALDatasetRasterIO( hDS, GF_Read,
1826 src_xoff, src_yoff, src_xsize, src_ysize,
1827 pBuffer, dst_xsize, dst_ysize, eDataType,
1828 image->format->bands, band_list,
1829 0, 0, 0 );
1830 free( band_list );
1831
1832 if( eErr != CE_None ) {
1833 msSetError( MS_IOERR, "GDALRasterIO() failed: %s",
1834 "msDrawRasterLayerGDAL_RawMode()", CPLGetLastErrorMsg() );
1835 free( pBuffer );
1836 free( f_nodatas );
1837 return -1;
1838 }
1839
1840 /* -------------------------------------------------------------------- */
1841 /* Transfer the data to the imageObj. */
1842 /* -------------------------------------------------------------------- */
1843 k = 0;
1844 for( band = 0; band < image->format->bands; band++ ) {
1845 for( i = dst_yoff; i < dst_yoff + dst_ysize; i++ ) {
1846 if( image->format->imagemode == MS_IMAGEMODE_INT16 ) {
1847 for( j = dst_xoff; j < dst_xoff + dst_xsize; j++ ) {
1848 int off = j + i * image->width
1849 + band*image->width*image->height;
1850 int off_mask = j + i * image->width;
1851
1852 if( ( i_nodatas && ((GInt16 *) pBuffer)[k] == i_nodatas[band] )
1853 || SKIP_MASK(j,i)) {
1854 k++;
1855 continue;
1856 }
1857
1858 image->img.raw_16bit[off] = ((GInt16 *) pBuffer)[k++];
1859 MS_SET_BIT(image->img_mask,off_mask);
1860 }
1861 } else if( image->format->imagemode == MS_IMAGEMODE_FLOAT32 ) {
1862 for( j = dst_xoff; j < dst_xoff + dst_xsize; j++ ) {
1863 int off = j + i * image->width
1864 + band*image->width*image->height;
1865 int off_mask = j + i * image->width;
1866
1867 if( ( f_nodatas && ((float *) pBuffer)[k] == f_nodatas[band] )
1868 || SKIP_MASK(j,i)) {
1869 k++;
1870 continue;
1871 }
1872
1873 image->img.raw_float[off] = ((float *) pBuffer)[k++];
1874 MS_SET_BIT(image->img_mask,off_mask);
1875 }
1876 } else if( image->format->imagemode == MS_IMAGEMODE_BYTE ) {
1877 for( j = dst_xoff; j < dst_xoff + dst_xsize; j++ ) {
1878 int off = j + i * image->width
1879 + band*image->width*image->height;
1880 int off_mask = j + i * image->width;
1881
1882 if( ( b_nodatas && ((unsigned char *) pBuffer)[k] == b_nodatas[band] )
1883 || SKIP_MASK(j,i)) {
1884 k++;
1885 continue;
1886 }
1887
1888 image->img.raw_byte[off] = ((unsigned char *) pBuffer)[k++];
1889 MS_SET_BIT(image->img_mask,off_mask);
1890 }
1891 }
1892 }
1893 }
1894
1895 free( pBuffer );
1896 free( f_nodatas );
1897
1898 return 0;
1899 }
1900
1901 /************************************************************************/
1902 /* msDrawRasterLayerGDAL_16BitClassifcation() */
1903 /* */
1904 /* Handle the rendering of rasters going through a 16bit */
1905 /* classification lookup instead of the more common 8bit one. */
1906 /* */
1907 /* Currently we are using one data path where we load the */
1908 /* raster into a floating point buffer, and then scale to */
1909 /* 16bit. Eventually we could add optimizations for some of */
1910 /* the 16bit cases at the cost of some complication. */
1911 /************************************************************************/
1912
1913 static int
msDrawRasterLayerGDAL_16BitClassification(mapObj * map,layerObj * layer,rasterBufferObj * rb,GDALDatasetH hDS,GDALRasterBandH hBand,int src_xoff,int src_yoff,int src_xsize,int src_ysize,int dst_xoff,int dst_yoff,int dst_xsize,int dst_ysize)1914 msDrawRasterLayerGDAL_16BitClassification(
1915 mapObj *map, layerObj *layer, rasterBufferObj *rb,
1916 GDALDatasetH hDS, GDALRasterBandH hBand,
1917 int src_xoff, int src_yoff, int src_xsize, int src_ysize,
1918 int dst_xoff, int dst_yoff, int dst_xsize, int dst_ysize )
1919
1920 {
1921 float *pafRawData;
1922 double dfScaleMin=0.0, dfScaleMax=0.0, dfScaleRatio;
1923 int nPixelCount = dst_xsize * dst_ysize, i, nBucketCount=0;
1924 GDALDataType eDataType;
1925 float fDataMin=0.0, fDataMax=255.0, fNoDataValue;
1926 const char *pszScaleInfo;
1927 const char *pszBuckets;
1928 int *cmap, c, j, k, bGotNoData = FALSE, bGotFirstValue;
1929 unsigned char *rb_cmap[4];
1930 CPLErr eErr;
1931 rasterBufferObj *mask_rb = NULL;
1932 rasterBufferObj s_mask_rb;
1933 int lastC;
1934 struct mstimeval starttime={0}, endtime={0};
1935
1936 const char *pszClassifyScaled;
1937 int bClassifyScaled = FALSE;
1938
1939 if(layer->mask) {
1940 int ret;
1941 layerObj *maskLayer = GET_LAYER(map, msGetLayerIndex(map,layer->mask));
1942 mask_rb = &s_mask_rb;
1943 memset(mask_rb, 0, sizeof(s_mask_rb));
1944 ret = MS_IMAGE_RENDERER(maskLayer->maskimage)->getRasterBufferHandle(maskLayer->maskimage,mask_rb);
1945 if(ret != MS_SUCCESS)
1946 return -1;
1947 }
1948
1949 assert( rb->type == MS_BUFFER_BYTE_RGBA );
1950
1951 /* ==================================================================== */
1952 /* Read the requested data in one gulp into a floating point */
1953 /* buffer. */
1954 /* ==================================================================== */
1955 pafRawData = (float *) malloc(sizeof(float) * dst_xsize * dst_ysize );
1956 if( pafRawData == NULL ) {
1957 msSetError( MS_MEMERR, "Out of memory allocating working buffer.",
1958 "msDrawRasterLayerGDAL_16BitClassification()" );
1959 return -1;
1960 }
1961
1962 eErr = GDALRasterIO( hBand, GF_Read,
1963 src_xoff, src_yoff, src_xsize, src_ysize,
1964 pafRawData, dst_xsize, dst_ysize, GDT_Float32, 0, 0 );
1965
1966 if( eErr != CE_None ) {
1967 free( pafRawData );
1968 msSetError( MS_IOERR, "GDALRasterIO() failed: %s",
1969 "msDrawRasterLayerGDAL_16BitClassification()",
1970 CPLGetLastErrorMsg() );
1971 return -1;
1972 }
1973
1974 fNoDataValue = (float) msGetGDALNoDataValue( layer, hBand, &bGotNoData );
1975
1976 /* ==================================================================== */
1977 /* Determine scaling. */
1978 /* ==================================================================== */
1979 eDataType = GDALGetRasterDataType( hBand );
1980
1981 /* -------------------------------------------------------------------- */
1982 /* Scan for absolute min/max of this block. */
1983 /* -------------------------------------------------------------------- */
1984 bGotFirstValue = FALSE;
1985
1986 for( i = 1; i < nPixelCount; i++ ) {
1987 if( bGotNoData && pafRawData[i] == fNoDataValue )
1988 continue;
1989
1990 if( CPLIsNan(pafRawData[i]) )
1991 continue;
1992
1993 if( !bGotFirstValue ) {
1994 fDataMin = fDataMax = pafRawData[i];
1995 bGotFirstValue = TRUE;
1996 } else {
1997 fDataMin = MS_MIN(fDataMin,pafRawData[i]);
1998 fDataMax = MS_MAX(fDataMax,pafRawData[i]);
1999 }
2000 }
2001
2002 /* -------------------------------------------------------------------- */
2003 /* Fetch the scale classification option. */
2004 /* -------------------------------------------------------------------- */
2005 pszClassifyScaled = CSLFetchNameValue( layer->processing, "CLASSIFY_SCALED" );
2006 if( pszClassifyScaled != NULL )
2007 bClassifyScaled = CSLTestBoolean(pszClassifyScaled);
2008
2009 /* -------------------------------------------------------------------- */
2010 /* Fetch the scale processing option. */
2011 /* -------------------------------------------------------------------- */
2012 pszBuckets = CSLFetchNameValue( layer->processing, "SCALE_BUCKETS" );
2013 pszScaleInfo = CSLFetchNameValue( layer->processing, "SCALE" );
2014
2015 if( pszScaleInfo != NULL ) {
2016 char **papszTokens;
2017
2018 papszTokens = CSLTokenizeStringComplex( pszScaleInfo, " ,",
2019 FALSE, FALSE );
2020 if( CSLCount(papszTokens) == 1
2021 && EQUAL(papszTokens[0],"AUTO") ) {
2022 dfScaleMin = dfScaleMax = 0.0;
2023 } else if( CSLCount(papszTokens) != 2 ) {
2024 free( pafRawData );
2025 CSLDestroy( papszTokens );
2026 msSetError( MS_MISCERR,
2027 "SCALE PROCESSING option unparsable for layer %s.",
2028 "msDrawGDAL()",
2029 layer->name );
2030 return -1;
2031 } else {
2032 dfScaleMin = atof(papszTokens[0]);
2033 dfScaleMax = atof(papszTokens[1]);
2034 }
2035 CSLDestroy( papszTokens );
2036 }
2037
2038 /* -------------------------------------------------------------------- */
2039 /* Special integer cases for scaling. */
2040 /* */
2041 /* TODO: Treat Int32 and UInt32 case the same way *if* the min */
2042 /* and max are less than 65536 apart. */
2043 /* -------------------------------------------------------------------- */
2044 if( eDataType == GDT_Byte || eDataType == GDT_Int16
2045 || eDataType == GDT_UInt16 ) {
2046 if( pszScaleInfo == NULL ) {
2047 dfScaleMin = fDataMin - 0.5;
2048 dfScaleMax = fDataMax + 0.5;
2049 }
2050
2051 if( pszBuckets == NULL ) {
2052 nBucketCount = (int) floor(fDataMax - fDataMin + 1.1);
2053 }
2054 }
2055
2056 /* -------------------------------------------------------------------- */
2057 /* General case if no scaling values provided in mapfile. */
2058 /* -------------------------------------------------------------------- */
2059 else if( dfScaleMin == 0.0 && dfScaleMax == 0.0 ) {
2060 double dfEpsilon = (fDataMax - fDataMin) / (65536*2);
2061 dfScaleMin = fDataMin - dfEpsilon;
2062 dfScaleMax = fDataMax + dfEpsilon;
2063 }
2064
2065 /* -------------------------------------------------------------------- */
2066 /* How many buckets? Only set if we don't already have a value. */
2067 /* -------------------------------------------------------------------- */
2068 if( pszBuckets == NULL ) {
2069 if( nBucketCount == 0 )
2070 nBucketCount = 65536;
2071 } else {
2072 nBucketCount = atoi(pszBuckets);
2073 if( nBucketCount < 2 ) {
2074 free( pafRawData );
2075 msSetError( MS_MISCERR,
2076 "SCALE_BUCKETS PROCESSING option is not a value of 2 or more: %s.",
2077 "msDrawRasterLayerGDAL_16BitClassification()",
2078 pszBuckets );
2079 return -1;
2080 }
2081 }
2082
2083 /* -------------------------------------------------------------------- */
2084 /* Compute scaling ratio. */
2085 /* -------------------------------------------------------------------- */
2086 if( dfScaleMax == dfScaleMin )
2087 dfScaleMax = dfScaleMin + 1.0;
2088
2089 dfScaleRatio = nBucketCount / (dfScaleMax - dfScaleMin);
2090
2091 if( layer->debug > 0 )
2092 msDebug( "msDrawRasterGDAL_16BitClassification(%s):\n"
2093 " scaling to %d buckets from range=%g,%g.\n",
2094 layer->name, nBucketCount, dfScaleMin, dfScaleMax );
2095
2096 /* ==================================================================== */
2097 /* Compute classification lookup table. */
2098 /* ==================================================================== */
2099
2100 cmap = (int *) msSmallCalloc(sizeof(int),nBucketCount);
2101 rb_cmap[0] = (unsigned char *) msSmallCalloc(1,nBucketCount);
2102 rb_cmap[1] = (unsigned char *) msSmallCalloc(1,nBucketCount);
2103 rb_cmap[2] = (unsigned char *) msSmallCalloc(1,nBucketCount);
2104 rb_cmap[3] = (unsigned char *) msSmallCalloc(1,nBucketCount);
2105
2106
2107 if(layer->debug >= MS_DEBUGLEVEL_TUNING) {
2108 msGettimeofday(&starttime, NULL);
2109 }
2110
2111 lastC = -1;
2112 for(i=0; i < nBucketCount; i++) {
2113 double dfOriginalValue;
2114
2115 cmap[i] = -1;
2116
2117 // i = (int) ((dfOriginalValue - dfScaleMin) * dfScaleRatio+1)-1;
2118 dfOriginalValue = (i+0.5) / dfScaleRatio + dfScaleMin;
2119
2120 /* The creation of buckets takes a significant time when they are many, and many classes
2121 as well. When iterating over buckets, a faster strategy is to reuse first the last used
2122 class index. */
2123 if(bClassifyScaled == TRUE)
2124 c = msGetClass_FloatRGB_WithFirstClassToTry(layer, (float) i, -1, -1, -1, lastC);
2125 else
2126 c = msGetClass_FloatRGB_WithFirstClassToTry(layer, (float) dfOriginalValue, -1, -1, -1, lastC);
2127
2128 lastC = c;
2129 if( c != -1 ) {
2130 int s;
2131
2132 /* change colour based on colour range? */
2133 for(s=0; s<layer->class[c]->numstyles; s++) {
2134 if( MS_VALID_COLOR(layer->class[c]->styles[s]->mincolor)
2135 && MS_VALID_COLOR(layer->class[c]->styles[s]->maxcolor) )
2136 msValueToRange(layer->class[c]->styles[s],dfOriginalValue, MS_COLORSPACE_RGB);
2137 }
2138 if( MS_TRANSPARENT_COLOR(layer->class[c]->styles[0]->color) ) {
2139 /* leave it transparent */
2140 } else if( MS_VALID_COLOR(layer->class[c]->styles[0]->color)) {
2141 /* use class color */
2142 rb_cmap[0][i] = layer->class[c]->styles[0]->color.red;
2143 rb_cmap[1][i] = layer->class[c]->styles[0]->color.green;
2144 rb_cmap[2][i] = layer->class[c]->styles[0]->color.blue;
2145 rb_cmap[3][i] = (255*layer->class[c]->styles[0]->opacity / 100);
2146 }
2147 }
2148 }
2149
2150 if(layer->debug >= MS_DEBUGLEVEL_TUNING) {
2151 msGettimeofday(&endtime, NULL);
2152 msDebug("msDrawRasterGDAL_16BitClassification() bucket creation time: %.3fs\n",
2153 (endtime.tv_sec+endtime.tv_usec/1.0e6)-
2154 (starttime.tv_sec+starttime.tv_usec/1.0e6) );
2155 }
2156
2157 /* ==================================================================== */
2158 /* Now process the data, applying to the working imageObj. */
2159 /* ==================================================================== */
2160 k = 0;
2161
2162 for( i = dst_yoff; i < dst_yoff + dst_ysize; i++ ) {
2163 for( j = dst_xoff; j < dst_xoff + dst_xsize; j++ ) {
2164 float fRawValue = pafRawData[k++];
2165 int iMapIndex;
2166
2167 /*
2168 * Skip nodata pixels ... no processing.
2169 */
2170 if( bGotNoData && fRawValue == fNoDataValue ) {
2171 continue;
2172 }
2173
2174 if( CPLIsNan(fRawValue) )
2175 continue;
2176
2177 if(SKIP_MASK(j,i))
2178 continue;
2179
2180 /*
2181 * The funny +1/-1 is to avoid odd rounding around zero.
2182 * We could use floor() but sometimes it is expensive.
2183 */
2184 iMapIndex = (int) ((fRawValue - dfScaleMin) * dfScaleRatio+1)-1;
2185
2186 if( iMapIndex >= nBucketCount || iMapIndex < 0 ) {
2187 continue;
2188 }
2189 /* currently we never have partial alpha so keep simple */
2190 if( rb_cmap[3][iMapIndex] > 0 )
2191 RB_SET_PIXEL( rb, j, i,
2192 rb_cmap[0][iMapIndex],
2193 rb_cmap[1][iMapIndex],
2194 rb_cmap[2][iMapIndex],
2195 rb_cmap[3][iMapIndex] );
2196 }
2197 }
2198
2199 /* -------------------------------------------------------------------- */
2200 /* Cleanup */
2201 /* -------------------------------------------------------------------- */
2202 free( pafRawData );
2203 free( cmap );
2204 free( rb_cmap[0] );
2205 free( rb_cmap[1] );
2206 free( rb_cmap[2] );
2207 free( rb_cmap[3] );
2208
2209 assert( k == dst_xsize * dst_ysize );
2210
2211 return 0;
2212 }
2213
2214 /************************************************************************/
2215 /* msGetGDALNoDataValue() */
2216 /************************************************************************/
2217 double
msGetGDALNoDataValue(layerObj * layer,void * hBand,int * pbGotNoData)2218 msGetGDALNoDataValue( layerObj *layer, void *hBand, int *pbGotNoData )
2219
2220 {
2221 const char *pszNODATAOpt;
2222
2223 *pbGotNoData = FALSE;
2224
2225 /* -------------------------------------------------------------------- */
2226 /* Check for a NODATA setting in the .map file. */
2227 /* -------------------------------------------------------------------- */
2228 pszNODATAOpt = CSLFetchNameValue(layer->processing,"NODATA");
2229 if( pszNODATAOpt != NULL ) {
2230 if( EQUAL(pszNODATAOpt,"OFF") || strlen(pszNODATAOpt) == 0 )
2231 return -1234567.0;
2232 if( !EQUAL(pszNODATAOpt,"AUTO") ) {
2233 *pbGotNoData = TRUE;
2234 return atof(pszNODATAOpt);
2235 }
2236 }
2237
2238 /* -------------------------------------------------------------------- */
2239 /* Check for a NODATA setting on the raw file. */
2240 /* -------------------------------------------------------------------- */
2241 if( hBand == NULL )
2242 return -1234567.0;
2243 else
2244 return GDALGetRasterNoDataValue( hBand, pbGotNoData );
2245 }
2246
2247 /************************************************************************/
2248 /* msGetGDALBandList() */
2249 /* */
2250 /* max_bands - pass zero for no limit. */
2251 /* band_count - location in which length of the return band */
2252 /* list is placed. */
2253 /* */
2254 /* returns malloc() list of band numbers (*band_count long). */
2255 /************************************************************************/
2256
msGetGDALBandList(layerObj * layer,void * hDS,int max_bands,int * band_count)2257 int *msGetGDALBandList( layerObj *layer, void *hDS,
2258 int max_bands, int *band_count )
2259
2260 {
2261 int i, file_bands;
2262 int *band_list = NULL;
2263
2264 file_bands = GDALGetRasterCount( hDS );
2265 if( file_bands == 0 ) {
2266 msSetError( MS_IMGERR,
2267 "Attempt to operate on GDAL file with no bands, layer=%s.",
2268 "msGetGDALBandList()",
2269 layer->name );
2270
2271 return NULL;
2272 }
2273
2274 /* -------------------------------------------------------------------- */
2275 /* Use all (or first) bands. */
2276 /* -------------------------------------------------------------------- */
2277 if( CSLFetchNameValue( layer->processing, "BANDS" ) == NULL ) {
2278 if( max_bands > 0 )
2279 *band_count = MS_MIN(file_bands,max_bands);
2280 else
2281 *band_count = file_bands;
2282
2283 band_list = (int *) malloc(sizeof(int) * *band_count );
2284
2285 /* FIXME MS_CHECK_ALLOC leaks papszItems */
2286 MS_CHECK_ALLOC(band_list, sizeof(int) * *band_count, NULL);
2287
2288 for( i = 0; i < *band_count; i++ )
2289 band_list[i] = i+1;
2290 return band_list;
2291 }
2292
2293 /* -------------------------------------------------------------------- */
2294 /* get bands from BANDS processing directive. */
2295 /* -------------------------------------------------------------------- */
2296 else {
2297 char **papszItems = CSLTokenizeStringComplex(
2298 CSLFetchNameValue(layer->processing,"BANDS"), " ,", FALSE, FALSE );
2299
2300 if( CSLCount(papszItems) == 0 ) {
2301 CSLDestroy( papszItems );
2302 msSetError( MS_IMGERR, "BANDS PROCESSING directive has no items.",
2303 "msGetGDALBandList()" );
2304 return NULL;
2305 } else if( max_bands != 0 && CSLCount(papszItems) > max_bands ) {
2306 msSetError( MS_IMGERR, "BANDS PROCESSING directive has wrong number of bands, expected at most %d, got %d.",
2307 "msGetGDALBandList()",
2308 max_bands, CSLCount(papszItems) );
2309 CSLDestroy( papszItems );
2310 return NULL;
2311 }
2312
2313 *band_count = CSLCount(papszItems);
2314 band_list = (int *) malloc(sizeof(int) * *band_count);
2315 MS_CHECK_ALLOC(band_list, sizeof(int) * *band_count, NULL);
2316
2317 for( i = 0; i < *band_count; i++ ) {
2318 band_list[i] = atoi(papszItems[i]);
2319 if( band_list[i] < 1 || band_list[i] > GDALGetRasterCount(hDS) ) {
2320 msSetError( MS_IMGERR,
2321 "BANDS PROCESSING directive includes illegal band '%s', should be from 1 to %d.",
2322 "msGetGDALBandList()",
2323 papszItems[i], GDALGetRasterCount(hDS) );
2324 CSLDestroy( papszItems );
2325 free( band_list );
2326 return NULL;
2327 }
2328 }
2329 CSLDestroy( papszItems );
2330
2331 return band_list;
2332 }
2333 }
2334