NorMIT-nav  16.5
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["VTK type min value"] = string_cast(image->getVTKMinValue());
266  retval["VTK type max value"] = string_cast(image->getVTKMaxValue());
267  retval["Modality"] = image->getModality().toStdString();
268  retval["Name"] = image->getName().toStdString();
269  retval["Parent space"] = image->getParentSpace().toStdString();
270  retval["Shading"] = image->getShadingOn() ? "on" : "off";
271  retval["Space"] = image->getSpace().toStdString();
272  retval["Type"] = image->getType().toStdString();
273  retval["Uid"] = image->getUid().toStdString();
274  retval["Acquisition time"] = string_cast(image->getAcquisitionTime().toString(timestampSecondsFormatNice()));
275  retval["Voxels with min value"] = string_cast(calculateNumVoxelsWithMinValue(image));
276  retval["Voxels with max value"] = string_cast(calculateNumVoxelsWithMaxValue(image));
277  retval["rMd"] = string_cast(image->get_rMd());
278 
279  std::map<std::string, std::string> volumeMap = getDisplayFriendlyInfo(image->getBaseVtkImageData());
280  retval.insert(volumeMap.begin(), volumeMap.end());
281 
282  return retval;
283 }
284 
285 std::map<std::string, std::string> getDisplayFriendlyInfo(vtkImageDataPtr image)
286 {
287  std::map<std::string, std::string> retval;
288  if(!image)
289  return retval;
290 
291  double spacing_x, spacing_y, spacing_z;
292  image->GetSpacing(spacing_x, spacing_y, spacing_z);
293  retval["Spacing"] = string_cast(spacing_x)+" mm , "+string_cast(spacing_y)+" mm , "+string_cast(spacing_z)+" mm ";
294  int dims[3];
295  image->GetDimensions(dims);
296  retval["Dimensions"] = string_cast(dims[0])+" , "+string_cast(dims[1])+" , "+string_cast(dims[2]);
297  retval["Size"] = string_cast(dims[0]*spacing_x)+" mm , "+string_cast(dims[1]*spacing_y)+" mm, "+string_cast(dims[2]*spacing_z)+" mm";
298  float actualMemorySizeKB = (float)image->GetActualMemorySize();
299  retval["Actual memory size"] = string_cast(actualMemorySizeKB/(1024*1024))+" GB, "+string_cast(actualMemorySizeKB/1024)+" MB, "+string_cast(actualMemorySizeKB)+" kB"+string_cast(actualMemorySizeKB*1024)+" bytes";
300  retval["Scalar components"] = string_cast(image->GetNumberOfScalarComponents());
301  retval["Number of components for points"] = string_cast(image->GetPointData()->GetNumberOfComponents());
302  retval["Scalar type"] = string_cast(image->GetScalarTypeAsString());
303  retval["Scalar size"] = string_cast(image->GetScalarSize());
304  int extent[6];
305  image->GetExtent(extent);
306  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]);
307 
308  return retval;
309 }
310 
311 void printDisplayFriendlyInfo(std::map<std::string, std::string> map)
312 {
313  report("----- DisplayFriendlyInfo -----");
314  std::map<std::string, std::string>::iterator it;
315  for(it = map.begin(); it != map.end(); ++it)
316  {
317  QString message((it->first+": "+it->second).c_str());
318  report(message);
319  }
320  report("-------------------------------");
321 }
322 
324 {
325  return static_cast<int*>(image->getHistogram()->GetOutput()->GetScalarPointer(image->getRange(), 0, 0))[0];
326 }
328 {
329  return static_cast<int*>(image->getHistogram()->GetOutput()->GetScalarPointer(0,0,0))[0];
330 }
331 
333 {
334  if (data.empty())
335  return DoubleBoundingBox3D(0, 0, 0, 0, 0, 0);
336 
337  std::vector<Vector3D> corners_r;
338 
339  for (unsigned i = 0; i < data.size(); ++i)
340  {
341  Transform3D qMd = qMr * data[i]->get_rMd();
342  DoubleBoundingBox3D bb = data[i]->boundingBox();
343 
344  corners_r.push_back(qMd.coord(bb.corner(0, 0, 0)));
345  corners_r.push_back(qMd.coord(bb.corner(0, 0, 1)));
346  corners_r.push_back(qMd.coord(bb.corner(0, 1, 0)));
347  corners_r.push_back(qMd.coord(bb.corner(0, 1, 1)));
348  corners_r.push_back(qMd.coord(bb.corner(1, 0, 0)));
349  corners_r.push_back(qMd.coord(bb.corner(1, 0, 1)));
350  corners_r.push_back(qMd.coord(bb.corner(1, 1, 0)));
351  corners_r.push_back(qMd.coord(bb.corner(1, 1, 1)));
352  }
353 
355  return bb_sigma;
356 }
357 
358 DoubleBoundingBox3D findEnclosingBoundingBox(std::vector<ImagePtr> images, Transform3D qMr)
359 {
360  std::vector<DataPtr> datas(images.size());
361  for (unsigned i = 0; i < images.size(); ++i)
362  datas[i] = images[i];
363  return findEnclosingBoundingBox(datas, qMr);
364 }
365 
367 {
368  vtkImageDataPtr retval = image;
369  if (image->GetNumberOfScalarComponents() > 2)
370  {
371  vtkSmartPointer<vtkImageLuminance> luminance = vtkSmartPointer<vtkImageLuminance>::New();
372  luminance->SetInputData(image);
373  luminance->Update();
374  retval = luminance->GetOutput();
375 // retval->Update();
376  }
377  return retval;
378 }
379 
380 vtkImageDataPtr convertImageDataTo8Bit(vtkImageDataPtr image, double windowWidth, double windowLevel)
381 {
382  vtkImageDataPtr retval = image;
383  if (image->GetScalarSize() > 1)
384  {
385  vtkImageShiftScalePtr imageCast = vtkImageShiftScalePtr::New();
386  imageCast->SetInputData(image);
387 
388 // double scalarMax = windowWidth/2.0 + windowLevel;
389  double scalarMin = windowWidth/2.0 - windowLevel;
390 
391  double addToScalarValue = -scalarMin;
392  double multiplyToScalarValue = 255/windowWidth;
393 
394  imageCast->SetShift(addToScalarValue);
395  imageCast->SetScale(multiplyToScalarValue);
396  imageCast->SetOutputScalarTypeToUnsignedChar();
397  imageCast->ClampOverflowOn();
398  imageCast->Update();
399  retval = imageCast->GetOutput();
400 // retval->Update();
401  }
402  return retval;
403 }
404 
406 {
407  image->Modified();
408  image->GetPointData()->Modified();
409  image->GetPointData()->GetScalars()->Modified();
410 }
411 
412 } // 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)