From d8136668b1fb70e620f3be97787382752e4ee31d Mon Sep 17 00:00:00 2001 From: James Hogan Date: Tue, 17 Mar 2009 23:36:26 +0000 Subject: [PATCH] Final changes before project hand-in --- geo/tcChannel.cpp | 2 +- geo/tcChannel.h | 4 + geo/tcCustomSunDirection.cpp | 6 +- geo/tcCustomSunDirection.h | 2 +- geo/tcElevationData.cpp | 7 + geo/tcElevationData.h | 3 + geo/tcElevationOptimization.cpp | 805 ++++++++++++++++++++++++++++++---------- geo/tcElevationOptimization.h | 82 ++-- geo/tcGeo.h | 7 + geo/tcGlobe.cpp | 21 +- geo/tcNegativeProduct.cpp | 117 ++++-- geo/tcNegativeProduct.h | 15 +- geo/tcSrtmModel.cpp | 20 +- maths/Matrix.cpp | 30 +- maths/Matrix.h | 6 +- tcViewportWidget.cpp | 4 +- 16 files changed, 837 insertions(+), 294 deletions(-) diff --git a/geo/tcChannel.cpp b/geo/tcChannel.cpp index 77b203b..0fe871b 100644 --- a/geo/tcChannel.cpp +++ b/geo/tcChannel.cpp @@ -106,7 +106,7 @@ void tcChannel::newCrossSection(const tcGeo& p1, const tcGeo& p2) { m_crossSectioned = true; m_crossSection[0] = p1; - m_crossSection[1] = p1; + m_crossSection[1] = p2; } /* diff --git a/geo/tcChannel.h b/geo/tcChannel.h index 36f9ee7..cc87108 100644 --- a/geo/tcChannel.h +++ b/geo/tcChannel.h @@ -200,12 +200,16 @@ class tcChannel : public QObject /// Radiance scaling information. float m_radianceOffset, m_radianceGain; + protected: + /// Whether cross sectioned. bool m_crossSectioned; /// Current cross section. tcGeo m_crossSection[2]; + private: + /* * Static variables */ diff --git a/geo/tcCustomSunDirection.cpp b/geo/tcCustomSunDirection.cpp index 70e2097..d66d4f3 100644 --- a/geo/tcCustomSunDirection.cpp +++ b/geo/tcCustomSunDirection.cpp @@ -35,7 +35,7 @@ */ /// Primary constructor. -tcCustomSunDirection::tcCustomSunDirection(tcSrtmModel* dem, const tcGeo& sunDirection, const tcGeo& sunReference) +tcCustomSunDirection::tcCustomSunDirection(tcSrtmModel* dem, const tcGeo& sunDirection, const tcGeo& sunReference, bool shading) : tcShadowClassifyingData() { setName(QObject::tr("Custom sun direction")); @@ -54,6 +54,10 @@ tcCustomSunDirection::tcCustomSunDirection(tcSrtmModel* dem, const tcGeo& sunDir tcRaytracedShadowMap* raytrace = new tcRaytracedShadowMap(tcGeo(pixelRes, pixelRes), dem, this); channelManager()->addChannel(raytrace); setShadowClassification(raytrace); + + tcLambertianShading* diffuse = new tcLambertianShading(channelManager()->channel(0), dem, this); + channelManager()->addChannel(diffuse); + setShading(diffuse); } /// Destructor. diff --git a/geo/tcCustomSunDirection.h b/geo/tcCustomSunDirection.h index 1b3a34f..a6118f9 100644 --- a/geo/tcCustomSunDirection.h +++ b/geo/tcCustomSunDirection.h @@ -39,7 +39,7 @@ class tcCustomSunDirection : public tcShadowClassifyingData */ /// Primary constructor. - tcCustomSunDirection(tcSrtmModel* dem, const tcGeo& sunDirection, const tcGeo& sunReference); + tcCustomSunDirection(tcSrtmModel* dem, const tcGeo& sunDirection, const tcGeo& sunReference, bool shading = false); /// Destructor. virtual ~tcCustomSunDirection(); diff --git a/geo/tcElevationData.cpp b/geo/tcElevationData.cpp index ae3d91b..cef5020 100644 --- a/geo/tcElevationData.cpp +++ b/geo/tcElevationData.cpp @@ -51,6 +51,13 @@ tcElevationData::~tcElevationData() * Elevation interface */ +/// Sample elevation at a coordinate. +float tcElevationData::elevationAt(const tcGeo& coord) const +{ + maths::Vector<2,double> pt = m_geoToPixels * coord; + return sampleFloat(pt[0]/(width()-1), pt[1]/(height()-1)); +} + /// Get the max elevation over a range of pixels. float tcElevationData::maxElevation(const tcGeo& swCorner, const tcGeo& neCorner, bool* inRange) const { diff --git a/geo/tcElevationData.h b/geo/tcElevationData.h index c0bb95e..16ba2b9 100644 --- a/geo/tcElevationData.h +++ b/geo/tcElevationData.h @@ -61,6 +61,9 @@ class tcElevationData : public tcPixelData * Elevation interface */ + /// Sample elevation at a coordinate. + float elevationAt(const tcGeo& coord) const; + /// Get the max elevation over a range of pixels. float maxElevation(const tcGeo& swCorner, const tcGeo& neCorner, bool* inRange) const; diff --git a/geo/tcElevationOptimization.cpp b/geo/tcElevationOptimization.cpp index 0f547b2..5c9f213 100644 --- a/geo/tcElevationOptimization.cpp +++ b/geo/tcElevationOptimization.cpp @@ -107,8 +107,8 @@ tcElevationOptimization::tcElevationOptimization(const QList*[m_numDatasets]) -, m_data(0) +, m_data(new tcTypedPixelData*[m_numDatasets]) +, m_datum(0) , m_shadowData(new Reference >[m_numDatasets]) , m_shadingData(new Reference >[m_numDatasets]) , m_voidsOnly(false) @@ -127,6 +127,7 @@ tcElevationOptimization::tcElevationOptimization(const QListsetWindowFlags(Qt::Tool | Qt::WindowStaysOnTopHint); m_plot->setTitle("Standard Deviations"); @@ -153,7 +154,7 @@ tcElevationOptimization::tcElevationOptimization(const QListaddDerivitive(this); } - m_datum[i] = 0; + m_data[i] = 0; } } @@ -179,13 +180,13 @@ tcElevationOptimization::~tcElevationOptimization() { m_shadingChannels[i]->removeDerivitive(this); } - delete m_datum[i]; + delete m_data[i]; } delete [] m_shadowChannels; delete [] m_shadingChannels; delete m_elevationData; - delete [] m_datum; - delete m_data; + delete [] m_data; + delete m_datum; delete [] m_shadowData; delete [] m_shadingData; delete m_configWidget; @@ -307,6 +308,29 @@ tcChannelConfigWidget* tcElevationOptimization::configWidget() btnLoadSnapshots->setToolTip(tr("Load a sequence of snapshots for processing")); connect(btnLoadSnapshots, SIGNAL(clicked(bool)), this, SLOT(loadSnapShots())); + // Samples resolution + m_crossResolution = new QSpinBox(m_configWidget); + buttons4Layout->addWidget(m_crossResolution); + m_crossResolution->setToolTip(tr("Cross sectioning resolution")); + m_crossResolution->setRange(64, 1024); + m_crossResolution->setValue(1024); + + // Frequency of samples + m_crossFrequency = new QSpinBox(m_configWidget); + buttons4Layout->addWidget(m_crossFrequency); + m_crossFrequency->setToolTip(tr("Cross sectioning frequency")); + m_crossFrequency->setRange(1, 10); + m_crossFrequency->setValue(1); + + // Button to save cross sectioning data + QPushButton* btnCrossSectionSave = new QPushButton(tr("Save cross section data"), m_configWidget); + buttons4Layout->addWidget(btnCrossSectionSave); + connect(btnCrossSectionSave, SIGNAL(clicked(bool)), this, SLOT(saveCrossSections())); + + // Current frame + m_snapshotCounter = new QLabel("0", m_configWidget); + buttons4Layout->addWidget(m_snapshotCounter); + // Snapshot slider m_snapshotSlider = new QSlider(Qt::Horizontal, m_configWidget); layout->addWidget(m_snapshotSlider); @@ -314,56 +338,81 @@ tcChannelConfigWidget* tcElevationOptimization::configWidget() connect(m_snapshotSlider, SIGNAL(sliderMoved(int)), this, SLOT(loadSnapShot(int))); // Rows of weight sliders -#define NUM_WEIGHTS 9 - QString weightNames[NUM_WEIGHTS] = { + QString weightNames[WeightsMax] = { tr("Variance"), tr("Roughness"), tr("Steepness"), tr("Lit facing sun"), tr("Shape from shading"), + tr("Photometric Stereo"), tr("Below shadow line"), - tr("Transition elevation"), + tr("Transition elevation (top)"), + tr("Transition elevation (bottom)"), tr("Transition slope"), - tr("Transition curving down") + tr("Transition curving down"), + tr("Elevation Step"), + tr("Range") + }; + float defaults[WeightsMax][3] = { + {0.0f, 10.0f, 1.0f}, // variance + {0.0f, 5.0f, 0.1f}, // roughenss + {0.0f, 5.0f, 0.1f}, // steepness + {0.0f, 1000.0f, 255.0f}, // lit facing + {0.0f, 100.0f, 0.0f}, // shading + {0.0f, 100.0f, 0.0f}, // photometric stereo + {0.0f, 10.0f, 1.0f}, // belowShadowLine + {0.0f, 10.0f, 1.0f}, // transitionElevTop + {0.0f, 10.0f, 1.0f}, // transitionElevBottom + {0.0f, 1000.0f, 100.0f}, // transitionTangential + {0.0f, 10.0f, 1.0f}, // transitionCurvingDown + {0.0f, 20.0f, 4.0f}, // elevation step + {0.0f, 4.0f, 0.0f}, // range }; QWidget* sliders = new QWidget(m_configWidget); layout->addWidget(sliders); QGridLayout* slidersLayout = new QGridLayout(sliders); int i = 0; - for (; i < NUM_WEIGHTS; ++i) + for (; i < WeightsMax; ++i) { QLabel* label = new QLabel(weightNames[i], sliders); slidersLayout->addWidget(label, i, 0); - QSlider* slider1 = new QSlider(Qt::Horizontal, sliders); - slidersLayout->addWidget(slider1, i, 1); - - QSlider* slider2 = new QSlider(Qt::Horizontal, sliders); - slidersLayout->addWidget(slider2, i, 2); - } - - QLabel* deltaHLabel = new QLabel(tr("Elevation Step"), sliders); - slidersLayout->addWidget(deltaHLabel, i, 0); - for (int s = 0; s < 2; ++s) - { - m_deltaHSlider[s] = new QSlider(Qt::Horizontal, sliders); - slidersLayout->addWidget(m_deltaHSlider[s], i, 1+s); - m_deltaHSlider[s]->setRange(1,100); - m_deltaHSlider[s]->setValue(40); - } - ++i; + int steps = 1000; + m_weightLimits[0][i] = defaults[i][0]; + m_weightLimits[1][i] = (defaults[i][1] - defaults[i][0])/steps; + int val = (int)((defaults[i][2] - m_weightLimits[0][i]) / m_weightLimits[1][i]); - QLabel* rangeLabel = new QLabel(tr("Range"), sliders); - slidersLayout->addWidget(rangeLabel, i, 0); - for (int s = 0; s < 2; ++s) + for (int sl = 0; sl < 2; ++sl) { - m_rangeSlider[s] = new QSlider(Qt::Horizontal, sliders); - slidersLayout->addWidget(m_rangeSlider[s], i, 1+s); - m_rangeSlider[s]->setRange(0,4); - m_rangeSlider[s]->setValue(0); + m_weightSliders[sl][i] = new QSlider(Qt::Horizontal, sliders); + slidersLayout->addWidget(m_weightSliders[sl][i], i, sl+1); + m_weightSliders[sl][i]->setRange(0, steps); + m_weightSliders[sl][i]->setValue(val); + connect(m_weightSliders[sl][i], SIGNAL(sliderMoved(int)), this, SLOT(updateWeightLabels())); + connect(m_weightSliders[sl][i], SIGNAL(valueChanged(int)), this, SLOT(updateWeightLabels())); } - ++i; + m_weightLabels[i] = new QLabel("-", sliders); + slidersLayout->addWidget(m_weightLabels[i], i, 3); + } + updateWeightLabels(); + + // Another row of controls, for weight saving, loading + QWidget* buttons5 = new QWidget(m_configWidget); + layout->addWidget(buttons5); + QHBoxLayout* buttons5Layout = new QHBoxLayout(buttons5); + + // Save weights + QPushButton* btnSaveWeights = new QPushButton(tr("Save weights"), m_configWidget); + buttons5Layout->addWidget(btnSaveWeights); + btnInitSnapshotting->setToolTip(tr("Save the weights to a file")); + connect(btnSaveWeights, SIGNAL(clicked(bool)), this, SLOT(saveWeights())); + + // Load weights + QPushButton* btnLoadWeights = new QPushButton(tr("Load weights"), m_configWidget); + buttons5Layout->addWidget(btnLoadWeights); + btnLoadWeights->setToolTip(tr("Load the weights from a file")); + connect(btnLoadWeights, SIGNAL(clicked(bool)), this, SLOT(loadWeights())); } return m_configWidget; } @@ -383,13 +432,13 @@ void tcElevationOptimization::loadFromDem() delete m_elevationData; m_elevationData = 0; - delete m_data; - m_data = 0; + delete m_datum; + m_datum = 0; for (int i = 0; i < m_numDatasets; ++i) { - delete m_datum[i]; - m_datum[i] = 0; + delete m_data[i]; + m_data[i] = 0; } m_loadingFromDem = true; @@ -415,7 +464,7 @@ void tcElevationOptimization::applyHardConstraints() int index = j*width + i; // Only apply constraints to unreliable pixels - PerPixelCache* cache = &m_data->buffer()[index]; + PerPixelCache* cache = &m_datum->buffer()[index]; if (cache->reliability != 0) { continue; @@ -426,7 +475,7 @@ void tcElevationOptimization::applyHardConstraints() // Per dataset for (int ds = 0; ds < m_numDatasets; ++ds) { - PerDatasetPixelCache* dsCache = &m_datum[ds]->buffer()[index]; + PerDatasetPixelCache* dsCache = &m_data[ds]->buffer()[index]; bool isShadowed = (m_shadowData[ds]->sampleFloat((float)i/(width-1), (float)j/(height-1)) < 0.5f); bool isShadowEntrance = isShadowed && @@ -447,7 +496,7 @@ void tcElevationOptimization::applyHardConstraints() dsCache->shadowTransition[TransitionExit][1] < height) { int exitIndex = width*dsCache->shadowTransition[TransitionExit][1] + dsCache->shadowTransition[TransitionExit][0]; - PerPixelCache* exitCache = &m_data->buffer()[exitIndex]; + PerPixelCache* exitCache = &m_datum->buffer()[exitIndex]; // Great, we can say for definite the elevation of this pixel because the corresponding pixel is reliable if (exitCache->reliability != 0) { @@ -456,7 +505,7 @@ void tcElevationOptimization::applyHardConstraints() exitDistance = dsCache->shadowTransitionDistance[TransitionExit]; if (isShadowEntrance) { - float newElevation = exitElevation + exitDistance*dsCache->sinLightElevation; + float newElevation = exitElevation + exitDistance*dsCache->dlz_dlxy; m_elevationData->buffer()[index] = m_elevationData->buffer()[index]*elevationConstraintCount + newElevation; m_elevationData->buffer()[index] /= ++elevationConstraintCount; cache->originalElevation = m_elevationData->buffer()[index]; @@ -473,7 +522,7 @@ void tcElevationOptimization::applyHardConstraints() dsCache->shadowTransition[TransitionEntrance][1] < height) { int entranceIndex = width*dsCache->shadowTransition[TransitionEntrance][1] + dsCache->shadowTransition[TransitionEntrance][0]; - PerPixelCache* entranceCache = &m_data->buffer()[entranceIndex]; + PerPixelCache* entranceCache = &m_datum->buffer()[entranceIndex]; // Great, we can say for definite the elevation of this pixel because the corresponding pixel is reliable if (entranceCache->reliability != 0) { @@ -482,7 +531,7 @@ void tcElevationOptimization::applyHardConstraints() entranceDistance = dsCache->shadowTransitionDistance[TransitionEntrance]; if (isShadowExit) { - float newElevation = entranceElevation - entranceDistance*dsCache->sinLightElevation; + float newElevation = entranceElevation - entranceDistance*dsCache->dlz_dlxy; m_elevationData->buffer()[index] = m_elevationData->buffer()[index]*elevationConstraintCount + newElevation; m_elevationData->buffer()[index] /= ++elevationConstraintCount; cache->originalElevation = m_elevationData->buffer()[index]; @@ -534,7 +583,7 @@ void tcElevationOptimization::voidFillBilinear() for (int i = 0; i < width; ++i) { int index = j*width + i; - if (m_data->buffer()[index].reliability != 0) + if (m_datum->buffer()[index].reliability != 0) { float value = m_elevationData->buffer()[index]; if (lastNonVoid != -1) @@ -549,19 +598,19 @@ void tcElevationOptimization::voidFillBilinear() float alt = startAlt + dAlt * (ii - lastNonVoid); m_elevationData->buffer()[iIndex] = alt; // Keep track of the gap - m_data->buffer()[iIndex].temporary.i = gap; + m_datum->buffer()[iIndex].temporary.i = gap; } } else { - m_data->buffer()[index].temporary.i = -1; + m_datum->buffer()[index].temporary.i = -1; } lastNonVoid = i; lastNonVoidValue = value; } else { - m_data->buffer()[index].temporary.i = -1; + m_datum->buffer()[index].temporary.i = -1; } } } @@ -574,7 +623,7 @@ void tcElevationOptimization::voidFillBilinear() { int index = j*width + i; float value = m_elevationData->buffer()[index]; - if (m_data->buffer()[index].reliability != 0) + if (m_datum->buffer()[index].reliability != 0) { if (lastNonVoid != -1) { @@ -587,7 +636,7 @@ void tcElevationOptimization::voidFillBilinear() // Average with previous value if non zero, otherwise overwrite int iIndex = jj*width + i; float alt = startAlt + dAlt * (jj - lastNonVoid); - int otherGap = m_data->buffer()[iIndex].temporary.i; + int otherGap = m_datum->buffer()[iIndex].temporary.i; if (otherGap != -1) { float prevAlt = m_elevationData->buffer()[iIndex]; @@ -632,7 +681,7 @@ void tcElevationOptimization::voidFillBicubic() for (int i = 0; i < width; ++i) { int index = j*width + i; - if (m_data->buffer()[index].reliability != 0) + if (m_datum->buffer()[index].reliability != 0) { float value = m_elevationData->buffer()[index]; if (lastNonVoid != -1) @@ -646,7 +695,7 @@ void tcElevationOptimization::voidFillBicubic() float beforeStartAlt = startAlt; if (lastNonVoid > 0) { - if (0 != m_data->buffer()[j*width + lastNonVoid - 1].reliability) + if (0 != m_datum->buffer()[j*width + lastNonVoid - 1].reliability) { beforeStartAlt = startAlt + (m_elevationData->buffer()[j*width + lastNonVoid - 1] - startAlt) * gap; } @@ -655,7 +704,7 @@ void tcElevationOptimization::voidFillBicubic() float afterEndAlt = endAlt; if (i < width-1) { - if (0 != m_data->buffer()[j*width + i + 1].reliability) + if (0 != m_datum->buffer()[j*width + i + 1].reliability) { afterEndAlt = endAlt + (m_elevationData->buffer()[j*width + i + 1] - endAlt) * gap; } @@ -670,19 +719,19 @@ void tcElevationOptimization::voidFillBicubic() float alt = maths::catmullRomSpline.Spline(u, u2, u3, beforeStartAlt, startAlt, endAlt, afterEndAlt); m_elevationData->buffer()[iIndex] = alt; // Keep track of the gap - m_data->buffer()[iIndex].temporary.i = gap; + m_datum->buffer()[iIndex].temporary.i = gap; } } else { - m_data->buffer()[index].temporary.i = -1; + m_datum->buffer()[index].temporary.i = -1; } lastNonVoid = i; lastNonVoidValue = value; } else { - m_data->buffer()[index].temporary.i = -1; + m_datum->buffer()[index].temporary.i = -1; } } } @@ -695,7 +744,7 @@ void tcElevationOptimization::voidFillBicubic() { int index = j*width + i; float value = m_elevationData->buffer()[index]; - if (m_data->buffer()[index].reliability != 0) + if (m_datum->buffer()[index].reliability != 0) { if (lastNonVoid != -1) { @@ -708,7 +757,7 @@ void tcElevationOptimization::voidFillBicubic() float beforeStartAlt = startAlt; if (lastNonVoid > 0) { - if (0 != m_data->buffer()[(lastNonVoid-1)*width + i].reliability) + if (0 != m_datum->buffer()[(lastNonVoid-1)*width + i].reliability) { beforeStartAlt = startAlt + (m_elevationData->buffer()[(lastNonVoid-1)*width + i] - startAlt) * gap; } @@ -717,7 +766,7 @@ void tcElevationOptimization::voidFillBicubic() float afterEndAlt = endAlt; if (j < height-1) { - if (0 != m_data->buffer()[(j+1)*width + i].reliability) + if (0 != m_datum->buffer()[(j+1)*width + i].reliability) { afterEndAlt = endAlt + (m_elevationData->buffer()[(j+1)*width + i] - endAlt) * gap; } @@ -731,7 +780,7 @@ void tcElevationOptimization::voidFillBicubic() float u2 = u*u; float u3 = u2*u; float alt = maths::catmullRomSpline.Spline(u, u2, u3, beforeStartAlt, startAlt, endAlt, afterEndAlt); - int otherGap = m_data->buffer()[iIndex].temporary.i; + int otherGap = m_datum->buffer()[iIndex].temporary.i; if (otherGap != -1) { float prevAlt = m_elevationData->buffer()[iIndex]; @@ -792,20 +841,17 @@ void tcElevationOptimization::optimize() for (int i = 0; i < passes; ++i) { startProcessing("Optimizing elevation data"); - progress((float)i/passes); - m_weights.variance = 0.1f;//1.0f - m_weights.roughness = 1.0f;//1.0f; - m_weights.steepness = 1.0f;//1.0f; - m_weights.litFacing = 255.0f;//255.0f; - m_weights.shading = 0.0f;//100.0f; - m_weights.belowShadowLine = 1.0f;//1.0f; - m_weights.transitionElev = 1e-1f;//1e-1f; - m_weights.transitionTangential = 100.0f;//100.0f; - m_weights.transitionCurvingDown = 1.0f;//1.0f; float passage = (float)i/(passes-1); float passageNeg = 1.0f - passageNeg; - float deltaH = 0.1f*(passageNeg*m_deltaHSlider[0]->value() + passage*m_deltaHSlider[1]->value());; - float range = passageNeg*m_rangeSlider[0]->value() + passage*m_rangeSlider[1]->value();; + progress(passage); + for (int w = 0; w < WeightsMax; ++w) + { + float start = m_weightLimits[0][w] + m_weightLimits[1][w]*m_weightSliders[0][w]->value(); + float end = m_weightLimits[0][w] + m_weightLimits[1][w]*m_weightSliders[1][w]->value(); + m_weights[w] = passageNeg*start + passage*end; + } + float deltaH = m_weights[ElevationStep]; + float range = m_weights[Range]; optimizeElevations(0, 0, width-1, height-1, deltaH, range); if (m_stopProcessing) { @@ -855,18 +901,21 @@ void tcElevationOptimization::sampleStdDeviations() float stdDevVoid = 0.0f; float stdDevNonVoid = 0.0f; float stdDevAll = 0.0f; - int index = 0; + int total = 0; int voids = 0; - for (int j = 0; j < height; ++j) + int marginx = height>>3; + int marginy = width>>3; + for (int j = marginy; j < height-marginy; ++j) { - for (int i = 0; i < width; ++i) + for (int i = marginx; i < width-marginx; ++i) { - float altitude = m_data->buffer()[index].accurateElevation; + int index = j*width + i; + float altitude = m_datum->buffer()[index].accurateElevation; float square = m_elevationData->buffer()[index] - altitude; square *= square; stdDevAll += square; - if (m_data->buffer()[index].reliability) + if (m_datum->buffer()[index].reliability) { stdDevNonVoid += square; } @@ -875,16 +924,18 @@ void tcElevationOptimization::sampleStdDeviations() stdDevVoid += square; ++voids; } - ++index; + ++total; } } - stdDevNonVoid = sqrt(stdDevNonVoid / (index-voids)); + stdDevNonVoid = sqrt(stdDevNonVoid / (total-voids)); stdDevVoid = sqrt(stdDevVoid / voids); - stdDevAll = sqrt(stdDevAll / index); + stdDevAll = sqrt(stdDevAll / total); emit statistics(tr("Overall std deviation: %1 m\n" - "Void std deviation: %2 m\n" - "Non-void std deviation: %3 m") - .arg(stdDevAll).arg(stdDevVoid).arg(stdDevNonVoid)); + "Void std deviation: %2 (%4\%) m\n" + "Non-void std deviation: %3 (%5\%) m") + .arg(stdDevAll) + .arg(stdDevVoid).arg(stdDevNonVoid) + .arg(100.0f*voids/total).arg(100.0f*(total-voids)/total)); //m_configWidget->update(); if (m_snapshotReadFile.isEmpty() || m_snapshotReading) @@ -914,7 +965,13 @@ void tcElevationOptimization::sampleStdDeviations() { for (int i = 0; i < width; ++i) { - stream << m_elevationData->buffer()[index]; + uint16_t alt = floor(0.5f + m_elevationData->buffer()[index]); + uint16_t orig = floor(0.5f + m_datum->buffer()[index].originalElevation); + if (!m_datum->buffer()[index].reliability) + { + alt |= 0x8000; + } + stream << alt << orig; ++index; } } @@ -955,6 +1012,67 @@ void tcElevationOptimization::saveStdDeviations() } } +/// Save cross section data. +void tcElevationOptimization::saveCrossSections() +{ + if (!m_crossSectioned) + { + QMessageBox::warning(m_configWidget, + tr("Cross section save failed."), + tr("You need to select a cross section with middle mouse button drag.")); + return; + } + + QString filename = QFileDialog::getSaveFileName(m_configWidget, + tr("Save cross section data"), + QString(), + tr("Crosssection Files (*.xsec)")); + + if (!filename.isEmpty()) + { + // Open the file ready to write the data + QFile file(filename); + if (!file.open(QIODevice::WriteOnly | QIODevice::Text)) + { + QMessageBox::warning(m_configWidget, + tr("Open failed."), + tr("Could not open '%1' for writing.").arg(filename)); + return; + } + + QTextStream out(&file); + + int resolution = m_crossResolution->value(); + int frequency = m_crossFrequency->value(); + int snaps = m_snapshotSlider->maximum(); + tcGeo delta = (m_crossSection[1] - m_crossSection[0]) / (resolution-1); + // save distances + { + float deltaDist = tcGeo::angleBetween(m_crossSection[1], m_crossSection[0])/(resolution-1)*6378.137e3; + float distance = 0.0f; + for (int p = 0; p < resolution; ++p) + { + out << distance << ((p == resolution-1) ? "\n" : "\t"); + distance += deltaDist; + } + } + for (int snap = 1; snap <= snaps; snap += frequency) + { + // Load snapshot + loadSnapShot(snap); + + tcGeo tmp = m_crossSection[0]; + for (int p = 0; p < resolution; ++p) + { + float elev = m_elevationData->elevationAt(tmp); + out << elev << ((p == resolution-1) ? "\n" : "\t"); + tmp.setLon(tmp.lon() + delta.lon()); + tmp.setLat(tmp.lat() + delta.lat()); + } + } + } +} + /// Invalidate and redraw. void tcElevationOptimization::invalidateAndRedrawSlot() { @@ -965,6 +1083,8 @@ void tcElevationOptimization::invalidateAndRedrawSlot() /// Initialise snapshotting. void tcElevationOptimization::initSnapShotting() { + m_snapshotReadFile = QString(); + m_snapshotSlider->setEnabled(false); m_snapshotFile = QFileDialog::getSaveFileName(m_configWidget, tr("Save snapshotting config file"), QString(), @@ -1063,9 +1183,108 @@ void tcElevationOptimization::loadSnapShots() } } +/// Update weight labels. +void tcElevationOptimization::updateWeightLabels() +{ + for (int i = 0; i < WeightsMax; ++i) + { + m_weightLabels[i]->setText(QString("%1 -> %2") + .arg(m_weightLimits[0][i] + m_weightLimits[1][i]*m_weightSliders[0][i]->value()) + .arg(m_weightLimits[0][i] + m_weightLimits[1][i]*m_weightSliders[1][i]->value()) + ); + } +} + +/// Save weights. +void tcElevationOptimization::saveWeights() +{ + QString filename = QFileDialog::getSaveFileName(m_configWidget, + tr("Save weights"), + QString(), + tr("Elevation optimization weights (*.eow)")); + + if (!filename.isEmpty()) + { + // Open the file ready to write the data + QFile file(filename); + if (!file.open(QIODevice::WriteOnly | QIODevice::Text)) + { + QMessageBox::warning(m_configWidget, + tr("Open failed."), + tr("Could not open '%1' for writing.").arg(filename)); + return; + } + + QTextStream stream(&file); + for (int w = 0; w < WeightsMax; ++w) + { + float start = m_weightLimits[0][w] + m_weightLimits[1][w]*m_weightSliders[0][w]->value(); + float end = m_weightLimits[0][w] + m_weightLimits[1][w]*m_weightSliders[1][w]->value(); + stream << start << "\t" << end << "\n"; + } + } +} + +/// Load weights. +void tcElevationOptimization::loadWeights() +{ + QString filename = QFileDialog::getOpenFileName(m_configWidget, + tr("Load weights"), + QString(), + tr("Elevation optimization weights (*.eow)")); + + if (!filename.isEmpty()) + { + // Open the file ready to write the data + QFile file(filename); + if (!file.open(QIODevice::ReadOnly | QIODevice::Text)) + { + QMessageBox::warning(m_configWidget, + tr("Load failed."), + tr("Could not open '%1' for reading.").arg(filename)); + return; + } + + QTextStream stream(&file); + float starts[WeightsMax], ends[WeightsMax]; + for (int w = 0; w < WeightsMax; ++w) + { + stream >> starts[w] >> ends[w]; + if (stream.atEnd()) + { + if (w == WeightsMax-1) + { + // Go back to shading shifting the values + while (w > PhotometricStereo) + { + starts[w] = starts[w-1]; + ends[w] = ends[w-1]; + --w; + } + starts[w] = ends[w] = 0.0f; + break; + } + else + { + starts[w] = ends[w] = 0.0f; + QMessageBox::warning(m_configWidget, + tr("Panic!"), + tr("Not enough weights")); + } + } + } + for (int w = 0; w < WeightsMax; ++w) + { + m_weightSliders[0][w]->setValue((starts[w] - m_weightLimits[0][w]) / m_weightLimits[1][w]); + m_weightSliders[1][w]->setValue((ends[w] - m_weightLimits[0][w]) / m_weightLimits[1][w]); + } + } +} + /// Load a single snapshot. void tcElevationOptimization::loadSnapShot(int id) { + QString counter = "-"; // Load a snapshot if (!m_snapshotReadFile.isEmpty()) { @@ -1090,7 +1309,13 @@ void tcElevationOptimization::loadSnapShot(int id) { for (int i = 0; i < width; ++i) { - stream >> m_elevationData->buffer()[index]; + uint16_t alt; + uint16_t orig; + stream >> alt >> orig; + m_datum->buffer()[index].reliability = ((0 != (alt & 0x8000)) ? 0 : 1); + m_datum->buffer()[index].originalElevation = orig; + alt &= 0x7FFF; + m_elevationData->buffer()[index] = alt; ++index; } } @@ -1100,9 +1325,11 @@ void tcElevationOptimization::loadSnapShot(int id) invalidate(); m_configWidget->requestRedraw(); } + counter = QString("%1").arg(id); } } } + m_snapshotCounter->setText(counter); } /* @@ -1146,14 +1373,14 @@ tcAbstractPixelData* tcElevationOptimization::loadPortion(double x1, double y1, } /// @todo What if its a different size because the portion has changed? - if (0 == m_elevationData || 0 == m_data || changed) + if (0 == m_elevationData || 0 == m_datum || changed) { delete m_elevationData; for (int i = 0; i < m_numDatasets; ++i) { - delete m_datum[i]; + delete m_data[i]; } - delete m_data; + delete m_datum; // Create a new pixel buffer, initialised to altitude @@ -1163,10 +1390,10 @@ tcAbstractPixelData* tcElevationOptimization::loadPortion(double x1, double y1, tcGeo geoDiagonal = geoAt(maths::Vector<2,float>(x1 + (x2-x1)/(width-1), y1 + (y2-y1)/(height-1))) - geoBase; m_elevationData->setGeoToPixels(tcAffineTransform2(geoBase, geoDiagonal).inverse()); - m_data = new tcTypedPixelData(width, height); + m_datum = new tcTypedPixelData(width, height); for (int i = 0; i < m_numDatasets; ++i) { - m_datum[i] = new tcTypedPixelData(width, height); + m_data[i] = new tcTypedPixelData(width, height); } int maxWidthHeight = qMax(width, height); @@ -1184,7 +1411,7 @@ tcAbstractPixelData* tcElevationOptimization::loadPortion(double x1, double y1, y1 + (y2-y1)*j / (height - 1)); tcGeo geoCoord = geoAt(coord); - PerPixelCache* cache = &m_data->buffer()[index]; + PerPixelCache* cache = &m_datum->buffer()[index]; // Find resolution maths::Vector<2,float> coordx (x1 + (x2-x1)*(i+1) / (width - 1), @@ -1195,7 +1422,7 @@ tcAbstractPixelData* tcElevationOptimization::loadPortion(double x1, double y1, tcGeo geoCoordy = geoAt(coordy) - geoCoord; // Find approx float res = geoCoordx.angle(); - cache->resolution[0] = res; + cache->resolution[0] = res*cos(geoCoord.lat()); cache->resolution[1] = geoCoordy.angle(); cache->resolution *= 6378.137e3; @@ -1213,9 +1440,10 @@ tcAbstractPixelData* tcElevationOptimization::loadPortion(double x1, double y1, tcGeoImageData* imagery = m_datasets[ds]; maths::Vector<3,float> light = lightDirectionAt(geoCoord, imagery); - PerDatasetPixelCache* dsCache = &m_datum[ds]->buffer()[index]; + PerDatasetPixelCache* dsCache = &m_data[ds]->buffer()[index]; dsCache->light = light; - dsCache->sinLightElevation = light[2]/light.slice<0,2>().mag(); + dsCache->lxy = light.slice<0,2>().mag(); + dsCache->dlz_dlxy = light[2]/dsCache->lxy; tcGeo lightScaled(light[0]*res/sin(geoCoord.lat()), light[1]*res); @@ -1250,9 +1478,63 @@ tcAbstractPixelData* tcElevationOptimization::loadPortion(double x1, double y1, } dsCache->shadowTransition[depthDirection][0] = transitionI; dsCache->shadowTransition[depthDirection][1] = transitionJ; - dsCache->shadowTransitionDistance[depthDirection] = (pos-geoCoord).angle()*6378.137e3; + dsCache->shadowTransitionDistance[depthDirection] = tcGeo::angleBetween(pos,geoCoord)*6378.137e3; } } + + { + /* + * matrix L={Light row vec of shading channel i;..} + * invert L + * multiply by shading values for shading channels + * should be mostly about the same magnitude + * normalize + */ + // Estimate a normal from the shading channels + maths::Vector<3,float> totalN(0.0f, 0.0f, 0.01f); + + int ds[3]; + for (ds[0] = 0; ds[0] < m_numDatasets; ++ds[0]) + { + if (0 == m_shadingData[ds[0]]) + { + continue; + } + for (ds[1] = ds[0]+1; ds[1] < m_numDatasets; ++ds[1]) + { + if (0 == m_shadingData[ds[1]]) + { + continue; + } + for (ds[2] = ds[1]+1; ds[2] < m_numDatasets; ++ds[2]) + { + if (0 == m_shadingData[ds[2]]) + { + continue; + } + maths::Matrix<3,float> L(m_data[ds[0]]->buffer()[index].light, + m_data[ds[1]]->buffer()[index].light, + m_data[ds[2]]->buffer()[index].light); + maths::Vector<3,float> S; + bool singular = false; + L.invert(&singular); + if (singular) + { + continue; + } + for (int dsi = 0; dsi < 3; ++dsi) + { + S[dsi] = m_shadingData[ds[dsi]]->sampleFloat((float)i/(width-1), (float)j/(height-1)); + } + maths::Vector<3,float> N = S * L; + N.normalize(); + totalN += N; + } + } + } + cache->normalEstimate = totalN; + cache->normalEstimate.normalize(); + } } } } @@ -1264,6 +1546,9 @@ tcAbstractPixelData* tcElevationOptimization::loadPortion(double x1, double y1, dem()->updateFromElevationData(m_elevationData); } + // Don't bother outputting anything + //return new tcTypedPixelData(1,1); + // Downscale the elevation for viewing startProcessing("Downscaling for viewing"); tcPixelData* result = new tcPixelData(width, height); @@ -1276,72 +1561,143 @@ tcAbstractPixelData* tcElevationOptimization::loadPortion(double x1, double y1, progress((float)j/(height-1)); for (int i = 0; i < width; ++i) { - PerDatasetPixelCache* dsCache = &m_datum[0]->buffer()[index]; - bool isShadowed = (m_shadowData[0]->sampleFloat((float)i/(width-1), (float)j/(height-1)) < 0.5f); - bool isShadowEntrance = isShadowed && - dsCache->shadowTransition[TransitionEntrance][0] == i && - dsCache->shadowTransition[TransitionEntrance][1] == j; - bool isShadowExit = isShadowed && - dsCache->shadowTransition[TransitionExit][0] == i && - dsCache->shadowTransition[TransitionExit][1] == j; - int xdist = abs(dsCache->shadowTransition[TransitionEntrance][0] - i); - int ydist = abs(dsCache->shadowTransition[TransitionEntrance][1] - j); - bool bothIn = isShadowed && (dsCache->shadowTransition[TransitionExit][0] >= 0 && + //* + PerPixelCache* cache = &m_datum->buffer()[index]; + + result->buffer()[index] = 255 * cache->reliability; + ++index; + continue; + +#if 0 + if (0) + { + float dh_dx = 0.0f; + float dh_dy = 0.0f; + if (i > 0 && i < width-1) + { + float hp1 = m_elevationData->buffer()[index+1]; + float hm1 = m_elevationData->buffer()[index-1]; + dh_dx = (hp1-hm1)/2 / cache->resolution[0]; + } + if (j > 0 && j < height-1) + { + float hp1 = m_elevationData->buffer()[index+width]; + float hm1 = m_elevationData->buffer()[index-width]; + dh_dy = (hp1-hm1)/2 / cache->resolution[1]; + } + // (0.0f, -1.0f, dh_dy) x (1.0f, 0.0f, dh_dx) + maths::Vector<3,float> norm(-dh_dx, dh_dy, 1.0f); + norm.normalize(); + result->buffer()[index] = norm; + } +#endif + +#if 0 + for (int ds = 0; ds < m_numDatasets; ++ds) + { + result->buffer()[index] = 0; + PerDatasetPixelCache* dsCache = &m_data[ds]->buffer()[index]; + + float hwm1 = m_elevationData->sampleFloat( + ((float)i + dsCache->light[0]/cache->resolution[0])/(width-1), + ((float)j - dsCache->light[1]/cache->resolution[1])/(height-1) + ); + float hwp1 = m_elevationData->sampleFloat( + ((float)i - dsCache->light[0]/cache->resolution[0])/(width-1), + ((float)j + dsCache->light[1]/cache->resolution[1])/(height-1) + ); + float dh_dw = (hwm1 - hwp1)/(2*dsCache->lxy); + float facingness = dsCache->dlz_dlxy-dh_dw; +// if (facingness < 0.0f) +// result->buffer()[index] = 128.0f; + + result->buffer()[index] = 128.0f * facingness*facingness; + +result->buffer()[index] = 0; + + bool isShadowed = (m_shadowData[ds]->sampleFloat((float)i/(width-1), (float)j/(height-1)) < 0.5f); + bool isShadowEntrance = isShadowed && + dsCache->shadowTransition[TransitionEntrance][0] == i && + dsCache->shadowTransition[TransitionEntrance][1] == j; + bool isShadowExit = isShadowed && + dsCache->shadowTransition[TransitionExit][0] == i && + dsCache->shadowTransition[TransitionExit][1] == j; + float elevation = m_elevationData->buffer()[index]; + result->buffer()[index] = (elevation - 2500.0f)/1500.0f*255.0f; + + if (false && isShadowed) + { + // maximum elevation is between elevation at entrance and elevation at exit + if (m_weights[BelowShadowLine] != 0.0f && + dsCache->shadowTransition[TransitionExit][0] >= 0 && dsCache->shadowTransition[TransitionExit][0] < width && dsCache->shadowTransition[TransitionExit][1] >= 0 && dsCache->shadowTransition[TransitionExit][1] < height && dsCache->shadowTransition[TransitionEntrance][0] >= 0 && dsCache->shadowTransition[TransitionEntrance][0] < width && dsCache->shadowTransition[TransitionEntrance][1] >= 0 && - dsCache->shadowTransition[TransitionEntrance][1] < height); + dsCache->shadowTransition[TransitionEntrance][1] < height) + { + int exitIndex = width*dsCache->shadowTransition[TransitionExit][1] + dsCache->shadowTransition[TransitionExit][0]; + int entranceIndex = width*dsCache->shadowTransition[TransitionEntrance][1] + dsCache->shadowTransition[TransitionEntrance][0]; + float exitElevation = m_elevationData->buffer()[exitIndex]; + float entranceElevation = m_elevationData->buffer()[entranceIndex]; + + result->buffer()[index] = qMax(0.0f, elevation - (exitElevation * dsCache->shadowTransitionDistance[TransitionEntrance] + +entranceElevation * dsCache->shadowTransitionDistance[TransitionExit]) + / (dsCache->shadowTransitionDistance[TransitionEntrance] + +dsCache->shadowTransitionDistance[TransitionExit])); + } + } - //* - result->buffer()[index] = isShadowed ? qMin(255, 100 * (xdist + ydist)) : 255; - result->buffer()[index] = bothIn?255:0; - result->buffer()[index] = cost(i, j, i, j); - //*/ + result->buffer()[index] = 128.0f*m_shadingData[ds]->sampleFloat((float)i/(width-1), (float)j/(height-1)); #if 0 - /* - * matrix L={Light row vec of shading channel i;..} - * invert L - * multiply by shading values for shading channels - * should be mostly about the same magnitude - * normalize - */ - -#define SHAPEFROMSHADING_COMBOS 1 - int chans[SHAPEFROMSHADING_COMBOS][3] = { - /*{0,1,2}, - {0,1,3}, - {0,1,4}, - {0,2,3}, - {0,2,4}, - {0,3,4}, - {1,2,3}, - {1,2,4}, - {1,3,4},*/ - {2,3,4}, - }; - maths::Vector<3,float> totalN(0.0f); - int comb = 0; - for (; comb < SHAPEFROMSHADING_COMBOS; ++comb) - { - maths::Matrix<3,float> L(m_datum[chans[comb][0]]->buffer()[index].light, - m_datum[chans[comb][1]]->buffer()[index].light, - m_datum[chans[comb][2]]->buffer()[index].light); - maths::Vector<3,float> S; - L.invert(); - for (int ds = 0; ds < 3; ++ds) + int xdist = abs(dsCache->shadowTransition[TransitionEntrance][0] - i); + int ydist = abs(dsCache->shadowTransition[TransitionEntrance][1] - j); + bool bothIn = isShadowed && (dsCache->shadowTransition[TransitionExit][0] >= 0 && + dsCache->shadowTransition[TransitionExit][0] < width && + dsCache->shadowTransition[TransitionExit][1] >= 0 && + dsCache->shadowTransition[TransitionExit][1] < height && + dsCache->shadowTransition[TransitionEntrance][0] >= 0 && + dsCache->shadowTransition[TransitionEntrance][0] < width && + dsCache->shadowTransition[TransitionEntrance][1] >= 0 && + dsCache->shadowTransition[TransitionEntrance][1] < height); + + if (isShadowEntrance && (dsCache->shadowTransition[TransitionExit][0] >= 0 && + dsCache->shadowTransition[TransitionExit][0] < width && + dsCache->shadowTransition[TransitionExit][1] >= 0 && + dsCache->shadowTransition[TransitionExit][1] < height)) + { + int exitIndex = width*dsCache->shadowTransition[TransitionExit][1] + dsCache->shadowTransition[TransitionExit][0]; + + float exitElevation = m_elevationData->buffer()[exitIndex]; + float diff = exitElevation + dsCache->shadowTransitionDistance[TransitionExit]*dsCache->dlz_dlxy - elevation; + result->buffer()[index] += diff/5000.0f*255.0f; + } + else { - S[ds] = m_shadingData[chans[comb][ds]]->sampleFloat((float)i/(width-1), (float)j/(height-1)); + result->buffer()[index] = 0; } - maths::Vector<3,float> N = S * L; - N.normalize(); - totalN += N; +#endif } - result->buffer()[index] = totalN / comb; - result->buffer()[index].normalize(); + // */ + + //result->buffer()[index] = isShadowed ? qMin(255, 100 * (xdist + ydist)) : 255; + /*result->buffer()[index] = (isShadowed && (dsCache->shadowTransition[TransitionExit][0] >= 0 && + dsCache->shadowTransition[TransitionExit][0] < width && + dsCache->shadowTransition[TransitionExit][1] >= 0 && + dsCache->shadowTransition[TransitionExit][1] < height))?255:0;*/ + //float mycostChange = 100000.0f*255.0f*costChange(i, j, +4.0f, 0); + //result->buffer()[index] = 128.0f-mycostChange; + //result->buffer()[index] = m_datum->buffer()[index].reliability*255; + +#else + +#if 0 + result->buffer()[index] = cache->normalEstimate; +#endif + #endif ++index; @@ -1363,6 +1719,7 @@ T MySqr(T t) * Private functions */ +#include /// Calculate the cost for a set of pixels. float tcElevationOptimization::cost(int x1, int y1, int x2, int y2) const { @@ -1398,10 +1755,14 @@ float tcElevationOptimization::cost(int x1, int y1, int x2, int y2) const int countLitFacing = 0; float costShading = 0.0f; int countShading = 0; + float costStereo = 0.0f; + int countStereo = 0; float costBelowShadowLine = 0.0f; int countBelowShadowLine = 0; - float costTransitionElev = 0.0f; - int countTransitionElev = 0; + float costTransitionElevTop = 0.0f; + int countTransitionElevTop = 0; + float costTransitionElevBottom = 0.0f; + int countTransitionElevBottom = 0; float costTransitionTangential = 0.0f; int countTransitionTangential = 0; float costTransitionCurvingDown = 0.0f; @@ -1420,7 +1781,7 @@ float tcElevationOptimization::cost(int x1, int y1, int x2, int y2) const // Basic data retrieval float elevation = *elevationPtr; - PerPixelCache* cache = &m_data->buffer()[index]; + PerPixelCache* cache = &m_datum->buffer()[index]; if (m_voidsOnly && cache->reliability == 0) { continue; @@ -1448,22 +1809,49 @@ float tcElevationOptimization::cost(int x1, int y1, int x2, int y2) const } // Minimise variance from original value - costVariance += cache->reliability * (elevation-cache->originalElevation)*(elevation-cache->originalElevation); - ++countVariance; + if (m_weights[Variance] != 0.0f) + { + // treat borderline pixels as reliable + float reliability = cache->reliability; + if (i < 5 || i >= width-5 || j < 5 || j >= height-5) + { + reliability = 10.0f; + } + costVariance += reliability * (elevation-cache->originalElevation)*(elevation-cache->originalElevation); + ++countVariance; + } // Minimise roughness - costRoughness += (d2h_dx2*d2h_dx2 + d2h_dy2*d2h_dy2); - ++countRoughness; + if (m_weights[Roughness] != 0.0f) + { + costRoughness += (d2h_dx2*d2h_dx2 + d2h_dy2*d2h_dy2); + ++countRoughness; + } // Also minimise slope // This works really well at getting a nice slope - costSteepness += (dh_dx*dh_dx + dh_dy*dh_dy); - ++countSteepness; + if (m_weights[Steepness] != 0.0f) + { + costSteepness += (dh_dx*dh_dx + dh_dy*dh_dy); + ++countSteepness; + } + + // Photometric stereo + if (cache->normalEstimate[2] > 0.2f && m_weights[PhotometricStereo] != 0.0f) + { + // (0.0f, -1.0f, dh_dy) x (1.0f, 0.0f, dh_dx) + maths::Vector<3,float> norm(-dh_dx, dh_dy, 1.0f); + norm.normalize(); + costStereo += MySqr(acos(cache->normalEstimate*norm)); + //costStereo += MySqr(dh_dx - cache->normalEstimate[0]/cache->normalEstimate[2]); + //costStereo += MySqr(dh_dy - cache->normalEstimate[1]/cache->normalEstimate[2]); + ++countStereo; + } // Per dataset for (int ds = 0; ds < m_numDatasets; ++ds) { - PerDatasetPixelCache* dsCache = &m_datum[ds]->buffer()[index]; + PerDatasetPixelCache* dsCache = &m_data[ds]->buffer()[index]; bool isShadowed = (m_shadowData[ds]->sampleFloat((float)i/(width-1), (float)j/(height-1)) < 0.5f); bool isShadowEntrance = isShadowed && @@ -1475,18 +1863,24 @@ float tcElevationOptimization::cost(int x1, int y1, int x2, int y2) const // Calculate measures relating to the light direction // Gradient - float hwm1 = m_elevationData->sampleFloat(((float)i - dsCache->light[0])/(width-1), - ((float)j - dsCache->light[1])/(height-1)); - float hwp1 = m_elevationData->sampleFloat(((float)i + dsCache->light[0])/(width-1), - ((float)j + dsCache->light[1])/(height-1)); - /// @todo Ensure units are right (this is in terms of h metres, w pixels) - float dh_dw = (hwm1 - hwp1)/2 / cache->resolution[0]; // 30 metre resolution? + // Two sample points are 2*dsCache->lxy apart + float hwm1 = m_elevationData->sampleFloat( + ((float)i + dsCache->light[0]/cache->resolution[0])/(width-1), + ((float)j - dsCache->light[1]/cache->resolution[1])/(height-1) + ); + float hwp1 = m_elevationData->sampleFloat( + ((float)i - dsCache->light[0]/cache->resolution[0])/(width-1), + ((float)j + dsCache->light[1]/cache->resolution[1])/(height-1) + ); + float dh_dw = (hwm1 - hwp1)/(2*dsCache->lxy); float d2h_dw2 = (hwm1+hwp1)/2 - elevation; + float facingness = dsCache->dlz_dlxy-dh_dw; if (isShadowed) { // maximum elevation is between elevation at entrance and elevation at exit - if (dsCache->shadowTransition[TransitionExit][0] >= 0 && + if (m_weights[BelowShadowLine] != 0.0f && + dsCache->shadowTransition[TransitionExit][0] >= 0 && dsCache->shadowTransition[TransitionExit][0] < width && dsCache->shadowTransition[TransitionExit][1] >= 0 && dsCache->shadowTransition[TransitionExit][1] < height && @@ -1509,46 +1903,57 @@ float tcElevationOptimization::cost(int x1, int y1, int x2, int y2) const } else { - float facingness = -dh_dw+dsCache->light[0]; - if (facingness < 0.0f) + if (facingness < 0.0f && m_weights[LitFacing] != 0.0f) { // Lit pixels must face light costLitFacing += facingness*facingness; ++countLitFacing; } // Simple shape from shading - if (0 != m_shadingData[ds]) + if (0 != m_shadingData[ds] && m_weights[Shading] != 0.0f) { float intensity = m_shadingData[ds]->sampleFloat((float)i/(width-1), (float)j/(height-1)); // intensity = cos(theta) = L.N - maths::Vector<3,float> xTan(1.0f, 0.0f, dh_dx); - maths::Vector<3,float> yTan(0.0f, 1.0f, dh_dy); - costShading += MySqr(acos(qMin(1.0f, intensity)) - acos(qMin(1.0f, dsCache->light*maths::cross(xTan, yTan).normalized()))); + // (0.0f, -1.0f, dh_dy) x (1.0f, 0.0f, dh_dx) + maths::Vector<3,float> norm(-dh_dx, dh_dy, 1.0f); + norm.normalize(); + costShading += MySqr(acos(qMin(1.0f, intensity)) - acos(qMin(1.0f, dsCache->light*norm))); ++countShading; } } if (isShadowEntrance) { - if (dsCache->shadowTransition[TransitionExit][0] >= 0 && + if (m_weights[TransitionElevTop] != 0.0f && + dsCache->shadowTransition[TransitionExit][0] >= 0 && dsCache->shadowTransition[TransitionExit][0] < width && dsCache->shadowTransition[TransitionExit][1] >= 0 && dsCache->shadowTransition[TransitionExit][1] < height) { int exitIndex = width*dsCache->shadowTransition[TransitionExit][1] + dsCache->shadowTransition[TransitionExit][0]; float exitElevation = m_elevationData->buffer()[exitIndex]; - costTransitionElev += MySqr(exitElevation + dsCache->shadowTransitionDistance[TransitionExit]*dsCache->sinLightElevation - elevation); - ++countTransitionElev; + //if (dsCache->shadowTransitionDistance[TransitionExit]*dsCache->dlz_dlxy > 100.0f) + { + costTransitionElevTop += MySqr(exitElevation + dsCache->shadowTransitionDistance[TransitionExit]*dsCache->dlz_dlxy - elevation); + ++countTransitionElevTop; + } } // gradient tangential to sun direction - costTransitionTangential += MySqr(dh_dw - dsCache->sinLightElevation); - ++countTransitionTangential; + if (m_weights[TransitionTangential] != 0.0f) + { + costTransitionTangential += MySqr(facingness); + ++countTransitionTangential; + } // second derivitive should be negative at entrance - costTransitionCurvingDown += MySqr(qMax(0.0f, d2h_dw2)); - ++countTransitionCurvingDown; + if (m_weights[TransitionCurvingDown] != 0.0f) + { + costTransitionCurvingDown += MySqr(qMax(0.0f, d2h_dw2)); + ++countTransitionCurvingDown; + } } - if (0 && isShadowExit) + if (isShadowExit) { - if (dsCache->shadowTransition[TransitionEntrance][0] >= 0 && + if (m_weights[TransitionElevBottom] != 0.0f && + dsCache->shadowTransition[TransitionEntrance][0] >= 0 && dsCache->shadowTransition[TransitionEntrance][0] < width && dsCache->shadowTransition[TransitionEntrance][1] >= 0 && dsCache->shadowTransition[TransitionEntrance][1] < height) @@ -1557,23 +1962,25 @@ float tcElevationOptimization::cost(int x1, int y1, int x2, int y2) const float entranceElevation = m_elevationData->buffer()[entranceIndex]; // this has less effect, its usually "more" right /// @todo Smooth high reliability into low reliability - costTransitionElev += 0.1f * MySqr(entranceElevation - dsCache->shadowTransitionDistance[TransitionEntrance]*dsCache->sinLightElevation - elevation); - ++countTransitionElev; + costTransitionElevBottom += 1.0f * MySqr(entranceElevation - dsCache->shadowTransitionDistance[TransitionEntrance]*dsCache->dlz_dlxy - elevation); + ++countTransitionElevBottom; } } } } } - return m_weights.variance*costVariance//countVariance - + m_weights.roughness*costRoughness//countRoughness - + m_weights.steepness*costSteepness//countSteepness - + m_weights.litFacing*costLitFacing//countLitFacing - + m_weights.shading*costShading//countShading - + m_weights.belowShadowLine*costBelowShadowLine//countBelowShadowLine - + m_weights.transitionElev*costTransitionElev//countTransitionElev - + m_weights.transitionTangential*costTransitionTangential//countTransitionTangential - + m_weights.transitionCurvingDown*costTransitionCurvingDown//countTransitionCurvingDown + return m_weights[Variance]*costVariance//countVariance + + m_weights[Roughness]*costRoughness//countRoughness + + m_weights[Steepness]*costSteepness//countSteepness + + m_weights[LitFacing]*costLitFacing//countLitFacing + + m_weights[Shading]*costShading//countShading + + m_weights[PhotometricStereo]*costStereo//countStereo + + m_weights[BelowShadowLine]*costBelowShadowLine//countBelowShadowLine + + m_weights[TransitionElevTop]*costTransitionElevTop//countTransitionElevTop + + m_weights[TransitionElevBottom]*costTransitionElevBottom//countTransitionElevBottom + + m_weights[TransitionTangential]*costTransitionTangential//countTransitionTangential + + m_weights[TransitionCurvingDown]*costTransitionCurvingDown//countTransitionCurvingDown ; } @@ -1610,10 +2017,10 @@ float tcElevationOptimization::costChange(int x, int y, float dh, int range) { int dj = j - range; float* elev = pixelElevation(x+di, y+dj); - if (0 != elev && (!m_voidsOnly || 0 != m_data->buffer()[j*m_data->width() + i].reliability)) + if (0 != elev && (!m_voidsOnly || 0 != m_datum->buffer()[j*m_datum->width() + i].reliability)) { oldElevations[i][j] = *elev; - *elev += dh*elevationFalloffFunction(float(di*di + dj*dj) / ((range-1)*(range-1))); + *elev += dh*elevationFalloffFunction(float(di*di + dj*dj) / ((range+1)*(range+1))); } } } @@ -1625,7 +2032,7 @@ float tcElevationOptimization::costChange(int x, int y, float dh, int range) { int dj = j - range; float* elev = pixelElevation(x+di, y+dj); - if (0 != elev && (!m_voidsOnly || 0 != m_data->buffer()[j*m_data->width() + i].reliability)) + if (0 != elev && (!m_voidsOnly || 0 != m_datum->buffer()[j*m_datum->width() + i].reliability)) { *elev = oldElevations[i][j]; } @@ -1645,9 +2052,9 @@ void tcElevationOptimization::optimizeElevations(int x1, int y1, int x2, int y2, { return; } - int index = j*m_data->width() + i; + int index = j*m_datum->width() + i; //while (true) - if (!m_voidsOnly || 0 != m_data->buffer()[j*m_data->width() + i].reliability) + if (!m_voidsOnly || 0 != m_datum->buffer()[j*m_datum->width() + i].reliability) { // Find the cost of adjusting the elevation by dh in each direction float costInc = costChange(i,j, dh, range); diff --git a/geo/tcElevationOptimization.h b/geo/tcElevationOptimization.h index dc79766..b898850 100644 --- a/geo/tcElevationOptimization.h +++ b/geo/tcElevationOptimization.h @@ -42,6 +42,7 @@ class QSpinBox; class QSlider; class QThread; class QPushButton; +class QLabel; /// Optimized elevation. class tcElevationOptimization : public tcChannelDem @@ -114,6 +115,9 @@ class tcElevationOptimization : public tcChannelDem /// Sample standard deviations. void saveStdDeviations(); + /// Save cross section data. + void saveCrossSections(); + /// Invalidate and redraw. void invalidateAndRedrawSlot(); @@ -126,6 +130,15 @@ class tcElevationOptimization : public tcChannelDem /// Load a single snapshot. void loadSnapShot(int id); + /// Update weight labels. + void updateWeightLabels(); + + /// Save weights. + void saveWeights(); + + /// Load weights. + void loadWeights(); + protected: /* @@ -190,10 +203,10 @@ class tcElevationOptimization : public tcChannelDem maths::Vector<2,int> shadowTransition[2]; /// Light direction. maths::Vector<3,float> light; - /// Length of light vector trace along ground. - float sinLightElevation; - /// How lit pixel is. - float lit; + /// Length of light direction in x y. + float lxy; + /// Gradiant of light direction. + float dlz_dlxy; /// Distance to shadow entrance and exit. float shadowTransitionDistance[2]; }; @@ -209,6 +222,8 @@ class tcElevationOptimization : public tcChannelDem unsigned char reliability; /// Accurate elevation. float accurateElevation; + /// Estimated normal. + maths::Vector<3,float> normalEstimate; /// Temporary values. union { float f; @@ -239,10 +254,10 @@ class tcElevationOptimization : public tcChannelDem tcElevationData* m_elevationData; /// Per-dataset-pixel data. - tcTypedPixelData** m_datum; + tcTypedPixelData** m_data; /// Per-pixel data. - tcTypedPixelData* m_data; + tcTypedPixelData* m_datum; /// Temporary shadow data. Reference >* m_shadowData; @@ -250,19 +265,35 @@ class tcElevationOptimization : public tcChannelDem /// Temporary shading data. Reference >* m_shadingData; - /// Weights. - struct + enum Weights { - float variance; - float roughness; - float steepness; - float litFacing; - float shading; - float belowShadowLine; - float transitionElev; - float transitionTangential; - float transitionCurvingDown; - } m_weights; + Variance, + Roughness, + Steepness, + LitFacing, + Shading, + PhotometricStereo, + BelowShadowLine, + TransitionElevTop, + TransitionElevBottom, + TransitionTangential, + TransitionCurvingDown, + ElevationStep, + Range, + WeightsMax + }; + + /// Weights. + float m_weights[WeightsMax]; + + /// Weight sliders. + QSlider* m_weightSliders[2][WeightsMax]; + + /// Weight labels. + QLabel* m_weightLabels[WeightsMax]; + + /// Weight limits. + float m_weightLimits[2][WeightsMax]; /// Update voids only. bool m_voidsOnly; @@ -292,12 +323,6 @@ class tcElevationOptimization : public tcChannelDem /// Number of iterations to perform. QSpinBox* m_spinIterations; - /// Slider for delta h. - QSlider* m_deltaHSlider[2]; - - /// Slider for range. - QSlider* m_rangeSlider[2]; - /// Plot of standard deviations. QwtPlot* m_plot; @@ -322,6 +347,12 @@ class tcElevationOptimization : public tcChannelDem /// Standard deviation data. QwtArray m_stdDevNonVoidData; + /// Cross sectioning resolution. + QSpinBox* m_crossResolution; + + /// Cross sectioning frequency. + QSpinBox* m_crossFrequency; + /// Snapshotting config file. QString m_snapshotFile; @@ -333,6 +364,9 @@ class tcElevationOptimization : public tcChannelDem /// Snapshot time slider. QSlider* m_snapshotSlider; + + /// Counter of snapshot frames. + QLabel* m_snapshotCounter; }; #endif diff --git a/geo/tcGeo.h b/geo/tcGeo.h index 2c35884..aac3b10 100644 --- a/geo/tcGeo.h +++ b/geo/tcGeo.h @@ -296,6 +296,13 @@ class tcGeo return sqrt(m_longitude*m_longitude + m_latitude*m_latitude); } + static double angleBetween(const tcGeo& lhs, const tcGeo& rhs) + { + //return acos((maths::Vector<3,double>)lhs* (maths::Vector<3,double>)rhs); + maths::Vector<3,double> x = cross((maths::Vector<3,double>)lhs, (maths::Vector<3,double>)rhs); + return asin(x.mag()); + } + private: /* diff --git a/geo/tcGlobe.cpp b/geo/tcGlobe.cpp index 13060aa..af184b9 100644 --- a/geo/tcGlobe.cpp +++ b/geo/tcGlobe.cpp @@ -70,23 +70,28 @@ tcGlobe::tcGlobe(double meanRadius) // high azimuth, low elevation datasets << new tcLandsatData("/home/james/cs/pro/data/LE71950282001307EDC00/", m_elevation); datasets << new tcLandsatData("/home/james/cs/pro/data/LE71950282002006EDC00/", m_elevation); - datasets << new tcLandsatData("/home/james/cs/pro/data/LE71950282002310EDC00/", m_elevation); + //datasets << new tcLandsatData("/home/james/cs/pro/data/LE71950282002310EDC00/", m_elevation); datasets << new tcLandsatData("/home/james/cs/pro/data/LE71950282000081EDC00/", m_elevation); - datasets << new tcLandsatData("/home/james/cs/pro/data/etp195r28_5t19900910/", m_elevation); + //datasets << new tcLandsatData("/home/james/cs/pro/data/etp195r28_5t19900910/", m_elevation); -#define NEW_DATASET(AZ, EL) \ +// Relative to 0N, 0E +#define NEW_DATASET(AZ, EL, TILT) \ + datasets << new tcCustomSunDirection(m_elevation+1, tcGeo((180.0-AZ)*M_PI/180, (EL)*M_PI/180), \ + tcGeo(7.0*M_PI/180, (TILT)*M_PI/180), true) +// Relative to our region of interest +#define NEW_DATASET_REL(AZ, EL) \ datasets << new tcCustomSunDirection(m_elevation+1, tcGeo((180.0-AZ)*M_PI/180, (EL)*M_PI/180), \ tcGeo(7.0*M_PI/180, 46.0*M_PI/180)) - /* - for (int i = 0; i < 360; i += 5) + + // Mustn't be co-planar + for (int i = -10; i <= 10; i += 20) { - for (int j = 30; j <= 30; j += 2) + for (int j = 30; j <= 150; j += 10) { - NEW_DATASET((double)i, (double)j); + //NEW_DATASET((double)90, (double)j, (double)i); } } - */ foreach (tcShadowClassifyingData* dataset, datasets) { diff --git a/geo/tcNegativeProduct.cpp b/geo/tcNegativeProduct.cpp index b12a6ba..8f49d15 100644 --- a/geo/tcNegativeProduct.cpp +++ b/geo/tcNegativeProduct.cpp @@ -32,6 +32,7 @@ #include #include #include +#include #include #include #include @@ -143,35 +144,42 @@ tcChannelConfigWidget* tcNegativeProduct::configWidget() m_powerLabels += labelPowerValue; } - QWidget* widgetOptimise = new QWidget(m_configWidget); - layout->addWidget(widgetOptimise); - QHBoxLayout* layoutOptimise = new QHBoxLayout(widgetOptimise); - m_comboClassifier = new QComboBox(widgetOptimise); - layoutOptimise->addWidget(m_comboClassifier); + QWidget* widgetOptimize = new QWidget(m_configWidget); + layout->addWidget(widgetOptimize); + QHBoxLayout* layoutOptimize = new QHBoxLayout(widgetOptimize); + m_comboClassifier = new QComboBox(widgetOptimize); + layoutOptimize->addWidget(m_comboClassifier); for (int i = 0; i < channelManager()->numChannels(); ++i) { tcChannel* channel = channelManager()->channel(i); m_comboClassifier->addItem(channel->name()); } - QPushButton* btnOptimiseThreshold = new QPushButton(tr("Optimise Threshold"), widgetOptimise); - layoutOptimise->addWidget(btnOptimiseThreshold); - connect(btnOptimiseThreshold, SIGNAL(clicked(bool)), this, SLOT(optimiseThreshold())); - QPushButton* btnSaveThresholdErrors = new QPushButton(tr("Save Threshold Errors"), widgetOptimise); - layoutOptimise->addWidget(btnSaveThresholdErrors); + QPushButton* btnCalculateError = new QPushButton(tr("Calculate error"), widgetOptimize); + layoutOptimize->addWidget(btnCalculateError); + connect(btnCalculateError, SIGNAL(clicked(bool)), this, SLOT(calculateError())); + QPushButton* btnOptimizeThreshold = new QPushButton(tr("Optimize Threshold"), widgetOptimize); + layoutOptimize->addWidget(btnOptimizeThreshold); + connect(btnOptimizeThreshold, SIGNAL(clicked(bool)), this, SLOT(optimizeThreshold())); + m_subsample = new QSpinBox(widgetOptimize); + layoutOptimize->addWidget(m_subsample); + m_subsample->setRange(1,20); + m_subsample->setValue(1); + QPushButton* btnSaveThresholdErrors = new QPushButton(tr("Save Threshold Errors"), widgetOptimize); + layoutOptimize->addWidget(btnSaveThresholdErrors); connect(btnSaveThresholdErrors, SIGNAL(clicked(bool)), this, SLOT(saveThresholdErrors())); - QWidget* widgetOptimiseParameters = new QWidget(m_configWidget); - layout->addWidget(widgetOptimiseParameters); - QHBoxLayout* layoutOptimiseParameters = new QHBoxLayout(widgetOptimiseParameters); - m_sliderLearningRate = new QSlider(Qt::Horizontal, widgetOptimiseParameters); - layoutOptimiseParameters->addWidget(m_sliderLearningRate); + QWidget* widgetOptimizeParameters = new QWidget(m_configWidget); + layout->addWidget(widgetOptimizeParameters); + QHBoxLayout* layoutOptimizeParameters = new QHBoxLayout(widgetOptimizeParameters); + m_sliderLearningRate = new QSlider(Qt::Horizontal, widgetOptimizeParameters); + layoutOptimizeParameters->addWidget(m_sliderLearningRate); m_sliderLearningRate->setRange(0, 5); m_sliderLearningRate->setValue(5); - QPushButton* btnOptimiseParams = new QPushButton(tr("Optimise Parameters"), widgetOptimiseParameters); - layoutOptimiseParameters->addWidget(btnOptimiseParams); - connect(btnOptimiseParams, SIGNAL(clicked(bool)), this, SLOT(optimiseParameters())); + QPushButton* btnOptimizeParams = new QPushButton(tr("Optimize Parameters"), widgetOptimizeParameters); + layoutOptimizeParameters->addWidget(btnOptimizeParams); + connect(btnOptimizeParams, SIGNAL(clicked(bool)), this, SLOT(optimizeParameters())); - m_labelError = new QLabel(tr("Error: Unknown"), widgetOptimise); + m_labelError = new QLabel(tr("Error: Unknown"), widgetOptimize); layout->addWidget(m_labelError); QWidget* widgetSaveLoad = new QWidget(m_configWidget); @@ -341,8 +349,8 @@ void tcNegativeProduct::changePower(int channel) m_configWidget->requestRedraw(); } -/// Optimise the threshold for the parameters using a classification channel. -void tcNegativeProduct::optimiseThreshold() +/// Optimize the threshold for the parameters using a classification channel. +void tcNegativeProduct::optimizeThreshold() { tcChannel* classifier = channelManager()->channel(m_comboClassifier->currentIndex()); if (0 != classifier) @@ -371,9 +379,10 @@ void tcNegativeProduct::optimiseThreshold() int height = data->height(); int pixels = width*height; int baseError = 0; - for (int j = 0; j < height; ++j) + int subsample = m_subsample->value(); + for (int j = 0; j < height; j += subsample) { - for (int i = 0; i < width; ++i) + for (int i = 0; i < width; i += subsample) { int index = j*width + i; float val = data->buffer()[index]; @@ -449,7 +458,7 @@ void tcNegativeProduct::optimiseThreshold() delete [] bins; m_optimalThresholdError = minError; - m_labelError->setText(tr("Error: %1% (%2/%3 pixels)").arg(100.0f*minError/pixels).arg(minError).arg(pixels)); + calculateError(); m_sliderThreshold->setValue(min + minBin); setThreshold(min + minBin); } @@ -474,7 +483,7 @@ void tcNegativeProduct::saveThresholdErrors() setThresholdEnabled(false); } - optimiseThreshold(); + optimizeThreshold(); m_saveThresholdError = QString(); if (previousThresholding = m_thresholdEnabled) @@ -483,8 +492,8 @@ void tcNegativeProduct::saveThresholdErrors() } } -/// Optimise the parameters using a classification channel. -void tcNegativeProduct::optimiseParameters() +/// Optimize the parameters using a classification channel. +void tcNegativeProduct::optimizeParameters() { tcChannel* classifier = channelManager()->channel(m_comboClassifier->currentIndex()); if (0 != classifier) @@ -507,19 +516,19 @@ void tcNegativeProduct::optimiseParameters() for (int c = 0; c < m_inputChannels.size(); ++c) { // Vary the sliders a little bit to see which direction is best - optimiseThreshold(); + optimizeThreshold(); progress(tr("learning rate %1, step: %2, channel: %3, %4").arg(learningRate).arg(i).arg(m_inputChannels[c]->name()).arg(m_labelError->text())); int baseError = m_optimalThresholdError; int power = m_powerSliders[c]->value(); m_powerSliders[c]->setValue(power-learningRate); changePower(c); - optimiseThreshold(); + optimizeThreshold(); int reduceError = m_optimalThresholdError; m_powerSliders[c]->setValue(power+learningRate); changePower(c); - optimiseThreshold(); + optimizeThreshold(); int increaseError = m_optimalThresholdError; if (baseError < increaseError && baseError < reduceError) @@ -547,7 +556,7 @@ void tcNegativeProduct::optimiseParameters() learningRate /= 2; perLRate *= 2; } - optimiseThreshold(); + optimizeThreshold(); if (previousThresholding = m_thresholdEnabled) { setThresholdEnabled(previousThresholding); @@ -556,6 +565,52 @@ void tcNegativeProduct::optimiseParameters() } } +/// Calculate errors. +void tcNegativeProduct::calculateError() +{ + m_labelError->setText(tr("Error: Unknown")); + tcChannel* classifier = channelManager()->channel(m_comboClassifier->currentIndex()); + if (0 != classifier) + { + Reference abstractData = portion(); + Reference classification = classifier->portion(); + if (0 != abstractData && 0 != classification) + { + Reference > data = dynamicCast*>(abstractData); + Q_ASSERT(0 != data); + + // indexed first by reference shadow map litness, then by negative product litness + int counts[2][2] = { {0, 0}, {0, 0} }; + + int width = data->width(); + int height = data->height(); + for (int j = 0; j < height; ++j) + { + for (int i = 0; i < width; ++i) + { + int index = j*width + i; + int val = (data->buffer()[index] > 0.5f ? 1 : 0); + int classified = (classification->sampleFloat((float)i/(width-1), (float)j/(height-1)) > 0.5f ? 1 : 0); + ++counts[classified][val]; + } + } + + int total = counts[0][0]+counts[0][1]+counts[1][0]+counts[1][1]; + int correct = counts[0][0]+counts[1][1]; + int incorrect = counts[0][1]+counts[1][0]; + m_labelError->setText(tr("Error: %1% (%2/%3 pixels) [%4,%5,%6,%7]") + .arg(100.0f*incorrect/total) + .arg(incorrect) + .arg(total) + .arg(counts[0][0]) + .arg(counts[0][1]) + .arg(counts[1][0]) + .arg(counts[1][1])); + } + } + +} + /// Save parameters to a file. void tcNegativeProduct::saveParams() { diff --git a/geo/tcNegativeProduct.h b/geo/tcNegativeProduct.h index 9d7ef4f..345e3a5 100644 --- a/geo/tcNegativeProduct.h +++ b/geo/tcNegativeProduct.h @@ -30,6 +30,7 @@ #include class QSlider; +class QSpinBox; class QLabel; class QComboBox; @@ -106,14 +107,17 @@ class tcNegativeProduct : public tcChannel /// Indicates that the power of a channel has changed. void changePower(int channel); - /// Optimise the threshold for the parameters using a classification channel. - void optimiseThreshold(); + /// Optimize the threshold for the parameters using a classification channel. + void optimizeThreshold(); /// Save threshold errors. void saveThresholdErrors(); - /// Optimise the parameters using a classification channel. - void optimiseParameters(); + /// Optimize the parameters using a classification channel. + void optimizeParameters(); + + /// Calculate errors. + void calculateError(); /// Save parameters to a file. void saveParams(); @@ -175,6 +179,9 @@ class tcNegativeProduct : public tcChannel /// Slider for log learning rate. QSlider* m_sliderLearningRate; + /// Spinbox for subsampling. + QSpinBox* m_subsample; + /// Label to show error. QLabel* m_labelError; diff --git a/geo/tcSrtmModel.cpp b/geo/tcSrtmModel.cpp index a43b96f..593f393 100644 --- a/geo/tcSrtmModel.cpp +++ b/geo/tcSrtmModel.cpp @@ -187,22 +187,10 @@ maths::Vector<3,float> tcSrtmModel::normalAt(const tcGeo& coord, bool corrected, } maths::Vector<2,float> res((radiusEarth+alt) * angRes.lon() * cos(coord.lat()), (radiusEarth+alt) * angRes.lat()); - return ( maths::Vector<3,float>( - -dx, - 0.0f, - res[0] - ).normalized() + - maths::Vector<3,float>( - 0.0f, - -dy, - res[1] - ).normalized() - ).normalized(); - /*return maths::Vector<3,float>( - -dx, - -dy, - res[0] + res[1] - ).normalized();*/ + // (1.0f, 0.0f, dx/res[0]) x (0.0f, -1.0f, dy/res[1]) + maths::Vector<3,float> norm(-dx/res[0], -dy/res[1], 1.0f); + norm.normalize(); + return norm; } } if (0 != accurate) diff --git a/maths/Matrix.cpp b/maths/Matrix.cpp index dc62f7e..ba2681a 100644 --- a/maths/Matrix.cpp +++ b/maths/Matrix.cpp @@ -66,7 +66,7 @@ namespace maths { // Invert the 3x3 matrix /// not written by me template <> - Matrix<3,float> & Matrix<3,float>::invert() + Matrix<3,float> & Matrix<3,float>::invert(bool* singular) { Matrix a = *this; Matrix b; @@ -77,6 +77,11 @@ namespace maths { unsigned int rowMax; // Points to max abs value row in this column unsigned int row; float tmp; + + if (0 != singular) + { + *singular = false; + } // Go through columns for (c=0; c<3; c++) { @@ -90,7 +95,13 @@ namespace maths { // If the max value here is 0, we can't invert. Return identity. if (a[rowMax][c] == 0.0F) + { + if (0 != singular) + { + *singular = true; + } return(identity()); + } // Swap row "rowMax" with row "c" for (cc=0; cc<3; cc++) @@ -131,7 +142,7 @@ namespace maths { // Invert the Matrix<4,float> template <> - Matrix<4,float> & Matrix<4,float>::invert() + Matrix<4,float> & Matrix<4,float>::invert(bool* singular) { Matrix<4,float> a(*this); Matrix<4,float> b; @@ -142,6 +153,11 @@ namespace maths { unsigned int rowMax; // Points to max abs value row in this column unsigned int row; float tmp; + + if (0 != singular) + { + *singular = false; + } // Go through columns for (c=0; c<4; c++) { @@ -154,8 +170,14 @@ namespace maths { } // If the max value here is 0, we can't invert. Return identity. - //if (a[rowMax][c] == 0.0F) - // return(identity()); + if (a[rowMax][c] == 0.0F) + { + if (0 != singular) + { + *singular = true; + } + return(identity()); + } // Swap row "rowMax" with row "c" for (cc=0; cc<4; cc++) { diff --git a/maths/Matrix.h b/maths/Matrix.h index b73c5fc..9597f06 100644 --- a/maths/Matrix.h +++ b/maths/Matrix.h @@ -366,13 +366,13 @@ namespace maths } /// Invert the matrix. - Matrix & invert(); + Matrix & invert(bool* singular = 0); /// Get the inverted matrix. - inline Matrix inverted() const + inline Matrix inverted(bool* singular = 0) const { Matrix ret = *this; - return ret.invert(); + return ret.invert(singular); } }; // ::maths::Matrix diff --git a/tcViewportWidget.cpp b/tcViewportWidget.cpp index 81ff984..3c3bc5d 100644 --- a/tcViewportWidget.cpp +++ b/tcViewportWidget.cpp @@ -571,7 +571,7 @@ void tcViewportWidget::paintGL() } tcGeo delta = ending - coord; - int nseg = 2 + fabs(delta.angle())*100*meanRadius / m_observer->range(); + int nseg = 2 + tcGeo::angleBetween(ending, coord)*100*meanRadius / m_observer->range(); double dlon = delta.lon() / (nseg-1); double dlat = delta.lat() / (nseg-1); @@ -666,7 +666,7 @@ void tcViewportWidget::mouseMoveEvent(QMouseEvent* event) { mouseDescription = m_mouseCoord.describe() + " " + tr("%1 m, range: %2 m") .arg(m_globe->altitudeAt(m_mouseCoord)) - .arg((m_mouseCoord-m_observer->focus()).angle() * m_globe->meanRadius()); + .arg(tcGeo::angleBetween(m_mouseCoord,m_observer->focus()) * m_globe->meanRadius()); emit mouseGeoChanged(m_mouseCoord); } else -- 2.11.4.GIT