diff --git a/Plugins/org.mitk.gui.qt.igt.app.hummelprotocolmeasurements/src/internal/mitkHummelProtocolEvaluation.cpp b/Plugins/org.mitk.gui.qt.igt.app.hummelprotocolmeasurements/src/internal/mitkHummelProtocolEvaluation.cpp index d37fb8f074..53d7be616d 100644 --- a/Plugins/org.mitk.gui.qt.igt.app.hummelprotocolmeasurements/src/internal/mitkHummelProtocolEvaluation.cpp +++ b/Plugins/org.mitk.gui.qt.igt.app.hummelprotocolmeasurements/src/internal/mitkHummelProtocolEvaluation.cpp @@ -1,384 +1,388 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #define _USE_MATH_DEFINES #include "mitkHummelProtocolEvaluation.h" #include #include #include #include #include #include #include #include bool mitk::HummelProtocolEvaluation::Evaluate15cmDistances(mitk::PointSet::Pointer p, HummelProtocolMeasurementVolume m, std::vector &Results) { if (m != mitk::HummelProtocolEvaluation::standard) { MITK_WARN << "15 cm distances are only evaluated for standard volumes, aborting!"; return false; } MITK_INFO << "########### 15 cm distance errors #############"; //convert measurements to matrix - itk::Matrix, 9, 10> matrix = ParseMatrixStandardVolume(p); + std::array ,9> matrix = ParseMatrixStandardVolume(p); //these are the variables for the results: std::vector distances; std::vector descriptions; //evaluation of rows int distanceCounter = 0; for (int row = 0; row < 9; row++) //rows for (int distance = 0; distance < 7; distance++) { distanceCounter++; mitk::Point3D point1 = p->GetPoint(row * 10 + distance); mitk::Point3D point2 = p->GetPoint(row * 10 + distance + 3); distances.push_back(point1.EuclideanDistanceTo(point2)); std::stringstream description; description << "Distance(" << distanceCounter << ") " << (row + 1) << "/" << (distance + 1) << " to " << (row + 1) << "/" << (distance + 4); descriptions.push_back(description.str()); } //evaluation of columns for (int column = 0; column < 10; column++) for (int row = 0; row < 6; row++) { distanceCounter++; mitk::Point3D point1 = matrix[row][column]; mitk::Point3D point2 = matrix[row + 3][column]; distances.push_back(point1.EuclideanDistanceTo(point2)); std::stringstream description; description << "Distance(" << distanceCounter << ") " << (row + 1) << "/" << (column + 1) << " to " << (row + 4) << "/" << (column + 1); descriptions.push_back(description.str()); } //compute all errors for (int i = 0; i < distances.size(); i++) { HummelProtocolDistanceError currentError; - currentError.distanceError = abs(distances.at(i) - (double)150.0); + currentError.distanceError = fabs(distances.at(i) - (double)150.0); currentError.description = descriptions.at(i); Results.push_back(currentError); MITK_INFO << "Error " << currentError.description << " : " << currentError.distanceError; } //compute statistics std::vector statistics = mitk::HummelProtocolEvaluation::ComputeStatistics(Results); for (auto currentError : statistics) { Results.push_back(currentError); MITK_INFO << currentError.description << " : " << currentError.distanceError; } return true; } bool mitk::HummelProtocolEvaluation::Evaluate30cmDistances(mitk::PointSet::Pointer p, HummelProtocolMeasurementVolume m, std::vector &Results) { if (m != mitk::HummelProtocolEvaluation::standard) { MITK_WARN << "30 cm distances are only evaluated for standard volumes, aborting!"; return false; } MITK_INFO << "########### 30 cm distance errors #############"; //convert measurements to matrix - itk::Matrix, 9, 10> matrix = ParseMatrixStandardVolume(p); + std::array ,9> matrix = ParseMatrixStandardVolume(p); //these are the variables for the results: std::vector distances; std::vector descriptions; //evaluation of rows int distanceCounter = 0; for (int row = 0; row < 9; row++) //rows for (int distance = 0; distance < 4; distance++) { distanceCounter++; mitk::Point3D point1 = p->GetPoint(row * 10 + distance); mitk::Point3D point2 = p->GetPoint(row * 10 + distance + 6); distances.push_back(point1.EuclideanDistanceTo(point2)); std::stringstream description; description << "Distance(" << distanceCounter << ") " << (row + 1) << "/" << (distance + 1) << " to " << (row + 1) << "/" << (distance + 7); descriptions.push_back(description.str()); } //evaluation of columns for (int column = 0; column < 10; column++) for (int row = 0; row < 3; row++) { distanceCounter++; mitk::Point3D point1 = matrix[row][column]; mitk::Point3D point2 = matrix[row + 6][column]; distances.push_back(point1.EuclideanDistanceTo(point2)); std::stringstream description; description << "Distance(" << distanceCounter << ") " << (row + 1) << "/" << (column + 1) << " to " << (row + 7) << "/" << (column + 1); descriptions.push_back(description.str()); } //compute all errors for (int i = 0; i < distances.size(); i++) { HummelProtocolDistanceError currentError; - currentError.distanceError = abs(distances.at(i) - (double)300.0); + currentError.distanceError = fabs(distances.at(i) - (double)300.0); currentError.description = descriptions.at(i); Results.push_back(currentError); MITK_INFO << "Error " << currentError.description << " : " << currentError.distanceError; } //compute statistics std::vector statistics = mitk::HummelProtocolEvaluation::ComputeStatistics(Results); for (auto currentError : statistics) { Results.push_back(currentError); MITK_INFO << currentError.description << " : " << currentError.distanceError; } return true; } bool mitk::HummelProtocolEvaluation::EvaluateAccumulatedDistances(mitk::PointSet::Pointer p, HummelProtocolMeasurementVolume m, std::vector &Results) { if (m != mitk::HummelProtocolEvaluation::standard) { MITK_WARN << "Accumulated distances are only evaluated for standard volumes, aborting!"; return false; } //convert measurements to matrix - itk::Matrix, 9, 10> matrix = ParseMatrixStandardVolume(p); + std::array ,9> matrix = ParseMatrixStandardVolume(p); MITK_INFO << "########### accumulated distance errors #############"; int distanceCounter = 0; //evaluation of rows for (int row = 0; row < 9; row++) //rows for (int distance = 0; distance < 9; distance++) { distanceCounter++; mitk::Point3D point1 = p->GetPoint(row * 10); mitk::Point3D point2 = p->GetPoint(row * 10 + distance + 1); std::stringstream description; description << "Distance(" << distanceCounter << ") " << (row + 1) << "/1 to " << (row + 1) << "/" << (distance + 2); //compute error HummelProtocolDistanceError currentError; - currentError.distanceError = abs(point1.EuclideanDistanceTo(point2) - (double)(50.0*(distance+1))); + currentError.distanceError = fabs(point1.EuclideanDistanceTo(point2) - (double)(50.0*(distance+1))); currentError.description = description.str(); Results.push_back(currentError); MITK_INFO << "Error " << currentError.description << " : " << currentError.distanceError; } //compute statistics std::vector statistics = mitk::HummelProtocolEvaluation::ComputeStatistics(Results); for (auto currentError : statistics) { Results.push_back(currentError); MITK_INFO << currentError.description << " : " << currentError.distanceError; } return true; } bool mitk::HummelProtocolEvaluation::Evaluate5cmDistances(mitk::PointSet::Pointer p, HummelProtocolMeasurementVolume m, std::vector &Results) { MITK_INFO << "########### 5 cm distance errors #############"; std::vector distances; std::vector descriptions; switch (m) { case small: if (p->GetSize() != 12) { MITK_WARN << "Wrong number of points: " << p->GetSize() << " (expected 12), aborting"; return false; } MITK_INFO << "Computing Hummel protocol distance errors for small measurement volumes (12 points)..."; //row 1 distances.push_back(p->GetPoint(0).EuclideanDistanceTo(p->GetPoint(1))); //0 descriptions.push_back("Distance 4/4 to 4/5"); distances.push_back(p->GetPoint(1).EuclideanDistanceTo(p->GetPoint(2))); //1 descriptions.push_back("Distance 4/5 to 4/6"); distances.push_back(p->GetPoint(2).EuclideanDistanceTo(p->GetPoint(3))); //2 descriptions.push_back("Distance 4/6 to 4/7"); //row 2 distances.push_back(p->GetPoint(4).EuclideanDistanceTo(p->GetPoint(5))); //3 descriptions.push_back("Distance 5/4 to 5/5"); distances.push_back(p->GetPoint(5).EuclideanDistanceTo(p->GetPoint(6))); //4 descriptions.push_back("Distance 5/5 to 5/6"); distances.push_back(p->GetPoint(6).EuclideanDistanceTo(p->GetPoint(7))); //5 descriptions.push_back("Distance 5/6 to 5/7"); //row 3 distances.push_back(p->GetPoint(8).EuclideanDistanceTo(p->GetPoint(9))); //6 descriptions.push_back("Distance 6/4 to 6/5"); distances.push_back(p->GetPoint(9).EuclideanDistanceTo(p->GetPoint(10))); //7 descriptions.push_back("Distance 6/5 to 6/6"); distances.push_back(p->GetPoint(10).EuclideanDistanceTo(p->GetPoint(11))); //8 descriptions.push_back("Distance 6/6 to 6/7"); //column 1 distances.push_back(p->GetPoint(0).EuclideanDistanceTo(p->GetPoint(4))); //9 descriptions.push_back("Distance 4/4 to 5/4"); distances.push_back(p->GetPoint(4).EuclideanDistanceTo(p->GetPoint(8))); //10 descriptions.push_back("Distance 5/4 to 6/4"); //column 2 distances.push_back(p->GetPoint(1).EuclideanDistanceTo(p->GetPoint(5))); //11 descriptions.push_back("Distance 4/5 to 5/5"); distances.push_back(p->GetPoint(5).EuclideanDistanceTo(p->GetPoint(9))); //12 descriptions.push_back("Distance 5/5 to 6/5"); //column 3 distances.push_back(p->GetPoint(2).EuclideanDistanceTo(p->GetPoint(6))); //13 descriptions.push_back("Distance 4/6 to 5/6"); distances.push_back(p->GetPoint(6).EuclideanDistanceTo(p->GetPoint(10))); //14 descriptions.push_back("Distance 5/6 to 6/6"); //column 4 distances.push_back(p->GetPoint(3).EuclideanDistanceTo(p->GetPoint(7))); //15 descriptions.push_back("Distance 4/7 to 5/7"); distances.push_back(p->GetPoint(7).EuclideanDistanceTo(p->GetPoint(11))); //16 descriptions.push_back("Distance 5/7 to 6/7"); break; case standard: if (p->GetSize() != 90) { MITK_WARN << "Wrong number of points (expected 90), aborting"; return false; } MITK_INFO << "Computing Hummel protocol distance errors for standard measurement volumes (90 points)..."; int distanceCounter = 0; //convert measurements to matrix -itk::Matrix, 9, 10> matrix = ParseMatrixStandardVolume(p); +std::array ,9> matrix = ParseMatrixStandardVolume(p); //evaluation of rows for (int row = 0; row < 9; row++) //rows for (int distance = 0; distance < 9; distance++) { distanceCounter++; mitk::Point3D point1 = p->GetPoint(row*10 + distance); mitk::Point3D point2 = p->GetPoint(row*10 + distance+1); distances.push_back(point1.EuclideanDistanceTo(point2)); std::stringstream description; description << "Distance(" << distanceCounter << ") " << (row + 1) << "/" << (distance + 1) << " to " << (row + 1) << "/" << (distance + 2); descriptions.push_back(description.str()); } //evaluation of columns for (int column = 0; column < 10; column++) for (int row = 0; row < 8; row++) { distanceCounter++; mitk::Point3D point1 = matrix[row][column]; mitk::Point3D point2 = matrix[row+1][column]; distances.push_back(point1.EuclideanDistanceTo(point2)); std::stringstream description; description << "Distance(" << distanceCounter << ") " << (row+1 )<< "/" << (column+1) << " to " << (row + 2) << "/" << (column + 1); descriptions.push_back(description.str()); } break; } //compute all errors for (int i = 0; i < distances.size(); i++) { HummelProtocolDistanceError currentError; -currentError.distanceError = abs(distances.at(i) - (double)50.0); +currentError.distanceError = fabs(distances.at(i) - (double)50.0); currentError.description = descriptions.at(i); Results.push_back(currentError); MITK_INFO << "Error " << currentError.description << " : " << currentError.distanceError; } //compute statistics std::vector statistics = mitk::HummelProtocolEvaluation::ComputeStatistics(Results); for (auto currentError : statistics) { Results.push_back(currentError); MITK_INFO << currentError.description << " : " << currentError.distanceError; } return true; } -itk::Matrix, 9, 10> mitk::HummelProtocolEvaluation::ParseMatrixStandardVolume(mitk::PointSet::Pointer p) +std::array, 9> mitk::HummelProtocolEvaluation::ParseMatrixStandardVolume(mitk::PointSet::Pointer p) { - itk::Matrix, 9, 10> returnValue = itk::Matrix, 9, 10>(); + + std::array ,9> returnValue; + if (p->GetSize() != 90) { MITK_WARN << "PointSet does not have the right size. Expected 90 got " << p->GetSize() << " ... aborting!"; return returnValue; } for (int row = 0; row < 9; row++) for (int column = 0; column < 10; column++) returnValue[row][column] = p->GetPoint(row * 10 + column); + } std::vector mitk::HummelProtocolEvaluation::ComputeStatistics(std::vector values) { std::vector returnValue; //convert input values to boost / using boost accumulators for statistics boost::accumulators::accumulator_set > > acc; for (mitk::HummelProtocolEvaluation::HummelProtocolDistanceError each : values) { acc(each.distanceError); } returnValue.push_back({ values.size(), "N" }); returnValue.push_back({ boost::accumulators::mean(acc), "Mean" }); //double quantile25th = boost::accumulators::quantile(acc, boost::accumulators::quantile_probability = 0.25); //returnValue.push_back({ boost::accumulators::median(acc), "Median" }); //returnValue.push_back({ boost::accumulators::variance(acc), "Variance" }); returnValue.push_back({ boost::accumulators::min(acc), "Min" }); returnValue.push_back({ boost::accumulators::max(acc), "Max" }); //don't get the boost stuff working correctly, so computing the quantiles, median and standard deviation by myself: std::vector quantile; for (mitk::HummelProtocolEvaluation::HummelProtocolDistanceError each : values) {quantile.push_back(each.distanceError);} auto const Q1 = quantile.size() / 4; auto const Q2 = quantile.size() / 2; auto const Q3 = Q1 + Q2; std::sort(quantile.begin(),quantile.end()); returnValue.push_back({ quantile[Q1], "Quartile 1" }); returnValue.push_back({ quantile[Q2], "Median" }); returnValue.push_back({ quantile[Q3], "Quartile 3" }); double mean = boost::accumulators::mean(acc); double errorSum = 0; for (mitk::HummelProtocolEvaluation::HummelProtocolDistanceError each : values) { double error = pow((each.distanceError - mean),2); errorSum += error; } double variance = errorSum / values.size(); double sampleVariance = errorSum / (values.size()-1); double standardDev = sqrt(variance); double sampleStandardDev = sqrt(sampleVariance); returnValue.push_back({ variance, "Variance" }); returnValue.push_back({ sampleVariance, "Sample Variance" }); returnValue.push_back({ standardDev, "Standard Deviation" }); returnValue.push_back({ sampleStandardDev, "Sample Standard Deviation" }); return returnValue; } + diff --git a/Plugins/org.mitk.gui.qt.igt.app.hummelprotocolmeasurements/src/internal/mitkHummelProtocolEvaluation.h b/Plugins/org.mitk.gui.qt.igt.app.hummelprotocolmeasurements/src/internal/mitkHummelProtocolEvaluation.h index 9f20434e89..08d2d8ac9f 100644 --- a/Plugins/org.mitk.gui.qt.igt.app.hummelprotocolmeasurements/src/internal/mitkHummelProtocolEvaluation.h +++ b/Plugins/org.mitk.gui.qt.igt.app.hummelprotocolmeasurements/src/internal/mitkHummelProtocolEvaluation.h @@ -1,86 +1,96 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef MITKHummelProtocolEvaluation_H_HEADER_INCLUDED_ #define MITKHummelProtocolEvaluation_H_HEADER_INCLUDED_ #include -#include +#include +#include +#include +#include namespace mitk { /**Documentation * \brief Static methods for evaluations according to the assessment protocol * for EM trackers published by Hummel et al. 2005 [1]. * * [1] Hummel, J. et al. - Design and application of an assessment protocol for electromagnetic tracking systems. Med Phys 32(7), July 2005 * * \ingroup IGT */ class HummelProtocolEvaluation { public: + /** Distance error with description. */ struct HummelProtocolDistanceError {double distanceError; std::string description;}; /** Tracking volumes for evaluation. * standard: The standard volume of 9 x 10 measurment points as described in [1] * small: A small volume in the center 3 x 4 measurement points, for smaller field generators [2] * [2] Maier-Hein, L. et al. - Standardized assessment of new electromagnetic field generators in an interventional radiology setting. Med Phys 39(6), June 2012 */ enum HummelProtocolMeasurementVolume { small, standard }; /** Evaluates the 5 cm distances as defined by the Hummel protocol [1,2]. * @return Returns true if evaluation was successfull, false if not. * @param[out] Results Please give an empty vector. The results will be added to this vector. */ static bool Evaluate5cmDistances(mitk::PointSet::Pointer p, HummelProtocolMeasurementVolume m, std::vector &Results); /** Evaluates the 15 cm distances as defined by the Hummel protocol [1,2]. * @return Returns true if evaluation was successfull, false if not. * @param[out] Results Please give an empty vector. The results will be added to this vector. */ static bool Evaluate15cmDistances(mitk::PointSet::Pointer p, HummelProtocolMeasurementVolume m, std::vector &Results); /** Evaluates the 30 cm distances as defined by the Hummel protocol [1,2]. * @return Returns true if evaluation was successfull, false if not. * @param[out] Results Please give an empty vector. The results will be added to this vector. */ static bool Evaluate30cmDistances(mitk::PointSet::Pointer p, HummelProtocolMeasurementVolume m, std::vector &Results); /** Evaluates the accumulated distances as defined by the Hummel protocol [1,2]. * @return Returns true if evaluation was successfull, false if not. * @param[out] Results Please give an empty vector. The results will be added to this vector. */ static bool EvaluateAccumulatedDistances(mitk::PointSet::Pointer p, HummelProtocolMeasurementVolume m, std::vector &Results); /** Computes statistics (as mean, standard deviation, quantiles, min, max, etc.) on the given values. * The results are stored inside the return value. */ static std::vector ComputeStatistics(std::vector values); protected: /** Converts a pointset holding all measurement points of the hummel protocol in line-by-line order * to an array representing the hummel board. */ - static itk::Matrix, 9, 10> ParseMatrixStandardVolume(mitk::PointSet::Pointer p); + static std::array ,9> ParseMatrixStandardVolume(mitk::PointSet::Pointer p); + //It would be really wonderfull if we could replace std::array ,9> by mitk::Matrix< mitk::Point3D, 9, 10 > but + //unfortunatly this version does not compile under Linux. To be precise under Linux only matrices like this: mitk::Matriy compile + //even the usage of a double pointer (eg mitk::Matrix) does not compile. We always got an error message saying: + //vnl_c_vector.h:42:49: error: invalid use of incomplete type ‘class vnl_numeric_traits >’ + //Under Windows this error does not appear there everything compiles fine. + }; } // namespace mitk #endif /* MITKHummelProtocolEvaluation_H_HEADER_INCLUDED_ */ diff --git a/Plugins/org.mitk.gui.qt.igt.app.hummelprotocolmeasurements/src/internal/mitkNavigationDataCSVSequentialPlayer.h b/Plugins/org.mitk.gui.qt.igt.app.hummelprotocolmeasurements/src/internal/mitkNavigationDataCSVSequentialPlayer.h index 011e75f06b..11ba570b96 100644 --- a/Plugins/org.mitk.gui.qt.igt.app.hummelprotocolmeasurements/src/internal/mitkNavigationDataCSVSequentialPlayer.h +++ b/Plugins/org.mitk.gui.qt.igt.app.hummelprotocolmeasurements/src/internal/mitkNavigationDataCSVSequentialPlayer.h @@ -1,114 +1,114 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef MITKNavigationDataCSVSequentialPlayer_H_HEADER_INCLUDED_ #define MITKNavigationDataCSVSequentialPlayer_H_HEADER_INCLUDED_ #include #include "tinyxml.h" namespace mitk { /**Documentation * \brief This class is a NavigationDataPlayer which can play CSV formatted * files in sequential order, which means it doesn't care about timestamps and just * outputs the navigationdatas in their sequential order. * * It is thought to interpret custom csv files. To do so please adapt the column * numbers of position and orientation in the internal method GetNavigationDataOutOfOneLine(). * * So far only one (the first) tool is read in as required for the hummel protocol measurements. * * This class can also interpret MITK style csv files (set filetype to NavigationDataCSV), but * you can also use the MITK navigation data player class inside the MITK-IGT module which * is newer and better maintained. * * \ingroup IGT */ class NavigationDataCSVSequentialPlayer : public NavigationDataPlayerBase { public: mitkClassMacro(NavigationDataCSVSequentialPlayer, NavigationDataPlayerBase); itkNewMacro(Self); /** * \brief sets the file name and path (if XMLString is set, this is neglected) */ void SetFileName(const std::string& _FileName); /** * \brief returns the file name and path */ itkGetStringMacro(FileName); enum Filetype { NavigationDataCSV, //for csv files from the MITK navigation data player ManualLoggingCSV //for custum csv files }; /** * \brief Sets the file type. ManualLoggingCSV is default and is thought for your * custom csv files. You can also set it to NavigationDataCSV, then this * player interprets MITK style csv files. */ itkSetMacro(Filetype, Filetype); /** * \return Returns true if the player reached the end of the file. */ bool IsAtEnd(); /** * \brief Used for pipeline update just to tell the pipeline * that we always have to update */ virtual void UpdateOutputInformation(); - int mitk::NavigationDataCSVSequentialPlayer::GetNumberOfSnapshots(); + int GetNumberOfSnapshots(); protected: NavigationDataCSVSequentialPlayer(); virtual ~NavigationDataCSVSequentialPlayer(); /// /// do the work here /// virtual void GenerateData(); std::string m_FileName; int m_CurrentPos; Filetype m_Filetype; //member for the navigation datas which were read (only one output is supported at the moment) std::vector m_NavigationDatas; std::vector GetNavigationDatasFromFile(std::string filename); std::vector GetFileContentLineByLine(std::string filename); mitk::NavigationData::Pointer GetNavigationDataOutOfOneLine(std::string line); void FillOutputEmpty(int number); mitk::NavigationData::Pointer GetEmptyNavigationData(); }; } // namespace mitk #endif /* MITKNavigationDataCSVSequentialPlayer_H_HEADER_INCLUDED_ */