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