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                false)
  39 , m_litChannel(lit)
  40 , m_exitDepth(exitDepth)
  41 {
  42   m_litChannel->addDerivitive(this);
  43 }
  44
  45 /// Destructor.
  46 tcShadowDepth::~tcShadowDepth()
  47 {
  48   m_litChannel->removeDerivitive(this);
  49 }
  50
  51 /*
  52  * Interface for derived class to implement
  53  */
  54
  55 #include <QtDebug>
  56 void tcShadowDepth::roundPortion(double* x1, double* y1, double* x2, double* y2)
  57 {
  58   m_litChannel->roundPortion(x1,y1,x2,y2);
  59 }
  60
  61 tcAbstractPixelData* tcShadowDepth::loadPortion(double x1, double y1, double x2, double y2, bool changed)
  62 {
  63   Reference<tcAbstractPixelData> channelData = m_litChannel->loadPortion(x1, y1, x2, y2, changed);
  64   int width = channelData->width();
  65   int height = channelData->height();
  66
  67   // Create a new pixel buffer
  68   tcPixelData<float>* data = new tcPixelData<float>(width, height);
  69   startProcessing(tr("Depthing shadows"));
  70   /*for (int i = 0; i < width*height; ++i)
  71   {
  72     data->buffer()[i] = -2.0f;
  73   }*/
  74
  75   tempData temp;
  76   temp.x1 = x1;
  77   temp.y1 = y1;
  78   temp.x2 = x2;
  79   temp.y2 = y2;
  80   temp.width = width;
  81   temp.height = height;
  82   temp.lit = channelData;
  83   temp.shadowDepths = data;
  84
  85   int maxWidthHeight = qMax(width, height);
  86   double minXY12 = qMin(x2-x1, y2-y1);
  87
  88   Reference<tcPixelData<float> > litData = dynamicCast<tcPixelData<float>*>(channelData);
  89   if (litData != 0)
  90   for (int j = 0; j < height; ++j)
  91   {
  92     progress((float)j/(height-1));
  93     for (int i = 0; i < width; ++i)
  94     {
  95       //pixelShadowDepth(&data, i, j);
  96
  97       int index = j*width + i;
  98
  99       // Not actually in shadow
 100       if (litData->buffer()[index] > 0.5f)
 101       {
 102         data->buffer()[index] = 0.0f;
 103       }
 104       // In shadow
 105       else
 106       {
 107         maths::Vector<2,float> coord      (x1 + (x2-x1)*i / (width  - 1),
 108                                            y1 + (y2-y1)*j / (height - 1));
 109         /*
 110          * Cx=x1+(x2-x1)*Tx
 111          * Cy=1-y1-(y2-y1)*Ty
 112          * Tx=(Cx-x1)/(x2-x1)
 113          * Ty=(1-y1-Cy)/(y2-y1)
 114          */
 115
 116         // Find local light direction
 117         tcGeo geoCoord = geoAt(coord);
 118         maths::Vector<3,float> light = lightDirectionAt(geoCoord);
 119
 120         // This appears to be almost accurate: the light vector in texture space
 121         // Perhaps this is the size of the texture not being taken into account
 122         tcGeo offset2(geoAt(coord + maths::Vector<2,float>(-light[0]/(maxWidthHeight-1)*minXY12,
 123                                                            light[1]/(maxWidthHeight-1)*minXY12)));
 124
 125         // This doesn't appear to be accurate: the light vector in geographical space
 126         //float magn = maths::Vector<2,float>(offset1.lon()-geoCoord.lon(), offset1.lat()-geoCoord.lat()).mag();
 127         //tcGeo offset2(geoCoord + tcGeo(light[0]*magn,
 128         //                               light[1]*magn));
 129
 130         // Trace in sun direction until we reach non shadow.
 131         tcGeo pos = geoCoord;
 132         tcGeo delta;
 133         if (m_exitDepth)
 134         {
 135           delta = offset2 - geoCoord;
 136         }
 137         else
 138         {
 139           delta = geoCoord - offset2;
 140         }
 141         pos = pos + delta;
 142         maths::Vector<2,float> tex = coord;
 143         while (tex[0] >= x1 && tex[0] <= x2 && tex[1] >= y1 && tex[1] <= y2)
 144         {
 145           if (channelData->sampleFloat((tex[0]-x1)/(x2-x1), (tex[1]-y1)/(y2-y1)) > 0.5f)
 146           {
 147             break;
 148           }
 149           pos = pos + delta;
 150           tex = textureAt(pos);
 151         }
 152         data->buffer()[index] = (tex-coord).mag()*200;
 153       }
 154     }
 155   }
 156   endProcessing();
 157   return data;
 158 }
 159
 160 /*
 161  * Private functions
 162  */
 163
 164 /// Calculate the shadow depth of a pixel, recursively as necessary.
 165 float tcShadowDepth::pixelShadowDepth(const tempData* data, int i, int j)
 166 {
 167   int index = j*data->width + i;
 168   float& depth = data->shadowDepths->buffer()[index];
 169
 170   // Already set, just return that
 171   if (depth >= 1.0f)
 172   {
 173     return depth;
 174   }
 175
 176   // Not actually in shadow
 177   if (data->lit->sampleFloat(i, j) > 0.5f)
 178   {
 179     return 0.0f;
 180   }
 181
 182   // In shadow, not yet calculated depth
 183
 184   // Find local light direction
 185   maths::Vector<2,float> coord      (data->x1 + (data->x2-data->x1)*i / (data->width  - 1),
 186                                      data->y1 + (data->y2-data->y1)*j / (data->height - 1));
 187   tcGeo geoCoord = geoAt(coord);
 188   maths::Vector<3,float> light = lightDirectionAt(geoCoord);
 189
 190   int coords[2][2];
 191   float weights[2];
 192
 193   // Pick an axis direction
 194   if (light[1] < light[0]) // SE
 195   {
 196     if (-light[1] < light[0]) // E
 197     {
 198       coords[0][0] = i+1; coords[0][1] = j;
 199       weights[1] = fabs(light[1]/light[0]);
 200     }
 201     else // S
 202     {
 203       coords[0][0] = i; coords[0][1] = j+1;
 204       weights[1] = fabs(light[0]/light[1]);
 205     }
 206   }
 207   else // NW
 208   {
 209     if (-light[1] < light[0]) // W
 210     {
 211       coords[0][0] = i-1; coords[0][1] = j;
 212       weights[1] = fabs(light[1]/light[0]);
 213     }
 214     else // N
 215     {
 216       coords[0][0] = i; coords[0][1] = j-1;
 217       weights[1] = fabs(light[0]/light[1]);
 218     }
 219   }
 220
 221   // Pick a diagonal direction
 222   if (light[0] > 0.0f) // E (NE,SE)
 223   {
 224     if (light[1] > 0.0f) // NE
 225     {
 226       coords[1][0] = i+1; coords[1][1] = j-1;
 227     }
 228     else // SE
 229     {
 230       coords[1][0] = i+1; coords[1][1] = j+1;
 231     }
 232   }
 233   else // W (NW,SW)
 234   {
 235     if (light[1] > 0.0f) // NW
 236     {
 237       coords[1][0] = i-1; coords[1][1] = j+1;
 238     }
 239     else // SW
 240     {
 241       coords[1][0] = i-1; coords[1][1] = j+1;
 242     }
 243   }
 244   weights[0] = 1.0f - weights[1];
 245
 246   // Get depths of two points.
 247   float depths[2];
 248   for (int ii = 0; ii < 2; ++ii)
 249   {
 250     if (coords[ii][0] >= 0 && coords[ii][0] < data->width &&
 251         coords[ii][1] >= 0 && coords[ii][1] < data->height)
 252     {
 253       depths[ii] = pixelShadowDepth(data, coords[ii][0], coords[ii][1]);
 254     }
 255     else
 256     {
 257       depths[ii] = -1.0f;
 258     }
 259   }
 260
 261   // Each depth can be:
 262   // -1 (out of range or unknown)
 263   // zero (lit)
 264   // positive (unlit, x 2d distance from boundary)
 265
 266   // If both are unknown, result is unknown
 267   if (depths[0] == -1.0f && depths[1] == -1.0f)
 268   {
 269     depth = -1.0f;
 270   }
 271   // If one is unknown, use whichever has highest weight
 272   else if (depths[0] == -1.0f || depths[1] == -1.0f)
 273   {
 274
 275   }
 276   // Normal weighting
 277   else
 278   {
 279     depth = depths[0]*weights[0] + depths[1]*weights[1];
 280   }
 281   return depth;
 282 }