Fraxinus  17.12-rc2
An IGT application
cxVolumeHelpers.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 #include "cxVolumeHelpers.h"
34 
35 #include <vtkUnsignedCharArray.h>
36 #include <vtkImageData.h>
37 #include <vtkPointData.h>
38 #include <vtkDoubleArray.h>
39 #include <vtkImageResample.h>
40 #include <vtkImageClip.h>
41 #include <vtkImageShiftScale.h>
42 #include <vtkImageAccumulate.h>
43 #include <vtkImageLuminance.h>
44 #include <vtkImageExtractComponents.h>
45 #include <vtkImageAppendComponents.h>
46 
47 #include "cxImage.h"
48 
49 #include "cxUtilHelpers.h"
50 #include "cxImageTF3D.h"
51 #include "cxImageLUT2D.h"
53 #include "cxLogger.h"
54 #include "cxEnumConverter.h"
55 #include "cxTime.h"
57 #include "cxPatientModelService.h"
58 
59 typedef vtkSmartPointer<vtkDoubleArray> vtkDoubleArrayPtr;
60 
61 namespace cx
62 {
63 
64 namespace {
65 template<class TYPE, int TYPE_FROM_VTK>
66 vtkImageDataPtr generateVtkImageDataGeneric(Eigen::Array3i dim,
67  Vector3D spacing,
68  const TYPE initValue,
69  int components)
70 {
71  vtkImageDataPtr data = vtkImageDataPtr::New();
72  data->SetSpacing(spacing[0], spacing[1], spacing[2]);
73  data->SetExtent(0, dim[0]-1, 0, dim[1]-1, 0, dim[2]-1);
74 // data->SetScalarType(TYPE_FROM_VTK);
75 // data->SetNumberOfScalarComponents(components);
76 // data->AllocateScalars();
77  data->AllocateScalars(TYPE_FROM_VTK, components);
78 
79  int scalarSize = dim[0]*dim[1]*dim[2]*components;
80 
81  TYPE* ptr = reinterpret_cast<TYPE*>(data->GetScalarPointer());
82  std::fill(ptr, ptr+scalarSize, initValue);
83 
84  // FIXED: replace with setDeepModified(image)
85  // A trick to get a full LUT in Image (automatic LUT generation)
86  // Can't seem to fix this by calling Image::resetTransferFunctions() after volume is modified
87 /* if (scalarSize > 0)
88  {
89  ptr[0] = 150;
90 // ptr[0] = 255;
91  if (scalarSize > 1)
92  ptr[1] = 50;
93  data->GetScalarRange();// Update internal data in vtkImageData. Seems like it is not possible to update this data after the volume has been changed.
94  ptr[0] = initValue;
95  if (scalarSize > 1)
96  ptr[1] = initValue;
97  }*/
98 // data->UpdateInformation(); // update extent etc
99  setDeepModified(data);
100 
101  return data;
102 }
103 } // namespace
104 
106  Vector3D spacing,
107  const unsigned char initValue,
108  int components)
109 {
110  return generateVtkImageDataGeneric<unsigned char, VTK_UNSIGNED_CHAR>(dim, spacing, initValue, components);
111 }
112 
114  Vector3D spacing,
115  const unsigned short initValue,
116  int components)
117 {
118  return generateVtkImageDataGeneric<unsigned short, VTK_UNSIGNED_SHORT>(dim, spacing, initValue, components);
119 }
120 
122  Vector3D spacing,
123  const short initValue,
124  int components)
125 {
126  return generateVtkImageDataGeneric<short, VTK_SHORT>(dim, spacing, initValue, components);
127 }
128 
130  Vector3D spacing,
131  double initValue)
132 {
133  return generateVtkImageDataGeneric<double, VTK_DOUBLE>(dim, spacing, initValue, 1);
134 }
135 
137 {
138  Eigen::Array3i dim(data->GetDimensions());
139 
140  unsigned short* ptr = reinterpret_cast<unsigned short*>(data->GetScalarPointer());
141 
142 // int scalarSize = dim[0]*dim[1]*dim[2]*1;
143  for (int z=0; z<dim[2]; ++z)
144  for (int y=0; y<dim[1]; ++y)
145  for (int x=0; x<dim[0]; ++x)
146  {
147  int mod = maxValue;
148  int val = int((double(z)/dim[2]*mod*6))%mod;// + y%255 + x/255;
149  if (val < mod/3)
150  val = 0;
151  ptr[z*dim[0]*dim[1] + y*dim[0] + x] = val;
152  }
153 
154  // unsigned char* ptr = reinterpret_cast<unsigned char*>(imageData->GetScalarPointer());
155 
156  // int scalarSize = dim[0]*dim[1]*dim[2]*1;
157  // for (int z=0; z<dim[2]; ++z)
158  // for (int y=0; y<dim[1]; ++y)
159  // for (int x=0; x<dim[0]; ++x)
160  // {
161  // int val = int((double(z)/dim[2]*255*6))%255;// + y%255 + x/255;
162  // val = val/3;
163  // ptr[z*dim[0]*dim[1] + y*dim[0] + x] = val;
164  // }
165  setDeepModified(data);
166 }
167 
172 ImagePtr createDerivedImage(PatientModelServicePtr dataManager, QString uid, QString name, vtkImageDataPtr raw, ImagePtr parent)
173 {
174  ImagePtr retval = dataManager->createSpecificData<Image>(uid, name);
175  retval->setVtkImageData(raw);
176  retval->intitializeFromParentImage(parent);
177  return retval;
178 }
179 
190 ImagePtr convertImageToUnsigned(PatientModelServicePtr dataManager, ImagePtr image, vtkImageDataPtr suggestedConvertedVolume, bool verbose)
191 {
192  vtkImageDataPtr input = image->getBaseVtkImageData();
193 
194  if (input->GetScalarTypeMin() >= 0)
195  return image;
196 
197  // start by shifting up to zero
198  int shift = -input->GetScalarRange()[0];
199  // if CT: always shift by 1024 (houndsfield units definition)
200  if (image->getModality().contains("CT", Qt::CaseInsensitive))
201  shift = 1024;
202 
203  vtkImageDataPtr convertedImageData = suggestedConvertedVolume; // use input if given
204 
205  // convert volume
206  if (!convertedImageData)
207  {
208  vtkImageShiftScalePtr cast = vtkImageShiftScalePtr::New();
209  cast->SetInputData(input);
210  cast->ClampOverflowOn();
211 
212  cast->SetShift(shift);
213 
214  // total intensity range of voxels:
215  double range = input->GetScalarRange()[1] - input->GetScalarRange()[0];
216 
217  // to to fit within smallest type
218  if (range <= VTK_UNSIGNED_SHORT_MAX-VTK_UNSIGNED_SHORT_MIN)
219  cast->SetOutputScalarType(VTK_UNSIGNED_SHORT);
220  else if (range <= VTK_UNSIGNED_INT_MAX-VTK_UNSIGNED_INT_MIN)
221  cast->SetOutputScalarType(VTK_UNSIGNED_INT);
222  // else if (range <= VTK_UNSIGNED_LONG_MAX-VTK_UNSIGNED_LONG_MIN) // not supported by vtk - it seems (crash in rendering)
223  // cast->SetOutputScalarType(VTK_UNSIGNED_LONG);
224  else
225  cast->SetOutputScalarType(VTK_UNSIGNED_INT);
226 
227  cast->Update();
228  if (verbose)
229  report(QString("Converting image %1 from %2 to %3, shift=%4")
230  .arg(image->getName())
231  .arg(input->GetScalarTypeAsString())
232  .arg(cast->GetOutput()->GetScalarTypeAsString())
233  .arg(shift));
234  convertedImageData = cast->GetOutput();
235  }
236 
237  ImagePtr retval = createDerivedImage(dataManager,
238  image->getUid()+"_u", image->getName()+" u",
239  convertedImageData, image);
240 
241  ImageTF3DPtr TF3D = retval->getTransferFunctions3D()->createCopy();
242  ImageLUT2DPtr LUT2D = retval->getLookupTable2D()->createCopy();
243 
244  TF3D->shift(shift);
245  LUT2D->shift(shift);
246 
247  retval->setLookupTable2D(LUT2D);
248  retval->setTransferFunctions3D(TF3D);
249 
250  return retval;
251 }
252 
253 std::map<std::string, std::string> getDisplayFriendlyInfo(ImagePtr image)
254 {
255  std::map<std::string, std::string> retval;
256  if(!image)
257  return retval;
258 
259  //image
260  retval["Filename"] = image->getFilename().toStdString();
261  retval["Coordinate system"] = image->getCoordinateSystem().toString().toStdString();
262  retval["Image type"] = image->getImageType().toStdString();
263  retval["Scalar minimum"] = string_cast(image->getMin());
264  retval["Scalar maximum"] = string_cast(image->getMax());
265  retval["Range (max - min)"] = string_cast(image->getRange());
266  retval["Maximum alpha value"] = string_cast(image->getMaxAlphaValue());
267  retval["VTK type min value"] = string_cast(image->getVTKMinValue());
268  retval["VTK type max value"] = string_cast(image->getVTKMaxValue());
269  retval["Modality"] = image->getModality().toStdString();
270  retval["Name"] = image->getName().toStdString();
271  retval["Parent space"] = image->getParentSpace().toStdString();
272  retval["Shading"] = image->getShadingOn() ? "on" : "off";
273  retval["Space"] = image->getSpace().toStdString();
274  retval["Type"] = image->getType().toStdString();
275  retval["Uid"] = image->getUid().toStdString();
276  retval["Acquisition time"] = string_cast(image->getAcquisitionTime().toString(timestampSecondsFormatNice()));
277  retval["Voxels with min value"] = string_cast(calculateNumVoxelsWithMinValue(image));
278  retval["Voxels with max value"] = string_cast(calculateNumVoxelsWithMaxValue(image));
279  retval["rMd"] = matrixAsSingleLineString(image->get_rMd());
280 
281  std::map<std::string, std::string> volumeMap = getDisplayFriendlyInfo(image->getBaseVtkImageData());
282  retval.insert(volumeMap.begin(), volumeMap.end());
283 
284  return retval;
285 }
286 
287 std::map<std::string, std::string> getDisplayFriendlyInfo(vtkImageDataPtr image)
288 {
289  std::map<std::string, std::string> retval;
290  if(!image)
291  return retval;
292 
293  double spacing_x, spacing_y, spacing_z;
294  image->GetSpacing(spacing_x, spacing_y, spacing_z);
295  retval["Spacing"] = string_cast(spacing_x)+" mm , "+string_cast(spacing_y)+" mm , "+string_cast(spacing_z)+" mm ";
296  int dims[3];
297  image->GetDimensions(dims);
298  retval["Dimensions"] = string_cast(dims[0])+" , "+string_cast(dims[1])+" , "+string_cast(dims[2]);
299  retval["Size"] = string_cast(dims[0]*spacing_x)+" mm , "+string_cast(dims[1]*spacing_y)+" mm, "+string_cast(dims[2]*spacing_z)+" mm";
300  float actualMemorySizeKB = (float)image->GetActualMemorySize();
301  retval["Actual memory size"] = string_cast(actualMemorySizeKB/(1024*1024))+" GB, "+string_cast(actualMemorySizeKB/1024)+" MB, "+string_cast(actualMemorySizeKB)+" kB"+string_cast(actualMemorySizeKB*1024)+" bytes";
302  retval["Scalar components"] = string_cast(image->GetNumberOfScalarComponents());
303  retval["Number of components for points"] = string_cast(image->GetPointData()->GetNumberOfComponents());
304  retval["Scalar type"] = string_cast(image->GetScalarTypeAsString());
305  retval["Scalar size"] = string_cast(image->GetScalarSize());
306  int extent[6];
307  image->GetExtent(extent);
308  retval["Extent"] = string_cast(extent[0])+" , "+string_cast(extent[1])+" , "+string_cast(extent[2])+" , "+string_cast(extent[3])+" , "+string_cast(extent[4])+" , "+string_cast(extent[5]);
309 
310  return retval;
311 }
312 
313 void printDisplayFriendlyInfo(std::map<std::string, std::string> map)
314 {
315  report("----- DisplayFriendlyInfo -----");
316  std::map<std::string, std::string>::iterator it;
317  for(it = map.begin(); it != map.end(); ++it)
318  {
319  QString message((it->first+": "+it->second).c_str());
320  report(message);
321  }
322  report("-------------------------------");
323 }
324 
326 {
327  return static_cast<int*>(image->getHistogram()->GetOutput()->GetScalarPointer(image->getRange(), 0, 0))[0];
328 }
330 {
331  return static_cast<int*>(image->getHistogram()->GetOutput()->GetScalarPointer(0,0,0))[0];
332 }
333 
335 {
336  if (data.empty())
337  return DoubleBoundingBox3D(0, 0, 0, 0, 0, 0);
338 
339  std::vector<Vector3D> corners_r;
340 
341  for (unsigned i = 0; i < data.size(); ++i)
342  {
343  Transform3D qMd = qMr * data[i]->get_rMd();
344  DoubleBoundingBox3D bb = data[i]->boundingBox();
345 
346  corners_r.push_back(qMd.coord(bb.corner(0, 0, 0)));
347  corners_r.push_back(qMd.coord(bb.corner(0, 0, 1)));
348  corners_r.push_back(qMd.coord(bb.corner(0, 1, 0)));
349  corners_r.push_back(qMd.coord(bb.corner(0, 1, 1)));
350  corners_r.push_back(qMd.coord(bb.corner(1, 0, 0)));
351  corners_r.push_back(qMd.coord(bb.corner(1, 0, 1)));
352  corners_r.push_back(qMd.coord(bb.corner(1, 1, 0)));
353  corners_r.push_back(qMd.coord(bb.corner(1, 1, 1)));
354  }
355 
357  return bb_sigma;
358 }
359 
360 DoubleBoundingBox3D findEnclosingBoundingBox(std::vector<ImagePtr> images, Transform3D qMr)
361 {
362  std::vector<DataPtr> datas(images.size());
363  for (unsigned i = 0; i < images.size(); ++i)
364  datas[i] = images[i];
365  return findEnclosingBoundingBox(datas, qMr);
366 }
367 
369 {
370  vtkImageDataPtr retval = image;
371 
372  //vtkImageLuminance demands 3 components
374 
375  if (input->GetNumberOfScalarComponents() > 2)
376  {
377  vtkSmartPointer<vtkImageLuminance> luminance = vtkSmartPointer<vtkImageLuminance>::New();
378  luminance->SetInputData(input);
379  luminance->Update();
380  retval = luminance->GetOutput();
381 // retval->Update();
382  }
383  return retval;
384 }
385 
387 {
388  vtkImageDataPtr retval = image;
389 
390  if (image->GetNumberOfScalarComponents() >= 4)
391  {
392  vtkImageAppendComponentsPtr merger = vtkImageAppendComponentsPtr::New();
393  vtkImageExtractComponentsPtr splitterRGBA = vtkImageExtractComponentsPtr::New();
394  splitterRGBA->SetInputData(image);
395  splitterRGBA->SetComponents(0, 1, 2);
396  merger->AddInputConnection(splitterRGBA->GetOutputPort());
397 
398  merger->Update();
399  retval = merger->GetOutput();
400  }
401  return retval;
402 }
403 
404 vtkImageDataPtr convertImageDataTo8Bit(vtkImageDataPtr image, double windowWidth, double windowLevel)
405 {
406  vtkImageDataPtr retval = image;
407  if (image->GetScalarSize() > 1)
408  {
409  vtkImageShiftScalePtr imageCast = vtkImageShiftScalePtr::New();
410  imageCast->SetInputData(image);
411 
412 // double scalarMax = windowWidth/2.0 + windowLevel;
413  double scalarMin = windowWidth/2.0 - windowLevel;
414 
415  double addToScalarValue = -scalarMin;
416  double multiplyToScalarValue = 255/windowWidth;
417 
418  imageCast->SetShift(addToScalarValue);
419  imageCast->SetScale(multiplyToScalarValue);
420  imageCast->SetOutputScalarTypeToUnsignedChar();
421  imageCast->ClampOverflowOn();
422  imageCast->Update();
423  retval = imageCast->GetOutput();
424 // retval->Update();
425  }
426  return retval;
427 }
428 
430 {
431  image->Modified();
432  image->GetPointData()->Modified();
433  image->GetPointData()->GetScalars()->Modified();
434 }
435 
436 } // namespace cx
std::map< std::string, std::string > getDisplayFriendlyInfo(MeshPtr mesh)
Vector3D corner(int x, int y, int z) const
vtkSmartPointer< class vtkImageShiftScale > vtkImageShiftScalePtr
vtkImageDataPtr generateVtkImageDataUnsignedShort(Eigen::Array3i dim, Vector3D spacing, const unsigned short initValue, int components)
virtual void setVtkImageData(const vtkImageDataPtr &data, bool resetTransferFunctions=true)
Definition: cxImage.cpp:289
Transform3D Transform3D
Transform3D is a representation of an affine 3D transform.
std::string matrixAsSingleLineString(cx::Transform3D transform)
boost::shared_ptr< class Image > ImagePtr
Definition: cxDicomWidget.h:48
QString timestampSecondsFormatNice()
Definition: cxTime.cpp:47
int calculateNumVoxelsWithMinValue(ImagePtr image)
Find number of voxels containing min scalar value.
DoubleBoundingBox3D findEnclosingBoundingBox(std::vector< DataPtr > data, Transform3D qMr)
vtkImageDataPtr generateVtkImageDataSignedShort(Eigen::Array3i dim, Vector3D spacing, const short initValue, int components)
static DoubleBoundingBox3D fromCloud(std::vector< Vector3D > cloud)
std::string string_cast(const T &val)
ImagePtr convertImageToUnsigned(PatientModelServicePtr dataManager, ImagePtr image, vtkImageDataPtr suggestedConvertedVolume, bool verbose)
ImagePtr createDerivedImage(PatientModelServicePtr dataManager, QString uid, QString name, vtkImageDataPtr raw, ImagePtr parent)
boost::shared_ptr< class ImageLUT2D > ImageLUT2DPtr
vtkSmartPointer< class vtkImageAppendComponents > vtkImageAppendComponentsPtr
vtkImageDataPtr convertImageDataTo8Bit(vtkImageDataPtr image, double windowWidth, double windowLevel)
Have never been used or tested. Create a test for it.
vtkImageDataPtr convertFrom4To3Components(vtkImageDataPtr image)
A volumetric data set.
Definition: cxImage.h:66
boost::shared_ptr< class PatientModelService > PatientModelServicePtr
void fillShortImageDataWithGradient(vtkImageDataPtr data, int maxValue)
int calculateNumVoxelsWithMaxValue(ImagePtr image)
Find number of voxels containing max scalar value.
vtkSmartPointer< vtkDoubleArray > vtkDoubleArrayPtr
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.
vtkImageDataPtr generateVtkImageData(Eigen::Array3i dim, Vector3D spacing, const unsigned char initValue, int components)
Eigen::Vector3d Vector3D
Vector3D is a representation of a point or vector in 3D.
Definition: cxVector3D.h:63
void fill(Eigen::Affine3d *self, vtkMatrix4x4Ptr m)
vtkSmartPointer< class vtkImageExtractComponents > vtkImageExtractComponentsPtr
vtkImageDataPtr convertImageDataToGrayScale(vtkImageDataPtr image)
void report(QString msg)
Definition: cxLogger.cpp:90
void setDeepModified(vtkImageDataPtr image)
vtkSmartPointer< class vtkImageData > vtkImageDataPtr
void printDisplayFriendlyInfo(std::map< std::string, std::string > map)
boost::shared_ptr< class ImageTF3D > ImageTF3DPtr
vtkImageDataPtr generateVtkImageDataDouble(Eigen::Array3i dim, Vector3D spacing, double initValue)
Namespace for all CustusX production code.