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 <qwt_plot.h>
  32 #include <qwt_plot_curve.h>
  33
  34 #include <QObject>
  35 #include <QWidget>
  36 #include <QVBoxLayout>
  37 #include <QHBoxLayout>
  38 #include <QGridLayout>
  39 #include <QPushButton>
  40 #include <QSpinBox>
  41 #include <QLabel>
  42 #include <QSlider>
  43 #include <QMessageBox>
  44 #include <QFileDialog>
  45 #include <QTextStream>
  46 #include <QThread>
  47
  48 /*
  49  * Private types
  50  */
  51
  52 /// Worker thread class.
  53 class tcElevationOptimization::WorkerThread : public QThread
  54 {
  55   public:
  56     /*
  57      * Constructors + destructor
  58      */
  59
  60     /// Constructor.
  61     WorkerThread(tcElevationOptimization* self)
  62     : QThread()
  63     , m_this(self)
  64     {
  65     }
  66
  67     /// Destructor.
  68     virtual ~WorkerThread()
  69     {
  70     }
  71
  72   protected:
  73     /// Reimplemented.
  74     virtual void run()
  75     {
  76       m_this->optimize();
  77
  78       m_this->m_stopProcessing = false;
  79       m_this->m_btnOptimize->setText(tr("Optimize Iteratively"));
  80       m_this->m_btnOptimize->update();
  81       m_this->m_processingMutex.unlock();
  82     }
  83
  84   private:
  85
  86     /*
  87      * Variables
  88      */
  89
  90     /// Main optimization object.
  91     tcElevationOptimization* m_this;
  92 };
  93
  94 /*
  95  * Constructors + destructor
  96  */
  97
  98 /// Primary constructor.
  99 tcElevationOptimization::tcElevationOptimization(const QList<tcShadowClassifyingData*>& datasets, tcSrtmModel* dem, tcGeoImageData* imagery)
 100 : tcChannelDem(dem, imagery,
 101                tr("Elevation Optimization"),
 102                tr("Optimizes elevation using image data."),
 103                false, true)
 104 , m_numDatasets(datasets.size())
 105 , m_datasets(new tcGeoImageData*[m_numDatasets])
 106 , m_texToDatasetTex(new tcAffineTransform2<double>[m_numDatasets])
 107 , m_shadowChannels(new tcChannel*[m_numDatasets])
 108 , m_shadingChannels(new tcChannel*[m_numDatasets])
 109 , m_elevationData(0)
 110 , m_datum(new tcTypedPixelData<PerDatasetPixelCache>*[m_numDatasets])
 111 , m_data(0)
 112 , m_shadowData(new Reference<tcPixelData<float> >[m_numDatasets])
 113 , m_shadingData(new Reference<tcPixelData<float> >[m_numDatasets])
 114 , m_voidsOnly(false)
 115 , m_loadingFromDem(false)
 116 , m_thread(0)
 117 , m_processingMutex()
 118 , m_stopProcessing(false)
 119 , m_btnOptimize(0)
 120 , m_configWidget(0)
 121 , m_spinIterations(0)
 122 , m_plot(new QwtPlot())
 123 , m_stdDevCurve(new QwtPlotCurve())
 124 , m_stdDevVoidCurve(new QwtPlotCurve())
 125 , m_stdDevNonVoidCurve(new QwtPlotCurve())
 126 , m_snapshotFile()
 127 , m_snapshotReadFile()
 128 , m_snapshotReading(false)
 129 , m_snapshotSlider(0)
 130 {
 131   m_plot->setWindowFlags(Qt::Tool | Qt::WindowStaysOnTopHint);
 132   m_plot->setTitle("Standard Deviations");
 133   connect(this, SIGNAL(replotStdDeviations()), m_plot, SLOT(replot()), Qt::QueuedConnection);
 134   connect(this, SIGNAL(invalidateAndRedraw()), this, SLOT(invalidateAndRedrawSlot()), Qt::BlockingQueuedConnection);
 135   m_stdDevVoidCurve->setPen(QPen(Qt::black));
 136   m_stdDevVoidCurve->attach(m_plot);
 137   m_stdDevNonVoidCurve->setPen(QPen(Qt::green));
 138   m_stdDevNonVoidCurve->attach(m_plot);
 139   m_stdDevCurve->setPen(QPen(Qt::blue));
 140   m_stdDevCurve->attach(m_plot);
 141   for (int i = 0; i < m_numDatasets; ++i)
 142   {
 143     tcShadowClassifyingData* landsat = datasets[i];
 144     m_datasets[i] = landsat;
 145     m_texToDatasetTex[i] = imagery->texToGeo() * landsat->geoToTex();
 146     m_shadowChannels[i] = landsat->shadowClassification();
 147     if (0 != m_shadowChannels[i])
 148     {
 149       m_shadowChannels[i]->addDerivitive(this);
 150     }
 151     m_shadingChannels[i] = landsat->shading();
 152     if (0 != m_shadingChannels[i])
 153     {
 154       m_shadingChannels[i]->addDerivitive(this);
 155     }
 156     m_datum[i] = 0;
 157   }
 158 }
 159
 160 /// Destructor.
 161 tcElevationOptimization::~tcElevationOptimization()
 162 {
 163   if (!m_processingMutex.tryLock())
 164   {
 165     m_stopProcessing = true;
 166     m_thread->wait();
 167   }
 168   delete m_thread;
 169   delete m_plot;
 170   delete [] m_datasets;
 171   delete [] m_texToDatasetTex;
 172   for (int i = 0; i < m_numDatasets; ++i)
 173   {
 174     if (0 != m_shadowChannels[i])
 175     {
 176       m_shadowChannels[i]->removeDerivitive(this);
 177     }
 178     if (0 != m_shadingChannels[i])
 179     {
 180       m_shadingChannels[i]->removeDerivitive(this);
 181     }
 182     delete m_datum[i];
 183   }
 184   delete [] m_shadowChannels;
 185   delete [] m_shadingChannels;
 186   delete m_elevationData;
 187   delete [] m_datum;
 188   delete m_data;
 189   delete [] m_shadowData;
 190   delete [] m_shadingData;
 191   delete m_configWidget;
 192 }
 193
 194 /*
 195  * Main image interface
 196  */
 197
 198 tcChannelConfigWidget* tcElevationOptimization::configWidget()
 199 {
 200   if (0 == m_configWidget)
 201   {
 202     m_configWidget = new tcChannelConfigWidget();
 203     m_configWidget->setWindowTitle(tr("%1 Configuration").arg(name()));
 204
 205     QVBoxLayout* layout = new QVBoxLayout(m_configWidget);
 206
 207     QLabel* stats = new QLabel(m_configWidget);
 208       layout->addWidget(stats);
 209       connect(this, SIGNAL(statistics(const QString&)), stats, SLOT(setText(const QString&)));
 210
 211     // Row of optimisation controls
 212     QWidget* buttons1 = new QWidget(m_configWidget);
 213       layout->addWidget(buttons1);
 214       QHBoxLayout* buttons1Layout = new QHBoxLayout(buttons1);
 215
 216       // reset DEM (sets all corrections to the original (interpolated) values
 217       QPushButton* btnResetDem = new QPushButton(tr("Reset DEM"), m_configWidget);
 218         buttons1Layout->addWidget(btnResetDem);
 219         btnResetDem->setToolTip(tr("Reset DEM corrected values to interpolated values"));
 220         connect(btnResetDem, SIGNAL(clicked(bool)), this, SLOT(resetDem()));
 221
 222       // load base from DEM (loads DEM elevations into elevation buffer)
 223       QPushButton* btnLoadDem = new QPushButton(tr("Load DEM"), m_configWidget);
 224         buttons1Layout->addWidget(btnLoadDem);
 225         btnLoadDem->setToolTip(tr("Load corrected values from DEM into elevation field"));
 226         connect(btnLoadDem, SIGNAL(clicked(bool)), this, SLOT(loadFromDem()));
 227
 228       // number of iterations (select number of iterations to perform)
 229       m_spinIterations = new QSpinBox(m_configWidget);
 230         buttons1Layout->addWidget(m_spinIterations);
 231         m_spinIterations->setToolTip(tr("Number of iterations of optimization to perform"));
 232         m_spinIterations->setRange(1, 1000);
 233         m_spinIterations->setValue(10);
 234
 235       // optimize (iterates on elevation buffer and write back to DEM at intervals and end)
 236       m_btnOptimize = new QPushButton(tr("Optimize Iteratively"), m_configWidget);
 237         buttons1Layout->addWidget(m_btnOptimize);
 238         m_btnOptimize->setToolTip(tr("Iteratively optimize elevation field and write results back to DEM"));
 239         connect(m_btnOptimize, SIGNAL(clicked(bool)), this, SLOT(startStopOptimize()));
 240
 241     // Another row of optimisation controls
 242     QWidget* buttons2 = new QWidget(m_configWidget);
 243       layout->addWidget(buttons2);
 244       QHBoxLayout* buttons2Layout = new QHBoxLayout(buttons2);
 245
 246       // apply hard constraints, adjusting reliabilities
 247       QPushButton* btnHardConstrain = new QPushButton(tr("Hard Constrain"), m_configWidget);
 248         buttons2Layout->addWidget(btnHardConstrain);
 249         btnLoadDem->setToolTip(tr("Apply hard shadow constraints"));
 250         connect(btnHardConstrain, SIGNAL(clicked(bool)), this, SLOT(applyHardConstraints()));
 251
 252       // bilinear interpolation of voids
 253       QPushButton* btnBilinear = new QPushButton(tr("Bilinear"), m_configWidget);
 254         buttons2Layout->addWidget(btnBilinear);
 255         btnLoadDem->setToolTip(tr("Void-fill using bilinear interpolation"));
 256         connect(btnBilinear, SIGNAL(clicked(bool)), this, SLOT(voidFillBilinear()));
 257
 258       // bicubic interpolation of voids
 259       QPushButton* btnBicubic = new QPushButton(tr("Bicubic"), m_configWidget);
 260         buttons2Layout->addWidget(btnBicubic);
 261         btnLoadDem->setToolTip(tr("Void-fill using bicubic interpolation"));
 262         connect(btnBicubic, SIGNAL(clicked(bool)), this, SLOT(voidFillBicubic()));
 263
 264     // Another row of controls, for graphs
 265     QWidget* buttons3 = new QWidget(m_configWidget);
 266       layout->addWidget(buttons3);
 267       QHBoxLayout* buttons3Layout = new QHBoxLayout(buttons3);
 268
 269       // Show standard deviation curve
 270       QPushButton* btnStdDev = new QPushButton(tr("Std Dev"), m_configWidget);
 271         buttons3Layout->addWidget(btnStdDev);
 272         btnStdDev->setToolTip(tr("Show standard deviation graphs"));
 273         connect(btnStdDev, SIGNAL(clicked(bool)), m_plot, SLOT(show()));
 274
 275       // Clear standard deviation
 276       QPushButton* btnStdDevClear = new QPushButton(tr("Clear"), m_configWidget);
 277         buttons3Layout->addWidget(btnStdDevClear);
 278         btnStdDevClear->setToolTip(tr("Clear standard deviation graphs"));
 279         connect(btnStdDevClear, SIGNAL(clicked(bool)), this, SLOT(clearStdDeviations()));
 280
 281       // Sample standard deviation
 282       QPushButton* btnStdDevSample = new QPushButton(tr("Sample"), m_configWidget);
 283         buttons3Layout->addWidget(btnStdDevSample);
 284         btnStdDevSample->setToolTip(tr("Sample standard deviation"));
 285         connect(btnStdDevSample, SIGNAL(clicked(bool)), this, SLOT(sampleStdDeviations()));
 286
 287       // Save standard deviation
 288       QPushButton* btnStdDevSave = new QPushButton(tr("Save"), m_configWidget);
 289         buttons3Layout->addWidget(btnStdDevSave);
 290         btnStdDevSave->setToolTip(tr("Save standard deviation samples to a text file for external processing"));
 291         connect(btnStdDevSave, SIGNAL(clicked(bool)), this, SLOT(saveStdDeviations()));
 292
 293     // Another row of controls, for terrain snapshotting
 294     QWidget* buttons4 = new QWidget(m_configWidget);
 295       layout->addWidget(buttons4);
 296       QHBoxLayout* buttons4Layout = new QHBoxLayout(buttons4);
 297
 298       // Initialise snapshotting
 299       QPushButton* btnInitSnapshotting = new QPushButton(tr("Start snapshotting"), m_configWidget);
 300         buttons4Layout->addWidget(btnInitSnapshotting);
 301         btnInitSnapshotting->setToolTip(tr("Start snapshotting the terrain at intervals for later analysis"));
 302         connect(btnInitSnapshotting, SIGNAL(clicked(bool)), this, SLOT(initSnapShotting()));
 303
 304       // Load snapshots
 305       QPushButton* btnLoadSnapshots = new QPushButton(tr("Load snapshots"), m_configWidget);
 306         buttons4Layout->addWidget(btnLoadSnapshots);
 307         btnLoadSnapshots->setToolTip(tr("Load a sequence of snapshots for processing"));
 308         connect(btnLoadSnapshots, SIGNAL(clicked(bool)), this, SLOT(loadSnapShots()));
 309
 310     // Snapshot slider
 311     m_snapshotSlider = new QSlider(Qt::Horizontal, m_configWidget);
 312       layout->addWidget(m_snapshotSlider);
 313       m_snapshotSlider->setEnabled(false);
 314       connect(m_snapshotSlider, SIGNAL(sliderMoved(int)), this, SLOT(loadSnapShot(int)));
 315
 316     // Rows of weight sliders
 317 #define NUM_WEIGHTS 9
 318     QString weightNames[NUM_WEIGHTS] = {
 319       tr("Variance"),
 320       tr("Roughness"),
 321       tr("Steepness"),
 322       tr("Lit facing sun"),
 323       tr("Shape from shading"),
 324       tr("Below shadow line"),
 325       tr("Transition elevation"),
 326       tr("Transition slope"),
 327       tr("Transition curving down")
 328     };
 329     QWidget* sliders = new QWidget(m_configWidget);
 330       layout->addWidget(sliders);
 331       QGridLayout* slidersLayout = new QGridLayout(sliders);
 332
 333       int i = 0;
 334       for (; i < NUM_WEIGHTS; ++i)
 335       {
 336         QLabel* label = new QLabel(weightNames[i], sliders);
 337           slidersLayout->addWidget(label, i, 0);
 338
 339         QSlider* slider1 = new QSlider(Qt::Horizontal, sliders);
 340           slidersLayout->addWidget(slider1, i, 1);
 341
 342         QSlider* slider2 = new QSlider(Qt::Horizontal, sliders);
 343           slidersLayout->addWidget(slider2, i, 2);
 344       }
 345
 346       QLabel* deltaHLabel = new QLabel(tr("Elevation Step"), sliders);
 347         slidersLayout->addWidget(deltaHLabel, i, 0);
 348         for (int s = 0; s < 2; ++s)
 349         {
 350           m_deltaHSlider[s] = new QSlider(Qt::Horizontal, sliders);
 351             slidersLayout->addWidget(m_deltaHSlider[s], i, 1+s);
 352             m_deltaHSlider[s]->setRange(1,100);
 353             m_deltaHSlider[s]->setValue(40);
 354         }
 355       ++i;
 356
 357       QLabel* rangeLabel = new QLabel(tr("Range"), sliders);
 358         slidersLayout->addWidget(rangeLabel, i, 0);
 359         for (int s = 0; s < 2; ++s)
 360         {
 361           m_rangeSlider[s] = new QSlider(Qt::Horizontal, sliders);
 362             slidersLayout->addWidget(m_rangeSlider[s], i, 1+s);
 363             m_rangeSlider[s]->setRange(0,4);
 364             m_rangeSlider[s]->setValue(0);
 365         }
 366       ++i;
 367   }
 368   return m_configWidget;
 369 }
 370
 371 /*
 372  * Slots
 373  */
 374
 375 /// Reset DEM.
 376 void tcElevationOptimization::resetDem()
 377 {
 378 }
 379
 380 /// Load elevation data from DEM.
 381 void tcElevationOptimization::loadFromDem()
 382 {
 383   delete m_elevationData;
 384   m_elevationData = 0;
 385
 386   delete m_data;
 387   m_data = 0;
 388
 389   for (int i = 0; i < m_numDatasets; ++i)
 390   {
 391     delete m_datum[i];
 392     m_datum[i] = 0;
 393   }
 394
 395   m_loadingFromDem = true;
 396   invalidate();
 397   m_configWidget->requestRedraw();
 398   m_loadingFromDem = false;
 399 }
 400
 401 /// Apply hard constraints to elevation data.
 402 void tcElevationOptimization::applyHardConstraints()
 403 {
 404   if (0 != m_elevationData)
 405   {
 406     const int width = m_elevationData->width();
 407     const int height = m_elevationData->height();
 408     // Now start optimising
 409     startProcessing("Applying hard constraints");
 410     for (int j = 0; j <= height; ++j)
 411     {
 412       progress((float)j/(height-1));
 413       for (int i = 0; i <= width; ++i)
 414       {
 415         int index = j*width + i;
 416
 417         // Only apply constraints to unreliable pixels
 418         PerPixelCache* cache = &m_data->buffer()[index];
 419         if (cache->reliability != 0)
 420         {
 421           continue;
 422         }
 423
 424         int elevationConstraintCount = 0;
 425
 426         // Per dataset
 427         for (int ds = 0; ds < m_numDatasets; ++ds)
 428         {
 429           PerDatasetPixelCache* dsCache = &m_datum[ds]->buffer()[index];
 430
 431           bool isShadowed = (m_shadowData[ds]->sampleFloat((float)i/(width-1), (float)j/(height-1)) < 0.5f);
 432           bool isShadowEntrance = isShadowed &&
 433                                   dsCache->shadowTransition[TransitionEntrance][0] == i &&
 434                                   dsCache->shadowTransition[TransitionEntrance][1] == j;
 435           bool isShadowExit = isShadowed &&
 436                               dsCache->shadowTransition[TransitionExit][0] == i &&
 437                               dsCache->shadowTransition[TransitionExit][1] == j;
 438
 439           if (isShadowed)
 440           {
 441             bool exitReliable = false;
 442             float exitElevation = 0.0f;
 443             float exitDistance = 0.0f;
 444             if (dsCache->shadowTransition[TransitionExit][0] >= 0 &&
 445                 dsCache->shadowTransition[TransitionExit][0] < width &&
 446                 dsCache->shadowTransition[TransitionExit][1] >= 0 &&
 447                 dsCache->shadowTransition[TransitionExit][1] < height)
 448             {
 449               int exitIndex = width*dsCache->shadowTransition[TransitionExit][1] + dsCache->shadowTransition[TransitionExit][0];
 450               PerPixelCache* exitCache = &m_data->buffer()[exitIndex];
 451               // Great, we can say for definite the elevation of this pixel because the corresponding pixel is reliable
 452               if (exitCache->reliability != 0)
 453               {
 454                 exitReliable = true;
 455                 exitElevation = m_elevationData->buffer()[exitIndex];
 456                 exitDistance = dsCache->shadowTransitionDistance[TransitionExit];
 457                 if (isShadowEntrance)
 458                 {
 459                   float newElevation = exitElevation + exitDistance*dsCache->sinLightElevation;
 460                   m_elevationData->buffer()[index] = m_elevationData->buffer()[index]*elevationConstraintCount + newElevation;
 461                   m_elevationData->buffer()[index] /= ++elevationConstraintCount;
 462                   cache->originalElevation = m_elevationData->buffer()[index];
 463                   cache->reliability = 1;
 464                 }
 465               }
 466             }
 467             bool entranceReliable = false;
 468             float entranceElevation = 0.0f;
 469             float entranceDistance = 0.0f;
 470             if (dsCache->shadowTransition[TransitionEntrance][0] >= 0 &&
 471                 dsCache->shadowTransition[TransitionEntrance][0] < width &&
 472                 dsCache->shadowTransition[TransitionEntrance][1] >= 0 &&
 473                 dsCache->shadowTransition[TransitionEntrance][1] < height)
 474             {
 475               int entranceIndex = width*dsCache->shadowTransition[TransitionEntrance][1] + dsCache->shadowTransition[TransitionEntrance][0];
 476               PerPixelCache* entranceCache = &m_data->buffer()[entranceIndex];
 477               // Great, we can say for definite the elevation of this pixel because the corresponding pixel is reliable
 478               if (entranceCache->reliability != 0)
 479               {
 480                 entranceReliable = true;
 481                 entranceElevation = m_elevationData->buffer()[entranceIndex];
 482                 entranceDistance = dsCache->shadowTransitionDistance[TransitionEntrance];
 483                 if (isShadowExit)
 484                 {
 485                   float newElevation = entranceElevation - entranceDistance*dsCache->sinLightElevation;
 486                   m_elevationData->buffer()[index] = m_elevationData->buffer()[index]*elevationConstraintCount + newElevation;
 487                   m_elevationData->buffer()[index] /= ++elevationConstraintCount;
 488                   cache->originalElevation = m_elevationData->buffer()[index];
 489                   cache->reliability = 1;
 490                 }
 491               }
 492             }
 493
 494             if (exitReliable && entranceReliable)
 495             {
 496               // maximum elevation is between elevation at entrance and elevation at exit
 497               float maxElevation = (exitElevation*entranceDistance + entranceElevation*exitDistance) / (entranceDistance + exitDistance);
 498               if (m_elevationData->buffer()[index] > maxElevation)
 499               {
 500                 m_elevationData->buffer()[index] = maxElevation;
 501               }
 502             }
 503           }
 504         }
 505       }
 506     }
 507     endProcessing();
 508
 509     invalidate();
 510     m_configWidget->requestRedraw();
 511   }
 512   else
 513   {
 514     QMessageBox::warning(m_configWidget,
 515                          tr("Elevation field not initialised."),
 516                          tr("Please ensure the optimisation channel is displayed."));
 517   }
 518 }
 519
 520 /// Void-fill bilinear.
 521 void tcElevationOptimization::voidFillBilinear()
 522 {
 523   if (0 != m_elevationData)
 524   {
 525     const int width = m_elevationData->width();
 526     const int height = m_elevationData->height();
 527     // Now start optimising
 528     startProcessing("Bilinear void filling");
 529     for (int j = 0; j < height; ++j)
 530     {
 531       progress(0.5f*j/(height-1));
 532       int lastNonVoid = -1;
 533       float lastNonVoidValue;
 534       for (int i = 0; i < width; ++i)
 535       {
 536         int index = j*width + i;
 537         if (m_data->buffer()[index].reliability != 0)
 538         {
 539           float value = m_elevationData->buffer()[index];
 540           if (lastNonVoid != -1)
 541           {
 542             int gap = i - lastNonVoid;
 543             float startAlt = lastNonVoidValue;
 544             float dAlt = (value - lastNonVoidValue) / gap;
 545             // Lovely, between lastNonVoid and i is a void
 546             for (int ii = lastNonVoid+1; ii < i; ++ii)
 547             {
 548               int iIndex = j*width + ii;
 549               float alt = startAlt + dAlt * (ii - lastNonVoid);
 550               m_elevationData->buffer()[iIndex] = alt;
 551               // Keep track of the gap
 552               m_data->buffer()[iIndex].temporary.i = gap;
 553             }
 554           }
 555           else
 556           {
 557             m_data->buffer()[index].temporary.i = -1;
 558           }
 559           lastNonVoid = i;
 560           lastNonVoidValue = value;
 561         }
 562         else
 563         {
 564           m_data->buffer()[index].temporary.i = -1;
 565         }
 566       }
 567     }
 568     for (int i = 0; i < width; ++i)
 569     {
 570       progress(0.5f + 0.5f*i/(width-1));
 571       int lastNonVoid = -1;
 572       float lastNonVoidValue;
 573       for (int j = 0; j < height; ++j)
 574       {
 575         int index = j*width + i;
 576         float value = m_elevationData->buffer()[index];
 577         if (m_data->buffer()[index].reliability != 0)
 578         {
 579           if (lastNonVoid != -1)
 580           {
 581             int gap = j - lastNonVoid;
 582             float startAlt = lastNonVoidValue;
 583             float dAlt = (value - lastNonVoidValue) / gap;
 584             // Lovely, between lastNonVoid and j is a void
 585             for (int jj = lastNonVoid+1; jj < j; ++jj)
 586             {
 587               // Average with previous value if non zero, otherwise overwrite
 588               int iIndex = jj*width + i;
 589               float alt = startAlt + dAlt * (jj - lastNonVoid);
 590               int otherGap = m_data->buffer()[iIndex].temporary.i;
 591               if (otherGap != -1)
 592               {
 593                 float prevAlt = m_elevationData->buffer()[iIndex];
 594                 // Prefer the value of the smaller gap
 595                 alt = (alt * otherGap + prevAlt * gap) / (gap+otherGap);
 596               }
 597               m_elevationData->buffer()[iIndex] = alt;
 598             }
 599           }
 600           lastNonVoid = j;
 601           lastNonVoidValue = value;
 602         }
 603       }
 604     }
 605     endProcessing();
 606
 607     invalidate();
 608     m_configWidget->requestRedraw();
 609   }
 610   else
 611   {
 612     QMessageBox::warning(m_configWidget,
 613                          tr("Elevation field not initialised."),
 614                          tr("Please ensure the optimisation channel is displayed."));
 615   }
 616 }
 617
 618 /// Void-fill bicubic.
 619 void tcElevationOptimization::voidFillBicubic()
 620 {
 621   if (0 != m_elevationData)
 622   {
 623     const int width = m_elevationData->width();
 624     const int height = m_elevationData->height();
 625     // Now start optimising
 626     startProcessing("Bicubic void filling");
 627     for (int j = 0; j < height; ++j)
 628     {
 629       progress(0.5f*j/(height-1));
 630       int lastNonVoid = -1;
 631       float lastNonVoidValue;
 632       for (int i = 0; i < width; ++i)
 633       {
 634         int index = j*width + i;
 635         if (m_data->buffer()[index].reliability != 0)
 636         {
 637           float value = m_elevationData->buffer()[index];
 638           if (lastNonVoid != -1)
 639           {
 640             int gap = i - lastNonVoid;
 641             // Obtain some control points for the splines
 642             // use the two samples either side of the gap as control points 1 and 2
 643             // use the ones before and after these samples, extended to the size of
 644             // the gap as control points 0 and 3
 645             float startAlt = lastNonVoidValue;
 646             float beforeStartAlt = startAlt;
 647             if (lastNonVoid > 0)
 648             {
 649               if (0 != m_data->buffer()[j*width + lastNonVoid - 1].reliability)
 650               {
 651                 beforeStartAlt = startAlt + (m_elevationData->buffer()[j*width + lastNonVoid - 1] - startAlt) * gap;
 652               }
 653             }
 654             float endAlt = value;
 655             float afterEndAlt = endAlt;
 656             if (i < width-1)
 657             {
 658               if (0 != m_data->buffer()[j*width + i + 1].reliability)
 659               {
 660                 afterEndAlt = endAlt + (m_elevationData->buffer()[j*width + i + 1] - endAlt) * gap;
 661               }
 662             }
 663             // Lovely, between lastNonVoid and i is a void
 664             for (int ii = lastNonVoid+1; ii < i; ++ii)
 665             {
 666               int iIndex = j*width + ii;
 667               float u = (float)(ii-lastNonVoid) / gap;
 668               float u2 = u*u;
 669               float u3 = u2*u;
 670               float alt = maths::catmullRomSpline.Spline(u, u2, u3, beforeStartAlt, startAlt, endAlt, afterEndAlt);
 671               m_elevationData->buffer()[iIndex] = alt;
 672               // Keep track of the gap
 673               m_data->buffer()[iIndex].temporary.i = gap;
 674             }
 675           }
 676           else
 677           {
 678             m_data->buffer()[index].temporary.i = -1;
 679           }
 680           lastNonVoid = i;
 681           lastNonVoidValue = value;
 682         }
 683         else
 684         {
 685           m_data->buffer()[index].temporary.i = -1;
 686         }
 687       }
 688     }
 689     for (int i = 0; i < width; ++i)
 690     {
 691       progress(0.5f + 0.5f*i/(width-1));
 692       int lastNonVoid = -1;
 693       float lastNonVoidValue;
 694       for (int j = 0; j < height; ++j)
 695       {
 696         int index = j*width + i;
 697         float value = m_elevationData->buffer()[index];
 698         if (m_data->buffer()[index].reliability != 0)
 699         {
 700           if (lastNonVoid != -1)
 701           {
 702             int gap = j - lastNonVoid;
 703             // Obtain some control points for the splines
 704             // use the two samples either side of the gap as control points 1 and 2
 705             // use the ones before and after these samples, extended to the size of
 706             // the gap as control points 0 and 3
 707             float startAlt = lastNonVoidValue;
 708             float beforeStartAlt = startAlt;
 709             if (lastNonVoid > 0)
 710             {
 711               if (0 != m_data->buffer()[(lastNonVoid-1)*width + i].reliability)
 712               {
 713                 beforeStartAlt = startAlt + (m_elevationData->buffer()[(lastNonVoid-1)*width + i] - startAlt) * gap;
 714               }
 715             }
 716             float endAlt = value;
 717             float afterEndAlt = endAlt;
 718             if (j < height-1)
 719             {
 720               if (0 != m_data->buffer()[(j+1)*width + i].reliability)
 721               {
 722                 afterEndAlt = endAlt + (m_elevationData->buffer()[(j+1)*width + i] - endAlt) * gap;
 723               }
 724             }
 725             // Lovely, between lastNonVoid and j is a void
 726             for (int jj = lastNonVoid+1; jj < j; ++jj)
 727             {
 728               // Average with previous value if non zero, otherwise overwrite
 729               int iIndex = jj*width + i;
 730               float u = (float)(jj-lastNonVoid) / gap;
 731               float u2 = u*u;
 732               float u3 = u2*u;
 733               float alt = maths::catmullRomSpline.Spline(u, u2, u3, beforeStartAlt, startAlt, endAlt, afterEndAlt);
 734               int otherGap = m_data->buffer()[iIndex].temporary.i;
 735               if (otherGap != -1)
 736               {
 737                 float prevAlt = m_elevationData->buffer()[iIndex];
 738                 // Prefer the value of the smaller gap
 739                 alt = (alt * otherGap + prevAlt * gap) / (gap+otherGap);
 740               }
 741               m_elevationData->buffer()[iIndex] = alt;
 742             }
 743           }
 744           lastNonVoid = j;
 745           lastNonVoidValue = value;
 746         }
 747       }
 748     }
 749     endProcessing();
 750
 751     invalidate();
 752     m_configWidget->requestRedraw();
 753   }
 754   else
 755   {
 756     QMessageBox::warning(m_configWidget,
 757                          tr("Elevation field not initialised."),
 758                          tr("Please ensure the optimisation channel is displayed."));
 759   }
 760 }
 761
 762 /// Start/stop optimization process.
 763 void tcElevationOptimization::startStopOptimize()
 764 {
 765   if (m_processingMutex.tryLock())
 766   {
 767     m_btnOptimize->setText(tr("Stop Optimization"));
 768     m_btnOptimize->update();
 769     if (!m_thread)
 770     {
 771       m_thread = new WorkerThread(this);
 772     }
 773     m_thread->start();
 774   }
 775   else
 776   {
 777     // Warn the worker thread that its time to stop.
 778     m_stopProcessing = true;
 779   }
 780 }
 781
 782 /// Optimize a number of times.
 783 void tcElevationOptimization::optimize()
 784 {
 785   if (0 != m_elevationData)
 786   {
 787     const int width = m_elevationData->width();
 788     const int height = m_elevationData->height();
 789     // Now start optimising
 790     startProcessing("Optimizing elevation data");
 791     int passes = m_spinIterations->value();
 792     for (int i = 0; i < passes; ++i)
 793     {
 794       startProcessing("Optimizing elevation data");
 795       progress((float)i/passes);
 796       m_weights.variance = 0.1f;//1.0f
 797       m_weights.roughness = 1.0f;//1.0f;
 798       m_weights.steepness = 1.0f;//1.0f;
 799       m_weights.litFacing = 255.0f;//255.0f;
 800       m_weights.shading = 0.0f;//100.0f;
 801       m_weights.belowShadowLine = 1.0f;//1.0f;
 802       m_weights.transitionElev = 1e-1f;//1e-1f;
 803       m_weights.transitionTangential = 100.0f;//100.0f;
 804       m_weights.transitionCurvingDown = 1.0f;//1.0f;
 805       float passage = (float)i/(passes-1);
 806       float passageNeg = 1.0f - passageNeg;
 807       float deltaH = 0.1f*(passageNeg*m_deltaHSlider[0]->value() + passage*m_deltaHSlider[1]->value());;
 808       float range = passageNeg*m_rangeSlider[0]->value() + passage*m_rangeSlider[1]->value();;
 809       optimizeElevations(0, 0, width-1, height-1, deltaH, range);
 810       if (m_stopProcessing)
 811       {
 812         emit invalidateAndRedraw();
 813         endProcessing();
 814         break;
 815       }
 816
 817       if ((i+1) % 5 == 0 || i == passes-1)
 818       {
 819         emit invalidateAndRedraw();
 820         startProcessing("Optimizing elevation data");
 821         progress((float)i/passes);
 822       }
 823
 824       sampleStdDeviations();
 825     }
 826     endProcessing();
 827   }
 828   else
 829   {
 830     QMessageBox::warning(m_configWidget,
 831                          tr("Elevation field not initialised."),
 832                          tr("Please ensure the optimisation channel is displayed."));
 833   }
 834 }
 835
 836 /// Clear standard deviations.
 837 void tcElevationOptimization::clearStdDeviations()
 838 {
 839   m_plotX.clear();
 840   m_stdDevData.clear();
 841   m_stdDevCurve->setData(m_plotX, m_stdDevData);
 842   m_stdDevVoidData.clear();
 843   m_stdDevVoidCurve->setData(m_plotX, m_stdDevVoidData);
 844   m_stdDevNonVoidData.clear();
 845   m_stdDevNonVoidCurve->setData(m_plotX, m_stdDevNonVoidData);
 846   m_plot->replot();
 847   m_plot->update();
 848 }
 849
 850 /// Sample standard deviations.
 851 void tcElevationOptimization::sampleStdDeviations()
 852 {
 853   int width = m_elevationData->width();
 854   int height = m_elevationData->height();
 855   float stdDevVoid = 0.0f;
 856   float stdDevNonVoid = 0.0f;
 857   float stdDevAll = 0.0f;
 858   int index = 0;
 859   int voids = 0;
 860   for (int j = 0; j < height; ++j)
 861   {
 862     for (int i = 0; i < width; ++i)
 863     {
 864       float altitude = m_data->buffer()[index].accurateElevation;
 865
 866       float square = m_elevationData->buffer()[index] - altitude;
 867       square *= square;
 868       stdDevAll += square;
 869       if (m_data->buffer()[index].reliability)
 870       {
 871         stdDevNonVoid += square;
 872       }
 873       else
 874       {
 875         stdDevVoid += square;
 876         ++voids;
 877       }
 878       ++index;
 879     }
 880   }
 881   stdDevNonVoid = sqrt(stdDevNonVoid / (index-voids));
 882   stdDevVoid = sqrt(stdDevVoid / voids);
 883   stdDevAll = sqrt(stdDevAll / index);
 884   emit statistics(tr("Overall std deviation: %1 m\n"
 885                      "Void std deviation: %2 m\n"
 886                      "Non-void std deviation: %3 m")
 887                   .arg(stdDevAll).arg(stdDevVoid).arg(stdDevNonVoid));
 888   //m_configWidget->update();
 889
 890   if (m_snapshotReadFile.isEmpty() || m_snapshotReading)
 891   {
 892     m_plotX.append(m_plotX.size()+1);
 893     m_stdDevData.append(stdDevAll);
 894     m_stdDevCurve->setData(m_plotX, m_stdDevData);
 895     m_stdDevVoidData.append(stdDevVoid);
 896     m_stdDevVoidCurve->setData(m_plotX, m_stdDevVoidData);
 897     m_stdDevNonVoidData.append(stdDevNonVoid);
 898     m_stdDevNonVoidCurve->setData(m_plotX, m_stdDevNonVoidData);
 899     // this will get queued for event loop
 900     emit replotStdDeviations();
 901
 902     // Take a snapshot
 903     if (!m_snapshotFile.isEmpty())
 904     {
 905       QString fname = QString("%1.%2").arg(m_snapshotFile).arg(m_plotX.size(), 4, 10, QLatin1Char('0'));
 906       QFile snap(fname);
 907       if (snap.open(QIODevice::WriteOnly))
 908       {
 909         QDataStream stream(&snap);
 910         stream.setByteOrder(QDataStream::LittleEndian);
 911         stream << width << height;
 912         int index = 0;
 913         for (int j = 0; j < height; ++j)
 914         {
 915           for (int i = 0; i < width; ++i)
 916           {
 917             stream << m_elevationData->buffer()[index];
 918             ++index;
 919           }
 920         }
 921       }
 922     }
 923   }
 924 }
 925
 926 /// Sample standard deviations.
 927 void tcElevationOptimization::saveStdDeviations()
 928 {
 929   QString filename = QFileDialog::getSaveFileName(m_configWidget,
 930                         tr("Save standard deviation data"),
 931                         QString(),
 932                         tr("ASCII Data Files (*.adf)"));
 933
 934   if (!filename.isEmpty())
 935   {
 936     // Open the file ready to write the data
 937     QFile file(filename);
 938     if (!file.open(QIODevice::WriteOnly | QIODevice::Text))
 939     {
 940       QMessageBox::warning(m_configWidget,
 941           tr("Open failed."),
 942           tr("Could not open '%1' for writing.").arg(filename));
 943       return;
 944     }
 945
 946     QTextStream out(&file);
 947
 948     for (int i = 0; i < m_plotX.size(); ++i)
 949     {
 950       out << m_plotX[i] << "\t"
 951           << m_stdDevData[i] << "\t"
 952           << m_stdDevVoidData[i] << "\t"
 953           << m_stdDevNonVoidData[i] << "\n";
 954     }
 955   }
 956 }
 957
 958 /// Invalidate and redraw.
 959 void tcElevationOptimization::invalidateAndRedrawSlot()
 960 {
 961   invalidate();
 962   m_configWidget->requestRedraw();
 963 }
 964
 965 /// Initialise snapshotting.
 966 void tcElevationOptimization::initSnapShotting()
 967 {
 968   m_snapshotFile = QFileDialog::getSaveFileName(m_configWidget,
 969                         tr("Save snapshotting config file"),
 970                         QString(),
 971                         tr("Snapshot Config Files (*.sdf)"));
 972
 973   if (!m_snapshotFile.isEmpty())
 974   {
 975     // Open the file ready to write the data
 976     QFile file(m_snapshotFile);
 977     if (!file.open(QIODevice::WriteOnly))
 978     {
 979       QMessageBox::warning(m_configWidget,
 980           tr("Open failed."),
 981           tr("Could not open '%1' for writing.").arg(m_snapshotFile));
 982       return;
 983     }
 984
 985     const double* portion = portionPosition();
 986     tcGeo p1(imagery()->texToGeo()*maths::Vector<2,double>(portion[0], portion[1]));
 987     tcGeo p2(imagery()->texToGeo()*maths::Vector<2,double>(portion[2], portion[3]));
 988
 989     QDataStream stream(&file);
 990     stream.setByteOrder(QDataStream::LittleEndian);
 991     stream << p1.lon()*180.0/M_PI
 992            << p1.lat()*180.0/M_PI
 993            << p2.lon()*180.0/M_PI
 994            << p2.lat()*180.0/M_PI;
 995   }
 996 }
 997
 998 #include <QtDebug>
 999 /// Load a sequence of snapshots.
