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