EVOLUTION-MANAGER

Edit File: gdaldem_lib.cpp

/****************************************************************************** * * Project: GDAL DEM Utilities * Purpose: * Authors: Matthew Perry, perrygeo at gmail.com * Even Rouault, even dot rouault at mines dash paris dot org * Howard Butler, hobu.inc at gmail.com * Chris Yesson, chris dot yesson at ioz dot ac dot uk * ****************************************************************************** * Copyright (c) 2006, 2009 Matthew Perry * Copyright (c) 2009-2013, Even Rouault <even dot rouault at mines-paris dot org> * Portions derived from GRASS 4.1 (public domain) See * http://trac.osgeo.org/gdal/ticket/2975 for more information regarding * history of this code * * 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 or substantial portions of the 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. **************************************************************************** * * Slope and aspect calculations based on original method for GRASS GIS 4.1 * by Michael Shapiro, U.S.Army Construction Engineering Research Laboratory * Olga Waupotitsch, U.S.Army Construction Engineering Research Laboratory * Marjorie Larson, U.S.Army Construction Engineering Research Laboratory * as found in GRASS's r.slope.aspect module. * * Horn's formula is used to find the first order derivatives in x and y directions * for slope and aspect calculations: Horn, B. K. P. (1981). * "Hill Shading and the Reflectance Map", Proceedings of the IEEE, 69(1):14-47. * * Other reference : * Burrough, P.A. and McDonell, R.A., 1998. Principles of Geographical Information * Systems. p. 190. * * Shaded relief based on original method for GRASS GIS 4.1 by Jim Westervelt, * U.S. Army Construction Engineering Research Laboratory * as found in GRASS's r.shaded.relief (formerly shade.rel.sh) module. * ref: "r.mapcalc: An Algebra for GIS and Image Processing", * by Michael Shapiro and Jim Westervelt, U.S. Army Construction Engineering * Research Laboratory (March/1991) * * Color table of named colors and lookup code derived from src/libes/gis/named_colr.c * of GRASS 4.1 * * TRI - Terrain Ruggedness Index is as described in Wilson et al. (2007) * this is based on the method of Valentine et al. (2004) * * TPI - Topographic Position Index follows the description in * Wilson et al. (2007), following Weiss (2001). The radius is fixed * at 1 cell width/height * * Roughness - follows the definition in Wilson et al. (2007), which follows * Dartnell (2000). * * References for TRI/TPI/Roughness: * Dartnell, P. 2000. Applying Remote Sensing Techniques to map Seafloor * Geology/Habitat Relationships. Masters Thesis, San Francisco State * University, pp. 108. * Valentine, P. C., S. J. Fuller, L. A. Scully. 2004. Terrain Ruggedness * Analysis and Distribution of Boulder Ridges in the Stellwagen Bank National * Marine Sanctuary Region (poster). Galway, Ireland: 5th International * Symposium on Marine Geological and Biological Habitat Mapping (GeoHAB), * May 2004. * Weiss, A. D. 2001. Topographic Positions and Landforms Analysis (poster), * ESRI International User Conference, July 2001. San Diego, CA: ESRI. * Wilson, M. F. J.; O'Connell, B.; Brown, C.; Guinan, J. C. & Grehan, A. J. * Multiscale terrain analysis of multibeam bathymetry data for habitat mapping * on the continental slope Marine Geodesy, 2007, 30, 3-35 ****************************************************************************/ // Include before others for mingw for VSIStatBufL #include "cpl_conv.h" #include "cpl_port.h" #include "gdal_utils.h" #include "gdal_utils_priv.h" #include "commonutils.h" #include <cfloat> #include <cmath> #include <cstdio> #include <cstdlib> #include <cstring> #if HAVE_SYS_STAT_H #include <sys/stat.h> #endif #include <algorithm> #include <limits> #include "cpl_error.h" #include "cpl_progress.h" #include "cpl_string.h" #include "cpl_vsi.h" #include "gdal.h" #include "gdal_priv.h" #if defined(__SSE2__) || defined(_M_X64) #define HAVE_16_SSE_REG #define HAVE_SSE2 #include "emmintrin.h" #endif CPL_CVSID("$Id: gdaldem_lib.cpp 4db55e60ad36bf21a576d1179c0f1788d2f8dbf2 2019-03-04 22:16:21 +0100 Martin Ždila $") static const double kdfDegreesToRadians = M_PI / 180.0; static const double kdfRadiansToDegrees = 180.0 / M_PI; typedef enum { COLOR_SELECTION_INTERPOLATE, COLOR_SELECTION_NEAREST_ENTRY, COLOR_SELECTION_EXACT_ENTRY } ColorSelectionMode; struct GDALDEMProcessingOptions { /*! output format. Use the short format name. */ char *pszFormat; /*! the progress function to use */ GDALProgressFunc pfnProgress; /*! pointer to the progress data variable */ void *pProgressData; double z; double scale; double az; double alt; int slopeFormat; bool bAddAlpha; bool bZeroForFlat; bool bAngleAsAzimuth; ColorSelectionMode eColorSelectionMode; bool bComputeAtEdges; bool bZevenbergenThorne; bool bCombined; bool bIgor; bool bMultiDirectional; char** papszCreateOptions; int nBand; }; /************************************************************************/ /* ComputeVal() */ /************************************************************************/ template<class T> struct GDALGeneric3x3ProcessingAlg { typedef float (*type) (const T* pafWindow, float fDstNoDataValue, void* pData); }; template<class T> struct GDALGeneric3x3ProcessingAlg_multisample { typedef int (*type) (const T* pafThreeLineWin, int nLine1Off, int nLine2Off, int nLine3Off, int nXSize, void* pData, float* pafOutputBuf); }; template<class T> static float ComputeVal( bool bSrcHasNoData, T fSrcNoDataValue, bool bIsSrcNoDataNan, T* afWin, float fDstNoDataValue, typename GDALGeneric3x3ProcessingAlg<T>::type pfnAlg, void* pData, bool bComputeAtEdges ); template<> float ComputeVal( bool bSrcHasNoData, float fSrcNoDataValue, bool bIsSrcNoDataNan, float* afWin, float fDstNoDataValue, GDALGeneric3x3ProcessingAlg<float>::type pfnAlg, void* pData, bool bComputeAtEdges ) { if( bSrcHasNoData && ((!bIsSrcNoDataNan && ARE_REAL_EQUAL(afWin[4], fSrcNoDataValue)) || (bIsSrcNoDataNan && CPLIsNan(afWin[4]))) ) { return fDstNoDataValue; } else if( bSrcHasNoData ) { for( int k = 0; k < 9; k++ ) { if( (!bIsSrcNoDataNan && ARE_REAL_EQUAL(afWin[k], fSrcNoDataValue)) || (bIsSrcNoDataNan && CPLIsNan(afWin[k])) ) { if( bComputeAtEdges ) afWin[k] = afWin[4]; else return fDstNoDataValue; } } } return pfnAlg(afWin, fDstNoDataValue, pData); } template<> float ComputeVal( bool bSrcHasNoData, GInt32 fSrcNoDataValue, bool /* bIsSrcNoDataNan */, GInt32* afWin, float fDstNoDataValue, GDALGeneric3x3ProcessingAlg<GInt32>::type pfnAlg, void* pData, bool bComputeAtEdges ) { if( bSrcHasNoData && afWin[4] == fSrcNoDataValue ) { return fDstNoDataValue; } else if( bSrcHasNoData ) { for( int k = 0; k < 9; k++ ) { if( afWin[k] == fSrcNoDataValue ) { if( bComputeAtEdges ) afWin[k] = afWin[4]; else return fDstNoDataValue; } } } return pfnAlg(afWin, fDstNoDataValue, pData); } /************************************************************************/ /* INTERPOL() */ /************************************************************************/ template<class T> static T INTERPOL(T a, T b, int bSrcHasNodata, T fSrcNoDataValue); template<> float INTERPOL(float a, float b, int bSrcHasNoData, float fSrcNoDataValue) { return ((bSrcHasNoData && (ARE_REAL_EQUAL(a, fSrcNoDataValue) || ARE_REAL_EQUAL(b, fSrcNoDataValue))) ? fSrcNoDataValue : 2 * (a) - (b)); } template<> GInt32 INTERPOL(GInt32 a, GInt32 b, int bSrcHasNoData, GInt32 fSrcNoDataValue) { if( bSrcHasNoData && ((a == fSrcNoDataValue) || (b == fSrcNoDataValue)) ) return fSrcNoDataValue; int nVal = 2 * a - b; if( bSrcHasNoData && fSrcNoDataValue == nVal ) return nVal + 1; return nVal; } /************************************************************************/ /* GDALGeneric3x3Processing() */ /************************************************************************/ template<class T> static CPLErr GDALGeneric3x3Processing( GDALRasterBandH hSrcBand, GDALRasterBandH hDstBand, typename GDALGeneric3x3ProcessingAlg<T>::type pfnAlg, typename GDALGeneric3x3ProcessingAlg_multisample<T>::type pfnAlg_multisample, void *pData, bool bComputeAtEdges, GDALProgressFunc pfnProgress, void *pProgressData ) { if( pfnProgress == nullptr ) pfnProgress = GDALDummyProgress; /* -------------------------------------------------------------------- */ /* Initialize progress counter. */ /* -------------------------------------------------------------------- */ if( !pfnProgress( 0.0, nullptr, pProgressData ) ) { CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" ); return CE_Failure; } const int nXSize = GDALGetRasterBandXSize(hSrcBand); const int nYSize = GDALGetRasterBandYSize(hSrcBand); // 1 line destination buffer. float *pafOutputBuf = static_cast<float *>( VSI_MALLOC2_VERBOSE(sizeof(float), nXSize)); // 3 line rotating source buffer. T *pafThreeLineWin = static_cast<T *>( VSI_MALLOC2_VERBOSE(3 * sizeof(T), nXSize + 1)); if( pafOutputBuf == nullptr || pafThreeLineWin == nullptr ) { VSIFree(pafOutputBuf); VSIFree(pafThreeLineWin); return CE_Failure; } GDALDataType eReadDT; int bSrcHasNoData = FALSE; const double dfNoDataValue = GDALGetRasterNoDataValue(hSrcBand, &bSrcHasNoData); int bIsSrcNoDataNan = FALSE; T fSrcNoDataValue = 0; if( std::numeric_limits<T>::is_integer ) { eReadDT = GDT_Int32; if( bSrcHasNoData ) { GDALDataType eSrcDT = GDALGetRasterDataType( hSrcBand ); CPLAssert( eSrcDT == GDT_Byte || eSrcDT == GDT_UInt16 || eSrcDT == GDT_Int16 ); const int nMinVal = (eSrcDT == GDT_Byte ) ? 0 : (eSrcDT == GDT_UInt16) ? 0 : -32768; const int nMaxVal = (eSrcDT == GDT_Byte ) ? 255 : (eSrcDT == GDT_UInt16) ? 65535 : 32767; if( fabs(dfNoDataValue - floor(dfNoDataValue + 0.5)) < 1e-2 && dfNoDataValue >= nMinVal && dfNoDataValue <= nMaxVal ) { fSrcNoDataValue = static_cast<T>(floor(dfNoDataValue + 0.5)); } else { bSrcHasNoData = FALSE; } } } else { eReadDT = GDT_Float32; fSrcNoDataValue = static_cast<T>(dfNoDataValue); bIsSrcNoDataNan = bSrcHasNoData && CPLIsNan(dfNoDataValue); } int bDstHasNoData = FALSE; float fDstNoDataValue = static_cast<float>(GDALGetRasterNoDataValue(hDstBand, &bDstHasNoData)); if( !bDstHasNoData ) fDstNoDataValue = 0.0; int nLine1Off = 0; int nLine2Off = nXSize; int nLine3Off = 2*nXSize; // Move a 3x3 pafWindow over each cell // (where the cell in question is #4) // // 0 1 2 // 3 4 5 // 6 7 8 /* Preload the first 2 lines */ bool abLineHasNoDataValue[3] = { CPL_TO_BOOL(bSrcHasNoData), CPL_TO_BOOL(bSrcHasNoData), CPL_TO_BOOL(bSrcHasNoData) }; // Create an extra scope for VC12 to ignore i. { for( int i = 0; i < 2 && i < nYSize; i++ ) { if( GDALRasterIO( hSrcBand, GF_Read, 0, i, nXSize, 1, pafThreeLineWin + i * nXSize, nXSize, 1, eReadDT, 0, 0) != CE_None ) { CPLFree(pafOutputBuf); CPLFree(pafThreeLineWin); return CE_Failure; } if( std::numeric_limits<T>::is_integer && bSrcHasNoData ) { abLineHasNoDataValue[i] = false; for( int iX = 0; iX < nXSize; iX++ ) { if( pafThreeLineWin[i * nXSize + iX] == fSrcNoDataValue ) { abLineHasNoDataValue[i] = true; break; } } } } } // End extra scope for VC12 CPLErr eErr = CE_None; if( bComputeAtEdges && nXSize >= 2 && nYSize >= 2 ) { for( int j = 0; j < nXSize; j++ ) { int jmin = (j == 0) ? j : j - 1; int jmax = (j == nXSize - 1) ? j : j + 1; T afWin[9] = { INTERPOL(pafThreeLineWin[jmin], pafThreeLineWin[nXSize + jmin], bSrcHasNoData, fSrcNoDataValue), INTERPOL(pafThreeLineWin[j], pafThreeLineWin[nXSize + j], bSrcHasNoData, fSrcNoDataValue), INTERPOL(pafThreeLineWin[jmax], pafThreeLineWin[nXSize + jmax], bSrcHasNoData, fSrcNoDataValue), pafThreeLineWin[jmin], pafThreeLineWin[j], pafThreeLineWin[jmax], pafThreeLineWin[nXSize + jmin], pafThreeLineWin[nXSize + j], pafThreeLineWin[nXSize + jmax] }; pafOutputBuf[j] = ComputeVal( CPL_TO_BOOL(bSrcHasNoData), fSrcNoDataValue, CPL_TO_BOOL(bIsSrcNoDataNan), afWin, fDstNoDataValue, pfnAlg, pData, bComputeAtEdges); } eErr = GDALRasterIO(hDstBand, GF_Write, 0, 0, nXSize, 1, pafOutputBuf, nXSize, 1, GDT_Float32, 0, 0); } else { // Exclude the edges for( int j = 0; j < nXSize; j++ ) { pafOutputBuf[j] = fDstNoDataValue; } eErr = GDALRasterIO(hDstBand, GF_Write, 0, 0, nXSize, 1, pafOutputBuf, nXSize, 1, GDT_Float32, 0, 0); if( eErr == CE_None && nYSize > 1 ) { eErr = GDALRasterIO(hDstBand, GF_Write, 0, nYSize - 1, nXSize, 1, pafOutputBuf, nXSize, 1, GDT_Float32, 0, 0); } } if( eErr != CE_None ) { CPLFree(pafOutputBuf); CPLFree(pafThreeLineWin); return eErr; } int i = 1; // Used after for. for( ; i < nYSize-1; i++ ) { /* Read third line of the line buffer */ eErr = GDALRasterIO( hSrcBand, GF_Read, 0, i+1, nXSize, 1, pafThreeLineWin + nLine3Off, nXSize, 1, eReadDT, 0, 0); if( eErr != CE_None ) { CPLFree(pafOutputBuf); CPLFree(pafThreeLineWin); return eErr; } // In case none of the 3 lines have nodata values, then no need to // check it in ComputeVal() bool bOneOfThreeLinesHasNoData = CPL_TO_BOOL(bSrcHasNoData); if( std::numeric_limits<T>::is_integer && bSrcHasNoData ) { bool bLastLineHasNoDataValue = false; int iX = 0; for( ; iX + 3 < nXSize; iX +=4 ) { if( pafThreeLineWin[nLine3Off + iX] == fSrcNoDataValue || pafThreeLineWin[nLine3Off + iX + 1] == fSrcNoDataValue || pafThreeLineWin[nLine3Off + iX + 2] == fSrcNoDataValue || pafThreeLineWin[nLine3Off + iX + 3] == fSrcNoDataValue ) { bLastLineHasNoDataValue = true; break; } } if( !bLastLineHasNoDataValue ) { for( ; iX < nXSize; iX++ ) { if( pafThreeLineWin[nLine3Off + iX] == fSrcNoDataValue ) { bLastLineHasNoDataValue = true; } } } abLineHasNoDataValue[nLine3Off / nXSize] = bLastLineHasNoDataValue; bOneOfThreeLinesHasNoData = abLineHasNoDataValue[0] || abLineHasNoDataValue[1] || abLineHasNoDataValue[2]; } if( bComputeAtEdges && nXSize >= 2 ) { int j = 0; T afWin[9] = { INTERPOL(pafThreeLineWin[nLine1Off + j], pafThreeLineWin[nLine1Off + j+1], bSrcHasNoData, fSrcNoDataValue), pafThreeLineWin[nLine1Off + j], pafThreeLineWin[nLine1Off + j+1], INTERPOL(pafThreeLineWin[nLine2Off + j], pafThreeLineWin[nLine2Off + j+1], bSrcHasNoData, fSrcNoDataValue), pafThreeLineWin[nLine2Off + j], pafThreeLineWin[nLine2Off + j+1], INTERPOL(pafThreeLineWin[nLine3Off + j], pafThreeLineWin[nLine3Off + j+1], bSrcHasNoData, fSrcNoDataValue), pafThreeLineWin[nLine3Off + j], pafThreeLineWin[nLine3Off + j+1] }; pafOutputBuf[j] = ComputeVal( CPL_TO_BOOL(bOneOfThreeLinesHasNoData), fSrcNoDataValue, CPL_TO_BOOL(bIsSrcNoDataNan), afWin, fDstNoDataValue, pfnAlg, pData, bComputeAtEdges); } else { // Exclude the edges pafOutputBuf[0] = fDstNoDataValue; } int j = 1; if( pfnAlg_multisample && !bOneOfThreeLinesHasNoData ) { j = pfnAlg_multisample(pafThreeLineWin, nLine1Off, nLine2Off, nLine3Off, nXSize, pData, pafOutputBuf); } for( ; j < nXSize - 1; j++ ) { T afWin[9] = { pafThreeLineWin[nLine1Off + j-1], pafThreeLineWin[nLine1Off + j], pafThreeLineWin[nLine1Off + j+1], pafThreeLineWin[nLine2Off + j-1], pafThreeLineWin[nLine2Off + j], pafThreeLineWin[nLine2Off + j+1], pafThreeLineWin[nLine3Off + j-1], pafThreeLineWin[nLine3Off + j], pafThreeLineWin[nLine3Off + j+1] }; pafOutputBuf[j] = ComputeVal( CPL_TO_BOOL(bOneOfThreeLinesHasNoData), fSrcNoDataValue, CPL_TO_BOOL(bIsSrcNoDataNan), afWin, fDstNoDataValue, pfnAlg, pData, bComputeAtEdges); } if( bComputeAtEdges && nXSize >= 2 ) { j = nXSize - 1; T afWin[9] = { pafThreeLineWin[nLine1Off + j-1], pafThreeLineWin[nLine1Off + j], INTERPOL(pafThreeLineWin[nLine1Off + j], pafThreeLineWin[nLine1Off + j-1], bSrcHasNoData, fSrcNoDataValue), pafThreeLineWin[nLine2Off + j-1], pafThreeLineWin[nLine2Off + j], INTERPOL(pafThreeLineWin[nLine2Off + j], pafThreeLineWin[nLine2Off + j-1], bSrcHasNoData, fSrcNoDataValue), pafThreeLineWin[nLine3Off + j-1], pafThreeLineWin[nLine3Off + j], INTERPOL(pafThreeLineWin[nLine3Off + j], pafThreeLineWin[nLine3Off + j-1], bSrcHasNoData, fSrcNoDataValue) }; pafOutputBuf[j] = ComputeVal( CPL_TO_BOOL(bOneOfThreeLinesHasNoData), fSrcNoDataValue, CPL_TO_BOOL(bIsSrcNoDataNan), afWin, fDstNoDataValue, pfnAlg, pData, bComputeAtEdges); } else { // Exclude the edges if( nXSize > 1 ) pafOutputBuf[nXSize - 1] = fDstNoDataValue; } /* ----------------------------------------- * Write Line to Raster */ eErr = GDALRasterIO(hDstBand, GF_Write, 0, i, nXSize, 1, pafOutputBuf, nXSize, 1, GDT_Float32, 0, 0); if( eErr != CE_None ) { CPLFree(pafOutputBuf); CPLFree(pafThreeLineWin); return eErr; } if( !pfnProgress( 1.0 * (i+1) / nYSize, nullptr, pProgressData ) ) { CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" ); eErr = CE_Failure; CPLFree(pafOutputBuf); CPLFree(pafThreeLineWin); return eErr; } const int nTemp = nLine1Off; nLine1Off = nLine2Off; nLine2Off = nLine3Off; nLine3Off = nTemp; } if( bComputeAtEdges && nXSize >= 2 && nYSize >= 2 ) { for( int j = 0; j < nXSize; j++ ) { int jmin = (j == 0) ? j : j - 1; int jmax = (j == nXSize - 1) ? j : j + 1; T afWin[9] = { pafThreeLineWin[nLine1Off + jmin], pafThreeLineWin[nLine1Off + j], pafThreeLineWin[nLine1Off + jmax], pafThreeLineWin[nLine2Off + jmin], pafThreeLineWin[nLine2Off + j], pafThreeLineWin[nLine2Off + jmax], INTERPOL(pafThreeLineWin[nLine2Off + jmin], pafThreeLineWin[nLine1Off + jmin], bSrcHasNoData, fSrcNoDataValue), INTERPOL(pafThreeLineWin[nLine2Off + j], pafThreeLineWin[nLine1Off + j], bSrcHasNoData, fSrcNoDataValue), INTERPOL(pafThreeLineWin[nLine2Off + jmax], pafThreeLineWin[nLine1Off + jmax], bSrcHasNoData, fSrcNoDataValue), }; pafOutputBuf[j] = ComputeVal( CPL_TO_BOOL(bSrcHasNoData), fSrcNoDataValue, CPL_TO_BOOL(bIsSrcNoDataNan), afWin, fDstNoDataValue, pfnAlg, pData, bComputeAtEdges); } eErr = GDALRasterIO(hDstBand, GF_Write, 0, i, nXSize, 1, pafOutputBuf, nXSize, 1, GDT_Float32, 0, 0); if( eErr != CE_None ) { CPLFree(pafOutputBuf); CPLFree(pafThreeLineWin); return eErr; } } pfnProgress( 1.0, nullptr, pProgressData ); eErr = CE_None; CPLFree(pafOutputBuf); CPLFree(pafThreeLineWin); return eErr; } /************************************************************************/ /* GradientAlg */ /************************************************************************/ typedef enum { HORN, ZEVENBERGEN_THORNE } GradientAlg; template<class T, GradientAlg alg> struct Gradient { static void inline calc(const T* afWin, double inv_ewres, double inv_nsres, double&x, double&y); }; template<class T> struct Gradient<T, HORN> { static void calc(const T* afWin, double inv_ewres, double inv_nsres, double&x, double&y) { x = ((afWin[0] + afWin[3] + afWin[3] + afWin[6]) - (afWin[2] + afWin[5] + afWin[5] + afWin[8])) * inv_ewres; y = ((afWin[6] + afWin[7] + afWin[7] + afWin[8]) - (afWin[0] + afWin[1] + afWin[1] + afWin[2])) * inv_nsres; } }; template<class T> struct Gradient<T, ZEVENBERGEN_THORNE> { static void calc(const T* afWin, double inv_ewres, double inv_nsres, double&x, double&y) { x = (afWin[3] - afWin[5]) * inv_ewres; y = (afWin[7] - afWin[1]) * inv_nsres; } }; /************************************************************************/ /* GDALHillshade() */ /************************************************************************/ typedef struct { double inv_nsres; double inv_ewres; double sin_altRadians; double cos_alt_mul_z; double azRadians; double cos_az_mul_cos_alt_mul_z; double sin_az_mul_cos_alt_mul_z; double square_z; double sin_altRadians_mul_254; double cos_az_mul_cos_alt_mul_z_mul_254; double sin_az_mul_cos_alt_mul_z_mul_254; double square_z_mul_square_inv_res; double cos_az_mul_cos_alt_mul_z_mul_254_mul_inv_res; double sin_az_mul_cos_alt_mul_z_mul_254_mul_inv_res; double z_scaled; } GDALHillshadeAlgData; /* Unoptimized formulas are : x = psData->z*((afWin[0] + afWin[3] + afWin[3] + afWin[6]) - (afWin[2] + afWin[5] + afWin[5] + afWin[8])) / (8.0 * psData->ewres * psData->scale); y = psData->z*((afWin[6] + afWin[7] + afWin[7] + afWin[8]) - (afWin[0] + afWin[1] + afWin[1] + afWin[2])) / (8.0 * psData->nsres * psData->scale); slope = atan(sqrt(x*x + y*y)); aspect = atan2(y,x); cang = sin(alt) * cos(slope) + cos(alt) * sin(slope) * cos(az - M_PI/2 - aspect); We can avoid a lot of trigonometric computations: since cos(atan(x)) = 1 / sqrt(1+x^2) ==> cos(slope) = 1 / sqrt(1+ x*x+y*y) and sin(atan(x)) = x / sqrt(1+x^2) ==> sin(slope) = sqrt(x*x + y*y) / sqrt(1+ x*x+y*y) and cos(az - M_PI/2 - aspect) = cos(-az + M_PI/2 + aspect) = cos(M_PI/2 - (az - aspect)) = sin(az - aspect) = -sin(aspect-az) ==> cang = (sin(alt) - cos(alt) * sqrt(x*x + y*y) * sin(aspect-as)) / sqrt(1+ x*x+y*y) But: sin(aspect - az) = sin(aspect)*cos(az) - cos(aspect)*sin(az)) and as sin(aspect)=sin(atan2(y,x)) = y / sqrt(xx_plus_yy) and cos(aspect)=cos(atan2(y,x)) = x / sqrt(xx_plus_yy) sin(aspect - az) = (y * cos(az) - x * sin(az)) / sqrt(xx_plus_yy) so we get a final formula with just one transcendental function (reciprocal of square root): cang = (psData->sin_altRadians - (y * psData->cos_az_mul_cos_alt_mul_z - x * psData->sin_az_mul_cos_alt_mul_z)) / sqrt(1 + psData->square_z * xx_plus_yy); */ #ifdef HAVE_SSE2 inline double ApproxADivByInvSqrtB( double a, double b ) { __m128d regB = _mm_load_sd( &b ); __m128d regB_half = _mm_mul_sd( regB, _mm_set1_pd( 0.5 ) ); // Compute rough approximation of 1 / sqrt(b) with _mm_rsqrt_ss regB = _mm_cvtss_sd( regB, _mm_rsqrt_ss( _mm_cvtsd_ss( _mm_setzero_ps(), regB ) ) ); // And perform one step of Newton-Raphson approximation to improve it // approx_inv_sqrt_x = approx_inv_sqrt_x*(1.5 - // 0.5*x*approx_inv_sqrt_x*approx_inv_sqrt_x); regB = _mm_mul_sd(regB, _mm_sub_sd( _mm_set1_pd( 1.5 ), _mm_mul_sd(regB_half, _mm_mul_sd(regB, regB)) ) ); double dOut; _mm_store_sd( &dOut, regB ); return a * dOut; } #else inline double ApproxADivByInvSqrtB( double a, double b ) { return a / sqrt(b); } #endif static double NormalizeAngle (double angle, double normalizer) { angle = std::fmod(angle, normalizer); if (angle < 0) angle = normalizer + angle; return angle; } static double DifferenceBetweenAngles (double angle1, double angle2, double normalizer) { double diff = NormalizeAngle (angle1, normalizer) - NormalizeAngle (angle2, normalizer); diff = std::abs(diff); if (diff > normalizer / 2) diff = normalizer - diff; return diff; } template<class T, GradientAlg alg> static float GDALHillshadeIgorAlg (const T* afWin, float /*fDstNoDataValue*/, void* pData) { GDALHillshadeAlgData* psData = static_cast<GDALHillshadeAlgData*>(pData); double slopeDegrees; if (alg == HORN) { const double dx = ((afWin[0] + afWin[3] + afWin[3] + afWin[6]) - (afWin[2] + afWin[5] + afWin[5] + afWin[8])) * psData->inv_ewres; const double dy = ((afWin[6] + afWin[7] + afWin[7] + afWin[8]) - (afWin[0] + afWin[1] + afWin[1] + afWin[2])) * psData->inv_nsres; const double key = (dx * dx + dy * dy); slopeDegrees = atan(sqrt(key) * psData->z_scaled) * kdfRadiansToDegrees; } else // ZEVENBERGEN_THORNE { const double dx = (afWin[3] - afWin[5]) * psData->inv_ewres; const double dy = (afWin[7] - afWin[1]) * psData->inv_nsres; const double key = dx * dx + dy * dy; slopeDegrees = atan(sqrt(key) * psData->z_scaled) * kdfRadiansToDegrees; } double aspect; if (alg == HORN) { const double dx = ((afWin[2] + afWin[5] + afWin[5] + afWin[8]) - (afWin[0] + afWin[3] + afWin[3] + afWin[6])); const double dy2 = ((afWin[6] + afWin[7] + afWin[7] + afWin[8]) - (afWin[0] + afWin[1] + afWin[1] + afWin[2])); aspect = atan2(dy2, -dx); } else // ZEVENBERGEN_THORNE { const double dx = afWin[5] - afWin[3]; const double dy = afWin[7] - afWin[1]; aspect = atan2(dy, -dx); } double slopeStrength = slopeDegrees / 90; double aspectDiff = DifferenceBetweenAngles(aspect, M_PI * 3 / 2 - psData->azRadians, M_PI * 2); double aspectStrength = 1 - aspectDiff / M_PI; double shadowness = 1.0 - slopeStrength * aspectStrength; return static_cast<float>(255.0 * shadowness); } template<class T, GradientAlg alg> static float GDALHillshadeAlg (const T* afWin, float /*fDstNoDataValue*/, void* pData) { GDALHillshadeAlgData* psData = static_cast<GDALHillshadeAlgData*>(pData); // First Slope ... double x, y; Gradient<T, alg>::calc(afWin, psData->inv_ewres, psData->inv_nsres, x, y); const double xx_plus_yy = x * x + y * y; // ... then the shade value const double cang_mul_254 = ApproxADivByInvSqrtB( psData->sin_altRadians_mul_254 - (y * psData->cos_az_mul_cos_alt_mul_z_mul_254 - x * psData->sin_az_mul_cos_alt_mul_z_mul_254), 1 + psData->square_z * xx_plus_yy); const double cang = cang_mul_254 <= 0.0 ? 1.0 : 1.0 + cang_mul_254; return static_cast<float>(cang); } template<class T> static float GDALHillshadeAlg_same_res (const T* afWin, float /*fDstNoDataValue*/, void* pData) { GDALHillshadeAlgData* psData = static_cast<GDALHillshadeAlgData*>(pData); // First Slope ... /*x = (afWin[0] + afWin[3] + afWin[3] + afWin[6]) - (afWin[2] + afWin[5] + afWin[5] + afWin[8]); y = (afWin[0] + afWin[1] + afWin[1] + afWin[2]) - (afWin[6] + afWin[7] + afWin[7] + afWin[8]);*/ T accX = afWin[0] - afWin[8]; const T six_minus_two = afWin[6] - afWin[2]; T accY = accX; const T three_minus_five = afWin[3] - afWin[5]; const T one_minus_seven = afWin[1] - afWin[7]; accX += three_minus_five; accY += one_minus_seven; accX += three_minus_five; accY += one_minus_seven; accX += six_minus_two; accY -= six_minus_two; const double x = accX; const double y = accY; const double xx_plus_yy = x * x + y * y; // ... then the shade value const double cang_mul_254 = ApproxADivByInvSqrtB( psData->sin_altRadians_mul_254 + (x * psData->sin_az_mul_cos_alt_mul_z_mul_254_mul_inv_res + y * psData->cos_az_mul_cos_alt_mul_z_mul_254_mul_inv_res), 1 + psData->square_z_mul_square_inv_res * xx_plus_yy); const double cang = cang_mul_254 <= 0.0 ? 1.0 : 1.0 + cang_mul_254; return static_cast<float>(cang); } #ifdef HAVE_16_SSE_REG template<class T> static int GDALHillshadeAlg_same_res_multisample( const T* pafThreeLineWin, int nLine1Off, int nLine2Off, int nLine3Off, int nXSize, void* pData, float* pafOutputBuf ) { // Only valid for T == int GDALHillshadeAlgData* psData = static_cast<GDALHillshadeAlgData*>(pData); const __m128d reg_fact_x = _mm_load1_pd( &(psData->sin_az_mul_cos_alt_mul_z_mul_254_mul_inv_res)); const __m128d reg_fact_y = _mm_load1_pd ( &(psData->cos_az_mul_cos_alt_mul_z_mul_254_mul_inv_res)); const __m128d reg_constant_num = _mm_load1_pd( &(psData->sin_altRadians_mul_254)); const __m128d reg_constant_denom = _mm_load1_pd( &(psData->square_z_mul_square_inv_res)); const __m128d reg_half = _mm_set1_pd(0.5); const __m128d reg_one = _mm_add_pd(reg_half, reg_half); const __m128 reg_one_float = _mm_set1_ps(1); int j = 1; // Used after for. for( ; j < nXSize - 4; j+= 4 ) { const T* firstLine = pafThreeLineWin + nLine1Off + j-1; const T* secondLine = pafThreeLineWin + nLine2Off + j-1; const T* thirdLine = pafThreeLineWin + nLine3Off + j-1; __m128i firstLine0 = _mm_loadu_si128( reinterpret_cast<__m128i const*>(firstLine) ); __m128i firstLine1 = _mm_loadu_si128( reinterpret_cast<__m128i const*>(firstLine + 1) ); __m128i firstLine2 = _mm_loadu_si128( reinterpret_cast<__m128i const*>(firstLine + 2) ); __m128i thirdLine0 = _mm_loadu_si128( reinterpret_cast<__m128i const*>(thirdLine) ); __m128i thirdLine1 = _mm_loadu_si128( reinterpret_cast<__m128i const*>(thirdLine + 1) ); __m128i thirdLine2 = _mm_loadu_si128( reinterpret_cast<__m128i const*>(thirdLine + 2) ); __m128i accX = _mm_sub_epi32( firstLine0, thirdLine2); const __m128i six_minus_two = _mm_sub_epi32( thirdLine0, firstLine2 ); __m128i accY = accX; const __m128i three_minus_five = _mm_sub_epi32( _mm_loadu_si128( reinterpret_cast<__m128i const*>(secondLine) ), _mm_loadu_si128( reinterpret_cast<__m128i const*>(secondLine+2) ) ); const __m128i one_minus_seven = _mm_sub_epi32( firstLine1, thirdLine1 ); accX = _mm_add_epi32(accX, three_minus_five); accY = _mm_add_epi32(accY, one_minus_seven); accX = _mm_add_epi32(accX, three_minus_five); accY = _mm_add_epi32(accY, one_minus_seven); accX = _mm_add_epi32(accX, six_minus_two); accY = _mm_sub_epi32(accY, six_minus_two); __m128d reg_x0 = _mm_cvtepi32_pd(accX); __m128d reg_x1 = _mm_cvtepi32_pd(_mm_srli_si128(accX, 8)); __m128d reg_y0 = _mm_cvtepi32_pd(accY); __m128d reg_y1 = _mm_cvtepi32_pd(_mm_srli_si128(accY, 8)); __m128d reg_xx_plus_yy0 = _mm_add_pd( _mm_mul_pd(reg_x0, reg_x0), _mm_mul_pd(reg_y0, reg_y0) ); __m128d reg_xx_plus_yy1 = _mm_add_pd( _mm_mul_pd(reg_x1, reg_x1), _mm_mul_pd(reg_y1, reg_y1) ); __m128d reg_numerator0 = _mm_add_pd(reg_constant_num, _mm_add_pd( _mm_mul_pd(reg_fact_x, reg_x0), _mm_mul_pd(reg_fact_y, reg_y0) ) ); __m128d reg_numerator1 = _mm_add_pd(reg_constant_num, _mm_add_pd( _mm_mul_pd(reg_fact_x, reg_x1), _mm_mul_pd(reg_fact_y, reg_y1) ) ); __m128d reg_denominator0 = _mm_add_pd(reg_one, _mm_mul_pd(reg_constant_denom, reg_xx_plus_yy0)); __m128d reg_denominator1 = _mm_add_pd(reg_one, _mm_mul_pd(reg_constant_denom, reg_xx_plus_yy1)); __m128d regB0 = reg_denominator0; __m128d regB1 = reg_denominator1; __m128d regB0_half = _mm_mul_pd( regB0, reg_half ); __m128d regB1_half = _mm_mul_pd( regB1, reg_half ); // Compute rough approximation of 1 / sqrt(b) with _mm_rsqrt_ps regB0 = _mm_cvtps_pd( _mm_rsqrt_ps( _mm_cvtpd_ps( regB0 ) ) ); regB1 = _mm_cvtps_pd( _mm_rsqrt_ps( _mm_cvtpd_ps( regB1 ) ) ); // And perform one step of Newton-Raphson approximation to improve it // approx_inv_sqrt_x = approx_inv_sqrt_x*(1.5 - // 0.5*x*approx_inv_sqrt_x*approx_inv_sqrt_x); const __m128d reg_one_and_a_half = _mm_add_pd(reg_one, reg_half); regB0 = _mm_mul_pd(regB0, _mm_sub_pd( reg_one_and_a_half, _mm_mul_pd(regB0_half, _mm_mul_pd(regB0, regB0)) ) ); regB1 = _mm_mul_pd(regB1, _mm_sub_pd( reg_one_and_a_half, _mm_mul_pd(regB1_half, _mm_mul_pd(regB1, regB1)) ) ); reg_numerator0 = _mm_mul_pd(reg_numerator0, regB0); reg_numerator1 = _mm_mul_pd(reg_numerator1, regB1); __m128 res = _mm_castsi128_ps( _mm_unpacklo_epi64 (_mm_castps_si128(_mm_cvtpd_ps(reg_numerator0)), _mm_castps_si128(_mm_cvtpd_ps(reg_numerator1)))); res = _mm_add_ps(res, reg_one_float); res = _mm_max_ps(res, reg_one_float); _mm_storeu_ps( pafOutputBuf + j, res); } return j; } #endif static const double INV_SQUARE_OF_HALF_PI = 1.0 / ((M_PI*M_PI)/4); template<class T, GradientAlg alg> static float GDALHillshadeCombinedAlg (const T* afWin, float /*fDstNoDataValue*/, void* pData) { GDALHillshadeAlgData* psData = static_cast<GDALHillshadeAlgData*>(pData); // First Slope ... double x, y; Gradient<T, alg>::calc(afWin, psData->inv_ewres, psData->inv_nsres, x, y); const double xx_plus_yy = x * x + y * y; const double slope = xx_plus_yy * psData->square_z; // ... then the shade value double cang = acos( ApproxADivByInvSqrtB( psData->sin_altRadians - (y * psData->cos_az_mul_cos_alt_mul_z - x * psData->sin_az_mul_cos_alt_mul_z), 1 + slope)); // combined shading cang = 1 - cang * atan(sqrt(slope)) * INV_SQUARE_OF_HALF_PI; const float fcang = cang <= 0.0 ? 1.0f : static_cast<float>(1.0 + (254.0 * cang)); return fcang; } static void* GDALCreateHillshadeData( double* adfGeoTransform, double z, double scale, double alt, double az, bool bZevenbergenThorne ) { GDALHillshadeAlgData* pData = static_cast<GDALHillshadeAlgData *>( CPLCalloc(1, sizeof(GDALHillshadeAlgData))); pData->inv_nsres = 1.0 / adfGeoTransform[5]; pData->inv_ewres = 1.0 / adfGeoTransform[1]; pData->sin_altRadians = sin(alt * kdfDegreesToRadians); pData->azRadians = az * kdfDegreesToRadians; pData->z_scaled = z / ((bZevenbergenThorne ? 2 : 8) * scale); pData->cos_alt_mul_z = cos(alt * kdfDegreesToRadians) * pData->z_scaled; pData->cos_az_mul_cos_alt_mul_z = cos(pData->azRadians) * pData->cos_alt_mul_z; pData->sin_az_mul_cos_alt_mul_z = sin(pData->azRadians) * pData->cos_alt_mul_z; pData->square_z = pData->z_scaled * pData->z_scaled; pData->sin_altRadians_mul_254 = 254.0 * pData->sin_altRadians; pData->cos_az_mul_cos_alt_mul_z_mul_254 = 254.0 * pData->cos_az_mul_cos_alt_mul_z; pData->sin_az_mul_cos_alt_mul_z_mul_254 = 254.0 * pData->sin_az_mul_cos_alt_mul_z; if( adfGeoTransform[1] == -adfGeoTransform[5] ) { pData->square_z_mul_square_inv_res = pData->square_z * pData->inv_ewres * pData->inv_ewres; pData->cos_az_mul_cos_alt_mul_z_mul_254_mul_inv_res = pData->cos_az_mul_cos_alt_mul_z_mul_254 * -pData->inv_ewres; pData->sin_az_mul_cos_alt_mul_z_mul_254_mul_inv_res = pData->sin_az_mul_cos_alt_mul_z_mul_254 * pData->inv_ewres; } return pData; } /************************************************************************/ /* GDALHillshadeMultiDirectional() */ /************************************************************************/ typedef struct { double inv_nsres; double inv_ewres; double square_z; double sin_altRadians_mul_127; double sin_altRadians_mul_254; double cos_alt_mul_z_mul_127; double cos225_az_mul_cos_alt_mul_z_mul_127; } GDALHillshadeMultiDirectionalAlgData; template<class T, GradientAlg alg> static float GDALHillshadeMultiDirectionalAlg (const T* afWin, float /*fDstNoDataValue*/, void* pData) { const GDALHillshadeMultiDirectionalAlgData* psData = static_cast<const GDALHillshadeMultiDirectionalAlgData*>(pData); // First Slope ... double x, y; Gradient<T, alg>::calc(afWin, psData->inv_ewres, psData->inv_nsres, x, y); // See http://pubs.usgs.gov/of/1992/of92-422/of92-422.pdf // W225 = sin^2(aspect - 225) = 0.5 * (1 - 2 * sin(aspect) * cos(aspect)) // W270 = sin^2(aspect - 270) = cos^2(aspect) // W315 = sin^2(aspect - 315) = 0.5 * (1 + 2 * sin(aspect) * cos(aspect)) // W360 = sin^2(aspect - 360) = sin^2(aspect) // hillshade= 0.5 * (W225 * hillshade(az=225) + // W270 * hillshade(az=270) + // W315 * hillshade(az=315) + // W360 * hillshade(az=360)) const double xx = x * x; const double yy = y * y; const double xx_plus_yy = xx + yy; if( xx_plus_yy == 0.0 ) return static_cast<float>(1.0 + psData->sin_altRadians_mul_254); // ... then the shade value from different azimuth double val225_mul_127 = psData->sin_altRadians_mul_127 + (x-y) * psData->cos225_az_mul_cos_alt_mul_z_mul_127; val225_mul_127 = ( val225_mul_127 <= 0.0) ? 0.0 : val225_mul_127; double val270_mul_127 = psData->sin_altRadians_mul_127 - x * psData->cos_alt_mul_z_mul_127; val270_mul_127 = ( val270_mul_127 <= 0.0) ? 0.0 : val270_mul_127; double val315_mul_127 = psData->sin_altRadians_mul_127 + (x+y) * psData->cos225_az_mul_cos_alt_mul_z_mul_127; val315_mul_127 = ( val315_mul_127 <= 0.0) ? 0.0 : val315_mul_127; double val360_mul_127 = psData->sin_altRadians_mul_127 - y * psData->cos_alt_mul_z_mul_127; val360_mul_127 = ( val360_mul_127 <= 0.0) ? 0.0 : val360_mul_127; // ... then the weighted shading const double weight_225 = 0.5 * xx_plus_yy - x * y; const double weight_270 = xx; const double weight_315 = xx_plus_yy - weight_225; const double weight_360 = yy; const double cang_mul_127 = ApproxADivByInvSqrtB( (weight_225 * val225_mul_127 + weight_270 * val270_mul_127 + weight_315 * val315_mul_127 + weight_360 * val360_mul_127) / xx_plus_yy, 1 + psData->square_z * xx_plus_yy); const double cang = 1.0 + cang_mul_127; return static_cast<float>(cang); } static void* GDALCreateHillshadeMultiDirectionalData( double* adfGeoTransform, double z, double scale, double alt, bool bZevenbergenThorne ) { GDALHillshadeMultiDirectionalAlgData* pData = static_cast<GDALHillshadeMultiDirectionalAlgData *>( CPLCalloc(1, sizeof(GDALHillshadeMultiDirectionalAlgData))); pData->inv_nsres = 1.0 / adfGeoTransform[5]; pData->inv_ewres = 1.0 / adfGeoTransform[1]; const double z_scaled = z / ((bZevenbergenThorne ? 2 : 8) * scale); const double cos_alt_mul_z = cos(alt * kdfDegreesToRadians) * z_scaled; pData->square_z = z_scaled * z_scaled; pData->sin_altRadians_mul_127 = 127.0 * sin(alt * kdfDegreesToRadians); pData->sin_altRadians_mul_254 = 254.0 * sin(alt * kdfDegreesToRadians); pData->cos_alt_mul_z_mul_127 = 127.0 * cos_alt_mul_z; pData->cos225_az_mul_cos_alt_mul_z_mul_127 = 127.0 * cos(225 * kdfDegreesToRadians) * cos_alt_mul_z; return pData; } /************************************************************************/ /* GDALSlope() */ /************************************************************************/ typedef struct { double nsres; double ewres; double scale; int slopeFormat; } GDALSlopeAlgData; template<class T> static float GDALSlopeHornAlg( const T* afWin, float /*fDstNoDataValue*/, void* pData ) { const GDALSlopeAlgData* psData = static_cast<const GDALSlopeAlgData*>(pData); const double dx = ((afWin[0] + afWin[3] + afWin[3] + afWin[6]) - (afWin[2] + afWin[5] + afWin[5] + afWin[8]))/psData->ewres; const double dy = ((afWin[6] + afWin[7] + afWin[7] + afWin[8]) - (afWin[0] + afWin[1] + afWin[1] + afWin[2]))/psData->nsres; const double key = (dx * dx + dy * dy); if( psData->slopeFormat == 1 ) return static_cast<float>( atan(sqrt(key) / (8*psData->scale)) * kdfRadiansToDegrees); return static_cast<float>(100*(sqrt(key) / (8*psData->scale))); } template<class T> static float GDALSlopeZevenbergenThorneAlg( const T* afWin, float /*fDstNoDataValue*/, void* pData ) { const GDALSlopeAlgData* psData = static_cast<const GDALSlopeAlgData*>(pData); const double dx = (afWin[3] - afWin[5]) / psData->ewres; const double dy = (afWin[7] - afWin[1]) / psData->nsres; const double key = dx * dx + dy * dy; if( psData->slopeFormat == 1 ) return static_cast<float>( atan(sqrt(key) / (2*psData->scale)) * kdfRadiansToDegrees); return static_cast<float>(100*(sqrt(key) / (2*psData->scale))); } static void* GDALCreateSlopeData( double* adfGeoTransform, double scale, int slopeFormat ) { GDALSlopeAlgData* pData = static_cast<GDALSlopeAlgData*>(CPLMalloc(sizeof(GDALSlopeAlgData))); pData->nsres = adfGeoTransform[5]; pData->ewres = adfGeoTransform[1]; pData->scale = scale; pData->slopeFormat = slopeFormat; return pData; } /************************************************************************/ /* GDALAspect() */ /************************************************************************/ typedef struct { bool bAngleAsAzimuth; } GDALAspectAlgData; template<class T> static float GDALAspectAlg( const T* afWin, float fDstNoDataValue, void* pData ) { const GDALAspectAlgData* psData = static_cast<const GDALAspectAlgData*>(pData); const double dx = ((afWin[2] + afWin[5] + afWin[5] + afWin[8]) - (afWin[0] + afWin[3] + afWin[3] + afWin[6])); const double dy = ((afWin[6] + afWin[7] + afWin[7] + afWin[8]) - (afWin[0] + afWin[1] + afWin[1] + afWin[2])); float aspect = static_cast<float>(atan2(dy,-dx) / kdfDegreesToRadians); if( dx == 0 && dy == 0 ) { /* Flat area */ aspect = fDstNoDataValue; } else if( psData->bAngleAsAzimuth ) { if( aspect > 90.0f ) aspect = 450.0f - aspect; else aspect = 90.0f - aspect; } else { if( aspect < 0 ) aspect += 360.0f; } if( aspect == 360.0f ) aspect = 0.0; return aspect; } template<class T> static float GDALAspectZevenbergenThorneAlg( const T* afWin, float fDstNoDataValue, void* pData ) { const GDALAspectAlgData* psData = static_cast<const GDALAspectAlgData*>(pData); const double dx = afWin[5] - afWin[3]; const double dy = afWin[7] - afWin[1]; float aspect = static_cast<float>(atan2(dy,-dx) / kdfDegreesToRadians); if( dx == 0 && dy == 0 ) { /* Flat area */ aspect = fDstNoDataValue; } else if( psData->bAngleAsAzimuth ) { if( aspect > 90.0f ) aspect = 450.0f - aspect; else aspect = 90.0f - aspect; } else { if( aspect < 0 ) aspect += 360.0f; } if( aspect == 360.0f ) aspect = 0.0; return aspect; } static void *GDALCreateAspectData( bool bAngleAsAzimuth ) { GDALAspectAlgData* pData = static_cast<GDALAspectAlgData *>(CPLMalloc(sizeof(GDALAspectAlgData))); pData->bAngleAsAzimuth = bAngleAsAzimuth; return pData; } /************************************************************************/ /* GDALColorRelief() */ /************************************************************************/ typedef struct { double dfVal; int nR; int nG; int nB; int nA; } ColorAssociation; static int GDALColorReliefSortColors(const ColorAssociation &pA, const ColorAssociation &pB) { /* Sort NaN in first position */ return (CPLIsNan(pA.dfVal) && !CPLIsNan(pB.dfVal)) || pA.dfVal < pB.dfVal; } static void GDALColorReliefProcessColors(ColorAssociation **ppasColorAssociation, int *pnColorAssociation, int bSrcHasNoData, double dfSrcNoDataValue) { ColorAssociation *pasColorAssociation = *ppasColorAssociation; int nColorAssociation = *pnColorAssociation; std::stable_sort(pasColorAssociation, pasColorAssociation + nColorAssociation, GDALColorReliefSortColors); ColorAssociation *pPrevious = (nColorAssociation > 0) ? &pasColorAssociation[0] : nullptr; int nAdded = 0; int nRepeatedEntryIndex = 0; for( int i = 1; i < nColorAssociation; ++i ) { ColorAssociation *pCurrent = &pasColorAssociation[i]; // NaN comparison is always false, so it handles itself if( bSrcHasNoData && pCurrent->dfVal == dfSrcNoDataValue ) { // Check if there is enough distance between the nodata value and // its predecessor. const double dfNewValue = pCurrent->dfVal - std::abs(pCurrent->dfVal) * DBL_EPSILON; if( dfNewValue > pPrevious->dfVal ) { // add one just below the nodata value ++nAdded; ColorAssociation sPrevious = *pPrevious; pasColorAssociation = static_cast<ColorAssociation *>( CPLRealloc(pasColorAssociation, (nColorAssociation + nAdded) * sizeof(ColorAssociation))); pCurrent = &pasColorAssociation[i]; pPrevious = nullptr; pasColorAssociation[nColorAssociation + nAdded - 1] = sPrevious; pasColorAssociation[nColorAssociation + nAdded - 1].dfVal = dfNewValue; } } else if( bSrcHasNoData && pPrevious->dfVal == dfSrcNoDataValue ) { // Check if there is enough distance between the nodata value and // its successor. const double dfNewValue = pPrevious->dfVal + std::abs(pPrevious->dfVal) * DBL_EPSILON; if( dfNewValue < pCurrent->dfVal ) { // add one just above the nodata value ++nAdded; ColorAssociation sCurrent = *pCurrent; pasColorAssociation = static_cast<ColorAssociation *>( CPLRealloc(pasColorAssociation, (nColorAssociation + nAdded) * sizeof(ColorAssociation))); pCurrent = &pasColorAssociation[i]; pPrevious = nullptr; pasColorAssociation[nColorAssociation + nAdded - 1] = sCurrent; pasColorAssociation[nColorAssociation + nAdded - 1].dfVal = dfNewValue; } } else if( nRepeatedEntryIndex == 0 && pCurrent->dfVal == pPrevious->dfVal ) { // second of a series of equivalent entries nRepeatedEntryIndex = i; } else if( nRepeatedEntryIndex != 0 && pCurrent->dfVal != pPrevious->dfVal ) { // Get the distance between the predecessor and successor of the // equivalent entries. double dfTotalDist = 0.0; double dfLeftDist = 0.0; if( nRepeatedEntryIndex >= 2 ) { ColorAssociation *pLower = &pasColorAssociation[nRepeatedEntryIndex - 2]; dfTotalDist = pCurrent->dfVal - pLower->dfVal; dfLeftDist = pPrevious->dfVal - pLower->dfVal; } else { dfTotalDist = pCurrent->dfVal - pPrevious->dfVal; } // check if this distance is enough const int nEquivalentCount = i - nRepeatedEntryIndex + 1; if( dfTotalDist > std::abs(pPrevious->dfVal) * nEquivalentCount * DBL_EPSILON ) { // balance the alterations double dfMultiplier = 0.5 - nEquivalentCount * dfLeftDist / dfTotalDist; for( int j = nRepeatedEntryIndex - 1; j < i; ++j ) { pasColorAssociation[j].dfVal += (std::abs(pPrevious->dfVal) * dfMultiplier) * DBL_EPSILON; dfMultiplier += 1.0; } } else { // Fallback to the old behaviour: keep equivalent entries as // they are. } nRepeatedEntryIndex = 0; } pPrevious = pCurrent; } if( nAdded ) { std::stable_sort(pasColorAssociation, pasColorAssociation + nColorAssociation + nAdded, GDALColorReliefSortColors); *ppasColorAssociation = pasColorAssociation; *pnColorAssociation = nColorAssociation + nAdded; } } static bool GDALColorReliefGetRGBA( ColorAssociation* pasColorAssociation, int nColorAssociation, double dfVal, ColorSelectionMode eColorSelectionMode, int* pnR, int* pnG, int* pnB, int* pnA) { int lower = 0; // Special case for NaN if( CPLIsNan(pasColorAssociation[0].dfVal) ) { if( CPLIsNan(dfVal) ) { *pnR = pasColorAssociation[0].nR; *pnG = pasColorAssociation[0].nG; *pnB = pasColorAssociation[0].nB; *pnA = pasColorAssociation[0].nA; return true; } else { lower = 1; } } // Find the index of the first element in the LUT input array that // is not smaller than the dfVal value. int i = 0; int upper = nColorAssociation - 1; while( true ) { const int mid = (lower + upper) / 2; if( upper - lower <= 1 ) { if( dfVal <= pasColorAssociation[lower].dfVal ) i = lower; else if( dfVal <= pasColorAssociation[upper].dfVal ) i = upper; else i = upper + 1; break; } else if( pasColorAssociation[mid].dfVal >= dfVal ) { upper = mid; } else { lower = mid; } } if( i == 0 ) { if( eColorSelectionMode == COLOR_SELECTION_EXACT_ENTRY && pasColorAssociation[0].dfVal != dfVal ) { *pnR = 0; *pnG = 0; *pnB = 0; *pnA = 0; return false; } else { *pnR = pasColorAssociation[0].nR; *pnG = pasColorAssociation[0].nG; *pnB = pasColorAssociation[0].nB; *pnA = pasColorAssociation[0].nA; return true; } } else if( i == nColorAssociation ) { if( eColorSelectionMode == COLOR_SELECTION_EXACT_ENTRY && pasColorAssociation[i-1].dfVal != dfVal ) { *pnR = 0; *pnG = 0; *pnB = 0; *pnA = 0; return false; } else { *pnR = pasColorAssociation[i-1].nR; *pnG = pasColorAssociation[i-1].nG; *pnB = pasColorAssociation[i-1].nB; *pnA = pasColorAssociation[i-1].nA; return true; } } else { if( eColorSelectionMode == COLOR_SELECTION_EXACT_ENTRY && pasColorAssociation[i-1].dfVal != dfVal ) { *pnR = 0; *pnG = 0; *pnB = 0; *pnA = 0; return false; } if( eColorSelectionMode == COLOR_SELECTION_NEAREST_ENTRY && pasColorAssociation[i-1].dfVal != dfVal ) { int index = i; if( dfVal - pasColorAssociation[i-1].dfVal < pasColorAssociation[i].dfVal - dfVal ) { --index; } *pnR = pasColorAssociation[index].nR; *pnG = pasColorAssociation[index].nG; *pnB = pasColorAssociation[index].nB; *pnA = pasColorAssociation[index].nA; return true; } if( pasColorAssociation[i-1].dfVal == dfVal ) { *pnR = pasColorAssociation[i-1].nR; *pnG = pasColorAssociation[i-1].nG; *pnB = pasColorAssociation[i-1].nB; *pnA = pasColorAssociation[i-1].nA; return true; } if( CPLIsNan(pasColorAssociation[i-1].dfVal) ) { *pnR = pasColorAssociation[i].nR; *pnG = pasColorAssociation[i].nG; *pnB = pasColorAssociation[i].nB; *pnA = pasColorAssociation[i].nA; return true; } const double dfRatio = (dfVal - pasColorAssociation[i-1].dfVal) / (pasColorAssociation[i].dfVal - pasColorAssociation[i-1].dfVal); *pnR = static_cast<int>( 0.45 + pasColorAssociation[i-1].nR + dfRatio * (pasColorAssociation[i].nR - pasColorAssociation[i-1].nR)); if( *pnR < 0 ) *pnR = 0; else if( *pnR > 255 ) *pnR = 255; *pnG = static_cast<int>( 0.45 + pasColorAssociation[i-1].nG + dfRatio * (pasColorAssociation[i].nG - pasColorAssociation[i-1].nG)); if( *pnG < 0 ) *pnG = 0; else if( *pnG > 255 ) *pnG = 255; *pnB = static_cast<int>( 0.45 + pasColorAssociation[i-1].nB + dfRatio * (pasColorAssociation[i].nB - pasColorAssociation[i-1].nB)); if( *pnB < 0 ) *pnB = 0; else if( *pnB > 255 ) *pnB = 255; *pnA = static_cast<int>( 0.45 + pasColorAssociation[i-1].nA + dfRatio * (pasColorAssociation[i].nA - pasColorAssociation[i-1].nA)); if( *pnA < 0 ) *pnA = 0; else if( *pnA > 255 ) *pnA = 255; return true; } } /* dfPct : percentage between 0 and 1 */ static double GDALColorReliefGetAbsoluteValFromPct( GDALRasterBandH hSrcBand, double dfPct ) { int bSuccessMin = FALSE; int bSuccessMax = FALSE; double dfMin = GDALGetRasterMinimum(hSrcBand, &bSuccessMin); double dfMax = GDALGetRasterMaximum(hSrcBand, &bSuccessMax); if( !bSuccessMin || !bSuccessMax ) { double dfMean = 0.0; double dfStdDev = 0.0; fprintf(stderr, "Computing source raster statistics...\n"); GDALComputeRasterStatistics(hSrcBand, FALSE, &dfMin, &dfMax, &dfMean, &dfStdDev, nullptr, nullptr); } return dfMin + dfPct * (dfMax - dfMin); } typedef struct { const char *name; float r, g, b; } NamedColor; static const NamedColor namedColors[] = { { "white", 1.00, 1.00, 1.00 }, { "black", 0.00, 0.00, 0.00 }, { "red", 1.00, 0.00, 0.00 }, { "green", 0.00, 1.00, 0.00 }, { "blue", 0.00, 0.00, 1.00 }, { "yellow", 1.00, 1.00, 0.00 }, { "magenta",1.00, 0.00, 1.00 }, { "cyan", 0.00, 1.00, 1.00 }, { "aqua", 0.00, 0.75, 0.75 }, { "grey", 0.75, 0.75, 0.75 }, { "gray", 0.75, 0.75, 0.75 }, { "orange", 1.00, 0.50, 0.00 }, { "brown", 0.75, 0.50, 0.25 }, { "purple", 0.50, 0.00, 1.00 }, { "violet", 0.50, 0.00, 1.00 }, { "indigo", 0.00, 0.50, 1.00 }, }; static bool GDALColorReliefFindNamedColor( const char *pszColorName, int *pnR, int *pnG, int *pnB ) { *pnR = 0; *pnG = 0; *pnB = 0; for( unsigned int i = 0; i < sizeof(namedColors) / sizeof(namedColors[0]); i++ ) { if( EQUAL(pszColorName, namedColors[i].name) ) { *pnR = static_cast<int>(255.0 * namedColors[i].r); *pnG = static_cast<int>(255.0 * namedColors[i].g); *pnB = static_cast<int>(255.0 * namedColors[i].b); return true; } } return false; } static ColorAssociation* GDALColorReliefParseColorFile( GDALRasterBandH hSrcBand, const char* pszColorFilename, int* pnColors ) { VSILFILE* fpColorFile = VSIFOpenL(pszColorFilename, "rt"); if( fpColorFile == nullptr ) { CPLError(CE_Failure, CPLE_AppDefined, "Cannot find %s", pszColorFilename); *pnColors = 0; return nullptr; } ColorAssociation* pasColorAssociation = nullptr; int nColorAssociation = 0; int bSrcHasNoData = FALSE; double dfSrcNoDataValue = GDALGetRasterNoDataValue(hSrcBand, &bSrcHasNoData); const char* pszLine = nullptr; bool bIsGMT_CPT = false; while( (pszLine = CPLReadLineL(fpColorFile)) != nullptr ) { if( pszLine[0] == '#' && strstr(pszLine, "COLOR_MODEL") ) { if( strstr(pszLine, "COLOR_MODEL = RGB") == nullptr ) { CPLError(CE_Failure, CPLE_AppDefined, "Only COLOR_MODEL = RGB is supported"); CPLFree(pasColorAssociation); *pnColors = 0; return nullptr; } bIsGMT_CPT = true; } char** papszFields = CSLTokenizeStringComplex(pszLine, " ,\t:", FALSE, FALSE ); /* Skip comment lines */ const int nTokens = CSLCount(papszFields); if( nTokens >= 1 && (papszFields[0][0] == '#' || papszFields[0][0] == '/') ) { CSLDestroy(papszFields); continue; } if( bIsGMT_CPT && nTokens == 8 ) { pasColorAssociation = static_cast<ColorAssociation *>( CPLRealloc(pasColorAssociation, (nColorAssociation + 2) * sizeof(ColorAssociation))); pasColorAssociation[nColorAssociation].dfVal = CPLAtof(papszFields[0]); pasColorAssociation[nColorAssociation].nR = atoi(papszFields[1]); pasColorAssociation[nColorAssociation].nG = atoi(papszFields[2]); pasColorAssociation[nColorAssociation].nB = atoi(papszFields[3]); pasColorAssociation[nColorAssociation].nA = 255; nColorAssociation++; pasColorAssociation[nColorAssociation].dfVal = CPLAtof(papszFields[4]); pasColorAssociation[nColorAssociation].nR = atoi(papszFields[5]); pasColorAssociation[nColorAssociation].nG = atoi(papszFields[6]); pasColorAssociation[nColorAssociation].nB = atoi(papszFields[7]); pasColorAssociation[nColorAssociation].nA = 255; nColorAssociation++; } else if( bIsGMT_CPT && nTokens == 4 ) { // The first token might be B (background), F (foreground) or N // (nodata) Just interested in N. if( EQUAL(papszFields[0], "N") && bSrcHasNoData ) { pasColorAssociation = static_cast<ColorAssociation *>( CPLRealloc(pasColorAssociation, (nColorAssociation + 1) * sizeof(ColorAssociation))); pasColorAssociation[nColorAssociation].dfVal = dfSrcNoDataValue; pasColorAssociation[nColorAssociation].nR = atoi(papszFields[1]); pasColorAssociation[nColorAssociation].nG = atoi(papszFields[2]); pasColorAssociation[nColorAssociation].nB = atoi(papszFields[3]); pasColorAssociation[nColorAssociation].nA = 255; nColorAssociation++; } } else if( !bIsGMT_CPT && nTokens >= 2 ) { pasColorAssociation = static_cast<ColorAssociation *>( CPLRealloc(pasColorAssociation, (nColorAssociation + 1) * sizeof(ColorAssociation))); if( EQUAL(papszFields[0], "nv") && bSrcHasNoData ) { pasColorAssociation[nColorAssociation].dfVal = dfSrcNoDataValue; } else if( strlen(papszFields[0]) > 1 && papszFields[0][strlen(papszFields[0])-1] == '%' ) { const double dfPct = CPLAtof(papszFields[0]) / 100.0; if( dfPct < 0.0 || dfPct > 1.0 ) { CPLError(CE_Failure, CPLE_AppDefined, "Wrong value for a percentage : %s", papszFields[0]); CSLDestroy(papszFields); VSIFCloseL(fpColorFile); CPLFree(pasColorAssociation); *pnColors = 0; return nullptr; } pasColorAssociation[nColorAssociation].dfVal = GDALColorReliefGetAbsoluteValFromPct(hSrcBand, dfPct); } else { pasColorAssociation[nColorAssociation].dfVal = CPLAtof(papszFields[0]); } if( nTokens >= 4 ) { pasColorAssociation[nColorAssociation].nR = atoi(papszFields[1]); pasColorAssociation[nColorAssociation].nG = atoi(papszFields[2]); pasColorAssociation[nColorAssociation].nB = atoi(papszFields[3]); pasColorAssociation[nColorAssociation].nA = (CSLCount(papszFields) >= 5 ) ? atoi(papszFields[4]) : 255; } else { int nR = 0; int nG = 0; int nB = 0; if( !GDALColorReliefFindNamedColor(papszFields[1], &nR, &nG, &nB) ) { CPLError(CE_Failure, CPLE_AppDefined, "Unknown color : %s", papszFields[1]); CSLDestroy(papszFields); VSIFCloseL(fpColorFile); CPLFree(pasColorAssociation); *pnColors = 0; return nullptr; } pasColorAssociation[nColorAssociation].nR = nR; pasColorAssociation[nColorAssociation].nG = nG; pasColorAssociation[nColorAssociation].nB = nB; pasColorAssociation[nColorAssociation].nA = (CSLCount(papszFields) >= 3 ) ? atoi(papszFields[2]) : 255; } nColorAssociation++; } CSLDestroy(papszFields); } VSIFCloseL(fpColorFile); if( nColorAssociation == 0 ) { CPLError(CE_Failure, CPLE_AppDefined, "No color association found in %s", pszColorFilename); *pnColors = 0; return nullptr; } GDALColorReliefProcessColors(&pasColorAssociation, &nColorAssociation, bSrcHasNoData, dfSrcNoDataValue); *pnColors = nColorAssociation; return pasColorAssociation; } static GByte* GDALColorReliefPrecompute(GDALRasterBandH hSrcBand, ColorAssociation* pasColorAssociation, int nColorAssociation, ColorSelectionMode eColorSelectionMode, int* pnIndexOffset) { const GDALDataType eDT = GDALGetRasterDataType(hSrcBand); GByte* pabyPrecomputed = nullptr; const int nIndexOffset = (eDT == GDT_Int16) ? 32768 : 0; *pnIndexOffset = nIndexOffset; const int nXSize = GDALGetRasterBandXSize(hSrcBand); const int nYSize = GDALGetRasterBandXSize(hSrcBand); if( eDT == GDT_Byte || ((eDT == GDT_Int16 || eDT == GDT_UInt16) && nXSize * nYSize > 65536) ) { const int iMax = (eDT == GDT_Byte) ? 256: 65536; pabyPrecomputed = static_cast<GByte *>(VSI_MALLOC2_VERBOSE(4, iMax)); if( pabyPrecomputed ) { for( int i = 0; i < iMax; i++ ) { int nR = 0; int nG = 0; int nB = 0; int nA = 0; GDALColorReliefGetRGBA (pasColorAssociation, nColorAssociation, i - nIndexOffset, eColorSelectionMode, &nR, &nG, &nB, &nA); pabyPrecomputed[4 * i] = static_cast<GByte>(nR); pabyPrecomputed[4 * i + 1] = static_cast<GByte>(nG); pabyPrecomputed[4 * i + 2] = static_cast<GByte>(nB); pabyPrecomputed[4 * i + 3] = static_cast<GByte>(nA); } } } return pabyPrecomputed; } /************************************************************************/ /* ==================================================================== */ /* GDALColorReliefDataset */ /* ==================================================================== */ /************************************************************************/ class GDALColorReliefRasterBand; class GDALColorReliefDataset : public GDALDataset { friend class GDALColorReliefRasterBand; GDALDatasetH hSrcDS; GDALRasterBandH hSrcBand; int nColorAssociation; ColorAssociation* pasColorAssociation; ColorSelectionMode eColorSelectionMode; GByte* pabyPrecomputed; int nIndexOffset; float* pafSourceBuf; int* panSourceBuf; int nCurBlockXOff; int nCurBlockYOff; public: GDALColorReliefDataset( GDALDatasetH hSrcDS, GDALRasterBandH hSrcBand, const char* pszColorFilename, ColorSelectionMode eColorSelectionMode, int bAlpha); ~GDALColorReliefDataset(); bool InitOK() const { return pafSourceBuf != nullptr || panSourceBuf != nullptr; } CPLErr GetGeoTransform( double * padfGeoTransform ) override; const OGRSpatialReference* GetSpatialRef() const override; }; /************************************************************************/ /* ==================================================================== */ /* GDALColorReliefRasterBand */ /* ==================================================================== */ /************************************************************************/ class GDALColorReliefRasterBand : public GDALRasterBand { friend class GDALColorReliefDataset; public: GDALColorReliefRasterBand( GDALColorReliefDataset *, int ); virtual CPLErr IReadBlock( int, int, void * ) override; virtual GDALColorInterp GetColorInterpretation() override; }; GDALColorReliefDataset::GDALColorReliefDataset( GDALDatasetH hSrcDSIn, GDALRasterBandH hSrcBandIn, const char* pszColorFilename, ColorSelectionMode eColorSelectionModeIn, int bAlpha ) : hSrcDS(hSrcDSIn), hSrcBand(hSrcBandIn), nColorAssociation(0), pasColorAssociation(nullptr), eColorSelectionMode(eColorSelectionModeIn), pabyPrecomputed(nullptr), nIndexOffset(0), pafSourceBuf(nullptr), panSourceBuf(nullptr), nCurBlockXOff(-1), nCurBlockYOff(-1) { pasColorAssociation = GDALColorReliefParseColorFile(hSrcBand, pszColorFilename, &nColorAssociation); nRasterXSize = GDALGetRasterXSize(hSrcDS); nRasterYSize = GDALGetRasterYSize(hSrcDS); int nBlockXSize = 0; int nBlockYSize = 0; GDALGetBlockSize( hSrcBand, &nBlockXSize, &nBlockYSize); pabyPrecomputed = GDALColorReliefPrecompute(hSrcBand, pasColorAssociation, nColorAssociation, eColorSelectionMode, &nIndexOffset); for( int i = 0; i < ((bAlpha) ? 4 : 3); i++ ) { SetBand(i + 1, new GDALColorReliefRasterBand(this, i+1)); } if( pabyPrecomputed ) panSourceBuf = static_cast<int *>( VSI_MALLOC3_VERBOSE(sizeof(int), nBlockXSize, nBlockYSize)); else pafSourceBuf = static_cast<float *>( VSI_MALLOC3_VERBOSE(sizeof(float), nBlockXSize, nBlockYSize)); } GDALColorReliefDataset::~GDALColorReliefDataset() { CPLFree(pasColorAssociation); CPLFree(pabyPrecomputed); CPLFree(panSourceBuf); CPLFree(pafSourceBuf); } CPLErr GDALColorReliefDataset::GetGeoTransform( double * padfGeoTransform ) { return GDALGetGeoTransform(hSrcDS, padfGeoTransform); } const OGRSpatialReference *GDALColorReliefDataset::GetSpatialRef() const { return GDALDataset::FromHandle(hSrcDS)->GetSpatialRef(); } GDALColorReliefRasterBand::GDALColorReliefRasterBand( GDALColorReliefDataset * poDSIn, int nBandIn) { poDS = poDSIn; nBand = nBandIn; eDataType = GDT_Byte; GDALGetBlockSize( poDSIn->hSrcBand, &nBlockXSize, &nBlockYSize); } CPLErr GDALColorReliefRasterBand::IReadBlock( int nBlockXOff, int nBlockYOff, void *pImage ) { GDALColorReliefDataset * poGDS = cpl::down_cast<GDALColorReliefDataset *>(poDS); const int nReqXSize = (nBlockXOff + 1) * nBlockXSize >= nRasterXSize ? nRasterXSize - nBlockXOff * nBlockXSize : nBlockXSize; const int nReqYSize = (nBlockYOff + 1) * nBlockYSize >= nRasterYSize ? nRasterYSize - nBlockYOff * nBlockYSize : nBlockYSize; if( poGDS->nCurBlockXOff != nBlockXOff || poGDS->nCurBlockYOff != nBlockYOff ) { poGDS->nCurBlockXOff = nBlockXOff; poGDS->nCurBlockYOff = nBlockYOff; const CPLErr eErr = GDALRasterIO( poGDS->hSrcBand, GF_Read, nBlockXOff * nBlockXSize, nBlockYOff * nBlockYSize, nReqXSize, nReqYSize, (poGDS->panSourceBuf) ? static_cast<void*>(poGDS->panSourceBuf) : static_cast<void*>(poGDS->pafSourceBuf), nReqXSize, nReqYSize, (poGDS->panSourceBuf) ? GDT_Int32 : GDT_Float32, 0, 0); if( eErr != CE_None ) { memset(pImage, 0, nBlockXSize * nBlockYSize); return eErr; } } int j = 0; if( poGDS->panSourceBuf ) { for( int y = 0; y < nReqYSize; y++ ) { for( int x = 0; x < nReqXSize; x++ ) { const int nIndex = poGDS->panSourceBuf[j] + poGDS->nIndexOffset; static_cast<GByte*>(pImage)[y * nBlockXSize + x] = poGDS->pabyPrecomputed[4*nIndex + nBand-1]; j++; } } } else { int anComponents[4] = { 0, 0, 0, 0 }; for( int y = 0; y < nReqYSize; y++ ) { for( int x = 0; x < nReqXSize; x++ ) { GDALColorReliefGetRGBA (poGDS->pasColorAssociation, poGDS->nColorAssociation, poGDS->pafSourceBuf[j], poGDS->eColorSelectionMode, &anComponents[0], &anComponents[1], &anComponents[2], &anComponents[3]); static_cast<GByte*>(pImage)[y * nBlockXSize + x] = static_cast<GByte>(anComponents[nBand-1]); j++; } } } return CE_None; } GDALColorInterp GDALColorReliefRasterBand::GetColorInterpretation() { return static_cast<GDALColorInterp>(GCI_RedBand + nBand - 1); } static CPLErr GDALColorRelief( GDALRasterBandH hSrcBand, GDALRasterBandH hDstBand1, GDALRasterBandH hDstBand2, GDALRasterBandH hDstBand3, GDALRasterBandH hDstBand4, const char* pszColorFilename, ColorSelectionMode eColorSelectionMode, GDALProgressFunc pfnProgress, void * pProgressData ) { if( hSrcBand == nullptr || hDstBand1 == nullptr || hDstBand2 == nullptr || hDstBand3 == nullptr ) return CE_Failure; int nColorAssociation = 0; ColorAssociation* pasColorAssociation = GDALColorReliefParseColorFile(hSrcBand, pszColorFilename, &nColorAssociation); if( pasColorAssociation == nullptr ) return CE_Failure; if( pfnProgress == nullptr ) pfnProgress = GDALDummyProgress; /* -------------------------------------------------------------------- */ /* Precompute the map from values to RGBA quadruplets */ /* for GDT_Byte, GDT_Int16 or GDT_UInt16 */ /* -------------------------------------------------------------------- */ int nIndexOffset = 0; GByte* pabyPrecomputed = GDALColorReliefPrecompute(hSrcBand, pasColorAssociation, nColorAssociation, eColorSelectionMode, &nIndexOffset); /* -------------------------------------------------------------------- */ /* Initialize progress counter. */ /* -------------------------------------------------------------------- */ const int nXSize = GDALGetRasterBandXSize(hSrcBand); const int nYSize = GDALGetRasterBandYSize(hSrcBand); float* pafSourceBuf = nullptr; int* panSourceBuf = nullptr; if( pabyPrecomputed ) panSourceBuf = static_cast<int *>( VSI_MALLOC2_VERBOSE(sizeof(int), nXSize)); else pafSourceBuf = static_cast<float *>( VSI_MALLOC2_VERBOSE(sizeof(float), nXSize)); GByte* pabyDestBuf1 = static_cast<GByte *>(VSI_MALLOC2_VERBOSE(4, nXSize)); GByte* pabyDestBuf2 = pabyDestBuf1 ? pabyDestBuf1 + nXSize : nullptr; GByte* pabyDestBuf3 = pabyDestBuf2 ? pabyDestBuf2 + nXSize : nullptr; GByte* pabyDestBuf4 = pabyDestBuf3 ? pabyDestBuf3 + nXSize : nullptr; if( (pabyPrecomputed != nullptr && panSourceBuf == nullptr) || (pabyPrecomputed == nullptr && pafSourceBuf == nullptr) || pabyDestBuf1 == nullptr ) { VSIFree(pabyPrecomputed); CPLFree(pafSourceBuf); CPLFree(panSourceBuf); CPLFree(pabyDestBuf1); CPLFree(pasColorAssociation); return CE_Failure; } if( !pfnProgress( 0.0, nullptr, pProgressData ) ) { CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" ); VSIFree(pabyPrecomputed); CPLFree(pafSourceBuf); CPLFree(panSourceBuf); CPLFree(pabyDestBuf1); CPLFree(pasColorAssociation); return CE_Failure; } int nR = 0; int nG = 0; int nB = 0; int nA = 0; for( int i = 0; i < nYSize; i++ ) { /* Read source buffer */ CPLErr eErr = GDALRasterIO( hSrcBand, GF_Read, 0, i, nXSize, 1, panSourceBuf ? static_cast<void*>(panSourceBuf) : static_cast<void*>(pafSourceBuf), nXSize, 1, panSourceBuf ? GDT_Int32 : GDT_Float32, 0, 0); if( eErr != CE_None ) { VSIFree(pabyPrecomputed); CPLFree(pafSourceBuf); CPLFree(panSourceBuf); CPLFree(pabyDestBuf1); CPLFree(pasColorAssociation); return eErr; } if( pabyPrecomputed ) { for( int j = 0; j < nXSize; j++ ) { int nIndex = panSourceBuf[j] + nIndexOffset; pabyDestBuf1[j] = pabyPrecomputed[4 * nIndex]; pabyDestBuf2[j] = pabyPrecomputed[4 * nIndex + 1]; pabyDestBuf3[j] = pabyPrecomputed[4 * nIndex + 2]; pabyDestBuf4[j] = pabyPrecomputed[4 * nIndex + 3]; } } else { for( int j = 0; j < nXSize; j++ ) { GDALColorReliefGetRGBA (pasColorAssociation, nColorAssociation, pafSourceBuf[j], eColorSelectionMode, &nR, &nG, &nB, &nA); pabyDestBuf1[j] = static_cast<GByte>(nR); pabyDestBuf2[j] = static_cast<GByte>(nG); pabyDestBuf3[j] = static_cast<GByte>(nB); pabyDestBuf4[j] = static_cast<GByte>(nA); } } /* ----------------------------------------- * Write Line to Raster */ eErr = GDALRasterIO(hDstBand1, GF_Write, 0, i, nXSize, 1, pabyDestBuf1, nXSize, 1, GDT_Byte, 0, 0); if( eErr != CE_None ) { VSIFree(pabyPrecomputed); CPLFree(pafSourceBuf); CPLFree(panSourceBuf); CPLFree(pabyDestBuf1); CPLFree(pasColorAssociation); return eErr; } eErr = GDALRasterIO(hDstBand2, GF_Write, 0, i, nXSize, 1, pabyDestBuf2, nXSize, 1, GDT_Byte, 0, 0); if( eErr != CE_None ) { VSIFree(pabyPrecomputed); CPLFree(pafSourceBuf); CPLFree(panSourceBuf); CPLFree(pabyDestBuf1); CPLFree(pasColorAssociation); return eErr; } eErr = GDALRasterIO(hDstBand3, GF_Write, 0, i, nXSize, 1, pabyDestBuf3, nXSize, 1, GDT_Byte, 0, 0); if( eErr != CE_None ) { VSIFree(pabyPrecomputed); CPLFree(pafSourceBuf); CPLFree(panSourceBuf); CPLFree(pabyDestBuf1); CPLFree(pasColorAssociation); return eErr; } if( hDstBand4 ) { eErr = GDALRasterIO(hDstBand4, GF_Write, 0, i, nXSize, 1, pabyDestBuf4, nXSize, 1, GDT_Byte, 0, 0); if( eErr != CE_None ) { VSIFree(pabyPrecomputed); CPLFree(pafSourceBuf); CPLFree(panSourceBuf); CPLFree(pabyDestBuf1); CPLFree(pasColorAssociation); return eErr; } } if( !pfnProgress( 1.0 * (i+1) / nYSize, nullptr, pProgressData ) ) { CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" ); VSIFree(pabyPrecomputed); CPLFree(pafSourceBuf); CPLFree(panSourceBuf); CPLFree(pabyDestBuf1); CPLFree(pasColorAssociation); return CE_Failure; } } pfnProgress( 1.0, nullptr, pProgressData ); VSIFree(pabyPrecomputed); CPLFree(pafSourceBuf); CPLFree(panSourceBuf); CPLFree(pabyDestBuf1); CPLFree(pasColorAssociation); return CE_None; } /************************************************************************/ /* GDALGenerateVRTColorRelief() */ /************************************************************************/ static CPLErr GDALGenerateVRTColorRelief( const char* pszDstFilename, GDALDatasetH hSrcDataset, GDALRasterBandH hSrcBand, const char* pszColorFilename, ColorSelectionMode eColorSelectionMode, bool bAddAlpha ) { int nColorAssociation = 0; ColorAssociation* pasColorAssociation = GDALColorReliefParseColorFile(hSrcBand, pszColorFilename, &nColorAssociation); if( pasColorAssociation == nullptr ) return CE_Failure; VSILFILE* fp = VSIFOpenL(pszDstFilename, "wt"); if( fp == nullptr ) { CPLFree(pasColorAssociation); return CE_Failure; } const int nXSize = GDALGetRasterBandXSize(hSrcBand); const int nYSize = GDALGetRasterBandYSize(hSrcBand); bool bOK = VSIFPrintfL(fp, "<VRTDataset rasterXSize=\"%d\" rasterYSize=\"%d\">\n", nXSize, nYSize) > 0; const char* pszProjectionRef = GDALGetProjectionRef(hSrcDataset); if( pszProjectionRef && pszProjectionRef[0] != '\0' ) { char* pszEscapedString = CPLEscapeString(pszProjectionRef, -1, CPLES_XML); bOK &= VSIFPrintfL(fp, " <SRS>%s</SRS>\n", pszEscapedString) > 0; VSIFree(pszEscapedString); } double adfGT[6] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 }; if( GDALGetGeoTransform(hSrcDataset, adfGT) == CE_None ) { bOK &= VSIFPrintfL(fp, " <GeoTransform> %.16g, %.16g, %.16g, " "%.16g, %.16g, %.16g</GeoTransform>\n", adfGT[0], adfGT[1], adfGT[2], adfGT[3], adfGT[4], adfGT[5]) > 0; } int nBlockXSize = 0; int nBlockYSize = 0; GDALGetBlockSize(hSrcBand, &nBlockXSize, &nBlockYSize); int bRelativeToVRT = FALSE; CPLString osPath = CPLGetPath(pszDstFilename); char* pszSourceFilename = CPLStrdup( CPLExtractRelativePath( osPath.c_str(), GDALGetDescription(hSrcDataset), &bRelativeToVRT )); const int nBands = 3 + (bAddAlpha ? 1 : 0); for( int iBand = 0; iBand < nBands; iBand++ ) { bOK &= VSIFPrintfL( fp, " <VRTRasterBand dataType=\"Byte\" band=\"%d\">\n", iBand + 1) > 0; bOK &= VSIFPrintfL( fp, " <ColorInterp>%s</ColorInterp>\n", GDALGetColorInterpretationName(static_cast<GDALColorInterp>(GCI_RedBand + iBand))) > 0; bOK &= VSIFPrintfL( fp, " <ComplexSource>\n") > 0; bOK &= VSIFPrintfL( fp, " <SourceFilename relativeToVRT=\"%d\">%s</SourceFilename>\n", bRelativeToVRT, pszSourceFilename) > 0; bOK &= VSIFPrintfL( fp, " <SourceBand>%d</SourceBand>\n", GDALGetBandNumber(hSrcBand)) > 0; bOK &= VSIFPrintfL( fp, " <SourceProperties RasterXSize=\"%d\" " "RasterYSize=\"%d\" DataType=\"%s\" " "BlockXSize=\"%d\" BlockYSize=\"%d\"/>\n", nXSize, nYSize, GDALGetDataTypeName(GDALGetRasterDataType(hSrcBand)), nBlockXSize, nBlockYSize) > 0; bOK &= VSIFPrintfL( fp, " " "<SrcRect xOff=\"0\" yOff=\"0\" xSize=\"%d\" ySize=\"%d\"/>\n", nXSize, nYSize) > 0; bOK &= VSIFPrintfL( fp, " " "<DstRect xOff=\"0\" yOff=\"0\" xSize=\"%d\" ySize=\"%d\"/>\n", nXSize, nYSize) > 0; bOK &= VSIFPrintfL(fp, " <LUT>") > 0; for( int iColor = 0; iColor < nColorAssociation; iColor++ ) { if( eColorSelectionMode == COLOR_SELECTION_NEAREST_ENTRY ) { if( iColor > 1 ) bOK &= VSIFPrintfL(fp, ",") > 0; } else if( iColor > 0 ) bOK &= VSIFPrintfL(fp, ",") > 0; const double dfVal = pasColorAssociation[iColor].dfVal; if( eColorSelectionMode == COLOR_SELECTION_EXACT_ENTRY ) { bOK &= VSIFPrintfL(fp, "%.18g:0,", dfVal - fabs(dfVal) * DBL_EPSILON) > 0; } else if( iColor > 0 && eColorSelectionMode == COLOR_SELECTION_NEAREST_ENTRY ) { const double dfMidVal = (dfVal + pasColorAssociation[iColor-1].dfVal) / 2.0; bOK &= VSIFPrintfL( fp, "%.18g:%d", dfMidVal - fabs(dfMidVal) * DBL_EPSILON, (iBand == 0) ? pasColorAssociation[iColor-1].nR : (iBand == 1) ? pasColorAssociation[iColor-1].nG : (iBand == 2) ? pasColorAssociation[iColor-1].nB : pasColorAssociation[iColor-1].nA) > 0; bOK &= VSIFPrintfL( fp, ",%.18g:%d", dfMidVal, (iBand == 0) ? pasColorAssociation[iColor].nR : (iBand == 1) ? pasColorAssociation[iColor].nG : (iBand == 2) ? pasColorAssociation[iColor].nB : pasColorAssociation[iColor].nA) > 0; } if( eColorSelectionMode != COLOR_SELECTION_NEAREST_ENTRY ) { if( dfVal != static_cast<double>(static_cast<int>(dfVal)) ) bOK &= VSIFPrintfL(fp, "%.18g", dfVal) > 0; else bOK &= VSIFPrintfL(fp, "%d", static_cast<int>(dfVal)) > 0; bOK &= VSIFPrintfL(fp, ":%d", (iBand == 0) ? pasColorAssociation[iColor].nR : (iBand == 1) ? pasColorAssociation[iColor].nG : (iBand == 2) ? pasColorAssociation[iColor].nB : pasColorAssociation[iColor].nA) > 0; } if( eColorSelectionMode == COLOR_SELECTION_EXACT_ENTRY ) { bOK &= VSIFPrintfL(fp, ",%.18g:0", dfVal + fabs(dfVal) * DBL_EPSILON) > 0; } } bOK &= VSIFPrintfL(fp, "</LUT>\n") > 0; bOK &= VSIFPrintfL(fp, " </ComplexSource>\n") > 0; bOK &= VSIFPrintfL(fp, " </VRTRasterBand>\n") > 0; } CPLFree(pszSourceFilename); bOK &= VSIFPrintfL(fp, "</VRTDataset>\n") > 0; VSIFCloseL(fp); CPLFree(pasColorAssociation); return (bOK) ? CE_None : CE_Failure; } /************************************************************************/ /* GDALTRIAlg() */ /************************************************************************/ template<class T> static T MyAbs(T x); template<> float MyAbs( float x ) { return fabsf(x); } template<> int MyAbs( int x ) { return x >= 0 ? x : -x; } template<class T> static float GDALTRIAlg( const T* afWin, float /*fDstNoDataValue*/, void* /*pData*/ ) { // Terrain Ruggedness is average difference in height return (MyAbs(afWin[0]-afWin[4]) + MyAbs(afWin[1]-afWin[4]) + MyAbs(afWin[2]-afWin[4]) + MyAbs(afWin[3]-afWin[4]) + MyAbs(afWin[5]-afWin[4]) + MyAbs(afWin[6]-afWin[4]) + MyAbs(afWin[7]-afWin[4]) + MyAbs(afWin[8]-afWin[4])) * 0.125f; } /************************************************************************/ /* GDALTPIAlg() */ /************************************************************************/ template<class T> static float GDALTPIAlg( const T* afWin, float /*fDstNoDataValue*/, void* /*pData*/ ) { // Terrain Position is the difference between // The central cell and the mean of the surrounding cells return afWin[4] - ((afWin[0]+ afWin[1]+ afWin[2]+ afWin[3]+ afWin[5]+ afWin[6]+ afWin[7]+ afWin[8]) * 0.125f ); } /************************************************************************/ /* GDALRoughnessAlg() */ /************************************************************************/ template<class T> static float GDALRoughnessAlg( const T* afWin, float /*fDstNoDataValue*/, void* /*pData*/ ) { // Roughness is the largest difference // between any two cells T pafRoughnessMin = afWin[0]; T pafRoughnessMax = afWin[0]; for( int k = 1; k < 9; k++ ) { if( afWin[k] > pafRoughnessMax ) { pafRoughnessMax=afWin[k]; } if( afWin[k] < pafRoughnessMin ) { pafRoughnessMin=afWin[k]; } } return static_cast<float>(pafRoughnessMax - pafRoughnessMin); } /************************************************************************/ /* ==================================================================== */ /* GDALGeneric3x3Dataset */ /* ==================================================================== */ /************************************************************************/ template<class T> class GDALGeneric3x3RasterBand; template<class T> class GDALGeneric3x3Dataset : public GDALDataset { friend class GDALGeneric3x3RasterBand<T>; typename GDALGeneric3x3ProcessingAlg<T>::type pfnAlg; void* pAlgData; GDALDatasetH hSrcDS; GDALRasterBandH hSrcBand; T* apafSourceBuf[3]; int bDstHasNoData; double dfDstNoDataValue; int nCurLine; bool bComputeAtEdges; public: GDALGeneric3x3Dataset( GDALDatasetH hSrcDS, GDALRasterBandH hSrcBand, GDALDataType eDstDataType, int bDstHasNoData, double dfDstNoDataValue, typename GDALGeneric3x3ProcessingAlg<T>::type pfnAlg, void* pAlgData, bool bComputeAtEdges ); ~GDALGeneric3x3Dataset(); bool InitOK() const { return apafSourceBuf[0] != nullptr && apafSourceBuf[1] != nullptr && apafSourceBuf[2] != nullptr; } CPLErr GetGeoTransform( double * padfGeoTransform ) override; const OGRSpatialReference* GetSpatialRef() const override; }; /************************************************************************/ /* ==================================================================== */ /* GDALGeneric3x3RasterBand */ /* ==================================================================== */ /************************************************************************/ template<class T> class GDALGeneric3x3RasterBand : public GDALRasterBand { friend class GDALGeneric3x3Dataset<T>; int bSrcHasNoData; T fSrcNoDataValue; int bIsSrcNoDataNan; GDALDataType eReadDT; void InitWithNoData( void* pImage ); public: GDALGeneric3x3RasterBand( GDALGeneric3x3Dataset<T> *poDS, GDALDataType eDstDataType ); virtual CPLErr IReadBlock( int, int, void * ) override; virtual double GetNoDataValue( int* pbHasNoData ) override; }; template<class T> GDALGeneric3x3Dataset<T>::GDALGeneric3x3Dataset( GDALDatasetH hSrcDSIn, GDALRasterBandH hSrcBandIn, GDALDataType eDstDataType, int bDstHasNoDataIn, double dfDstNoDataValueIn, typename GDALGeneric3x3ProcessingAlg<T>::type pfnAlgIn, void* pAlgDataIn, bool bComputeAtEdgesIn ) : pfnAlg(pfnAlgIn), pAlgData(pAlgDataIn), hSrcDS(hSrcDSIn), hSrcBand(hSrcBandIn), bDstHasNoData(bDstHasNoDataIn), dfDstNoDataValue(dfDstNoDataValueIn), nCurLine(-1), bComputeAtEdges(bComputeAtEdgesIn) { CPLAssert(eDstDataType == GDT_Byte || eDstDataType == GDT_Float32); nRasterXSize = GDALGetRasterXSize(hSrcDS); nRasterYSize = GDALGetRasterYSize(hSrcDS); SetBand(1, new GDALGeneric3x3RasterBand<T>(this, eDstDataType)); apafSourceBuf[0] = static_cast<T *>(VSI_MALLOC2_VERBOSE(sizeof(T), nRasterXSize)); apafSourceBuf[1] = static_cast<T *>(VSI_MALLOC2_VERBOSE(sizeof(T), nRasterXSize)); apafSourceBuf[2] = static_cast<T *>(VSI_MALLOC2_VERBOSE(sizeof(T), nRasterXSize)); } template<class T> GDALGeneric3x3Dataset<T>::~GDALGeneric3x3Dataset() { CPLFree(apafSourceBuf[0]); CPLFree(apafSourceBuf[1]); CPLFree(apafSourceBuf[2]); } template<class T> CPLErr GDALGeneric3x3Dataset<T>::GetGeoTransform( double * padfGeoTransform ) { return GDALGetGeoTransform(hSrcDS, padfGeoTransform); } template<class T> const OGRSpatialReference* GDALGeneric3x3Dataset<T>::GetSpatialRef() const { return GDALDataset::FromHandle(hSrcDS)->GetSpatialRef(); } template<class T> GDALGeneric3x3RasterBand<T>::GDALGeneric3x3RasterBand( GDALGeneric3x3Dataset<T> *poDSIn, GDALDataType eDstDataType ) : bSrcHasNoData(FALSE), fSrcNoDataValue(0), bIsSrcNoDataNan(FALSE), eReadDT(GDT_Unknown) { poDS = poDSIn; nBand = 1; eDataType = eDstDataType; nBlockXSize = poDS->GetRasterXSize(); nBlockYSize = 1; const double dfNoDataValue = GDALGetRasterNoDataValue(poDSIn->hSrcBand, &bSrcHasNoData); if( std::numeric_limits<T>::is_integer ) { eReadDT = GDT_Int32; if( bSrcHasNoData ) { GDALDataType eSrcDT = GDALGetRasterDataType(poDSIn->hSrcBand); CPLAssert( eSrcDT == GDT_Byte || eSrcDT == GDT_UInt16 || eSrcDT == GDT_Int16 ); const int nMinVal = (eSrcDT == GDT_Byte ) ? 0 : (eSrcDT == GDT_UInt16) ? 0 : -32768; const int nMaxVal = (eSrcDT == GDT_Byte ) ? 255 : (eSrcDT == GDT_UInt16) ? 65535 : 32767; if( fabs(dfNoDataValue - floor(dfNoDataValue + 0.5)) < 1e-2 && dfNoDataValue >= nMinVal && dfNoDataValue <= nMaxVal ) { fSrcNoDataValue = static_cast<T>(floor(dfNoDataValue + 0.5)); } else { bSrcHasNoData = FALSE; } } } else { eReadDT = GDT_Float32; fSrcNoDataValue = static_cast<T>(dfNoDataValue); bIsSrcNoDataNan = bSrcHasNoData && CPLIsNan(dfNoDataValue); } } template<class T> void GDALGeneric3x3RasterBand<T>::InitWithNoData(void* pImage) { auto poGDS = cpl::down_cast<GDALGeneric3x3Dataset<T> *>(poDS); if( eDataType == GDT_Byte ) { for( int j = 0; j < nBlockXSize; j++ ) static_cast<GByte*>(pImage)[j] = static_cast<GByte>(poGDS->dfDstNoDataValue); } else { for( int j = 0; j < nBlockXSize; j++ ) static_cast<float*>(pImage)[j] = static_cast<float>(poGDS->dfDstNoDataValue); } } template<class T> CPLErr GDALGeneric3x3RasterBand<T>::IReadBlock( int /*nBlockXOff*/, int nBlockYOff, void *pImage ) { auto poGDS = cpl::down_cast<GDALGeneric3x3Dataset<T> *>(poDS); if( poGDS->bComputeAtEdges && nRasterXSize >= 2 && nRasterYSize >= 2 ) { if( nBlockYOff == 0 ) { for( int i = 0; i < 2; i++ ) { CPLErr eErr = GDALRasterIO( poGDS->hSrcBand, GF_Read, 0, i, nBlockXSize, 1, poGDS->apafSourceBuf[i+1], nBlockXSize, 1, eReadDT, 0, 0); if( eErr != CE_None ) { InitWithNoData(pImage); return eErr; } } poGDS->nCurLine = 0; for( int j = 0; j < nRasterXSize; j++ ) { int jmin = (j == 0) ? j : j - 1; int jmax = (j == nRasterXSize - 1) ? j : j + 1; T afWin[9] = { INTERPOL(poGDS->apafSourceBuf[1][jmin], poGDS->apafSourceBuf[2][jmin], bSrcHasNoData, fSrcNoDataValue), INTERPOL(poGDS->apafSourceBuf[1][j], poGDS->apafSourceBuf[2][j], bSrcHasNoData, fSrcNoDataValue), INTERPOL(poGDS->apafSourceBuf[1][jmax], poGDS->apafSourceBuf[2][jmax], bSrcHasNoData, fSrcNoDataValue), poGDS->apafSourceBuf[1][jmin], poGDS->apafSourceBuf[1][j], poGDS->apafSourceBuf[1][jmax], poGDS->apafSourceBuf[2][jmin], poGDS->apafSourceBuf[2][j], poGDS->apafSourceBuf[2][jmax] }; const float fVal = ComputeVal( CPL_TO_BOOL(bSrcHasNoData), fSrcNoDataValue, CPL_TO_BOOL(bIsSrcNoDataNan), afWin, static_cast<float>(poGDS->dfDstNoDataValue), poGDS->pfnAlg, poGDS->pAlgData, poGDS->bComputeAtEdges); if( eDataType == GDT_Byte ) static_cast<GByte*>(pImage)[j] = static_cast<GByte>(fVal + 0.5); else static_cast<float*>(pImage)[j] = fVal; } return CE_None; } else if( nBlockYOff == nRasterYSize - 1 ) { if( poGDS->nCurLine != nRasterYSize - 2 ) { for( int i = 0; i < 2; i++ ) { CPLErr eErr = GDALRasterIO( poGDS->hSrcBand, GF_Read, 0, nRasterYSize - 2 + i, nBlockXSize, 1, poGDS->apafSourceBuf[i+1], nBlockXSize, 1, eReadDT, 0, 0); if( eErr != CE_None ) { InitWithNoData(pImage); return eErr; } } } for( int j = 0; j < nRasterXSize; j++ ) { int jmin = (j == 0) ? j : j - 1; int jmax = (j == nRasterXSize - 1) ? j : j + 1; T afWin[9] = { poGDS->apafSourceBuf[1][jmin], poGDS->apafSourceBuf[1][j], poGDS->apafSourceBuf[1][jmax], poGDS->apafSourceBuf[2][jmin], poGDS->apafSourceBuf[2][j], poGDS->apafSourceBuf[2][jmax], INTERPOL(poGDS->apafSourceBuf[2][jmin], poGDS->apafSourceBuf[1][jmin], bSrcHasNoData, fSrcNoDataValue), INTERPOL(poGDS->apafSourceBuf[2][j], poGDS->apafSourceBuf[1][j], bSrcHasNoData, fSrcNoDataValue), INTERPOL(poGDS->apafSourceBuf[2][jmax], poGDS->apafSourceBuf[1][jmax], bSrcHasNoData, fSrcNoDataValue) }; const float fVal = ComputeVal( CPL_TO_BOOL(bSrcHasNoData), fSrcNoDataValue, CPL_TO_BOOL(bIsSrcNoDataNan), afWin, static_cast<float>(poGDS->dfDstNoDataValue), poGDS->pfnAlg, poGDS->pAlgData, poGDS->bComputeAtEdges); if( eDataType == GDT_Byte ) static_cast<GByte*>(pImage)[j] = static_cast<GByte>(fVal + 0.5); else static_cast<float*>(pImage)[j] = fVal; } return CE_None; } } else if( nBlockYOff == 0 || nBlockYOff == nRasterYSize - 1 ) { InitWithNoData(pImage); return CE_None; } if( poGDS->nCurLine != nBlockYOff ) { if( poGDS->nCurLine + 1 == nBlockYOff ) { T* pafTmp = poGDS->apafSourceBuf[0]; poGDS->apafSourceBuf[0] = poGDS->apafSourceBuf[1]; poGDS->apafSourceBuf[1] = poGDS->apafSourceBuf[2]; poGDS->apafSourceBuf[2] = pafTmp; CPLErr eErr = GDALRasterIO( poGDS->hSrcBand, GF_Read, 0, nBlockYOff + 1, nBlockXSize, 1, poGDS->apafSourceBuf[2], nBlockXSize, 1, eReadDT, 0, 0); if( eErr != CE_None ) { InitWithNoData(pImage); return eErr; } } else { for( int i = 0; i < 3; i++ ) { const CPLErr eErr = GDALRasterIO(poGDS->hSrcBand, GF_Read, 0, nBlockYOff + i - 1, nBlockXSize, 1, poGDS->apafSourceBuf[i], nBlockXSize, 1, eReadDT, 0, 0); if( eErr != CE_None ) { InitWithNoData(pImage); return eErr; } } } poGDS->nCurLine = nBlockYOff; } if( poGDS->bComputeAtEdges && nRasterXSize >= 2 ) { int j = 0; T afWin[9] = { INTERPOL(poGDS->apafSourceBuf[0][j], poGDS->apafSourceBuf[0][j+1], bSrcHasNoData, fSrcNoDataValue), poGDS->apafSourceBuf[0][j], poGDS->apafSourceBuf[0][j+1], INTERPOL(poGDS->apafSourceBuf[1][j], poGDS->apafSourceBuf[1][j+1], bSrcHasNoData, fSrcNoDataValue), poGDS->apafSourceBuf[1][j], poGDS->apafSourceBuf[1][j+1], INTERPOL(poGDS->apafSourceBuf[2][j], poGDS->apafSourceBuf[2][j+1], bSrcHasNoData, fSrcNoDataValue), poGDS->apafSourceBuf[2][j], poGDS->apafSourceBuf[2][j+1] }; { const float fVal = ComputeVal( CPL_TO_BOOL(bSrcHasNoData), fSrcNoDataValue, CPL_TO_BOOL(bIsSrcNoDataNan), afWin, static_cast<float>(poGDS->dfDstNoDataValue), poGDS->pfnAlg, poGDS->pAlgData, poGDS->bComputeAtEdges); if( eDataType == GDT_Byte ) static_cast<GByte*>(pImage)[j] = static_cast<GByte>(fVal + 0.5); else static_cast<float*>(pImage)[j] = fVal; } j = nRasterXSize - 1; afWin[0] = poGDS->apafSourceBuf[0][j-1]; afWin[1] = poGDS->apafSourceBuf[0][j]; afWin[2] = INTERPOL(poGDS->apafSourceBuf[0][j], poGDS->apafSourceBuf[0][j-1], bSrcHasNoData, fSrcNoDataValue); afWin[3] = poGDS->apafSourceBuf[1][j-1]; afWin[4] = poGDS->apafSourceBuf[1][j]; afWin[5] = INTERPOL(poGDS->apafSourceBuf[1][j], poGDS->apafSourceBuf[1][j-1], bSrcHasNoData, fSrcNoDataValue); afWin[6] = poGDS->apafSourceBuf[2][j-1]; afWin[7] = poGDS->apafSourceBuf[2][j]; afWin[8] = INTERPOL(poGDS->apafSourceBuf[2][j], poGDS->apafSourceBuf[2][j-1], bSrcHasNoData, fSrcNoDataValue); const float fVal = ComputeVal( CPL_TO_BOOL(bSrcHasNoData), fSrcNoDataValue, CPL_TO_BOOL(bIsSrcNoDataNan), afWin, static_cast<float>(poGDS->dfDstNoDataValue), poGDS->pfnAlg, poGDS->pAlgData, poGDS->bComputeAtEdges); if( eDataType == GDT_Byte ) static_cast<GByte*>(pImage)[j] = static_cast<GByte>(fVal + 0.5); else static_cast<float*>(pImage)[j] = fVal; } else { if( eDataType == GDT_Byte ) { static_cast<GByte*>(pImage)[0] = static_cast<GByte>(poGDS->dfDstNoDataValue); if( nBlockXSize > 1 ) static_cast<GByte*>(pImage)[nBlockXSize - 1] = static_cast<GByte>(poGDS->dfDstNoDataValue); } else { static_cast<float*>(pImage)[0] = static_cast<float>(poGDS->dfDstNoDataValue); if( nBlockXSize > 1 ) static_cast<float*>(pImage)[nBlockXSize - 1] = static_cast<float>(poGDS->dfDstNoDataValue); } } for( int j = 1; j < nBlockXSize - 1; j++ ) { T afWin[9] = { poGDS->apafSourceBuf[0][j-1], poGDS->apafSourceBuf[0][j], poGDS->apafSourceBuf[0][j+1], poGDS->apafSourceBuf[1][j-1], poGDS->apafSourceBuf[1][j], poGDS->apafSourceBuf[1][j+1], poGDS->apafSourceBuf[2][j-1], poGDS->apafSourceBuf[2][j], poGDS->apafSourceBuf[2][j+1] }; const float fVal = ComputeVal( CPL_TO_BOOL(bSrcHasNoData), fSrcNoDataValue, CPL_TO_BOOL(bIsSrcNoDataNan), afWin, static_cast<float>(poGDS->dfDstNoDataValue), poGDS->pfnAlg, poGDS->pAlgData, poGDS->bComputeAtEdges); if( eDataType == GDT_Byte ) static_cast<GByte*>(pImage)[j] = static_cast<GByte>(fVal + 0.5); else static_cast<float*>(pImage)[j] = fVal; } return CE_None; } template<class T> double GDALGeneric3x3RasterBand<T>::GetNoDataValue( int* pbHasNoData ) { auto poGDS = cpl::down_cast<GDALGeneric3x3Dataset<T> *>(poDS); if( pbHasNoData ) *pbHasNoData = poGDS->bDstHasNoData; return poGDS->dfDstNoDataValue; } /************************************************************************/ /* ArgIsNumeric() */ /************************************************************************/ static int ArgIsNumeric( const char *pszArg ) { return CPLGetValueType(pszArg) != CPL_VALUE_STRING; } /************************************************************************/ /* GetAlgorithm() */ /************************************************************************/ typedef enum { INVALID, HILL_SHADE, SLOPE, ASPECT, COLOR_RELIEF, TRI, TPI, ROUGHNESS } Algorithm; static Algorithm GetAlgorithm(const char* pszProcessing) { if( EQUAL(pszProcessing, "shade") || EQUAL(pszProcessing, "hillshade") ) { return HILL_SHADE; } else if( EQUAL(pszProcessing, "slope") ) { return SLOPE; } else if( EQUAL(pszProcessing, "aspect") ) { return ASPECT; } else if( EQUAL(pszProcessing, "color-relief") ) { return COLOR_RELIEF; } else if( EQUAL(pszProcessing, "TRI") ) { return TRI; } else if( EQUAL(pszProcessing, "TPI") ) { return TPI; } else if( EQUAL(pszProcessing, "roughness") ) { return ROUGHNESS; } else { return INVALID; } } /************************************************************************/ /* GDALDEMProcessing() */ /************************************************************************/ /** * Apply a DEM processing. * * This is the equivalent of the <a href="gdaldem.html">gdaldem</a> utility. * * GDALDEMProcessingOptions* must be allocated and freed with * GDALDEMProcessingOptionsNew() and GDALDEMProcessingOptionsFree() * respectively. * * @param pszDest the destination dataset path. * @param hSrcDataset the source dataset handle. * @param pszProcessing the processing to apply (one of "hillshade", "slope", * "aspect", "color-relief", "TRI", "TPI", "Roughness") * @param pszColorFilename color file (mandatory for "color-relief" processing, * should be NULL otherwise) * @param psOptionsIn the options struct returned by * GDALDEMProcessingOptionsNew() or NULL. * @param pbUsageError the pointer to int variable to determine any usage * error has occurred or NULL. * @return the output dataset (new dataset that must be closed using * GDALClose()) or NULL in case of error. * * @since GDAL 2.1 */ GDALDatasetH GDALDEMProcessing( const char *pszDest, GDALDatasetH hSrcDataset, const char* pszProcessing, const char* pszColorFilename, const GDALDEMProcessingOptions *psOptionsIn, int *pbUsageError ) { if( hSrcDataset == nullptr ) { CPLError( CE_Failure, CPLE_AppDefined, "No source dataset specified."); if(pbUsageError) *pbUsageError = TRUE; return nullptr; } if( pszDest == nullptr ) { CPLError( CE_Failure, CPLE_AppDefined, "No target dataset specified."); if(pbUsageError) *pbUsageError = TRUE; return nullptr; } if( pszProcessing == nullptr ) { CPLError( CE_Failure, CPLE_AppDefined, "No target dataset specified."); if(pbUsageError) *pbUsageError = TRUE; return nullptr; } Algorithm eUtilityMode = GetAlgorithm(pszProcessing); if( eUtilityMode == INVALID ) { CPLError(CE_Failure, CPLE_IllegalArg, "Invalid processing"); if(pbUsageError) *pbUsageError = TRUE; return nullptr; } if( eUtilityMode == COLOR_RELIEF && pszColorFilename == nullptr ) { CPLError( CE_Failure, CPLE_AppDefined, "pszColorFilename == NULL."); if(pbUsageError) *pbUsageError = TRUE; return nullptr; } else if( eUtilityMode != COLOR_RELIEF && pszColorFilename != nullptr ) { CPLError( CE_Failure, CPLE_AppDefined, "pszColorFilename != NULL."); if(pbUsageError) *pbUsageError = TRUE; return nullptr; } if( psOptionsIn && psOptionsIn->bCombined && eUtilityMode != HILL_SHADE ) { CPLError(CE_Failure, CPLE_NotSupported, "-combined can only be used with hillshade"); if(pbUsageError) *pbUsageError = TRUE; return nullptr; } if( psOptionsIn && psOptionsIn->bIgor && eUtilityMode != HILL_SHADE ) { CPLError(CE_Failure, CPLE_NotSupported, "-igor can only be used with hillshade"); if(pbUsageError) *pbUsageError = TRUE; return nullptr; } if( psOptionsIn && psOptionsIn->bMultiDirectional && eUtilityMode != HILL_SHADE ) { CPLError(CE_Failure, CPLE_NotSupported, "-multidirectional can only be used with hillshade"); if(pbUsageError) *pbUsageError = TRUE; return nullptr; } GDALDEMProcessingOptions* psOptionsToFree = nullptr; const GDALDEMProcessingOptions* psOptions = psOptionsIn; if( !psOptions ) { psOptionsToFree = GDALDEMProcessingOptionsNew(nullptr, nullptr); psOptions = psOptionsToFree; } double adfGeoTransform[6] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 }; const int nXSize = GDALGetRasterXSize(hSrcDataset); const int nYSize = GDALGetRasterYSize(hSrcDataset); if( psOptions->nBand <= 0 || psOptions->nBand > GDALGetRasterCount(hSrcDataset) ) { CPLError(CE_Failure, CPLE_IllegalArg, "Unable to fetch band #%d", psOptions->nBand ); GDALDEMProcessingOptionsFree(psOptionsToFree); return nullptr; } GDALRasterBandH hSrcBand = GDALGetRasterBand( hSrcDataset, psOptions->nBand ); GDALGetGeoTransform(hSrcDataset, adfGeoTransform); CPLString osFormat; if( psOptions->pszFormat == nullptr ) { osFormat = GetOutputDriverForRaster(pszDest); if( osFormat.empty() ) { GDALDEMProcessingOptionsFree(psOptionsToFree); return nullptr; } } else { osFormat = psOptions->pszFormat; } GDALDriverH hDriver = GDALGetDriverByName(osFormat); if( hDriver == nullptr || (GDALGetMetadataItem( hDriver, GDAL_DCAP_CREATE, nullptr ) == nullptr && GDALGetMetadataItem( hDriver, GDAL_DCAP_CREATECOPY, nullptr ) == nullptr)) { CPLError(CE_Failure, CPLE_AppDefined, "Output driver `%s' not recognised to have output support.", osFormat.c_str() ); fprintf( stderr, "The following format drivers are configured\n" "and support output:\n" ); for( int iDr = 0; iDr < GDALGetDriverCount(); iDr++ ) { hDriver = GDALGetDriver(iDr); if( GDALGetMetadataItem( hDriver, GDAL_DCAP_RASTER, nullptr) != nullptr && (GDALGetMetadataItem( hDriver, GDAL_DCAP_CREATE, nullptr ) != nullptr || GDALGetMetadataItem( hDriver, GDAL_DCAP_CREATECOPY, nullptr ) != nullptr) ) { fprintf( stderr, " %s: %s\n", GDALGetDriverShortName( hDriver ), GDALGetDriverLongName( hDriver ) ); } } GDALDEMProcessingOptionsFree(psOptionsToFree); return nullptr; } double dfDstNoDataValue = 0.0; bool bDstHasNoData = false; void* pData = nullptr; GDALGeneric3x3ProcessingAlg<float>::type pfnAlgFloat = nullptr; GDALGeneric3x3ProcessingAlg<GInt32>::type pfnAlgInt32 = nullptr; GDALGeneric3x3ProcessingAlg_multisample<GInt32>::type pfnAlgInt32_multisample = nullptr; if( eUtilityMode == HILL_SHADE && psOptions->bMultiDirectional ) { dfDstNoDataValue = 0; bDstHasNoData = true; pData = GDALCreateHillshadeMultiDirectionalData(adfGeoTransform, psOptions->z, psOptions->scale, psOptions->alt, psOptions->bZevenbergenThorne); if( psOptions->bZevenbergenThorne ) { pfnAlgFloat = GDALHillshadeMultiDirectionalAlg<float, ZEVENBERGEN_THORNE>; pfnAlgInt32 = GDALHillshadeMultiDirectionalAlg<GInt32, ZEVENBERGEN_THORNE>; } else { pfnAlgFloat = GDALHillshadeMultiDirectionalAlg<float, HORN>; pfnAlgInt32 = GDALHillshadeMultiDirectionalAlg<GInt32, HORN>; } } else if( eUtilityMode == HILL_SHADE ) { dfDstNoDataValue = 0; bDstHasNoData = true; pData = GDALCreateHillshadeData(adfGeoTransform, psOptions->z, psOptions->scale, psOptions->alt, psOptions->az, psOptions->bZevenbergenThorne); if( psOptions->bZevenbergenThorne ) { if( psOptions->bCombined ) { pfnAlgFloat = GDALHillshadeCombinedAlg<float, ZEVENBERGEN_THORNE>; pfnAlgInt32 = GDALHillshadeCombinedAlg<GInt32, ZEVENBERGEN_THORNE>; } else if( psOptions->bIgor ) { pfnAlgFloat = GDALHillshadeIgorAlg<float, ZEVENBERGEN_THORNE>; pfnAlgInt32 = GDALHillshadeIgorAlg<GInt32, ZEVENBERGEN_THORNE>; } else { pfnAlgFloat = GDALHillshadeAlg<float, ZEVENBERGEN_THORNE>; pfnAlgInt32 = GDALHillshadeAlg<GInt32, ZEVENBERGEN_THORNE>; } } else { if( psOptions->bCombined ) { pfnAlgFloat = GDALHillshadeCombinedAlg<float, HORN>; pfnAlgInt32 = GDALHillshadeCombinedAlg<GInt32, HORN>; } else if( psOptions->bIgor ) { pfnAlgFloat = GDALHillshadeIgorAlg<float, HORN>; pfnAlgInt32 = GDALHillshadeIgorAlg<GInt32, HORN>; } else { if( adfGeoTransform[1] == -adfGeoTransform[5] ) { pfnAlgFloat = GDALHillshadeAlg_same_res<float>; pfnAlgInt32 = GDALHillshadeAlg_same_res<GInt32>; #ifdef HAVE_16_SSE_REG pfnAlgInt32_multisample = GDALHillshadeAlg_same_res_multisample<GInt32>; #endif } else { pfnAlgFloat = GDALHillshadeAlg<float, HORN>; pfnAlgInt32 = GDALHillshadeAlg<GInt32, HORN>; } } } } else if( eUtilityMode == SLOPE ) { dfDstNoDataValue = -9999; bDstHasNoData = true; pData = GDALCreateSlopeData(adfGeoTransform, psOptions->scale, psOptions->slopeFormat); if( psOptions->bZevenbergenThorne ) { pfnAlgFloat = GDALSlopeZevenbergenThorneAlg<float>; pfnAlgInt32 = GDALSlopeZevenbergenThorneAlg<GInt32>; } else { pfnAlgFloat = GDALSlopeHornAlg<float>; pfnAlgInt32 = GDALSlopeHornAlg<GInt32>; } } else if( eUtilityMode == ASPECT ) { if( !psOptions->bZeroForFlat ) { dfDstNoDataValue = -9999; bDstHasNoData = true; } pData = GDALCreateAspectData(psOptions->bAngleAsAzimuth); if( psOptions->bZevenbergenThorne ) { pfnAlgFloat = GDALAspectZevenbergenThorneAlg<float>; pfnAlgInt32 = GDALAspectZevenbergenThorneAlg<GInt32>; } else { pfnAlgFloat = GDALAspectAlg<float>; pfnAlgInt32 = GDALAspectAlg<GInt32>; } } else if( eUtilityMode == TRI ) { dfDstNoDataValue = -9999; bDstHasNoData = true; pfnAlgFloat = GDALTRIAlg<float>; pfnAlgInt32 = GDALTRIAlg<GInt32>; } else if( eUtilityMode == TPI ) { dfDstNoDataValue = -9999; bDstHasNoData = true; pfnAlgFloat = GDALTPIAlg<float>; pfnAlgInt32 = GDALTPIAlg<GInt32>; } else if( eUtilityMode == ROUGHNESS ) { dfDstNoDataValue = -9999; bDstHasNoData = true; pfnAlgFloat = GDALRoughnessAlg<float>; pfnAlgInt32 = GDALRoughnessAlg<GInt32>; } const GDALDataType eDstDataType = (eUtilityMode == HILL_SHADE || eUtilityMode == COLOR_RELIEF) ? GDT_Byte : GDT_Float32; if( EQUAL(osFormat, "VRT") ) { if( eUtilityMode == COLOR_RELIEF ) { GDALGenerateVRTColorRelief(pszDest, hSrcDataset, hSrcBand, pszColorFilename, psOptions->eColorSelectionMode, psOptions->bAddAlpha); CPLFree(pData); GDALDEMProcessingOptionsFree(psOptionsToFree); return GDALOpen(pszDest, GA_Update); } else { CPLError(CE_Failure, CPLE_NotSupported, "VRT driver can only be used with color-relief utility."); GDALDEMProcessingOptionsFree(psOptionsToFree); CPLFree(pData); return nullptr; } } // We might actually want to always go through the intermediate dataset bool bForceUseIntermediateDataset = false; GDALProgressFunc pfnProgress = psOptions->pfnProgress; void* pProgressData = psOptions->pProgressData; if( EQUAL(osFormat, "GTiff") ) { if( !EQUAL(CSLFetchNameValueDef(psOptions->papszCreateOptions, "COMPRESS", "NONE"), "NONE") && CPLTestBool(CSLFetchNameValueDef(psOptions->papszCreateOptions, "TILED", "NO")) ) { bForceUseIntermediateDataset = true; } else if( strcmp(pszDest, "/vsistdout/") == 0 ) { bForceUseIntermediateDataset = true; pfnProgress = GDALDummyProgress; pProgressData = nullptr; } #ifdef S_ISFIFO else { VSIStatBufL sStat; if( VSIStatExL(pszDest, &sStat, VSI_STAT_EXISTS_FLAG | VSI_STAT_NATURE_FLAG) == 0 && S_ISFIFO(sStat.st_mode) ) { bForceUseIntermediateDataset = true; } } #endif } const GDALDataType eSrcDT = GDALGetRasterDataType(hSrcBand); if( GDALGetMetadataItem( hDriver, GDAL_DCAP_RASTER, nullptr) != nullptr && ((bForceUseIntermediateDataset || GDALGetMetadataItem( hDriver, GDAL_DCAP_CREATE, nullptr ) == nullptr) && GDALGetMetadataItem( hDriver, GDAL_DCAP_CREATECOPY, nullptr ) != nullptr) ) { GDALDatasetH hIntermediateDataset = nullptr; if( eUtilityMode == COLOR_RELIEF ) { GDALColorReliefDataset* poDS = new GDALColorReliefDataset(hSrcDataset, hSrcBand, pszColorFilename, psOptions->eColorSelectionMode, psOptions->bAddAlpha); if( !(poDS->InitOK()) ) { delete poDS; CPLFree(pData); GDALDEMProcessingOptionsFree(psOptionsToFree); return nullptr; } hIntermediateDataset = static_cast<GDALDatasetH>(poDS); } else { if( eSrcDT == GDT_Byte || eSrcDT == GDT_Int16 || eSrcDT == GDT_UInt16 ) { GDALGeneric3x3Dataset<GInt32>* poDS = new GDALGeneric3x3Dataset<GInt32>(hSrcDataset, hSrcBand, eDstDataType, bDstHasNoData, dfDstNoDataValue, pfnAlgInt32, pData, psOptions->bComputeAtEdges); if( !(poDS->InitOK()) ) { delete poDS; CPLFree(pData); GDALDEMProcessingOptionsFree(psOptionsToFree); return nullptr; } hIntermediateDataset = static_cast<GDALDatasetH>(poDS); } else { GDALGeneric3x3Dataset<float>* poDS = new GDALGeneric3x3Dataset<float>(hSrcDataset, hSrcBand, eDstDataType, bDstHasNoData, dfDstNoDataValue, pfnAlgFloat, pData, psOptions->bComputeAtEdges); if( !(poDS->InitOK()) ) { delete poDS; CPLFree(pData); GDALDEMProcessingOptionsFree(psOptionsToFree); return nullptr; } hIntermediateDataset = static_cast<GDALDatasetH>(poDS); } } GDALDatasetH hOutDS = GDALCreateCopy( hDriver, pszDest, hIntermediateDataset, TRUE, psOptions->papszCreateOptions, pfnProgress, pProgressData ); GDALClose(hIntermediateDataset); CPLFree(pData); GDALDEMProcessingOptionsFree(psOptionsToFree); return hOutDS; } const int nDstBands = eUtilityMode == COLOR_RELIEF ? ((psOptions->bAddAlpha) ? 4 : 3) : 1; GDALDatasetH hDstDataset = GDALCreate( hDriver, pszDest, nXSize, nYSize, nDstBands, eDstDataType, psOptions->papszCreateOptions ); if( hDstDataset == nullptr ) { CPLError(CE_Failure, CPLE_AppDefined, "Unable to create dataset %s", pszDest ); GDALDEMProcessingOptionsFree(psOptionsToFree); CPLFree(pData); return nullptr; } GDALRasterBandH hDstBand = GDALGetRasterBand( hDstDataset, 1 ); GDALSetGeoTransform(hDstDataset, adfGeoTransform); GDALSetProjection(hDstDataset, GDALGetProjectionRef(hSrcDataset)); if( eUtilityMode == COLOR_RELIEF ) { GDALColorRelief (hSrcBand, GDALGetRasterBand(hDstDataset, 1), GDALGetRasterBand(hDstDataset, 2), GDALGetRasterBand(hDstDataset, 3), psOptions->bAddAlpha ? GDALGetRasterBand(hDstDataset, 4) : nullptr, pszColorFilename, psOptions->eColorSelectionMode, pfnProgress, pProgressData); } else { if( bDstHasNoData ) GDALSetRasterNoDataValue(hDstBand, dfDstNoDataValue); if( eSrcDT == GDT_Byte || eSrcDT == GDT_Int16 || eSrcDT == GDT_UInt16 ) { GDALGeneric3x3Processing<GInt32>(hSrcBand, hDstBand, pfnAlgInt32, pfnAlgInt32_multisample, pData, psOptions->bComputeAtEdges, pfnProgress, pProgressData); } else { GDALGeneric3x3Processing<float>(hSrcBand, hDstBand, pfnAlgFloat, nullptr, pData, psOptions->bComputeAtEdges, pfnProgress, pProgressData); } } CPLFree(pData); GDALDEMProcessingOptionsFree(psOptionsToFree); return hDstDataset; } /************************************************************************/ /* GDALDEMProcessingOptionsNew() */ /************************************************************************/ /** * Allocates a GDALDEMProcessingOptions struct. * * @param papszArgv NULL terminated list of options (potentially including filename and open options too), or NULL. * The accepted options are the ones of the <a href="gdaldem.html">gdaldem</a> utility. * @param psOptionsForBinary (output) may be NULL (and should generally be NULL), * otherwise (gdal_translate_bin.cpp use case) must be allocated with * GDALDEMProcessingOptionsForBinaryNew() prior to this function. Will be * filled with potentially present filename, open options,... * @return pointer to the allocated GDALDEMProcessingOptions struct. Must be freed with GDALDEMProcessingOptionsFree(). * * @since GDAL 2.1 */ GDALDEMProcessingOptions *GDALDEMProcessingOptionsNew( char** papszArgv, GDALDEMProcessingOptionsForBinary* psOptionsForBinary ) { GDALDEMProcessingOptions *psOptions = static_cast<GDALDEMProcessingOptions *>( CPLCalloc(1, sizeof(GDALDEMProcessingOptions))); Algorithm eUtilityMode = INVALID; psOptions->pszFormat = nullptr; psOptions->pfnProgress = GDALDummyProgress; psOptions->pProgressData = nullptr; psOptions->z = 1.0; psOptions->scale = 1.0; psOptions->az = 315.0; psOptions->alt = 45.0; // 0 = 'percent' or 1 = 'degrees' psOptions->slopeFormat = 1; psOptions->bAddAlpha = false; psOptions->bZeroForFlat = false; psOptions->bAngleAsAzimuth = true; psOptions->eColorSelectionMode = COLOR_SELECTION_INTERPOLATE; psOptions->bComputeAtEdges = false; psOptions->bZevenbergenThorne = false; psOptions->bCombined = false; psOptions->bIgor = false; psOptions->bMultiDirectional = false; psOptions->nBand = 1; psOptions->papszCreateOptions = nullptr; bool bAzimuthSpecified = false; bool bAltSpecified = false; /* -------------------------------------------------------------------- */ /* Handle command line arguments. */ /* -------------------------------------------------------------------- */ int argc = CSLCount(papszArgv); for( int i = 0; papszArgv != nullptr && i < argc; i++ ) { if( i == 0 && psOptionsForBinary ) { eUtilityMode = GetAlgorithm(papszArgv[0]); if(eUtilityMode == INVALID ) { CPLError(CE_Failure, CPLE_IllegalArg, "Invalid utility mode"); GDALDEMProcessingOptionsFree(psOptions); return nullptr; } psOptionsForBinary->pszProcessing = CPLStrdup(papszArgv[0]); continue; } if( i < argc-1 && (EQUAL(papszArgv[i],"-of") || EQUAL(papszArgv[i],"-f")) ) { ++i; CPLFree(psOptions->pszFormat); psOptions->pszFormat = CPLStrdup(papszArgv[i]); } else if( EQUAL(papszArgv[i],"-q") || EQUAL(papszArgv[i],"-quiet") ) { if( psOptionsForBinary ) psOptionsForBinary->bQuiet = TRUE; } else if( (EQUAL(papszArgv[i], "--z") || EQUAL(papszArgv[i], "-z")) && i + 1 < argc ) { ++i; if( !ArgIsNumeric(papszArgv[i]) ) { CPLError(CE_Failure, CPLE_IllegalArg, "Numeric value expected for -z"); GDALDEMProcessingOptionsFree(psOptions); return nullptr; } psOptions->z = CPLAtof(papszArgv[i]); } else if( EQUAL(papszArgv[i], "-p") ) { psOptions->slopeFormat = 0; } else if( EQUAL(papszArgv[i], "-alg") && i+1<argc ) { i++; if( EQUAL(papszArgv[i], "ZevenbergenThorne") ) { psOptions->bZevenbergenThorne = true; } else if( !EQUAL(papszArgv[i], "Horn") ) { CPLError(CE_Failure, CPLE_IllegalArg, "Numeric value expected for %s", papszArgv[i-1]); GDALDEMProcessingOptionsFree(psOptions); return nullptr; } } else if( EQUAL(papszArgv[i], "-trigonometric") ) { psOptions->bAngleAsAzimuth = false; } else if( EQUAL(papszArgv[i], "-zero_for_flat") ) { psOptions->bZeroForFlat = true; } else if( EQUAL(papszArgv[i], "-exact_color_entry") ) { psOptions->eColorSelectionMode = COLOR_SELECTION_EXACT_ENTRY; } else if( EQUAL(papszArgv[i], "-nearest_color_entry") ) { psOptions->eColorSelectionMode = COLOR_SELECTION_NEAREST_ENTRY; } else if( (EQUAL(papszArgv[i], "--s") || EQUAL(papszArgv[i], "-s") || EQUAL(papszArgv[i], "--scale") || EQUAL(papszArgv[i], "-scale")) && i+1<argc ) { ++i; if( !ArgIsNumeric(papszArgv[i]) ) { CPLError(CE_Failure, CPLE_IllegalArg, "Numeric value expected for %s", papszArgv[i-1]); GDALDEMProcessingOptionsFree(psOptions); return nullptr; } psOptions->scale = CPLAtof(papszArgv[i]); } else if( (EQUAL(papszArgv[i], "--az") || EQUAL(papszArgv[i], "-az") || EQUAL(papszArgv[i], "--azimuth") || EQUAL(papszArgv[i], "-azimuth")) && i+1<argc ) { ++i; if( !ArgIsNumeric(papszArgv[i]) ) { CPLError(CE_Failure, CPLE_IllegalArg, "Numeric value expected for %s", papszArgv[i-1]); GDALDEMProcessingOptionsFree(psOptions); return nullptr; } bAzimuthSpecified = true; psOptions->az = CPLAtof(papszArgv[i]); } else if( (EQUAL(papszArgv[i], "--alt") || EQUAL(papszArgv[i], "-alt") || EQUAL(papszArgv[i], "--altitude") || EQUAL(papszArgv[i], "-altitude")) && i+1<argc ) { ++i; if( !ArgIsNumeric(papszArgv[i]) ) { CPLError(CE_Failure, CPLE_IllegalArg, "Numeric value expected for %s", papszArgv[i-1]); GDALDEMProcessingOptionsFree(psOptions); return nullptr; } bAltSpecified = true; psOptions->alt = CPLAtof(papszArgv[i]); } else if( (EQUAL(papszArgv[i], "-combined") || EQUAL(papszArgv[i], "--combined")) ) { psOptions->bCombined = true; } else if( (EQUAL(papszArgv[i], "-igor") || EQUAL(papszArgv[i], "--igor")) ) { psOptions->bIgor = true; } else if( EQUAL(papszArgv[i], "-multidirectional") || EQUAL(papszArgv[i], "--multidirectional") ) { psOptions->bMultiDirectional = true; } else if( EQUAL(papszArgv[i], "-alpha")) { psOptions->bAddAlpha = true; } else if( EQUAL(papszArgv[i], "-compute_edges")) { psOptions->bComputeAtEdges = true; } else if( i + 1 < argc && (EQUAL(papszArgv[i], "--b") || EQUAL(papszArgv[i], "-b")) ) { psOptions->nBand = atoi(papszArgv[++i]); } else if( EQUAL(papszArgv[i],"-co") && i+1<argc ) { psOptions->papszCreateOptions = CSLAddString( psOptions->papszCreateOptions, papszArgv[++i] ); } else if( papszArgv[i][0] == '-' ) { CPLError(CE_Failure, CPLE_NotSupported, "Unknown option name '%s'", papszArgv[i]); GDALDEMProcessingOptionsFree(psOptions); return nullptr; } else if( psOptionsForBinary && psOptionsForBinary->pszSrcFilename == nullptr ) { psOptionsForBinary->pszSrcFilename = CPLStrdup(papszArgv[i]); } else if( psOptionsForBinary && eUtilityMode == COLOR_RELIEF && psOptionsForBinary->pszColorFilename == nullptr ) { psOptionsForBinary->pszColorFilename = CPLStrdup(papszArgv[i]); } else if( psOptionsForBinary && psOptionsForBinary->pszDstFilename == nullptr ) { psOptionsForBinary->pszDstFilename = CPLStrdup(papszArgv[i]); } else { CPLError(CE_Failure, CPLE_NotSupported, "Too many command options '%s'", papszArgv[i]); GDALDEMProcessingOptionsFree(psOptions); return nullptr; } } if( psOptions->bMultiDirectional + psOptions->bCombined + psOptions->bIgor > 1) { CPLError(CE_Failure, CPLE_NotSupported, "only one of -multidirectional, -combined or -igor can be used"); GDALDEMProcessingOptionsFree(psOptions); return nullptr; } if( psOptions->bMultiDirectional && bAzimuthSpecified ) { CPLError(CE_Failure, CPLE_NotSupported, "-multidirectional and -az cannot be used together"); GDALDEMProcessingOptionsFree(psOptions); return nullptr; } if( psOptions->bIgor && bAltSpecified ) { CPLError(CE_Failure, CPLE_NotSupported, "-igor and -alt cannot be used together"); GDALDEMProcessingOptionsFree(psOptions); return nullptr; } return psOptions; } /************************************************************************/ /* GDALDEMProcessingOptionsFree() */ /************************************************************************/ /** * Frees the GDALDEMProcessingOptions struct. * * @param psOptions the options struct for GDALDEMProcessing(). * * @since GDAL 2.1 */ void GDALDEMProcessingOptionsFree( GDALDEMProcessingOptions *psOptions ) { if( psOptions ) { CPLFree(psOptions->pszFormat); CSLDestroy(psOptions->papszCreateOptions); CPLFree(psOptions); } } /************************************************************************/ /* GDALDEMProcessingOptionsSetProgress() */ /************************************************************************/ /** * Set a progress function. * * @param psOptions the options struct for GDALDEMProcessing(). * @param pfnProgress the progress callback. * @param pProgressData the user data for the progress callback. * * @since GDAL 2.1 */ void GDALDEMProcessingOptionsSetProgress( GDALDEMProcessingOptions *psOptions, GDALProgressFunc pfnProgress, void *pProgressData ) { psOptions->pfnProgress = pfnProgress; psOptions->pProgressData = pProgressData; }