Fraxinus  2023.01.05-dev+develop.0da12
An IGT application
cxtestSyntheticVolumeComparer.cpp
Go to the documentation of this file.
1 /*=========================================================================
2 This file is part of CustusX, an Image Guided Therapy Application.
3 
4 Copyright (c) SINTEF Department of Medical Technology.
5 All rights reserved.
6 
7 CustusX is released under a BSD 3-Clause license.
8 
9 See Lisence.txt (https://github.com/SINTEFMedtek/CustusX/blob/master/License.txt) for details.
10 =========================================================================*/
11 
13 #include "cxVolumeHelpers.h"
14 #include "catch.hpp"
15 #include "cxImage.h"
16 #include "cxtestUtilities.h"
17 #include <vtkImageData.h>
19 #include <QDir>
20 #include <QString>
21 #include "cxTypeConversions.h"
22 #include "cxFileManagerService.h"
23 
24 namespace cxtest
25 {
26 
28 {
29 }
30 
32 {
33  mPhantom = phantom;
34  mNominalImage.reset();
35 }
36 
38 {
39  mTestImage = image;
40  mNominalImage.reset();
41 }
42 
43 //vtkImageDataPtr Utilities::create3DVtkImageData(Eigen::Array3i dim, const unsigned int voxelValue)
44 //{
45 // return cx::generateVtkImageData(dim, cx::Vector3D(1,1,1), voxelValue);
46 //}
47 
48 //cx::ImagePtr Utilities::create3DImage(Eigen::Array3i dim, const unsigned int voxelValue)
49 //{
50 // vtkImageDataPtr vtkImageData = create3DVtkImageData(dim, voxelValue);
51 // QString unique_string = qstring_cast(reinterpret_cast<long>(vtkImageData.GetPointer()));
52 // QString imagesUid = QString("TESTUID_%2_%1").arg(unique_string);
53 // cx::ImagePtr image(new cx::Image(imagesUid, vtkImageData));
54 
55 // return image;
56 //}
57 
58 
59 cx::ImagePtr SyntheticVolumeComparer::getNominalOutputImage() const
60 {
61  REQUIRE(mTestImage);
62 
63  if (!mNominalImage)
64  {
65  Eigen::Array3i dim(mTestImage->getBaseVtkImageData()->GetDimensions());
66  cx::Vector3D spacing(mTestImage->getBaseVtkImageData()->GetSpacing());
67  mNominalImage = cxtest::Utilities::create3DImage(dim, spacing, 0);
68  mNominalImage->get_rMd_History()->setRegistration(mTestImage->get_rMd());
69 
70  mPhantom->fillVolume(mNominalImage);
71  }
72 
73  return mNominalImage;
74 }
75 
77 {
78  float sse = this->getRMS();
79  if (this->getVerbose())
80  std::cout << "RMS value: " << sse << std::endl;
81  CHECK(sse < threshold);
82 }
83 
84 double SyntheticVolumeComparer::getRMS() const
85 {
86  double sse = cx::calculateRMSError(mTestImage->getBaseVtkImageData(), this->getNominalOutputImage()->getBaseVtkImageData());
87 // float sse = mPhantom->computeRMSError(mOutputData);
88 // std::cout << "RMS value: " << sse << std::endl;
89  return sse;
90 }
91 
93 {
94  cx::Vector3D c_n = calculateCentroid(this->getNominalOutputImage());
95  cx::Vector3D c_r = calculateCentroid(mTestImage);
96 
97  double difference = (c_n-c_r).norm();
98 
99  if (this->getVerbose())
100  std::cout << "c_n=[" << c_n << "] c_r=[" << c_r << "] diff=[" << difference << "]" << std::endl;
101 
102  CHECK(difference < val);
103 }
104 
106 {
107  double v_n = calculateMass(this->getNominalOutputImage());
108  double v_r = calculateMass(mTestImage);
109  double normalized_difference = fabs(v_n-v_r)/(v_n+v_r);
110 
111  if (this->getVerbose())
112  std::cout << "v_n=[" << v_n << "] v_r=[" << v_r << "] diff=[" << normalized_difference << "]" << std::endl;
113 
114  CHECK(normalized_difference<val);
115 }
116 
117 void SyntheticVolumeComparer::checkValueWithin(cx::Vector3D p_r, int lowerLimit, int upperLimit)
118 {
119  int val = this->getValue(mTestImage, p_r);
120  if (this->getVerbose())
121  std::cout << QString("p_r[%1]=%2, limit=[%3, %4]").arg(qstring_cast(p_r)).arg(val).arg(lowerLimit).arg(upperLimit) << std::endl;
122  CHECK(val>=lowerLimit);
123  CHECK(val<=upperLimit);
124 }
125 
126 double SyntheticVolumeComparer::getValue(cx::ImagePtr image, cx::Vector3D p_r)
127 {
128  vtkImageDataPtr raw = image->getBaseVtkImageData();
129  cx::Vector3D p_d = image->get_rMd().inv().coord(p_r);
130  Eigen::Array3d spacing(raw->GetSpacing());
131  Eigen::Array3i voxel = cx::round(p_r.array() / spacing).cast<int>();
132  double val = raw->GetScalarComponentAsDouble(voxel.data()[0],voxel.data()[1],voxel.data()[2], 0);
133  return val;
134 }
135 
137 {
138  port->save(this->getNominalOutputImage(), this->addFullPath(filename));
139 }
140 
142 {
143  port->save(mTestImage, this->addFullPath(filename));
144 }
145 
146 QString SyntheticVolumeComparer::addFullPath(QString filename)
147 {
148  QDir dir(QDir::homePath() + "/test/");
149  dir.mkdir(".");
150  QString fullFilename = dir.absoluteFilePath(filename);
151  return fullFilename;
152 }
153 
154 } // namespace cxtest
155 
156 
QString qstring_cast(const T &val)
boost::shared_ptr< class FileManagerService > FileManagerServicePtr
double calculateMass(cx::ImagePtr image)
void saveOutputToFile(QString filename, cx::FileManagerServicePtr port)
boost::shared_ptr< class Image > ImagePtr
Definition: cxDicomWidget.h:27
static cx::ImagePtr create3DImage(Eigen::Array3i dim=Eigen::Array3i(3, 3, 3), const unsigned int voxelValue=100)
void checkValueWithin(cx::Vector3D p_r, int lowerLimit, int upperLimit)
boost::shared_ptr< cxSyntheticVolume > cxSyntheticVolumePtr
cx::Vector3D calculateCentroid(cx::ImagePtr image)
void saveNominalOutputToFile(QString filename, cx::FileManagerServicePtr port)
void setPhantom(cx::cxSyntheticVolumePtr phantom)
double calculateRMSError(vtkImageDataPtr a, vtkImageDataPtr b)
Eigen::Vector3d Vector3D
Vector3D is a representation of a point or vector in 3D.
Definition: cxVector3D.h:42
vtkSmartPointer< class vtkImageData > vtkImageDataPtr
Vector3D round(const Vector3D &a)
Definition: cxVector3D.cpp:75