Fraxinus  16.5.0-fx-rc1
An IGT application
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
cxPNNReconstructionMethodService.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 
34 
35 #include <QFileInfo>
36 #include "cxLogger.h"
37 #include "cxTypeConversions.h"
38 #include "cxVolumeHelpers.h"
39 #include "cxTimeKeeper.h"
40 #include "cxUSFrameData.h"
41 #include <vtkImageData.h>
42 #include "cxImage.h"
43 #include "cxDoubleProperty.h"
44 
45 namespace cx
46 {
48 {
49 }
50 
52 {}
53 
55 {
56  return "PNN";
57 }
58 
59 std::vector<PropertyPtr> PNNReconstructionMethodService::getSettings(QDomElement root)
60 {
61  std::vector<PropertyPtr> retval;
62  retval.push_back(this->getInterpolationStepsOption(root));
63  return retval;
64 }
65 
66 DoublePropertyPtr PNNReconstructionMethodService::getInterpolationStepsOption(QDomElement root)
67 {
68  DoublePropertyPtr retval;
69  retval = DoubleProperty::initialize("interpolationSteps", "Distance (voxels)",
70  "Interpolation steps in voxels", 3, DoubleRange(1, 10, 1), 0, root);
71  return retval;
72 }
73 
74 void optimizedCoordTransform(Vector3D* p, boost::array<double, 16> tt)
75 {
76  double* t = tt.begin();
77  double x = (*p)[0];
78  double y = (*p)[1];
79  double z = (*p)[2];
80  (*p)[0] = t[0] * x + t[1] * y + t[2] * z + t[3];
81  (*p)[1] = t[4] * x + t[5] * y + t[6] * z + t[7];
82  (*p)[2] = t[8] * x + t[9] * y + t[10] * z + t[11];
83 }
84 
86  vtkImageDataPtr outputData, QDomElement settings)
87 {
88  input->validate();
89 
90  std::vector<TimedPosition> frameInfo = input->getFrames();
91  if (frameInfo.empty())
92  return false;
93  if (input->getDimensions()[2]==0)
94  return false;
95 
96  vtkImageDataPtr target = outputData;
97 
98  Eigen::Array3i inputDims = input->getDimensions();
99 
100  Eigen::Array3i targetDims(target->GetDimensions());
101  Vector3D targetSpacing(target->GetSpacing());
102 
103  //Create temporary volume
104  vtkImageDataPtr tempOutput = generateVtkImageData(targetDims, targetSpacing, 0);
105  ImagePtr tempOutputData = ImagePtr(new Image("tempOutput", tempOutput, "tempOutput"));
106 
107  int* outputDims = tempOutput->GetDimensions();
108  //int* outputDims = target->GetDimensions();
109 
110  if (inputDims[2] != static_cast<int> (frameInfo.size()))
111  reportWarning("inputDims[2] != frameInfo.size()" + qstring_cast(inputDims[2]) + " != "
112  + qstring_cast(frameInfo.size()));
113 
114  Vector3D inputSpacing(input->getSpacing());
115  Vector3D outputSpacing(tempOutput->GetSpacing());
116 
117  //Get raw data pointers
118  unsigned char *outputPointer = static_cast<unsigned char*> (tempOutput->GetScalarPointer());
119  unsigned char* maskPointer = static_cast<unsigned char*> (input->getMask()->GetScalarPointer());
120 
121  // Traverse all input pixels
122  for (int record = 0; record < inputDims[2]; record++)
123  {
124  unsigned char *inputPointer = input->getFrame(record);
125  boost::array<double, 16> recordTransform = frameInfo[record].mPos.flatten();
126 
127  for (int beam = 0; beam < inputDims[0]; beam++)
128  {
129  for (int sample = 0; sample < inputDims[1]; sample++)
130  {
131  if (!validPixel(beam, sample, inputDims, maskPointer))
132  continue;
133  Vector3D inputPoint(beam * inputSpacing[0], sample * inputSpacing[1], 0.0);
134  Vector3D outputPoint = inputPoint;
135  optimizedCoordTransform(&outputPoint, recordTransform);
136  int outputVoxelX = static_cast<int> ((outputPoint[0] / outputSpacing[0]) + 0.5);
137  int outputVoxelY = static_cast<int> ((outputPoint[1] / outputSpacing[1]) + 0.5);
138  int outputVoxelZ = static_cast<int> ((outputPoint[2] / outputSpacing[2]) + 0.5);
139 
140  if (validVoxel(outputVoxelX, outputVoxelY, outputVoxelZ, outputDims))
141  {
142  int outputIndex = outputVoxelX + outputVoxelY * outputDims[0] + outputVoxelZ * outputDims[0]
143  * outputDims[1];
144  int inputIndex = beam + sample * inputDims[0];
145 
146  // assign the max value found from all frames hitting this voxel. This removes black areas where (some of) multiple sweeps contains shadows.
147  outputPointer[outputIndex] = std::max<unsigned char>(inputPointer[inputIndex], outputPointer[outputIndex]);
148  // set minimum intensity value to 1. This separates "zero intensity" from "no intensity".
149  outputPointer[outputIndex] = std::max<unsigned char>(inputPointer[inputIndex], 1); //
150  }//validVoxel
151 
152  }//sample
153  }//beam
154  }//record
155 
156  // Fill holes
157  this->interpolate(tempOutputData, outputData, settings);
158 
159  setDeepModified(outputData);
160  return true;
161 }
162 
163 namespace
164 {
167 inline int getIndex_z_last(int x, int y, int z, const Eigen::Array3i& dim)
168 {
169  return x + y*dim[0] + z*dim[0]*dim[1];
170 }
173 inline int getIndex_x_last(int y, int z, int x, const Eigen::Array3i& dim)
174 {
175  return x + y*dim[0] + z*dim[0]*dim[1];
176 }
179 inline int getIndex_y_last(int z, int x, int y, const Eigen::Array3i& dim)
180 {
181  return x + y*dim[0] + z*dim[0]*dim[1];
182 }
183 
184 
190 template <class FUNCTION>
191 void maskAlongDim(int a_dim, int b_dim, int c_dim, const Eigen::Array3i& dim, unsigned char *inputPtr, unsigned char *maskPtr, FUNCTION getIndex)
192 {
193  for (int a = 0; a < a_dim; a++)
194  {
195  for (int b = 0; b < b_dim; b++)
196  {
197  int start = c_dim;
198  int stop = -1;
199  for (int c = 0; c < c_dim; c++)
200  {
201  int index = getIndex(a, b, c, dim);
202  if (inputPtr[index]>0)
203  {
204  start = c;
205  break;
206  }
207  }
208  for (int c = c_dim-1; c >=0; c--)
209  {
210  int index = getIndex(a, b, c, dim);
211  if (inputPtr[index]>0)
212  {
213  stop = c;
214  break;
215  }
216  }
217  for (int c = start; c <= stop; c++)
218  {
219  int index = getIndex(a, b, c, dim);
220  maskPtr[index] = 1;
221  }
222  }
223  }
224 }
225 } // unnamed namespace
226 
235 vtkImageDataPtr PNNReconstructionMethodService::createMask(vtkImageDataPtr inputData)
236 {
237  Eigen::Array3i dim(inputData->GetDimensions());
238  Vector3D spacing(inputData->GetSpacing());
239  vtkImageDataPtr mask = generateVtkImageData(dim, spacing, 0);
240  unsigned char *inputPtr = static_cast<unsigned char*> (inputData->GetScalarPointer());
241  unsigned char *maskPtr = static_cast<unsigned char*> (mask->GetScalarPointer());
242 
243  // mask along all 3 dimensions
244  maskAlongDim(dim[0], dim[1], dim[2], dim, inputPtr, maskPtr, &getIndex_z_last);
245  maskAlongDim(dim[1], dim[2], dim[0], dim, inputPtr, maskPtr, &getIndex_x_last);
246  maskAlongDim(dim[2], dim[0], dim[1], dim, inputPtr, maskPtr, &getIndex_y_last);
247 
248  return mask;
249 }
250 
251 void PNNReconstructionMethodService::interpolate(ImagePtr inputData, vtkImageDataPtr outputData, QDomElement settings)
252 {
253  TimeKeeper timer;
254  DoublePropertyPtr interpolationStepsOption = this->getInterpolationStepsOption(settings);
255  int interpolationSteps = static_cast<int> (interpolationStepsOption->getValue());
256 
257  vtkImageDataPtr input = inputData->getBaseVtkImageData();
258  vtkImageDataPtr output = outputData;
259  vtkImageDataPtr mask = this->createMask(input);
260 
261  Eigen::Array3i outputDims(output->GetDimensions());
262 
263  Eigen::Array3i inputDims(input->GetDimensions());
264 
265  unsigned char *inputPointer = static_cast<unsigned char*> (input->GetScalarPointer());
266  unsigned char *outputPointer = static_cast<unsigned char*> (output->GetScalarPointer());
267  unsigned char *maskPointer = static_cast<unsigned char*> (mask->GetScalarPointer());
268 
269  if ((outputDims[0] != inputDims[0]) || (outputDims[1] != inputDims[1]) || (outputDims[2] != inputDims[2]))
270  reportWarning("outputDims != inputDims. output: " + qstring_cast(outputDims[0]) + " "
271  + qstring_cast(outputDims[1]) + " " + qstring_cast(outputDims[2]) + " input: " + qstring_cast(inputDims[0])
272  + " " + qstring_cast(inputDims[1]) + " " + qstring_cast(inputDims[2]));
273 
274  int total = outputDims[0] * outputDims[1] * outputDims[2];
275  int removed = 0;
276  int ignored = 0;
277  // Traverse all voxels
278  for (int x = 0; x < outputDims[0]; x++)
279  {
280  for (int y = 0; y < outputDims[1]; y++)
281  {
282  for (int z = 0; z < outputDims[2]; z++)
283  {
284  int outputIndex = x + y * outputDims[0] + z * outputDims[0] * outputDims[1];
285 
286  // ignore if outside volume of interest
287  if (maskPointer[outputIndex]==0)
288  {
289  removed++;
290  }
291  // copy if value already exists
292  else if (inputPointer[outputIndex]>0)
293  {
294  outputPointer[outputIndex] = inputPointer[outputIndex];
295  ignored++;
296  }
297  // fill hole otherwise (empty space within the volume)
298  else
299  {
300  this->fillHole(inputPointer, outputPointer, x, y, z, outputDims, interpolationSteps);
301  }
302  }//z
303  }//y
304  }//x
305 
306  int valid = 100*double(ignored)/double(total);
307  int outside = 100*double(removed)/double(total);
308  int holes = 100*double(total-ignored-removed)/double(total);
309  reportDebug(
310  QString("PNN: Size: %1Mb, Valid voxels: %2\%, Outside mask: %3\% Filled holes [steps=%4, %5s]: %6\%")
311  .arg(total/1024/1024)
312  .arg(valid)
313  .arg(outside)
314  .arg(interpolationSteps)
315  .arg(timer.getElapsedSecondsAsString())
316  .arg(holes));
317 }
318 
323 void PNNReconstructionMethodService::fillHole(unsigned char *inputPointer, unsigned char *outputPointer, int x, int y, int z, const Eigen::Array3i& dim, int interpolationSteps)
324 {
325  int outputIndex = x + y * dim[0] + z * dim[0] * dim[1];
326  bool interpolated = false;
327  int localArea = 0;
328 
329  int count = 0;
330  double tempVal = 0;
331 
332  do
333  {
334  for (int i = -localArea; i < localArea + 1; i++)
335  {
336  for (int j = -localArea; j < localArea + 1; j++)
337  {
338  for (int k = -localArea; k < localArea + 1; k++)
339  {
340  int localIndex = outputIndex + i + j*dim[0] + k*dim[0]*dim[1];
341 
342  if (validVoxel(x + i, y + j, z + k, dim.data()) && inputPointer[localIndex] > 0.1)
343  {
344  tempVal += inputPointer[localIndex];
345  count++;
346  }
347  }//local voxel area
348  }
349  }
350 
351  if (count > 0)
352  {
353  interpolated = true;
354  if (tempVal == 0)
355  {
356  // keep noneness of index
357  }
358  else
359  {
360  outputPointer[outputIndex] = static_cast<int> ((tempVal / count) + 0.5);
361  outputPointer[outputIndex] = std::max<unsigned char>(1, outputPointer[outputIndex]);
362  }
363  }
364 
365  localArea++;
366 
367  } while (localArea <= interpolationSteps && !interpolated);
368 }
369 
370 }//namespace
QString qstring_cast(const T &val)
Utility class for describing a bounded numeric range.
Definition: cxDoubleRange.h:53
boost::shared_ptr< class Image > ImagePtr
Definition: cxDicomWidget.h:48
virtual bool reconstruct(ProcessedUSInputDataPtr input, vtkImageDataPtr outputData, QDomElement settings)
void reportWarning(QString msg)
Definition: cxLogger.cpp:91
A volumetric data set.
Definition: cxImage.h:66
Settings * settings()
Shortcut for accessing the settings instance.
Definition: cxSettings.cpp:42
vtkImageDataPtr generateVtkImageData(Eigen::Array3i dim, Vector3D spacing, const unsigned char initValue, int components)
virtual std::vector< PropertyPtr > getSettings(QDomElement root)
boost::shared_ptr< class DoubleProperty > DoublePropertyPtr
Eigen::Vector3d Vector3D
Vector3D is a representation of a point or vector in 3D.
Definition: cxVector3D.h:63
void setDeepModified(vtkImageDataPtr image)
static DoublePropertyPtr initialize(const QString &uid, QString name, QString help, double value, DoubleRange range, int decimals, QDomNode root=QDomNode())
boost::shared_ptr< class ProcessedUSInputData > ProcessedUSInputDataPtr
void optimizedCoordTransform(Vector3D *p, boost::array< double, 16 > tt)
vtkSmartPointer< class vtkImageData > vtkImageDataPtr
void reportDebug(QString msg)
Definition: cxLogger.cpp:89