CASToR  2.0
Tomographic Reconstruction (PET/SPECT/CT)
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Macros Groups Pages
iImageConvolverStationaryGaussian.cc
Go to the documentation of this file.
1 /*
2 This file is part of CASToR.
3 
4  CASToR is free software: you can redistribute it and/or modify it under the
5  terms of the GNU General Public License as published by the Free Software
6  Foundation, either version 3 of the License, or (at your option) any later
7  version.
8 
9  CASToR is distributed in the hope that it will be useful, but WITHOUT ANY
10  WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
11  FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
12  details.
13 
14  You should have received a copy of the GNU General Public License along with
15  CASToR (in file GNU_GPL.TXT). If not, see <http://www.gnu.org/licenses/>.
16 
17 Copyright 2017-2018 all CASToR contributors listed below:
18 
19  --> current contributors: Thibaut MERLIN, Simon STUTE, Didier BENOIT, Claude COMTAT, Marina FILIPOVIC, Mael MILLARDET
20  --> past contributors: Valentin VIELZEUF
21 
22 This is CASToR version 2.0.
23 */
24 
31 #include "vImageConvolver.hh"
33 #include "oImageSpace.hh"
35 
36 // =====================================================================
37 // ---------------------------------------------------------------------
38 // ---------------------------------------------------------------------
39 // =====================================================================
40 
42 {
43  m_nbSigmas = -1.;
44  m_transFWHM = -1.;
45  m_axialFWHM = -1.;
46  m_dimKernelXY = -1;
47  m_dimKernelXYZ = -1;
48 }
49 
50 // =====================================================================
51 // ---------------------------------------------------------------------
52 // ---------------------------------------------------------------------
53 // =====================================================================
54 
56 {
57 }
58 
59 // =====================================================================
60 // ---------------------------------------------------------------------
61 // ---------------------------------------------------------------------
62 // =====================================================================
63 
65 {
66  cout << "This convolver is based on a classic stationary Gaussian kernel." << endl;
67  cout << "It can be anisotropic so that the transaxial and axial FWHM can be different." << endl;
68  cout << "One also have to choose the number of sigmas of the Gaussian distribution that will be included in the kernel." << endl;
69  cout << "The following options can be used (in this particular order when provided as a list, separated by commas):" << endl;
70  cout << " transaxial FWHM: to set the transaxial FWHM (in mm)" << endl;
71  cout << " axial FWHM: to set the axial FWHM (in mm)" << endl;
72  cout << " number of sigmas: to set the number of sigmas included in the kernel (recommendations: at least 3. and maximum 5.)" << endl;
73 }
74 
75 // =====================================================================
76 // ---------------------------------------------------------------------
77 // ---------------------------------------------------------------------
78 // =====================================================================
79 
80 int iImageConvolverStationaryGaussian::ReadConfigurationFile(const string& a_configurationFile)
81 {
82  string key_word = "";
83  // Read the transaxial FWHM option
84  key_word = "transaxial FWHM";
85  if (ReadDataASCIIFile(a_configurationFile, key_word, &m_transFWHM, 1, KEYWORD_MANDATORY))
86  {
87  Cerr("***** iImageConvolverStationaryGaussian::ReadConfigurationFile() -> Failed to get the '" << key_word << "' keyword !" << endl);
88  return 1;
89  }
90  // Read the axial FWHM option
91  key_word = "axial FWHM";
92  if (ReadDataASCIIFile(a_configurationFile, key_word, &m_axialFWHM, 1, KEYWORD_MANDATORY))
93  {
94  Cerr("***** iImageConvolverStationaryGaussian::ReadConfigurationFile() -> Failed to get the '" << key_word << "' keyword !" << endl);
95  return 1;
96  }
97  // Read the number of sigma option
98  key_word = "number of sigmas";
99  if (ReadDataASCIIFile(a_configurationFile, key_word, &m_nbSigmas, 1, KEYWORD_MANDATORY))
100  {
101  Cerr("***** iImageConvolverStationaryGaussian::ReadConfigurationFile() -> Failed to get the '" << key_word << "' keyword !" << endl);
102  return 1;
103  }
104  // Normal end
105  return 0;
106 }
107 
108 // =====================================================================
109 // ---------------------------------------------------------------------
110 // ---------------------------------------------------------------------
111 // =====================================================================
112 
114 {
115  // There are 3 floating point variables as options
116  FLTNB options[3];
117  // Read them
118  if (ReadStringOption(a_optionsList, options, 3, ",", "Stationary Gaussian convolver configuration"))
119  {
120  Cerr("***** iImageConvolverStationaryGaussian::ReadConfigurationFile() -> Failed to correctly read the list of options !" << endl);
121  return 1;
122  }
123  // Affect options
124  m_transFWHM = options[0];
125  m_axialFWHM = options[1];
126  m_nbSigmas = options[2];
127  // Normal end
128  return 0;
129 }
130 
131 // =====================================================================
132 // ---------------------------------------------------------------------
133 // ---------------------------------------------------------------------
134 // =====================================================================
135 
137 {
138  // Check number of sigmas
139  if (m_nbSigmas<=0.)
140  {
141  Cerr("***** iImageConvolverStationaryGaussian::CheckSpecificParameters() -> Number of sigmas included in the kernel must be strictly positive !" << endl);
142  return 1;
143  }
144  // Check the transaxial FWHM
145  if (m_transFWHM<0.)
146  {
147  Cerr("***** iImageConvolverStationaryGaussian::CheckSpecificParameters() -> Transaxial FWHM is negative !" << endl);
148  return 1;
149  }
150  // Check the axial FWHM
151  if (m_axialFWHM<0.)
152  {
153  Cerr("***** iImageConvolverStationaryGaussian::CheckSpecificParameters() -> Axial FWHM is negative !" << endl);
154  return 1;
155  }
156  // Normal end
157  return 0;
158 }
159 
160 // =====================================================================
161 // ---------------------------------------------------------------------
162 // ---------------------------------------------------------------------
163 // =====================================================================
164 
166 {
167  // Verbose
168  if (m_verbose>=2)
169  {
170  Cout("iImageConvolverStationaryGaussian::BuildConvolutionKernel() -> Compute convolution kernel" << endl);
171  Cout(" --> Transaxial FWHM: " << m_transFWHM << " mm" << endl);
172  Cout(" --> Axial FWHM: " << m_axialFWHM << " mm" << endl);
173  Cout(" --> Included sigmas: " << m_nbSigmas << endl);
174  }
175 
176  // --------------------------------
177  // Stationary kernel
178  // --------------------------------
179 
180  // Set the number of kernels to 1
181  m_nbKernels = 1;
182  // Specify that it is stationary (we do not overload the Convolve() function here)
183  m_stationary = true;
184  // Compute kernel half size based on the number of sigmas included in the kernel
185  FLTNB two_sqrt_two_ln_two = 2.*sqrt(2.*log(2.));
186  FLTNB half_kern_floatX = m_transFWHM * m_nbSigmas / (two_sqrt_two_ln_two*mp_ImageDimensionsAndQuantification->GetVoxSizeX());
187  FLTNB half_kern_floatY = m_transFWHM * m_nbSigmas / (two_sqrt_two_ln_two*mp_ImageDimensionsAndQuantification->GetVoxSizeY());
188  FLTNB half_kern_floatZ = m_axialFWHM * m_nbSigmas / (two_sqrt_two_ln_two*mp_ImageDimensionsAndQuantification->GetVoxSizeZ());
189  // For each half size, if the decimal part is below 0.5, then it will be taken into account into the central kernel's voxel,
190  // otherwise, we had a voxel to the half size
191  INTNB half_kern_dimX = ((INTNB)(half_kern_floatX));
192  if (half_kern_floatX-((FLTNB)half_kern_dimX)>0.5) half_kern_dimX++;
193  INTNB half_kern_dimY = ((INTNB)(half_kern_floatY));
194  if (half_kern_floatY-((FLTNB)half_kern_dimY)>0.5) half_kern_dimY++;
195  INTNB half_kern_dimZ = ((INTNB)(half_kern_floatZ));
196  if (half_kern_floatZ-((FLTNB)half_kern_dimZ)>0.5) half_kern_dimZ++;
197  // Allocate kernel dimensions
198  mp_dimKernelX = (INTNB*)malloc(1*sizeof(INTNB));
199  mp_dimKernelY = (INTNB*)malloc(1*sizeof(INTNB));
200  mp_dimKernelZ = (INTNB*)malloc(1*sizeof(INTNB));
201  // Compute kernel dimensions
202  mp_dimKernelX[0] = half_kern_dimX * 2 + 1;
203  mp_dimKernelY[0] = half_kern_dimY * 2 + 1;
204  mp_dimKernelZ[0] = half_kern_dimZ * 2 + 1;
207  // Verbose
208  if (m_verbose>=2) Cout (" --> Kernel dimensions: [" << mp_dimKernelX[0] << ";" << mp_dimKernelY[0] << ";" << mp_dimKernelZ[0] << "]" << endl);
209  // Allocate kernel
210  m2p_kernel = (FLTNB**)malloc(1*sizeof(FLTNB*));
211  m2p_kernel[0] = (FLTNB*)malloc(m_dimKernelXYZ*sizeof(FLTNB));
212 
213  // --------------------------------
214  // Compute the kernel
215  // --------------------------------
216 
217  // Variables
218  FLTNB sigma_x = m_transFWHM / two_sqrt_two_ln_two;
219  FLTNB sigma_y = m_transFWHM / two_sqrt_two_ln_two;
220  FLTNB sigma_z = m_axialFWHM / two_sqrt_two_ln_two;
221  FLTNB mu_x = (FLTNB)(half_kern_dimX);
222  FLTNB mu_y = (FLTNB)(half_kern_dimY);
223  FLTNB mu_z = (FLTNB)(half_kern_dimZ);
224  // Compute kernel by integration of the Gaussian distribution (this method is not exactly exact... but is closer to it compared to single Gaussian points)
225  FLTNB sum_kernel = 0.;
226  for (int z=0, index=0; z<mp_dimKernelZ[0]; z++)
227  {
228  FLTNB kern_Z = 0.5 * ( erf( ( (((FLTNB)z)-mu_z+0.5)*mp_ImageDimensionsAndQuantification->GetVoxSizeZ()) / (sqrt(2.)*sigma_z) )
229  - erf( ( (((FLTNB)z)-mu_z-0.5)*mp_ImageDimensionsAndQuantification->GetVoxSizeZ()) / (sqrt(2.)*sigma_z) ) );
230  for (int y=0; y<mp_dimKernelY[0]; y++)
231  {
232  FLTNB kern_ZY = 0.5 * ( erf( ( (((FLTNB)y)-mu_y+0.5)*mp_ImageDimensionsAndQuantification->GetVoxSizeY()) / (sqrt(2.)*sigma_y) )
233  - erf( ( (((FLTNB)y)-mu_y-0.5)*mp_ImageDimensionsAndQuantification->GetVoxSizeY()) / (sqrt(2.)*sigma_y) ) )
234  * kern_Z;
235  for (int x=0; x<mp_dimKernelX[0]; x++, index++)
236  {
237  m2p_kernel[0][index] = 0.5 * ( erf( ( (((FLTNB)x)-mu_x+0.5)*mp_ImageDimensionsAndQuantification->GetVoxSizeX()) / (sqrt(2.)*sigma_x) )
238  - erf( ( (((FLTNB)x)-mu_x-0.5)*mp_ImageDimensionsAndQuantification->GetVoxSizeX()) / (sqrt(2.)*sigma_x) ) )
239  * kern_ZY;
240  sum_kernel += m2p_kernel[0][index];
241  }
242  }
243  }
244  // Normalize kernel
245  for (int v=0; v<m_dimKernelXYZ; v++) m2p_kernel[0][v] /= sum_kernel;
246  // Normal end
247  return 0;
248 }
249 
250 // =====================================================================
251 // ---------------------------------------------------------------------
252 // ---------------------------------------------------------------------
253 // =====================================================================
Declaration of class oImageDimensionsAndQuantification.
#define FLTNB
Definition: gVariables.hh:81
FLTNB GetVoxSizeX()
Get the voxel's size along the X axis, in mm.
int BuildConvolutionKernel()
A private function used to build the convolution kernel specific to the child convolver.
FLTNB GetVoxSizeZ()
Get the voxel's size along the Z axis, in mm.
int ReadConfigurationFile(const string &a_fileOptions)
A function used to read options from a configuration file.
int CheckSpecificParameters()
A private function used to check the parameters settings specific to the child module.
oImageDimensionsAndQuantification * mp_ImageDimensionsAndQuantification
FLTNB GetVoxSizeY()
Get the voxel's size along the Y axis, in mm.
#define Cerr(MESSAGE)
Declaration of class iImageConvolverStationaryGaussian.
int ReadDataASCIIFile(const string &a_file, const string &a_keyword, T *ap_return, int a_nbElts, bool a_mandatoryFlag)
Look for "a_nbElts" elts in the "a_file" file matching the "a_keyword" string passed as parameter a...
Definition: gOptions.cc:123
iImageConvolverStationaryGaussian()
The constructor of iImageConvolverStationaryGaussian.
#define KEYWORD_MANDATORY
Definition: gOptions.hh:48
int ReadOptionsList(const string &a_listOptions)
A function used to read options from a list of options.
#define INTNB
Definition: gVariables.hh:92
Declaration of class oImageSpace.
~iImageConvolverStationaryGaussian()
The destructor of iImageConvolverStationaryGaussian.
#define Cout(MESSAGE)
Declaration of class vImageConvolver.
int ReadStringOption(const string &a_input, T *ap_return, int a_nbElts, const string &sep, const string &a_option)
Parse the 'a_input' string corresponding to the 'a_option' into 'a_nbElts' elements, using the 'sep' separator. The results are returned in the templated 'ap_return' dynamic templated array. Call "ConvertFromString()" to perform the correct conversion depending on the type of the data to convert.
Definition: gOptions.cc:62
void ShowHelp()
A function used to show help about the child module.
This abstract class is the generic image convolver class used by the oImageConvolverManager.