1000 void tcElevationOptimization::loadSnapShots()
1001 {
1002   m_snapshotFile = QString();
1003   m_snapshotReadFile = QFileDialog::getOpenFileName(m_configWidget,
1004                         tr("Open snapshotting config file"),
1005                         QString(),
1006                         tr("Snapshot Config Files (*.sdf)"));
1007
1008   if (!m_snapshotReadFile.isEmpty())
1009   {
1010     // Open the file ready to write the data
1011     QFile file(m_snapshotReadFile);
1012     if (!file.open(QIODevice::ReadOnly | QIODevice::Text))
1013     {
1014       QMessageBox::warning(m_configWidget,
1015           tr("Open failed."),
1016           tr("Could not open '%1' for reading.").arg(m_snapshotReadFile));
1017       return;
1018     }
1019
1020     QDataStream stream(&file);
1021     stream.setByteOrder(QDataStream::LittleEndian);
1022
1023     double portion[2][2];
1024     stream >> portion[0][0]
1025            >> portion[0][1]
1026            >> portion[1][0]
1027            >> portion[1][1];
1028     tcGeo p1(portion[0][0]*M_PI/180.0, portion[0][1]*M_PI/180.0);
1029     tcGeo p2(portion[1][0]*M_PI/180.0, portion[1][1]*M_PI/180.0);
1030
1031     // force update of region
1032     m_configWidget->requestSlice(p1, p2);
1033
1034     // Find number of snaps
1035     int max = 0;
1036     bool exists = true;
1037     QFile snap;
1038     while (exists)
1039     {
1040       ++max;
1041       QString fname = QString("%1.%2").arg(m_snapshotReadFile).arg(max, 4, 10, QLatin1Char('0'));
1042       snap.setFileName(fname);
1043       exists = snap.exists();
1044     }
1045     --max;
1046     m_snapshotSlider->setRange(1, max);
1047     m_snapshotSlider->setEnabled(true);
1048
1049     // load first sample
1050     clearStdDeviations();
1051     m_snapshotReading = true;
1052     for (int i = 1; i <= max; ++i)
1053     {
1054       loadSnapShot(i);
1055     }
1056     m_snapshotReading = false;
1057     loadSnapShot(max);
1058     m_snapshotSlider->setValue(max);
1059   }
1060   else
1061   {
1062     m_snapshotSlider->setEnabled(false);
1063   }
1064 }
1065
1066 /// Load a single snapshot.
1067 void tcElevationOptimization::loadSnapShot(int id)
1068 {
1069   // Load a snapshot
1070   if (!m_snapshotReadFile.isEmpty())
1071   {
1072     QString fname = QString("%1.%2").arg(m_snapshotReadFile).arg(id, 4, 10, QLatin1Char('0'));
1073     QFile snap(fname);
1074     if (snap.open(QIODevice::ReadOnly))
1075     {
1076       QDataStream stream(&snap);
1077       stream.setByteOrder(QDataStream::LittleEndian);
1078       int width, height;
1079       stream >> width >> height;
1080       if (m_elevationData->width() != width || m_elevationData->height() != height)
1081       {
1082         QMessageBox::warning(m_configWidget,
1083                              tr("Snapshot dimentions do not match with current."),
1084                              tr("Panic!"));
1085       }
1086       else
1087       {
1088         int index = 0;
1089         for (int j = 0; j < height; ++j)
1090         {
1091           for (int i = 0; i < width; ++i)
1092           {
1093             stream >> m_elevationData->buffer()[index];
1094             ++index;
1095           }
1096         }
1097         sampleStdDeviations();
1098         if (!m_snapshotReading)
1099         {
1100           invalidate();
1101           m_configWidget->requestRedraw();
1102         }
1103       }
1104     }
1105   }
1106 }
1107
1108 /*
1109  * Interface for derived class to implement
1110  */
1111
1112 void tcElevationOptimization::roundPortion(double* x1, double* y1, double* x2, double* y2)
1113 {
1114   //m_shadowChannel->roundPortion(x1,y1,x2,y2);
1115 }
1116
1117 tcAbstractPixelData* tcElevationOptimization::loadPortion(double x1, double y1, double x2, double y2, bool changed)
1118 {
1119   int width = -1;
1120   int height = -1;
1121
1122   for (int i = 0; i < m_numDatasets; ++i)
1123   {
1124     maths::Vector<2,double> xy1 = m_texToDatasetTex[i] * maths::Vector<2,double>(x1,y1);
1125     maths::Vector<2,double> xy2 = m_texToDatasetTex[i] * maths::Vector<2,double>(x2,y2);
1126     Reference<tcAbstractPixelData> channelData = m_shadowChannels[i]->portion(xy1[0], xy1[1], xy2[0], xy2[1]);
1127     if (i == 0)
1128     {
1129       width = channelData->width();
1130       height = channelData->height();
1131     }
1132     m_shadowData[i] = dynamicCast<tcPixelData<float>*>(channelData);
1133     if (0 == m_shadowData[i])
1134     {
1135       return 0;
1136     }
1137     if (0 != m_shadingChannels[i])
1138     {
1139       Reference<tcAbstractPixelData> shadingData = m_shadingChannels[i]->portion(xy1[0], xy1[1], xy2[0], xy2[1]);
1140       m_shadingData[i] = dynamicCast<tcPixelData<float>*>(shadingData);
1141     }
1142     else
1143     {
1144       m_shadingData[i] = 0;
1145     }
1146   }
1147
1148   /// @todo What if its a different size because the portion has changed?
1149   if (0 == m_elevationData || 0 == m_data || changed)
1150   {
1151     delete m_elevationData;
1152     for (int i = 0; i < m_numDatasets; ++i)
1153     {
1154       delete m_datum[i];
1155     }
1156     delete m_data;
1157
1158     // Create a new pixel buffer, initialised to altitude
1159
1160     m_elevationData = new tcElevationData(width, height);
1161     // Find transform of elevation data
1162     tcGeo geoBase = geoAt(maths::Vector<2,float>(x1, y1));
1163     tcGeo geoDiagonal = geoAt(maths::Vector<2,float>(x1 + (x2-x1)/(width-1), y1 + (y2-y1)/(height-1))) - geoBase;
1164     m_elevationData->setGeoToPixels(tcAffineTransform2<double>(geoBase, geoDiagonal).inverse());
1165
1166     m_data = new tcTypedPixelData<PerPixelCache>(width, height);
1167     for (int i = 0; i < m_numDatasets; ++i)
1168     {
1169       m_datum[i] = new tcTypedPixelData<PerDatasetPixelCache>(width, height);
1170     }
1171
1172     int maxWidthHeight = qMax(width, height);
1173     double minXY12 = qMin(x2-x1, y2-y1);
1174
1175     startProcessing("Caching elevation characteristics");
1176     for (int j = 0; j < height; ++j)
1177     {
1178       progress((float)j/(height-1));
1179       for (int i = 0; i < width; ++i)
1180       {
1181         int index = j*width + i;
1182         // Transform coordinates
1183         maths::Vector<2,float> coord      (x1 + (x2-x1)*i / (width  - 1),
1184                                            y1 + (y2-y1)*j / (height - 1));
1185         tcGeo                  geoCoord = geoAt(coord);
1186
1187         PerPixelCache* cache = &m_data->buffer()[index];
1188
1189         // Find resolution
1190         maths::Vector<2,float> coordx      (x1 + (x2-x1)*(i+1) / (width  - 1),
1191                                             y1 + (y2-y1)*j     / (height - 1));
1192         maths::Vector<2,float> coordy      (x1 + (x2-x1)*i     / (width  - 1),
1193                                             y1 + (y2-y1)*(j+1) / (height - 1));
1194         tcGeo                  geoCoordx = geoAt(coordx) - geoCoord;
1195         tcGeo                  geoCoordy = geoAt(coordy) - geoCoord;
1196         // Find approx
1197         float res = geoCoordx.angle();
1198         cache->resolution[0] = res;
1199         cache->resolution[1] = geoCoordy.angle();
1200         cache->resolution *= 6378.137e3;
1201
1202         // Get some elevation data
1203         bool accurate;
1204         float altitude        = altitudeAt(geoCoord, true, &accurate);
1205         m_elevationData->buffer()[index] = altitude;
1206         cache->originalElevation = altitude;
1207         cache->reliability = accurate ? 1 : 0;
1208         cache->accurateElevation = dem()[1].altitudeAt(geoCoord, true, 0);
1209
1210         // Per dataset cache
1211         for (int ds = 0; ds < m_numDatasets; ++ds)
1212         {
1213           tcGeoImageData* imagery = m_datasets[ds];
1214           maths::Vector<3,float> light = lightDirectionAt(geoCoord, imagery);
1215
1216           PerDatasetPixelCache* dsCache = &m_datum[ds]->buffer()[index];
1217           dsCache->light = light;
1218           dsCache->sinLightElevation = light[2]/light.slice<0,2>().mag();
1219
1220           tcGeo lightScaled(light[0]*res/sin(geoCoord.lat()),
1221                             light[1]*res);
1222           tcGeo offset(geoCoord + lightScaled);
1223           for (int depthDirection = 0; depthDirection < 2; ++depthDirection)
1224           {
1225             tcGeo pos = geoCoord;
1226             tcGeo delta;
1227             if (depthDirection == TransitionEntrance)
1228             {
1229               delta = offset - geoCoord;
1230             }
1231             else
1232             {
1233               delta = geoCoord - offset;
1234             }
1235             delta = delta / 2;
1236             pos = pos + delta;
1237             maths::Vector<2,float> tex = coord;
1238             int transitionI = i;
1239             int transitionJ = j;
1240             while (transitionI >= 0 && transitionI < width && transitionJ >= 0 && transitionJ < height)
1241             {
1242               if (m_shadowData[ds]->sampleFloat((tex[0]-x1)/(x2-x1), (tex[1]-y1)/(y2-y1)) > 0.5f)
1243               {
1244                 break;
1245               }
1246               transitionI = 0.5f + (tex[0]-x1)/(x2-x1)*(width -1);
1247               transitionJ = 0.5f + (tex[1]-y1)/(y2-y1)*(height-1);
1248               pos = pos + delta;
1249               tex = textureAt(pos);
1250             }
1251             dsCache->shadowTransition[depthDirection][0] = transitionI;
1252             dsCache->shadowTransition[depthDirection][1] = transitionJ;
1253             dsCache->shadowTransitionDistance[depthDirection] = (pos-geoCoord).angle()*6378.137e3;
1254           }
1255         }
1256       }
1257     }
1258   }
1259
1260   if (!m_loadingFromDem)
1261   {
1262     // Update elevation model
1263     startProcessing("Updating elevation model");
1264     dem()->updateFromElevationData(m_elevationData);
1265   }
1266
1267   // Downscale the elevation for viewing
1268   startProcessing("Downscaling for viewing");
1269   tcPixelData<GLubyte>* result = new tcPixelData<GLubyte>(width, height);
1270   //tcTypedPixelData<maths::Vector<3,float> >* result = new tcTypedPixelData<maths::Vector<3,float> >(width, height);
1271   Q_ASSERT(0 != result);
1272
1273   int index = 0;
1274   for (int j = 0; j < height; ++j)
1275   {
1276     progress((float)j/(height-1));
1277     for (int i = 0; i < width; ++i)
1278     {
1279       PerDatasetPixelCache* dsCache = &m_datum[0]->buffer()[index];
1280       bool isShadowed = (m_shadowData[0]->sampleFloat((float)i/(width-1), (float)j/(height-1)) < 0.5f);
1281       bool isShadowEntrance = isShadowed &&
1282                               dsCache->shadowTransition[TransitionEntrance][0] == i &&
1283                               dsCache->shadowTransition[TransitionEntrance][1] == j;
1284       bool isShadowExit = isShadowed &&
1285                           dsCache->shadowTransition[TransitionExit][0] == i &&
1286                           dsCache->shadowTransition[TransitionExit][1] == j;
1287       int xdist = abs(dsCache->shadowTransition[TransitionEntrance][0] - i);
1288       int ydist = abs(dsCache->shadowTransition[TransitionEntrance][1] - j);
1289       bool bothIn = isShadowed && (dsCache->shadowTransition[TransitionExit][0] >= 0 &&
1290               dsCache->shadowTransition[TransitionExit][0] < width &&
1291               dsCache->shadowTransition[TransitionExit][1] >= 0 &&
1292               dsCache->shadowTransition[TransitionExit][1] < height &&
1293               dsCache->shadowTransition[TransitionEntrance][0] >= 0 &&
1294               dsCache->shadowTransition[TransitionEntrance][0] < width &&
1295               dsCache->shadowTransition[TransitionEntrance][1] >= 0 &&
1296               dsCache->shadowTransition[TransitionEntrance][1] < height);
1297
1298       //*
1299       result->buffer()[index] = isShadowed ? qMin(255, 100 * (xdist + ydist)) : 255;
1300       result->buffer()[index] = bothIn?255:0;
1301       result->buffer()[index] = cost(i, j, i, j);
1302       //*/
1303
1304 #if 0
1305       /*
1306        * matrix L={Light row vec of shading channel i;..}
1307        * invert L
1308        * multiply by shading values for shading channels
1309        * should be mostly about the same magnitude
1310        * normalize
1311        */
1312
1313 #define SHAPEFROMSHADING_COMBOS 1
1314       int chans[SHAPEFROMSHADING_COMBOS][3] = {
1315         /*{0,1,2},
1316         {0,1,3},
1317         {0,1,4},
1318         {0,2,3},
1319         {0,2,4},
1320         {0,3,4},
1321         {1,2,3},
1322         {1,2,4},
1323         {1,3,4},*/
1324         {2,3,4},
1325       };
1326       maths::Vector<3,float> totalN(0.0f);
1327       int comb = 0;
1328       for (; comb < SHAPEFROMSHADING_COMBOS; ++comb)
1329       {
1330         maths::Matrix<3,float> L(m_datum[chans[comb][0]]->buffer()[index].light,
1331                                  m_datum[chans[comb][1]]->buffer()[index].light,
1332                                  m_datum[chans[comb][2]]->buffer()[index].light);
1333         maths::Vector<3,float> S;
1334         L.invert();
1335         for (int ds = 0; ds < 3; ++ds)
1336         {
1337           S[ds] = m_shadingData[chans[comb][ds]]->sampleFloat((float)i/(width-1), (float)j/(height-1));
1338         }
1339         maths::Vector<3,float> N = S * L;
1340         N.normalize();
1341         totalN += N;
1342       }
1343       result->buffer()[index] = totalN / comb;
1344       result->buffer()[index].normalize();
1345 #endif
1346
1347       ++index;
1348     }
1349   }
1350
1351   endProcessing();
1352
1353   return result;
1354 }
1355
1356 template <typename T>
1357 T MySqr(T t)
1358 {
1359   return t*t;
1360 }
1361
1362 /*
1363  * Private functions
1364  */
1365
1366 /// Calculate the cost for a set of pixels.
1367 float tcElevationOptimization::cost(int x1, int y1, int x2, int y2) const
1368 {
1369   int width = m_elevationData->width();
1370   int height = m_elevationData->height();
1371   if (x1 < 0)
1372   {
1373     x1 = 0;
1374   }
1375   if (x2 >= width)
1376   {
1377     x2 = width-1;
1378   }
1379   if (y1 < 0)
1380   {
1381     y1 = 0;
1382   }
1383   if (y2 >= height)
1384   {
1385     y2 = height-1;
1386   }
1387
1388   int maxWidthHeight = qMax(width, height);
1389   //double minXY12 = qMin(m_window.x2-m_window.x1, m_window.y2-m_window.y1);
1390
1391   float costVariance = 0.0f;
1392   int countVariance = 0;
1393   float costRoughness = 0.0f;
1394   int countRoughness = 0;
1395   float costSteepness = 0.0f;
1396   int countSteepness = 0;
1397   float costLitFacing = 0.0f;
1398   int countLitFacing = 0;
1399   float costShading = 0.0f;
1400   int countShading = 0;
1401   float costBelowShadowLine = 0.0f;
1402   int countBelowShadowLine = 0;
1403   float costTransitionElev = 0.0f;
1404   int countTransitionElev = 0;
1405   float costTransitionTangential = 0.0f;
1406   int countTransitionTangential = 0;
1407   float costTransitionCurvingDown = 0.0f;
1408   int countTransitionCurvingDown = 0;
1409   for (int j = y1; j <= y2; ++j)
1410   {
1411     for (int i = x1; i <= x2; ++i)
1412     {
1413       float* elevationPtr = pixelElevation(i, j);
1414       if (0 == elevationPtr)
1415       {
1416         continue;
1417       }
1418
1419       int index = j*width + i;
1420
1421       // Basic data retrieval
1422       float elevation = *elevationPtr;
1423       PerPixelCache* cache = &m_data->buffer()[index];
1424       if (m_voidsOnly && cache->reliability == 0)
1425       {
1426         continue;
1427       }
1428
1429       // Gradient
1430       // Derivitive of gradient
1431       float dh_dx = 0.0f;
1432       float dh_dy = 0.0f;
1433       float d2h_dx2 = 0.0f;
1434       float d2h_dy2 = 0.0f;
1435       if (i > 0 && i < width-1)
1436       {
1437         float hp1 = m_elevationData->buffer()[index+1];
1438         float hm1 = m_elevationData->buffer()[index-1];
1439         d2h_dx2 = (hp1+hm1)/2 - elevation;
1440         dh_dx = (hp1-hm1)/2 / cache->resolution[0];
1441       }
1442       if (j > 0 && j < height-1)
1443       {
1444         float hp1 = m_elevationData->buffer()[index+width];
1445         float hm1 = m_elevationData->buffer()[index-width];
1446         d2h_dy2 = (hp1+hm1)/2 - elevation;
1447         dh_dy = (hp1-hm1)/2 / cache->resolution[1];
1448       }
1449
1450       // Minimise variance from original value
1451       costVariance += cache->reliability * (elevation-cache->originalElevation)*(elevation-cache->originalElevation);
1452       ++countVariance;
1453
1454       // Minimise roughness
1455       costRoughness += (d2h_dx2*d2h_dx2 + d2h_dy2*d2h_dy2);
1456       ++countRoughness;
1457
1458       // Also minimise slope
1459       // This works really well at getting a nice slope
1460       costSteepness += (dh_dx*dh_dx + dh_dy*dh_dy);
1461       ++countSteepness;
1462
1463       // Per dataset
1464       for (int ds = 0; ds < m_numDatasets; ++ds)
1465       {
1466         PerDatasetPixelCache* dsCache = &m_datum[ds]->buffer()[index];
1467
1468         bool isShadowed = (m_shadowData[ds]->sampleFloat((float)i/(width-1), (float)j/(height-1)) < 0.5f);
1469         bool isShadowEntrance = isShadowed &&
1470                                 dsCache->shadowTransition[TransitionEntrance][0] == i &&
1471                                 dsCache->shadowTransition[TransitionEntrance][1] == j;
1472         bool isShadowExit = isShadowed &&
1473                             dsCache->shadowTransition[TransitionExit][0] == i &&
1474                             dsCache->shadowTransition[TransitionExit][1] == j;
1475
1476         // Calculate measures relating to the light direction
1477         // Gradient
1478         float hwm1 = m_elevationData->sampleFloat(((float)i - dsCache->light[0])/(width-1),
1479             ((float)j - dsCache->light[1])/(height-1));
1480         float hwp1 = m_elevationData->sampleFloat(((float)i + dsCache->light[0])/(width-1),
1481             ((float)j + dsCache->light[1])/(height-1));
1482         /// @todo Ensure units are right (this is in terms of h metres, w pixels)
1483         float dh_dw = (hwm1 - hwp1)/2 / cache->resolution[0]; // 30 metre resolution?
1484         float d2h_dw2 = (hwm1+hwp1)/2 - elevation;
1485
1486         if (isShadowed)
1487         {
1488           // maximum elevation is between elevation at entrance and elevation at exit
1489           if (dsCache->shadowTransition[TransitionExit][0] >= 0 &&
1490               dsCache->shadowTransition[TransitionExit][0] < width &&
1491               dsCache->shadowTransition[TransitionExit][1] >= 0 &&
1492               dsCache->shadowTransition[TransitionExit][1] < height &&
1493               dsCache->shadowTransition[TransitionEntrance][0] >= 0 &&
1494               dsCache->shadowTransition[TransitionEntrance][0] < width &&
1495               dsCache->shadowTransition[TransitionEntrance][1] >= 0 &&
1496               dsCache->shadowTransition[TransitionEntrance][1] < height)
1497           {
1498             int exitIndex = width*dsCache->shadowTransition[TransitionExit][1] + dsCache->shadowTransition[TransitionExit][0];
1499             int entranceIndex = width*dsCache->shadowTransition[TransitionEntrance][1] + dsCache->shadowTransition[TransitionEntrance][0];
1500             float exitElevation = m_elevationData->buffer()[exitIndex];
1501             float entranceElevation = m_elevationData->buffer()[entranceIndex];
1502
1503             costBelowShadowLine += MySqr(qMax(0.0f, elevation - (exitElevation     * dsCache->shadowTransitionDistance[TransitionEntrance]
1504                                                                 +entranceElevation * dsCache->shadowTransitionDistance[TransitionExit])
1505                                                               / (dsCache->shadowTransitionDistance[TransitionEntrance]
1506                                                                 +dsCache->shadowTransitionDistance[TransitionExit])));
1507             ++countBelowShadowLine;
1508           }
1509         }
1510         else
1511         {
1512           float facingness = -dh_dw+dsCache->light[0];
1513           if (facingness < 0.0f)
1514           {
1515             // Lit pixels must face light
1516             costLitFacing += facingness*facingness;
1517             ++countLitFacing;
1518           }
1519           // Simple shape from shading
1520           if (0 != m_shadingData[ds])
1521           {
1522             float intensity = m_shadingData[ds]->sampleFloat((float)i/(width-1), (float)j/(height-1));
1523             // intensity = cos(theta) = L.N
1524             maths::Vector<3,float> xTan(1.0f, 0.0f, dh_dx);
1525             maths::Vector<3,float> yTan(0.0f, 1.0f, dh_dy);
1526             costShading += MySqr(acos(qMin(1.0f, intensity)) - acos(qMin(1.0f, dsCache->light*maths::cross(xTan, yTan).normalized())));
1527             ++countShading;
1528           }
1529         }
1530         if (isShadowEntrance)
1531         {
1532           if (dsCache->shadowTransition[TransitionExit][0] >= 0 &&
1533               dsCache->shadowTransition[TransitionExit][0] < width &&
1534               dsCache->shadowTransition[TransitionExit][1] >= 0 &&
1535               dsCache->shadowTransition[TransitionExit][1] < height)
1536           {
1537             int exitIndex = width*dsCache->shadowTransition[TransitionExit][1] + dsCache->shadowTransition[TransitionExit][0];
1538             float exitElevation = m_elevationData->buffer()[exitIndex];
1539             costTransitionElev += MySqr(exitElevation + dsCache->shadowTransitionDistance[TransitionExit]*dsCache->sinLightElevation - elevation);
1540             ++countTransitionElev;
1541           }
1542           // gradient tangential to sun direction
1543           costTransitionTangential += MySqr(dh_dw - dsCache->sinLightElevation);
1544           ++countTransitionTangential;
1545           // second derivitive should be negative at entrance
1546           costTransitionCurvingDown += MySqr(qMax(0.0f, d2h_dw2));
1547           ++countTransitionCurvingDown;
1548         }
1549         if (0 && isShadowExit)
1550         {
1551           if (dsCache->shadowTransition[TransitionEntrance][0] >= 0 &&
1552               dsCache->shadowTransition[TransitionEntrance][0] < width &&
1553               dsCache->shadowTransition[TransitionEntrance][1] >= 0 &&
1554               dsCache->shadowTransition[TransitionEntrance][1] < height)
1555           {
1556             int entranceIndex = width*dsCache->shadowTransition[TransitionEntrance][1] + dsCache->shadowTransition[TransitionEntrance][0];
1557             float entranceElevation = m_elevationData->buffer()[entranceIndex];
1558             // this has less effect, its usually "more" right
1559             /// @todo Smooth high reliability into low reliability
1560             costTransitionElev += 0.1f * MySqr(entranceElevation - dsCache->shadowTransitionDistance[TransitionEntrance]*dsCache->sinLightElevation - elevation);
1561             ++countTransitionElev;
1562           }
1563         }
1564       }
1565     }
1566   }
1567
1568   return m_weights.variance*costVariance//countVariance
1569        + m_weights.roughness*costRoughness//countRoughness
1570        + m_weights.steepness*costSteepness//countSteepness
1571        + m_weights.litFacing*costLitFacing//countLitFacing
1572        + m_weights.shading*costShading//countShading
1573        + m_weights.belowShadowLine*costBelowShadowLine//countBelowShadowLine
1574        + m_weights.transitionElev*costTransitionElev//countTransitionElev
1575        + m_weights.transitionTangential*costTransitionTangential//countTransitionTangential
1576        + m_weights.transitionCurvingDown*costTransitionCurvingDown//countTransitionCurvingDown
1577        ;
1578 }
1579
1580 /*
1581  * f(0) = 1
1582  * f(range) = small
1583  * f(x) = exp[ - (2*x/range)^2 ]
1584  */
1585 inline float elevationFalloffFunction(float xsq)
1586 {
1587   return qMax(0.0f, 1.0f - xsq);
1588   //return exp(-4*xsq);
1589 }
1590 /*
1591  */
1592
1593 /// Calculate the relative cost of changing a pixel's elevation.
1594 float tcElevationOptimization::costChange(int x, int y, float dh, int range)
1595 {
1596   // Override range for now
1597 #define ELEV_RANGE 4
1598   range = qMin(range, ELEV_RANGE);
1599   static const int effectiveRange = 1 + range;
1600   float oldElevations[1+2*ELEV_RANGE][1+2*ELEV_RANGE];
1601
1602   // Assume a single change in elevation only affects cost of pixels to effectiveRange from it.
1603   // Now technically this doesn't hold for shadow edges as all pixels across the shadow can be affected.
1604   // If the pixel is a shadow edge, take into account the cost of extra pixels.
1605   float currentCost = cost(x-effectiveRange, y-effectiveRange, x+effectiveRange, y+effectiveRange);
1606   for (int i = 0; i < 1+range*2; ++i)
1607   {
1608     int di = i - range;
1609     for (int j = 0; j < 1+range*2; ++j)
1610     {
1611       int dj = j - range;
1612       float* elev = pixelElevation(x+di, y+dj);
1613       if (0 != elev && (!m_voidsOnly || 0 != m_data->buffer()[j*m_data->width() + i].reliability))
1614       {
1615         oldElevations[i][j] = *elev;
1616         *elev += dh*elevationFalloffFunction(float(di*di + dj*dj) / ((range-1)*(range-1)));
1617       }
1618     }
1619   }
1620   float newCost = cost(x-effectiveRange, y-effectiveRange, x+effectiveRange, y+effectiveRange);
1621   for (int i = 0; i < 1+range*2; ++i)
1622   {
1623     int di = i - range;
1624     for (int j = 0; j < 1+range*2; ++j)
1625     {
1626       int dj = j - range;
1627       float* elev = pixelElevation(x+di, y+dj);
1628       if (0 != elev && (!m_voidsOnly || 0 != m_data->buffer()[j*m_data->width() + i].reliability))
1629       {
1630         *elev = oldElevations[i][j];
1631       }
1632     }
1633   }
1634   return newCost - currentCost;
1635 }
1636
1637 /// Optimize a whole range of elevations.
1638 void tcElevationOptimization::optimizeElevations(int x1, int y1, int x2, int y2, float dh, int range)
1639 {
1640   for (int j = y1; j <= y2; ++j)
1641   {
1642     for (int i = x1; i <= x2; ++i)
1643     {
1644       if (m_stopProcessing)
1645       {
1646         return;
1647       }
1648       int index = j*m_data->width() + i;
1649       //while (true)
1650       if (!m_voidsOnly || 0 != m_data->buffer()[j*m_data->width() + i].reliability)
1651       {
1652         // Find the cost of adjusting the elevation by dh in each direction
1653         float costInc = costChange(i,j, dh, range);
1654         float costDec = costChange(i,j,-dh, range);
1655         if (costInc < 0.0f && costInc < costDec)
1656         {
1657           adjustElevation(i,j,dh, range);
1658         }
1659         else if (costDec < 0.0f)
1660         {
1661           adjustElevation(i,j,-dh, range);
1662         }
1663         else
1664         {
1665           //break;
1666         }
1667       }
1668     }
1669   }
1670 }
1671
1672 /// Get pointer to elevation at a pixel.
1673 float* tcElevationOptimization::pixelElevation(int x, int y) const
1674 {
1675   if (x < 0 || x >= m_elevationData->width() ||
1676       y < 0 || y >= m_elevationData->height())
1677   {
1678     return 0;
1679   }
1680   int index = y*m_elevationData->width() + x;
1681   return &(m_elevationData->buffer()[index]);
1682 }
1683
1684 /// Adjust elevation at a pixel.
1685 void tcElevationOptimization::adjustElevation(int x, int y, float dh, int range)
1686 {
1687   range = qMin(range, ELEV_RANGE);
1688   for (int i = 0; i < 1+range*2; ++i)
1689   {
1690     int di = i - range;
1691     for (int j = 0; j < 1+range*2; ++j)
1692     {
1693       int dj = j - range;
1694       float* elev = pixelElevation(x+di, y+dj);
1695       if (0 != elev)
1696       {
1697         *elev += dh*elevationFalloffFunction(float(di*di + dj*dj) / ((range+1)*(range+1)));
1698       }
1699     }
1700   }
1701 }