Fraxinus  16.5.0-fx-rc9
An IGT application
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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) 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 
35 
36 #include "cxLogger.h"
37 #include "cxReconstructCore.h"
38 #include <vtkImageData.h>
39 #include <QFileInfo>
40 #include "cxTime.h"
41 #include "cxTypeConversions.h"
43 #include "cxVolumeHelpers.h"
45 #include "cxTimeKeeper.h"
46 #include "cxUSFrameData.h"
47 
49 #include "cxPatientModelService.h"
50 
51 namespace cx
52 {
53 
55  mInput(ReconstructCore::InputParams()),
56  mPatientModelService(patientModelService)
57 {
58  mMaxTimeDiff = 100; // TODO: Change default value for max allowed time difference between tracking and image time tags
59 }
60 
62 {
63 }
64 
65 std::vector<ProcessedUSInputDataPtr> ReconstructPreprocessor::createProcessedInput(std::vector<bool> angio)
66 {
67 
68  std::vector<std::vector<vtkImageDataPtr> > frames = mFileData.mUsRaw->initializeFrames(angio);
69 
70  std::vector<ProcessedUSInputDataPtr> retval;
71 
72  for (unsigned i=0; i<angio.size(); ++i)
73  {
75  input.reset(new ProcessedUSInputData(frames[i],
76  mFileData.mFrames,
77  mFileData.getMask(),
78  mFileData.mFilename,
79  QFileInfo(mFileData.mFilename).completeBaseName() ));
80  CX_ASSERT(Eigen::Array3i(frames[i][0]->GetDimensions()).isApprox(Eigen::Array3i(mFileData.getMask()->GetDimensions())));
81  retval.push_back(input);
82  }
83  return retval;
84 }
85 
86 namespace
87 {
88 bool within(int x, int min, int max)
89 {
90  return (x >= min) && (x <= max);
91 }
92 }
93 
103 void ReconstructPreprocessor::calibrateTimeStamps(double offset, double scale)
104 {
105  for (unsigned int i = 0; i < mFileData.mPositions.size(); i++)
106  {
107  mFileData.mPositions[i].mTime = scale * mFileData.mPositions[i].mTime + offset;
108  }
109 }
110 
117 void ReconstructPreprocessor::alignTimeSeries()
118 {
119  report("Generate time calibration based on input time stamps.");
120  double framesSpan = mFileData.mFrames.back().mTime - mFileData.mFrames.front().mTime;
121  double positionsSpan = mFileData.mPositions.back().mTime - mFileData.mPositions.front().mTime;
122  double scale = framesSpan / positionsSpan;
123 
124  double offset = mFileData.mFrames.front().mTime - scale * mFileData.mPositions.front().mTime;
125 
126  this->calibrateTimeStamps(offset, scale);
127 }
128 
132 void ReconstructPreprocessor::cropInputData()
133 {
134  //IntBoundingBox3D
135  ProbeDefinition sector = mFileData.mProbeDefinition.mData;
136  IntBoundingBox3D cropbox(sector.getClipRect_p().begin());
137  cropbox = this->reduceCropboxToImageSize(cropbox, sector.getSize());
138  Eigen::Vector3i shift = cropbox.corner(0,0,0).cast<int>();
139  Eigen::Vector3i size = cropbox.range().cast<int>() + Eigen::Vector3i(1,1,0); // convert from extent format to size format by adding 1
140  mFileData.mUsRaw->setCropBox(cropbox);
141 
142  DoubleBoundingBox3D clipRect_p = sector.getClipRect_p();
143  Vector3D origin_p = sector.getOrigin_p();
144 
145  for (unsigned i=0; i<3; ++i)
146  {
147  clipRect_p[2*i] -= shift[i];
148  clipRect_p[2*i+1] -= shift[i];
149  origin_p[i] -= shift[i];
150  }
151 
152  sector.setClipRect_p(clipRect_p);
153  sector.setOrigin_p(origin_p);
154  sector.setSize(QSize(size[0], size[1]));
155  mFileData.mProbeDefinition.setData(sector);
156 }
157 
158 IntBoundingBox3D ReconstructPreprocessor::reduceCropboxToImageSize(IntBoundingBox3D cropbox, QSize size)
159 {
160  cropbox[0] = std::max(cropbox[0], 0);
161  cropbox[2] = std::max(cropbox[2], 0);
162  cropbox[4] = std::max(cropbox[4], 0);
163 
164  cropbox[1] = std::min(cropbox[1], size.width() - 1);
165  cropbox[3] = std::min(cropbox[3], size.height() - 1);
166 
167  return cropbox;
168 }
169 
179 void ReconstructPreprocessor::applyTimeCalibration()
180 {
181  double timeshift = mInput.mExtraTimeCalibration;
182  // The shift is on frames. The calibrate function applies to tracker positions,
183  // hence the positive sign. (real use: subtract from frame data)
184  // std::cout << "TIMESHIFT " << timeshift << std::endl;
185  // timeshift = -timeshift;
186  if (!similar(0.0, timeshift))
187  report("Applying reconstruction-time calibration to tracking data: " + qstring_cast(
188  timeshift) + "ms");
189  this->calibrateTimeStamps(timeshift, 1.0);
190 
191  // ignore calibrations
192  if (mInput.mAlignTimestamps)
193  {
194  this->alignTimeSeries();
195  }
196 }
197 
198 
200 {
201  RemoveDataType() : count(0), err(-1) {}
202  void add(double _err) { ++count; err=((err<0)?_err:std::min(_err, err)); }
203  int count;
204  double err;
205 };
206 
213 void ReconstructPreprocessor::interpolatePositions()
214 {
215  mFileData.mUsRaw->resetRemovedFrames();
216  int startFrames = mFileData.mFrames.size();
217 
218  std::map<int,RemoveDataType> removedData;
219 
220  for (unsigned i_frame = 0; i_frame < mFileData.mFrames.size();)
221  {
222  std::vector<TimedPosition>::iterator posIter;
223  posIter = lower_bound(mFileData.mPositions.begin(), mFileData.mPositions.end(), mFileData.mFrames[i_frame]);
224 
225  unsigned i_pos = distance(mFileData.mPositions.begin(), posIter);
226  if (i_pos != 0)
227  i_pos--;
228 
229  if (i_pos >= mFileData.mPositions.size() - 1)
230  i_pos = mFileData.mPositions.size() - 2;
231 
232  // Remove frames too far from the positions
233  // Don't increment frame index since the index now points to the next element
234  double timeToPos1 = timeToPosition(i_frame, i_pos);
235  double timeToPos2 = timeToPosition(i_frame, i_pos+1);
236  if ((timeToPos1 > mMaxTimeDiff) || (timeToPos2 > mMaxTimeDiff))
237  {
238  removedData[i_frame].add(std::max(timeToPos1, timeToPos2));
239 
240  mFileData.mFrames.erase(mFileData.mFrames.begin() + i_frame);
241  mFileData.mUsRaw->removeFrame(i_frame);
242  }
243  else
244  {
245  double t_delta_tracking = mFileData.mPositions[i_pos + 1].mTime - mFileData.mPositions[i_pos].mTime;
246  double t = 0;
247  if (!similar(t_delta_tracking, 0))
248  t = (mFileData.mFrames[i_frame].mTime - mFileData.mPositions[i_pos].mTime) / t_delta_tracking;
249  // mFileData.mFrames[i_frame].mPos = interpolate(mFileData.mPositions[i_pos].mPos, mFileData.mPositions[i_pos + 1].mPos, t);
250  mFileData.mFrames[i_frame].mPos = cx::USReconstructInputDataAlgorithm::slerpInterpolate(mFileData.mPositions[i_pos].mPos, mFileData.mPositions[i_pos + 1].mPos, t);
251  i_frame++;// Only increment if we didn't delete the frame
252  }
253  }
254 
255  int removeCount=0;
256  for (std::map<int,RemoveDataType>::iterator iter=removedData.begin(); iter!=removedData.end(); ++iter)
257  {
258  int first = iter->first+removeCount;
259  int last = first + iter->second.count-1;
260  report(QString("Removed input frame [%1-%2]. Time diff=%3").arg(first).arg(last).arg(iter->second.err, 0, 'f', 1));
261  removeCount += iter->second.count;
262  }
263 
264  double removed = double(startFrames - mFileData.mFrames.size()) / double(startFrames);
265  if (removed > 0.02)
266  {
267  double percent = removed * 100;
268  if (percent > 1)
269  {
270  reportWarning("Removed " + QString::number(percent, 'f', 1) + "% of the "+ qstring_cast(startFrames) + " frames.");
271  }
272  else
273  {
274  report("Removed " + QString::number(percent, 'f', 1) + "% of the " + qstring_cast(startFrames) + " frames.");
275  }
276  }
277 }
278 
279 
280 double ReconstructPreprocessor::timeToPosition(unsigned i_frame, unsigned i_pos)
281 {
282  return fabs(mFileData.mFrames[i_frame].mTime - mFileData.mPositions[i_pos].mTime);
283 }
284 
288 std::vector<Vector3D> ReconstructPreprocessor::generateInputRectangle()
289 {
290  std::vector<Vector3D> retval(4);
291  vtkImageDataPtr mask = mFileData.getMask();
292  if (!mask)
293  {
294  reportError("Reconstructer::generateInputRectangle() + requires mask");
295  return retval;
296  }
297  Eigen::Array3i dims = mFileData.mUsRaw->getDimensions();
298  Vector3D spacing = mFileData.mUsRaw->getSpacing();
299 
300  Eigen::Array3i maskDims(mask->GetDimensions());
301 
302  if (( maskDims[0]<dims[0] )||( maskDims[1]<dims[1] ))
303  reportError(QString("input data (%1) and mask (%2) dim mimatch")
304  .arg(qstring_cast(dims))
305  .arg(qstring_cast(maskDims)));
306 
307  int xmin = maskDims[0];
308  int xmax = 0;
309  int ymin = maskDims[1];
310  int ymax = 0;
311 
312  unsigned char* ptr = static_cast<unsigned char*> (mask->GetScalarPointer());
313  for (int x = 0; x < maskDims[0]; x++)
314  for (int y = 0; y < maskDims[1]; y++)
315  {
316  unsigned char val = ptr[x + y * maskDims[0]];
317  if (val != 0)
318  {
319  xmin = std::min(xmin, x);
320  ymin = std::min(ymin, y);
321  xmax = std::max(xmax, x);
322  ymax = std::max(ymax, y);
323  }
324  }
325 
326  //Test: reduce the output volume by reducing the mask when determining
327  // output volume size
328  double red = mInput.mMaskReduce;
329  int reduceX = (xmax - xmin) * (red / 100);
330  int reduceY = (ymax - ymin) * (red / 100);
331 
332  xmin += reduceX;
333  xmax -= reduceX;
334  ymin += reduceY;
335  ymax -= reduceY;
336 
337  retval[0] = Vector3D(xmin * spacing[0], ymin * spacing[1], 0);
338  retval[1] = Vector3D(xmax * spacing[0], ymin * spacing[1], 0);
339  retval[2] = Vector3D(xmin * spacing[0], ymax * spacing[1], 0);
340  retval[3] = Vector3D(xmax * spacing[0], ymax * spacing[1], 0);
341 
342  return retval;
343 }
344 
351 Transform3D ReconstructPreprocessor::applyOutputOrientation()
352 {
353  Transform3D prMdd = Transform3D::Identity();
354 
355  if (mInput.mOrientation == "PatientReference")
356  {
357  // identity
358  }
359  else if (mInput.mOrientation == "MiddleFrame")
360  {
361  prMdd = mFileData.mFrames[mFileData.mFrames.size() / 2].mPos;
362  }
363  else
364  {
365  reportError("no orientation algorithm selected in reconstruction");
366  }
367 
368  // apply the selected orientation to the frames.
369  Transform3D ddMpr = prMdd.inv();
370  for (unsigned i = 0; i < mFileData.mFrames.size(); i++)
371  {
372  // mPos = prMu
373  mFileData.mFrames[i].mPos = ddMpr * mFileData.mFrames[i].mPos;
374  // mPos = ddMu
375  }
376 
377  return prMdd;
378 }
379 
389 void ReconstructPreprocessor::findExtentAndOutputTransform()
390 {
391  if (mFileData.mFrames.empty())
392  return;
393  // A first guess for d'Mu with correct orientation
394  Transform3D prMdd = this->applyOutputOrientation();
395  //mFrames[i].mPos = d'Mu, d' = only rotation
396 
397  // Find extent of all frames as a point cloud
398  std::vector<Vector3D> inputRect = this->generateInputRectangle();
399  std::vector<Vector3D> outputRect;
400  for (unsigned slice = 0; slice < mFileData.mFrames.size(); slice++)
401  {
402  Transform3D dMu = mFileData.mFrames[slice].mPos;
403  for (unsigned i = 0; i < inputRect.size(); i++)
404  {
405  outputRect.push_back(dMu.coord(inputRect[i]));
406  }
407  }
408 
409  DoubleBoundingBox3D extent = DoubleBoundingBox3D::fromCloud(outputRect);
410 
411  // Translate dMu to output volume origo
412  Transform3D T_origo = createTransformTranslate(extent.corner(0, 0, 0));
413  Transform3D prMd = prMdd * T_origo; // transform from output space to patref, use when storing volume.
414  for (unsigned i = 0; i < mFileData.mFrames.size(); i++)
415  {
416  mFileData.mFrames[i].mPos = T_origo.inv() * mFileData.mFrames[i].mPos;
417  }
418 
419  // Calculate optimal output image spacing and dimensions based on US frame spacing
420  double inputSpacing = std::min(mFileData.mUsRaw->getSpacing()[0], mFileData.mUsRaw->getSpacing()[1]);
421  mOutputVolumeParams = OutputVolumeParams(extent, inputSpacing, mInput.mMaxOutputVolumeSize);
422 
423  mOutputVolumeParams.set_rMd((mPatientModelService->get_rMpr()) * prMd);
424 }
425 
429 void ReconstructPreprocessor::updateFromOriginalFileData()
430 {
431  // uncomment to test cropping of data before reconstructing
432  this->cropInputData();
433 
434  // Only use this if the time stamps have different formats
435  // The function assumes that both lists of time stamps start at the same time,
436  // and that is not completely correct.
437  //this->calibrateTimeStamps();
438  // Use the time calibration from the acquisition module
439  //this->calibrateTimeStamps(0.0, 1.0);
440  this->applyTimeCalibration();
441 
443  //mPos (in mPositions) is now prMu
444 
445  this->interpolatePositions();
446  // mFrames: now mPos as prMu
447 
448  if (mFileData.mFrames.empty()) // if all positions are filtered out
449  return;
450 
451  this->findExtentAndOutputTransform();
452  //mPos in mFrames is now dMu
453 }
454 
456 {
457  mInput = input;
458  mFileData = fileData;
459  this->updateFromOriginalFileData();
460 }
461 
462 
463 } /* namespace cx */
QString qstring_cast(const T &val)
void reportError(QString msg)
Definition: cxLogger.cpp:92
#define CX_ASSERT(statement)
Definition: cxLogger.h:131
Transform3D Transform3D
Transform3D is a representation of an affine 3D transform.
std::vector< TimedPosition > mFrames
static void transformTrackingPositionsTo_prMu(USReconstructInputData *data)
static DoubleBoundingBox3D fromCloud(std::vector< Vector3D > cloud)
bool similar(const DoubleBoundingBox3D &a, const DoubleBoundingBox3D &b, double tol)
ReconstructPreprocessor(PatientModelServicePtr patientModelService)
void reportWarning(QString msg)
Definition: cxLogger.cpp:91
boost::shared_ptr< class PatientModelService > PatientModelServicePtr
Transform3D createTransformTranslate(const Vector3D &translation)
std::vector< ProcessedUSInputDataPtr > createProcessedInput(std::vector< bool > angio)
std::vector< TimedPosition > mPositions
ProbeDefinition mData
Definition: cxProbeSector.h:75
Eigen::Vector3d Vector3D
Vector3D is a representation of a point or vector in 3D.
Definition: cxVector3D.h:63
Eigen::Vector3i corner(int x, int y, int z) const
void report(QString msg)
Definition: cxLogger.cpp:90
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)
Algorithm part of reconstruction - no dependencies on parameter classes.
void setData(ProbeDefinition data)