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