Fraxinus  2023.01.05-dev+develop.0da12
An IGT application
cxReconstructPreprocessor.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 
14 
15 #include "cxLogger.h"
16 #include "cxReconstructCore.h"
17 #include <vtkImageData.h>
18 #include <QFileInfo>
19 #include "cxTime.h"
20 #include "cxTypeConversions.h"
22 #include "cxVolumeHelpers.h"
24 #include "cxTimeKeeper.h"
25 #include "cxUSFrameData.h"
26 
28 #include "cxPatientModelService.h"
29 #include "cxMathUtils.h"
30 #include "cxPositionFilter.h"
31 
32 #include <QThread>
33 
34 namespace cx
35 {
36 
38  mInput(ReconstructCore::InputParams()),
39  mPatientModelService(patientModelService)
40 {
41  mMaxTimeDiff = 250; // TODO: Change default value for max allowed time difference between tracking and image time tags
42  // TODO Legg inn andre rekonstruksjonsparametre?
43  //
44  //
45 }
46 
48 {
49 }
50 
51 std::vector<ProcessedUSInputDataPtr> ReconstructPreprocessor::createProcessedInput(std::vector<bool> angio)
52 {
53 
54  std::vector<std::vector<vtkImageDataPtr> > frames = mFileData.mUsRaw->initializeFrames(angio);
55 
56  std::vector<ProcessedUSInputDataPtr> retval;
57 
58  for (unsigned i=0; i<angio.size(); ++i)
59  {
61  input.reset(new ProcessedUSInputData(frames[i],
62  mFileData.mFrames,
63  mFileData.getMask(),
64  mFileData.mFilename,
65  QFileInfo(mFileData.mFilename).completeBaseName() ));
66  CX_ASSERT(Eigen::Array3i(frames[i][0]->GetDimensions()).isApprox(Eigen::Array3i(mFileData.getMask()->GetDimensions())));
67  retval.push_back(input);
68  }
69  return retval;
70 }
71 
81 void ReconstructPreprocessor::calibrateTimeStamps(double offset, double scale)
82 {
83  for (unsigned int i = 0; i < mFileData.mPositions.size(); i++)
84  {
85  mFileData.mPositions[i].mTime = scale * mFileData.mPositions[i].mTime + offset;
86  }
87 }
88 
95 void ReconstructPreprocessor::alignTimeSeries()
96 {
97  report("Generate time calibration based on input time stamps.");
98  double framesSpan = mFileData.mFrames.back().mTime - mFileData.mFrames.front().mTime;
99  double positionsSpan = mFileData.mPositions.back().mTime - mFileData.mPositions.front().mTime;
100  double scale = framesSpan / positionsSpan;
101 
102  double offset = mFileData.mFrames.front().mTime - scale * mFileData.mPositions.front().mTime;
103 
104  this->calibrateTimeStamps(offset, scale);
105 }
106 
110 void ReconstructPreprocessor::cropInputData()
111 {
112  //IntBoundingBox3D
113  ProbeDefinition sector = mFileData.mProbeDefinition.mData;
114  IntBoundingBox3D cropbox(sector.getClipRect_p().begin());
115  cropbox = this->reduceCropboxToImageSize(cropbox, sector.getSize());
116  Eigen::Vector3i shift = cropbox.corner(0,0,0).cast<int>();
117  Eigen::Vector3i size = cropbox.range().cast<int>() + Eigen::Vector3i(1,1,0); // convert from extent format to size format by adding 1
118  mFileData.mUsRaw->setCropBox(cropbox);
119 
120  DoubleBoundingBox3D clipRect_p = sector.getClipRect_p();
121  Vector3D origin_p = sector.getOrigin_p();
122 
123  for (unsigned i=0; i<3; ++i)
124  {
125  clipRect_p[2*i] -= shift[i];
126  clipRect_p[2*i+1] -= shift[i];
127  origin_p[i] -= shift[i];
128  }
129 
130  sector.setClipRect_p(clipRect_p);
131  sector.setOrigin_p(origin_p);
132  sector.setSize(QSize(size[0], size[1]));
133  mFileData.mProbeDefinition.setData(sector);
134 }
135 
136 IntBoundingBox3D ReconstructPreprocessor::reduceCropboxToImageSize(IntBoundingBox3D cropbox, QSize size)
137 {
138  cropbox[0] = std::max(cropbox[0], 0);
139  cropbox[2] = std::max(cropbox[2], 0);
140  cropbox[4] = std::max(cropbox[4], 0);
141 
142  cropbox[1] = std::min(cropbox[1], size.width() - 1);
143  cropbox[3] = std::min(cropbox[3], size.height() - 1);
144 
145  return cropbox;
146 }
147 
157 void ReconstructPreprocessor::applyTimeCalibration()
158 {
159  double timeshift = mInput.mExtraTimeCalibration;
160  // The shift is on frames. The calibrate function applies to tracker positions,
161  // hence the positive sign. (real use: subtract from frame data)
162  // std::cout << "TIMESHIFT " << timeshift << std::endl;
163  // timeshift = -timeshift;
164  if (!similar(0.0, timeshift))
165  report("Applying reconstruction-time calibration to tracking data: " + qstring_cast(
166  timeshift) + "ms");
167  this->calibrateTimeStamps(timeshift, 1.0);
168 
169  // ignore calibrations
170  if (mInput.mAlignTimestamps)
171  {
172  this->alignTimeSeries();
173  }
174 }
175 
176 
178 {
179  RemoveDataType() : count(0), err(-1) {}
180  void add(double _err) { ++count; err=((err<0)?_err:std::min(_err, err)); }
181  unsigned count;
182  double err;
183 };
184 
185 void ReconstructPreprocessor::filterPositions()
186 {
187  unsigned filterStrength = mInput.mPosFilterStrength;
188  CX_LOG_DEBUG() << "Filter length specified " << filterStrength;
189  //PositionsPtr positions = new PositionsPtr(mFileData.mPositions);
190  PositionFilter positionFilter(filterStrength, mFileData.mPositions);
191  positionFilter.filterPositions();
192 }
193 
194 void ReconstructPreprocessor::positionThinning()
195 {
196  // If enabled, try to remove "suspect" data (large jumps etc.)
197  // Replace tracking positions that deviate greatly from neighbours with an interpolated value
198 }
199 
200 
207 void ReconstructPreprocessor::interpolatePositions()
208 {
209  mFileData.mUsRaw->resetRemovedFrames();
210  unsigned long startFrames = mFileData.mFrames.size();
211 
212  std::map<unsigned,RemoveDataType> removedData;
213 
214  for (unsigned i_frame = 0; i_frame < mFileData.mFrames.size();)
215  {
216  std::vector<TimedPosition>::iterator posIter;
217  posIter = lower_bound(mFileData.mPositions.begin(), mFileData.mPositions.end(), mFileData.mFrames[i_frame]);
218 
219  unsigned i_pos = unsigned(distance(mFileData.mPositions.begin(), posIter));
220  if (i_pos != 0)
221  i_pos--;
222 
223  if (i_pos >= mFileData.mPositions.size() - 1)
224  i_pos = unsigned(mFileData.mPositions.size()) - 2;
225 
226  // Remove frames too far from the positions
227  // Don't increment frame index since the index now points to the next element
228  double timeToPos1 = timeToPosition(i_frame, unsigned(i_pos));
229  double timeToPos2 = timeToPosition(i_frame, i_pos+1);
230  if ((timeToPos1 > mMaxTimeDiff) || (timeToPos2 > mMaxTimeDiff))
231  {
232  removedData[i_frame].add(std::max(timeToPos1, timeToPos2));
233 
234  mFileData.mFrames.erase(mFileData.mFrames.begin() + i_frame);
235  mFileData.mUsRaw->removeFrame(i_frame);
236  }
237  else
238  {
239  double t_delta_tracking = mFileData.mPositions[i_pos + 1].mTime - mFileData.mPositions[i_pos].mTime;
240  double t = 0;
241  if (!similar(t_delta_tracking, 0))
242  t = (mFileData.mFrames[i_frame].mTime - mFileData.mPositions[i_pos].mTime) / t_delta_tracking;
243  // mFileData.mFrames[i_frame].mPos = interpolate(mFileData.mPositions[i_pos].mPos, mFileData.mPositions[i_pos + 1].mPos, t);
244  mFileData.mFrames[i_frame].mPos = cx::USReconstructInputDataAlgorithm::slerpInterpolate(mFileData.mPositions[i_pos].mPos, mFileData.mPositions[i_pos + 1].mPos, t);
245  i_frame++;// Only increment if we didn't delete the frame
246  }
247  }
248 
249  unsigned removeCount=0;
250  for (std::map<unsigned,RemoveDataType>::iterator iter=removedData.begin(); iter!=removedData.end(); ++iter)
251  {
252  unsigned first = iter->first+removeCount;
253  unsigned last = first + iter->second.count-1;
254  report(QString("Removed input frame [%1-%2]. Time diff=%3").arg(first).arg(last).arg(iter->second.err, 0, 'f', 1));
255  removeCount += iter->second.count;
256  }
257 
258  double removed = double(startFrames - mFileData.mFrames.size()) / double(startFrames);
259  if (removed > 0.02)
260  {
261  double percent = removed * 100;
262  if (percent > 1)
263  {
264  reportWarning("Removed " + QString::number(percent, 'f', 1) + "% of the "+ qstring_cast(startFrames) + " frames.");
265  }
266  else
267  {
268  report("Removed " + QString::number(percent, 'f', 1) + "% of the " + qstring_cast(startFrames) + " frames.");
269  }
270  }
271 }
272 
273 
274 double ReconstructPreprocessor::timeToPosition(unsigned i_frame, unsigned i_pos)
275 {
276  return fabs(mFileData.mFrames[i_frame].mTime - mFileData.mPositions[i_pos].mTime);
277 }
278 
282 std::vector<Vector3D> ReconstructPreprocessor::generateInputRectangle()
283 {
284  std::vector<Vector3D> retval(4);
285  vtkImageDataPtr mask = mFileData.getMask();
286  if (!mask)
287  {
288  reportError("Reconstructer::generateInputRectangle() + requires mask");
289  return retval;
290  }
291  Eigen::Array3i dims = mFileData.mUsRaw->getDimensions();
292  Vector3D spacing = mFileData.mUsRaw->getSpacing();
293 
294  Eigen::Array3i maskDims(mask->GetDimensions());
295 
296  if (( maskDims[0]<dims[0] )||( maskDims[1]<dims[1] ))
297  reportError(QString("input data (%1) and mask (%2) dim mimatch")
298  .arg(qstring_cast(dims))
299  .arg(qstring_cast(maskDims)));
300 
301  int xmin = maskDims[0];
302  int xmax = 0;
303  int ymin = maskDims[1];
304  int ymax = 0;
305 
306  unsigned char* ptr = static_cast<unsigned char*> (mask->GetScalarPointer());
307  for (int x = 0; x < maskDims[0]; x++)
308  for (int y = 0; y < maskDims[1]; y++)
309  {
310  unsigned char val = ptr[x + y * maskDims[0]];
311  if (val != 0)
312  {
313  xmin = std::min(xmin, x);
314  ymin = std::min(ymin, y);
315  xmax = std::max(xmax, x);
316  ymax = std::max(ymax, y);
317  }
318  }
319 
320  //Reduce the output volume by reducing the mask when determining output volume size
321  double red = mInput.mMaskReduce;
322  int reduceX = int((xmax - xmin) * (red / 100));
323  int reduceY = int((ymax - ymin) * (red / 100));
324 
325  xmin += reduceX;
326  xmax -= reduceX;
327  ymin += reduceY;
328  ymax -= reduceY;
329 
330  retval[0] = Vector3D(xmin * spacing[0], ymin * spacing[1], 0);
331  retval[1] = Vector3D(xmax * spacing[0], ymin * spacing[1], 0);
332  retval[2] = Vector3D(xmin * spacing[0], ymax * spacing[1], 0);
333  retval[3] = Vector3D(xmax * spacing[0], ymax * spacing[1], 0);
334 
335  return retval;
336 }
337 
344 Transform3D ReconstructPreprocessor::applyOutputOrientation()
345 {
346  Transform3D prMdd = Transform3D::Identity();
347 
348  if (mInput.mOrientation == "PatientReference")
349  {
350  // identity
351  }
352  else if (mInput.mOrientation == "MiddleFrame")
353  {
354  prMdd = mFileData.mFrames[mFileData.mFrames.size() / 2].mPos;
355  }
356  else
357  {
358  reportError("no orientation algorithm selected in reconstruction");
359  }
360 
361  // apply the selected orientation to the frames.
362  Transform3D ddMpr = prMdd.inv();
363  for (unsigned i = 0; i < mFileData.mFrames.size(); i++)
364  {
365  // mPos = prMu
366  mFileData.mFrames[i].mPos = ddMpr * mFileData.mFrames[i].mPos;
367  // mPos = ddMu
368  }
369 
370  return prMdd;
371 }
372 
382 void ReconstructPreprocessor::findExtentAndOutputTransform()
383 {
384  if (mFileData.mFrames.empty())
385  return;
386  // A first guess for d'Mu with correct orientation
387  Transform3D prMdd = this->applyOutputOrientation();
388  //mFrames[i].mPos = d'Mu, d' = only rotation
389 
390  // Find extent of all frames as a point cloud
391  std::vector<Vector3D> inputRect = this->generateInputRectangle();
392  std::vector<Vector3D> outputRect;
393  for (unsigned slice = 0; slice < mFileData.mFrames.size(); slice++)
394  {
395  Transform3D dMu = mFileData.mFrames[slice].mPos;
396  for (unsigned i = 0; i < inputRect.size(); i++)
397  {
398  outputRect.push_back(dMu.coord(inputRect[i]));
399  }
400  }
401 
403 
404  // Translate dMu to output volume origo
405  Transform3D T_origo = createTransformTranslate(extent.corner(0, 0, 0));
406  Transform3D prMd = prMdd * T_origo; // transform from output space to patref, use when storing volume.
407  for (unsigned i = 0; i < mFileData.mFrames.size(); i++)
408  {
409  mFileData.mFrames[i].mPos = T_origo.inv() * mFileData.mFrames[i].mPos;
410  }
411 
412  // Calculate optimal output image spacing and dimensions based on US frame spacing
413  double inputSpacing = std::min(mFileData.mUsRaw->getSpacing()[0], mFileData.mUsRaw->getSpacing()[1]);
414  mOutputVolumeParams = OutputVolumeParams(extent, inputSpacing, mInput.mMaxOutputVolumeSize);
415 
416  mOutputVolumeParams.set_rMd((mPatientModelService->get_rMpr()) * prMd);
417 }
418 
422 void ReconstructPreprocessor::updateFromOriginalFileData()
423 {
424  // uncomment to test cropping of data before reconstructing
425  this->cropInputData();
426 
427  // Only use this if the time stamps have different formats
428  // The function assumes that both lists of time stamps start at the same time,
429  // and that is not completely correct.
430  //this->calibrateTimeStamps();
431  // Use the time calibration from the acquisition module
432  //this->calibrateTimeStamps(0.0, 1.0);
433  this->applyTimeCalibration();
434 
435  // Smooth tracking data before further processing
436  // User preferences apply
437  //this->positionThinning(); //Do thinning before filtering if enabled
438  this->filterPositions();
439 
441  //mPos (in mPositions) is now prMu
442 
443  this->interpolatePositions();
444  // mFrames: now mPos as prMu
445 
446  if (mFileData.mFrames.empty()) // if all positions are filtered out
447  return;
448 
449  this->findExtentAndOutputTransform();
450  //mPos in mFrames is now dMu
451 }
452 
454 {
455  mInput = input;
456  mFileData = fileData;
457  this->updateFromOriginalFileData();
458 }
459 
460 
461 } /* namespace cx */
QString qstring_cast(const T &val)
void reportError(QString msg)
Definition: cxLogger.cpp:71
Vector3D corner(int x, int y, int z) const
#define CX_ASSERT(statement)
Definition: cxLogger.h:116
Transform3D Transform3D
Transform3D is a representation of an affine 3D transform.
std::vector< TimedPosition > mFrames
Helper struct for sending and controlling output volume properties.
static void transformTrackingPositionsTo_prMu(USReconstructInputData *data)
static DoubleBoundingBox3D fromCloud(std::vector< Vector3D > cloud)
Vector3D getOrigin_p() const
void setOrigin_p(Vector3D origin_p)
ReconstructPreprocessor(PatientModelServicePtr patientModelService)
void reportWarning(QString msg)
Definition: cxLogger.cpp:70
void setClipRect_p(DoubleBoundingBox3D clipRect_p)
boost::shared_ptr< class PatientModelService > PatientModelServicePtr
Transform3D createTransformTranslate(const Vector3D &translation)
std::vector< ProcessedUSInputDataPtr > createProcessedInput(std::vector< bool > angio)
Representation of an integer bounding box in 3D. The data are stored as {xmin,xmax,ymin,ymax,zmin,zmax}, in order to simplify communication with vtk.
Representation of a floating-point bounding box in 3D. The data are stored as {xmin,xmax,ymin,ymax,zmin,zmax}, in order to simplify communication with vtk.
std::vector< TimedPosition > mPositions
ProbeDefinition mData
Definition: cxProbeSector.h:54
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.
Eigen::Vector3i corner(int x, int y, int z) const
#define CX_LOG_DEBUG
Definition: cxLogger.h:95
void report(QString msg)
Definition: cxLogger.cpp:69
bool similar(const CameraInfo &lhs, const CameraInfo &rhs, double tol)
DoubleBoundingBox3D getClipRect_p() const
static Transform3D slerpInterpolate(const Transform3D &a, const Transform3D &b, double t)
QString mFilename
filename used for current data read
boost::shared_ptr< class ProcessedUSInputData > ProcessedUSInputDataPtr
vtkSmartPointer< class vtkImageData > vtkImageDataPtr
USFrameDataPtr mUsRaw
All imported US data frames with pointers to each frame.
void initialize(ReconstructCore::InputParams input, USReconstructInputData fileData)
void setSize(QSize size)
Algorithm part of reconstruction - no dependencies on parameter classes.
void setData(ProbeDefinition data)
Namespace for all CustusX production code.