//########################################################################## //# # //# CLOUDCOMPARE PLUGIN: q3DMASC # //# # //# This program is free software; you can redistribute it and/or modify # //# it under the terms of the GNU General Public License as published by # //# the Free Software Foundation; version 2 or later of the License. # //# # //# This program is distributed in the hope that it will be useful, # //# but WITHOUT ANY WARRANTY; without even the implied warranty of # //# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # //# GNU General Public License for more details. # //# # //# COPYRIGHT: Dimitri Lague / CNRS / UEB # //# # //########################################################################## #include "q3DMASCClassifier.h" //Local #include "ScalarFieldWrappers.h" #include "q3DMASCTools.h" //qCC_db #include #include #include #include //qPDALIO #include "../../../core/IO/qPDALIO/include/LASFields.h" //CCPluginAPI #include //Qt #include #include #include #include #include "qTrain3DMASCDialog.h" #include "confusionmatrix.h" #if defined(_OPENMP) #include #endif #if defined(CC_MAC_OS) || defined(CC_LINUX) #include #endif using namespace masc; Classifier::Classifier() { } bool Classifier::isValid() const { return (m_rtrees && m_rtrees->isClassifier() && m_rtrees->isTrained()); } static IScalarFieldWrapper::Shared GetSource(const Feature::Source& fs, const ccPointCloud* cloud) { IScalarFieldWrapper::Shared source(nullptr); switch (fs.type) { case Feature::Source::ScalarField: { assert(!fs.name.isEmpty()); int sfIdx = cloud->getScalarFieldIndexByName(qPrintable(fs.name)); if (sfIdx >= 0) { source.reset(new ScalarFieldWrapper(cloud->getScalarField(sfIdx))); } else { ccLog::Warning(QObject::tr("Internal error: unknown scalar field '%1'").arg(fs.name)); return IScalarFieldWrapper::Shared(nullptr); } } break; case Feature::Source::DimX: source.reset(new DimScalarFieldWrapper(cloud, DimScalarFieldWrapper::DimX)); break; case Feature::Source::DimY: source.reset(new DimScalarFieldWrapper(cloud, DimScalarFieldWrapper::DimY)); break; case Feature::Source::DimZ: source.reset(new DimScalarFieldWrapper(cloud, DimScalarFieldWrapper::DimZ)); break; case Feature::Source::Red: source.reset(new ColorScalarFieldWrapper(cloud, ColorScalarFieldWrapper::Red)); break; case Feature::Source::Green: source.reset(new ColorScalarFieldWrapper(cloud, ColorScalarFieldWrapper::Green)); break; case Feature::Source::Blue: source.reset(new ColorScalarFieldWrapper(cloud, ColorScalarFieldWrapper::Blue)); break; } return source; } bool Classifier::classify( const Feature::Source::Set& featureSources, ccPointCloud* cloud, QString& errorMessage, QWidget* parentWidget/*=nullptr*/, ccMainAppInterface* app/*nullptr*/ ) { if (!cloud) { assert(false); errorMessage = QObject::tr("Invalid input"); return false; } if (!isValid()) { errorMessage = QObject::tr("Invalid classifier"); return false; } if (featureSources.empty()) { errorMessage = QObject::tr("Training method called without any feature (source)?!"); return false; } // add a ccConfidence value if needed int cvConfidenceIdx = cloud->getScalarFieldIndexByName("Classification_confidence"); if (cvConfidenceIdx >= 0) // if the scalar field exists, delete it cloud->deleteScalarField(cvConfidenceIdx); cvConfidenceIdx = cloud->addScalarField("Classification_confidence"); ccScalarField* cvConfidenceSF = static_cast(cloud->getScalarField(cvConfidenceIdx)); //look for the classification field CCCoreLib::ScalarField* classificationSF = Tools::GetClassificationSF(cloud); ccScalarField* classifSFBackup = nullptr; if (classificationSF) //save classification field (if any) by renaming it "Classification_backup" { ccLog::Warning("Classification SF found: copy it in Classification_backup, a confusion matrix will be generated"); // delete Classification_backup field (if any) int sfIdx = cloud->getScalarFieldIndexByName("Classification_backup"); if (sfIdx >= 0) cloud->deleteScalarField(sfIdx); classificationSF->setName("Classification_backup"); // rename the classification field classifSFBackup = static_cast(classificationSF); } //create the classification SF ccScalarField* _classificationSF = new ccScalarField(LAS_FIELD_NAMES[LAS_CLASSIFICATION]); if (!_classificationSF->resizeSafe(cloud->size())) { _classificationSF->release(); errorMessage = QObject::tr("Not enough memory"); return false; } cloud->addScalarField(_classificationSF); classificationSF = _classificationSF; assert(classificationSF); classificationSF->fill(0); //0 = no classification? int sampleCount = static_cast(cloud->size()); int attributesPerSample = static_cast(featureSources.size()); ccLog::Print(QObject::tr("[3DMASC] Classifying %1 points with %2 feature(s)").arg(sampleCount).arg(attributesPerSample)); //create the field wrappers std::vector< IScalarFieldWrapper::Shared > wrappers; { wrappers.reserve(attributesPerSample); for (int fIndex = 0; fIndex < attributesPerSample; ++fIndex) { const Feature::Source& fs = featureSources[fIndex]; IScalarFieldWrapper::Shared source = GetSource(fs, cloud); if (!source || !source->isValid()) { assert(false); errorMessage = QObject::tr("Internal error: invalid source '%1'").arg(fs.name); return false; } wrappers.push_back(source); } } QScopedPointer pDlg; if (parentWidget) { pDlg.reset(new ccProgressDialog(parentWidget)); pDlg->setLabelText(QString("Classify (%1 points)").arg(sampleCount)); pDlg->show(); QCoreApplication::processEvents(); } CCCoreLib::NormalizedProgress nProgress(pDlg.data(), cloud->size()); bool success = true; int numberOfTrees = static_cast(m_rtrees->getRoots().size()); bool cancelled = false; #ifndef _DEBUG #if defined(_OPENMP) #pragma omp parallel for num_threads(omp_get_max_threads() - 2) #endif #endif for (int i = 0; i < static_cast(cloud->size()); ++i) { { //allocate the data matrix cv::Mat test_data; try { test_data.create(1, attributesPerSample, CV_32FC1); } catch (const cv::Exception& cvex) { errorMessage = cvex.msg.c_str(); success = false; cancelled = true; } if (!cancelled) { for (int fIndex = 0; fIndex < attributesPerSample; ++fIndex) { double value = wrappers[fIndex]->pointValue(i); test_data.at(0, fIndex) = static_cast(value); } float predictedClass = m_rtrees->predict(test_data.row(0), cv::noArray(), cv::ml::DTrees::PREDICT_MAX_VOTE); classificationSF->setValue(i, static_cast(predictedClass)); // compute the confidence cv::Mat result; m_rtrees->getVotes(test_data, result, cv::ml::DTrees::PREDICT_MAX_VOTE); int classIndex = -1; for (int col = 0; col < result.cols; col++) // look for the index of the predicted class if (predictedClass == result.at(0, col)) { classIndex = col; break; } if (classIndex != -1) { float nbVotes = result.at(1, classIndex); // get the number of votes cvConfidenceSF->setValue(i, static_cast(nbVotes / numberOfTrees)); // compute the confidence } else cvConfidenceSF->setValue(i, CCCoreLib::NAN_VALUE); if (pDlg && !nProgress.oneStep()) { //process cancelled by the user success = false; cancelled = true; } } } } classificationSF->computeMinAndMax(); cvConfidenceSF->computeMinAndMax(); //show the classification field by default { int classifSFIdx = cloud->getScalarFieldIndexByName(classificationSF->getName()); cloud->setCurrentDisplayedScalarField(classifSFIdx); cloud->showSF(true); } if (parentWidget && cloud->getDisplay()) { cloud->getDisplay()->redraw(); QCoreApplication::processEvents(); } if (classifSFBackup != nullptr) { if (app) { ConfusionMatrix *confusionMatrix = new ConfusionMatrix(*classifSFBackup, *classificationSF); } } return success; } bool Classifier::evaluate(const Feature::Source::Set& featureSources, ccPointCloud* testCloud, AccuracyMetrics& metrics, QString& errorMessage, Train3DMASCDialog& train3DMASCDialog, CCCoreLib::ReferenceCloud* testSubset/*=nullptr=*/, QString outputSFName/*=QString()*/, QWidget* parentWidget/*=nullptr*/, ccMainAppInterface *app/*=nullptr*/) { if (!testCloud) { //invalid input assert(false); errorMessage = QObject::tr("Invalid input cloud"); return false; } metrics.sampleCount = metrics.goodGuess = 0; metrics.ratio = 0.0f; if (!m_rtrees || !m_rtrees->isTrained()) { errorMessage = QObject::tr("Classifier hasn't been trained yet"); return false; } if (featureSources.empty()) { errorMessage = QObject::tr("Training method called without any feature (source)?!"); return false; } if (testSubset && testSubset->getAssociatedCloud() != testCloud) { errorMessage = QObject::tr("Invalid test subset (associated point cloud is different)"); return false; } //look for the classification field CCCoreLib::ScalarField* classifSF = Tools::GetClassificationSF(testCloud); if (!classifSF || classifSF->size() < testCloud->size()) { assert(false); errorMessage = QObject::tr("Missing/invalid 'Classification' field on input cloud"); return false; } ccScalarField* outSF = nullptr; ccScalarField* cvConfidenceSF = nullptr; if (!outputSFName.isEmpty()) { int outIdx = testCloud->getScalarFieldIndexByName(qPrintable(outputSFName)); if (outIdx >= 0) testCloud->deleteScalarField(outIdx); else ccLog::Print("add " + outputSFName + " to the TEST cloud"); outIdx = testCloud->addScalarField(qPrintable(outputSFName)); outSF = static_cast(testCloud->getScalarField(outIdx)); } if (outSF) // add a Classification_confidence value to the test cloud if needed { int cvConfidenceIdx = testCloud->getScalarFieldIndexByName("Classification_confidence"); if (cvConfidenceIdx >= 0) // if the scalar field exists, delete it testCloud->deleteScalarField(cvConfidenceIdx); else ccLog::Print("add Classification_confidence to the TEST cloud"); cvConfidenceIdx = testCloud->addScalarField("Classification_confidence"); cvConfidenceSF = static_cast(testCloud->getScalarField(cvConfidenceIdx)); } unsigned testSampleCount = (testSubset ? testSubset->size() : testCloud->size()); int attributesPerSample = static_cast(featureSources.size()); ccLog::Print(QObject::tr("[3DMASC] Testing data: %1 samples with %2 feature(s)").arg(testSampleCount).arg(attributesPerSample)); //allocate the data matrix cv::Mat test_data; try { test_data.create(static_cast(testSampleCount), attributesPerSample, CV_32FC1); } catch (const cv::Exception& cvex) { errorMessage = cvex.msg.c_str(); return false; } QScopedPointer pDlg; if (parentWidget) { pDlg.reset(new ccProgressDialog(parentWidget)); pDlg->setLabelText(QString("Evaluating the classifier on %1 points").arg(testSampleCount)); pDlg->show(); QCoreApplication::processEvents(); } CCCoreLib::NormalizedProgress nProgress(pDlg.data(), testSampleCount); //fill the data matrix for (int fIndex = 0; fIndex < attributesPerSample; ++fIndex) { const Feature::Source& fs = featureSources[fIndex]; IScalarFieldWrapper::Shared source = GetSource(fs, testCloud); if (!source || !source->isValid()) { assert(false); errorMessage = QObject::tr("Internal error: invalid source '%1'").arg(fs.name); return false; } for (unsigned i = 0; i < testSampleCount; ++i) { unsigned pointIndex = (testSubset ? testSubset->getPointGlobalIndex(i) : i); double value = source->pointValue(pointIndex); test_data.at(i, fIndex) = static_cast(value); } } int numberOfTrees = static_cast(m_rtrees->getRoots().size()); //estimate the efficiency of the classifier std::vector actualClass(testSampleCount); std::vector predictectedClass(testSampleCount); { metrics.sampleCount = testSampleCount; metrics.goodGuess = 0; for (unsigned i = 0; i < testSampleCount; ++i) { unsigned pointIndex = (testSubset ? testSubset->getPointGlobalIndex(i) : i); ScalarType pointClass = classifSF->getValue(pointIndex); int iClass = static_cast(pointClass); //if (iClass < 0 || iClass > 255) //{ // errorMessage = QObject::tr("Classification values out of range (0-255)"); // return false; //} float fPredictedClass = m_rtrees->predict(test_data.row(i), cv::noArray(), cv::ml::DTrees::PREDICT_MAX_VOTE); int iPredictedClass = static_cast(fPredictedClass); actualClass.at(i) = iClass; predictectedClass.at(i) = iPredictedClass; if (iPredictedClass == iClass) { ++metrics.goodGuess; } if (outSF) { outSF->setValue(pointIndex, static_cast(iPredictedClass)); if (cvConfidenceSF) { // compute the confidence cv::Mat result; m_rtrees->getVotes(test_data.row(i), result, cv::ml::DTrees::PREDICT_MAX_VOTE); int classIndex = -1; for (int col = 0; col < result.cols; col++) // look for the index of the predicted class if (iPredictedClass == result.at(0, col)) { classIndex = col; break; } if (classIndex != -1) { float nbVotes = result.at(1, classIndex); // get the number of votes cvConfidenceSF->setValue(i, static_cast(nbVotes / numberOfTrees)); // compute the confidence } else cvConfidenceSF->setValue(i, CCCoreLib::NAN_VALUE); } } if (pDlg && !nProgress.oneStep()) { //process cancelled by the user return false; } } if (outSF) outSF->computeMinAndMax(); if (cvConfidenceSF) cvConfidenceSF->computeMinAndMax(); metrics.ratio = static_cast(metrics.goodGuess) / metrics.sampleCount; } ConfusionMatrix* confusionMatrix = new ConfusionMatrix(actualClass, predictectedClass); train3DMASCDialog.addConfusionMatrixAndSaveTraces(confusionMatrix); if (app) { confusionMatrix->show(); } //show the Classification_prediction field by default if (outSF) { int classifSFIdx = testCloud->getScalarFieldIndexByName(outSF->getName()); testCloud->setCurrentDisplayedScalarField(classifSFIdx); testCloud->showSF(true); } if (parentWidget && testCloud->getDisplay()) { testCloud->getDisplay()->redraw(); QCoreApplication::processEvents(); } return true; } bool Classifier::train( const ccPointCloud* cloud, const RandomTreesParams& params, const Feature::Source::Set& featureSources, QString& errorMessage, CCCoreLib::ReferenceCloud* trainSubset/*=nullptr*/, ccMainAppInterface* app/*=nullptr*/, QWidget* parentWidget/*=nullptr*/) { if (featureSources.empty()) { errorMessage = QObject::tr("Training method called without any feature (source)?!"); return false; } if (!cloud) { errorMessage = QObject::tr("Invalid input cloud"); return false; } if (trainSubset && trainSubset->getAssociatedCloud() != cloud) { errorMessage = QObject::tr("Invalid train subset (associated point cloud is different)"); return false; } //look for the classification field CCCoreLib::ScalarField* classifSF = Tools::GetClassificationSF(cloud); if (!classifSF || classifSF->size() < cloud->size()) { assert(false); errorMessage = QObject::tr("Missing/invalid 'Classification' field on input cloud"); return false; } int sampleCount = static_cast(trainSubset ? trainSubset->size() : cloud->size()); int attributesPerSample = static_cast(featureSources.size()); if (app) { app->dispToConsole(QString("[3DMASC] Training data: %1 samples with %2 feature(s)").arg(sampleCount).arg(attributesPerSample)); } cv::Mat training_data, train_labels; try { training_data.create(sampleCount, attributesPerSample, CV_32FC1); train_labels.create(sampleCount, 1, CV_32FC1); } catch (const cv::Exception& cvex) { errorMessage = cvex.msg.c_str(); return false; } //fill the classification labels vector { for (int i = 0; i < sampleCount; ++i) { int pointIndex = (trainSubset ? static_cast(trainSubset->getPointGlobalIndex(i)) : i); ScalarType pointClass = classifSF->getValue(pointIndex); int iClass = static_cast(pointClass); //if (iClass < 0 || iClass > 255) //{ // errorMessage = QObject::tr("Classification values out of range (0-255)"); // return false; //} train_labels.at(i) = static_cast(iClass); } } //fill the training data matrix for (int fIndex = 0; fIndex < attributesPerSample; ++fIndex) { const Feature::Source& fs = featureSources[fIndex]; IScalarFieldWrapper::Shared source = GetSource(fs, cloud); if (!source || !source->isValid()) { assert(false); errorMessage = QObject::tr("Internal error: invalid source '%1'").arg(fs.name); return false; } for (int i = 0; i < sampleCount; ++i) { int pointIndex = (trainSubset ? static_cast(trainSubset->getPointGlobalIndex(i)) : i); double value = source->pointValue(pointIndex); training_data.at(i, fIndex) = static_cast(value); } } QScopedPointer pDlg; if (parentWidget) { pDlg.reset(new QProgressDialog(parentWidget)); pDlg->setRange(0, 0); //infinite loop pDlg->setLabelText("Training classifier"); pDlg->show(); QCoreApplication::processEvents(); } m_rtrees = cv::ml::RTrees::create(); m_rtrees->setMaxDepth(params.maxDepth); m_rtrees->setMinSampleCount(params.minSampleCount); m_rtrees->setRegressionAccuracy(0); // If true then surrogate splits will be built. These splits allow to work with missing data and compute variable importance correctly. Default value is false. m_rtrees->setUseSurrogates(false); m_rtrees->setPriors(cv::Mat()); //m_rtrees->setMaxCategories(params.maxCategories); //not important? m_rtrees->setCalculateVarImportance(true); m_rtrees->setActiveVarCount(params.activeVarCount); cv::TermCriteria terminationCriteria(cv::TermCriteria::MAX_ITER, params.maxTreeCount, std::numeric_limits::epsilon()); m_rtrees->setTermCriteria(terminationCriteria); QFuture future = QtConcurrent::run([&]() { // Code in this block will run in another thread try { cv::Mat sampleIndexes = cv::Mat::zeros(1, training_data.rows, CV_8U); // cv::Mat trainSamples = sampleIndexes.colRange(0, sampleCount); // trainSamples.setTo(cv::Scalar::all(1)); cv::Mat varTypes(training_data.cols + 1, 1, CV_8U); varTypes.setTo(cv::Scalar::all(cv::ml::VAR_ORDERED)); varTypes.at(training_data.cols) = cv::ml::VAR_CATEGORICAL; cv::Ptr trainData = cv::ml::TrainData::create(training_data, cv::ml::ROW_SAMPLE, train_labels, /* samples layout responses */ cv::noArray(), sampleIndexes, /* varIdx sampleIdx */ cv::noArray(), varTypes); // sampleWeights varType bool success = m_rtrees->train(trainData); if (!success || !m_rtrees->isClassifier()) { errorMessage = "Training failed"; return false; } } catch (const cv::Exception& cvex) { m_rtrees.release(); errorMessage = cvex.msg.c_str(); return false; } catch (const std::exception& stdex) { errorMessage = stdex.what(); return false; } catch (...) { errorMessage = QObject::tr("Unknown error"); return false; } return true; }); while (!future.isFinished()) { #if defined(CC_WINDOWS) ::Sleep(500); #else usleep(500 * 1000); #endif if (pDlg) { if (pDlg->wasCanceled()) { // future.cancel(); QMessageBox msgBox; msgBox.setText("The training is still in progress, not possible to cancel."); msgBox.exec(); // break; pDlg->reset(); pDlg->show(); } pDlg->setValue(pDlg->value() + 1); } QCoreApplication::processEvents(); } if (pDlg) { pDlg->close(); QCoreApplication::processEvents(); } if (future.isCanceled() || !future.result() || !m_rtrees->isTrained()) { errorMessage = QObject::tr("Training failed for an unknown reason..."); m_rtrees.release(); return false; } return true; } bool Classifier::toFile(QString filename, QWidget* parentWidget/*=nullptr*/) const { if (!m_rtrees) { ccLog::Warning(QObject::tr("Classifier hasn't been trained, can't save it")); return false; } //save the classifier QProgressDialog pDlg(parentWidget); pDlg.setRange(0, 0); //infinite loop pDlg.setLabelText(QObject::tr("Saving classifier")); pDlg.show(); QCoreApplication::processEvents(); cv::String cvFilename = filename.toStdString(); m_rtrees->save(cvFilename); pDlg.close(); QCoreApplication::processEvents(); ccLog::Print("Classifier file saved to: " + QString::fromStdString(cvFilename)); return true; } bool Classifier::fromFile(QString filename, QWidget* parentWidget/*=nullptr*/) { //load the classifier QScopedPointer pDlg; if (parentWidget) { pDlg.reset(new QProgressDialog(parentWidget)); pDlg->setRange(0, 0); //infinite loop pDlg->setLabelText(QObject::tr("Loading classifier")); pDlg->show(); QCoreApplication::processEvents(); } try { m_rtrees = cv::ml::RTrees::load(filename.toStdString()); } catch (const cv::Exception& cvex) { ccLog::Warning(cvex.msg.c_str()); ccLog::Error("Failed to load file: " + filename); return false; } if (pDlg) { pDlg->close(); QCoreApplication::processEvents(); } if (m_rtrees->empty() || !m_rtrees->isClassifier()) { ccLog::Error(QObject::tr("Loaded classifier is invalid")); return false; } else if (!m_rtrees->isTrained()) { ccLog::Warning(QObject::tr("Loaded classifier doesn't seem to be trained")); } return true; }