Fraxinus  2023.01.05-dev+develop.0da12
An IGT application
cxSyntheticVolume.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 
12 #include "cxSyntheticVolume.h"
13 #include "vtkImageData.h"
14 #include "cxImage.h"
15 #include <cstdlib>
16 #include <time.h>
17 #include "cxTypeConversions.h"
18 #include <QTime>
19 #include "cxLogger.h"
21 #include "cxVolumeHelpers.h"
22 
23 
24 double noiseValue(double noiseSigma,
25  double noiseMean)
26 {
27  double random_value_1 = (rand()+1.0)/(RAND_MAX+1.0);
28  double random_value_2 = (rand()+1.0)/(RAND_MAX+1.0);
29 
30  double random_normal = sqrt(-2*log(random_value_1)) * cos(2*M_PI*random_value_2);
31 
32  return random_normal*noiseSigma + noiseMean;
33 }
34 
35 namespace cx {
36 
37 
39 cxSyntheticVolume::sampleUsData(const std::vector<Transform3D>& planes_rMt,
40  const ProbeDefinition& probe,
41  const Transform3D& output_dMr,
42  const double noiseSigma,
43  const unsigned char noiseMean) const
44 {
45  cx::ProbeSector sector;
46  sector.setData(probe);
47 
48  std::vector<Transform3D> planes_rMf(planes_rMt.size());
49  for (unsigned i=0; i<planes_rMt.size(); ++i)
50  planes_rMf[i] = planes_rMt[i] * sector.get_tMu() * sector.get_uMv();
51 
52  Eigen::Array2f pixelSpacing = probe.getSpacing().block(0,0,2,1).cast<float>();
53  Eigen::Array2i sliceDimension(probe.getSize().width(), probe.getSize().height());
54 
55  return this->sampleUsData(planes_rMf,
56  pixelSpacing,
57  sliceDimension,
58  output_dMr,
59  noiseSigma, noiseMean);
60 }
61 
64  const ProbeDefinition& probe,
65  const double noiseSigma,
66  const unsigned char noiseMean) const
67 {
68  cx::ProbeSector sector;
69  sector.setData(probe);
70 
71  Transform3D rMf = plane_rMt * sector.get_tMu() * sector.get_uMv();
72  Eigen::Array2f pixelSpacing = probe.getSpacing().block(0,0,2,1).cast<float>();
73  Eigen::Array2i sliceDimension(probe.getSize().width(), probe.getSize().height());
74 
75  return this->sampleUsData(rMf,
76  pixelSpacing,
77  sliceDimension,
78  noiseSigma, noiseMean);
79 }
80 
82 cxSyntheticVolume::sampleUsData(const std::vector<Transform3D>& planes_rMf,
83  const Eigen::Array2f& pixelSpacing,
84  const Eigen::Array2i& sliceDimension,
85  const Transform3D& output_dMr,
86  const double noiseSigma, const unsigned char noiseMean) const
87 {
88  std::vector<TimedPosition> positions;
89  std::vector<vtkImageDataPtr> images;
90  // For each plane
91  for(std::vector<Transform3D>::const_iterator i = planes_rMf.begin();
92  planes_rMf.end() != i;
93  ++i)
94  {
95  const Transform3D rMf = *i;
96 
97  vtkImageDataPtr us_frame;
98  us_frame = this->sampleUsData(rMf, pixelSpacing, sliceDimension, noiseSigma, noiseMean);
99 
100  // Build the TimedPosition for this frame
101  TimedPosition t;
102  t.mTime = i - planes_rMf.begin();
103  t.mPos = output_dMr*rMf;
104 
105  positions.push_back(t);
106  images.push_back(us_frame);
107  }
108 
109  vtkImageDataPtr mask = this->createEmptyMask(sliceDimension);
110 // std::cout << "elapsed: " << time.elapsed() << std::endl;
111 
113  ret.reset(new ProcessedUSInputData(images, positions, mask, "VIRTUAL_DATA", "VIRTUAL_DATA_"));
114  return ret;
115 }
116 
119  const Eigen::Array2f& pixelSpacing,
120  const Eigen::Array2i& sliceDimension,
121  const double noiseSigma, const unsigned char noiseMean) const
122 {
123 
124 
125  const Vector3D p0 = rMf.coord(Vector3D(0,0,0));
126  const Vector3D e_x = rMf.vector(Vector3D(pixelSpacing[0],0,0));
127  const Vector3D e_y = rMf.vector(Vector3D(0,pixelSpacing[1],0));
128 
129  vtkImageDataPtr us_frame = vtkImageDataPtr::New();
130  us_frame->SetExtent(0, sliceDimension[0]-1, 0, sliceDimension[1]-1, 0, 0);
131  us_frame->SetSpacing(pixelSpacing[0], pixelSpacing[1], 0.0);
132  us_frame->AllocateScalars(VTK_UNSIGNED_CHAR, 1);
133 
134  unsigned char* us_data = (unsigned char*)us_frame->GetScalarPointer();
135  // For each pixel on that plane
136  for(unsigned int px = 0; px < sliceDimension[0]; px++)
137  {
138  // optimization: use transformed pixel vectors
139  const Vector3D px0_vol = p0 + e_x*px;
140 
141  for(unsigned int py = 0; py < sliceDimension[1]; py++)
142  {
143  const Vector3D volume_coords = px0_vol + e_y*py;
144 
145  // Evaluate volume at that position
146  const unsigned char val = this->evaluate(volume_coords);
147 
148  const double noise_val = noiseValue(noiseSigma, noiseMean);
149  const int noised_val = noise_val + val;
150  unsigned char final_val = this->constrainToUnsignedChar(noised_val);
151  // Store that value in the US slice
152  us_data[px + py*sliceDimension[0]] = final_val;
153 
154  }
155  }
156 
157  setDeepModified(us_frame);
158  return us_frame;
159 }
160 
161 vtkImageDataPtr cxSyntheticVolume::createEmptyMask(const Eigen::Array2i& sliceDimension) const
162 {
163  vtkImageDataPtr mask = vtkImageDataPtr::New();
164  mask->SetExtent(0, sliceDimension[0]-1, 0, sliceDimension[1]-1, 0, 0);
165  mask->AllocateScalars(VTK_UNSIGNED_CHAR, 1);
166  unsigned char* mask_data = (unsigned char*)mask->GetScalarPointer();
167  memset(mask_data, 1, sizeof(unsigned char)*sliceDimension[0]*sliceDimension[1]);
168  setDeepModified(mask);
169  return mask;
170 }
171 
172 unsigned char cxSyntheticVolume::constrainToUnsignedChar(const int val) const
173 {
174  if(val < 0)
175  {
176  return 0;
177  }
178  else if(val > 255)
179  {
180  return 255;
181  }
182  else
183  {
184  return (unsigned char)val;
185  }
186 }
187 
189 {
190  vtkImageDataPtr input = vol->getBaseVtkImageData();
191 
192  vtkImageDataPtr nominal = vtkImageDataPtr::New();
193  nominal->DeepCopy(input);
194  cx::ImagePtr nominal_img(new cx::Image("nominal", nominal));
195  nominal_img->get_rMd_History()->setRegistration(vol->get_rMd());
196  this->fillVolume(nominal_img);
197 
198  return calculateRMSError(input, nominal);
199 }
200 
202 {
203  vtkImageDataPtr raw = vol->getBaseVtkImageData();
204  cx::Transform3D rMd = vol->get_rMd();
205 
206  Eigen::Array3i dims = Eigen::Array3i(raw->GetDimensions());
207  Eigen::Array3d spacing = Eigen::Array3d(raw->GetSpacing());
208  unsigned char *pixels = (unsigned char*)raw->GetScalarPointer();
209 
210  for(int z = 0; z < dims[2]; ++z)
211  {
212  for(int y = 0; y < dims[1]; ++y)
213  {
214  for(int x = 0; x < dims[0]; ++x)
215  {
216  Vector3D p_d = Vector3D(x, y, z).array()*spacing;
217  int index = x + y*dims[0] + z*dims[1]*dims[0];
218  pixels[index] = this->evaluate(rMd.coord(p_d));
219  }
220  }
221  }
222  setDeepModified(raw);
223 }
224 
226 {
227  CX_ASSERT(Eigen::Array3i(a->GetDimensions()).isApprox(Eigen::Array3i(b->GetDimensions())));
228  CX_ASSERT(Eigen::Array3d(a->GetSpacing()).isApprox(Eigen::Array3d(b->GetSpacing())));
229 
230  float sse = 0.0f;
231  Eigen::Array3i dims = Eigen::Array3i(a->GetDimensions());
232  unsigned char *pa = (unsigned char*)a->GetScalarPointer();
233  unsigned char *pb = (unsigned char*)b->GetScalarPointer();
234 
235  for(int z = 0; z < dims[2]; ++z)
236  {
237  for(int y = 0; y < dims[1]; ++y)
238  {
239  for(int x = 0; x < dims[0]; ++x)
240  {
241  int index = x + y*dims[0] + z*dims[1]*dims[0];
242  float error = pa[index] - pb[index];
243  sse += error*error;
244  }
245  }
246  }
247 
248  return sqrt(sse/dims.prod());
249 }
250 
255 {
257  {
258  mSum = 0;
259  mThreshold = 2;
260  }
261 
262  void operator()(int x, int y, int z, unsigned char* val)
263  {
264  if (*val < mThreshold)
265  return;
266  mSum += *val;
267  }
268 
269  double getVolume()
270  {
271  return mSum;
272  }
273 
274 private:
275  double mSum;
276  unsigned char mThreshold;
277 };
278 
283 {
285  {
286  mSum = 0;
287  mWeight = Vector3D::Zero();
288  mThreshold = 2;
289  }
290 
291  void operator()(int x, int y, int z, unsigned char* val)
292  {
293  if (*val < mThreshold)
294  return;
295  mSum += *val;
296  mWeight += Vector3D(x,y,z) * (*val);
297  }
298 
300  {
301  return mWeight/mSum;
302  }
303 
304 private:
305  double mSum;
306  Vector3D mWeight;
307  unsigned char mThreshold;
308 };
309 
310 template<class FUNCTOR>
311 void applyFunctor(cx::ImagePtr image, FUNCTOR& func)
312 {
313  vtkImageDataPtr raw = image->getBaseVtkImageData();
314  Eigen::Array3i dims = Eigen::Array3i(raw->GetDimensions());
315  unsigned char *pixels = (unsigned char*)raw->GetScalarPointer();
316 
317  for(int z = 0; z < dims[2]; ++z)
318  {
319  for(int y = 0; y < dims[1]; ++y)
320  {
321  for(int x = 0; x < dims[0]; ++x)
322  {
323  int index = x + y*dims[0] + z*dims[1]*dims[0];
324  func(x,y,z, pixels+index);
325  }
326  }
327  }
328 }
329 
331 {
332  vtkImageDataPtr raw = image->getBaseVtkImageData();
333  cx::Transform3D rMd = image->get_rMd();
334  Eigen::Array3d spacing = Eigen::Array3d(raw->GetSpacing());
335 
336  CentroidFunctor func;
337  applyFunctor(image, func);
338  Vector3D ci = func.getCentroid().array() * spacing;
339  return rMd.coord(ci);
340 }
341 
343 {
344  MassFunctor func;
345  applyFunctor(image, func);
346  return func.getVolume();
347 }
348 
349 }
void fillVolume(cx::ImagePtr vol)
#define CX_ASSERT(statement)
Definition: cxLogger.h:116
One position with timestamp.
Transform3D Transform3D
Transform3D is a representation of an affine 3D transform.
double calculateMass(cx::ImagePtr image)
boost::shared_ptr< class Image > ImagePtr
Definition: cxDicomWidget.h:27
void operator()(int x, int y, int z, unsigned char *val)
virtual unsigned char evaluate(const Vector3D &p) const =0
void applyFunctor(cx::ImagePtr image, FUNCTOR &func)
cx::Vector3D calculateCentroid(cx::ImagePtr image)
Transform3D get_uMv() const
get transform from inverted image space v (origin in ul corner) to image space u. ...
A volumetric data set.
Definition: cxImage.h:45
double calculateRMSError(vtkImageDataPtr a, vtkImageDataPtr b)
Transform3D get_tMu() const
get transform from image space u to probe tool space t.
virtual ProcessedUSInputDataPtr sampleUsData(const std::vector< Transform3D > &planes_rMf, const Eigen::Array2f &pixelSpacing, const Eigen::Array2i &sliceDimension, const Transform3D &output_dMr, const double noiseSigma, const unsigned char noiseMean) const
unsigned char constrainToUnsignedChar(const int val) const
Eigen::Vector3d Vector3D
Vector3D is a representation of a point or vector in 3D.
Definition: cxVector3D.h:42
Definition of characteristics for an Ultrasound Probe Sector.
Vector3D getSpacing() const
double noiseValue(double noiseSigma, double noiseMean)
void setDeepModified(vtkImageDataPtr image)
vtkImageDataPtr createEmptyMask(const Eigen::Array2i &sliceDimension) const
Utility functions for drawing an US Probe sector.
Definition: cxProbeSector.h:38
virtual float computeRMSError(cx::ImagePtr vol)
boost::shared_ptr< class ProcessedUSInputData > ProcessedUSInputDataPtr
vtkSmartPointer< class vtkImageData > vtkImageDataPtr
#define M_PI
void setData(ProbeDefinition data)
Namespace for all CustusX production code.
void operator()(int x, int y, int z, unsigned char *val)