mirror of https://gitee.com/openkylin/gdal.git
4059 lines
143 KiB
C++
4059 lines
143 KiB
C++
/******************************************************************************
|
|
*
|
|
* Project: GDAL DEM Utilities
|
|
* Purpose:
|
|
* Authors: Matthew Perry, perrygeo at gmail.com
|
|
* Even Rouault, even dot rouault at spatialys.com
|
|
* 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 spatialys.com>
|
|
* 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 -
|
|
* For bathymetric use cases, implements
|
|
* Terrain Ruggedness Index is as described in Wilson et al. (2007)
|
|
* this is based on the method of Valentine et al. (2004)
|
|
*
|
|
* For terrestrial use cases, implements
|
|
* Riley, S.J., De Gloria, S.D., Elliot, R. (1999): A Terrain Ruggedness
|
|
* that Quantifies Topographic Heterogeneity. Intermountain Journal of Science,
|
|
*Vol.5, No.1-4, pp.23-27
|
|
*
|
|
*
|
|
* 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
|
|
|
|
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;
|
|
|
|
namespace
|
|
{
|
|
enum class GradientAlg
|
|
{
|
|
HORN,
|
|
ZEVENBERGEN_THORNE,
|
|
};
|
|
|
|
enum class TRIAlg
|
|
{
|
|
WILSON,
|
|
RILEY,
|
|
};
|
|
} // namespace
|
|
|
|
struct GDALDEMProcessingOptions
|
|
{
|
|
/*! output format. Use the short format name. */
|
|
char *pszFormat = nullptr;
|
|
|
|
/*! the progress function to use */
|
|
GDALProgressFunc pfnProgress = nullptr;
|
|
|
|
/*! pointer to the progress data variable */
|
|
void *pProgressData = nullptr;
|
|
|
|
double z = 1.0;
|
|
double scale = 1.0;
|
|
double az = 315.0;
|
|
double alt = 45.0;
|
|
int slopeFormat = 1; // 0 = 'percent' or 1 = 'degrees'
|
|
bool bAddAlpha = false;
|
|
bool bZeroForFlat = false;
|
|
bool bAngleAsAzimuth = true;
|
|
ColorSelectionMode eColorSelectionMode = COLOR_SELECTION_INTERPOLATE;
|
|
bool bComputeAtEdges = false;
|
|
bool bGradientAlgSpecified = false;
|
|
GradientAlg eGradientAlg = GradientAlg::HORN;
|
|
bool bTRIAlgSpecified = false;
|
|
TRIAlg eTRIAlg = TRIAlg::RILEY;
|
|
bool bCombined = false;
|
|
bool bIgor = false;
|
|
bool bMultiDirectional = false;
|
|
char **papszCreateOptions = nullptr;
|
|
int nBand = 1;
|
|
};
|
|
|
|
/************************************************************************/
|
|
/* 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 */
|
|
/************************************************************************/
|
|
|
|
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, GradientAlg::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, GradientAlg::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 == GradientAlg::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 == GradientAlg::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,
|
|
GradientAlg eAlg)
|
|
{
|
|
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 / ((eAlg == GradientAlg::ZEVENBERGEN_THORNE ? 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,
|
|
GradientAlg eAlg)
|
|
{
|
|
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 / ((eAlg == GradientAlg::ZEVENBERGEN_THORNE ? 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];
|
|
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];
|
|
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 behavior: 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"))
|
|
{
|
|
if (!bSrcHasNoData)
|
|
{
|
|
CPLError(CE_Warning, CPLE_AppDefined,
|
|
"Input dataset has no nodata value. "
|
|
"Ignoring 'nv' entry in color palette");
|
|
CSLDestroy(papszFields);
|
|
continue;
|
|
}
|
|
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);
|
|
CPLFree(pasColorAssociation);
|
|
*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 = GDALGetRasterBandYSize(hSrcBand);
|
|
if (eDT == GDT_Byte || ((eDT == GDT_Int16 || eDT == GDT_UInt16) &&
|
|
static_cast<GIntBig>(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;
|
|
}
|
|
|
|
// Implements Wilson et al. (2007), for bathymetric use cases
|
|
template <class T>
|
|
static float GDALTRIAlgWilson(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;
|
|
}
|
|
|
|
// Implements Riley, S.J., De Gloria, S.D., Elliot, R. (1999): A Terrain
|
|
// Ruggedness that Quantifies Topographic Heterogeneity. Intermountain Journal
|
|
// of Science, Vol.5, No.1-4, pp.23-27 for terrestrial use cases
|
|
template <class T>
|
|
static float GDALTRIAlgRiley(const T *afWin, float /*fDstNoDataValue*/,
|
|
void * /*pData*/)
|
|
{
|
|
const auto square = [](double x) { return x * x; };
|
|
|
|
return static_cast<float>(
|
|
std::sqrt(square(afWin[0] - afWin[4]) + square(afWin[1] - afWin[4]) +
|
|
square(afWin[2] - afWin[4]) + square(afWin[3] - afWin[4]) +
|
|
square(afWin[5] - afWin[4]) + square(afWin[6] - afWin[4]) +
|
|
square(afWin[7] - afWin[4]) + square(afWin[8] - afWin[4])));
|
|
}
|
|
|
|
/************************************************************************/
|
|
/* 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> *poDSIn,
|
|
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="/programs/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 pointer to a integer output variable to store if 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;
|
|
}
|
|
|
|
if (psOptions->bGradientAlgSpecified &&
|
|
!(eUtilityMode == HILL_SHADE || eUtilityMode == SLOPE ||
|
|
eUtilityMode == ASPECT))
|
|
{
|
|
CPLError(CE_Failure, CPLE_IllegalArg,
|
|
"This value of -alg is only value for hillshade, slope or "
|
|
"aspect processing");
|
|
GDALDEMProcessingOptionsFree(psOptionsToFree);
|
|
return nullptr;
|
|
}
|
|
|
|
if (psOptions->bTRIAlgSpecified && !(eUtilityMode == TRI))
|
|
{
|
|
CPLError(CE_Failure, CPLE_IllegalArg,
|
|
"This value of -alg is only value for TRI processing");
|
|
GDALDEMProcessingOptionsFree(psOptionsToFree);
|
|
return nullptr;
|
|
}
|
|
|
|
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->eGradientAlg);
|
|
if (psOptions->eGradientAlg == GradientAlg::ZEVENBERGEN_THORNE)
|
|
{
|
|
pfnAlgFloat = GDALHillshadeMultiDirectionalAlg<
|
|
float, GradientAlg::ZEVENBERGEN_THORNE>;
|
|
pfnAlgInt32 = GDALHillshadeMultiDirectionalAlg<
|
|
GInt32, GradientAlg::ZEVENBERGEN_THORNE>;
|
|
}
|
|
else
|
|
{
|
|
pfnAlgFloat =
|
|
GDALHillshadeMultiDirectionalAlg<float, GradientAlg::HORN>;
|
|
pfnAlgInt32 =
|
|
GDALHillshadeMultiDirectionalAlg<GInt32, GradientAlg::HORN>;
|
|
}
|
|
}
|
|
else if (eUtilityMode == HILL_SHADE)
|
|
{
|
|
dfDstNoDataValue = 0;
|
|
bDstHasNoData = true;
|
|
pData = GDALCreateHillshadeData(adfGeoTransform, psOptions->z,
|
|
psOptions->scale, psOptions->alt,
|
|
psOptions->az, psOptions->eGradientAlg);
|
|
if (psOptions->eGradientAlg == GradientAlg::ZEVENBERGEN_THORNE)
|
|
{
|
|
if (psOptions->bCombined)
|
|
{
|
|
pfnAlgFloat =
|
|
GDALHillshadeCombinedAlg<float,
|
|
GradientAlg::ZEVENBERGEN_THORNE>;
|
|
pfnAlgInt32 =
|
|
GDALHillshadeCombinedAlg<GInt32,
|
|
GradientAlg::ZEVENBERGEN_THORNE>;
|
|
}
|
|
else if (psOptions->bIgor)
|
|
{
|
|
pfnAlgFloat =
|
|
GDALHillshadeIgorAlg<float,
|
|
GradientAlg::ZEVENBERGEN_THORNE>;
|
|
pfnAlgInt32 =
|
|
GDALHillshadeIgorAlg<GInt32,
|
|
GradientAlg::ZEVENBERGEN_THORNE>;
|
|
}
|
|
else
|
|
{
|
|
pfnAlgFloat =
|
|
GDALHillshadeAlg<float, GradientAlg::ZEVENBERGEN_THORNE>;
|
|
pfnAlgInt32 =
|
|
GDALHillshadeAlg<GInt32, GradientAlg::ZEVENBERGEN_THORNE>;
|
|
}
|
|
}
|
|
else
|
|
{
|
|
if (psOptions->bCombined)
|
|
{
|
|
pfnAlgFloat =
|
|
GDALHillshadeCombinedAlg<float, GradientAlg::HORN>;
|
|
pfnAlgInt32 =
|
|
GDALHillshadeCombinedAlg<GInt32, GradientAlg::HORN>;
|
|
}
|
|
else if (psOptions->bIgor)
|
|
{
|
|
pfnAlgFloat = GDALHillshadeIgorAlg<float, GradientAlg::HORN>;
|
|
pfnAlgInt32 = GDALHillshadeIgorAlg<GInt32, GradientAlg::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, GradientAlg::HORN>;
|
|
pfnAlgInt32 = GDALHillshadeAlg<GInt32, GradientAlg::HORN>;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
else if (eUtilityMode == SLOPE)
|
|
{
|
|
dfDstNoDataValue = -9999;
|
|
bDstHasNoData = true;
|
|
|
|
pData = GDALCreateSlopeData(adfGeoTransform, psOptions->scale,
|
|
psOptions->slopeFormat);
|
|
if (psOptions->eGradientAlg == GradientAlg::ZEVENBERGEN_THORNE)
|
|
{
|
|
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->eGradientAlg == GradientAlg::ZEVENBERGEN_THORNE)
|
|
{
|
|
pfnAlgFloat = GDALAspectZevenbergenThorneAlg<float>;
|
|
pfnAlgInt32 = GDALAspectZevenbergenThorneAlg<GInt32>;
|
|
}
|
|
else
|
|
{
|
|
pfnAlgFloat = GDALAspectAlg<float>;
|
|
pfnAlgInt32 = GDALAspectAlg<GInt32>;
|
|
}
|
|
}
|
|
else if (eUtilityMode == TRI)
|
|
{
|
|
dfDstNoDataValue = -9999;
|
|
bDstHasNoData = true;
|
|
if (psOptions->eTRIAlg == TRIAlg::WILSON)
|
|
{
|
|
pfnAlgFloat = GDALTRIAlgWilson<float>;
|
|
pfnAlgInt32 = GDALTRIAlgWilson<GInt32>;
|
|
}
|
|
else
|
|
{
|
|
pfnAlgFloat = GDALTRIAlgRiley<float>;
|
|
pfnAlgInt32 = GDALTRIAlgRiley<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="/programs/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 = new GDALDEMProcessingOptions;
|
|
Algorithm eUtilityMode = INVALID;
|
|
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->bGradientAlgSpecified = true;
|
|
psOptions->eGradientAlg = GradientAlg::ZEVENBERGEN_THORNE;
|
|
}
|
|
else if (EQUAL(papszArgv[i], "Horn"))
|
|
{
|
|
psOptions->bGradientAlgSpecified = true;
|
|
psOptions->eGradientAlg = GradientAlg::HORN;
|
|
}
|
|
else if (EQUAL(papszArgv[i], "Riley"))
|
|
{
|
|
psOptions->bTRIAlgSpecified = true;
|
|
psOptions->eTRIAlg = TRIAlg::RILEY;
|
|
}
|
|
else if (EQUAL(papszArgv[i], "Wilson"))
|
|
{
|
|
psOptions->bTRIAlgSpecified = true;
|
|
psOptions->eTRIAlg = TRIAlg::WILSON;
|
|
}
|
|
else
|
|
{
|
|
CPLError(CE_Failure, CPLE_IllegalArg,
|
|
"Invalid value for -alg: %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);
|
|
|
|
delete 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;
|
|
}
|