geo/tcShadowDepth.cpp

   1 /***************************************************************************
   2  *   This file is part of Tecorrec.                                        *
   3  *   Copyright 2008 James Hogan <james@albanarts.com>                      *
   4  *                                                                         *
   5  *   Tecorrec is free software: you can redistribute it and/or modify      *
   6  *   it under the terms of the GNU General Public License as published by  *
   7  *   the Free Software Foundation, either version 2 of the License, or     *
   8  *   (at your option) any later version.                                   *
   9  *                                                                         *
  10  *   Tecorrec is distributed in the hope that it will be useful,           *
  11  *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
  12  *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
  13  *   GNU General Public License for more details.                          *
  14  *                                                                         *
  15  *   You should have received a copy of the GNU General Public License     *
  16  *   along with Tecorrec.  If not, write to the Free Software Foundation,  *
  17  *   Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.   *
  18  ***************************************************************************/
  19
  20 /**
  21  * @file tcShadowDepth.cpp
  22  * @brief 2D depth into shadow along light direction.
  23  */
  24
  25 #include "tcShadowDepth.h"
  26
  27 #include <QObject>
  28
  29 /*
  30  * Constructors + destructor
  31  */
  32
  33 /// Primary constructor.
  34 tcShadowDepth::tcShadowDepth(tcChannel* lit, tcSrtmModel* dem, tcGeoImageData* imagery, bool exitDepth)
  35 : tcChannelDem(dem, imagery,
  36                tr("Shadow depth"),
  37                tr("2D depth into shadows along light direction."))
  38 , m_litChannel(lit)
  39 , m_exitDepth(exitDepth)
  40 {
  41   m_litChannel->addDerivitive(this);
  42 }
  43
  44 /// Destructor.
  45 tcShadowDepth::~tcShadowDepth()
  46 {
  47   m_litChannel->removeDerivitive(this);
  48 }
  49
  50 /*
  51  * Interface for derived class to implement
  52  */
  53
  54 void tcShadowDepth::roundPortion(double* x1, double* y1, double* x2, double* y2)
  55 {
  56   m_litChannel->roundPortion(x1,y1,x2,y2);
  57 }
  58
  59 tcAbstractPixelData* tcShadowDepth::loadPortion(double x1, double y1, double x2, double y2)
  60 {
  61   Reference<tcAbstractPixelData> channelData = m_litChannel->loadPortion(x1, y1, x2, y2);
  62   int width = channelData->width();
  63   int height = channelData->height();
  64
  65   // Create a new pixel buffer
  66   tcPixelData<float>* data = new tcPixelData<float>(width, height);
  67   startProcessing(tr("Depthing shadows"));
  68   /*for (int i = 0; i < width*height; ++i)
  69   {
  70     data->buffer()[i] = -2.0f;
  71   }*/
  72
  73   tempData temp;
  74   temp.x1 = x1;
  75   temp.y1 = y1;
  76   temp.x2 = x2;
  77   temp.y2 = y2;
  78   temp.width = width;
  79   temp.height = height;
  80   temp.lit = channelData;
  81   temp.shadowDepths = data;
  82
  83   int maxWidthHeight = qMax(width, height);
  84   double minXY12 = qMin(x2-x1, y2-y1);
  85
  86   Reference<tcPixelData<float> > litData = dynamicCast<tcPixelData<float>*>(channelData);
  87   if (litData != 0)
  88   for (int j = 0; j < height; ++j)
  89   {
  90     progress((float)j/(height-1));
  91     for (int i = 0; i < width; ++i)
  92     {
  93       //pixelShadowDepth(&data, i, j);
  94
  95       int index = j*width + i;
  96
  97       // Not actually in shadow
  98       if (litData->buffer()[index] > 0.5f)
  99       {
 100         data->buffer()[index] = 0.0f;
 101       }
 102       // In shadow
 103       else
 104       {
 105         maths::Vector<2,float> coord      (       x1 + (x2-x1)*i / (width  - 1),
 106                                            1.0f - y1 - (y2-y1)*j / (height - 1));
 107         /*
 108          * Cx=x1+(x2-x1)*Tx
 109          * Cy=1-y1-(y2-y1)*Ty
 110          * Tx=(Cx-x1)/(x2-x1)
 111          * Ty=(1-y1-Cy)/(y2-y1)
 112          */
 113
 114         // Find local light direction
 115         tcGeo geoCoord = geoAt(coord);
 116         maths::Vector<3,float> light = lightDirectionAt(geoCoord);
 117
 118         // This appears to be almost accurate: the light vector in texture space
 119         // Perhaps this is the size of the texture not being taken into account
 120         tcGeo offset2(geoAt(coord + maths::Vector<2,float>(light[0]/(maxWidthHeight-1)*minXY12,
 121                                                            light[1]/(maxWidthHeight-1)*minXY12)));
 122
 123         // This doesn't appear to be accurate: the light vector in geographical space
 124         //float magn = maths::Vector<2,float>(offset1.lon()-geoCoord.lon(), offset1.lat()-geoCoord.lat()).mag();
 125         //tcGeo offset2(geoCoord + tcGeo(light[0]*magn,
 126         //                               light[1]*magn));
 127
 128         // Trace in sun direction until we reach non shadow.
 129         tcGeo pos = geoCoord;
 130         tcGeo delta;
 131         if (m_exitDepth)
 132         {
 133           delta = offset2 - geoCoord;
 134         }
 135         else
 136         {
 137           delta = geoCoord - offset2;
 138         }
 139         pos = pos + delta;
 140         maths::Vector<2,float> tex = coord;
 141         while (tex[0] >= x1 && tex[0] <= x2 && tex[1] >= 1.0f-y2 && tex[1] <= 1.0f-y1)
 142         {
 143           if (channelData->sampleFloat((tex[0]-x1)/(x2-x1), (1.0f-tex[1]-y1)/(y2-y1)) > 0.5f)
 144           {
 145             break;
 146           }
 147           pos = pos + delta;
 148           tex = textureAt(pos);
 149         }
 150         data->buffer()[index] = (tex-coord).mag()*200;
 151       }
 152     }
 153   }
 154   endProcessing();
 155   return data;
 156 }
 157
 158 /*
 159  * Private functions
 160  */
 161
 162 /// Calculate the shadow depth of a pixel, recursively as necessary.
 163 float tcShadowDepth::pixelShadowDepth(const tempData* data, int i, int j)
 164 {
 165   int index = j*data->width + i;
 166   float& depth = data->shadowDepths->buffer()[index];
 167
 168   // Already set, just return that
 169   if (depth >= 1.0f)
 170   {
 171     return depth;
 172   }
 173
 174   // Not actually in shadow
 175   if (data->lit->sampleFloat(i, j) > 0.5f)
 176   {
 177     return 0.0f;
 178   }
 179
 180   // In shadow, not yet calculated depth
 181
 182   // Find local light direction
 183   maths::Vector<2,float> coord      (       data->x1 + (data->x2-data->x1)*i / (data->width  - 1),
 184                                      1.0f - data->y1 - (data->y2-data->y1)*j / (data->height - 1));
 185   tcGeo geoCoord = geoAt(coord);
 186   maths::Vector<3,float> light = lightDirectionAt(geoCoord);
 187
 188   int coords[2][2];
 189   float weights[2];
 190
 191   // Pick an axis direction
 192   if (light[1] < light[0]) // SE
 193   {
 194     if (-light[1] < light[0]) // E
 195     {
 196       coords[0][0] = i+1; coords[0][1] = j;
 197       weights[1] = fabs(light[1]/light[0]);
 198     }
 199     else // S
 200     {
 201       coords[0][0] = i; coords[0][1] = j+1;
 202       weights[1] = fabs(light[0]/light[1]);
 203     }
 204   }
 205   else // NW
 206   {
 207     if (-light[1] < light[0]) // W
 208     {
 209       coords[0][0] = i-1; coords[0][1] = j;
 210       weights[1] = fabs(light[1]/light[0]);
 211     }
 212     else // N
 213     {
 214       coords[0][0] = i; coords[0][1] = j-1;
 215       weights[1] = fabs(light[0]/light[1]);
 216     }
 217   }
 218
 219   // Pick a diagonal direction
 220   if (light[0] > 0.0f) // E (NE,SE)
 221   {
 222     if (light[1] > 0.0f) // NE
 223     {
 224       coords[1][0] = i+1; coords[1][1] = j-1;
 225     }
 226     else // SE
 227     {
 228       coords[1][0] = i+1; coords[1][1] = j+1;
 229     }
 230   }
 231   else // W (NW,SW)
 232   {
 233     if (light[1] > 0.0f) // NW
 234     {
 235       coords[1][0] = i-1; coords[1][1] = j+1;
 236     }
 237     else // SW
 238     {
 239       coords[1][0] = i-1; coords[1][1] = j+1;
 240     }
 241   }
 242   weights[0] = 1.0f - weights[1];
 243
 244   // Get depths of two points.
 245   float depths[2];
 246   for (int ii = 0; ii < 2; ++ii)
 247   {
 248     if (coords[ii][0] >= 0 && coords[ii][0] < data->width &&
 249         coords[ii][1] >= 0 && coords[ii][1] < data->height)
 250     {
 251       depths[ii] = pixelShadowDepth(data, coords[ii][0], coords[ii][1]);
 252     }
 253     else
 254     {
 255       depths[ii] = -1.0f;
 256     }
 257   }
 258
 259   // Each depth can be:
 260   // -1 (out of range or unknown)
 261   // zero (lit)
 262   // positive (unlit, x 2d distance from boundary)
 263
 264   // If both are unknown, result is unknown
 265   if (depths[0] == -1.0f && depths[1] == -1.0f)
 266   {
 267     depth = -1.0f;
 268   }
 269   // If one is unknown, use whichever has highest weight
 270   else if (depths[0] == -1.0f || depths[1] == -1.0f)
 271   {
 272
 273   }
 274   // Normal weighting
 275   else
 276   {
 277     depth = depths[0]*weights[0] + depths[1]*weights[1];
 278   }
 279   return depth;
 280 }