From 6704f8e2f0bf004c1baa71e7b18171ef245bd0ed Mon Sep 17 00:00:00 2001 From: James Hogan Date: Sun, 7 Dec 2008 02:13:26 +0000 Subject: [PATCH] Raytracing to determine whether a pixel is lit or not --- geo/CMakeLists.txt | 1 + geo/tcChannelDem.cpp | 14 +++++ geo/tcChannelDem.h | 6 ++ geo/tcLambertianShading.cpp | 5 +- geo/tcLandsatData.cpp | 4 ++ geo/tcNegativeProduct.cpp | 16 +++++ ...bertianShading.cpp => tcRaytracedShadowMap.cpp} | 35 +++++------ geo/{tcChannelDem.h => tcRaytracedShadowMap.h} | 53 +++++------------ geo/tcSrtmModel.cpp | 68 +++++++++++++++++++++- geo/tcSrtmModel.h | 6 ++ 10 files changed, 146 insertions(+), 62 deletions(-) copy geo/{tcLambertianShading.cpp => tcRaytracedShadowMap.cpp} (60%) copy geo/{tcChannelDem.h => tcRaytracedShadowMap.h} (53%) diff --git a/geo/CMakeLists.txt b/geo/CMakeLists.txt index 1ef0d29..c700171 100644 --- a/geo/CMakeLists.txt +++ b/geo/CMakeLists.txt @@ -28,6 +28,7 @@ set(tecorrec_geo_SRCS tcShadowClassification.cpp tcNegativeProduct.cpp tcLambertianShading.cpp + tcRaytracedShadowMap.cpp tcGeoImageData.cpp tcSensor.cpp tcSpectrum.cpp diff --git a/geo/tcChannelDem.cpp b/geo/tcChannelDem.cpp index 88ba422..78bc582 100644 --- a/geo/tcChannelDem.cpp +++ b/geo/tcChannelDem.cpp @@ -92,3 +92,17 @@ maths::Vector<3,float> tcChannelDem::lightDirectionAt(const tcGeo& coordinate) c /// @todo Improve me to vary a bit across the image return m_imagery->sunDirection(); } + +/// Find whether a point is lit by the light. +bool tcChannelDem::isLitAt(const tcGeo& coordinate) const +{ + return m_dem->litAt(coordinate, lightDirectionAt(coordinate)) > 0.5f; +} + +/// Find how lit a point is by the light. +float tcChannelDem::litAt(const tcGeo& coordinate) const +{ + // Sun is half a degree in the sky + static const float sunSize = 0.25f * M_PI / 180.0f; + return m_dem->litAt(coordinate, lightDirectionAt(coordinate), sunSize); +} diff --git a/geo/tcChannelDem.h b/geo/tcChannelDem.h index baaedca..8d4f1e1 100644 --- a/geo/tcChannelDem.h +++ b/geo/tcChannelDem.h @@ -72,6 +72,12 @@ class tcChannelDem : public tcChannel /// Get the lighting direction vector at a geographical coordinate. maths::Vector<3,float> lightDirectionAt(const tcGeo& coordinate) const; + /// Find whether a point is lit by the light. + bool isLitAt(const tcGeo& coordinate) const; + + /// Find how lit a point is by the light. + float litAt(const tcGeo& coordinate) const; + private: /* diff --git a/geo/tcLambertianShading.cpp b/geo/tcLambertianShading.cpp index e75f476..b36d867 100644 --- a/geo/tcLambertianShading.cpp +++ b/geo/tcLambertianShading.cpp @@ -57,13 +57,12 @@ void tcLambertianShading::roundPortion(double* x1, double* y1, double* x2, doubl tcAbstractPixelData* tcLambertianShading::loadPortion(double x1, double y1, double x2, double y2) { - tcPixelData* data = 0; - Reference > channelData = dynamic_cast*>(m_referenceChannel->loadPortion(x1, y1, x2, y2)); + Reference channelData = m_referenceChannel->loadPortion(x1, y1, x2, y2); int width = channelData->width(); int height = channelData->height(); // Create a new pixel buffer - data = new tcPixelData(width, height); + tcPixelData* data = new tcPixelData(width, height); for (int j = 0; j < height; ++j) { for (int i = 0; i < width; ++i) diff --git a/geo/tcLandsatData.cpp b/geo/tcLandsatData.cpp index bf9b234..f8d8b8c 100644 --- a/geo/tcLandsatData.cpp +++ b/geo/tcLandsatData.cpp @@ -33,6 +33,7 @@ #include "tcShadowClassification.h" #include "tcNegativeProduct.h" #include "tcLambertianShading.h" +#include "tcRaytracedShadowMap.h" #include #include @@ -166,6 +167,9 @@ tcLandsatData::tcLandsatData(const QString& path, tcSrtmModel* dem) tcLambertianShading* diffuse = new tcLambertianShading(channelManager()->channel(0), dem, this); channelManager()->addChannel(diffuse); + tcRaytracedShadowMap* raytrace = new tcRaytracedShadowMap(channelManager()->channel(0), dem, this); + channelManager()->addChannel(raytrace); + // Only load one metafile break; } diff --git a/geo/tcNegativeProduct.cpp b/geo/tcNegativeProduct.cpp index 5ba2e34..2ccd0ff 100644 --- a/geo/tcNegativeProduct.cpp +++ b/geo/tcNegativeProduct.cpp @@ -94,5 +94,21 @@ tcAbstractPixelData* tcNegativeProduct::loadPortion(double x1, double y1, double data->buffer()[i] *= 1.0f - (float)channelData->buffer()[i] / 255.0f; } } +#if 0 + if (0 != data) + { + for (int i = 0; i < width*height; ++i) + { + if (data->buffer()[i] > 0.01f) + { + data->buffer()[i] = 0.0f; + } + else + { + data->buffer()[i] = 1.0f; + } + } + } +#endif return data; } diff --git a/geo/tcLambertianShading.cpp b/geo/tcRaytracedShadowMap.cpp similarity index 60% copy from geo/tcLambertianShading.cpp copy to geo/tcRaytracedShadowMap.cpp index e75f476..8b5a8ca 100644 --- a/geo/tcLambertianShading.cpp +++ b/geo/tcRaytracedShadowMap.cpp @@ -18,11 +18,11 @@ ***************************************************************************/ /** - * @file tcLambertianShading.cpp - * @brief Shading channel using lambertian reflectance. + * @file tcRaytracedShadowMap.cpp + * @brief Shadow map using ray tracing of DEM. */ -#include "tcLambertianShading.h" +#include "tcRaytracedShadowMap.h" #include @@ -31,18 +31,16 @@ */ /// Primary constructor. -tcLambertianShading::tcLambertianShading(tcChannel* reference, tcSrtmModel* dem, tcGeoImageData* imagery) +tcRaytracedShadowMap::tcRaytracedShadowMap(tcChannel* reference, tcSrtmModel* dem, tcGeoImageData* imagery) : tcChannelDem(dem, imagery, - QObject::tr("Lambertian Shading"), - QObject::tr("Uses lambertian shading and elevation data to shade pixels. " - "Cast shadows are not taken into account. " - "Self shadows are set to zero.")) + QObject::tr("Raytraced Shadow Map"), + QObject::tr("Uses ray tracing and elevation data to highlight shadows.")) , m_referenceChannel(reference) { } /// Destructor. -tcLambertianShading::~tcLambertianShading() +tcRaytracedShadowMap::~tcRaytracedShadowMap() { } @@ -50,20 +48,19 @@ tcLambertianShading::~tcLambertianShading() * Interface for derived class to implement */ -void tcLambertianShading::roundPortion(double* x1, double* y1, double* x2, double* y2) +void tcRaytracedShadowMap::roundPortion(double* x1, double* y1, double* x2, double* y2) { m_referenceChannel->roundPortion(x1,y1,x2,y2); } -tcAbstractPixelData* tcLambertianShading::loadPortion(double x1, double y1, double x2, double y2) +tcAbstractPixelData* tcRaytracedShadowMap::loadPortion(double x1, double y1, double x2, double y2) { - tcPixelData* data = 0; - Reference > channelData = dynamic_cast*>(m_referenceChannel->loadPortion(x1, y1, x2, y2)); + Reference channelData = m_referenceChannel->loadPortion(x1, y1, x2, y2); int width = channelData->width(); int height = channelData->height(); // Create a new pixel buffer - data = new tcPixelData(width, height); + tcPixelData* data = new tcPixelData(width, height); for (int j = 0; j < height; ++j) { for (int i = 0; i < width; ++i) @@ -72,14 +69,10 @@ tcAbstractPixelData* tcLambertianShading::loadPortion(double x1, double y1, doub // Transform coordinates maths::Vector<2,float> coord ( x1 + (x2-x1)*i / (width - 1), 1.0f - y1 - (y2-y1)*j / (height - 1)); - tcGeo geoCoord = geoAt(coord); + tcGeo geoCoord = geoAt(coord); + float lit = litAt(geoCoord); - // Get some elevation data - maths::Vector<3,float> normal = normalAt(geoCoord); - maths::Vector<3,float> light = lightDirectionAt(geoCoord); - // Find dot product between normal and light vectors and limit it to 0.0 - float lambert = normal * light; - data->buffer()[index] = (lambert >= 0.0f ? lambert : 0.0f); + data->buffer()[index] = lit; } } return data; diff --git a/geo/tcChannelDem.h b/geo/tcRaytracedShadowMap.h similarity index 53% copy from geo/tcChannelDem.h copy to geo/tcRaytracedShadowMap.h index baaedca..64890d8 100644 --- a/geo/tcChannelDem.h +++ b/geo/tcRaytracedShadowMap.h @@ -17,23 +17,18 @@ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. * ***************************************************************************/ -#ifndef _tcChannelDem_h_ -#define _tcChannelDem_h_ +#ifndef _tcRaytracedShadowMap_h_ +#define _tcRaytracedShadowMap_h_ /** - * @file tcChannelDem.h - * @brief An abstract image channel which needs access to the digital elevation model. + * @file tcRaytracedShadowMap.h + * @brief Shadow map using ray tracing of DEM. */ -#include "tcChannel.h" -#include "tcGeo.h" -#include +#include "tcChannelDem.h" -class tcSrtmModel; -class tcGeoImageData; - -/// An abstract image channel which needs access to the digital elevation model. -class tcChannelDem : public tcChannel +/// Shadow map using ray tracing of DEM. +class tcRaytracedShadowMap : public tcChannelDem { public: @@ -42,35 +37,22 @@ class tcChannelDem : public tcChannel */ /// Primary constructor. - tcChannelDem(tcSrtmModel* dem, tcGeoImageData* imagery, const QString& name, const QString& description); + tcRaytracedShadowMap(tcChannel* reference, tcSrtmModel* dem, tcGeoImageData* imagery); /// Destructor. - virtual ~tcChannelDem(); + virtual ~tcRaytracedShadowMap(); protected: /* - * Interface for derived classes to access elevation model data + * Interface for derived class to implement */ - // Main transformation - - /// Get the geographical coordinate at a pixel. - tcGeo geoAt(const maths::Vector<2,float>& textureCoordinate) const; - - /// Get the texture coordinate at a geographical coordinate. - maths::Vector<2,float> textureAt(const tcGeo& geoCoordinate) const; - - // Wrapper functions + // Reimplemented + virtual void roundPortion(double* x1, double* y1, double* x2, double* y2); - /// Get the altitude at a geographical coordinate. - float altitudeAt(const tcGeo& coordinate) const; - - /// Get the texture-space normal vector at a geographical coordinate. - maths::Vector<3,float> normalAt(const tcGeo& coordinate) const; - - /// Get the lighting direction vector at a geographical coordinate. - maths::Vector<3,float> lightDirectionAt(const tcGeo& coordinate) const; + // Reimplemented + virtual tcAbstractPixelData* loadPortion(double x1, double y1, double x2, double y2); private: @@ -78,11 +60,8 @@ class tcChannelDem : public tcChannel * Variables */ - /// Digital elevation model. - tcSrtmModel* m_dem; - - /// Imagery object. - tcGeoImageData* m_imagery; + /// We'll use the same resolution as this reference channel. + tcChannel* m_referenceChannel; }; #endif diff --git a/geo/tcSrtmModel.cpp b/geo/tcSrtmModel.cpp index 9a9a7ef..c09f936 100644 --- a/geo/tcSrtmModel.cpp +++ b/geo/tcSrtmModel.cpp @@ -153,13 +153,13 @@ double tcSrtmModel::altitudeAt(const tcGeo& coord, bool corrected, bool* accurat /// Get the normal at a coordinate. maths::Vector<3,float> tcSrtmModel::normalAt(const tcGeo& coord, bool corrected, bool* accurate) { - const double radiusEarth = 6378.137e3; if (0 != m_cache) { tcSrtmData* data = m_cache->loadData(coord); if (0 != data) { // Sample altitudes a fixed angle in each direction. + const double radiusEarth = 6378.137e3; tcGeo angRes = data->sampleResolution(); bool xpa, xma, ypa, yma; double xp = altitudeAt(coord + tcGeo(angRes.lon()*0.5, 0.0), corrected, &xpa); @@ -181,9 +181,75 @@ maths::Vector<3,float> tcSrtmModel::normalAt(const tcGeo& coord, bool corrected, ).normalized(); } } + if (0 != accurate) + { + *accurate = false; + } return maths::Vector<3,float>(0.0f, 0.0f, 1.0f); } +/// Use raytracing to find how lit a point is. +float tcSrtmModel::litAt(const tcGeo& coord, const maths::Vector<3,float>& light, float angularRadius, bool corrected, bool* accurate) +{ + if (0 != m_cache) + { + tcSrtmData* data = m_cache->loadData(coord); + if (0 != data) + { + // Sample altitudes a fixed angle in each direction. + const float radiusEarth = 6378.137e3; + tcGeo angRes = data->sampleResolution(); + float distance = 0.0f; + float step = (angRes.lon() + angRes.lat()) / 4.0; + bool ac = true; + float result = 1.0f; + bool xa; + tcGeo cur = coord; + float alt = altitudeAt(cur, corrected, &xa); + ac = ac && xa; + Q_ASSERT(light[2] > 0.0); + while (alt < 8000.0) + { + // Move towards light source + cur.setLon(cur.lon() + step*light[0] / cos(cur.lat())); + cur.setLat(cur.lat() + step*light[1]); + alt += radiusEarth*step*light[2]; + distance += radiusEarth*step; + // Check for collision + float x = altitudeAt(cur, corrected, &xa); + ac = ac && xa; + float flex = distance*angularRadius; + if (x >= alt + flex) + { + result = 0.0f; + break; + } + else if (x >= alt - flex) + { + float portion = (x - alt + flex) / (2.0f*flex); + if (result > portion) + { + result = portion; + } + } + // Can take bigger steps the further we are from the ground + // here we are assuming that the terrain isn't extremely steep + step = 0.5f * (alt - x + flex) / radiusEarth; + } + if (0 != accurate) + { + *accurate = ac; + } + return result; + } + } + if (0 != accurate) + { + *accurate = false; + } + return 1.0f; +} + /// Set the data set to use. void tcSrtmModel::setDataSet(const QString& name) { diff --git a/geo/tcSrtmModel.h b/geo/tcSrtmModel.h index 169d0b5..2246876 100644 --- a/geo/tcSrtmModel.h +++ b/geo/tcSrtmModel.h @@ -105,6 +105,12 @@ class tcSrtmModel /// Get the normal at a coordinate. maths::Vector<3,float> normalAt(const tcGeo& coord, bool corrected = false, bool* accurate = 0); + /// Use raytracing to find whether a point is lit. + bool isLitAt(const tcGeo& coord, const maths::Vector<3,float>& light, bool corrected = false, bool* accurate = 0); + + /// Use raytracing to find how lit a point is. + float litAt(const tcGeo& coord, const maths::Vector<3,float>& light, float angularRadius = 0.0f, bool corrected = false, bool* accurate = 0); + /// Set the data set to use. void setDataSet(const QString& name); -- 2.11.4.GIT