Medical Imaging Interaction Toolkit  2016.11.0
Medical Imaging Interaction Toolkit
itkShCoefficientImageImporter.cpp
Go to the documentation of this file.
1 /*===================================================================
2 
3 The Medical Imaging Interaction Toolkit (MITK)
4 
5 Copyright (c) German Cancer Research Center,
6 Division of Medical and Biological Informatics.
7 All rights reserved.
8 
9 This software is distributed WITHOUT ANY WARRANTY; without
10 even the implied warranty of MERCHANTABILITY or FITNESS FOR
11 A PARTICULAR PURPOSE.
12 
13 See LICENSE.txt or http://www.mitk.org for details.
14 
15 ===================================================================*/
16 #ifndef __itkShCoefficientImageImporter_cpp
17 #define __itkShCoefficientImageImporter_cpp
18 
19 #include <time.h>
20 #include <stdio.h>
21 #include <stdlib.h>
22 
24 #include <itkImageRegionIterator.h>
25 #include <boost/math/special_functions.hpp>
26 
27 using namespace boost::math;
28 
29 namespace itk {
30 
31 template< class PixelType, int ShOrder >
33  : m_Toolkit(FSL)
34 {
35  m_ShBasis.set_size(QBALL_ODFSIZE, (ShOrder+1)*(ShOrder+2)/2);
36 }
37 
38 template< class PixelType, int ShOrder >
41 {
42  CalcShBasis();
43  if (m_InputImage.IsNull())
44  return;
45 
46  Vector<float, 4> spacing4 = m_InputImage->GetSpacing();
47  Point<float, 4> origin4 = m_InputImage->GetOrigin();
48  Matrix<double, 4, 4> direction4 = m_InputImage->GetDirection();
49  ImageRegion<4> imageRegion4 = m_InputImage->GetLargestPossibleRegion();
50 
51  Vector<double, 3> spacing3;
52  Point<float, 3> origin3;
53  Matrix<double, 3, 3> direction3;
54  ImageRegion<3> imageRegion3;
55 
56  spacing3[0] = spacing4[0]; spacing3[1] = spacing4[1]; spacing3[2] = spacing4[2];
57  origin3[0] = origin4[0]; origin3[1] = origin4[1]; origin3[2] = origin4[2];
58  for (int r=0; r<3; r++)
59  for (int c=0; c<3; c++)
60  direction3[r][c] = direction4[r][c];
61  imageRegion3.SetSize(0, imageRegion4.GetSize()[0]);
62  imageRegion3.SetSize(1, imageRegion4.GetSize()[1]);
63  imageRegion3.SetSize(2, imageRegion4.GetSize()[2]);
64 
65  m_QballImage = QballImageType::New();
66  m_QballImage->SetSpacing( spacing3 );
67  m_QballImage->SetOrigin( origin3 );
68  m_QballImage->SetDirection( direction3 );
69  m_QballImage->SetRegions( imageRegion3 );
70  m_QballImage->Allocate();
71  Vector< PixelType, QBALL_ODFSIZE > nullVec1; nullVec1.Fill(0.0);
72  m_QballImage->FillBuffer(nullVec1);
73 
74  m_CoefficientImage = CoefficientImageType::New();
75  m_CoefficientImage->SetSpacing( spacing3 );
76  m_CoefficientImage->SetOrigin( origin3 );
77  m_CoefficientImage->SetDirection( direction3 );
78  m_CoefficientImage->SetRegions( imageRegion3 );
79  m_CoefficientImage->Allocate();
80  Vector< PixelType, (ShOrder*ShOrder + ShOrder + 2)/2 + ShOrder > nullVec2; nullVec2.Fill(0.0);
81  m_CoefficientImage->FillBuffer(nullVec2);
82 
83  typedef ImageRegionConstIterator< InputImageType > InputIteratorType;
84  int x = imageRegion4.GetSize(0);
85  int y = imageRegion4.GetSize(1);
86  int z = imageRegion4.GetSize(2);
87  int numCoeffs = imageRegion4.GetSize(3);
88 
89  for (int a=0; a<x; a++)
90  for (int b=0; b<y; b++)
91  for (int c=0; c<z; c++)
92  {
93  vnl_matrix<double> coeffs((ShOrder*ShOrder + ShOrder + 2)/2 + ShOrder,1);
94 
95  typename InputImageType::IndexType index;
96  index.SetElement(0,a);
97  index.SetElement(1,b);
98  index.SetElement(2,c);
100  for (int d=0; d<numCoeffs; d++)
101  {
102  index.SetElement(3,d);
103  pix[d] = m_InputImage->GetPixel(index);
104  coeffs[d][0] = pix[d];
105  }
106  typename CoefficientImageType::IndexType index2;
107  index2.SetElement(0,a);
108  index2.SetElement(1,b);
109  index2.SetElement(2,c);
110  m_CoefficientImage->SetPixel(index2, pix);
111 
112  typename QballImageType::PixelType pix2;
113  vnl_matrix<double> odf = m_ShBasis*coeffs;
114  for (int d=0; d<QBALL_ODFSIZE; d++)
115  pix2[d] = odf(d,0)*M_PI*4/QBALL_ODFSIZE;
116 
117  m_QballImage->SetPixel(index2,pix2);
118  }
119 }
120 
121 // generate spherical harmonic values of the desired order for each input direction
122 template< class PixelType, int ShOrder >
125 {
126  vnl_matrix_fixed<double, 2, QBALL_ODFSIZE> sphCoords = GetSphericalOdfDirections();
127  int j, m; double mag, plm;
128 
129  for (int p=0; p<QBALL_ODFSIZE; p++)
130  {
131  j=0;
132  for (int l=0; l<=ShOrder; l=l+2)
133  for (m=-l; m<=l; m++)
134  {
135  switch (m_Toolkit)
136  {
137  case FSL:
138  plm = legendre_p<double>(l,abs(m),cos(sphCoords(0,p)));
139  mag = sqrt((double)(2*l+1)/(4.0*M_PI)*factorial<double>(l-abs(m))/factorial<double>(l+abs(m)))*plm;
140  if (m<0)
141  m_ShBasis(p,j) = sqrt(2.0)*mag*cos(-m*sphCoords(1,p));
142  else if (m==0)
143  m_ShBasis(p,j) = mag;
144  else
145  m_ShBasis(p,j) = pow(-1.0, m)*sqrt(2.0)*mag*sin(m*sphCoords(1,p));
146  break;
147  case MRTRIX:
148  plm = legendre_p<double>(l,abs(m),-cos(sphCoords(0,p)));
149  mag = sqrt((double)(2*l+1)/(4.0*M_PI)*factorial<double>(l-abs(m))/factorial<double>(l+abs(m)))*plm;
150  if (m>0)
151  m_ShBasis(p,j) = mag*cos(m*sphCoords(1,p));
152  else if (m==0)
153  m_ShBasis(p,j) = mag;
154  else
155  m_ShBasis(p,j) = mag*sin(-m*sphCoords(1,p));
156  break;
157  }
158 
159  j++;
160  }
161  }
162 }
163 
164 // convert cartesian to spherical coordinates
165 template< class PixelType, int ShOrder >
166 vnl_matrix_fixed<double, 2, QBALL_ODFSIZE> ShCoefficientImageImporter< PixelType, ShOrder >
168 {
170  vnl_matrix_fixed<double, 3, QBALL_ODFSIZE>* dir = odf.GetDirections();
171  vnl_matrix_fixed<double, 2, QBALL_ODFSIZE> sphCoords;
172 
173  for (int i=0; i<QBALL_ODFSIZE; i++)
174  {
175  double mag = dir->get_column(i).magnitude();
176 
177  if( mag<mitk::eps )
178  {
179  sphCoords(0,i) = M_PI/2; // theta
180  sphCoords(1,i) = M_PI/2; // phi
181  }
182  else
183  {
184  sphCoords(0,i) = acos(dir->get(2,i)/mag); // theta
185  sphCoords(1,i) = atan2(dir->get(1,i), dir->get(0,i)); // phi
186  }
187  }
188  return sphCoords;
189 }
190 
191 }
192 
193 #endif // __itkShCoefficientImageImporter_cpp
Represents an ODF for Q-Ball imaging.
#define QBALL_ODFSIZE
MITKCORE_EXPORT const ScalarType eps
unsigned short PixelType
vnl_matrix_fixed< double, 2, QBALL_ODFSIZE > GetSphericalOdfDirections()
static itkEventMacro(BoundingShapeInteractionEvent, itk::AnyEvent) class MITKBOUNDINGSHAPE_EXPORT BoundingShapeInteractor Pointer New()
Basic interaction methods for mitk::GeometryData.