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