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 #include "tcShadowClassifyingData.h"
  29 #include <SplineDefinitions.h>
  30
  31 #include <QObject>
  32 #include <QWidget>
  33 #include <QVBoxLayout>
  34 #include <QHBoxLayout>
  35 #include <QGridLayout>
  36 #include <QPushButton>
  37 #include <QSpinBox>
  38 #include <QLabel>
  39 #include <QSlider>
  40 #include <QMessageBox>
  41
  42 /*
  43  * Constructors + destructor
  44  */
  45
  46 /// Primary constructor.
  47 tcElevationOptimization::tcElevationOptimization(const QList<tcShadowClassifyingData*>& datasets, tcSrtmModel* dem, tcGeoImageData* imagery)
  48 : tcChannelDem(dem, imagery,
  49                tr("Elevation Optimization"),
  50                tr("Optimizes elevation using image data."),
  51                false)
  52 , m_numDatasets(datasets.size())
  53 , m_datasets(new tcGeoImageData*[m_numDatasets])
  54 , m_texToDatasetTex(new tcAffineTransform2<double>[m_numDatasets])
  55 , m_shadowChannels(new tcChannel*[m_numDatasets])
  56 , m_shadingChannels(new tcChannel*[m_numDatasets])
  57 , m_elevationData(0)
  58 , m_datum(new tcTypedPixelData<PerDatasetPixelCache>*[m_numDatasets])
  59 , m_data(0)
  60 , m_shadowData(new Reference<tcPixelData<float> >[m_numDatasets])
  61 , m_shadingData(new Reference<tcPixelData<float> >[m_numDatasets])
  62 , m_loadingFromDem(false)
  63 , m_configWidget(0)
  64 , m_spinIterations(0)
  65 {
  66   for (int i = 0; i < m_numDatasets; ++i)
  67   {
  68     tcShadowClassifyingData* landsat = datasets[i];
  69     m_datasets[i] = landsat;
  70     m_texToDatasetTex[i] = imagery->texToGeo() * landsat->geoToTex();
  71     m_shadowChannels[i] = landsat->shadowClassification();
  72     if (0 != m_shadowChannels[i])
  73     {
  74       m_shadowChannels[i]->addDerivitive(this);
  75     }
  76     m_shadingChannels[i] = landsat->shading();
  77     if (0 != m_shadingChannels[i])
  78     {
  79       m_shadingChannels[i]->addDerivitive(this);
  80     }
  81     m_datum[i] = 0;
  82   }
  83 }
  84
  85 /// Destructor.
  86 tcElevationOptimization::~tcElevationOptimization()
  87 {
  88   delete [] m_datasets;
  89   delete [] m_texToDatasetTex;
  90   for (int i = 0; i < m_numDatasets; ++i)
  91   {
  92     if (0 != m_shadowChannels[i])
  93     {
  94       m_shadowChannels[i]->removeDerivitive(this);
  95     }
  96     if (0 != m_shadingChannels[i])
  97     {
  98       m_shadingChannels[i]->removeDerivitive(this);
  99     }
 100     delete m_datum[i];
 101   }
 102   delete [] m_shadowChannels;
 103   delete [] m_shadingChannels;
 104   delete m_elevationData;
 105   delete [] m_datum;
 106   delete m_data;
 107   delete [] m_shadowData;
 108   delete [] m_shadingData;
 109   delete m_configWidget;
 110 }
 111
 112 /*
 113  * Main image interface
 114  */
 115
 116 tcChannelConfigWidget* tcElevationOptimization::configWidget()
 117 {
 118   if (0 == m_configWidget)
 119   {
 120     m_configWidget = new tcChannelConfigWidget();
 121     m_configWidget->setWindowTitle(tr("%1 Configuration").arg(name()));
 122
 123     QVBoxLayout* layout = new QVBoxLayout(m_configWidget);
 124
 125     // Row of optimisation controls
 126     QWidget* buttons1 = new QWidget(m_configWidget);
 127       layout->addWidget(buttons1);
 128       QHBoxLayout* buttons1Layout = new QHBoxLayout(buttons1);
 129
 130       // reset DEM (sets all corrections to the original (interpolated) values
 131       QPushButton* btnResetDem = new QPushButton(tr("Reset DEM"), m_configWidget);
 132         buttons1Layout->addWidget(btnResetDem);
 133         btnResetDem->setToolTip(tr("Reset DEM corrected values to interpolated values"));
 134         connect(btnResetDem, SIGNAL(clicked(bool)), this, SLOT(resetDem()));
 135
 136       // load base from DEM (loads DEM elevations into elevation buffer)
 137       QPushButton* btnLoadDem = new QPushButton(tr("Load DEM"), m_configWidget);
 138         buttons1Layout->addWidget(btnLoadDem);
 139         btnLoadDem->setToolTip(tr("Load corrected values from DEM into elevation field"));
 140         connect(btnLoadDem, SIGNAL(clicked(bool)), this, SLOT(loadFromDem()));
 141
 142       // number of iterations (select number of iterations to perform)
 143       m_spinIterations = new QSpinBox(m_configWidget);
 144         buttons1Layout->addWidget(m_spinIterations);
 145         m_spinIterations->setToolTip(tr("Number of iterations of optimization to perform"));
 146         m_spinIterations->setRange(1, 1000);
 147         m_spinIterations->setValue(10);
 148
 149       // optimize (iterates on elevation buffer and write back to DEM at intervals and end)
 150       QPushButton* btnOptimize = new QPushButton(tr("Optimize Iteratively"), m_configWidget);
 151         buttons1Layout->addWidget(btnOptimize);
 152         btnOptimize->setToolTip(tr("Iteratively optimize elevation field and write results back to DEM"));
 153         connect(btnOptimize, SIGNAL(clicked(bool)), this, SLOT(optimize()));
 154
 155     // Another row of optimisation controls
 156     QWidget* buttons2 = new QWidget(m_configWidget);
 157       layout->addWidget(buttons2);
 158       QHBoxLayout* buttons2Layout = new QHBoxLayout(buttons2);
 159
 160       // apply hard constraints, adjusting reliabilities
 161       QPushButton* btnHardConstrain = new QPushButton(tr("Hard Constrain"), m_configWidget);
 162         buttons2Layout->addWidget(btnHardConstrain);
 163         btnLoadDem->setToolTip(tr("Apply hard shadow constraints"));
 164         connect(btnHardConstrain, SIGNAL(clicked(bool)), this, SLOT(applyHardConstraints()));
 165
 166       // bilinear interpolation of voids
 167       QPushButton* btnBilinear = new QPushButton(tr("Bilinear"), m_configWidget);
 168         buttons2Layout->addWidget(btnBilinear);
 169         btnLoadDem->setToolTip(tr("Void-fill using bilinear interpolation"));
 170         connect(btnBilinear, SIGNAL(clicked(bool)), this, SLOT(voidFillBilinear()));
 171
 172       // bicubic interpolation of voids
 173       QPushButton* btnBicubic = new QPushButton(tr("Bicubic"), m_configWidget);
 174         buttons2Layout->addWidget(btnBicubic);
 175         btnLoadDem->setToolTip(tr("Void-fill using bicubic interpolation"));
 176         connect(btnBicubic, SIGNAL(clicked(bool)), this, SLOT(voidFillBicubic()));
 177
 178     // Rows of weight sliders
 179 #define NUM_WEIGHTS 9
 180     QString weightNames[NUM_WEIGHTS] = {
 181       tr("Variance"),
 182       tr("Roughness"),
 183       tr("Steepness"),
 184       tr("Lit facing sun"),
 185       tr("Shape from shading"),
 186       tr("Below shadow line"),
 187       tr("Transition elevation"),
 188       tr("Transition slope"),
 189       tr("Transition curving down")
 190     };
 191     QWidget* sliders = new QWidget(m_configWidget);
 192       layout->addWidget(sliders);
 193       QGridLayout* slidersLayout = new QGridLayout(sliders);
 194
 195       for (int i = 0; i < NUM_WEIGHTS; ++i)
 196       {
 197         QLabel* label = new QLabel(weightNames[i], sliders);
 198           slidersLayout->addWidget(label, i, 0);
 199
 200         QSlider* slider = new QSlider(Qt::Horizontal, sliders);
 201           slidersLayout->addWidget(slider, i, 1);
 202       }
 203   }
 204   return m_configWidget;
 205 }
 206
 207 /*
 208  * Slots
 209  */
 210
 211 /// Reset DEM.
 212 void tcElevationOptimization::resetDem()
 213 {
 214 }
 215
 216 /// Load elevation data from DEM.
 217 void tcElevationOptimization::loadFromDem()
 218 {
 219   delete m_elevationData;
 220   m_elevationData = 0;
 221
 222   delete m_data;
 223   m_data = 0;
 224
 225   for (int i = 0; i < m_numDatasets; ++i)
 226   {
 227     delete m_datum[i];
 228     m_datum[i] = 0;
 229   }
 230
 231   m_loadingFromDem = true;
 232   invalidate();
 233   m_configWidget->requestRedraw();
 234   m_loadingFromDem = false;
 235 }
 236
 237 /// Apply hard constraints to elevation data.
 238 void tcElevationOptimization::applyHardConstraints()
 239 {
 240   if (0 != m_elevationData)
 241   {
 242     const int width = m_elevationData->width();
 243     const int height = m_elevationData->height();
 244     // Now start optimising
 245     startProcessing("Applying hard constraints");
 246     for (int j = 0; j <= height; ++j)
 247     {
 248       progress((float)j/(height-1));
 249       for (int i = 0; i <= width; ++i)
 250       {
 251         int index = j*width + i;
 252
 253         // Only apply constraints to unreliable pixels
 254         PerPixelCache* cache = &m_data->buffer()[index];
 255         if (cache->reliability != 0)
 256         {
 257           continue;
 258         }
 259
 260         // Per dataset
 261         for (int ds = 0; ds < m_numDatasets; ++ds)
 262         {
 263           PerDatasetPixelCache* dsCache = &m_datum[ds]->buffer()[index];
 264
 265           bool isShadowed = (m_shadowData[ds]->sampleFloat((float)i/(width-1), (float)j/(height-1)) < 0.5f);
 266           bool isShadowEntrance = isShadowed &&
 267                                   dsCache->shadowTransition[TransitionEntrance][0] == i &&
 268                                   dsCache->shadowTransition[TransitionEntrance][1] == j;
 269           bool isShadowExit = isShadowed &&
 270                               dsCache->shadowTransition[TransitionExit][0] == i &&
 271                               dsCache->shadowTransition[TransitionExit][1] == j;
 272
 273           if (isShadowEntrance)
 274           {
 275             if (dsCache->shadowTransition[TransitionExit][0] >= 0 &&
 276                 dsCache->shadowTransition[TransitionExit][0] < width &&
 277                 dsCache->shadowTransition[TransitionExit][1] >= 0 &&
 278                 dsCache->shadowTransition[TransitionExit][1] < height)
 279             {
 280               int exitIndex = width*dsCache->shadowTransition[TransitionExit][1] + dsCache->shadowTransition[TransitionExit][0];
 281               PerPixelCache* exitCache = &m_data->buffer()[exitIndex];
 282               // Great, we can say for definite the elevation of this pixel because the corresponding pixel is reliable
 283               if (exitCache->reliability != 0)
 284               {
 285                 float exitElevation = m_elevationData->buffer()[exitIndex];
 286                 m_elevationData->buffer()[index] = exitElevation + dsCache->shadowTransitionDistance[TransitionExit]*dsCache->sinLightElevation;
 287                 cache->reliability = 1;
 288                 break;
 289               }
 290             }
 291           }
 292           else if (isShadowExit)
 293           {
 294             if (dsCache->shadowTransition[TransitionEntrance][0] >= 0 &&
 295                 dsCache->shadowTransition[TransitionEntrance][0] < width &&
 296                 dsCache->shadowTransition[TransitionEntrance][1] >= 0 &&
 297                 dsCache->shadowTransition[TransitionEntrance][1] < height)
 298             {
 299               int entranceIndex = width*dsCache->shadowTransition[TransitionEntrance][1] + dsCache->shadowTransition[TransitionEntrance][0];
 300               PerPixelCache* entranceCache = &m_data->buffer()[entranceIndex];
 301               // Great, we can say for definite the elevation of this pixel because the corresponding pixel is reliable
 302               if (entranceCache->reliability != 0)
 303               {
 304                 float entranceElevation = m_elevationData->buffer()[entranceIndex];
 305                 m_elevationData->buffer()[index] = entranceElevation - dsCache->shadowTransitionDistance[TransitionEntrance]*dsCache->sinLightElevation;
 306                 cache->reliability = 1;
 307                 break;
 308               }
 309             }
 310           }
 311         }
 312       }
 313     }
 314     endProcessing();
 315
 316     invalidate();
 317     m_configWidget->requestRedraw();
 318   }
 319   else
 320   {
 321     QMessageBox::warning(m_configWidget,
 322                          tr("Elevation field not initialised."),
 323                          tr("Please ensure the optimisation channel is displayed."));
 324   }
 325 }
 326
 327 /// Void-fill bilinear.
 328 void tcElevationOptimization::voidFillBilinear()
 329 {
 330   if (0 != m_elevationData)
 331   {
 332     const int width = m_elevationData->width();
 333     const int height = m_elevationData->height();
 334     // Now start optimising
 335     startProcessing("Bilinear void filling");
 336     for (int j = 0; j < height; ++j)
 337     {
 338       progress(0.5f*j/(height-1));
 339       int lastNonVoid = -1;
 340       float lastNonVoidValue;
 341       for (int i = 0; i < width; ++i)
 342       {
 343         int index = j*width + i;
 344         if (m_data->buffer()[index].reliability != 0)
 345         {
 346           float value = m_elevationData->buffer()[index];
 347           if (lastNonVoid != -1)
 348           {
 349             int gap = i - lastNonVoid;
 350             float startAlt = lastNonVoidValue;
 351             float dAlt = (value - lastNonVoidValue) / gap;
 352             // Lovely, between lastNonVoid and i is a void
 353             for (int ii = lastNonVoid+1; ii < i; ++ii)
 354             {
 355               int iIndex = j*width + ii;
 356               float alt = startAlt + dAlt * (ii - lastNonVoid);
 357               m_elevationData->buffer()[iIndex] = alt;
 358               // Keep track of the gap
 359               m_data->buffer()[iIndex].temporary.i = gap;
 360             }
 361           }
 362           else
 363           {
 364             m_data->buffer()[index].temporary.i = -1;
 365           }
 366           lastNonVoid = i;
 367           lastNonVoidValue = value;
 368         }
 369         else
 370         {
 371           m_data->buffer()[index].temporary.i = -1;
 372         }
 373       }
 374     }
 375     for (int i = 0; i < width; ++i)
 376     {
 377       progress(0.5f + 0.5f*i/(width-1));
 378       int lastNonVoid = -1;
 379       float lastNonVoidValue;
 380       for (int j = 0; j < height; ++j)
 381       {
 382         int index = j*width + i;
 383         float value = m_elevationData->buffer()[index];
 384         if (m_data->buffer()[index].reliability != 0)
 385         {
 386           if (lastNonVoid != -1)
 387           {
 388             int gap = j - lastNonVoid;
 389             float startAlt = lastNonVoidValue;
 390             float dAlt = (value - lastNonVoidValue) / gap;
 391             // Lovely, between lastNonVoid and j is a void
 392             for (int jj = lastNonVoid+1; jj < j; ++jj)
 393             {
 394               // Average with previous value if non zero, otherwise overwrite
 395               int iIndex = jj*width + i;
 396               float alt = startAlt + dAlt * (jj - lastNonVoid);
 397               int otherGap = m_data->buffer()[iIndex].temporary.i;
 398               if (otherGap != -1)
 399               {
 400                 float prevAlt = m_elevationData->buffer()[iIndex];
 401                 // Prefer the value of the smaller gap
 402                 alt = (alt * otherGap + prevAlt * gap) / (gap+otherGap);
 403               }
 404               m_elevationData->buffer()[iIndex] = alt;
 405             }
 406           }
 407           lastNonVoid = j;
 408           lastNonVoidValue = value;
 409         }
 410       }
 411     }
 412     endProcessing();
 413
 414     invalidate();
 415     m_configWidget->requestRedraw();
 416   }
 417   else
 418   {
 419     QMessageBox::warning(m_configWidget,
 420                          tr("Elevation field not initialised."),
 421                          tr("Please ensure the optimisation channel is displayed."));
 422   }
 423 }
 424
 425 /// Void-fill bicubic.
 426 void tcElevationOptimization::voidFillBicubic()
 427 {
 428   if (0 != m_elevationData)
 429   {
 430     const int width = m_elevationData->width();
 431     const int height = m_elevationData->height();
 432     // Now start optimising
 433     startProcessing("Bicubic void filling");
 434     for (int j = 0; j < height; ++j)
 435     {
 436       progress(0.5f*j/(height-1));
 437       int lastNonVoid = -1;
 438       float lastNonVoidValue;
 439       for (int i = 0; i < width; ++i)
 440       {
 441         int index = j*width + i;
 442         if (m_data->buffer()[index].reliability != 0)
 443         {
 444           float value = m_elevationData->buffer()[index];
 445           if (lastNonVoid != -1)
 446           {
 447             int gap = i - lastNonVoid;
 448             // Obtain some control points for the splines
 449             // use the two samples either side of the gap as control points 1 and 2
 450             // use the ones before and after these samples, extended to the size of
 451             // the gap as control points 0 and 3
 452             float startAlt = lastNonVoidValue;
 453             float beforeStartAlt = startAlt;
 454             if (lastNonVoid > 0)
 455             {
 456               if (0 != m_data->buffer()[j*width + lastNonVoid - 1].reliability)
 457               {
 458                 beforeStartAlt = startAlt + (m_elevationData->buffer()[j*width + lastNonVoid - 1] - startAlt) * gap;
 459               }
 460             }
 461             float endAlt = value;
 462             float afterEndAlt = endAlt;
 463             if (i < width-1)
 464             {
 465               if (0 != m_data->buffer()[j*width + i + 1].reliability)
 466               {
 467                 afterEndAlt = endAlt + (m_elevationData->buffer()[j*width + i + 1] - endAlt) * gap;
 468               }
 469             }
 470             // Lovely, between lastNonVoid and i is a void
 471             for (int ii = lastNonVoid+1; ii < i; ++ii)
 472             {
 473               int iIndex = j*width + ii;
 474               float u = (float)(ii-lastNonVoid) / gap;
 475               float u2 = u*u;
 476               float u3 = u2*u;
 477               float alt = maths::catmullRomSpline.Spline(u, u2, u3, beforeStartAlt, startAlt, endAlt, afterEndAlt);
 478               m_elevationData->buffer()[iIndex] = alt;
 479               // Keep track of the gap
 480               m_data->buffer()[iIndex].temporary.i = gap;
 481             }
 482           }
 483           else
 484           {
 485             m_data->buffer()[index].temporary.i = -1;
 486           }
 487           lastNonVoid = i;
 488           lastNonVoidValue = value;
 489         }
 490         else
 491         {
 492           m_data->buffer()[index].temporary.i = -1;
 493         }
 494       }
 495     }
 496     for (int i = 0; i < width; ++i)
 497     {
 498       progress(0.5f + 0.5f*i/(width-1));
 499       int lastNonVoid = -1;
 500       float lastNonVoidValue;
 501       for (int j = 0; j < height; ++j)
 502       {
 503         int index = j*width + i;
 504         float value = m_elevationData->buffer()[index];
 505         if (m_data->buffer()[index].reliability != 0)
 506         {
 507           if (lastNonVoid != -1)
 508           {
 509             int gap = j - lastNonVoid;
 510             // Obtain some control points for the splines
 511             // use the two samples either side of the gap as control points 1 and 2
 512             // use the ones before and after these samples, extended to the size of
 513             // the gap as control points 0 and 3
 514             float startAlt = lastNonVoidValue;
 515             float beforeStartAlt = startAlt;
 516             if (lastNonVoid > 0)
 517             {
 518               if (0 != m_data->buffer()[(lastNonVoid-1)*width + i].reliability)
 519               {
 520                 beforeStartAlt = startAlt + (m_elevationData->buffer()[(lastNonVoid-1)*width + i] - startAlt) * gap;
 521               }
 522             }
 523             float endAlt = value;
 524             float afterEndAlt = endAlt;
 525             if (j < height-1)
 526             {
 527               if (0 != m_data->buffer()[(j+1)*width + i].reliability)
 528               {
 529                 afterEndAlt = endAlt + (m_elevationData->buffer()[(j+1)*width + i] - endAlt) * gap;
 530               }
 531             }
 532             // Lovely, between lastNonVoid and j is a void
 533             for (int jj = lastNonVoid+1; jj < j; ++jj)
 534             {
 535               // Average with previous value if non zero, otherwise overwrite
 536               int iIndex = jj*width + i;
 537               float u = (float)(jj-lastNonVoid) / gap;
 538               float u2 = u*u;
 539               float u3 = u2*u;
 540               float alt = maths::catmullRomSpline.Spline(u, u2, u3, beforeStartAlt, startAlt, endAlt, afterEndAlt);
 541               int otherGap = m_data->buffer()[iIndex].temporary.i;
 542               if (otherGap != -1)
 543               {
 544                 float prevAlt = m_elevationData->buffer()[iIndex];
 545                 // Prefer the value of the smaller gap
 546                 alt = (alt * otherGap + prevAlt * gap) / (gap+otherGap);
 547               }
 548               m_elevationData->buffer()[iIndex] = alt;
 549             }
 550           }
 551           lastNonVoid = j;
 552           lastNonVoidValue = value;
 553         }
 554       }
 555     }
 556     endProcessing();
 557
 558     invalidate();
 559     m_configWidget->requestRedraw();
 560   }
 561   else
 562   {
 563     QMessageBox::warning(m_configWidget,
 564                          tr("Elevation field not initialised."),
 565                          tr("Please ensure the optimisation channel is displayed."));
 566   }
 567 }
 568
 569 /// Optimize a number of times.
 570 void tcElevationOptimization::optimize()
 571 {
 572   if (0 != m_elevationData)
 573   {
 574     const int width = m_elevationData->width();
 575     const int height = m_elevationData->height();
 576     // Now start optimising
 577     startProcessing("Optimizing elevation data");
 578     int passes = m_spinIterations->value();
 579     for (int i = 0; i < passes; ++i)
 580     {
 581       progress((float)i/passes);
 582       m_weights.variance = 1.0f;
 583       m_weights.roughness = 50.0f *(1.0f - (float)i/passes);
 584       m_weights.steepness = 0;//5000;//200.0f;
 585       m_weights.litFacing = 1000.0f;
 586       m_weights.shading = 1.0f * i/passes;
 587       m_weights.belowShadowLine = 2000.0f*i/passes;
 588       m_weights.transitionElev = 1000.0f;
 589       m_weights.transitionTangential = 1000.0f*i/passes;
 590       m_weights.transitionCurvingDown = 10.0f*i/passes;
 591       optimizeElevations(0, 0, width-1, height-1, 8.0f - 7.8f*i/passes, (int)(4.5f - 4.0f*i/passes));
 592
 593       if ((i+1) % 2 == 0 || i == passes-1)
 594       {
 595         invalidate();
 596         m_configWidget->requestRedraw();
 597
 598         startProcessing("Optimizing elevation data");
 599         progress((float)i/passes);
 600       }
 601     }
 602     endProcessing();
 603   }
 604   else
 605   {
 606     QMessageBox::warning(m_configWidget,
 607                          tr("Elevation field not initialised."),
 608                          tr("Please ensure the optimisation channel is displayed."));
 609   }
 610 }
 611
 612 /*
 613  * Interface for derived class to implement
 614  */
 615
 616 void tcElevationOptimization::roundPortion(double* x1, double* y1, double* x2, double* y2)
 617 {
 618   //m_shadowChannel->roundPortion(x1,y1,x2,y2);
 619 }
 620
 621 tcAbstractPixelData* tcElevationOptimization::loadPortion(double x1, double y1, double x2, double y2)
 622 {
 623   int width = -1;
 624   int height = -1;
 625
 626   for (int i = 0; i < m_numDatasets; ++i)
 627   {
 628     maths::Vector<2,double> xy1 = m_texToDatasetTex[i] * maths::Vector<2,double>(x1,y1);
 629     maths::Vector<2,double> xy2 = m_texToDatasetTex[i] * maths::Vector<2,double>(x2,y2);
 630     Reference<tcAbstractPixelData> channelData = m_shadowChannels[i]->portion(xy1[0], xy1[1], xy2[0], xy2[1]);
 631     if (i == 0)
 632     {
 633       width = channelData->width();
 634       height = channelData->height();
 635     }
 636     m_shadowData[i] = dynamicCast<tcPixelData<float>*>(channelData);
 637     if (0 == m_shadowData[i])
 638     {
 639       return 0;
 640     }
 641     if (0 != m_shadingChannels[i])
 642     {
 643       Reference<tcAbstractPixelData> shadingData = m_shadingChannels[i]->portion(xy1[0], xy1[1], xy2[0], xy2[1]);
 644       m_shadingData[i] = dynamicCast<tcPixelData<float>*>(shadingData);
 645     }
 646     else
 647     {
 648       m_shadingData[i] = 0;
 649     }
 650   }
 651
 652   /// @todo What if its a different size because the portion has changed?
 653   if (0 == m_elevationData || 0 == m_data)
 654   {
 655     delete m_elevationData;
 656     for (int i = 0; i < m_numDatasets; ++i)
 657     {
 658       delete m_datum[i];
 659     }
 660     delete m_data;
 661
 662     // Create a new pixel buffer, initialised to altitude
 663
 664     m_elevationData = new tcElevationData(width, height);
 665     // Find transform of elevation data
 666     tcGeo geoBase = geoAt(maths::Vector<2,float>(x1, y1));
 667     tcGeo geoDiagonal = geoAt(maths::Vector<2,float>(x1 + (x2-x1)/(width-1), y1 + (y2-y1)/(height-1))) - geoBase;
 668     m_elevationData->setGeoToPixels(tcAffineTransform2<double>(geoBase, geoDiagonal).inverse());
 669
 670     m_data = new tcTypedPixelData<PerPixelCache>(width, height);
 671     for (int i = 0; i < m_numDatasets; ++i)
 672     {
 673       m_datum[i] = new tcTypedPixelData<PerDatasetPixelCache>(width, height);
 674     }
 675
 676     int maxWidthHeight = qMax(width, height);
 677     double minXY12 = qMin(x2-x1, y2-y1);
 678
 679     startProcessing("Caching elevation characteristics");
 680     for (int j = 0; j < height; ++j)
 681     {
 682       progress((float)j/(height-1));
 683       for (int i = 0; i < width; ++i)
 684       {
 685         int index = j*width + i;
 686         // Transform coordinates
 687         maths::Vector<2,float> coord      (x1 + (x2-x1)*i / (width  - 1),
 688                                            y1 + (y2-y1)*j / (height - 1));
 689         tcGeo                  geoCoord = geoAt(coord);
 690
 691         PerPixelCache* cache = &m_data->buffer()[index];
 692
 693         // Find resolution
 694         maths::Vector<2,float> coordx      (x1 + (x2-x1)*(i+1) / (width  - 1),
 695                                             y1 + (y2-y1)*j     / (height - 1));
 696         maths::Vector<2,float> coordy      (x1 + (x2-x1)*i     / (width  - 1),
 697                                             y1 + (y2-y1)*(j+1) / (height - 1));
 698         tcGeo                  geoCoordx = geoAt(coordx) - geoCoord;
 699         tcGeo                  geoCoordy = geoAt(coordy) - geoCoord;
 700         // Find approx
 701         float res = geoCoordx.angle();
 702         cache->resolution[0] = res;
 703         cache->resolution[1] = geoCoordy.angle();
 704         cache->resolution *= 6378.137e3;
 705
 706         // Get some elevation data
 707         bool accurate;
 708         float altitude        = altitudeAt(geoCoord, true, &accurate);
 709         m_elevationData->buffer()[index] = altitude;
 710         cache->originalElevation = altitude;
 711         cache->reliability = accurate ? 1 : 0;
 712
 713         // Per dataset cache
 714         for (int ds = 0; ds < m_numDatasets; ++ds)
 715         {
 716           tcGeoImageData* imagery = m_datasets[ds];
 717           maths::Vector<3,float> light = lightDirectionAt(geoCoord, imagery);
 718
 719           PerDatasetPixelCache* dsCache = &m_datum[ds]->buffer()[index];
 720           dsCache->light = light;
 721           dsCache->sinLightElevation = light[2]/light.slice<0,2>().mag();
 722
 723           tcGeo lightScaled(light[0]*res/sin(geoCoord.lat()),
 724                             light[1]*res);
 725           tcGeo offset(geoCoord + lightScaled);
 726           for (int depthDirection = 0; depthDirection < 2; ++depthDirection)
 727           {
 728             tcGeo pos = geoCoord;
 729             tcGeo delta;
 730             if (depthDirection == TransitionEntrance)
 731             {
 732               delta = offset - geoCoord;
 733             }
 734             else
 735             {
 736               delta = geoCoord - offset;
 737             }
 738             pos = pos + delta;
 739             maths::Vector<2,float> tex = coord;
 740             dsCache->shadowTransition[depthDirection][0] = i;
 741             dsCache->shadowTransition[depthDirection][1] = j;
 742             while (tex[0] >= x1 && tex[0] <= x2 && tex[1] >= y2 && tex[1] <= y1)
 743             {
 744               if (m_shadowData[ds]->sampleFloat((tex[0]-x1)/(x2-x1), (tex[1]-y1)/(y2-y1)) > 0.5f)
 745               {
 746                 break;
 747               }
 748               dsCache->shadowTransition[depthDirection][0] = 0.5f + (tex[0]-x1)/(x2-x1)*(width -1);
 749               dsCache->shadowTransition[depthDirection][1] = 0.5f + (tex[1]-y1)/(y2-y1)*(height-1);
 750               pos = pos + delta;
 751               tex = textureAt(pos);
 752             }
 753             dsCache->shadowTransitionDistance[depthDirection] = (pos-geoCoord).angle()*6378.137e3;
 754           }
 755         }
 756       }
 757     }
 758   }
 759
 760   if (!m_loadingFromDem)
 761   {
 762     // Update elevation model
 763     startProcessing("Updating elevation model");
 764     dem()->updateFromElevationData(m_elevationData);
 765   }
 766
 767   // Downscale the elevation for viewing
 768   startProcessing("Downscaling for viewing");
 769   tcPixelData<GLubyte>* result = new tcPixelData<GLubyte>(width, height);
 770   int index = 0;
 771   for (int j = 0; j < height; ++j)
 772   {
 773     progress((float)j/(height-1));
 774     for (int i = 0; i < width; ++i)
 775     {
 776       result->buffer()[index] = 255.0f * m_elevationData->buffer()[index] / 4000.0f;
 777       ++index;
 778     }
 779   }
 780
 781   endProcessing();
 782
 783   return result;
 784 }
 785
 786 template <typename T>
 787 T MySqr(T t)
 788 {
 789   return t*t;
 790 }
 791
 792 /*
 793  * Private functions
 794  */
 795
 796 /// Calculate the cost for a set of pixels.
 797 float tcElevationOptimization::cost(int x1, int y1, int x2, int y2) const
 798 {
 799   int width = m_elevationData->width();
 800   int height = m_elevationData->height();
 801   if (x1 < 0)
 802   {
 803     x1 = 0;
 804   }
 805   if (x2 >= width)
 806   {
 807     x2 = width-1;
 808   }
 809   if (y1 < 0)
 810   {
 811     y1 = 0;
 812   }
 813   if (y2 >= height)
 814   {
 815     y2 = height-1;
 816   }
 817
 818   int maxWidthHeight = qMax(width, height);
 819   //double minXY12 = qMin(m_window.x2-m_window.x1, m_window.y2-m_window.y1);
 820
 821   float costVariance = 0.0f;
 822   int countVariance = 0;
 823   float costRoughness = 0.0f;
 824   int countRoughness = 0;
 825   float costSteepness = 0.0f;
 826   int countSteepness = 0;
 827   float costLitFacing = 0.0f;
 828   int countLitFacing = 0;
 829   float costShading = 0.0f;
 830   int countShading = 0;
 831   float costBelowShadowLine = 0.0f;
 832   int countBelowShadowLine = 0;
 833   float costTransitionElev = 0.0f;
 834   int countTransitionElev = 0;
 835   float costTransitionTangential = 0.0f;
 836   int countTransitionTangential = 0;
 837   float costTransitionCurvingDown = 0.0f;
 838   int countTransitionCurvingDown = 0;
 839   for (int j = y1; j <= y2; ++j)
 840   {
 841     for (int i = x1; i <= x2; ++i)
 842     {
 843       float* elevationPtr = pixelElevation(i, j);
 844       if (0 == elevationPtr)
 845       {
 846         continue;
 847       }
 848
 849       int index = j*width + i;
 850
 851       // Basic data retrieval
 852       float elevation = *elevationPtr;
 853       PerPixelCache* cache = &m_data->buffer()[index];
 854
 855       // Gradient
 856       // Derivitive of gradient
 857       float dh_dx = 0.0f;
 858       float dh_dy = 0.0f;
 859       float d2h_dx2 = 0.0f;
 860       float d2h_dy2 = 0.0f;
 861       if (i > 0 && i < width-1)
 862       {
 863         float hp1 = m_elevationData->buffer()[index+1];
 864         float hm1 = m_elevationData->buffer()[index-1];
 865         d2h_dx2 = (hp1+hm1)/2 - elevation;
 866         dh_dx = (hp1-hm1)/2 / cache->resolution[0];
 867       }
 868       if (j > 0 && j < height-1)
 869       {
 870         float hp1 = m_elevationData->buffer()[index+width];
 871         float hm1 = m_elevationData->buffer()[index-width];
 872         d2h_dy2 = (hp1+hm1)/2 - elevation;
 873         dh_dy = (hp1-hm1)/2 / cache->resolution[1];
 874       }
 875
 876       // Minimise variance from original value
 877       costVariance += cache->reliability * (elevation-cache->originalElevation)*(elevation-cache->originalElevation);
 878       ++countVariance;
 879
 880       // Minimise roughness
 881       costRoughness += (d2h_dx2*d2h_dx2 + d2h_dy2*d2h_dy2);
 882       ++countRoughness;
 883
 884       // Also minimise slope
 885       // This works really well at getting a nice slope
 886       costSteepness += (dh_dx*dh_dx + dh_dy*dh_dy);
 887       ++countSteepness;
 888
 889       // Per dataset
 890       for (int ds = 0; ds < m_numDatasets; ++ds)
 891       {
 892         PerDatasetPixelCache* dsCache = &m_datum[ds]->buffer()[index];
 893
 894         bool isShadowed = (m_shadowData[ds]->sampleFloat((float)i/(width-1), (float)j/(height-1)) < 0.5f);
 895         bool isShadowEntrance = isShadowed &&
 896                                 dsCache->shadowTransition[TransitionEntrance][0] == i &&
 897                                 dsCache->shadowTransition[TransitionEntrance][1] == j;
 898         bool isShadowExit = isShadowed &&
 899                             dsCache->shadowTransition[TransitionExit][0] == i &&
 900                             dsCache->shadowTransition[TransitionExit][1] == j;
 901
 902         // Calculate measures relating to the light direction
 903         // Gradient
 904         float hwm1 = m_elevationData->sampleFloat(((float)i - dsCache->light[0])/(width-1),
 905             ((float)j - dsCache->light[1])/(height-1));
 906         float hwp1 = m_elevationData->sampleFloat(((float)i + dsCache->light[0])/(width-1),
 907             ((float)j + dsCache->light[1])/(height-1));
 908         /// @todo Ensure units are right (this is in terms of h metres, w pixels)
 909         float dh_dw = (hwm1 - hwp1)/2 / cache->resolution[0]; // 30 metre resolution?
 910         float d2h_dw2 = (hwm1+hwp1)/2 - elevation;
 911
 912         if (isShadowed)
 913         {
 914           // maximum elevation is between elevation at entrance and elevation at exit
 915           if (dsCache->shadowTransition[TransitionExit][0] >= 0 &&
 916               dsCache->shadowTransition[TransitionExit][0] < width &&
 917               dsCache->shadowTransition[TransitionExit][1] >= 0 &&
 918               dsCache->shadowTransition[TransitionExit][1] < height &&
 919               dsCache->shadowTransition[TransitionEntrance][0] >= 0 &&
 920               dsCache->shadowTransition[TransitionEntrance][0] < width &&
 921               dsCache->shadowTransition[TransitionEntrance][1] >= 0 &&
 922               dsCache->shadowTransition[TransitionEntrance][1] < height)
 923           {
 924             int exitIndex = width*dsCache->shadowTransition[TransitionExit][1] + dsCache->shadowTransition[TransitionExit][0];
 925             int entranceIndex = width*dsCache->shadowTransition[TransitionEntrance][1] + dsCache->shadowTransition[TransitionEntrance][0];
 926             float exitElevation = m_elevationData->buffer()[exitIndex];
 927             float entranceElevation = m_elevationData->buffer()[entranceIndex];
 928
 929             costBelowShadowLine += MySqr(qMax(0.0f, elevation - (exitElevation     * dsCache->shadowTransitionDistance[TransitionEntrance]
 930                                                                 +entranceElevation * dsCache->shadowTransitionDistance[TransitionExit])
 931                                                               / (dsCache->shadowTransitionDistance[TransitionEntrance]
 932                                                                 +dsCache->shadowTransitionDistance[TransitionExit])));
 933             ++countBelowShadowLine;
 934           }
 935         }
 936         else
 937         {
 938           float facingness = -dh_dw-dsCache->light[0];
 939           if (facingness < 0.0f)
 940           {
 941             // Lit pixels must face light
 942             costLitFacing += facingness*facingness;
 943             ++countLitFacing;
 944           }
 945           // Simple shape from shading
 946           if (0 != m_shadingData[ds])
 947           {
 948             float intensity = m_shadingData[ds]->sampleFloat((float)i/(width-1), (float)j/(height-1));
 949             // intensity = cos(theta) = L.N
 950             maths::Vector<3,float> xTan(1.0f, 0.0f, dh_dx);
 951             maths::Vector<3,float> yTan(0.0f, 1.0f, dh_dy);
 952             costShading += MySqr(acos(qMin(1.0f, intensity)) - acos(qMin(1.0f, dsCache->light*maths::cross(xTan, yTan).normalized())));
 953             ++countShading;
 954           }
 955         }
 956         if (isShadowEntrance)
 957         {
 958           if (dsCache->shadowTransition[TransitionExit][0] >= 0 &&
 959               dsCache->shadowTransition[TransitionExit][0] < width &&
 960               dsCache->shadowTransition[TransitionExit][1] >= 0 &&
 961               dsCache->shadowTransition[TransitionExit][1] < height)
 962           {
 963             int exitIndex = width*dsCache->shadowTransition[TransitionExit][1] + dsCache->shadowTransition[TransitionExit][0];
 964             float exitElevation = m_elevationData->buffer()[exitIndex];
 965             costTransitionElev += (1.0f - cache->reliability) * MySqr(exitElevation + dsCache->shadowTransitionDistance[TransitionExit]*dsCache->sinLightElevation - elevation);
 966             ++countTransitionElev;
 967           }
 968           // gradient tangential to sun direction
 969           costTransitionTangential += MySqr(dh_dw - dsCache->sinLightElevation);
 970           ++countTransitionTangential;
 971           // second derivitive should be negative at entrance
 972           costTransitionCurvingDown += MySqr(qMax(0.0f, d2h_dw2));
 973           ++countTransitionCurvingDown;
 974         }
 975         if (0 && isShadowExit)
 976         {
 977           if (dsCache->shadowTransition[TransitionEntrance][0] >= 0 &&
 978               dsCache->shadowTransition[TransitionEntrance][0] < width &&
 979               dsCache->shadowTransition[TransitionEntrance][1] >= 0 &&
 980               dsCache->shadowTransition[TransitionEntrance][1] < height)
 981           {
 982             int entranceIndex = width*dsCache->shadowTransition[TransitionEntrance][1] + dsCache->shadowTransition[TransitionEntrance][0];
 983             float entranceElevation = m_elevationData->buffer()[entranceIndex];
 984             // this has less effect, its usually "more" right
 985             /// @todo Smooth high reliability into low reliability
 986             costTransitionElev += 0.1f * (1.0f - cache->reliability) * MySqr(entranceElevation - dsCache->shadowTransitionDistance[TransitionEntrance]*dsCache->sinLightElevation - elevation);
 987             ++countTransitionElev;
 988           }
 989         }
 990       }
 991     }
 992   }
 993
 994   return m_weights.variance*costVariance//countVariance
 995        + m_weights.roughness*costRoughness//countRoughness
 996        + m_weights.steepness*costSteepness//countSteepness
 997        + m_weights.litFacing*costLitFacing//countLitFacing
 998        + m_weights.shading*costShading//countShading
 999        + m_weights.belowShadowLine*costBelowShadowLine//countBelowShadowLine
1000        + m_weights.transitionElev*costTransitionElev//countTransitionElev
1001        + m_weights.transitionTangential*costTransitionTangential//countTransitionTangential
1002        + m_weights.transitionCurvingDown*costTransitionCurvingDown//countTransitionCurvingDown
1003        ;
1004 }
1005
1006 /*
1007  * f(0) = 1
1008  * f(range) = small
1009  * f(x) = exp[ - (2*x/range)^2 ]
1010  */
1011 inline float elevationFalloffFunction(float xsq)
1012 {
1013   return qMax(0.0f, 1.0f - xsq);
1014   //return exp(-4*xsq);
1015 }
1016 /*
1017  */
1018
1019 /// Calculate the relative cost of changing a pixel's elevation.
1020 float tcElevationOptimization::costChange(int x, int y, float dh, int range)
1021 {
1022   // Override range for now
1023 #define ELEV_RANGE 4
1024   range = qMin(range, ELEV_RANGE);
1025   static const int effectiveRange = 1 + range;
1026   float oldElevations[1+2*ELEV_RANGE][1+2*ELEV_RANGE];
1027
1028   // Assume a single change in elevation only affects cost of pixels to effectiveRange from it.
1029   // Now technically this doesn't hold for shadow edges as all pixels across the shadow can be affected.
1030   // If the pixel is a shadow edge, take into account the cost of extra pixels.
1031   float currentCost = cost(x-effectiveRange, y-effectiveRange, x+effectiveRange, y+effectiveRange);
1032   for (int i = 0; i < 1+range*2; ++i)
1033   {
1034     int di = i - range;
1035     for (int j = 0; j < 1+range*2; ++j)
1036     {
1037       int dj = j - range;
1038       float* elev = pixelElevation(x+di, y+dj);
1039       if (0 != elev)
1040       {
1041         oldElevations[i][j] = *elev;
1042         *elev += dh*elevationFalloffFunction(float(di*di + dj*dj) / ((range-1)*(range-1)));
1043       }
1044     }
1045   }
1046   float newCost = cost(x-effectiveRange, y-effectiveRange, x+effectiveRange, y+effectiveRange);
1047   for (int i = 0; i < 1+range*2; ++i)
1048   {
1049     int di = i - range;
1050     for (int j = 0; j < 1+range*2; ++j)
1051     {
1052       int dj = j - range;
1053       float* elev = pixelElevation(x+di, y+dj);
1054       if (0 != elev)
1055       {
1056         *elev = oldElevations[i][j];
1057       }
1058     }
1059   }
1060   return newCost - currentCost;
1061 }
1062
1063 /// Optimize a whole range of elevations.
1064 void tcElevationOptimization::optimizeElevations(int x1, int y1, int x2, int y2, float dh, int range)
1065 {
1066   for (int j = y1; j <= y2; ++j)
1067   {
1068     for (int i = x1; i <= x2; ++i)
1069     {
1070       int index = j*m_data->width() + i;
1071       //while (true)
1072       {
1073         // Find the cost of adjusting the elevation by dh in each direction
1074         float costInc = costChange(i,j, dh, range);
1075         float costDec = costChange(i,j,-dh, range);
1076         if (costInc < 0.0f && costInc < costDec)
1077         {
1078           adjustElevation(i,j,dh, range);
1079         }
1080         else if (costDec < 0.0f)
1081         {
1082           adjustElevation(i,j,-dh, range);
1083         }
1084         else
1085         {
1086           //break;
1087         }
1088       }
1089     }
1090   }
1091 }
1092
1093 /// Get pointer to elevation at a pixel.
1094 float* tcElevationOptimization::pixelElevation(int x, int y) const
1095 {
1096   if (x < 0 || x >= m_elevationData->width() ||
1097       y < 0 || y >= m_elevationData->height())
1098   {
1099     return 0;
1100   }
1101   int index = y*m_elevationData->width() + x;
1102   return &(m_elevationData->buffer()[index]);
1103 }
1104
1105 /// Adjust elevation at a pixel.
1106 void tcElevationOptimization::adjustElevation(int x, int y, float dh, int range)
1107 {
1108   range = qMin(range, ELEV_RANGE);
1109   for (int i = 0; i < 1+range*2; ++i)
1110   {
1111     int di = i - range;
1112     for (int j = 0; j < 1+range*2; ++j)
1113     {
1114       int dj = j - range;
1115       float* elev = pixelElevation(x+di, y+dj);
1116       if (0 != elev)
1117       {
1118         *elev += dh*elevationFalloffFunction(float(di*di + dj*dj) / ((range+1)*(range+1)));
1119       }
1120     }
1121   }
1122 }