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