Fraxinus  17.12-rc4
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) 2008-2014, SINTEF Department of Medical Technology
5 All rights reserved.
6 
7 Redistribution and use in source and binary forms, with or without
8 modification, are permitted provided that the following conditions are met:
9 
10 1. Redistributions of source code must retain the above copyright notice,
11  this list of conditions and the following disclaimer.
12 
13 2. Redistributions in binary form must reproduce the above copyright notice,
14  this list of conditions and the following disclaimer in the documentation
15  and/or other materials provided with the distribution.
16 
17 3. Neither the name of the copyright holder nor the names of its contributors
18  may be used to endorse or promote products derived from this software
19  without specific prior written permission.
20 
21 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
22 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
23 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
24 DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
25 FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
26 DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
27 SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
28 CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
29 OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
30 OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
31 =========================================================================*/
32 
33 #include "cxSyntheticVolume.h"
34 #include "vtkImageData.h"
35 #include "cxImage.h"
36 #include <cstdlib>
37 #include <time.h>
38 #include "cxTypeConversions.h"
39 #include <QTime>
40 #include "cxLogger.h"
42 #include "cxVolumeHelpers.h"
43 
44 
45 double noiseValue(double noiseSigma,
46  double noiseMean)
47 {
48  double random_value_1 = (rand()+1.0)/(RAND_MAX+1.0);
49  double random_value_2 = (rand()+1.0)/(RAND_MAX+1.0);
50 
51  double random_normal = sqrt(-2*log(random_value_1)) * cos(2*M_PI*random_value_2);
52 
53  return random_normal*noiseSigma + noiseMean;
54 }
55 
56 namespace cx {
57 
58 
60 cxSyntheticVolume::sampleUsData(const std::vector<Transform3D>& planes_rMt,
61  const ProbeDefinition& probe,
62  const Transform3D& output_dMr,
63  const double noiseSigma,
64  const unsigned char noiseMean) const
65 {
66  cx::ProbeSector sector;
67  sector.setData(probe);
68 
69  std::vector<Transform3D> planes_rMf(planes_rMt.size());
70  for (unsigned i=0; i<planes_rMt.size(); ++i)
71  planes_rMf[i] = planes_rMt[i] * sector.get_tMu() * sector.get_uMv();
72 
73  Eigen::Array2f pixelSpacing = probe.getSpacing().block(0,0,2,1).cast<float>();
74  Eigen::Array2i sliceDimension(probe.getSize().width(), probe.getSize().height());
75 
76  return this->sampleUsData(planes_rMf,
77  pixelSpacing,
78  sliceDimension,
79  output_dMr,
80  noiseSigma, noiseMean);
81 }
82 
85  const ProbeDefinition& probe,
86  const double noiseSigma,
87  const unsigned char noiseMean) const
88 {
89  cx::ProbeSector sector;
90  sector.setData(probe);
91 
92  Transform3D rMf = plane_rMt * sector.get_tMu() * sector.get_uMv();
93  Eigen::Array2f pixelSpacing = probe.getSpacing().block(0,0,2,1).cast<float>();
94  Eigen::Array2i sliceDimension(probe.getSize().width(), probe.getSize().height());
95 
96  return this->sampleUsData(rMf,
97  pixelSpacing,
98  sliceDimension,
99  noiseSigma, noiseMean);
100 }
101 
103 cxSyntheticVolume::sampleUsData(const std::vector<Transform3D>& planes_rMf,
104  const Eigen::Array2f& pixelSpacing,
105  const Eigen::Array2i& sliceDimension,
106  const Transform3D& output_dMr,
107  const double noiseSigma, const unsigned char noiseMean) const
108 {
109  std::vector<TimedPosition> positions;
110  std::vector<vtkImageDataPtr> images;
111  // For each plane
112  for(std::vector<Transform3D>::const_iterator i = planes_rMf.begin();
113  planes_rMf.end() != i;
114  ++i)
115  {
116  const Transform3D rMf = *i;
117 
118  vtkImageDataPtr us_frame;
119  us_frame = this->sampleUsData(rMf, pixelSpacing, sliceDimension, noiseSigma, noiseMean);
120 
121  // Build the TimedPosition for this frame
122  TimedPosition t;
123  t.mTime = i - planes_rMf.begin();
124  t.mPos = output_dMr*rMf;
125 
126  positions.push_back(t);
127  images.push_back(us_frame);
128  }
129 
130  vtkImageDataPtr mask = this->createEmptyMask(sliceDimension);
131 // std::cout << "elapsed: " << time.elapsed() << std::endl;
132 
134  ret.reset(new ProcessedUSInputData(images, positions, mask, "VIRTUAL_DATA", "VIRTUAL_DATA_"));
135  return ret;
136 }
137 
140  const Eigen::Array2f& pixelSpacing,
141  const Eigen::Array2i& sliceDimension,
142  const double noiseSigma, const unsigned char noiseMean) const
143 {
144 
145 
146  const Vector3D p0 = rMf.coord(Vector3D(0,0,0));
147  const Vector3D e_x = rMf.vector(Vector3D(pixelSpacing[0],0,0));
148  const Vector3D e_y = rMf.vector(Vector3D(0,pixelSpacing[1],0));
149 
150  vtkImageDataPtr us_frame = vtkImageDataPtr::New();
151  us_frame->SetExtent(0, sliceDimension[0]-1, 0, sliceDimension[1]-1, 0, 0);
152  us_frame->SetSpacing(pixelSpacing[0], pixelSpacing[1], 0.0);
153  us_frame->AllocateScalars(VTK_UNSIGNED_CHAR, 1);
154 
155  unsigned char* us_data = (unsigned char*)us_frame->GetScalarPointer();
156  // For each pixel on that plane
157  for(unsigned int px = 0; px < sliceDimension[0]; px++)
158  {
159  // optimization: use transformed pixel vectors
160  const Vector3D px0_vol = p0 + e_x*px;
161 
162  for(unsigned int py = 0; py < sliceDimension[1]; py++)
163  {
164  const Vector3D volume_coords = px0_vol + e_y*py;
165 
166  // Evaluate volume at that position
167  const unsigned char val = this->evaluate(volume_coords);
168 
169  const double noise_val = noiseValue(noiseSigma, noiseMean);
170  const int noised_val = noise_val + val;
171  unsigned char final_val = this->constrainToUnsignedChar(noised_val);
172  // Store that value in the US slice
173  us_data[px + py*sliceDimension[0]] = final_val;
174 
175  }
176  }
177 
178  setDeepModified(us_frame);
179  return us_frame;
180 }
181 
182 vtkImageDataPtr cxSyntheticVolume::createEmptyMask(const Eigen::Array2i& sliceDimension) const
183 {
184  vtkImageDataPtr mask = vtkImageDataPtr::New();
185  mask->SetExtent(0, sliceDimension[0]-1, 0, sliceDimension[1]-1, 0, 0);
186  mask->AllocateScalars(VTK_UNSIGNED_CHAR, 1);
187  unsigned char* mask_data = (unsigned char*)mask->GetScalarPointer();
188  memset(mask_data, 1, sizeof(unsigned char)*sliceDimension[0]*sliceDimension[1]);
189  setDeepModified(mask);
190  return mask;
191 }
192 
193 unsigned char cxSyntheticVolume::constrainToUnsignedChar(const int val) const
194 {
195  if(val < 0)
196  {
197  return 0;
198  }
199  else if(val > 255)
200  {
201  return 255;
202  }
203  else
204  {
205  return (unsigned char)val;
206  }
207 }
208 
210 {
211  vtkImageDataPtr input = vol->getBaseVtkImageData();
212 
213  vtkImageDataPtr nominal = vtkImageDataPtr::New();
214  nominal->DeepCopy(input);
215  cx::ImagePtr nominal_img(new cx::Image("nominal", nominal));
216  nominal_img->get_rMd_History()->setRegistration(vol->get_rMd());
217  this->fillVolume(nominal_img);
218 
219  return calculateRMSError(input, nominal);
220 }
221 
223 {
224  vtkImageDataPtr raw = vol->getBaseVtkImageData();
225  cx::Transform3D rMd = vol->get_rMd();
226 
227  Eigen::Array3i dims = Eigen::Array3i(raw->GetDimensions());
228  Eigen::Array3d spacing = Eigen::Array3d(raw->GetSpacing());
229  unsigned char *pixels = (unsigned char*)raw->GetScalarPointer();
230 
231  for(int z = 0; z < dims[2]; ++z)
232  {
233  for(int y = 0; y < dims[1]; ++y)
234  {
235  for(int x = 0; x < dims[0]; ++x)
236  {
237  Vector3D p_d = Vector3D(x, y, z).array()*spacing;
238  int index = x + y*dims[0] + z*dims[1]*dims[0];
239  pixels[index] = this->evaluate(rMd.coord(p_d));
240  }
241  }
242  }
243  setDeepModified(raw);
244 }
245 
247 {
248  CX_ASSERT(Eigen::Array3i(a->GetDimensions()).isApprox(Eigen::Array3i(b->GetDimensions())));
249  CX_ASSERT(Eigen::Array3d(a->GetSpacing()).isApprox(Eigen::Array3d(b->GetSpacing())));
250 
251  float sse = 0.0f;
252  Eigen::Array3i dims = Eigen::Array3i(a->GetDimensions());
253  unsigned char *pa = (unsigned char*)a->GetScalarPointer();
254  unsigned char *pb = (unsigned char*)b->GetScalarPointer();
255 
256  for(int z = 0; z < dims[2]; ++z)
257  {
258  for(int y = 0; y < dims[1]; ++y)
259  {
260  for(int x = 0; x < dims[0]; ++x)
261  {
262  int index = x + y*dims[0] + z*dims[1]*dims[0];
263  float error = pa[index] - pb[index];
264  sse += error*error;
265  }
266  }
267  }
268 
269  return sqrt(sse/dims.prod());
270 }
271 
276 {
278  {
279  mSum = 0;
280  mThreshold = 2;
281  }
282 
283  void operator()(int x, int y, int z, unsigned char* val)
284  {
285  if (*val < mThreshold)
286  return;
287  mSum += *val;
288  }
289 
290  double getVolume()
291  {
292  return mSum;
293  }
294 
295 private:
296  double mSum;
297  unsigned char mThreshold;
298 };
299 
304 {
306  {
307  mSum = 0;
308  mWeight = Vector3D::Zero();
309  mThreshold = 2;
310  }
311 
312  void operator()(int x, int y, int z, unsigned char* val)
313  {
314  if (*val < mThreshold)
315  return;
316  mSum += *val;
317  mWeight += Vector3D(x,y,z) * (*val);
318  }
319 
321  {
322  return mWeight/mSum;
323  }
324 
325 private:
326  double mSum;
327  Vector3D mWeight;
328  unsigned char mThreshold;
329 };
330 
331 template<class FUNCTOR>
332 void applyFunctor(cx::ImagePtr image, FUNCTOR& func)
333 {
334  vtkImageDataPtr raw = image->getBaseVtkImageData();
335  Eigen::Array3i dims = Eigen::Array3i(raw->GetDimensions());
336  unsigned char *pixels = (unsigned char*)raw->GetScalarPointer();
337 
338  for(int z = 0; z < dims[2]; ++z)
339  {
340  for(int y = 0; y < dims[1]; ++y)
341  {
342  for(int x = 0; x < dims[0]; ++x)
343  {
344  int index = x + y*dims[0] + z*dims[1]*dims[0];
345  func(x,y,z, pixels+index);
346  }
347  }
348  }
349 }
350 
352 {
353  vtkImageDataPtr raw = image->getBaseVtkImageData();
354  cx::Transform3D rMd = image->get_rMd();
355  Eigen::Array3d spacing = Eigen::Array3d(raw->GetSpacing());
356 
357  CentroidFunctor func;
358  applyFunctor(image, func);
359  Vector3D ci = func.getCentroid().array() * spacing;
360  return rMd.coord(ci);
361 }
362 
364 {
365  MassFunctor func;
366  applyFunctor(image, func);
367  return func.getVolume();
368 }
369 
370 }
void fillVolume(cx::ImagePtr vol)
#define CX_ASSERT(statement)
Definition: cxLogger.h:137
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:48
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:66
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:63
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:59
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)