CASToR  2.0
Tomographic Reconstruction (PET/SPECT/CT)
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Macros Groups Pages
vImageConvolver.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"
34 
35 // =====================================================================
36 // ---------------------------------------------------------------------
37 // ---------------------------------------------------------------------
38 // =====================================================================
39 
41 {
42  // Standards
44  m_verbose = -1;
45  // Booleans
46  m_checked = false;
47  m_initialized = false;
48  m_stationary = false;
49  // Padded image
50  mp_paddedImage = NULL;
51  m_offsetX = -1;
52  m_offsetY = -1;
53  m_offsetZ = -1;
54  m_dimPadX = -1;
55  m_dimPadY = -1;
56  m_dimPadZ = -1;
57  m_dimPadXY = -1;
58  m_dimPadXYZ = -1;
59  // Convolution kernel
60  m_nbKernels = -1;
61  mp_dimKernelX = NULL;
62  mp_dimKernelY = NULL;
63  mp_dimKernelZ = NULL;
64  m2p_kernel = NULL;
65 }
66 
67 // =====================================================================
68 // ---------------------------------------------------------------------
69 // ---------------------------------------------------------------------
70 // =====================================================================
71 
73 {
74  // Simply free what should be
75  if (mp_paddedImage) free(mp_paddedImage);
76  if (mp_dimKernelX) free(mp_dimKernelX);
77  if (mp_dimKernelY) free(mp_dimKernelY);
78  if (mp_dimKernelZ) free(mp_dimKernelZ);
79  if (m2p_kernel)
80  {
81  for (int k=0; k<m_nbKernels; k++) if (m2p_kernel[k]) free(m2p_kernel[k]);
82  free(m2p_kernel);
83  }
84 }
85 
86 // =====================================================================
87 // ---------------------------------------------------------------------
88 // ---------------------------------------------------------------------
89 // =====================================================================
90 
91 // Check global parameters then check the specific parameters
93 {
94  // Check verbose
95  if (m_verbose<0)
96  {
97  Cerr("***** vImageConvolver::CheckParameters() -> The verbose level must be set positive !" << endl);
98  return 1;
99  }
100  // Check oImageDimensionsAndQuantification
102  {
103  Cerr("***** vImageConvolver::CheckParameters() -> The image dimensions must be set through a oImageDimensionsAndQuantification !" << endl);
104  return 1;
105  }
106  // Finally check the specific parameters from the child convolver
108  {
109  Cerr("***** vImageConvolver::CheckParameters() -> A problem occured while checking specific parameters of the child convolver !" << endl);
110  return 1;
111  }
112  // All set
113  m_checked = true;
114  // Normal end
115  return 0;
116 }
117 
118 // =====================================================================
119 // ---------------------------------------------------------------------
120 // ---------------------------------------------------------------------
121 // =====================================================================
122 
123 // Check that the kernel has been successfully built, then allocate the padded image
125 {
126  // Check if parameters have been checked
127  if (!m_checked)
128  {
129  Cerr("***** vImageConvolver::Initialize() -> Parameters have not been checked ! Please call CheckParameters() before." << endl);
130  return 1;
131  }
132 
133  // Build the convolution kernel
135  {
136  Cerr("***** vImageConvolver::Initialize() -> A problem occured while building the convolution kernel !" << endl);
137  return 1;
138  }
139 
140  // ---------------------------------------------------------------------------------------
141  // Check that all kernel dimensions are odd and find the maximum size along each dimension
142  // ---------------------------------------------------------------------------------------
143 
144  // Check that the number of kernels is not null or negative
145  if (m_nbKernels<1)
146  {
147  Cerr("***** vImageConvolver::Initialize() -> The number of kernels is negative or null !" << endl);
148  return 1;
149  }
150  // Check that the kernel dimensions have been allocated
151  if (mp_dimKernelX==NULL || mp_dimKernelY==NULL || mp_dimKernelZ==NULL)
152  {
153  Cerr("***** vImageConvolver::Initialize() -> The kernel dimensions are not allocated !" << endl);
154  return 1;
155  }
156  // Check that the kernel has been allocated
157  if (m2p_kernel==NULL)
158  {
159  Cerr("***** vImageConvolver::Initialize() -> The kernel is not allocated !" << endl);
160  return 1;
161  }
162  // In case the kernel is stationary and the user forgot to explicitely specify it, we force it here
163  if (m_nbKernels==1) m_stationary = true;
164 
165  // Maximum kernel sizes along each dimension
166  INTNB max_kern_size_Z = 0;
167  INTNB max_kern_size_Y = 0;
168  INTNB max_kern_size_X = 0;
169 
170  // Loop over all kernels
171  for (int k=0; k<m_nbKernels; k++)
172  {
173  // Check oddity along Z
174  if (mp_dimKernelZ[k]%2==0)
175  {
176  Cerr("***** vImageConvolver::Initialize() -> Dimension of " << k << "th convolution kernel along Z (" << mp_dimKernelZ[k] << " is not odd !" << endl);
177  return 1;
178  }
179  // Check oddity along Y
180  if (mp_dimKernelY[k]%2==0)
181  {
182  Cerr("***** vImageConvolver::Initialize() -> Dimension of " << k << "th convolution kernel along Y (" << mp_dimKernelY[k] << " is not odd !" << endl);
183  return 1;
184  }
185  // Check oddity along X
186  if (mp_dimKernelX[k]%2==0)
187  {
188  Cerr("***** vImageConvolver::Initialize() -> Dimension of " << k << "th convolution kernel along X (" << mp_dimKernelX[k] << " is not odd !" << endl);
189  return 1;
190  }
191  // Check maximum along Z
192  if (mp_dimKernelZ[k]>max_kern_size_Z) max_kern_size_Z = mp_dimKernelZ[k];
193  // Check maximum along Y
194  if (mp_dimKernelY[k]>max_kern_size_Y) max_kern_size_Y = mp_dimKernelY[k];
195  // Check maximum along X
196  if (mp_dimKernelX[k]>max_kern_size_X) max_kern_size_X = mp_dimKernelX[k];
197  }
198 
199  // Set the offsets
200  m_offsetZ = max_kern_size_Z / 2;
201  m_offsetY = max_kern_size_Y / 2;
202  m_offsetX = max_kern_size_X / 2;
203 
204  // In the convolution we assume an intensity conservation along the Z axis, so we check if the offset is higher or equal to the image dimension.
205  // In this case, the convolution will still work (it is implemented so), but computation time will be lost for nothing, so we throw a warning.
207  {
208  FLTNB over_computation_factor = ((FLTNB)(max_kern_size_Z)) / ((FLTNB)(mp_ImageDimensionsAndQuantification->GetNbVoxZ()*2-1));
209  Cerr("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl);
210  Cerr("!!!!! vImageConvolver::Initialize() -> The maximum half size of the convolution kernels is higher or equal to the number of slices of the image." << endl);
211  Cerr("!!!!! This won't affect the results, but a convolution running time factor of " << over_computation_factor << " will be lost for nothing." << endl);
212  Cerr("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl);
213  }
214 
215  // ---------------------------------------------------------------------------------------
216  // Build the padded image
217  // ---------------------------------------------------------------------------------------
218 
219  // Compute the padded image dimensions
220  m_dimPadZ = mp_ImageDimensionsAndQuantification->GetNbVoxZ() + max_kern_size_Z - 1;
221  m_dimPadY = mp_ImageDimensionsAndQuantification->GetNbVoxY() + max_kern_size_Y - 1;
222  m_dimPadX = mp_ImageDimensionsAndQuantification->GetNbVoxX() + max_kern_size_X - 1;
225 
226  // Allocate the padded image and fill it with 0.
227  mp_paddedImage = (FLTNB*)calloc(m_dimPadXYZ,sizeof(FLTNB));
228 
229  // All set
230  m_initialized = true;
231  // Normal end
232  return 0;
233 }
234 
235 // =====================================================================
236 // ---------------------------------------------------------------------
237 // ---------------------------------------------------------------------
238 // =====================================================================
239 
241 {
242  #ifdef CASTOR_DEBUG
243  if (!m_initialized)
244  {
245  Cerr("***** vImageConvolver::ApplyConvolution() -> Called while not initialized !" << endl);
246  Exit(EXIT_DEBUG);
247  }
248  #endif
249 
250  // Copy image into padded buffer image
251  CopyToPaddedImage(ap_image);
252  // Then convolve
253  if (Convolve(ap_image))
254  {
255  Cerr("***** vImageConvolver::ApplyConvolution() -> A problem occured while actually convolving !" << endl);
256  return 1;
257  }
258  // Normal end
259  return 0;
260 }
261 
262 // =====================================================================
263 // ---------------------------------------------------------------------
264 // ---------------------------------------------------------------------
265 // =====================================================================
266 
268 {
269  #ifdef CASTOR_DEBUG
270  if (!m_initialized)
271  {
272  Cerr("***** vImageConvolver::ApplyConvolutionTranspose() -> Called while not initialized !" << endl);
273  Exit(EXIT_DEBUG);
274  }
275  #endif
276 
277  // Copy image into padded buffer image
278  CopyToPaddedImage(ap_image);
279  // Then convolve
280  if (ConvolveTranspose(ap_image))
281  {
282  Cerr("***** vImageConvolver::ApplyConvolutionTranspose() -> A problem occured while actually convolving !" << endl);
283  return 1;
284  }
285  // Normal end
286  return 0;
287 }
288 
289 // =====================================================================
290 // ---------------------------------------------------------------------
291 // ---------------------------------------------------------------------
292 // =====================================================================
293 
295 {
296  #ifdef CASTOR_DEBUG
297  if (!m_initialized)
298  {
299  Cerr("***** vImageConvolver::CopyToPaddedImage() -> Called while not initialized !" << endl);
300  Exit(EXIT_DEBUG);
301  }
302  #endif
303 
304  // Copy the image into the padded buffer, centered with respect to the pad
305  INTNB z;
306  #pragma omp parallel for private(z) schedule(guided)
308  {
309  // Pad coordinate Z
310  INTNB padz = z + m_offsetZ;
312  {
313  // Pad coordinate Y
314  INTNB pady = y + m_offsetY;
316  {
317  // Pad coordinate X
318  INTNB padx = x + m_offsetX;
319  // Global original coordinate
321  // Global pad coordinate
322  INTNB coord_pad = padz*m_dimPadXY + pady*m_dimPadX + padx;
323  // Affect buffer
324  mp_paddedImage[coord_pad] = ap_inputImage[coord_orig];
325  // Zero image
326  ap_inputImage[coord_orig] = 0.;
327  }
328  }
329  }
330 }
331 
332 // =====================================================================
333 // ---------------------------------------------------------------------
334 // ---------------------------------------------------------------------
335 // =====================================================================
336 
337 int vImageConvolver::Convolve(FLTNB* ap_outputImage)
338 {
339  #ifdef CASTOR_DEBUG
340  if (!m_initialized)
341  {
342  Cerr("***** vImageConvolver::Convolve() -> Called while not initialized !" << endl);
343  Exit(EXIT_DEBUG);
344  }
345  #endif
346 
347  // This function is designed to work universally and exclusively on a stationary kernel.
348  // Note that we suppose that a stationary kernel is also symmetric (but not necessarily isotropic).
349  // For asymetric and non-stationary kernels, have to overload this function in the child convolver.
350  // Note also that we choose the out-to-in implementation (outer loop on the output image) so that the
351  // parallelism is operated on the output and is thus thread-safe.
352  // Finally note that transaxially, we do not implement an image intensity preservation policy as we suppose that
353  // the distance between the object and the image border is taller than the kernel half-width (it should be so).
354  if (m_stationary)
355  {
356  INTNB stationary_kernel = 0;
357  // Convolve (OpenMP on Z dimension)
358  INTNB z;
359  #pragma omp parallel for private(z) schedule(guided)
361  {
362  // As a first step, we virtually go through the kernel along Z to compute the actual kernel integral that will
363  // be really used in the convolution. For each virtual slice contributing to the convolved slice, we test if it
364  // is inside or outside the image; in the former case we count it into the integral and in the latter we don't.
365  FLTNB kernel_integral = 0.;
366  for (INTNB zk=0; zk<mp_dimKernelZ[stationary_kernel]; zk++)
367  {
368  // Index of the virtual slice which will (or not) contribute to the convolved slice
369  INTNB virtualZ = z + zk - m_offsetZ;
370  // Test if it is not in the image, then we move onto the next virtual slice
371  if (virtualZ<0 || virtualZ>=mp_ImageDimensionsAndQuantification->GetNbVoxZ()) continue;
372  // Otherwise, we update the kernel integral
373  INTNB kernZ_base = zk * mp_dimKernelX[stationary_kernel] * mp_dimKernelY[stationary_kernel];
374  for (INTNB xyk=0; xyk<mp_dimKernelX[stationary_kernel]*mp_dimKernelY[stationary_kernel]; xyk++)
375  kernel_integral += m2p_kernel[stationary_kernel][kernZ_base+xyk];
376  }
377  // Some precomputation for maximum efficiency
380  {
381  // Some precomputation for maximum efficiency
382  INTNB indexZY_base = indexZ_base + y*mp_ImageDimensionsAndQuantification->GetNbVoxX();
384  {
385  // Output image index
386  INTNB index_img = indexZY_base + x;
387  // Inner loops on convolution kernel dimensions
388  for (INTNB zk=0, index_kern=0; zk<mp_dimKernelZ[stationary_kernel]; zk++)
389  {
390  // Some precomputation for maximum efficiency
391  INTNB indexZ_pad = (z + zk) * m_dimPadXY;
392  for (INTNB yk=0; yk<mp_dimKernelY[stationary_kernel]; yk++)
393  {
394  // Some precomputation for maximum efficiency
395  INTNB indexZY_pad = indexZ_pad + (y + yk) * m_dimPadX;
396  for (INTNB xk=0; xk<mp_dimKernelX[stationary_kernel]; xk++, index_kern++)
397  {
398  // Padded image index
399  INTNB index_pad = indexZY_pad + x + xk;
400  // Apply contribution
401  ap_outputImage[index_img] += mp_paddedImage[index_pad] * m2p_kernel[stationary_kernel][index_kern] / kernel_integral;
402  }
403  }
404  }
405  }
406  }
407  }
408  // Normal end
409  return 0;
410  }
411  // Otherwise the convolve function has to be implemented by the child non-stationary convolver
412  else
413  {
414  Cerr("***** vImageConvolver::Convolve() -> To use a non-stationary kernel, you should overload this function to implement your own in your child convolver !" << endl);
415  return 1;
416  }
417 }
418 
419 // =====================================================================
420 // ---------------------------------------------------------------------
421 // ---------------------------------------------------------------------
422 // =====================================================================
423 
425 {
426  #ifdef CASTOR_DEBUG
427  if (!m_initialized)
428  {
429  Cerr("***** vImageConvolver::ConvolveTranspose() -> Called while not initialized !" << endl);
430  Exit(EXIT_DEBUG);
431  }
432  #endif
433 
434  // This function is designed to work universally on a stationary kernel.
435  // Note that we suppose a stationary kernel is also symmetric (otherwise it is absurd)
436  if (m_stationary)
437  {
438  return Convolve(ap_outputImage);
439  }
440  // Otherwise the convolve function has to be implemented by the child non-stationary convolver
441  else
442  {
443  Cerr("***** vImageConvolver::ConvolveTranspose() -> To use a non-stationary kernel, you should overload this function to implement your own in your child convolver !" << endl);
444  return 1;
445  }
446 }
447 
448 // =====================================================================
449 // ---------------------------------------------------------------------
450 // ---------------------------------------------------------------------
451 // =====================================================================
Declaration of class oImageDimensionsAndQuantification.
#define FLTNB
Definition: gVariables.hh:81
virtual int ConvolveTranspose(FLTNB *ap_outputImage)
A private function used to apply the transpose convolution on the padded image to the provided output...
int Initialize()
A public function used to initialize the module.
virtual int BuildConvolutionKernel()=0
A private function used to build the convolution kernel specific to the child convolver.
void Exit(int code)
oImageDimensionsAndQuantification * mp_ImageDimensionsAndQuantification
#define Cerr(MESSAGE)
int CheckParameters()
A public function used to check the parameters settings.
INTNB GetNbVoxXY()
Get the number of voxels in a slice.
int ApplyConvolutionTranspose(FLTNB *ap_image)
A public function used to apply the transpose convolution module on the provided image.
virtual ~vImageConvolver()
The destructor of vImageConvolver.
vImageConvolver()
The constructor of vImageConvolver.
#define INTNB
Definition: gVariables.hh:92
int ApplyConvolution(FLTNB *ap_image)
A public function used to apply the convolution module on the provided image.
Declaration of class oImageSpace.
virtual int CheckSpecificParameters()=0
A private function used to check the parameters settings specific to the child convolver.
#define EXIT_DEBUG
Definition: gVariables.hh:97
virtual int Convolve(FLTNB *ap_outputImage)
A private function used to apply the convolution on the padded image to the provided output image...
INTNB GetNbVoxX()
Get the number of voxels along the X axis.
INTNB GetNbVoxZ()
Get the number of voxels along the Z axis.
Declaration of class vImageConvolver.
void CopyToPaddedImage(FLTNB *ap_inputImage)
A private function used to copy the provided image into the padded buffer.
INTNB GetNbVoxY()
Get the number of voxels along the Y axis.