CASToR  2.0
Tomographic Reconstruction (PET/SPECT/CT)
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Macros Groups Pages
vDynamicModel.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 "vDynamicModel.hh"
32 
33 // =====================================================================
34 // ---------------------------------------------------------------------
35 // ---------------------------------------------------------------------
36 // =====================================================================
37 /*
38  \fn vDynamicModel
39  \brief Constructor of vDynamicModel. Simply set all data members to default values.
40 */
42 {
43  mp_ID = NULL;
44  m_nbTimeBF = -1;
45  m_nbModelParam = -1;
46  m_verbose = -1;
47  m_checked = false;
48  m_initialized = false;
49  m2p_parametricImages = NULL;
50  m2p_modelTACs = NULL;
51  m2p_outputParImages = NULL;
52  m_saveParImageFlag = true;
53 }
54 
55 
56 
57 
58 // =====================================================================
59 // ---------------------------------------------------------------------
60 // ---------------------------------------------------------------------
61 // =====================================================================
67 
68 
69 
70 
71 // =====================================================================
72 // ---------------------------------------------------------------------
73 // ---------------------------------------------------------------------
74 // =====================================================================
82 {
83  if(m_verbose>=2) Cout("vDynamicModel::CheckParameters ..."<< endl);
84 
85  // Check image dimensions
86  if (mp_ID==NULL)
87  {
88  Cerr("***** vDynamicModel::CheckParameters() -> No image dimensions provided !" << endl);
89  return 1;
90  }
91 
92  // Check verbosity
93  if (m_verbose<0)
94  {
95  Cerr("***** vDynamicModel::CheckParameters() -> Wrong verbosity level provided !" << endl);
96  return 1;
97  }
98 
99  // Check number of basis functions
100  if (m_nbTimeBF <0)
101  {
102  Cerr("***** vDynamicModel::CheckParameters() -> Basis functions number has not been initialized !" << endl);
103  return 1;
104  }
105 
106  // Check number of parameters in the model
107  if (m_nbModelParam <0)
108  {
109  Cerr("***** vDynamicModel::CheckParameters() -> Basis functions number has not been initialized !" << endl);
110  return 1;
111  }
112 
113  // Check parameters of the child class (if this function is overloaded)
115  {
116  Cerr("***** vDynamicModel::CheckParameters() -> An error occurred while checking parameters of the child dynamic class !" << endl);
117  return 1;
118  }
119 
120  // Normal end
121  m_checked = true;
122  return 0;
123 }
124 
125 
126 
127 
128 // =====================================================================
129 // ---------------------------------------------------------------------
130 // ---------------------------------------------------------------------
131 // =====================================================================
137 {
138  if(m_verbose >=2) Cout("vDynamicModel::ComputeOutputParImage ..." <<endl);
139 
140  // If we save parametric image
142  for(int p=0 ; p<m_nbModelParam ; p++)
143  {
144  // If output image matrix is allocated, then copy the current parametric images
145  if(m2p_outputParImages != NULL
147  {
148  for(int v=0 ; v<mp_ID->GetNbVoxXYZ() ; v++)
150  }
151  // Point to parametric images otherwise
152  else
154  }
155 
156 }
157 
158 
159 
160 
161 // =====================================================================
162 // ---------------------------------------------------------------------
163 // ---------------------------------------------------------------------
164 // =====================================================================
165 /*
166  \fn vDynamicModel::ApplyOutputFOVMaskingOnParametricImages
167  \brief Mask the outside of the transaxial FOV based on the m_fovOutPercent
168  \details Similar to the eponym function in ImageSpace, but on parametric images
169 */
171 {
172  // If the output FOV percent is under 0 (the default value) and number of axial slices to be removed is 0, we do not mask anything
173  if ( mp_ID->GetFOVOutPercent()<=0. &&
174  mp_ID->GetNbSliceOutMask()==0 ) return 0;
175 
176  // Verbose
177  if (m_verbose>=2)
178  {
179  Cout("vDynamicModel::ApplyOutputFOVMaskingOnParametricImages() -> Mask output image" << endl);
180  if (mp_ID->GetFOVOutPercent()>0.) Cout(" --> Mask transaxial FOV outside " << mp_ID->GetFOVOutPercent() << " %" << endl);
181  if (mp_ID->GetNbSliceOutMask()>0) Cout(" --> Mask " << mp_ID->GetNbSliceOutMask() << " slices from both axial FOV limits" << endl);
182  }
183  // -----------------------------------------------
184  // Transaxial FOV masking
185  // -----------------------------------------------
186  if (mp_ID->GetFOVOutPercent()>0.)
187  {
188  // Precast half the number of voxels over X and Y minus 1 (for efficiency)
189  FLTNB flt_base_x = 0.5*((FLTNB)(mp_ID->GetNbVoxX()-1));
190  FLTNB flt_base_y = 0.5*((FLTNB)(mp_ID->GetNbVoxY()-1));
191 
192  // Compute FOV elipse radius over X and Y, then squared
193  FLTNB squared_radius_x = 0.5 * ((FLTNB)(mp_ID->GetNbVoxX())) * mp_ID->GetVoxSizeX()
194  * mp_ID->GetFOVOutPercent() / 100.;
195  squared_radius_x *= squared_radius_x;
196  FLTNB squared_radius_y = 0.5 * ((FLTNB)(mp_ID->GetNbVoxY())) * mp_ID->GetVoxSizeY()
197  * mp_ID->GetFOVOutPercent() / 100.;
198  squared_radius_y *= squared_radius_y;
199 
200  // We assume that the computation of the distance from the center for a given
201  // voxel and comparing it with the output FOV percentage costs more than performing
202  // the loops in an inverse order compared to how the image is stored in memory.
203  // Thus we begin the loops over X, then Y, then we test and if test passes, we
204  // do the remaining loop over Z and over all dynamic dimensions.
205  int x;
206  #pragma omp parallel for private(x) schedule(guided)
207  for (x=0; x<mp_ID->GetNbVoxX(); x++)
208  {
209  // Compute X distance from image center, then squared
210  FLTNB squared_distance_x = (((FLTNB)x)-flt_base_x) * mp_ID->GetVoxSizeX();
211  squared_distance_x *= squared_distance_x;
212  // Loop over Y
213  for (int y=0; y<mp_ID->GetNbVoxY(); y++)
214  {
215  // Compute Y distance from image center, then squared
216  FLTNB squared_distance_y = (((FLTNB)y)-flt_base_y) * mp_ID->GetVoxSizeY();
217  squared_distance_y *= squared_distance_y;
218  // Test if the voxel is inside the FOV elipse, then we skip this voxel
219  if ( squared_distance_x/squared_radius_x + squared_distance_y/squared_radius_y <= 1. ) continue;
220  // Loop over Z
221  for (int z=0; z<mp_ID->GetNbVoxZ(); z++)
222  {
223  // Compute global voxel index
224  INTNB index = z*mp_ID->GetNbVoxXY() + y*mp_ID->GetNbVoxX() + x;
225 
226  for (int b=0 ; b<m_nbTimeBF ; b++)
227  m2p_outputParImages[b][index] = 0.;
228 
229  }
230  }
231  }
232  }
233 
234  // -----------------------------------------------
235  // Axial FOV masking
236  // -----------------------------------------------
237  if (mp_ID->GetNbSliceOutMask()>0)
238  {
239  INTNB removed_slices = mp_ID->GetNbSliceOutMask();
240 
241  // Mask slices
242  for (int b=0 ; b<m_nbTimeBF ; b++)
243  for (int z=0; z<removed_slices; z++)
244  {
245  // First slices
246  INTNB base_z_first = z*mp_ID->GetNbVoxXY();
247  // Loop over Y and X
248  for (int i=0; i<mp_ID->GetNbVoxXY(); i++)
249  {
250  INTNB index = base_z_first + i;
251  m2p_outputParImages[b][index] = 0.;
252  }
253  // Last slices
254  INTNB base_z_last = (mp_ID->GetNbVoxZ()-1-z)*mp_ID->GetNbVoxXY();
255  // Loop over Y and X
256  for (int i=0; i<mp_ID->GetNbVoxXY(); i++)
257  {
258  INTNB index = base_z_last + i;
259  m2p_outputParImages[b][index] = 0.;
260  }
261  }
262  }
263 
264  // End
265  return 0;
266 }
267 
268 
269 
270 
271 // =====================================================================
272 // ---------------------------------------------------------------------
273 // ---------------------------------------------------------------------
274 // =====================================================================
275 /*
276  \fn SaveParametricImages
277  \param a_iteration : current iteration index
278  \param a_subset : current number of subsets (or -1 by default)
279  \brief This function is pure virtual so must be implemented by children \n
280  Call SaveParametricImages() function of the dynamic model object,
281  in order to write on disk any parametric image of the model
282  \return 0 if success, positive value otherwise
283 */
284 int vDynamicModel::SaveParametricImages(int a_iteration, int a_subset)
285 {
286  if(m_verbose >=2) Cout("vDynamicModel::SaveParametricImages ..." <<endl);
287 
288 
290  {
291  // Get the output manager
292  sOutputManager* p_output_manager = sOutputManager::GetInstance();
293 
294  // Interfile
295  string path_to_image = p_output_manager->GetPathName() + p_output_manager->GetBaseName();
296 
297  // Add a suffix for iteration
298  if (a_iteration >= 0)
299  {
300  stringstream ss; ss << a_iteration + 1;
301  path_to_image.append("par_it").append(ss.str());
302  }
303 
304  // Add a suffix for subset (if not negative by default), this means that we save a 'subset' image
305  if (a_subset >= 0)
306  {
307  stringstream ss; ss << a_subset + 1;
308  path_to_image.append("_ss").append(ss.str());
309  }
310 
311  // Write interfile parametric image
312  if(IntfWriteImgDynCoeffFile(path_to_image,
314  mp_ID,
316  m_verbose) )
317  {
318  Cerr("***** iPatlakModel::SaveParametricImages()-> Error writing Interfile of output image !" << endl);
319  return 1;
320  }
321 
322  }
323 
324  return 0;
325 }
virtual int CheckSpecificParameters()=0
This function is used to check the parameters of the child functions before initialization if require...
virtual int CheckParameters()
This function is used to check parameters after the latter have been all set using Set functions...
FLTNB ** m2p_modelTACs
#define FLTNB
Definition: gVariables.hh:81
FLTNB GetVoxSizeX()
Get the voxel's size along the X axis, in mm.
FLTNB GetFOVOutPercent()
Get the percentage of transaxial unmasked FOV.
oImageDimensionsAndQuantification * mp_ID
virtual int ApplyOutputFOVMaskingOnParametricImages()
Mask the outside of the transaxial FOV based on the m_fovOutPercent.
static sOutputManager * GetInstance()
Instanciate the singleton object and Initialize member variables if not already done, return a pointer to this object otherwise.
uint16_t m_nbModelParam
FLTNB GetVoxSizeY()
Get the voxel's size along the Y axis, in mm.
#define Cerr(MESSAGE)
const string & GetPathName()
Singleton class that manages output writing on disk (images, sinograms, etc). It also manages loggi...
FLTNB ** m2p_outputParImages
INTNB GetNbVoxXY()
Get the number of voxels in a slice.
virtual ~vDynamicModel()
Destructor of vDynamicModel.
vDynamicModel()
Constructor of vDynamicModel. Simply set all data members to default values.
#define INTNB
Definition: gVariables.hh:92
const string & GetBaseName()
int IntfWriteImgDynCoeffFile(const string &a_pathToImg, FLTNB **a2p_ImgMatrix, oImageDimensionsAndQuantification *ap_ID, int a_nbParImgs, int vb)
Function dedicated to Interfile image writing for dynamic coefficients images.
FLTNB ** m2p_parametricImages
INTNB GetNbVoxXYZ()
Get the total number of voxels.
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 vDynamicModel.
#define Cout(MESSAGE)
virtual void ComputeOutputParImage()
Compute output image using the m2p_parametricImages matrix Store the result in the m2p_outputParImage...
virtual int SaveParametricImages(int a_iteration, int a_subset=-1)
This function is pure virtual so must be implemented by children Call SaveParametricImages() functi...
INTNB GetNbSliceOutMask()
Get the number of extrem slices that will be masked at each side.
INTNB GetNbVoxY()
Get the number of voxels along the Y axis.