Commit 2131cc9a authored by Emmanuel Christophe's avatar Emmanuel Christophe

ENH: improve dem rendering

parent d430085e
......@@ -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
......
......@@ -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 );
......
......@@ -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();
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment