Fraxinus  18.10
An IGT application
cxUSFrameData.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 
13 #include "cxUSFrameData.h"
14 
15 #include <vtkImageData.h>
16 #include <vtkImageLuminance.h>
17 #include <vtkImageClip.h>
18 #include <vtkImageAppend.h>
19 #include <vtkMetaImageWriter.h>
20 #include <vtkImageImport.h>
21 #include "cxTypeConversions.h"
22 #include "cxDataReaderWriter.h"
23 #include <QFileInfo>
24 #include "cxTimeKeeper.h"
25 #include "cxImageDataContainer.h"
26 #include "cxVolumeHelpers.h"
27 #include "cxLogger.h"
28 
29 
30 typedef vtkSmartPointer<vtkImageAppend> vtkImageAppendPtr;
31 
32 namespace cx
33 {
34 
35 ProcessedUSInputData::ProcessedUSInputData(std::vector<vtkImageDataPtr> frames, std::vector<TimedPosition> pos, vtkImageDataPtr mask, QString path, QString uid) :
36  mProcessedImage(frames),
37  mFrames(pos),
38  mMask(mask),
39  mPath(path),
40  mUid(uid)
41 {
42  this->validate();
43 }
44 
46 {
47  std::vector<TimedPosition> frameInfo = this->getFrames();
48  Eigen::Array3i inputDims = this->getDimensions();
49 
50  if (inputDims[2] != static_cast<int>(frameInfo.size()))
51  {
52  QString msg("Mismatch in US input data: inputDims[2] != frameInfo.size() : %1 != %2");
53  reportWarning(msg.arg(inputDims[2]).arg(frameInfo.size()));
54  return false;
55  }
56 
57  if (inputDims[2]==0)
58  {
59  reportWarning("Empty US input data");
60  return false;
61  }
62 
63  return true;
64 }
65 
66 unsigned char* ProcessedUSInputData::getFrame(unsigned int index) const
67 {
68  CX_ASSERT(index < mProcessedImage.size());
69 
70  // Raw data pointer
71  unsigned char *inputPointer = static_cast<unsigned char*> (mProcessedImage[index]->GetScalarPointer());
72  return inputPointer;
73 }
74 
76 {
77  Eigen::Array3i retval;
78  retval[0] = mProcessedImage[0]->GetDimensions()[0];
79  retval[1] = mProcessedImage[0]->GetDimensions()[1];
80  retval[2] = mProcessedImage.size();
81  return retval;
82 }
83 
85 {
86  Vector3D retval = Vector3D(mProcessedImage[0]->GetSpacing());
87  retval[2] = retval[0]; // set z-spacing to arbitrary value.
88  return retval;
89 }
90 
91 std::vector<TimedPosition> ProcessedUSInputData::getFrames() const
92 {
93  return mFrames;
94 }
95 
97 {
98  return mMask;
99 }
100 
102 {
103  return mPath;
104 }
105 
107 {
108  return mUid;
109 }
110 
111 
115 
117 {
118  USFrameDataPtr retval(new USFrameData());
119 
120  retval->mName = QFileInfo(inputFrameData->getName()).completeBaseName();
121  retval->mImageContainer.reset(new cx::SplitFramesContainer(inputFrameData->getBaseVtkImageData()));
122  retval->resetRemovedFrames();
123 
124  return retval;
125 }
126 
133 USFrameDataPtr USFrameData::create(QString inputFilename)
134 {
135  QFileInfo info(inputFilename);
136 
137  TimeKeeper timer;
138  QString mhdSingleFile = info.absolutePath()+"/"+info.completeBaseName()+".mhd";
139 
140  if (QFileInfo(mhdSingleFile).exists())
141  {
142  vtkImageDataPtr image = MetaImageReader().loadVtkImageData(mhdSingleFile);
143  // load from single file
144  USFrameDataPtr retval = USFrameData::create(ImagePtr(new Image(mhdSingleFile, image)));
145  retval->mName = QFileInfo(mhdSingleFile).completeBaseName();
146  timer.printElapsedms(QString("Loading single %1").arg(inputFilename));
147  return retval;
148  }
149  else
150  {
151  USFrameDataPtr retval(new USFrameData());
152  retval->mName = QFileInfo(inputFilename).completeBaseName();
153  retval->mImageContainer.reset(new cx::CachedImageDataContainer(inputFilename, -1));
154  retval->resetRemovedFrames();
155  return retval;
156  }
157 }
158 
160 {
161  if (!images)
162  return USFrameDataPtr();
163  USFrameDataPtr retval(new USFrameData());
164  retval->mName = name;
165  retval->mImageContainer = images;
166  retval->resetRemovedFrames();
167 
168  return retval;
169 }
170 
171 USFrameDataPtr USFrameData::create(QString name, std::vector<vtkImageDataPtr> frames)
172 {
173  cx::ImageDataContainerPtr container;
174  container.reset(new cx::FramesDataContainer(frames));
175  return create(name, container);
176 }
177 
178 
180  mCropbox(0,0,0,0,0,0), mPurgeInput(true)
181 {
182 }
183 
185 {
186  if (mImageContainer->empty())
187  return;
188 
189  mReducedToFull.clear();
190  for (unsigned i=0; i<mImageContainer->size(); ++i)
191  mReducedToFull.push_back(i);
192 }
193 
195 {
196  mImageContainer.reset();
197 }
198 
202 void USFrameData::removeFrame(unsigned int index)
203 {
204  CX_ASSERT(index < mReducedToFull.size());
205  mReducedToFull.erase(mReducedToFull.begin()+index);
206 }
207 
208 Eigen::Array3i USFrameData::getDimensions() const
209 {
210  Eigen::Array3i retval(mImageContainer->get(0)->GetDimensions());
211 
212  if (mCropbox.range()[0]!=0)
213  {
214  // WARNING: cannot use dimensions from the cropImage function - it uses extent,
215  // which may be larger. We need the real size, given by wholeextent/updateextent.
216  vtkImageDataPtr sample = this->cropImageExtent(mImageContainer->get(0), mCropbox);
217 // Eigen::Array<int, 6, 1> extent(sample->GetWholeExtent());
218  Eigen::Array<int, 6, 1> extent(sample->GetExtent());
219  retval[0] = extent[1]-extent[0]+1;
220  retval[1] = extent[3]-extent[2]+1;
221 // std::cout << "extent dim: " << retval[0] << " " << retval[1] << std::endl;
222  }
223 
224  retval[2] = mReducedToFull.size();
225  return retval;
226 }
227 
228 
230 {
231  Vector3D retval(1,1,1);
232  if (!mImageContainer->empty())
233  retval = Vector3D(mImageContainer->get(0)->GetSpacing());
234  retval[2] = retval[0]; // set z-spacing to arbitrary value.
235  return retval;
236 }
237 
239 {
240  // ensure clip never happens in z dir.
241  cropbox[4] = -100000;
242  cropbox[5] = 100000;
243  mCropbox = cropbox;
244 }
245 
253 {
254  vtkImageClipPtr clip = vtkImageClipPtr::New();
255  clip->SetInputData(input);
256  clip->SetOutputWholeExtent(cropbox.begin());
257 // clip->ClipDataOn(); // this line sets wholeextent==extent, but uses lots of time (~6ms/frame)
258  clip->Update();
259  vtkImageDataPtr rawResult = clip->GetOutput();
260 // rawResult->Update();
261  rawResult->Crop(cropbox.begin()); // moved in from toGrayscaleAndEffectuateCropping when VTK5->6
262  return rawResult;
263 }
264 
270 {
271  vtkImageDataPtr grayScaleData;
272 
273  if (input->GetNumberOfScalarComponents() == 1) // already gray
274  {
275  // crop (slow)
276  grayScaleData = input;
277 // outData->Crop();
278  }
279  else
280  {
281  // convert and crop as side effect (optimization)
282  grayScaleData = convertImageDataToGrayScale(input);
283  }
284 
285  vtkImageDataPtr outData = this->convertTo8bit(grayScaleData);
286 
287  vtkImageDataPtr copy = vtkImageDataPtr::New();
288  copy->DeepCopy(outData);
289  return copy;
290 }
291 
292 vtkImageDataPtr USFrameData::convertTo8bit(vtkImageDataPtr input) const
293 {
294  vtkImageDataPtr retval = input;
295  if (input->GetScalarSize() > 1)
296  {
297  ImagePtr tempImage = cx::ImagePtr(new cx::Image("tempImage", input, "tempImage"));
298  tempImage->resetTransferFunctions();
299  retval = tempImage->get8bitGrayScaleVtkImageData();
300  }
301  return retval;
302 }
303 
304 // for testing
305 //void printStuff(QString text, vtkImageDataPtr data)
306 //{
307 // std::cout << text << std::endl;
308 // Eigen::Array<int, 6, 1> wholeExtent(data->GetWholeExtent());
309 // Eigen::Array<int, 6, 1> extent(data->GetExtent());
310 // Eigen::Array<int, 6, 1> updateExtent(data->GetUpdateExtent());
311 // std::cout << "\t whole ext " << wholeExtent << std::endl;
312 // std::cout << "\tupdate ext " << updateExtent << std::endl;
313 // std::cout << "\t ext " << extent << std::endl;
315 
316 // Eigen::Array<vtkIdType,3,1> hinc;
317 // data->GetContinuousIncrements(data->GetWholeExtent(), hinc[0], hinc[1], hinc[2]);
318 // std::cout << "\t whole inc " << hinc << std::endl;
319 
320 // Eigen::Array<vtkIdType,3,1> uinc;
321 // data->GetContinuousIncrements(data->GetUpdateExtent(), uinc[0], uinc[1], uinc[2]);
322 // std::cout << "\t uinc " << uinc << std::endl;
323 
324 // Eigen::Array<vtkIdType,3,1> inc;
325 // data->GetContinuousIncrements(data->GetExtent(), inc[0], inc[1], inc[2]);
326 // std::cout << "\t inc " << inc << std::endl;
327 
328 // Eigen::Array3i dim(data->GetDimensions());
329 // std::cout << "\t dim " << dim << std::endl;
330 
331 // Eigen::Array3i max(wholeExtent[1] - wholeExtent[0], wholeExtent[3] - wholeExtent[2], wholeExtent[5] - wholeExtent[4]);
332 // std::cout << "\t max " << max << std::endl;
333 //}
334 
336 {
337  // Some of the code here is borrowed from the vtk examples:
338  // http://public.kitware.com/cgi-bin/viewcvs.cgi/*checkout*/Examples/Build/vtkMy/Imaging/vtkImageFoo.cxx?root=VTK&content-type=text/plain
339 
340  if (inData->GetNumberOfScalarComponents() != 3)
341  {
342  if(frameNum == 0) //Only report warning once
343  reportWarning("Angio requested for grayscale ultrasound");
344  return grayFrame;
345  }
346 
347  vtkImageDataPtr outData = vtkImageDataPtr::New();
348  outData->DeepCopy(grayFrame);
349 // outData->Update(); // updates whole extent.
350 
351 // printStuff("Clipped color in", inData);
352 // printStuff("Grayscale in", outData);
353 
354 // int* inExt = inData->GetWholeExtent();
355  int* outExt = outData->GetExtent();
356 
357  // Remember that the input might (and do due to vtkImageClip) contain leaps.
358  // This means that the wholeextent might be larger than the extent, thus
359  // we must use a startoffset and leaps between lines.
360 
361  unsigned char *inPtr = static_cast<unsigned char*> (inData->GetScalarPointerForExtent(inData->GetExtent()));
362  unsigned char *outPtr = static_cast<unsigned char*> (outData->GetScalarPointerForExtent(outData->GetExtent()));
363 
364  int maxX = outExt[1] - outExt[0];
365  int maxY = outExt[3] - outExt[2];
366  int maxZ = outExt[5] - outExt[4];
367 
368  Eigen::Array<vtkIdType,3,1> inInc;
369  inData->GetContinuousIncrements(inData->GetExtent(), inInc[0], inInc[1], inInc[2]);
370  CX_ASSERT(inInc[0]==0);
371  // we assume (wholeextent == extent) for the outData in the algorithm below. assert here.
372  Eigen::Array<vtkIdType,3,1> outInc;
373  outData->GetContinuousIncrements(outData->GetExtent(), outInc[0], outInc[1], outInc[2]);
374  CX_ASSERT(outInc[0]==0);
375  CX_ASSERT(outInc[1]==0);
376  CX_ASSERT(outInc[2]==0);
377 
378  for (int z=0; z<=maxZ; z++)
379  {
380  for (int y=0; y<=maxY; y++)
381  {
382  for (int x=0; x<=maxX; x++)
383  {
384  //Look at 3 scalar components at the same time (RGB),
385  //if color is grayscale or close to grayscale: set to zero.
386 
387  if (((*inPtr) == (*(inPtr + 1))) && ((*inPtr) == (*(inPtr + 2))))
388  {
389  (*outPtr) = 0;
390  }
391  else
392  {
393  // strong check: look for near-gray values and remove them.
394  double r = inPtr[0];
395  double g = inPtr[1];
396  double b = inPtr[2];
397  int metric = (fabs(r-g) + fabs(r-b) + fabs(g-b)) / 3; // average absolute diff must be less than or equal to this
398 // std::cout << QString(" %1,%2,%3 \t %4, %5 -- %6").arg(int(inPtr[0])).arg(int(inPtr[1])).arg(int(inPtr[2])).arg(idxR).arg(idxY).arg(metric) << std::endl;
399  if (metric <= 3)
400  (*outPtr) = 0;
401 
402  }
403  //Assume the outVolume is treated with the luminance filter first
404  outPtr++;
405  inPtr += 3;
406  }
407  inPtr += inInc[1];
408  }
409  inPtr += inInc[2];
410  }
411 
412  return outData;
413 }
414 
416 {
417  mPurgeInput = value;
418 }
419 
420 std::vector<std::vector<vtkImageDataPtr> > USFrameData::initializeFrames(std::vector<bool> angio)
421 {
422 
423  std::vector<std::vector<vtkImageDataPtr> > raw(angio.size());
424 
425  for (unsigned i=0; i<raw.size(); ++i)
426  {
427  raw[i].resize(mReducedToFull.size());
428  }
429 
430  // apply cropping and angio
431  for (unsigned i=0; i<mReducedToFull.size(); ++i)
432  {
435 
436  if (mCropbox.range()[0]!=0)
437  current = this->cropImageExtent(current, mCropbox);
438 
439  // optimization: grayFrame is used in both calculations: compute once
440  vtkImageDataPtr grayFrame = this->to8bitGrayscaleAndEffectuateCropping(current);
441 
442  for (unsigned j=0; j<angio.size(); ++j)
443  {
444 
445  if (angio[j])
446  {
447  vtkImageDataPtr angioFrame = this->useAngio(current, grayFrame, i);
448  raw[j][i] = angioFrame;
449  }
450  else
451  {
452  raw[j][i] = grayFrame;
453  }
454  }
455 
456  if (mPurgeInput)
457  mImageContainer->purge(mReducedToFull[i]);
458  }
459 
460  if (mPurgeInput)
461  mImageContainer->purgeAll();
462 
463  return raw;
464 }
465 
467 {
468  mImageContainer->purgeAll();
469 }
470 
472 {
473  USFrameDataPtr retval(new USFrameData(*this));
474  this->resetRemovedFrames();
475  return retval;
476 }
477 
478 QString USFrameData::getName() const
479 {
480  return QFileInfo(mName).completeBaseName();
481 }
482 
484 {
485  return mImageContainer->size();
486 }
487 
489 {
490  TimeKeeper timer;
491  vtkImageDataPtr source = mImageContainer->get(index);
492 
493 // use this test code to display angio output:
494 // vtkImageDataPtr current = mImageContainer->get(index);
495 // current = this->cropImage(current, IntBoundingBox3D(157, 747, 68, 680, 0, 0));
496 // vtkImageDataPtr grayFrame = this->toGrayscale(current);
497 // static vtkImageDataPtr source;
498 // source = this->useAngio(current, grayFrame);
499 
500 // source->Update(); // VTK5->6
501 
502  import->SetImportVoidPointer(source->GetScalarPointer());
503  import->SetDataScalarType(source->GetScalarType());
504  import->SetDataSpacing(source->GetSpacing());
505  import->SetNumberOfScalarComponents(source->GetNumberOfScalarComponents());
506  import->SetWholeExtent(source->GetExtent());
507  import->SetDataExtentToWholeExtent();
508 }
509 
511 {
512  int numberOfFrames = mReducedToFull.size();
513 
514  if (numberOfFrames > 0)
515  if(mImageContainer->get(0)->GetDataDimension() == 3)
516  return true;
517  return false;
518 }
519 
521 {
522  if (mImageContainer->get(0)->GetNumberOfScalarComponents() == 1)
523  return true;
524  return false;
525 }
526 
527 }//namespace cx
vtkImageDataPtr cropImageExtent(vtkImageDataPtr input, IntBoundingBox3D cropbox) const
Use only US angio data as input. Removes grayscale from the US data and converts the remaining color ...
bool is8bit() const
unsigned char * getFrame(unsigned int index) const
vtkSmartPointer< vtkImageAppend > vtkImageAppendPtr
std::vector< std::vector< vtkImageDataPtr > > initializeFrames(std::vector< bool > angio)
cx::ImageDataContainerPtr mImageContainer
#define CX_ASSERT(statement)
Definition: cxLogger.h:116
void setPurgeInputDataAfterInitialize(bool value)
boost::shared_ptr< class Image > ImagePtr
Definition: cxDicomWidget.h:27
boost::shared_ptr< class USFrameData > USFrameDataPtr
void printElapsedms(QString text="Elapsed") const
std::vector< int > mReducedToFull
map from indexes in the reduced volume to the full (original) volume.
virtual vtkImageDataPtr loadVtkImageData(QString filename)
Eigen::Array3i getDimensions() const
IntBoundingBox3D mCropbox
unsigned getNumImages()
QString getName() const
std::vector< TimedPosition > getFrames() const
Helper class encapsulating a 2S US data set.
Definition: cxUSFrameData.h:87
virtual USFrameDataPtr copy()
Eigen::Vector3i range() const
void reportWarning(QString msg)
Definition: cxLogger.cpp:70
vtkImageDataPtr getMask()
Eigen::Array3i getDimensions() const
A volumetric data set.
Definition: cxImage.h:45
vtkImageDataPtr useAngio(vtkImageDataPtr inData, vtkImageDataPtr grayFrame, int frameNum) const
void fillImageImport(vtkImageImportPtr import, int index)
fill import with a single frame
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.
vtkSmartPointer< class vtkImageClip > vtkImageClipPtr
Eigen::Vector3d Vector3D
Vector3D is a representation of a point or vector in 3D.
Definition: cxVector3D.h:42
vtkImageDataPtr convertImageDataToGrayScale(vtkImageDataPtr image)
vtkSmartPointer< class vtkImageImport > vtkImageImportPtr
Vector3D getSpacing() const
void setCropBox(IntBoundingBox3D mCropbox)
Vector3D getSpacing() const
vtkSmartPointer< class vtkImageData > vtkImageDataPtr
Reader for metaheader .mhd files.
vtkImageDataPtr to8bitGrayscaleAndEffectuateCropping(vtkImageDataPtr input) const
ProcessedUSInputData(std::vector< vtkImageDataPtr > frames, std::vector< TimedPosition > pos, vtkImageDataPtr mask, QString path, QString uid)
static USFrameDataPtr create(ImagePtr inputFrameData)
void removeFrame(unsigned int index)
boost::shared_ptr< class ImageDataContainer > ImageDataContainerPtr
Definition: cxUSFrameData.h:24
Namespace for all CustusX production code.