geo/tcElevationOptimization.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 tcElevationOptimization.cpp
  22  * @brief Optimized elevation.
  23  */
  24
  25 #include "tcElevationOptimization.h"
  26 #include "tcChannelConfigWidget.h"
  27 #include "tcSrtmModel.h"
  28
  29 #include <QObject>
  30 #include <QWidget>
  31 #include <QVBoxLayout>
  32 #include <QHBoxLayout>
  33 #include <QPushButton>
  34 #include <QSpinBox>
  35 #include <QMessageBox>
  36
  37 /*
  38  * Constructors + destructor
  39  */
  40
  41 /// Primary constructor.
  42 tcElevationOptimization::tcElevationOptimization(tcChannel* shadowChannel, tcChannel* shadingChannel, tcSrtmModel* dem, tcGeoImageData* imagery)
  43 : tcChannelDem(dem, imagery,
  44                tr("Elevation Optimization"),
  45                tr("Optimizes elevation using image data."))
  46 , m_shadowChannel(shadowChannel)
  47 , m_shadingChannel(shadingChannel)
  48 , m_elevationData(0)
  49 , m_data(0)
  50 , m_shadowData(0)
  51 , m_shadingData(0)
  52 , m_configWidget(0)
  53 , m_spinIterations(0)
  54 {
  55   m_shadowChannel->addDerivitive(this);
  56   m_shadingChannel->addDerivitive(this);
  57 }
  58
  59 /// Destructor.
  60 tcElevationOptimization::~tcElevationOptimization()
  61 {
  62   m_shadingChannel->removeDerivitive(this);
  63   delete m_configWidget;
  64   delete m_elevationData;
  65   delete m_data;
  66 }
  67
  68 /*
  69  * Main image interface
  70  */
  71
  72 tcChannelConfigWidget* tcElevationOptimization::configWidget()
  73 {
  74   if (0 == m_configWidget)
  75   {
  76     m_configWidget = new tcChannelConfigWidget();
  77     m_configWidget->setWindowTitle(tr("%1 Configuration").arg(name()));
  78
  79     QVBoxLayout* layout = new QVBoxLayout(m_configWidget);
  80
  81     // Row of optimisation controls
  82     QWidget* buttons = new QWidget(m_configWidget);
  83       layout->addWidget(buttons);
  84       QHBoxLayout* buttonsLayout = new QHBoxLayout(buttons);
  85
  86       // reset DEM (sets all corrections to the original (interpolated) values
  87       QPushButton* btnResetDem = new QPushButton(tr("Reset DEM"), m_configWidget);
  88         buttonsLayout->addWidget(btnResetDem);
  89         btnResetDem->setToolTip(tr("Reset DEM corrected values to interpolated values"));
  90         connect(btnResetDem, SIGNAL(clicked(bool)), this, SLOT(resetDem()));
  91
  92       // load base from DEM (loads DEM elevations into elevation buffer)
  93       QPushButton* btnLoadDem = new QPushButton(tr("Load DEM"), m_configWidget);
  94         buttonsLayout->addWidget(btnLoadDem);
  95         btnLoadDem->setToolTip(tr("Load corrected values from DEM into elevation field"));
  96         connect(btnLoadDem, SIGNAL(clicked(bool)), this, SLOT(loadFromDem()));
  97
  98       // number of iterations (select number of iterations to perform)
  99       m_spinIterations = new QSpinBox(m_configWidget);
 100         buttonsLayout->addWidget(m_spinIterations);
 101         m_spinIterations->setToolTip(tr("Number of iterations of optimization to perform"));
 102         m_spinIterations->setRange(1, 1000);
 103         m_spinIterations->setValue(10);
 104
 105       // optimize (iterates on elevation buffer and write back to DEM at intervals and end)
 106       QPushButton* btnOptimize = new QPushButton(tr("Optimize Iteratively"), m_configWidget);
 107         buttonsLayout->addWidget(btnOptimize);
 108         btnOptimize->setToolTip(tr("Iteratively optimize elevation field and write results back to DEM"));
 109         connect(btnOptimize, SIGNAL(clicked(bool)), this, SLOT(optimize()));
 110
 111     // Rows of weight sliders
 112   }
 113   return m_configWidget;
 114 }
 115
 116 /*
 117  * Slots
 118  */
 119
 120 /// Reset DEM.
 121 void tcElevationOptimization::resetDem()
 122 {
 123 }
 124
 125 /// Load elevation data from DEM.
 126 void tcElevationOptimization::loadFromDem()
 127 {
 128   delete m_elevationData;
 129   delete m_data;
 130   m_elevationData = 0;
 131   m_data = 0;
 132
 133   invalidate();
 134   m_configWidget->requestRedraw();
 135 }
 136
 137 /// Optimize a number of times.
 138 void tcElevationOptimization::optimize()
 139 {
 140   if (0 != m_elevationData)
 141   {
 142     const int width = m_elevationData->width();
 143     const int height = m_elevationData->height();
 144     // Now start optimising
 145     startProcessing("Optimizing elevation data");
 146     int passes = m_spinIterations->value();
 147     for (int i = 0; i < passes; ++i)
 148     {
 149       progress((float)i/passes);
 150       m_weights.variance = 0.0;//0001f;
 151       m_weights.roughness = 0.5f; //*i/passes;
 152       m_weights.litFacing = 0.01f;//30.0f;
 153       m_weights.shading = 300.0f;
 154       m_weights.transitionElev = 100.0f;
 155       m_weights.transitionTangential = 100.0f;
 156       m_weights.transitionCurvingDown = 100.0f;
 157       optimizeElevations(0, 0, width-1, height-1, 5.0f-4.5f*i/(passes));
 158
 159       if (i % 2 == 0 || i == passes-1)
 160       {
 161         invalidate();
 162         m_configWidget->requestRedraw();
 163       }
 164     }
 165     endProcessing();
 166   }
 167   else
 168   {
 169     QMessageBox::warning(m_configWidget,
 170                          tr("Elevation field not initialised."),
 171                          tr("Please ensure the optimisation channel is displayed."));
 172   }
 173 }
 174
 175 /*
 176  * Interface for derived class to implement
 177  */
 178
 179 void tcElevationOptimization::roundPortion(double* x1, double* y1, double* x2, double* y2)
 180 {
 181   m_shadowChannel->roundPortion(x1,y1,x2,y2);
 182 }
 183
 184 tcAbstractPixelData* tcElevationOptimization::loadPortion(double x1, double y1, double x2, double y2)
 185 {
 186   Reference<tcAbstractPixelData> channelData = m_shadowChannel->loadPortion(x1, y1, x2, y2);
 187   int width = channelData->width();
 188   int height = channelData->height();
 189   m_shadowData = dynamicCast<tcPixelData<float>*>(channelData);
 190   if (0 == m_shadowData)
 191   {
 192     return 0;
 193   }
 194   Reference<tcAbstractPixelData> shadingData = m_shadingChannel->loadPortion(x1, y1, x2, y2);
 195   m_shadingData = dynamicCast<tcPixelData<float>*>(shadingData);
 196
 197   /// @todo What if its a different size because the portion has changed?
 198   if (0 == m_elevationData || 0 == m_data)
 199   {
 200     delete m_elevationData;
 201     delete m_data;
 202
 203     // Create a new pixel buffer, initialised to altitude
 204
 205     m_elevationData = new tcElevationData(width, height);
 206     // Find transform of elevation data
 207     tcGeo geoBase = geoAt(maths::Vector<2,float>(x1, y1));
 208     tcGeo geoDiagonal = geoAt(maths::Vector<2,float>(x1 + (x2-x1)/(width-1), y1 + (y2-y1)/(height-1))) - geoBase;
 209     m_elevationData->setGeoToPixels(tcAffineTransform2<double>(geoBase, geoDiagonal).inverse());
 210
 211     m_data = new tcTypedPixelData<PerPixelCache>(width, height);
 212
 213     int maxWidthHeight = qMax(width, height);
 214     double minXY12 = qMin(x2-x1, y2-y1);
 215
 216     startProcessing("Caching elevation characteristics");
 217     for (int j = 0; j < height; ++j)
 218     {
 219       progress((float)j/(height-1));
 220       for (int i = 0; i < width; ++i)
 221       {
 222         int index = j*width + i;
 223         // Transform coordinates
 224         maths::Vector<2,float> coord      (x1 + (x2-x1)*i / (width  - 1),
 225                                            y1 + (y2-y1)*j / (height - 1));
 226         tcGeo                  geoCoord = geoAt(coord);
 227         maths::Vector<3,float> light = lightDirectionAt(geoCoord);
 228         PerPixelCache& cache = m_data->buffer()[index];
 229         cache.light = light;
 230
 231         // Find resolution
 232         maths::Vector<2,float> coordx      (x1 + (x2-x1)*(i+1) / (width  - 1),
 233                                             y1 + (y2-y1)*j     / (height - 1));
 234         maths::Vector<2,float> coordy      (x1 + (x2-x1)*i     / (width  - 1),
 235                                             y1 + (y2-y1)*(j+1) / (height - 1));
 236         tcGeo                  geoCoordx = geoAt(coordx) - geoCoord;
 237         tcGeo                  geoCoordy = geoAt(coordy) - geoCoord;
 238         // Find approx
 239         float res = geoCoordx.angle();
 240         cache.resolution[0] = res;
 241         cache.resolution[1] = geoCoordy.angle();
 242         cache.resolution *= 6378.137e3;
 243
 244         // Get some elevation data
 245         bool accurate;
 246         float altitude        = altitudeAt(geoCoord, true, &accurate);
 247         m_elevationData->buffer()[index] = altitude;
 248         cache.originalElevation = altitude;
 249         cache.reliability = accurate ? 1 : 0;
 250
 251         tcGeo lightScaled(light[0]*res/sin(geoCoord.lat()),
 252                           light[1]*res);
 253         tcGeo offset(geoCoord + lightScaled);
 254         for (int depthDirection = 0; depthDirection < 2; ++depthDirection)
 255         {
 256           tcGeo pos = geoCoord;
 257           tcGeo delta;
 258           if (depthDirection == TransitionEntrance)
 259           {
 260             delta = offset - geoCoord;
 261           }
 262           else
 263           {
 264             delta = geoCoord - offset;
 265           }
 266           pos = pos + delta;
 267           maths::Vector<2,float> tex = coord;
 268           cache.shadowTransition[depthDirection][0] = i;
 269           cache.shadowTransition[depthDirection][1] = j;
 270           while (tex[0] >= x1 && tex[0] <= x2 && tex[1] >= y1 && tex[1] <= y2)
 271           {
 272             if (channelData->sampleFloat((tex[0]-x1)/(x2-x1), (tex[1]-y1)/(y2-y1)) > 0.5f)
 273             {
 274               break;
 275             }
 276             cache.shadowTransition[depthDirection][0] = 0.5f + (tex[0]-x1)/(x2-x1)*(width -1);
 277             cache.shadowTransition[depthDirection][1] = 0.5f + (tex[1]-y1)/(y2-y1)*(height-1);
 278             pos = pos + delta;
 279             tex = textureAt(pos);
 280           }
 281           cache.shadowTransitionDistance[depthDirection] = (pos-geoCoord).angle()*6378.137e3;
 282         }
 283       }
 284     }
 285   }
 286
 287   // Update elevation model
 288   startProcessing("Updating elevation model");
 289   dem()->updateFromElevationData(m_elevationData);
 290
 291   // Downscale the elevation for viewing
 292   startProcessing("Downscaling for viewing");
 293   tcPixelData<GLubyte>* result = new tcPixelData<GLubyte>(width, height);
 294   int index = 0;
 295   for (int j = 0; j < height; ++j)
 296   {
 297     progress((float)j/(height-1));
 298     for (int i = 0; i < width; ++i)
 299     {
 300       result->buffer()[index] = 255.0f * m_elevationData->buffer()[index] / 4000.0f;
 301       ++index;
 302     }
 303   }
 304
 305   endProcessing();
 306
 307   return result;
 308 }
 309
 310 template <typename T>
 311 T MySqr(T t)
 312 {
 313   return t*t;
 314 }
 315
 316 /*
 317  * Private functions
 318  */
 319
 320 /// Calculate the cost for a set of pixels.
 321 float tcElevationOptimization::cost(int x1, int y1, int x2, int y2) const
 322 {
 323   int width = m_elevationData->width();
 324   int height = m_elevationData->height();
 325   if (x1 < 0)
 326   {
 327     x1 = 0;
 328   }
 329   if (x2 >= width)
 330   {
 331     x2 = width-1;
 332   }
 333   if (y1 < 0)
 334   {
 335     y1 = 0;
 336   }
 337   if (y2 >= height)
 338   {
 339     y2 = height-1;
 340   }
 341
 342   int maxWidthHeight = qMax(width, height);
 343   //double minXY12 = qMin(m_window.x2-m_window.x1, m_window.y2-m_window.y1);
 344
 345   float costVariance = 0.0f;
 346   int countVariance = 0;
 347   float costRoughness = 0.0f;
 348   int countRoughness = 0;
 349   float costLitFacing = 0.0f;
 350   int countLitFacing = 0;
 351   float costShading = 0.0f;
 352   int countShading = 0;
 353   float costTransitionElev = 0.0f;
 354   int countTransitionElev = 0;
 355   float costTransitionTangential = 0.0f;
 356   int countTransitionTangential = 0;
 357   float costTransitionCurvingDown = 0.0f;
 358   int countTransitionCurvingDown = 0;
 359   for (int j = y1; j <= y2; ++j)
 360   {
 361     for (int i = x1; i <= x2; ++i)
 362     {
 363       float* elevationPtr = pixelElevation(i, j);
 364       if (0 == elevationPtr)
 365       {
 366         continue;
 367       }
 368
 369       int index = j*width + i;
 370
 371       // Basic data retrieval
 372       float elevation = *elevationPtr;
 373       PerPixelCache* cache = &m_data->buffer()[index];
 374       bool isShadowed = (m_shadowData->buffer()[index] < 0.5f);
 375       bool isShadowEntrance = isShadowed &&
 376                               cache->shadowTransition[TransitionEntrance][0] == i &&
 377                               cache->shadowTransition[TransitionEntrance][1] == j;
 378       bool isShadowExit = isShadowed &&
 379                           cache->shadowTransition[TransitionExit][0] == i &&
 380                           cache->shadowTransition[TransitionExit][1] == j;
 381
 382       // Gradient
 383       // Derivitive of gradient
 384       float dh_dx = 0.0f;
 385       float dh_dy = 0.0f;
 386       float d2h_dx2 = 0.0f;
 387       float d2h_dy2 = 0.0f;
 388       if (i > 0 && i < width-1)
 389       {
 390         float hp1 = m_elevationData->buffer()[index+1];
 391         float hm1 = m_elevationData->buffer()[index-1];
 392         d2h_dx2 = (hp1+hm1)/2 - elevation;
 393         dh_dx = (hp1-hm1)/2 / cache->resolution[0];
 394       }
 395       if (j > 0 && j < height-1)
 396       {
 397         float hp1 = m_elevationData->buffer()[index+width];
 398         float hm1 = m_elevationData->buffer()[index-width];
 399         d2h_dy2 = (hp1+hm1)/2 - elevation;
 400         dh_dy = (hp1-hm1)/2 / cache->resolution[1];
 401       }
 402
 403       // Calculate measures relating to the light direction
 404       // Gradient
 405       float hwm1 = m_elevationData->sampleFloat(((float)i - cache->light[0])/(width-1),
 406                                                 ((float)j - cache->light[1])/(height-1));
 407       float hwp1 = m_elevationData->sampleFloat(((float)i + cache->light[0])/(width-1),
 408                                                 ((float)j + cache->light[1])/(height-1));
 409       /// @todo Ensure units are right (this is in terms of h metres, w pixels)
 410       float dh_dw = (hwm1 - hwp1)/2 / cache->resolution[0]; // 30 metre resolution?
 411       float d2h_dw2 = (hwm1+hwp1)/2 - elevation;
 412
 413       // Minimise variance from original value
 414       costVariance += cache->reliability * (elevation-cache->originalElevation)*(elevation-cache->originalElevation);
 415       ++countVariance;
 416
 417       // Minimise roughness
 418       costRoughness += (d2h_dx2*d2h_dx2 + d2h_dy2*d2h_dy2);
 419       ++countRoughness;
 420
 421       if (!isShadowed)
 422       {
 423         /*
 424         float facingness = -dh_dw-cache->light[0];
 425         if (facingness < 0.0f)
 426         {
 427           // Lit pixels must face light
 428           costLitFacing += facingness*facingness;
 429           ++countLitFacing;
 430         }
 431         */
 432         // Simple shape from shading
 433         float intensity = m_shadingData->buffer()[index];
 434         // intensity = cos(theta) = L.N
 435         maths::Vector<3,float> xTan(1.0f, 0.0f, dh_dx);
 436         maths::Vector<3,float> yTan(0.0f, 1.0f, dh_dy);
 437         costShading += MySqr(acos(qMin(1.0f, intensity)) - acos(qMin(1.0f, cache->light*maths::cross(xTan, yTan).normalized())));
 438         ++countShading;
 439       }
 440       if (isShadowEntrance)
 441       {
 442         if (cache->shadowTransition[TransitionExit][0] >= 0 &&
 443             cache->shadowTransition[TransitionExit][0] < width &&
 444             cache->shadowTransition[TransitionExit][1] >= 0 &&
 445             cache->shadowTransition[TransitionExit][1] < height)
 446         {
 447           int exitIndex = width*cache->shadowTransition[TransitionExit][1] + cache->shadowTransition[TransitionExit][0];
 448           float exitElevation = m_elevationData->buffer()[exitIndex];
 449           costTransitionElev += MySqr(exitElevation + cache->shadowTransitionDistance[TransitionExit]*cache->light[2] - elevation);
 450           ++countTransitionElev;
 451         }
 452         // gradient tangential to sun direction
 453         costTransitionTangential += MySqr(dh_dw - cache->light[2]);
 454         ++countTransitionTangential;
 455         // second derivitive should be negative at entrance
 456         costTransitionCurvingDown += MySqr(qMax(0.0f, d2h_dw2));
 457         ++countTransitionCurvingDown;
 458       }
 459       if (0 && isShadowExit)
 460       {
 461         if (cache->shadowTransition[TransitionEntrance][0] >= 0 &&
 462             cache->shadowTransition[TransitionEntrance][0] < width &&
 463             cache->shadowTransition[TransitionEntrance][1] >= 0 &&
 464             cache->shadowTransition[TransitionEntrance][1] < height)
 465         {
 466           int entranceIndex = width*cache->shadowTransition[TransitionEntrance][1] + cache->shadowTransition[TransitionEntrance][0];
 467           float entranceElevation = m_elevationData->buffer()[entranceIndex];
 468           costTransitionElev += MySqr(entranceElevation - cache->shadowTransitionDistance[TransitionEntrance]*cache->light[2] - elevation);
 469           ++countTransitionElev;
 470         }
 471       }
 472     }
 473   }
 474
 475   return m_weights.variance*costVariance//countVariance
 476        + m_weights.roughness*costRoughness//countRoughness
 477        + m_weights.litFacing*costLitFacing//countLitFacing
 478        + m_weights.shading*costShading//countShading
 479        + m_weights.transitionElev*costTransitionElev//countTransitionElev
 480        + m_weights.transitionTangential*costTransitionTangential//countTransitionTangential
 481        + m_weights.transitionCurvingDown*costTransitionCurvingDown//countTransitionCurvingDown
 482        ;
 483 }
 484
 485 /// Calculate the relative cost of changing a pixel's elevation.
 486 float tcElevationOptimization::costChange(int x, int y, float dh)
 487 {
 488   float* elev = pixelElevation(x, y);
 489   if (0 == elev)
 490   {
 491     return 0;
 492   }
 493
 494   // Assume a single change in elevation only affects cost of pixels to effectiveRange from it.
 495   // Now technically this doesn't hold for shadow edges as all pixels across the shadow can be affected.
 496   // If the pixel is a shadow edge, take into account the cost of extra pixels.
 497   static const int effectiveRange = 1;
 498   float currentCost = cost(x-effectiveRange, y-effectiveRange, x+effectiveRange, y+effectiveRange);
 499   float oldElevation = *elev;
 500   *elev += dh;
 501   float newCost = cost(x-effectiveRange, y-effectiveRange, x+effectiveRange, y+effectiveRange);
 502   *elev = oldElevation;
 503   return newCost - currentCost;
 504 }
 505
 506 /// Optimize a whole range of elevations.
 507 void tcElevationOptimization::optimizeElevations(int x1, int y1, int x2, int y2, float dh)
 508 {
 509   for (int j = y1; j <= y2; ++j)
 510   {
 511     for (int i = x1; i <= x2; ++i)
 512     {
 513       int index = j*m_data->width() + i;
 514       while (true)
 515       {
 516         // Find the cost of adjusting the elevation by dh in each direction
 517         float costInc = costChange(i,j, dh);
 518         float costDec = costChange(i,j,-dh);
 519         if (costInc < 0.0f && costInc < costDec)
 520         {
 521           adjustElevation(i,j,dh,true);
 522         }
 523         else if (costDec < 0.0f)
 524         {
 525           adjustElevation(i,j,-dh,true);
 526         }
 527         else
 528         {
 529           break;
 530         }
 531       }
 532     }
 533   }
 534 }
 535
 536 /// Get pointer to elevation at a pixel.
 537 float* tcElevationOptimization::pixelElevation(int x, int y) const
 538 {
 539   if (x < 0 || x >= m_elevationData->width() ||
 540       y < 0 || y >= m_elevationData->height())
 541   {
 542     return 0;
 543   }
 544   int index = y*m_elevationData->width() + x;
 545   return &(m_elevationData->buffer()[index]);
 546 }
 547
 548 /// Adjust elevation at a pixel.
 549 float tcElevationOptimization::adjustElevation(int x, int y, float dh, bool relative)
 550 {
 551   float* elev = pixelElevation(x,y);
 552   if (0 != elev)
 553   {
 554     float old = *elev;
 555     if (relative)
 556     {
 557       dh += old;
 558     }
 559     *elev = dh;
 560     return old;
 561   }
 562   else
 563   {
 564     return 0.0f;
 565   }
 566 }