MRG

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

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 );

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();