From 2131cc9a0f7a5790c3cbe983236f7336aaa98640 Mon Sep 17 00:00:00 2001
From: Emmanuel Christophe <emmanuel.christophe@orfeo-toolbox.org>
Date: Wed, 3 Mar 2010 22:49:57 +0800
Subject: [PATCH] ENH: improve dem rendering
---
Code/BasicFilters/otbHillShadingFunctor.h | 18 ++++++++++++------
Code/BasicFilters/otbReliefColormapFunctor.txx | 8 ++++----
Examples/BasicFilters/HillShadingExample.cxx | 13 +++++++++++++
3 files changed, 29 insertions(+), 10 deletions(-)
diff --git a/Code/BasicFilters/otbHillShadingFunctor.h b/Code/BasicFilters/otbHillShadingFunctor.h
index 35c826c345..9974b4667b 100644
--- a/Code/BasicFilters/otbHillShadingFunctor.h
+++ b/Code/BasicFilters/otbHillShadingFunctor.h
@@ -40,9 +40,10 @@ public:
typedef HillShadingFunctor Self;
typedef TNeighIter IteratorType;
+ typedef typename IteratorType::PixelType PixelType;
HillShadingFunctor(): m_AzimuthLight(30.0/180.0*CONST_PI), m_ElevationLight(45.0/180.0*CONST_PI),
- m_XRes(100.0), m_YRes(100.0), m_Scale(1.0)
+ m_XRes(100.0), m_YRes(100.0), m_Scale(0.1)
{
m_SinElev = vcl_sin(m_ElevationLight);
m_CosElev = vcl_cos(m_ElevationLight);
@@ -93,16 +94,16 @@ public:
const typename IteratorType::OffsetType LEFTDOWN ={{-1,1}};
// const typename IteratorType::OffsetType CENTER ={{0,0}};
- double xSlope = ((it.GetPixel(LEFTUP)+2*it.GetPixel(LEFT)+it.GetPixel(RIGHTDOWN))
- -(it.GetPixel(RIGHTUP)+2*it.GetPixel(RIGHT)+it.GetPixel(RIGHTDOWN)))
+ float xSlope = ((makeValid(it.GetPixel(LEFTUP))+2*makeValid(it.GetPixel(LEFT))+makeValid(it.GetPixel(RIGHTDOWN)))
+ -(makeValid(it.GetPixel(RIGHTUP))+2*makeValid(it.GetPixel(RIGHT))+makeValid(it.GetPixel(RIGHTDOWN))))
/(m_XRes*m_Scale);
// - as the azimuth is given compared to y axis pointing up
- double ySlope = -((it.GetPixel(LEFTUP)+2*it.GetPixel(UP)+it.GetPixel(RIGHTUP))
- -(it.GetPixel(LEFTDOWN)+2*it.GetPixel(DOWN)+it.GetPixel(RIGHTDOWN)))
+ float ySlope = -((makeValid(it.GetPixel(LEFTUP))+2*makeValid(it.GetPixel(UP))+makeValid(it.GetPixel(RIGHTUP)))
+ -(makeValid(it.GetPixel(LEFTDOWN))+2*makeValid(it.GetPixel(DOWN))+makeValid(it.GetPixel(RIGHTDOWN))))
/(m_YRes*m_Scale);
// permutation between x and y as the azimuth angle is given compared to the north-south axis
- double lambertian = ((m_SinElev*m_CosAz*ySlope)+(m_SinElev*m_SinAz*xSlope)+m_CosElev)
+ float lambertian = ((m_SinElev*m_CosAz*ySlope)+(m_SinElev*m_SinAz*xSlope)+m_CosElev)
/ vcl_sqrt(xSlope*xSlope+ySlope*ySlope+1);
return (lambertian+1)/2; //normalize between 0 and 1
@@ -113,6 +114,11 @@ private:
HillShadingFunctor(const Self&); //purposely not implemented
void operator=(const Self&); //purposely not implemented
+ inline PixelType makeValid(PixelType v) const
+ {
+ return v < itk::NumericTraits<PixelType>::Zero ? itk::NumericTraits<PixelType>::Zero : v;
+ }
+
double m_AzimuthLight; // in radian
double m_ElevationLight; // in radian
double m_XRes; // assumed to be positive provided in m
diff --git a/Code/BasicFilters/otbReliefColormapFunctor.txx b/Code/BasicFilters/otbReliefColormapFunctor.txx
index 2a88e05f6e..6584fcbb60 100644
--- a/Code/BasicFilters/otbReliefColormapFunctor.txx
+++ b/Code/BasicFilters/otbReliefColormapFunctor.txx
@@ -29,10 +29,10 @@ typename ReliefColormapFunctor<TScalar, TRGBPixel>::RGBPixelType
ReliefColormapFunctor<TScalar, TRGBPixel>
::operator()( const TScalar & v ) const
{
- float m_Borders[] = {0.0, 0.2, 0.43, 0.71, 1.0};
- float m_RedValues[] = {0.2, 0.94, 0.74, 0.92, 1.0};
- float m_GreenValues[] = {0.7, 0.98, 0.72, 0.86, 1.0};
- float m_BlueValues[] = {0.2, 0.59, 0.53, 0.69, 1.0};
+ float m_Borders[] = {0.0, 0.0001, 0.2, 0.43, 0.71, 1.0};
+ float m_RedValues[] = {0.7, 0.2, 0.94, 0.74, 0.92, 1.0};
+ float m_GreenValues[] = {0.7, 0.7, 0.98, 0.72, 0.86, 1.0};
+ float m_BlueValues[] = {1.0, 0.2, 0.59, 0.53, 0.69, 1.0};
// Map the input scalar between [0, 1].
RealType value = this->RescaleInputValue( v );
diff --git a/Examples/BasicFilters/HillShadingExample.cxx b/Examples/BasicFilters/HillShadingExample.cxx
index 8c033d5c45..84cf049dda 100644
--- a/Examples/BasicFilters/HillShadingExample.cxx
+++ b/Examples/BasicFilters/HillShadingExample.cxx
@@ -139,6 +139,16 @@ int main(int argc, char * argv[])
demToImage->SetOutputSpacing(spacing);
+ //Compute the resolution (Vincenty formula)
+ double lon1 = origin[0];
+ double lon2 = origin[0]+size[0]*spacing[0];
+ double lat1 = origin[1];
+ double lat2 = origin[1]+size[1]*spacing[1];
+ double R = 6371; // km
+ double d = vcl_acos(vcl_sin(lat1)*vcl_sin(lat2) +
+ vcl_cos(lat1)*vcl_cos(lat2) * vcl_cos(lon2-lon1)) * R;
+ double res = d / vcl_sqrt(2.0);
+
// Software Guide : BeginLatex
//
// After generating the dem image as in the DEMToImageGenerator example, you can declare
@@ -158,6 +168,9 @@ int main(int argc, char * argv[])
hillShading->SetInput(demToImage->GetOutput());
// Software Guide : EndCodeSnippet
+ hillShading->GetFunctor().SetXRes(res);
+ hillShading->GetFunctor().SetYRes(res);
+
typedef itk::RescaleIntensityImageFilter<ImageType, ScalarImageType> RescalerType;
RescalerType::Pointer rescaler = RescalerType::New();
--
GitLab