/****************************************************************************** * $Id$ * * Project: MapServer * Purpose: Assorted code related to resampling rasters. * Author: Frank Warmerdam, warmerdam@pobox.com * ****************************************************************************** * Copyright (c) 1996-2005 Regents of the University of Minnesota. * * Permission is hereby granted, free of charge, to any person obtaining a * copy of this software and associated documentation files (the "Software"), * to deal in the Software without restriction, including without limitation * the rights to use, copy, modify, merge, publish, distribute, sublicense, * and/or sell copies of the Software, and to permit persons to whom the * Software is furnished to do so, subject to the following conditions: * * The above copyright notice and this permission notice shall be included in * all copies of this Software or works derived from this Software. * * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER * DEALINGS IN THE SOFTWARE. *****************************************************************************/ #include #include "mapresample.h" #include "mapthread.h" MS_CVSID("$Id$") #ifndef MAX # define MIN(a,b) ((ab) ? a : b) #endif /************************************************************************/ /* InvGeoTransform() */ /* */ /* Invert a standard 3x2 "GeoTransform" style matrix with an */ /* implicit [1 0 0] final row. */ /************************************************************************/ int InvGeoTransform( double *gt_in, double *gt_out ) { double det, inv_det; /* we assume a 3rd row that is [1 0 0] */ /* Compute determinate */ det = gt_in[1] * gt_in[5] - gt_in[2] * gt_in[4]; if( fabs(det) < 0.000000000000001 ) return 0; inv_det = 1.0 / det; /* compute adjoint, and devide by determinate */ gt_out[1] = gt_in[5] * inv_det; gt_out[4] = -gt_in[4] * inv_det; gt_out[2] = -gt_in[2] * inv_det; gt_out[5] = gt_in[1] * inv_det; gt_out[0] = ( gt_in[2] * gt_in[3] - gt_in[0] * gt_in[5]) * inv_det; gt_out[3] = (-gt_in[1] * gt_in[3] + gt_in[0] * gt_in[4]) * inv_det; return 1; } #if defined(USE_PROJ) && defined(USE_GDAL) /************************************************************************/ /* msNearestRasterResample() */ /************************************************************************/ static int msNearestRasterResampler( imageObj *psSrcImage, colorObj offsite, imageObj *psDstImage, int *panCMap, SimpleTransformer pfnTransform, void *pCBData, int debug ) { double *x, *y; int nDstX, nDstY; int *panSuccess; int nDstXSize = psDstImage->width; int nDstYSize = psDstImage->height; int nSrcXSize = psSrcImage->width; int nSrcYSize = psSrcImage->height; int nFailedPoints = 0, nSetPoints = 0; gdImagePtr srcImg, dstImg; srcImg = psSrcImage->img.gd; dstImg = psDstImage->img.gd; x = (double *) malloc( sizeof(double) * nDstXSize ); y = (double *) malloc( sizeof(double) * nDstXSize ); panSuccess = (int *) malloc( sizeof(int) * nDstXSize ); for( nDstY = 0; nDstY < nDstYSize; nDstY++ ) { for( nDstX = 0; nDstX < nDstXSize; nDstX++ ) { x[nDstX] = nDstX + 0.5; y[nDstX] = nDstY + 0.5; } pfnTransform( pCBData, nDstXSize, x, y, panSuccess ); for( nDstX = 0; nDstX < nDstXSize; nDstX++ ) { int nSrcX, nSrcY; int nValue = 0; if( !panSuccess[nDstX] ) { nFailedPoints++; continue; } nSrcX = (int) x[nDstX]; nSrcY = (int) y[nDstX]; /* * We test the original floating point values to * avoid errors related to asymmetric rounding around zero. */ if( x[nDstX] < 0.0 || y[nDstX] < 0.0 || nSrcX >= nSrcXSize || nSrcY >= nSrcYSize ) { continue; } if( MS_RENDERER_GD(psSrcImage->format) ) { if( !gdImageTrueColor(psSrcImage->img.gd) ) { nValue = panCMap[srcImg->pixels[nSrcY][nSrcX]]; if( nValue == -1 ) continue; nSetPoints++; dstImg->pixels[nDstY][nDstX] = nValue; } else { int nValue = srcImg->tpixels[nSrcY][nSrcX]; int gd_alpha = gdTrueColorGetAlpha(nValue); if( gd_alpha == 0 ) { nSetPoints++; dstImg->tpixels[nDstY][nDstX] = nValue; } else if( gd_alpha == 127 ) /* overlay is transparent, do nothing */; else { nSetPoints++; dstImg->tpixels[nDstY][nDstX] = gdAlphaBlend( dstImg->tpixels[nDstY][nDstX], nValue ); } } } else if( MS_RENDERER_RAWDATA(psSrcImage->format) ) { int band; for( band = 0; band < psSrcImage->format->bands; band++ ) { if( psSrcImage->format->imagemode == MS_IMAGEMODE_INT16 ) { int nValue; nValue = psSrcImage->img.raw_16bit[ nSrcX + nSrcY * psSrcImage->width + band*psSrcImage->width*psSrcImage->height]; if( nValue == offsite.red ) continue; nSetPoints++; psDstImage->img.raw_16bit[ nDstX + nDstY * psDstImage->width + band*psDstImage->width*psDstImage->height] = nValue; } else if( psSrcImage->format->imagemode == MS_IMAGEMODE_FLOAT32) { float fValue; fValue = psSrcImage->img.raw_float[ nSrcX + nSrcY * psSrcImage->width + band*psSrcImage->width*psSrcImage->height]; if( fValue == offsite.red ) continue; nSetPoints++; psDstImage->img.raw_float[ nDstX + nDstY * psDstImage->width + band*psDstImage->width*psDstImage->height] = fValue; } else if(psSrcImage->format->imagemode == MS_IMAGEMODE_BYTE) { nValue = psSrcImage->img.raw_byte[ nSrcX + nSrcY * psSrcImage->width + band*psSrcImage->width*psSrcImage->height]; if( nValue == offsite.red ) continue; nSetPoints++; psDstImage->img.raw_byte[ nDstX + nDstY * psDstImage->width + band*psDstImage->width*psDstImage->height] = (unsigned char) nValue; } else { assert( 0 ); } } } } } free( panSuccess ); free( x ); free( y ); /* -------------------------------------------------------------------- */ /* Some debugging output. */ /* -------------------------------------------------------------------- */ if( nFailedPoints > 0 && debug ) { char szMsg[256]; sprintf( szMsg, "msNearestRasterResampler: " "%d failed to transform, %d actually set.\n", nFailedPoints, nSetPoints ); msDebug( szMsg ); } return 0; } /************************************************************************/ /* msSourceSample() */ /************************************************************************/ static void msSourceSample( imageObj *psSrcImage, int iSrcX, int iSrcY, double *padfPixelSum, double dfWeight, double *pdfWeightSum, colorObj *offsite ) { if( MS_RENDERER_GD(psSrcImage->format) ) { if( !gdImageTrueColor(psSrcImage->img.gd) ) { padfPixelSum[0] += (dfWeight * psSrcImage->img.gd->pixels[iSrcY][iSrcX]); *pdfWeightSum += dfWeight; } else { int nValue = psSrcImage->img.gd->tpixels[iSrcY][iSrcX]; int gd_alpha = gdTrueColorGetAlpha(nValue); if( gd_alpha != 127 ) { padfPixelSum[0] += dfWeight * gdTrueColorGetRed(nValue); padfPixelSum[1] += dfWeight * gdTrueColorGetGreen(nValue); padfPixelSum[2] += dfWeight * gdTrueColorGetBlue(nValue); *pdfWeightSum += dfWeight; } } } else { int band; for( band = 0; band < psSrcImage->format->bands; band++ ) { if( psSrcImage->format->imagemode == MS_IMAGEMODE_INT16 ) { int nValue; nValue = psSrcImage->img.raw_16bit[ iSrcX + iSrcY * psSrcImage->width + band*psSrcImage->width*psSrcImage->height]; /* if band 1 is nodata, skip the rest */ if( band == 0 && nValue == offsite->red ) return; padfPixelSum[band] += dfWeight * nValue; } else if( psSrcImage->format->imagemode == MS_IMAGEMODE_FLOAT32) { float fValue; fValue = psSrcImage->img.raw_float[ iSrcX + iSrcY * psSrcImage->width + band*psSrcImage->width*psSrcImage->height]; if( band == 0 && fValue == offsite->red ) return; padfPixelSum[band] += fValue * dfWeight; } else if(psSrcImage->format->imagemode == MS_IMAGEMODE_BYTE) { int nValue; nValue = psSrcImage->img.raw_byte[ iSrcX + iSrcY * psSrcImage->width + band*psSrcImage->width*psSrcImage->height]; if( band == 0 && nValue == offsite->red ) continue; padfPixelSum[band] += nValue * dfWeight; } else { assert( 0 ); return; } } *pdfWeightSum += dfWeight; } } /************************************************************************/ /* msBilinearRasterResample() */ /************************************************************************/ static int msBilinearRasterResampler( imageObj *psSrcImage, colorObj offsite, imageObj *psDstImage, int *panCMap, SimpleTransformer pfnTransform, void *pCBData, int debug ) { double *x, *y; int nDstX, nDstY, i; int *panSuccess; int nDstXSize = psDstImage->width; int nDstYSize = psDstImage->height; int nSrcXSize = psSrcImage->width; int nSrcYSize = psSrcImage->height; int nFailedPoints = 0, nSetPoints = 0; double *padfPixelSum; gdImagePtr srcImg, dstImg; int bandCount = MAX(4,psSrcImage->format->bands); padfPixelSum = (double *) malloc(sizeof(double) * bandCount); srcImg = psSrcImage->img.gd; dstImg = psDstImage->img.gd; x = (double *) malloc( sizeof(double) * nDstXSize ); y = (double *) malloc( sizeof(double) * nDstXSize ); panSuccess = (int *) malloc( sizeof(int) * nDstXSize ); for( nDstY = 0; nDstY < nDstYSize; nDstY++ ) { for( nDstX = 0; nDstX < nDstXSize; nDstX++ ) { x[nDstX] = nDstX + 0.5; y[nDstX] = nDstY + 0.5; } pfnTransform( pCBData, nDstXSize, x, y, panSuccess ); for( nDstX = 0; nDstX < nDstXSize; nDstX++ ) { int nSrcX, nSrcY, nSrcX2, nSrcY2; double dfRatioX2, dfRatioY2, dfWeightSum = 0.0; if( !panSuccess[nDstX] ) { nFailedPoints++; continue; } /* ** Offset to treat TL pixel corners as pixel location instead ** of the center. */ x[nDstX] -= 0.5; y[nDstX] -= 0.5; nSrcX = (int) floor(x[nDstX]); nSrcY = (int) floor(y[nDstX]); nSrcX2 = nSrcX+1; nSrcY2 = nSrcY+1; dfRatioX2 = x[nDstX] - nSrcX; dfRatioY2 = y[nDstX] - nSrcY; /* If we are right off the source, skip this pixel */ if( nSrcX2 < 0 || nSrcX >= nSrcXSize || nSrcY2 < 0 || nSrcY >= nSrcYSize ) continue; /* Trim in stuff one pixel off the edge */ nSrcX = MAX(nSrcX,0); nSrcY = MAX(nSrcY,0); nSrcX2 = MIN(nSrcX2,nSrcXSize-1); nSrcY2 = MIN(nSrcY2,nSrcYSize-1); memset( padfPixelSum, 0, sizeof(double) * bandCount); msSourceSample( psSrcImage, nSrcX, nSrcY, padfPixelSum, (1.0 - dfRatioX2) * (1.0 - dfRatioY2), &dfWeightSum, &offsite ); msSourceSample( psSrcImage, nSrcX2, nSrcY, padfPixelSum, (dfRatioX2) * (1.0 - dfRatioY2), &dfWeightSum, &offsite ); msSourceSample( psSrcImage, nSrcX, nSrcY2, padfPixelSum, (1.0 - dfRatioX2) * (dfRatioY2), &dfWeightSum, &offsite ); msSourceSample( psSrcImage, nSrcX2, nSrcY2, padfPixelSum, (dfRatioX2) * (dfRatioY2), &dfWeightSum, &offsite ); if( dfWeightSum == 0.0 ) continue; for( i = 0; i < bandCount; i++ ) padfPixelSum[i] /= dfWeightSum; if( MS_RENDERER_GD(psSrcImage->format) ) { if( !gdImageTrueColor(psSrcImage->img.gd) ) { int nResult = panCMap[(int) padfPixelSum[0]]; if( nResult != -1 ) { nSetPoints++; dstImg->pixels[nDstY][nDstX] = nResult; } } else { nSetPoints++; if( dfWeightSum > 0.99 ) dstImg->tpixels[nDstY][nDstX] = gdTrueColor( (int) padfPixelSum[0], (int) padfPixelSum[1], (int) padfPixelSum[2] ); else { int gd_color; int gd_alpha = (int) (127 - 127.9 * dfWeightSum); gd_alpha = MAX(0,MIN(127,gd_alpha)); gd_color = gdTrueColorAlpha( (int) padfPixelSum[0], (int) padfPixelSum[1], (int) padfPixelSum[2], gd_alpha ); dstImg->tpixels[nDstY][nDstX] = msAlphaBlend( dstImg->tpixels[nDstY][nDstX], gd_color ); } } } else if( MS_RENDERER_RAWDATA(psSrcImage->format) ) { int band; for( band = 0; band < psSrcImage->format->bands; band++ ) { if( psSrcImage->format->imagemode == MS_IMAGEMODE_INT16 ) { psDstImage->img.raw_16bit[ nDstX + nDstY * psDstImage->width + band*psDstImage->width*psDstImage->height] = (short) padfPixelSum[band]; } else if( psSrcImage->format->imagemode == MS_IMAGEMODE_FLOAT32) { psDstImage->img.raw_float[ nDstX + nDstY * psDstImage->width + band*psDstImage->width*psDstImage->height] = (float) padfPixelSum[band]; } else if( psSrcImage->format->imagemode == MS_IMAGEMODE_BYTE ) { psDstImage->img.raw_byte[ nDstX + nDstY * psDstImage->width + band*psDstImage->width*psDstImage->height] = (unsigned char) padfPixelSum[band]; } } } } } free( padfPixelSum ); free( panSuccess ); free( x ); free( y ); /* -------------------------------------------------------------------- */ /* Some debugging output. */ /* -------------------------------------------------------------------- */ if( nFailedPoints > 0 && debug ) { char szMsg[256]; sprintf( szMsg, "msBilinearRasterResampler: " "%d failed to transform, %d actually set.\n", nFailedPoints, nSetPoints ); msDebug( szMsg ); } return 0; } /************************************************************************/ /* msAverageSample() */ /************************************************************************/ static int msAverageSample( imageObj *psSrcImage, double dfXMin, double dfYMin, double dfXMax, double dfYMax, colorObj *offsite, double *padfPixelSum, double *pdfAlpha01 ) { int nXMin, nXMax, nYMin, nYMax, iX, iY; double dfWeightSum = 0.0; double dfMaxWeight = 0.0; nXMin = (int) dfXMin; nYMin = (int) dfYMin; nXMax = (int) ceil(dfXMax); nYMax = (int) ceil(dfYMax); *pdfAlpha01 = 0.0; for( iY = nYMin; iY < nYMax; iY++ ) { double dfYCellMin, dfYCellMax; dfYCellMin = MAX(iY,dfYMin); dfYCellMax = MIN(iY+1,dfYMax); for( iX = nXMin; iX < nXMax; iX++ ) { double dfXCellMin, dfXCellMax, dfWeight; dfXCellMin = MAX(iX,dfXMin); dfXCellMax = MIN(iX+1,dfXMax); dfWeight = (dfXCellMax-dfXCellMin) * (dfYCellMax-dfYCellMin); msSourceSample( psSrcImage, iX, iY, padfPixelSum, dfWeight, &dfWeightSum, offsite ); dfMaxWeight += dfWeight; } } if( dfWeightSum == 0.0 ) return MS_FALSE; for( iX = 0; iX < 4; iX++ ) padfPixelSum[iX] /= dfWeightSum; *pdfAlpha01 = dfWeightSum / dfMaxWeight; return MS_TRUE; } /************************************************************************/ /* msAverageRasterResample() */ /************************************************************************/ static int msAverageRasterResampler( imageObj *psSrcImage, colorObj offsite, imageObj *psDstImage, int *panCMap, SimpleTransformer pfnTransform, void *pCBData, int debug ) { double *x1, *y1, *x2, *y2; int nDstX, nDstY; int *panSuccess1, *panSuccess2; int nDstXSize = psDstImage->width; int nDstYSize = psDstImage->height; int nFailedPoints = 0, nSetPoints = 0; double *padfPixelSum; gdImagePtr srcImg, dstImg; int bandCount = MAX(4,psSrcImage->format->bands); padfPixelSum = (double *) malloc(sizeof(double) * bandCount); srcImg = psSrcImage->img.gd; dstImg = psDstImage->img.gd; x1 = (double *) malloc( sizeof(double) * (nDstXSize+1) ); y1 = (double *) malloc( sizeof(double) * (nDstXSize+1) ); x2 = (double *) malloc( sizeof(double) * (nDstXSize+1) ); y2 = (double *) malloc( sizeof(double) * (nDstXSize+1) ); panSuccess1 = (int *) malloc( sizeof(int) * (nDstXSize+1) ); panSuccess2 = (int *) malloc( sizeof(int) * (nDstXSize+1) ); for( nDstY = 0; nDstY < nDstYSize; nDstY++ ) { for( nDstX = 0; nDstX <= nDstXSize; nDstX++ ) { x1[nDstX] = nDstX; y1[nDstX] = nDstY; x2[nDstX] = nDstX; y2[nDstX] = nDstY+1; } pfnTransform( pCBData, nDstXSize+1, x1, y1, panSuccess1 ); pfnTransform( pCBData, nDstXSize+1, x2, y2, panSuccess2 ); for( nDstX = 0; nDstX < nDstXSize; nDstX++ ) { double dfXMin, dfYMin, dfXMax, dfYMax; double dfAlpha01; /* Do not generate a pixel unless all four corners transformed */ if( !panSuccess1[nDstX] || !panSuccess1[nDstX+1] || !panSuccess2[nDstX] || !panSuccess2[nDstX+1] ) { nFailedPoints++; continue; } dfXMin = MIN(MIN(x1[nDstX],x1[nDstX+1]), MIN(x2[nDstX],x2[nDstX+1])); dfYMin = MIN(MIN(y1[nDstX],y1[nDstX+1]), MIN(y2[nDstX],y2[nDstX+1])); dfXMax = MAX(MAX(x1[nDstX],x1[nDstX+1]), MAX(x2[nDstX],x2[nDstX+1])); dfYMax = MAX(MAX(y1[nDstX],y1[nDstX+1]), MAX(y2[nDstX],y2[nDstX+1])); dfXMin = MAX(dfXMin,0); dfYMin = MAX(dfYMin,0); dfXMax = MIN(dfXMax,psSrcImage->width); dfYMax = MIN(dfYMax,psSrcImage->height); memset( padfPixelSum, 0, sizeof(double)*bandCount ); if( !msAverageSample( psSrcImage, dfXMin, dfYMin, dfXMax, dfYMax, &offsite, padfPixelSum, &dfAlpha01 ) ) continue; if( MS_RENDERER_GD(psSrcImage->format) ) { if( !gdImageTrueColor(psSrcImage->img.gd) ) { int nResult = panCMap[(int) padfPixelSum[0]]; if( nResult != -1 ) { nSetPoints++; dstImg->pixels[nDstY][nDstX] = nResult; } } else { nSetPoints++; if( dfAlpha01 > 0.99 ) dstImg->tpixels[nDstY][nDstX] = gdTrueColor( (int) padfPixelSum[0], (int) padfPixelSum[1], (int) padfPixelSum[2] ); else { int gd_color; int gd_alpha = (int) (127 - 127.9 * dfAlpha01); gd_alpha = MAX(0,MIN(127,gd_alpha)); gd_color = gdTrueColorAlpha( (int) padfPixelSum[0], (int) padfPixelSum[1], (int) padfPixelSum[2], gd_alpha ); dstImg->tpixels[nDstY][nDstX] = msAlphaBlend( dstImg->tpixels[nDstY][nDstX], gd_color ); } } } else if( MS_RENDERER_RAWDATA(psSrcImage->format) ) { int band; for( band = 0; band < psSrcImage->format->bands; band++ ) { if( psSrcImage->format->imagemode == MS_IMAGEMODE_INT16 ) { psDstImage->img.raw_16bit[ nDstX + nDstY * psDstImage->width + band*psDstImage->width*psDstImage->height] = (short) padfPixelSum[band]; } else if( psSrcImage->format->imagemode == MS_IMAGEMODE_FLOAT32) { psDstImage->img.raw_float[ nDstX + nDstY * psDstImage->width + band*psDstImage->width*psDstImage->height] = (float) padfPixelSum[band]; } else if( psSrcImage->format->imagemode == MS_IMAGEMODE_BYTE ) { psDstImage->img.raw_byte[ nDstX + nDstY * psDstImage->width + band*psDstImage->width*psDstImage->height] = (unsigned char) padfPixelSum[band]; } } } } } free( padfPixelSum ); free( panSuccess1 ); free( x1 ); free( y1 ); free( panSuccess2 ); free( x2 ); free( y2 ); /* -------------------------------------------------------------------- */ /* Some debugging output. */ /* -------------------------------------------------------------------- */ if( nFailedPoints > 0 && debug ) { char szMsg[256]; sprintf( szMsg, "msAverageRasterResampler: " "%d failed to transform, %d actually set.\n", nFailedPoints, nSetPoints ); msDebug( szMsg ); } return 0; } /************************************************************************/ /* ==================================================================== */ /* PROJ.4 based transformer. */ /* ==================================================================== */ /************************************************************************/ typedef struct { projectionObj *psSrcProjObj; projPJ psSrcProj; int bSrcIsGeographic; double adfInvSrcGeoTransform[6]; projectionObj *psDstProjObj; projPJ psDstProj; int bDstIsGeographic; double adfDstGeoTransform[6]; int bUseProj; } msProjTransformInfo; /************************************************************************/ /* msInitProjTransformer() */ /************************************************************************/ void *msInitProjTransformer( projectionObj *psSrc, double *padfSrcGeoTransform, projectionObj *psDst, double *padfDstGeoTransform ) { msProjTransformInfo *psPTInfo; psPTInfo = (msProjTransformInfo *) calloc(1,sizeof(msProjTransformInfo)); /* -------------------------------------------------------------------- */ /* We won't even use PROJ.4 if either coordinate system is */ /* NULL. */ /* -------------------------------------------------------------------- */ psPTInfo->bUseProj = (psSrc->proj != NULL && psDst->proj != NULL && msProjectionsDiffer( psSrc, psDst ) ); /* -------------------------------------------------------------------- */ /* Record source image information. We invert the source */ /* transformation for more convenient inverse application in */ /* the transformer. */ /* -------------------------------------------------------------------- */ psPTInfo->psSrcProj = psSrc->proj; if( psPTInfo->bUseProj ) psPTInfo->bSrcIsGeographic = pj_is_latlong(psSrc->proj); else psPTInfo->bSrcIsGeographic = MS_FALSE; if( !InvGeoTransform(padfSrcGeoTransform, psPTInfo->adfInvSrcGeoTransform) ) return NULL; /* -------------------------------------------------------------------- */ /* Record destination image information. */ /* -------------------------------------------------------------------- */ psPTInfo->psDstProj = psDst->proj; if( psPTInfo->bUseProj ) psPTInfo->bDstIsGeographic = pj_is_latlong(psDst->proj); else psPTInfo->bDstIsGeographic = MS_FALSE; memcpy( psPTInfo->adfDstGeoTransform, padfDstGeoTransform, sizeof(double) * 6 ); return psPTInfo; } /************************************************************************/ /* msFreeProjTransformer() */ /************************************************************************/ void msFreeProjTransformer( void * pCBData ) { free( pCBData ); } /************************************************************************/ /* msProjTransformer */ /************************************************************************/ int msProjTransformer( void *pCBData, int nPoints, double *x, double *y, int *panSuccess ) { int i; msProjTransformInfo *psPTInfo = (msProjTransformInfo*) pCBData; double x_out; /* -------------------------------------------------------------------- */ /* Transform into destination georeferenced space. */ /* -------------------------------------------------------------------- */ for( i = 0; i < nPoints; i++ ) { x_out = psPTInfo->adfDstGeoTransform[0] + psPTInfo->adfDstGeoTransform[1] * x[i] + psPTInfo->adfDstGeoTransform[2] * y[i]; y[i] = psPTInfo->adfDstGeoTransform[3] + psPTInfo->adfDstGeoTransform[4] * x[i] + psPTInfo->adfDstGeoTransform[5] * y[i]; x[i] = x_out; panSuccess[i] = 1; } /* -------------------------------------------------------------------- */ /* Transform from degrees to radians if geographic. */ /* -------------------------------------------------------------------- */ if( psPTInfo->bDstIsGeographic ) { for( i = 0; i < nPoints; i++ ) { x[i] = x[i] * DEG_TO_RAD; y[i] = y[i] * DEG_TO_RAD; } } /* -------------------------------------------------------------------- */ /* Transform back to source projection space. */ /* -------------------------------------------------------------------- */ if( psPTInfo->bUseProj ) { double *z; int tr_result; z = (double *) calloc(sizeof(double),nPoints); msAcquireLock( TLOCK_PROJ ); tr_result = pj_transform( psPTInfo->psDstProj, psPTInfo->psSrcProj, nPoints, 1, x, y, z); msReleaseLock( TLOCK_PROJ ); if( tr_result != 0 ) { free( z ); for( i = 0; i < nPoints; i++ ) panSuccess[i] = 0; return MS_FALSE; } free( z ); for( i = 0; i < nPoints; i++ ) { if( x[i] == HUGE_VAL || y[i] == HUGE_VAL ) panSuccess[i] = 0; } } /* -------------------------------------------------------------------- */ /* Transform back to degrees if source is geographic. */ /* -------------------------------------------------------------------- */ if( psPTInfo->bSrcIsGeographic ) { for( i = 0; i < nPoints; i++ ) { if( panSuccess[i] ) { x[i] = x[i] * RAD_TO_DEG; y[i] = y[i] * RAD_TO_DEG; } } } /* -------------------------------------------------------------------- */ /* Transform to source raster space. */ /* -------------------------------------------------------------------- */ for( i = 0; i < nPoints; i++ ) { if( panSuccess[i] ) { x_out = psPTInfo->adfInvSrcGeoTransform[0] + psPTInfo->adfInvSrcGeoTransform[1] * x[i] + psPTInfo->adfInvSrcGeoTransform[2] * y[i]; y[i] = psPTInfo->adfInvSrcGeoTransform[3] + psPTInfo->adfInvSrcGeoTransform[4] * x[i] + psPTInfo->adfInvSrcGeoTransform[5] * y[i]; x[i] = x_out; } else { x[i] = -1; y[i] = -1; } } return 1; } /************************************************************************/ /* ==================================================================== */ /* Approximate transformer. */ /* ==================================================================== */ /************************************************************************/ typedef struct { SimpleTransformer pfnBaseTransformer; void *pBaseCBData; double dfMaxError; } msApproxTransformInfo; /************************************************************************/ /* msInitApproxTransformer() */ /************************************************************************/ static void *msInitApproxTransformer( SimpleTransformer pfnBaseTransformer, void *pBaseCBData, double dfMaxError ) { msApproxTransformInfo *psATInfo; psATInfo = (msApproxTransformInfo *) malloc(sizeof(msApproxTransformInfo)); psATInfo->pfnBaseTransformer = pfnBaseTransformer; psATInfo->pBaseCBData = pBaseCBData; psATInfo->dfMaxError = dfMaxError; return psATInfo; } /************************************************************************/ /* msFreeApproxTransformer() */ /************************************************************************/ static void msFreeApproxTransformer( void * pCBData ) { free( pCBData ); } /************************************************************************/ /* msApproxTransformer */ /************************************************************************/ static int msApproxTransformer( void *pCBData, int nPoints, double *x, double *y, int *panSuccess ) { msApproxTransformInfo *psATInfo = (msApproxTransformInfo *) pCBData; double x2[3], y2[3], dfDeltaX, dfDeltaY, dfError, dfDist; int nMiddle, anSuccess2[3], i, bSuccess; nMiddle = (nPoints-1)/2; /* -------------------------------------------------------------------- */ /* Bail if our preconditions are not met, or if error is not */ /* acceptable. */ /* -------------------------------------------------------------------- */ if( y[0] != y[nPoints-1] || y[0] != y[nMiddle] || x[0] == x[nPoints-1] || x[0] == x[nMiddle] || psATInfo->dfMaxError == 0.0 || nPoints <= 5 ) { return psATInfo->pfnBaseTransformer( psATInfo->pBaseCBData, nPoints, x, y, panSuccess ); } /* -------------------------------------------------------------------- */ /* Transform first, last and middle point. */ /* -------------------------------------------------------------------- */ x2[0] = x[0]; y2[0] = y[0]; x2[1] = x[nMiddle]; y2[1] = y[nMiddle]; x2[2] = x[nPoints-1]; y2[2] = y[nPoints-1]; bSuccess = psATInfo->pfnBaseTransformer( psATInfo->pBaseCBData, 3, x2, y2, anSuccess2 ); if( !bSuccess || !anSuccess2[0] || !anSuccess2[1] || !anSuccess2[2] ) return psATInfo->pfnBaseTransformer( psATInfo->pBaseCBData, nPoints, x, y, panSuccess ); /* -------------------------------------------------------------------- */ /* Is the error at the middle acceptable relative to an */ /* interpolation of the middle position? */ /* -------------------------------------------------------------------- */ dfDeltaX = (x2[2] - x2[0]) / (x[nPoints-1] - x[0]); dfDeltaY = (y2[2] - y2[0]) / (x[nPoints-1] - x[0]); dfError = fabs((x2[0] + dfDeltaX * (x[nMiddle] - x[0])) - x2[1]) + fabs((y2[0] + dfDeltaY * (x[nMiddle] - x[0])) - y2[1]); if( dfError > psATInfo->dfMaxError ) { bSuccess = msApproxTransformer( psATInfo, nMiddle, x, y, panSuccess ); if( !bSuccess ) { return psATInfo->pfnBaseTransformer( psATInfo->pBaseCBData, nPoints, x, y, panSuccess ); } bSuccess = msApproxTransformer( psATInfo, nPoints - nMiddle, x+nMiddle, y+nMiddle, panSuccess+nMiddle ); if( !bSuccess ) { return psATInfo->pfnBaseTransformer( psATInfo->pBaseCBData, nPoints, x, y, panSuccess ); } return 1; } /* -------------------------------------------------------------------- */ /* Error is OK, linearly interpolate all points along line. */ /* -------------------------------------------------------------------- */ for( i = nPoints-1; i >= 0; i-- ) { dfDist = (x[i] - x[0]); y[i] = y2[0] + dfDeltaY * dfDist; x[i] = x2[0] + dfDeltaX * dfDist; panSuccess[i] = 1; } return 1; } /************************************************************************/ /* msTransformMapToSource() */ /* */ /* Compute the extents of the current map view if transformed */ /* onto the source raster. */ /************************************************************************/ static int msTransformMapToSource( int nDstXSize, int nDstYSize, double * adfDstGeoTransform, projectionObj *psDstProj, int nSrcXSize, int nSrcYSize, double * adfInvSrcGeoTransform, projectionObj *psSrcProj, rectObj *psSrcExtent, int bUseGrid ) { int nFailures = 0; #define EDGE_STEPS 10 #define MAX_SIZE ((EDGE_STEPS+1)*(EDGE_STEPS+1)) int i, nSamples = 0, bOutInit = 0; double dfRatio; double x[MAX_SIZE], y[MAX_SIZE], z[MAX_SIZE]; /* -------------------------------------------------------------------- */ /* Collect edges in map image pixel/line coordinates */ /* -------------------------------------------------------------------- */ if( !bUseGrid ) { for( dfRatio = 0.0; dfRatio <= 1.001; dfRatio += (1.0/EDGE_STEPS) ) { assert( nSamples < MAX_SIZE ); x[nSamples ] = dfRatio * nDstXSize; y[nSamples++] = 0.0; x[nSamples ] = dfRatio * nDstXSize; y[nSamples++] = nDstYSize; x[nSamples ] = 0.0; y[nSamples++] = dfRatio * nDstYSize; x[nSamples ] = nDstXSize; y[nSamples++] = dfRatio * nDstYSize; } } /* -------------------------------------------------------------------- */ /* Collect a grid in the hopes of a more accurate region. */ /* -------------------------------------------------------------------- */ else { double dfRatio2; for( dfRatio = 0.0; dfRatio <= 1.001; dfRatio += (1.0/EDGE_STEPS) ) { for( dfRatio2=0.0; dfRatio2 <= 1.001; dfRatio2 += (1.0/EDGE_STEPS)) { assert( nSamples < MAX_SIZE ); x[nSamples ] = dfRatio2 * nDstXSize; y[nSamples++] = dfRatio * nDstYSize; } } } /* -------------------------------------------------------------------- */ /* transform to map georeferenced units */ /* -------------------------------------------------------------------- */ for( i = 0; i < nSamples; i++ ) { double x_out, y_out; x_out = adfDstGeoTransform[0] + x[i] * adfDstGeoTransform[1] + y[i] * adfDstGeoTransform[2]; y_out = adfDstGeoTransform[3] + x[i] * adfDstGeoTransform[4] + y[i] * adfDstGeoTransform[5]; x[i] = x_out; y[i] = y_out; z[i] = 0.0; } /* -------------------------------------------------------------------- */ /* Transform to layer georeferenced coordinates. */ /* -------------------------------------------------------------------- */ if( psDstProj->proj && psSrcProj->proj ) { int tr_result; if( pj_is_latlong(psDstProj->proj) ) { for( i = 0; i < nSamples; i++ ) { x[i] = x[i] * DEG_TO_RAD; y[i] = y[i] * DEG_TO_RAD; } } msAcquireLock( TLOCK_PROJ ); tr_result = pj_transform( psDstProj->proj, psSrcProj->proj, nSamples, 1, x, y, z ); msReleaseLock( TLOCK_PROJ ); if( tr_result != 0 ) return MS_FALSE; if( pj_is_latlong(psSrcProj->proj) ) { for( i = 0; i < nSamples; i++ ) { if( x[i] != HUGE_VAL && y[i] != HUGE_VAL ) { x[i] = x[i] * RAD_TO_DEG; y[i] = y[i] * RAD_TO_DEG; } } } } /* -------------------------------------------------------------------- */ /* If we just using the edges (not a grid) and we go some */ /* errors, then we need to restart using a grid pattern. */ /* -------------------------------------------------------------------- */ if( !bUseGrid ) { for( i = 0; i < nSamples; i++ ) { if( x[i] == HUGE_VAL || y[i] == HUGE_VAL ) { return msTransformMapToSource( nDstXSize, nDstYSize, adfDstGeoTransform, psDstProj, nSrcXSize, nSrcYSize, adfInvSrcGeoTransform,psSrcProj, psSrcExtent, 1 ); } } } /* -------------------------------------------------------------------- */ /* transform to layer raster coordinates, and collect bounds. */ /* -------------------------------------------------------------------- */ for( i = 0; i < nSamples; i++ ) { double x_out, y_out; if( x[i] == HUGE_VAL || y[i] == HUGE_VAL ) { nFailures++; continue; } x_out = adfInvSrcGeoTransform[0] + x[i]*adfInvSrcGeoTransform[1] + y[i]*adfInvSrcGeoTransform[2]; y_out = adfInvSrcGeoTransform[3] + x[i]*adfInvSrcGeoTransform[4] + y[i]*adfInvSrcGeoTransform[5]; if( !bOutInit ) { psSrcExtent->minx = psSrcExtent->maxx = x_out; psSrcExtent->miny = psSrcExtent->maxy = y_out; bOutInit = 1; } else { psSrcExtent->minx = MIN(psSrcExtent->minx, x_out); psSrcExtent->maxx = MAX(psSrcExtent->maxx, x_out); psSrcExtent->miny = MIN(psSrcExtent->miny, y_out); psSrcExtent->maxy = MAX(psSrcExtent->maxy, y_out); } } if( !bOutInit ) return MS_FALSE; /* -------------------------------------------------------------------- */ /* If we had some failures, we need to expand the region to */ /* represent our very coarse sampling grid. */ /* -------------------------------------------------------------------- */ if( nFailures > 0 ) { int nGrowAmountX = (int) (psSrcExtent->maxx - psSrcExtent->minx)/EDGE_STEPS + 1; int nGrowAmountY = (int) (psSrcExtent->maxy - psSrcExtent->miny)/EDGE_STEPS + 1; psSrcExtent->minx = MAX(psSrcExtent->minx - nGrowAmountX,0); psSrcExtent->miny = MAX(psSrcExtent->miny - nGrowAmountY,0); psSrcExtent->maxx = MIN(psSrcExtent->maxx + nGrowAmountX,nSrcXSize); psSrcExtent->maxy = MIN(psSrcExtent->maxy + nGrowAmountY,nSrcYSize); } return MS_TRUE; } #endif /* def USE_PROJ */ #ifdef USE_GDAL /************************************************************************/ /* msResampleGDALToMap() */ /************************************************************************/ int msResampleGDALToMap( mapObj *map, layerObj *layer, imageObj *image, GDALDatasetH hDS ) { /* -------------------------------------------------------------------- */ /* We require PROJ.4 4.4.2 or later. Earlier versions don't */ /* have PJD_GRIDSHIFT. */ /* -------------------------------------------------------------------- */ #if !defined(PJD_GRIDSHIFT) && !defined(PJ_VERSION) msSetError(MS_PROJERR, "Projection support is not available, so msResampleGDALToMap() fails.", "msProjectRect()"); return(MS_FAILURE); #else int nSrcXSize, nSrcYSize, nDstXSize, nDstYSize; int result, bSuccess; double adfSrcGeoTransform[6], adfDstGeoTransform[6]; double adfInvSrcGeoTransform[6], dfNominalCellSize; rectObj sSrcExtent, sOrigSrcExtent; mapObj sDummyMap; imageObj *srcImage; void *pTCBData; void *pACBData; int anCMap[256]; char **papszAlteredProcessing = NULL; int nLoadImgXSize, nLoadImgYSize; double dfOversampleRatio; const char *resampleMode = CSLFetchNameValue( layer->processing, "RESAMPLE" ); if( resampleMode == NULL ) resampleMode = "NEAREST"; /* -------------------------------------------------------------------- */ /* We will require source and destination to have a valid */ /* projection object. */ /* -------------------------------------------------------------------- */ if( map->projection.proj == NULL || layer->projection.proj == NULL ) { if( layer->debug ) msDebug( "msResampleGDALToMap(): " "Either map or layer projection is NULL, assuming compatible.\n" ); } /* -------------------------------------------------------------------- */ /* Initialize some information. */ /* -------------------------------------------------------------------- */ nDstXSize = image->width; nDstYSize = image->height; memcpy( adfDstGeoTransform, map->gt.geotransform, sizeof(double)*6 ); msGetGDALGeoTransform( hDS, map, layer, adfSrcGeoTransform ); nSrcXSize = GDALGetRasterXSize( hDS ); nSrcYSize = GDALGetRasterYSize( hDS ); InvGeoTransform( adfSrcGeoTransform, adfInvSrcGeoTransform ); /* -------------------------------------------------------------------- */ /* We need to find the extents in the source layer projection */ /* of the output requested region. We will accomplish this by */ /* collecting the extents of a region around the edge of the */ /* destination chunk. */ /* -------------------------------------------------------------------- */ if( CSLFetchBoolean( layer->processing, "LOAD_WHOLE_IMAGE", FALSE ) ) bSuccess = FALSE; else bSuccess = msTransformMapToSource( nDstXSize, nDstYSize, adfDstGeoTransform, &(map->projection), nSrcXSize, nSrcYSize,adfInvSrcGeoTransform, &(layer->projection), &sSrcExtent, FALSE ); /* -------------------------------------------------------------------- */ /* If the transformation failed, it is likely that we have such */ /* broad extents that the projection transformation failed at */ /* points around the extents. If so, we will assume we need */ /* the whole raster. This and later assumptions are likely to */ /* result in the raster being loaded at a higher resolution */ /* than really needed but should give decent results. */ /* -------------------------------------------------------------------- */ if( !bSuccess ) { if( layer->debug ) msDebug( "msTransformMapToSource(): " "pj_transform() failed. Out of bounds? Loading whole image.\n" ); sSrcExtent.minx = 0; sSrcExtent.maxx = nSrcXSize; sSrcExtent.miny = 0; sSrcExtent.maxy = nSrcYSize; } /* -------------------------------------------------------------------- */ /* Project desired extents out by 2 pixels, and then strip to */ /* available data. */ /* -------------------------------------------------------------------- */ memcpy( &sOrigSrcExtent, &sSrcExtent, sizeof(sSrcExtent) ); sSrcExtent.minx = floor(sSrcExtent.minx-1.0); sSrcExtent.maxx = ceil (sSrcExtent.maxx+1.0); sSrcExtent.miny = floor(sSrcExtent.miny-1.0); sSrcExtent.maxy = ceil (sSrcExtent.maxy+1.0); sSrcExtent.minx = MAX(0,sSrcExtent.minx); sSrcExtent.maxx = MIN(sSrcExtent.maxx, nSrcXSize ); sSrcExtent.miny = MAX(sSrcExtent.miny, 0 ); sSrcExtent.maxy = MIN(sSrcExtent.maxy, nSrcYSize ); if( sSrcExtent.maxx <= sSrcExtent.minx || sSrcExtent.maxy <= sSrcExtent.miny ) { if( layer->debug ) msDebug( "msResampleGDALToMap(): no overlap ... no result.\n" ); return 0; } /* -------------------------------------------------------------------- */ /* Determine desired oversampling ratio. Default to 2.0 if not */ /* otherwise set. */ /* -------------------------------------------------------------------- */ dfOversampleRatio = 2.0; if( CSLFetchNameValue( layer->processing, "OVERSAMPLE_RATIO" ) != NULL ) { dfOversampleRatio = atof(CSLFetchNameValue( layer->processing, "OVERSAMPLE_RATIO" )); } /* -------------------------------------------------------------------- */ /* Decide on a resolution to read from the source image at. We */ /* will operate from full resolution data, if we are requesting */ /* at near to full resolution. Otherwise we will read the data */ /* at twice the resolution of the eventual map. */ /* -------------------------------------------------------------------- */ dfNominalCellSize = sqrt(adfSrcGeoTransform[1] * adfSrcGeoTransform[1] + adfSrcGeoTransform[2] * adfSrcGeoTransform[2]); if( (sOrigSrcExtent.maxx - sOrigSrcExtent.minx) > dfOversampleRatio * nDstXSize && !CSLFetchBoolean( layer->processing, "LOAD_FULL_RES_IMAGE", FALSE )) sDummyMap.cellsize = (dfNominalCellSize * (sOrigSrcExtent.maxx - sOrigSrcExtent.minx)) / (dfOversampleRatio * nDstXSize); else sDummyMap.cellsize = dfNominalCellSize; nLoadImgXSize = (int) MAX(1,(sSrcExtent.maxx - sSrcExtent.minx) * (dfNominalCellSize / sDummyMap.cellsize)); nLoadImgYSize = (int) MAX(1,(sSrcExtent.maxy - sSrcExtent.miny) * (dfNominalCellSize / sDummyMap.cellsize)); /* ** Because the previous calculation involved some round off, we need ** to fixup the cellsize to ensure the map region represents the whole ** RAW_WINDOW (at least in X). Re: bug 1715. */ sDummyMap.cellsize = ((sSrcExtent.maxx - sSrcExtent.minx) * dfNominalCellSize) / nLoadImgXSize; if( layer->debug ) msDebug( "msResampleGDALToMap in effect: cellsize = %f\n", sDummyMap.cellsize ); adfSrcGeoTransform[0] += + adfSrcGeoTransform[1] * sSrcExtent.minx + adfSrcGeoTransform[2] * sSrcExtent.miny; adfSrcGeoTransform[1] *= (sDummyMap.cellsize / dfNominalCellSize); adfSrcGeoTransform[2] *= (sDummyMap.cellsize / dfNominalCellSize); adfSrcGeoTransform[3] += + adfSrcGeoTransform[4] * sSrcExtent.minx + adfSrcGeoTransform[5] * sSrcExtent.miny; adfSrcGeoTransform[4] *= (sDummyMap.cellsize / dfNominalCellSize); adfSrcGeoTransform[5] *= (sDummyMap.cellsize / dfNominalCellSize); papszAlteredProcessing = CSLDuplicate( layer->processing ); papszAlteredProcessing = CSLSetNameValue( papszAlteredProcessing, "RAW_WINDOW", CPLSPrintf( "%d %d %d %d", (int) sSrcExtent.minx, (int) sSrcExtent.miny, (int) (sSrcExtent.maxx-sSrcExtent.minx), (int) (sSrcExtent.maxy-sSrcExtent.miny))); /* -------------------------------------------------------------------- */ /* We clone this without referencing it knowing that the */ /* srcImage will take a reference on it. The sDummyMap is */ /* destroyed off the stack, so the missing map reference is */ /* never a problem. The image's dereference of the */ /* outputformat during the msFreeImage() calls will result in */ /* the output format being cleaned up. */ /* */ /* We make a copy so we can easily modify the outputformat used */ /* for the temporary image to include transparentency support. */ /* -------------------------------------------------------------------- */ sDummyMap.outputformat = msCloneOutputFormat( image->format ); #ifdef USE_AGG if( MS_RENDERER_AGG(sDummyMap.outputformat)) { /* Fix for bug 2540 : The resampling methods used later assume that the pixel format of the dummy image is in gd format (alpha between 127 and 0). If we keep an AGG outpuformat, the image is initialised with its alpha=0, which is interpreted as fully opaque instead of fully transparent. The fix is to completely forget about AGG here, and create a GD imageObj */ sDummyMap.outputformat->renderer=MS_RENDER_WITH_GD; } #endif sDummyMap.width = nLoadImgXSize; sDummyMap.height = nLoadImgYSize; /* -------------------------------------------------------------------- */ /* If we are working in 256 color GD mode, allocate 0 as the */ /* transparent color on the temporary image so it will be */ /* initialized to see-through. We pick an arbitrary rgb tuple */ /* as our transparent color, but ensure it is initalized in the */ /* map so that normal transparent avoidance will apply. */ /* -------------------------------------------------------------------- */ if( MS_RENDERER_GD(sDummyMap.outputformat) && !gdImageTrueColor( image->img.gd ) ) { sDummyMap.outputformat->transparent = MS_TRUE; sDummyMap.imagecolor.red = 117; sDummyMap.imagecolor.green = 17; sDummyMap.imagecolor.blue = 191; } /* -------------------------------------------------------------------- */ /* If we are working in RGB mode ensure we produce an RGBA */ /* image so the transparency can be preserved. */ /* -------------------------------------------------------------------- */ else if( MS_RENDERER_GD(sDummyMap.outputformat) && gdImageTrueColor( image->img.gd ) ) { assert( sDummyMap.outputformat->imagemode == MS_IMAGEMODE_RGB || sDummyMap.outputformat->imagemode == MS_IMAGEMODE_RGBA ); sDummyMap.outputformat->transparent = MS_TRUE; sDummyMap.outputformat->imagemode = MS_IMAGEMODE_RGBA; sDummyMap.imagecolor.red = map->imagecolor.red; sDummyMap.imagecolor.green = map->imagecolor.green; sDummyMap.imagecolor.blue = map->imagecolor.blue; } /* -------------------------------------------------------------------- */ /* Setup a dummy map object we can use to read from the source */ /* raster, with the newly established extents, and resolution. */ /* -------------------------------------------------------------------- */ srcImage = msImageCreate( nLoadImgXSize, nLoadImgYSize, sDummyMap.outputformat, NULL, NULL, &sDummyMap); if (srcImage == NULL) return -1; /* msSetError() should have been called already */ /* -------------------------------------------------------------------- */ /* Draw into the temporary image. Temporarily replace the */ /* layer processing directive so that we use our RAW_WINDOW. */ /* -------------------------------------------------------------------- */ { char **papszSavedProcessing = layer->processing; layer->processing = papszAlteredProcessing; result = msDrawRasterLayerGDAL( &sDummyMap, layer, srcImage, hDS ); layer->processing = papszSavedProcessing; CSLDestroy( papszAlteredProcessing ); if( result ) { msFreeImage( srcImage ); return result; } } /* -------------------------------------------------------------------- */ /* Do we need to generate a colormap remapping, potentially */ /* allocating new colors on the destination color map? */ /* -------------------------------------------------------------------- */ if( MS_RENDERER_GD(srcImage->format) && !gdImageTrueColor( srcImage->img.gd ) ) { int iColor, nColorCount; anCMap[0] = -1; /* color zero is always transparent */ nColorCount = gdImageColorsTotal( srcImage->img.gd ); for( iColor = 1; iColor < nColorCount; iColor++ ) { anCMap[iColor] = msAddColorGD( map, image->img.gd, 0, gdImageRed( srcImage->img.gd, iColor ), gdImageGreen( srcImage->img.gd, iColor ), gdImageBlue( srcImage->img.gd, iColor ) ); } for( iColor = nColorCount; iColor < 256; iColor++ ) anCMap[iColor] = -1; } /* -------------------------------------------------------------------- */ /* Setup transformations between our source image, and the */ /* target map image. */ /* -------------------------------------------------------------------- */ pTCBData = msInitProjTransformer( &(layer->projection), adfSrcGeoTransform, &(map->projection), adfDstGeoTransform ); if( pTCBData == NULL ) { if( layer->debug ) msDebug( "msInitProjTransformer() returned NULL.\n" ); msFreeImage( srcImage ); return MS_PROJERR; } /* -------------------------------------------------------------------- */ /* It is cheaper to use linear approximations as long as our */ /* error is modest (less than 0.333 pixels). */ /* -------------------------------------------------------------------- */ pACBData = msInitApproxTransformer( msProjTransformer, pTCBData, 0.333 ); /* -------------------------------------------------------------------- */ /* Perform the resampling. */ /* -------------------------------------------------------------------- */ if( EQUAL(resampleMode,"AVERAGE") ) result = msAverageRasterResampler( srcImage, layer->offsite, image, anCMap, msApproxTransformer, pACBData, layer->debug ); else if( EQUAL(resampleMode,"BILINEAR") ) result = msBilinearRasterResampler( srcImage, layer->offsite, image, anCMap, msApproxTransformer, pACBData, layer->debug ); else result = msNearestRasterResampler( srcImage, layer->offsite, image, anCMap, msApproxTransformer, pACBData, layer->debug ); /* -------------------------------------------------------------------- */ /* cleanup */ /* -------------------------------------------------------------------- */ msFreeImage( srcImage ); msFreeProjTransformer( pTCBData ); msFreeApproxTransformer( pACBData ); return result; #endif } #endif /* def USE_GDAL */