17 #include <vtkImageData.h>
39 mPatientModelService(patientModelService)
54 std::vector<std::vector<vtkImageDataPtr> > frames = mFileData.
mUsRaw->initializeFrames(angio);
56 std::vector<ProcessedUSInputDataPtr> retval;
58 for (
unsigned i=0; i<angio.size(); ++i)
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);
81 void ReconstructPreprocessor::calibrateTimeStamps(
double offset,
double scale)
83 for (
unsigned int i = 0; i < mFileData.
mPositions.size(); i++)
95 void ReconstructPreprocessor::alignTimeSeries()
97 report(
"Generate time calibration based on input time stamps.");
98 double framesSpan = mFileData.
mFrames.back().mTime - mFileData.
mFrames.front().mTime;
100 double scale = framesSpan / positionsSpan;
102 double offset = mFileData.
mFrames.front().mTime - scale * mFileData.
mPositions.front().mTime;
104 this->calibrateTimeStamps(offset, scale);
110 void ReconstructPreprocessor::cropInputData()
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);
118 mFileData.
mUsRaw->setCropBox(cropbox);
120 DoubleBoundingBox3D clipRect_p = sector.getClipRect_p();
121 Vector3D origin_p = sector.getOrigin_p();
123 for (
unsigned i=0; i<3; ++i)
125 clipRect_p[2*i] -= shift[i];
126 clipRect_p[2*i+1] -= shift[i];
127 origin_p[i] -= shift[i];
130 sector.setClipRect_p(clipRect_p);
131 sector.setOrigin_p(origin_p);
132 sector.setSize(QSize(size[0], size[1]));
136 IntBoundingBox3D ReconstructPreprocessor::reduceCropboxToImageSize(IntBoundingBox3D cropbox, QSize size)
138 cropbox[0] = std::max(cropbox[0], 0);
139 cropbox[2] = std::max(cropbox[2], 0);
140 cropbox[4] = std::max(cropbox[4], 0);
142 cropbox[1] = std::min(cropbox[1], size.width() - 1);
143 cropbox[3] = std::min(cropbox[3], size.height() - 1);
157 void ReconstructPreprocessor::applyTimeCalibration()
167 this->calibrateTimeStamps(timeshift, 1.0);
172 this->alignTimeSeries();
185 void ReconstructPreprocessor::filterPositions()
188 CX_LOG_DEBUG() <<
"Filter length specified " << filterStrength;
194 void ReconstructPreprocessor::positionThinning()
207 void ReconstructPreprocessor::interpolatePositions()
209 mFileData.
mUsRaw->resetRemovedFrames();
210 unsigned long startFrames = mFileData.
mFrames.size();
212 std::map<unsigned,RemoveDataType> removedData;
214 for (
unsigned i_frame = 0; i_frame < mFileData.
mFrames.size();)
216 std::vector<TimedPosition>::iterator posIter;
219 unsigned i_pos = unsigned(distance(mFileData.
mPositions.begin(), posIter));
224 i_pos = unsigned(mFileData.
mPositions.size()) - 2;
228 double timeToPos1 = timeToPosition(i_frame,
unsigned(i_pos));
229 double timeToPos2 = timeToPosition(i_frame, i_pos+1);
230 if ((timeToPos1 > mMaxTimeDiff) || (timeToPos2 > mMaxTimeDiff))
232 removedData[i_frame].add(std::max(timeToPos1, timeToPos2));
235 mFileData.
mUsRaw->removeFrame(i_frame);
239 double t_delta_tracking = mFileData.
mPositions[i_pos + 1].mTime - mFileData.
mPositions[i_pos].mTime;
241 if (!
similar(t_delta_tracking, 0))
242 t = (mFileData.
mFrames[i_frame].mTime - mFileData.
mPositions[i_pos].mTime) / t_delta_tracking;
249 unsigned removeCount=0;
250 for (std::map<unsigned,RemoveDataType>::iterator iter=removedData.begin(); iter!=removedData.end(); ++iter)
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;
258 double removed = double(startFrames - mFileData.
mFrames.size()) / double(startFrames);
261 double percent = removed * 100;
268 report(
"Removed " + QString::number(percent,
'f', 1) +
"% of the " +
qstring_cast(startFrames) +
" frames.");
274 double ReconstructPreprocessor::timeToPosition(
unsigned i_frame,
unsigned i_pos)
276 return fabs(mFileData.
mFrames[i_frame].mTime - mFileData.
mPositions[i_pos].mTime);
282 std::vector<Vector3D> ReconstructPreprocessor::generateInputRectangle()
284 std::vector<Vector3D> retval(4);
288 reportError(
"Reconstructer::generateInputRectangle() + requires mask");
291 Eigen::Array3i dims = mFileData.
mUsRaw->getDimensions();
294 Eigen::Array3i maskDims(mask->GetDimensions());
296 if (( maskDims[0]<dims[0] )||( maskDims[1]<dims[1] ))
297 reportError(QString(
"input data (%1) and mask (%2) dim mimatch")
301 int xmin = maskDims[0];
303 int ymin = maskDims[1];
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++)
310 unsigned char val = ptr[x + y * maskDims[0]];
313 xmin = std::min(xmin, x);
314 ymin = std::min(ymin, y);
315 xmax = std::max(xmax, x);
316 ymax = std::max(ymax, y);
322 int reduceX = int((xmax - xmin) * (red / 100));
323 int reduceY = int((ymax - ymin) * (red / 100));
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);
344 Transform3D ReconstructPreprocessor::applyOutputOrientation()
358 reportError(
"no orientation algorithm selected in reconstruction");
363 for (
unsigned i = 0; i < mFileData.
mFrames.size(); i++)
382 void ReconstructPreprocessor::findExtentAndOutputTransform()
387 Transform3D prMdd = this->applyOutputOrientation();
391 std::vector<Vector3D> inputRect = this->generateInputRectangle();
392 std::vector<Vector3D> outputRect;
393 for (
unsigned slice = 0; slice < mFileData.
mFrames.size(); slice++)
396 for (
unsigned i = 0; i < inputRect.size(); i++)
398 outputRect.push_back(dMu.coord(inputRect[i]));
407 for (
unsigned i = 0; i < mFileData.
mFrames.size(); i++)
409 mFileData.
mFrames[i].mPos = T_origo.inv() * mFileData.
mFrames[i].mPos;
413 double inputSpacing = std::min(mFileData.
mUsRaw->getSpacing()[0], mFileData.
mUsRaw->getSpacing()[1]);
416 mOutputVolumeParams.
set_rMd((mPatientModelService->get_rMpr()) * prMd);
422 void ReconstructPreprocessor::updateFromOriginalFileData()
425 this->cropInputData();
433 this->applyTimeCalibration();
438 this->filterPositions();
443 this->interpolatePositions();
449 this->findExtentAndOutputTransform();
456 mFileData = fileData;
457 this->updateFromOriginalFileData();