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