Fraxinus  2023.01.05-dev+develop.0da12
An IGT application
cxIGTLinkConversionImage.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 
14  Portions (c) Copyright 2008-2009 Brigham and Women's Hospital (BWH) All Rights Reserved.
15 
16  See Doc/copyright/copyright.txt
17  or http://www.slicer.org/copyright/copyright.txt for details.
18 
19  Program: 3D Slicer
20  Module: $HeadURL: http://svn.slicer.org/Slicer3/trunk/Modules/OpenIGTLinkIF/vtkIGTLToMRMLImage.h $
21  Date: $Date: 2010-11-23 00:58:13 -0500 (Tue, 23 Nov 2010) $
22  Version: $Revision: 15552 $
23 
24 ==========================================================================*/
25 #include "cxIGTLinkConversionImage.h"
26 #include "vtkImageData.h"
27 
28 #include <igtl_util.h>
29 #include "cxLogger.h"
30 #include "cxRegistrationTransform.h"
31 #include "cxIGTLinkConversionBase.h"
32 
33 
34 namespace cx
35 {
36 
37 igtl::ImageMessage::Pointer IGTLinkConversionImage::encode(ImagePtr image, PATIENT_COORDINATE_SYSTEM externalSpace)
38 {
39  igtl::ImageMessage::Pointer retval = igtl::ImageMessage::New();
40 
41  retval->SetDeviceName(cstring_cast(image->getName()));
42  IGTLinkConversionBase().encode_timestamp(image->getAcquisitionTime(), retval);
43  this->encode_vtkImageData(image->getBaseVtkImageData(), retval);
44  this->encode_rMd(image, retval, externalSpace);
45 
46  return retval;
47 }
48 
49 ImagePtr IGTLinkConversionImage::decode(igtl::ImageMessage *msg)
50 {
51  vtkImageDataPtr vtkImage = this->decode_vtkImageData(msg);
52  QDateTime timestamp = IGTLinkConversionBase().decode_timestamp(msg);
53  QString deviceName = msg->GetDeviceName();
54 
55  ImagePtr retval(new Image(deviceName, vtkImage));
56  retval->setAcquisitionTime(timestamp);
57  this->decode_rMd(msg, retval);
58 
59  return retval;
60 }
61 
62 
63 
64 namespace { // unnamed namespace
65 
66 //---------------------------------------------------------------------------
67 // Stream copy + byte swap
68 //---------------------------------------------------------------------------
69 int swapCopy16(igtlUint16 * dst, igtlUint16 * src, int n)
70 {
71  igtlUint16 * esrc = src + n;
72  while (src < esrc)
73  {
74  *dst = BYTE_SWAP_INT16(*src);
75  dst ++;
76  src ++;
77  }
78  return 1;
79 }
80 
81 int swapCopy32(igtlUint32 * dst, igtlUint32 * src, int n)
82 {
83  igtlUint32 * esrc = src + n;
84  while (src < esrc)
85  {
86  *dst = BYTE_SWAP_INT32(*src);
87  dst ++;
88  src ++;
89  }
90  return 1;
91 }
92 
93 int swapCopy64(igtlUint64 * dst, igtlUint64 * src, int n)
94 {
95  igtlUint64 * esrc = src + n;
96  while (src < esrc)
97  {
98  *dst = BYTE_SWAP_INT64(*src);
99  dst ++;
100  src ++;
101  }
102  return 1;
103 }
104 } // unnamed namespace
105 
106 vtkImageDataPtr IGTLinkConversionImage::decode_vtkImageData(igtl::ImageMessage *imgMsg)
107 {
108  // NOTE: This method is mostly a copy-paste from Slicer.
109  // MRML and coordinate stuff are removed.
110  // Avoid refactoring the internals, as it is important to keep the code similar to the origin.
111 
112  // Retrieve the image data
113  int size[3]; // image dimension
114  float spacing[3]; // spacing (mm/pixel)
115  int svsize[3]; // sub-volume size
116  int svoffset[3]; // sub-volume offset
117  int scalarType; // VTK scalar type
118  int numComponents; // number of scalar components
119  int endian;
120 
121  scalarType = IGTLToVTKScalarType( imgMsg->GetScalarType() );
122  endian = imgMsg->GetEndian();
123  imgMsg->GetDimensions(size);
124  imgMsg->GetSpacing(spacing);
125  numComponents = imgMsg->GetNumComponents();
126  imgMsg->GetSubVolume(svsize, svoffset);
127  // imgMsg->GetMatrix(matrix);
128 
129  // check if the IGTL data fits to the current MRML node
130  int sizeInNode[3]={0,0,0};
131  int scalarTypeInNode=VTK_VOID;
132  int numComponentsInNode=0;
133  // Get vtk image from MRML node
134  vtkSmartPointer<vtkImageData> imageData = vtkSmartPointer<vtkImageData>::New();
135  imageData->SetDimensions(size[0], size[1], size[2]);
136  imageData->SetExtent(0, size[0]-1, 0, size[1]-1, 0, size[2]-1);
137  imageData->SetOrigin(0.0, 0.0, 0.0);
138 // imageData->SetSpacing(1.0, 1.0, 1.0); // Slicer inserts spacing into its IKTtoRAS matrix, we dont.
139  imageData->SetSpacing(spacing[0], spacing[1], spacing[2]);
140  imageData->AllocateScalars(scalarType, numComponents);
141 
142  // Check scalar size
143  int scalarSize = imgMsg->GetScalarSize();
144 
145  int fByteSwap = 0;
146  // Check if bytes-swap is required
147  if (scalarSize > 1 &&
148  ((igtl_is_little_endian() && endian == igtl::ImageMessage::ENDIAN_BIG) ||
149  (!igtl_is_little_endian() && endian == igtl::ImageMessage::ENDIAN_LITTLE)))
150  {
151  // Needs byte swap
152  fByteSwap = 1;
153  }
154 
155  if (imgMsg->GetImageSize() == imgMsg->GetSubVolumeImageSize())
156  {
157  // In case that volume size == sub-volume size,
158  // image is read directly to the memory area of vtkImageData
159  // for better performance.
160  if (fByteSwap)
161  {
162  switch (scalarSize)
163  {
164  case 2:
165  swapCopy16((igtlUint16 *)imageData->GetScalarPointer(),
166  (igtlUint16 *)imgMsg->GetScalarPointer(),
167  imgMsg->GetSubVolumeImageSize() / 2);
168  break;
169  case 4:
170  swapCopy32((igtlUint32 *)imageData->GetScalarPointer(),
171  (igtlUint32 *)imgMsg->GetScalarPointer(),
172  imgMsg->GetSubVolumeImageSize() / 4);
173  break;
174  case 8:
175  swapCopy64((igtlUint64 *)imageData->GetScalarPointer(),
176  (igtlUint64 *)imgMsg->GetScalarPointer(),
177  imgMsg->GetSubVolumeImageSize() / 8);
178  break;
179  default:
180  break;
181  }
182  }
183  else
184  {
185  memcpy(imageData->GetScalarPointer(),
186  imgMsg->GetScalarPointer(), imgMsg->GetSubVolumeImageSize());
187  }
188  }
189  else
190  {
191  // In case of volume size != sub-volume size,
192  // image is loaded into ImageReadBuffer, then copied to
193  // the memory area of vtkImageData.
194  char* imgPtr = (char*) imageData->GetScalarPointer();
195  char* bufPtr = (char*) imgMsg->GetScalarPointer();
196  int sizei = size[0];
197  int sizej = size[1];
198  //int sizek = size[2];
199  int subsizei = svsize[0];
200 
201  int bg_i = svoffset[0];
202  //int ed_i = bg_i + svsize[0];
203  int bg_j = svoffset[1];
204  int ed_j = bg_j + svsize[1];
205  int bg_k = svoffset[2];
206  int ed_k = bg_k + svsize[2];
207 
208  if (fByteSwap)
209  {
210  switch (scalarSize)
211  {
212  case 2:
213  for (int k = bg_k; k < ed_k; k ++)
214  {
215  for (int j = bg_j; j < ed_j; j ++)
216  {
217  swapCopy16((igtlUint16 *)&imgPtr[(sizei*sizej*k + sizei*j + bg_i)*scalarSize],
218  (igtlUint16 *)bufPtr,
219  subsizei);
220  bufPtr += subsizei*scalarSize;
221  }
222  }
223  break;
224  case 4:
225  for (int k = bg_k; k < ed_k; k ++)
226  {
227  for (int j = bg_j; j < ed_j; j ++)
228  {
229  swapCopy32((igtlUint32 *)&imgPtr[(sizei*sizej*k + sizei*j + bg_i)*scalarSize],
230  (igtlUint32 *)bufPtr,
231  subsizei);
232  bufPtr += subsizei*scalarSize;
233  }
234  }
235  break;
236  case 8:
237  for (int k = bg_k; k < ed_k; k ++)
238  {
239  for (int j = bg_j; j < ed_j; j ++)
240  {
241  swapCopy64((igtlUint64 *)&imgPtr[(sizei*sizej*k + sizei*j + bg_i)*scalarSize],
242  (igtlUint64 *)bufPtr,
243  subsizei);
244  bufPtr += subsizei*scalarSize;
245  }
246  }
247  break;
248  default:
249  break;
250  }
251  }
252  else
253  {
254  for (int k = bg_k; k < ed_k; k ++)
255  {
256  for (int j = bg_j; j < ed_j; j ++)
257  {
258  memcpy(&imgPtr[(sizei*sizej*k + sizei*j + bg_i)*scalarSize],
259  bufPtr, subsizei*scalarSize);
260  bufPtr += subsizei*scalarSize;
261  }
262  }
263  }
264 
265  }
266 
267  imageData->Modified();
268  return imageData;
269 }
270 
271 void IGTLinkConversionImage::encode_vtkImageData(vtkImageDataPtr in, igtl::ImageMessage *outmsg)
272 {
273  // NOTE: This method is mostly a copy-paste from Slicer.
274  // MRML and coordinate stuff are removed.
275  // Avoid refactoring the internals, as it is important to keep the code similar to the origin.
276 
277  vtkImageDataPtr imageData = in;
278  int isize[3]; // image dimension
279  int scalarType; // scalar type
280  double *spacing; // spacing (mm/pixel)
281  int svoffset[] = {0, 0, 0}; // sub-volume offset
282  int endian;
283 
284  scalarType = imageData->GetScalarType();
285  imageData->GetDimensions(isize);
286  spacing = imageData->GetSpacing();
287  int numComponents = imageData->GetNumberOfScalarComponents();
288 
289  // Check endianness of the machine
290  endian = igtl::ImageMessage::ENDIAN_BIG;
291  if (igtl_is_little_endian())
292  {
293  endian = igtl::ImageMessage::ENDIAN_LITTLE;
294  }
295 
296  outmsg->SetDimensions(isize);
297  outmsg->SetSpacing((float)spacing[0], (float)spacing[1], (float)spacing[2]);
298  outmsg->SetScalarType(scalarType);
299  outmsg->SetEndian(endian);
300  outmsg->SetSubVolume(isize, svoffset);
301  outmsg->SetNumComponents(numComponents);
302  outmsg->AllocateScalars();
303 
304  memcpy(outmsg->GetScalarPointer(),
305  imageData->GetScalarPointer(),
306  outmsg->GetImageSize());
307 }
308 
309 
310 void IGTLinkConversionImage::decode_rMd(igtl::ImageMessage *msg, ImagePtr out)
311 {
312  Transform3D sMigtl = this->getMatrix(msg);
313 
314  Vector3D c = out->boundingBox().center();
315  Transform3D igtlMd = createTransformTranslate(-c);
316 
317  // s is the inbound reference system, i.e. LPS or RAS.
318  PATIENT_COORDINATE_SYSTEM s = this->getPatientCoordinateSystem(msg->GetCoordinateSystem());
319  Transform3D rMs = createTransformFromReferenceToExternal(s).inv();
320 
321  Transform3D rMd = rMs * sMigtl * igtlMd;
322 
323  out->get_rMd_History()->setRegistration(rMd);
324 }
325 
326 int IGTLinkConversionImage::getIgtlCoordinateSystem(PATIENT_COORDINATE_SYSTEM space) const
327 {
328  if (space==pcsRAS)
329  return igtl::ImageMessage::COORDINATE_RAS;
330  if (space==pcsLPS)
331  return igtl::ImageMessage::COORDINATE_LPS;
332 
333  return igtl::ImageMessage::COORDINATE_RAS; // default
334 }
335 
336 PATIENT_COORDINATE_SYSTEM IGTLinkConversionImage::getPatientCoordinateSystem(int igtlSpace) const
337 {
338  if (igtlSpace==igtl::ImageMessage::COORDINATE_RAS)
339  return pcsRAS;
340  if (igtlSpace==igtl::ImageMessage::COORDINATE_LPS)
341  return pcsLPS;
342 
343  return pcsRAS; // default
344 }
345 
346 void IGTLinkConversionImage::encode_rMd(ImagePtr image, igtl::ImageMessage *outmsg, PATIENT_COORDINATE_SYSTEM externalSpace)
347 {
348  Transform3D rMd = image->get_rMd();
349 
350  // NOTE: there seems to be a bug in Slicer3D: LPS not supported in OpenIGTLinkIF,
351  // thus using LPS will fail against Slicer3D.
352  // This doesn't matter much, as POLYDATA anyway must use RAS.
353 
354  // s is the outbound reference system, i.e. LPS or RAS.
355  Transform3D sMr = createTransformFromReferenceToExternal(externalSpace);
356  outmsg->SetCoordinateSystem(this->getIgtlCoordinateSystem(externalSpace));
357 
358 // if (coordinateSystem=="RAS")
359 // {
360 // // NOTE: there seems to be a bug in Slicer: LPS not supported in OpenIGTLinkIF.
361 // // This doesn't matter much, as POLYDATA anyway must use RAS.
362 // sMr = createTransformLPS2RAS();
363 // outmsg->SetCoordinateSystem(igtl::ImageMessage::COORDINATE_RAS);
364 // }
365 // else // i.e. LPS
366 // {
367 // sMr = Transform3D::Identity();
368 // outmsg->SetCoordinateSystem(igtl::ImageMessage::COORDINATE_LPS);
369 // }
370 
371  // shift origin to image center
372  Vector3D c = image->boundingBox().center();
373  Transform3D dMigtl = createTransformTranslate(c);
374  Transform3D sMigtl = sMr * rMd * dMigtl;
375 
376 // CX_LOG_CHANNEL_DEBUG("ca") << "rMd encode\n" << rMd;
377 // CX_LOG_CHANNEL_DEBUG("ca") << "sMr encode\n" << sMr;
378 // CX_LOG_CHANNEL_DEBUG("ca") << "sMigtl encode\n" << sMigtl;
379 
380  this->setMatrix(outmsg, sMigtl);
381 }
382 
383 Transform3D IGTLinkConversionImage::getMatrix(igtl::ImageMessage *msg)
384 {
385  igtl::Matrix4x4 matrix;
386  msg->GetMatrix(matrix);
387  return Transform3D::fromFloatArray(matrix);
388 }
389 
390 void IGTLinkConversionImage::setMatrix(igtl::ImageMessage *msg, Transform3D matrix)
391 {
392  igtl::Matrix4x4 m;
393  for (int r = 0; r < 4; ++r)
394  for (int c = 0; c < 4; ++c)
395  m[r][c] = matrix(r,c);
396 
397  msg->SetMatrix(m);
398 }
399 
400 int IGTLinkConversionImage::IGTLToVTKScalarType(int igtlType)
401 {
402  switch (igtlType)
403  {
404  case igtl::ImageMessage::TYPE_INT8: return VTK_CHAR;
405  case igtl::ImageMessage::TYPE_UINT8: return VTK_UNSIGNED_CHAR;
406  case igtl::ImageMessage::TYPE_INT16: return VTK_SHORT;
407  case igtl::ImageMessage::TYPE_UINT16: return VTK_UNSIGNED_SHORT;
408  case igtl::ImageMessage::TYPE_INT32: return VTK_UNSIGNED_LONG;
409  case igtl::ImageMessage::TYPE_UINT32: return VTK_UNSIGNED_LONG;
410  case igtl::ImageMessage::TYPE_FLOAT32: return VTK_FLOAT;
411  case igtl::ImageMessage::TYPE_FLOAT64: return VTK_DOUBLE;
412  default:
413  CX_LOG_CHANNEL_ERROR("igtl") << "Invalid IGTL scalar Type: "<< igtlType;
414  return VTK_VOID;
415  }
416 }
417 
418 
419 } //namespace cx
420 
421 
pcsRAS
pcsRAS
Right-Anterior-Superior, used by Slicer3D, ITK-Snap, nifti, MINC.
Definition: cxDefinitions.h:104
cx::Transform3D
Transform3D Transform3D
Transform3D is a representation of an affine 3D transform.
Definition: cxLandmarkPatientRegistrationWidget.h:33
cx::createTransformFromReferenceToExternal
cxResource_EXPORT Transform3D createTransformFromReferenceToExternal(PATIENT_COORDINATE_SYSTEM external)
Definition: cxTransform3D.cpp:271
cx::ImagePtr
boost::shared_ptr< class Image > ImagePtr
Definition: cxDicomWidget.h:27
cstring_cast
cstring_cast_Placeholder cstring_cast(const T &val)
Definition: cxTypeConversions.h:69
cx::Image
A volumetric data set.
Definition: cxImage.h:45
cx::createTransformTranslate
Transform3D createTransformTranslate(const Vector3D &translation)
Definition: cxTransform3D.cpp:164
pcsLPS
pcsLPS
Left-Posterior-Superior, used internally by CustusX, also DICOM, ITK.
Definition: cxDefinitions.h:104
cx::Vector3D
Eigen::Vector3d Vector3D
Vector3D is a representation of a point or vector in 3D.
Definition: cxVector3D.h:42
vtkImageDataPtr
vtkSmartPointer< class vtkImageData > vtkImageDataPtr
Definition: cxVideoConnectionWidget.h:30
CX_LOG_CHANNEL_ERROR
#define CX_LOG_CHANNEL_ERROR(channel)
Definition: cxLogger.h:111
cx
Namespace for all CustusX production code.
Definition: cx_dev_group_definitions.h:13