CASToR  3.0
Tomographic Reconstruction (PET/SPECT/CT)
iLinearModel.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-2019 all CASToR contributors listed below:
18 
19  --> Didier BENOIT, Claude COMTAT, Marina FILIPOVIC, Thibaut MERLIN, Mael MILLARDET, Simon STUTE, Valentin VIELZEUF
20 
21 This is CASToR version 3.0.
22 */
23 
31 #include "iLinearModel.hh"
32 
33 // =====================================================================
34 // ---------------------------------------------------------------------
35 // ---------------------------------------------------------------------
36 // =====================================================================
37 /*
38  \fn iLinearModel
39  \brief Constructor of iLinearModel. Simply set all data members to default values.
40 */
41 iLinearModel::iLinearModel() : vDynamicModel()
42 {
43  m_nbRgateBF = -1;
44  m_nbCgateBF = -1;
45  m_nbRGModelParam = -1;
46  m_nbCGModelParam = -1;
47 
48  m_OptimisationMethod = -1;
49  m_nbWeightFactors=-1;
50 
51  m2p_respBasisFunctions = NULL;
52  m2p_cardBasisFunctions = NULL;
53 
54  m_fileOptions = "";
55  m_listOptions = "";
56 
57  mp_corrBasisCoeffs = NULL;
58  mp_corrBasisFunctions = NULL;
59  m_nbLinearModelCycles = 1;
60  m_basisFunctionsUpdStartIte =-1;
61  m_basisFunctionsUpdRatio = 1;
62  m_basisFunctionsUpdIdx = 0;
63 }
64 
65 
66 
67 
68 // =====================================================================
69 // ---------------------------------------------------------------------
70 // ---------------------------------------------------------------------
71 // =====================================================================
72 /*
73  \fn ~iLinearModel
74  \brief Destructor of iLinearModel
75 */
76 iLinearModel::~iLinearModel()
77 {
78  if(m_initialized)
79  {
80  for(int rb=0 ; rb<m_nbRgateBF ; rb++)
81  {
82  if (m2p_respBasisFunctions[rb]) delete[] m2p_respBasisFunctions[rb];
83  if (m2p_RGParametricImages[rb]) delete[] m2p_RGParametricImages[rb];
84  }
85 
86  for(int cb=0 ; cb<m_nbCgateBF ; cb++)
87  {
88  if (m2p_cardBasisFunctions[cb]) delete[] m2p_cardBasisFunctions[cb];
89  if (m2p_CGParametricImages[cb]) delete[] m2p_CGParametricImages[cb];
90  }
91 
92  if (mp_corrBasisCoeffs != NULL) delete[] mp_corrBasisCoeffs;
93  if (mp_corrBasisFunctions != NULL) delete[] mp_corrBasisFunctions;
94 
95 
96  // Free NNLS variables
97  if(m_OptimisationMethod == OPTIMISATION_METHOD_NNLS)
98  {
99  for( int th=0 ; th<mp_ID->GetNbThreadsForImageComputation() ; th++ )
100  {
101  for(int n=0 ; n<m_nnlsN ; n++)
102  if ( m3p_nnlsA[ th ][ n ]) delete[] m3p_nnlsA[ th ][ n ];
103 
104  // Init 2D coefficient matrix for NNLS estimation
105  if( m3p_nnlsA[ th ] ) delete[] m3p_nnlsA[ th ];
106 
107  // Init solution vector and working matrix for NNLS estimation
108 
109  if( m2p_nnlsB[ th ] ) delete[] m2p_nnlsB[ th ];
110  if( m2p_nnlsMat[ th ] ) delete[] m2p_nnlsMat[ th ];
111  if( m2p_nnlsX[ th ] ) delete[] m2p_nnlsX[ th ];
112  if( m2p_nnlsWp[ th ] ) delete[] m2p_nnlsWp[ th ];
113  if( m2p_nnlsIdx[ th ] ) delete[] m2p_nnlsIdx[ th ];
114  }
115 
116  if( m3p_nnlsA ) delete[] m3p_nnlsA;
117  if( m2p_nnlsB ) delete[] m2p_nnlsB;
118  if( m2p_nnlsMat ) delete[] m2p_nnlsMat;
119  if( m2p_nnlsX ) delete[] m2p_nnlsX;
120  if( m2p_nnlsWp ) delete[] m2p_nnlsWp;
121  if( m2p_nnlsIdx ) delete[] m2p_nnlsIdx;
122  if( mp_w ) delete[] mp_w;
123  }
124  }
125 }
126 
127 
128 
129 
130 // =====================================================================
131 // ---------------------------------------------------------------------
132 // ---------------------------------------------------------------------
133 // =====================================================================
134 /*
135  \fn ShowHelp
136  \brief Print out specific help about the implementation of this model
137  and its initialization
138 */
139 void iLinearModel::ShowHelp()
140 {
141  cout << "-- This class implements a general linear dynamic model applied between the images of a dynamic acquisition." << endl;
142  cout << "-- The model is applied on a voxel-by-voxel basis between the images of the frames and/or respiratory/cardiac gates. " << endl;
143  cout << "-- The main keywords 'DYNAMIC FRAMING', 'RESPIRATORY GATING' and 'CARDIAC GATING' ahead of the parameters allow to define at which level the model parameters must be applied. " << endl;
144  cout << "-- Main parameters to define are:" << endl;
145  cout << " -> The number of basis functions / parametric images defined in the model" << endl;
146  cout << " -> Basis function initial values" << endl;
147  cout << " -> Parametric images initialization (optional)" << endl;
148  cout << " -> Optimisation_method : x (mandatory) optimization method for voxelwise parameters estimation." << endl;
149  cout << endl;
150  cout << " It can be initialized using a configuration text file with the following keywords and information :" << endl;
151  cout << " - Mandatory keywords :" << endl;
152  cout << " The following keywords are mandatory for at least one dynamic image level (Dynamic frame, Respiratory gating or Cardiac gating) :" << endl;
153  cout << " 'Number_basis_functions:' Enter the number of basis function for each image of time frame or respiratory/cardiac gate " << endl;
154  cout << " 'Basis_functions:' Enter the basis function (bf) coefficients for each image (im) of time frame or respiratory/cardiac gate " << endl;
155  cout << " on successive lines, separated by ',' :" << endl;
156  cout << " -> basis_functions: " << endl;
157  cout << " -> coeff_bf1_im1,coeff_bf1_im2,...,coeff_bf1_imn" << endl;
158  cout << " -> coeff_bf2_im1,coeff_bf2_im2,...,coeff_bf2_imn" << endl;
159  cout << " -> etc..." << endl;
160  cout << " Optimisation_method : x (mandatory) optimization method available options: " << endl;
161  cout << " x=0: Direct ( Implementation of basis functions side by system matrix in each tomographic iterative loop " << endl;
162  cout << " (! Currently not compatible with motion correction)" << endl;
163  cout << " x=1: Nested EM " << endl;
164  cout << " x=2: Iterative non-negative Least-Square " << endl;
165  cout << " (C.L. Lawson and R.J. Hanson, Solving Least Squares Problems)" << endl;
166  cout << endl;
167  cout << " - Optional keywords :" << endl;
168  cout << " 'Parametric_image_init: image_file' Set an image file to initialize the parametric images of each dynamic model. Default initialization: '1.0 for each voxel." << endl;
169  cout << " " << endl;
170  cout << " 'Number_weight_values:' Enter the number of weights to be applied for each dynamic frame for performing WLS optimisation (Optimisation method=2) " << endl;
171  cout << " 'Weight_values:' Enter the weight values to be applied for each dynamic frame ( within DYNAMIC FRAMING/ENDDF) " << endl;
172  cout << " on one single line, separated by ',' " << endl;
173  cout << " The previous parameters must be declared inside the couple of the following specific tags: " << endl;
174  cout << " - DYNAMIC FRAMING/ENDDF for chronological frame model (dynamic model applied to chronological frames of a dynamic aquisition)" << endl;
175  cout << " - RESPIRATORY GATING/ENDRG for respiratory model (dynamic model applied to respiratory gates of a dynamic aquisition)" << endl;
176  cout << " - CARDIAC GATING/ENDCG for cardiac model (dynamic model applied to cardiac gates of a dynamic aquisition)" << endl;
177  cout << " Different levels of dynamic model can be enabled simultaneously (i.e a dynamic frame model can be used simultaneously with respiratory and/or cardiac gating model) " << endl;
178  cout << " " << endl;
179  cout << " " << endl;
180  cout << " The following keywords are optional and common to each dynamic model (Dynamic frame, Respiratory gating or Cardiac gating) :" << endl;
181  cout << " 'Number_model_iterations: x' Number of iterations of the model parameters and basis functions updates in one cycle of Nested EM" << endl;
182  cout << " (Default x == 1) (one cycle consists in x iterations in which either the parametric images or the basis functions are updated) " << endl;
183  cout << " The ratio of parametric images / basis functions updates depends on the following parameter:" << endl;
184  cout << " 'Basis_function_update_ratio: x' Ratio of update between parametric images and basis functions updates cycle " << endl;
185  cout << " (Default x == 0) Cycles consist in x iterations of the parametric images, following by x iterations of the basis functions" << endl;
186  cout << " If x == 0, only the parametric images are updated" << endl;
187  cout << " 'Basis_function_start_ite : x' Starting iteration for the update of basis functions " << endl;
188  cout << " (Default x == -1) If negative, no update of the basis functions is performed (only parametric images are updated) " << endl;
189  cout << " " << endl;
190  cout << " " << endl;
191  cout << " ---------------------------" << endl;
192  cout << " Example of initialization: " << endl;
193  cout << " DYNAMIC FRAMING" << endl;
194  cout << " Number_basis_functions : 2" << endl;
195  cout << " Basis_functions :" << endl;
196  cout << " 23682.79, 25228.74, 26636.99, 27923.61, 29101.16" << endl;
197  cout << " 5.4, 4.91, 4.48, 4.1, 3.75" << endl;
198  cout << " ENDDF" << endl;
199  cout << " " << endl;
200  cout << " RESPIRATORY GATING" << endl;
201  cout << " Number_basis_functions : 6" << endl;
202  cout << " Basis_functions :" << endl;
203  cout << " 1, 0.8, 0.6, 0.4, 0.2, 0.01" << endl;
204  cout << " 0.7, 0.9, 0.7, 0.5, 0.3, 0.1" << endl;
205  cout << " 0.4, 0.6, 0.8, 0.6, 0.4, 0.2" << endl;
206  cout << " 0.2, 0.4, 0.6, 0.8, 0.6, 0.4" << endl;
207  cout << " 0.1, 0.3, 0.5, 0.7, 0.9, 0.7" << endl;
208  cout << " 0.01, 0.2, 0.4, 0.6, 0.8, 1" << endl;
209  cout << " ENDRG" << endl;
210  cout << " " << endl;
211  cout << " Number_model_iterations : 1" << endl;
212  cout << " Basis_function_start_ite : -1" << endl;
213  cout << " Basis_function_update_ratio : 0" << endl;
214  cout << " ---------------------------" << endl;
215  cout << " " << endl;
216  cout << " The following keywords are common to all dynamic models :" << endl;
217  cout << endl;
218  cout << " 'AIC_input_file: path/to/file' : This option will enable the Creation of Patlak Basis functions for direct Patlak Reconstruction ( For nested recon look in -help-dynamic-model ) . " << endl;
219  cout << " Provide text file with the Arterial Input Curve data points and time points, in two different horizontal lines starting with " << endl;
220  cout << " 'AIC_time_points:'" << endl;
221  cout << " and " << endl;
222  cout << " 'AIC_data_points:' " << endl;
223  cout << " to indicate which dataset corresponds to each line. Values must be separated by commas." << endl;
224  cout << " Also provide a value of the total number of data points, on a new line starting with " << endl;
225  cout << " 'AIC_number_of_points:'" << endl;
226  cout << " 'Number of iterations before image update: x' Set a number 'x' of iteration to reach before using the model to generate the images at each frames/gates" << endl;
227  cout << " (Default x == 0) " << endl;
228  cout << " 'No image update: x' If set to 1, the reconstructed images for the next iteration/subset are not reestimated using the model" << endl;
229  cout << " (Default x == 0) (the code just performs standard independent reconstruction of each frames/gates) " << endl;
230  cout << " 'No parameters update: x' If set to 1, the parameters / functions of the model are not estimated with the image" << endl;
231  cout << " (Default x == 0) (this could be used to test The EstimateImageWithModel() function with specific user-provided parametric images) " << endl;
232  cout << " 'Save parametric images : x' Enable (1)/Disable(0) saving parametric images on disk" << endl;
233  cout << " (Default x == 1) " << endl;
234  cout << " 'Save blacklisted voxels images : x' Enable (1)/Disable(0) saving blacklisted voxels images on disk" << endl;
235  cout << " (Default x == 0) " << endl;
236  cout << " " << endl;
237  cout << " " << endl;
238 }
239 
240 
241 
242 // =====================================================================
243 // ---------------------------------------------------------------------
244 // ---------------------------------------------------------------------
245 // =====================================================================
246 /*
247  \fn ReadAndCheckConfigurationFileSpecific
248  \brief This function is used to read options from a configuration file when the generic Linear Model is requested.
249  \return 0 if success, other value otherwise.
250 */
251 
252 int iLinearModel::ReadAndCheckConfigurationFileSpecific()
253 {
254  if(m_verbose >=3) Cout("iLinearModel::ReadAndCheckConfigurationFileSpecific ..."<< endl);
255 
256  // Apply the Generic linear Checks for all Linear Models
257  if( ReadAndCheckConfigurationFileSpecificToAllLinearModels()==1)
258  {
259  Cerr("***** iLinearModel::ReadAndCheckConfigurationFileSpecific -> Error while trying to read and check specific configuration " << m_fileOptions << endl);
260  return 1;
261  }
262 
263  // Normal End
264  return 0 ;
265 }
266 
267 // =====================================================================
268 // ---------------------------------------------------------------------
269 // ---------------------------------------------------------------------
270 // =====================================================================
271 /*
272  \fn ReadAndCheckConfigurationFileSpecificToAllLinearModels
273  \brief This function is used to read options from a configuration file that are generic to all linear dynamic models.
274  \return 0 if success, other value otherwise.
275 */
276 int iLinearModel::ReadAndCheckConfigurationFileSpecificToAllLinearModels()
277 {
278  if(m_verbose >=3) Cout("iLinearModel::ReadAndCheckConfigurationFileSpecificToAllLinearModels ..."<< endl);
279 
280  // The file will be fully processed in the Initialize() function
281  ifstream in_file(m_fileOptions.c_str(), ios::in);
282 
283  if(in_file)
284  {
285  // Check first the file contains the mandatory keyword(s)
286  bool file_is_good = false;
287 
288  string line="", kword_to_search="";
289  while(!in_file.eof())
290  {
291  getline(in_file, line);
292 
293  //remove comment
294  if (line.find("#") != string::npos) line = line.substr(0, line.find_first_of("#")) ;
295 
296  if (line.find("DYNAMIC FRAMING") != string::npos)
297  {
298  if(kword_to_search != "")
299  {
300  Cerr("***** iLinearModel::ReadAndCheckConfigurationFileSpecificToAllLinearModels -> Error, found an end tag 'END**' before 'DYNAMIC FRAMING' in configuration file: " << m_fileOptions << endl);
301  return 1;
302  }
303  else
304  kword_to_search = "ENDDF";
305  }
306 
307  else if (line.find("RESPIRATORY GATING") != string::npos)
308  {
309  if(kword_to_search != "")
310  {
311  Cerr("***** iLinearModel::ReadAndCheckConfigurationFileSpecificToAllLinearModels -> Error, found an end tag 'END**' before 'RESPIRATORY GATING' in configuration file: " << m_fileOptions << endl);
312  return 1;
313  }
314  else
315  kword_to_search = "ENDRG";
316  }
317 
318  else if (line.find("CARDIAC GATING") != string::npos)
319  {
320  if(kword_to_search != "")
321  {
322  Cerr("***** iLinearModel::ReadAndCheckConfigurationFileSpecificToAllLinearModels -> Error, found an end tag 'END**' before 'CARDIAC GATING' in configuration file: " << m_fileOptions << endl);
323  return 1;
324  }
325  else
326  kword_to_search = "ENDCG";
327  }
328 
329  else if (kword_to_search != ""
330  && line.find(kword_to_search) != string::npos)
331  {
332  kword_to_search = "";
333  file_is_good = true;
334  }
335  }
336 
337  if(!file_is_good)
338  {
339  Cerr("***** iLinearModel::ReadAndCheckConfigurationFileSpecificToAllLinearModels -> Error, no mandatory tags ('DYNAMIC FRAMING'/'ENDDF', 'RESPIRATORY GATING'/'ENDRG', 'CARDIAC GATING'/'ENDCG' detected in configuration file: "
340  << m_fileOptions << ". Check help for more info" << endl);
341  return 1;
342  }
343 
344  if( ReadDataASCIIFile(m_fileOptions, "Number_basis_functions", &m_nbTimeBF, 1, KEYWORD_OPTIONAL, "DYNAMIC FRAMING", "ENDDF") == 1)
345  {
346  Cerr("***** iLinearModel::ReadAndCheckConfigurationFileSpecificToAllLinearModels -> Error while trying to read 'Number_basis_functions' flag in " << m_fileOptions << endl);
347  return 1;
348  }
349 
350  if( ReadDataASCIIFile(m_fileOptions, "Number_basis_functions", &m_nbRgateBF, 1, KEYWORD_OPTIONAL, "RESPIRATORY GATING", "ENDRG") == 1)
351  {
352  Cerr("***** iLinearModel::ReadAndCheckConfigurationFileSpecificToAllLinearModels -> Error while trying to read 'Number_basis_functions' flag in " << m_fileOptions << endl);
353  return 1;
354  }
355 
356  if( ReadDataASCIIFile(m_fileOptions, "Number_basis_functions", &m_nbCgateBF, 1, KEYWORD_OPTIONAL, "CARDIAC GATING", "ENDCG") == 1)
357  {
358  Cerr("***** iLinearModel::ReadAndCheckConfigurationFileSpecificToAllLinearModels -> Error while trying to read 'Number_basis_functions' flag in " << m_fileOptions << endl);
359  return 1;
360  }
361 
362  if( ReadDataASCIIFile(m_fileOptions, "Number_weight_values", &m_nbWeightFactors, 1, KEYWORD_OPTIONAL, "DYNAMIC FRAMING", "ENDDF") == 1)
363  {
364  Cerr("***** iLinearModel::ReadAndCheckConfigurationFileSpecificToAllLinearModels -> Error while trying to read 'Number_basis_functions' flag in " << m_fileOptions << endl);
365  return 1;
366  }
367 
368  if( ReadDataASCIIFile(m_fileOptions, "Number_model_iterations", &m_nbLinearModelCycles, 1, KEYWORD_OPTIONAL) == 1)
369  {
370  Cerr("***** iLinearModel::ReadAndCheckConfigurationFileSpecificToAllLinearModels -> Error while trying to read 'Number_model_iterations' flag in " << m_fileOptions << endl);
371  return 1;
372  }
373 
374  if( ReadDataASCIIFile(m_fileOptions, "Basis_function_start_ite", &m_basisFunctionsUpdStartIte, 1, KEYWORD_OPTIONAL) == 1)
375  {
376  Cerr("***** iLinearModel::ReadAndCheckConfigurationFileSpecificToAllLinearModels -> Error while trying to read 'Basis_function_start_ite' flag in " << m_fileOptions << endl);
377  return 1;
378  }
379 
380  if( ReadDataASCIIFile(m_fileOptions, "Basis_function_update_ratio", &m_basisFunctionsUpdRatio, 1, KEYWORD_OPTIONAL) == 1)
381  {
382  Cerr("***** iLinearModel::ReadAndCheckConfigurationFileSpecificToAllLinearModels -> Error while trying to read 'Basis_function_update_ratio' flag in " << m_fileOptions << endl);
383  return 1;
384  }
385  // Optimisation method
386  if( ReadDataASCIIFile(m_fileOptions, "Optimisation_method", &m_OptimisationMethod, 1, KEYWORD_MANDATORY))
387  {
388  Cerr("***** ***** iLinearModel::ReadAndCheckConfigurationFileSpecificToAllLinearModels -> Error while trying to read 'Optimisation_method' keyword in " << m_fileOptions << endl);
389  return 1;
390  }
391 
392  // Print out information on optimisation method set
393  if(m_verbose >=2)
394  {
395  if (m_OptimisationMethod == OPTIMISATION_METHOD_NNLS)
396  Cout("iLinearModel::ReadAndCheckConfigurationFileSpecificToAllLinearModels() -> Selected optimization method : NNLS (2)" << endl);
397  else if (m_OptimisationMethod == OPTIMISATION_METHOD_NESTEM)
398  Cout("iLinearModel::ReadAndCheckConfigurationFileSpecificToAllLinearModels() -> Selected optimization method : Nested-EM (1)" << endl);
399  else if (m_OptimisationMethod == OPTIMISATION_METHOD_DR)
400  Cout("iLinearModel::ReadAndCheckConfigurationFileSpecificToAllLinearModels() -> Selected optimization method : Direct Dynamic (0)" << endl);
401  }
402  }
403  else
404  {
405  Cerr("***** iLinearModel::ReadAndCheckConfigurationFileSpecificToAllLinearModels -> Error while trying to read configuration file at: " << m_fileOptions << endl);
406  return 1;
407  }
408 
409 
410  // Initialize the number of parameters of the linear model as the number of (frame) basis functions
411  // (Variable used during parametric images writing on disk)
412  m_nbModelParam = m_nbTimeBF;
413  m_nbRGModelParam = m_nbRgateBF;
414  m_nbCGModelParam = m_nbCgateBF;
415 
416  return 0;
417 }
418 
419 
420 
421 
422 // =====================================================================
423 // ---------------------------------------------------------------------
424 // ---------------------------------------------------------------------
425 // =====================================================================
426 /*
427  \fn ReadAndCheckOptionsList
428  \brief This function is used to read parameters from a string.
429  \return 0 if success, other value otherwise.
430 */
431 int iLinearModel::ReadAndCheckOptionsList(string a_listOptions)
432 {
433  if(m_verbose >=3) Cout("iLinearModel::ReadAndCheckOptionsList ..."<< endl);
434 
435  // Just recover the string here, it will be processed in the Initialize() function
436  m_listOptions = a_listOptions;
437 
438 
439  // For now, just restricts the initialization using a configuration file as there are quite a lot of parameters to initialize for a use with command line options
440  Cerr("***** iLinearModel::ReadAndCheckOptionsList() -> Initialization with command line options is not implemented for this class. Please use a configuration file instead" << endl);
441  return 1;
442 
443  // Normal end
444  //return 0;
445 }
446 
447 
448 
449 
450 // =====================================================================
451 // ---------------------------------------------------------------------
452 // ---------------------------------------------------------------------
453 // =====================================================================
454 /*
455  \fn CheckSpecificParameters
456  \brief This function is used to check whether all member variables
457  have been correctly initialized or not.
458  \return 0 if success, positive value otherwise.
459 */
460 int iLinearModel::CheckSpecificParameters()
461 {
462 
463  // Perform generic checks for the Linear Model
464  if(CheckSpecificParametersForAllLinearModels())
465  {
466  Cerr("***** iLinearModel::CheckSpecificParameters -> A problem occurred while checking specific parameters ! " << endl);
467  return 1;
468  }
469 
470  // Normal end
471  return 0;
472 }
473 
474 
475 
476 // =====================================================================
477 // ---------------------------------------------------------------------
478 // ---------------------------------------------------------------------
479 // =====================================================================
480 /*
481  \fn CheckSpecificParametersForAllLinearModels
482  \brief This function is used to check whether all member variables of a
483  generic linear model have been correctly initialized or not.
484  \return 0 if success, positive value otherwise.
485 */
486 
487 int iLinearModel::CheckSpecificParametersForAllLinearModels()
488 {
489 
490  if(m_verbose >=2) Cout("iLinearModel::CheckSpecificParametersForAllLinearModels ..."<< endl);
491 
492  // Check at least one basis function number has been initialized
493  if (m_nbTimeBF<=0 && m_nbRgateBF<=0 && m_nbCgateBF<=0)
494  {
495  Cerr("***** iLinearModel::CheckSpecificParametersForAllLinearModels() -> Error, the variables corresponding to the number of basis function has not been initialized. There might be an error in the configuration process/file !" << endl);
496  return 1;
497  }
498  else
499  {
500  // Set the other variables to 1 if not initialized
501  if(m_nbTimeBF<0) m_nbTimeBF =1;
502  if(m_nbRgateBF<0) m_nbRgateBF=1;
503  if(m_nbCgateBF<0) m_nbCgateBF=1;
504  }
505 
506  // if weights have been provided for WLS check if their lenght is equal to number of time basis functions
507  if (m_nbWeightFactors>1)
508  {
509  if (m_nbWeightFactors!=mp_ID->GetNbTimeFrames())
510  {
511  Cerr("***** iLinearModel::CheckSpecificParametersForAllLinearModels() -> Error, the number of weight factors doesn't match the number of input frames. There might be an error in the configuration process/file !" << endl);
512  return 1;
513  }
514  }
515 
516 
517  // Check if we have somehow both a file and a list of options for init...
518  if(m_listOptions != "" && m_fileOptions != "")
519  {
520  Cerr("***** iLinearModel::CheckSpecificParametersForAllLinearModels -> Either a file or a list of options have to be selected to initialize the model, but not both ! " << endl);
521  return 1;
522  }
523 
524  // Check if we have no file not list of options for some reason...
525  if(m_listOptions == "" && m_fileOptions == "")
526  {
527  Cerr("***** iLinearModel::CheckSpecificParametersForAllLinearModels -> Either a file or a list of options should have been provided at this point ! " << endl);
528  return 1;
529  }
530 
531  // In case of regular linear regression selected make sure a dynamic model with 2 basis functions is set
532  if (m_OptimisationMethod==OPTIMISATION_METHOD_LS && m_nbTimeBF!=2 )
533  {
534  Cerr("***** iLinearModel::CheckSpecificParametersForAllLinearModels -> Optimization method Linear Regression (=3) only available when using two time basis functions ! " << endl);
535  return 1;
536  }
537  // Normal end
538  return 0;
539 }
540 
541 
542 // =====================================================================
543 // ---------------------------------------------------------------------
544 // ---------------------------------------------------------------------
545 // =====================================================================
546 /*
547  \fn InitializeSpecific
548  \brief This function is used to initialize the parametric images and basis functions
549  \return 0 if success, other value otherwise.
550 */
551 
552 int iLinearModel::InitializeSpecificToAllLinearModels()
553 {
554  if(m_verbose >=2) Cout("iLinearModel::InitializeSpecificToAllLinearModels ..."<< endl);
555 
556  // Forbid initialization without check
557  if (!m_checked)
558  {
559  Cerr("***** iLinearModel::InitializeSpecificToAllLinearModels() -> Must call CheckParameters functions before Initialize() !" << endl);
560  return 1;
561  }
562 
563  // Initialize blacklisted voxels image to zero;
564  mp_blackListedvoxelsImage = new FLTNB [mp_ID->GetNbVoxXYZ()];
565  for (int v = 0; v < mp_ID->GetNbVoxXYZ(); v++)
566  {
567  mp_blackListedvoxelsImage[v] = 0.0;
568  }
569 
570  // Allocate memory for Parametric images and time basis functions of the model
571  m2p_parametricImages = new FLTNB*[m_nbTimeBF];
572  m2p_outputParImages = new FLTNB*[m_nbTimeBF];
573  m2p_nestedModelTimeBasisFunctions = new FLTNB*[m_nbTimeBF];
574 
575  for(int b=0 ; b<m_nbTimeBF ; b++)
576  {
577  m2p_nestedModelTimeBasisFunctions[b] = new FLTNB[mp_ID->GetNbTimeFrames()];
578  m2p_parametricImages[b] = new FLTNB[mp_ID->GetNbVoxXYZ()];
579  }
580 
581  m2p_respBasisFunctions = new FLTNB*[m_nbRgateBF];
582  m2p_RGParametricImages = new FLTNB*[m_nbRgateBF];
583 
584  for(int rb=0 ; rb<m_nbRgateBF ; rb++)
585  {
586  m2p_respBasisFunctions[rb] = new FLTNB[mp_ID->GetNbRespGates()];
587  m2p_RGParametricImages[rb] = new FLTNB[mp_ID->GetNbVoxXYZ()];
588  }
589 
590  m2p_cardBasisFunctions = new FLTNB*[m_nbCgateBF];
591  m2p_CGParametricImages = new FLTNB*[m_nbCgateBF];
592 
593 
594  for(int cb=0 ; cb<m_nbCgateBF ; cb++)
595  {
596  m2p_cardBasisFunctions[cb] = new FLTNB[mp_ID->GetNbCardGates()];
597  m2p_CGParametricImages[cb] = new FLTNB[mp_ID->GetNbVoxXYZ()];
598  }
599 
600  // Memory allocation for correction images
601  mp_corrBasisCoeffs = new FLTNB[mp_ID->GetNbVoxXYZ()];
602  mp_corrBasisFunctions = new FLTNB[mp_ID->GetNbThreadsForProjection()];
603 
604  // Allocate output image matrices
605  for(int b=0 ; b<m_nbTimeBF ; b++)
606  m2p_outputParImages[b] = new FLTNB[mp_ID->GetNbVoxXYZ()];
607 
608 
609  // TODO: Consider other cases where linear models will be applied for Cardiac / Respiratory reconstruction.
610 
611  // --- Default Initialization basis functions --- //
612  for(int b=0 ; b<m_nbTimeBF ; b++)
613  for (int fr=0 ; fr<mp_ID->GetNbTimeFrames() ; fr++)
614  m2p_nestedModelTimeBasisFunctions[b][fr] = 1.;
615 
616  // --- Default Initialization of respiratory and cardiac basis functions --- //
617  for(int rb=0 ; rb<m_nbRgateBF ; rb++)
618  for (int rg=0 ; rg<mp_ID->GetNbRespGates() ; rg++)
619  m2p_respBasisFunctions[rb][rg] = 1.;
620 
621  for(int cb=0 ; cb<m_nbCgateBF ; cb++)
622  for (int cg=0 ; cg<mp_ID->GetNbCardGates() ; cg++)
623  m2p_cardBasisFunctions[cb][cg] = 1.;
624 
625 
626  // Initialize gate parametric images of linear model mother class
627  for(int v=0 ; v<mp_ID->GetNbVoxXYZ() ; v++)
628  {
629  m2p_RGParametricImages[ 0 ][ v ] = 1.;
630  m2p_CGParametricImages[ 0 ][ v ] = 1.;
631  }
632 
633  // Allocate NNLS variables and weights.
634  if(m_OptimisationMethod == OPTIMISATION_METHOD_NNLS)
635  {
636  m3p_nnlsA = new FLTNB** [mp_ID->GetNbThreadsForImageComputation()];
637  m2p_nnlsB = new FLTNB* [mp_ID->GetNbThreadsForImageComputation()];
638  m2p_nnlsMat = new FLTNB* [mp_ID->GetNbThreadsForImageComputation()];
639  m2p_nnlsX = new FLTNB* [mp_ID->GetNbThreadsForImageComputation()];
640  m2p_nnlsWp = new FLTNB* [mp_ID->GetNbThreadsForImageComputation()];
641  m2p_nnlsIdx = new int* [mp_ID->GetNbThreadsForImageComputation()];
642  mp_w = new FLTNB[ mp_ID->GetNbTimeFrames() ];
643 
644  // Get the largest nb of parameters and samples between time frames, resp or cardiac gates
645  // as all could be used as nnls samples depending on the linear model
646  int nnls_max_params = (m_nbTimeBF > m_nbRgateBF) ?
647  ( (m_nbTimeBF > m_nbCgateBF) ? m_nbTimeBF : m_nbCgateBF) :
648  ( (m_nbRgateBF > m_nbCgateBF) ? m_nbRgateBF : m_nbCgateBF) ;
649 
650  m_nnlsN = nnls_max_params;
651 
652  int nnls_max_samples = (mp_ID->GetNbTimeFrames() > mp_ID->GetNbRespGates()) ?
653  ( (mp_ID->GetNbTimeFrames() > mp_ID->GetNbCardGates()) ? mp_ID->GetNbTimeFrames() : mp_ID->GetNbCardGates()) :
654  ( (mp_ID->GetNbRespGates() > mp_ID->GetNbCardGates()) ? mp_ID->GetNbRespGates() : mp_ID->GetNbCardGates()) ;
655 
656 
657  for( int th=0 ; th<mp_ID->GetNbThreadsForImageComputation() ; th++ )
658  {
659  // Init 2D coefficient matrix for NNLS estimation
660  m3p_nnlsA[ th ] = new FLTNB*[ m_nnlsN ];
661 
662  for(int n=0 ; n<m_nnlsN ; n++)
663  m3p_nnlsA[ th ][n] = new FLTNB[ nnls_max_samples ];
664 
665  // Init solution vector and working matrix for NNLS estimation
666 
667  m2p_nnlsB[ th ] = new FLTNB[ nnls_max_samples ];
668  m2p_nnlsMat[ th ] = new FLTNB[ (m_nnlsN+2) * nnls_max_samples ];
669  m2p_nnlsX[ th ] = new FLTNB[ m_nnlsN ];
670  m2p_nnlsWp[ th ] = new FLTNB[ m_nnlsN ];
671  m2p_nnlsIdx[ th ] = new int[ m_nnlsN ];
672  }
673 
674  // Fill weights for NNLS to default values
675  // todo: specific weights for gate-related models (gate duration) - Also NECR related weights or user set weights (ZC)
676 
677  // Case of no weights input - set all weights to 1
678  if (m_nbWeightFactors<=1)
679  for(int t=0; t<mp_ID->GetNbTimeFrames(); t++)
680  mp_w[t] = 1;
681  }
682 
683 
684  // --- Data Initialization with a configuration file --- //
685  if(m_fileOptions != "")
686  {
687  ifstream in_file(m_fileOptions.c_str(), ios::in);
688 
689  if(in_file)
690  {
691  // TODO: Implement a check of length of input basis functions
692  // Frame basis functions Initialization
693  if(m_nbTimeBF>1
694  && ReadDataASCIIFile(m_fileOptions,
695  "Basis_functions",
696  m2p_nestedModelTimeBasisFunctions,
697  mp_ID->GetNbTimeFrames(),
698  m_nbTimeBF,
699  KEYWORD_OPTIONAL,
700  "DYNAMIC FRAMING",
701  "ENDDF")==1)
702  {
703  Cerr("***** iLinearModel::Initialize -> Error while trying to read frame basis functions coefficients !" << endl);
704  Cerr(" 'Basis_functions' keyword inside DYNAMIC FRAMING / ENDDF paragraph in " << m_fileOptions << endl);
705  return 1;
706  }
707 
708  // Resp gates basis functions Initialization
709  if(m_nbRgateBF>1
710  && ReadDataASCIIFile(m_fileOptions,
711  "Basis_functions",
712  m2p_respBasisFunctions,
713  mp_ID->GetNbRespGates(),
714  m_nbRgateBF,
715  KEYWORD_OPTIONAL,
716  "RESPIRATORY GATING",
717  "ENDRG") ==1 )
718  {
719  Cerr("***** iLinearModel::Initialize -> Error while trying to read respiratory gates basis functions coefficients !" << endl);
720  Cerr(" 'Basis_functions' keyword inside RESPIRATORY GATING / ENDRG paragraph in " << m_fileOptions << endl);
721  return 1;
722  }
723 
724  // Card gates basis functions Initialization
725  if(m_nbCgateBF>1
726  && ReadDataASCIIFile(m_fileOptions,
727  "Basis_functions",
728  m2p_cardBasisFunctions,
729  mp_ID->GetNbCardGates(),
730  m_nbCgateBF,
731  KEYWORD_OPTIONAL,
732  "CARDIAC GATING",
733  "ENDCG") ==1 )
734  {
735  Cerr("***** iLinearModel::Initialize -> Error while trying to read cardiac gates basis functions coefficients !" << endl);
736  Cerr(" 'Basis_functions' keyword inside CARDIAC GATING / ENDCG paragraph in " << m_fileOptions << endl);
737  return 1;
738  }
739 
740  // Optimisation method
741  if( ReadDataASCIIFile(m_fileOptions,
742  "Optimisation_method",
743  &m_OptimisationMethod,
744  1,
745  KEYWORD_MANDATORY))
746  {
747  Cerr("***** iLinearModel::Initialize -> Error while trying to read 'Optimisation_method' keyword in " << m_fileOptions << endl);
748  return 1;
749  }
750 
751  //
752 
753  if(m_nbWeightFactors>1
754  && ReadDataASCIIFile(m_fileOptions,
755  "Weight_values",
756  mp_w,
757  m_nbWeightFactors,
758  KEYWORD_OPTIONAL,
759  "DYNAMIC FRAMING",
760  "ENDDF")==1)
761  {
762  Cerr("***** iLinearModel::Initialize -> Error while trying to read frame weight factor coefficients !" << endl);
763  Cerr(" 'Weight_values' keyword inside DYNAMIC FRAMING / ENDDF paragraph in " << m_fileOptions << endl);
764  return 1;
765  }
766 
767 
768  // --- Parametric image initialization --- //
769 
770  // Frame model
771  string input_image = "";
772  int return_value = 0;
773 
774  return_value = ReadDataASCIIFile(m_fileOptions,
775  "Parametric_images_init",
776  &input_image,
777  1,
778  KEYWORD_OPTIONAL,
779  "DYNAMIC FRAMING",
780  "ENDDF");
781 
782  if( return_value == 0) // Image have been provided
783  {
784  // Read image // INTF_LERP_DISABLED = interpolation disabled for input image reading
785  if( IntfReadImgDynCoeffFile(input_image,
786  m2p_parametricImages,
787  mp_ID,
788  m_nbModelParam,
789  m_verbose,
790  INTF_LERP_DISABLED) ) // Image have been provided
791  {
792  Cerr("***** iLinearModel::Initialize -> Error while trying to read the provided initialization parametric images : " << input_image << endl);
793  return 1;
794  }
795  }
796  else if( return_value == 1) // Error during reading
797  {
798  Cerr("***** iLinearModel::Initialize -> Error while trying to read dynamic frame model parametric images !" << endl);
799  Cerr(" 'Parametric_image_init' keyword in " << m_fileOptions << endl);
800  return 1;
801  }
802  else //(return_value >= 1 ) // Keyword not found : no initialization provided
803  {
804  // Standard initialization
805  for(int b=0 ; b<m_nbTimeBF ; b++)
806  for(int v=0 ; v<mp_ID->GetNbVoxXYZ() ; v++)
807  m2p_parametricImages[b][v] = 1.;
808  }
809 
810 
811  // Resp gate model
812  input_image = "";
813  return_value = 0;
814 
815  return_value = ReadDataASCIIFile(m_fileOptions,
816  "Parametric_images_init",
817  &input_image,
818  1,
819  KEYWORD_OPTIONAL,
820  "RESPIRATORY GATING",
821  "ENDRG");
822 
823  if( return_value == 0) // Image have been provided
824  {
825  // Read image // INTF_LERP_DISABLED = interpolation disabled for input image reading
826  if( IntfReadImgDynCoeffFile(input_image,
827  m2p_RGParametricImages,
828  mp_ID,
829  m_nbRGModelParam,
830  m_verbose,
831  INTF_LERP_DISABLED) ) // Image have been provided
832  {
833  Cerr("***** iLinearModel::Initialize -> Error while trying to read the provided initialization parametric images : " << input_image << endl);
834  return 1;
835  }
836  }
837  else if( return_value == 1) // Error during reading
838  {
839  Cerr("***** iLinearModel::Initialize -> Error while trying to read respiratory gate model parametric images !" << endl);
840  Cerr(" 'Parametric_image_init' keyword in " << m_fileOptions << endl);
841  return 1;
842  }
843  else //(return_value >= 1 ) // Keyword not found : no initialization provided
844  {
845  // Standard initialization
846  for(int rb=0 ; rb<m_nbRgateBF ; rb++)
847  for (int v=0 ; v<mp_ID->GetNbVoxXYZ() ; v++)
848  m2p_RGParametricImages[rb][v] = 1.;
849  }
850 
851 
852 
853  // Card gate model
854  input_image = "";
855  return_value = 0;
856 
857  return_value = ReadDataASCIIFile(m_fileOptions,
858  "Parametric_images_init",
859  &input_image,
860  1,
861  KEYWORD_OPTIONAL,
862  "CARDIAC GATING",
863  "ENDCG");
864 
865  if( return_value == 0) // Image have been provided
866  {
867  // Read image // INTF_LERP_DISABLED = interpolation disabled for input image reading
868  if( IntfReadImgDynCoeffFile(input_image,
869  m2p_CGParametricImages,
870  mp_ID,
871  m_nbCGModelParam,
872  m_verbose,
873  INTF_LERP_DISABLED) ) // Image have been provided
874  {
875  Cerr("***** iLinearModel::Initialize -> Error while trying to read the provided initialization parametric images : " << input_image << endl);
876  return 1;
877  }
878  }
879  else if( return_value == 1) // Error during reading
880  {
881  Cerr("***** iLinearModel::Initialize -> Error while trying to read cardiac gate model parametric images !" << endl);
882  Cerr(" 'Parametric_image_init' keyword in " << m_fileOptions << endl);
883  return 1;
884  }
885  else //(return_value >= 1 ) // Keyword not found : no initialization provided
886  {
887  // Standard initialization
888  for(int cb=0 ; cb<m_nbCgateBF ; cb++)
889  for (int v=0 ; v<mp_ID->GetNbVoxXYZ() ; v++)
890  m2p_CGParametricImages[cb][v] = 1.;
891  }
892  }
893 
894  else
895  {
896  Cerr("***** iLinearModel::Initialize() -> Error while trying to read configuration file at: " << m_fileOptions << endl);
897  return 1;
898  }
899  }
900 
901  // If weight values have been provided check they are not negative
902  if (m_nbWeightFactors>1)
903  for(int t=0; t<mp_ID->GetNbTimeFrames(); t++)
904  {
905  if (mp_w[t]<=0)
906  {
907  Cerr("***** iLinearModel::CheckSpecificParametersForAllLinearModels() -> Error, negative weight factors found. Only positive non-zero factors allowed !" << endl);
908  return 1;
909  }
910  if(m_verbose >=4) Cout("iLinearModel::NNLS optimisation weight for frame: "<< t << " set at "<< mp_w[t]<< endl);
911  }
912 
913 
914 
915  // --- Data Initialization with a list of options --- //
916 
917  if(m_listOptions != "")
918  {
919  // TODO
920  }
921 
922  // Allocate output image matrices
923  for(int b=0 ; b<m_nbTimeBF ; b++)
924  m2p_outputParImages[b] = new FLTNB[mp_ID->GetNbVoxXYZ()];
925 
926  // If optimisation method requires an input of basis functions set m_ModelSpecificBasisFunctionsRequired flat
927  if (m_OptimisationMethod==OPTIMISATION_METHOD_DR)
928  {
929  m_ModelSpecificBasisFunctionsRequired = true;
930  m_noImageUpdateFlag = true;
931  }
932 
933  // TODO : print output parametric images rg et cg
934 
935 
936  // Normal End
937  return 0;
938 }
939 
940 
941 // =====================================================================
942 // ---------------------------------------------------------------------
943 // ---------------------------------------------------------------------
944 // =====================================================================
945 /*
946  \fn ShowBasisFunctions
947  \brief This function is used to initialize the parametric images and basis functions
948  \return 0 if success, other value otherwise.
949 */
950 
951 void iLinearModel::ShowBasisFunctions()
952 {
953 
954  // Display Basis Functions set for the model
955  if ((m_verbose>=2) & (m_nbTimeBF>1))
956  {
957  Cout("***** iLinearModel::InitializeSpecificToAllLinearModels() -> Time Basis Function coefficients :" << endl);
958  for(int b=0 ; b<m_nbTimeBF ; b++)
959  {
960  Cout(" ");
961  Cout("Basis function["<<b+1<<"] : ");
962  for(int fr=0 ; fr<mp_ID->GetNbTimeFrames()-1 ; fr++)
963  {
964  Cout(m2p_nestedModelTimeBasisFunctions[b][fr] << ", ");
965  }
966  // Print last value
967  Cout(m2p_nestedModelTimeBasisFunctions[b][mp_ID->GetNbTimeFrames()-1] << endl);
968  }
969  }
970 
971 
972 }
973 
974 // =====================================================================
975 // ---------------------------------------------------------------------
976 // ---------------------------------------------------------------------
977 // =====================================================================
978 /*
979  \fn InitializeSpecific
980  \brief This function is used to initialize the parametric images and basis functions
981  \return 0 if success, other value otherwise.
982 */
983 int iLinearModel::InitializeSpecific()
984 {
985  if(m_verbose >=2) Cout("iLinearModel::InitializeSpecific ..."<< endl);
986 
987  // Forbid initialization without check
988  if (!m_checked)
989  {
990  Cerr("***** oDynamicModelManager::InitializeSpecific() -> Must call CheckParameters functions before Initialize() !" << endl);
991  return 1;
992  }
993 
994  // Run generic Initialization for all Linear Models
995  if (InitializeSpecificToAllLinearModels())
996  {
997  Cerr("***** iLinearModel::InitializeSpecific() -> Error while performing generic initialisations for linear models !" << endl);
998  return 1;
999  }
1000 
1001 
1002  // Normal end
1003  m_initialized = true;
1004  return 0;
1005 }
1006 
1007 
1008 
1009 
1010 // =====================================================================
1011 // ---------------------------------------------------------------------
1012 // ---------------------------------------------------------------------
1013 // =====================================================================
1014 /*
1015  \fn EstimateModelParameters
1016  \param ap_ImageS : pointer to the ImageSpace
1017  \param a_ite : index of the actual iteration (not used)
1018  \param a_sset : index of the actual subset (not used)
1019  \brief Estimate model parameters (parametric images and basis functions)
1020  \return 0 if success, other value otherwise.
1021 */
1022 int iLinearModel::EstimateModelParameters(oImageSpace* ap_ImageS, int a_ite, int a_sset)
1023 {
1024 
1025  #ifdef CASTOR_DEBUG
1026  if (!m_initialized)
1027  {
1028  Cerr("***** iLinearModel::EstimateModelParameters() -> Called while not initialized !" << endl);
1029  Exit(EXIT_DEBUG);
1030  }
1031  #endif
1032 
1033  if(m_verbose >=3) Cout(" iLinearModel::EstimateModelParameters -> Optimisation method to use: " << m_OptimisationMethod << endl) ;
1034 
1035  // Use the nested EM parametric image estimation method
1036  if(m_OptimisationMethod == OPTIMISATION_METHOD_NESTEM)
1037  {
1038  if(NestedEM(ap_ImageS, a_ite) )
1039  {
1040  Cerr("***** iLinearModel::EstimateModelParameters() -> An error occured while using the nested EM parametric image estimation method !" << endl);
1041  return 1;
1042  }
1043  }
1044  // Direct reconstruction ( non nested ) - no update required
1045  else if(m_OptimisationMethod == OPTIMISATION_METHOD_DR)
1046  {
1047  Cout("iLinearModel::EstimateModelParameters() -> Skipping nested calculation, not required for direct method " << endl);
1048  }
1049 
1050  // Least Square for Post Reconsutrction Patlak Analysis
1051  else if(m_OptimisationMethod == OPTIMISATION_METHOD_NNLS )
1052  {
1053  if(EstimateParametersWithNNLS(ap_ImageS, a_ite) )
1054  {
1055  Cerr("***** iLinearModel::EstimateModelParameters() -> An error occured while using the Patlak LS Estimation Method !" << endl);
1056  return 1;
1057  }
1058  }
1059  // least-square linear regression
1060  else if(m_OptimisationMethod == OPTIMISATION_METHOD_LS )
1061  {
1062  if(Patlak_LS(ap_ImageS, a_ite) )
1063  {
1064  Cerr("***** iLinearModel::EstimateModelParameters() -> An error occured while using the Least-Squares linear regression Estimation Method !" << endl);
1065  return 1;
1066  }
1067  }
1068 
1069  else
1070  {
1071  Cerr("***** iLinearModel::EstimateModelParameters() -> Error : unknown method to estimate images ! !" << endl);
1072  return 1;
1073  }
1074 
1075  if(m_verbose >=3)
1076  {
1077  for (int fb=0 ; fb<m_nbTimeBF ; fb++)
1078  {
1079  FLTNB avg = 0;
1080  for (int v=0 ; v<mp_ID->GetNbVoxXYZ() ; v++)
1081  if(mp_maskModel != NULL && mp_maskModel[ v ] == 1)
1082  avg += m2p_parametricImages[fb][v];
1083 
1084  avg /= m_nbVoxelsMask;
1085  Cout( "iLinearModel::EstimateModelParameters() -> Frame parametric image["<<fb<<"] avg value:" << avg << endl);
1086  }
1087  }
1088 
1089  return 0;
1090 }
1091 
1092 
1093 // =====================================================================
1094 // ---------------------------------------------------------------------
1095 // ---------------------------------------------------------------------
1096 // =====================================================================
1097 /*
1098  \fn NestedEM
1099  \param ap_ImageS : pointer to the ImageSpace
1100  \param a_ite : index of the actual iteration (not used)
1101  \brief Estimate parametric images and basis functions (if enabled) using the nested EM method
1102  \return 0 if success, other value otherwise.
1103 */
1104 int iLinearModel::NestedEM(oImageSpace* ap_ImageS, int a_ite)
1105 {
1106  if(m_verbose >=3) Cout("iLinearModel::NestedEM() ..." <<endl);
1107 
1108  for (uint32_t it=0 ; it<m_nbLinearModelCycles ; it++)
1109  {
1110  if(m_verbose >=3) Cout("iLinearModel::NestedEM() cycle "<< it+1 << "/" << m_nbLinearModelCycles <<endl);
1111 
1112  // Step 1 : Generate a difference image from voxel-by-voxel division of the current estimation of the image and the model image generated from the current parametric images / basis functions of the model
1113  // The backward image matrix (which is useless at this point of the reconstruction) is used as provisional image to gather the model image
1114 
1115  int v;
1116  #pragma omp parallel for private(v) schedule(static, 1)
1117  for (int fr=0 ; fr<mp_ID->GetNbTimeFrames() ; fr++)
1118  for (int rg=0 ; rg<mp_ID->GetNbRespGates() ; rg++)
1119  for (int cg=0 ; cg<mp_ID->GetNbCardGates() ; cg++)
1120  for (v=0 ; v<mp_ID->GetNbVoxXYZ() ; v++)
1121  {
1122  // If a mask has been provided, check if the model applies in this voxel
1123  if(mp_maskModel != NULL && mp_maskModel[ v ] == 0)
1124  continue;
1125 
1126  // Reset this voxel to 0
1127  ap_ImageS->m6p_backwardImage[0][0][fr][rg][cg][v] = 0;
1128 
1129  for (int fb=0 ; fb<m_nbTimeBF ; fb++)
1130  for (int rb=0 ; rb<m_nbRgateBF ; rb++)
1131  for (int cb=0 ; cb<m_nbCgateBF ; cb++)
1132  {
1133  // Retrieve current estimation of image according to coeffs/basis functions
1134  ap_ImageS->m6p_backwardImage[0][0][fr][rg][cg][v] += m2p_parametricImages[fb][v]
1135  * m2p_nestedModelTimeBasisFunctions[fb][fr]
1136  * m2p_RGParametricImages[rb][v]
1137  * m2p_respBasisFunctions[rb][rg]
1138  * m2p_CGParametricImages[cb][v]
1139  * m2p_cardBasisFunctions[cb][cg] ;
1140  }
1141 
1142  // Recover correction images using the ratio of the current estimation of the image (m4p_image) and the image generated with coffs/basis functions
1143  // (Again, use backward image as provisional image to recover the result)
1144 
1145  // TODO: Find where NaN is generated
1146  if (ap_ImageS->m6p_backwardImage[0][0][fr][rg][cg][v] == 0 )
1147  {
1148  mp_blackListedvoxelsImage[v] = 1.0;
1149  ap_ImageS->m6p_backwardImage[0][0][fr][rg][cg][v] = 1;
1150  }
1151  else
1152  {
1153  ap_ImageS->m6p_backwardImage[0][0][fr][rg][cg][v] =
1154  ap_ImageS->m4p_image[fr][rg][cg][v] / ap_ImageS->m6p_backwardImage[0][0][fr][rg][cg][v];
1155  }
1156  }
1157 
1158  // Step 2 : Estimate either basis functions or parametric images
1159 
1160  // Step 2a : Basis functions estimation // TODO: shouldn't this condition be a single & as we want all conditions ??
1161  if( m_basisFunctionsUpdStartIte >= 0 // if (m_basisFunctionsUpdStartIte < 0) --> no update of basis functions, just update the parametric images
1162  && m_basisFunctionsUpdStartIte <= a_ite // if (m_basisFunctionsUpdStartIte >= 0) --> Check condition of minimal iteration before starting to update the basis functions
1163  && !(int((m_basisFunctionsUpdIdx-1)/m_basisFunctionsUpdRatio)&1) ) // Check if we must estimate basis functions or coefficients at this stage (depends on m_basisFunctionsUpdRatio)
1164  {
1165 
1166  FLTNB parametric_image_norm = 0;
1167 
1168  if( m_nbTimeBF>1 )
1169  {
1170  if(m_verbose >=3) Cout("iLinearModel::NestedEM() -> Estimate Basis Functions - Frame basis functions estimation step" <<endl);
1171 
1172  for (int fb=0 ; fb<m_nbTimeBF ; fb++)
1173  {
1174  // Compute normalization related to the vowelwise time basis functions coefficients.
1175  parametric_image_norm = 0;
1176 
1177  // Compute normalization related to the parametric images.
1178  for (int rb=0 ; rb<m_nbRgateBF ; rb++)
1179  for (int cb=0 ; cb<m_nbRgateBF ; cb++)
1180  for (int v=0 ; v<mp_ID->GetNbVoxXYZ() ; v++)
1181  if(mp_maskModel != NULL && mp_maskModel[ v ] == 1)
1182  parametric_image_norm += m2p_parametricImages[fb][v]
1183  * m2p_RGParametricImages[rb][v]
1184  * m2p_CGParametricImages[cb][v];
1185 
1186  parametric_image_norm *= mp_ID->GetNbRespGates()
1187  * mp_ID->GetNbCardGates();
1188 
1189  for (int fr=0 ; fr<mp_ID->GetNbTimeFrames() ; fr++)
1190  {
1191  // Initialization of corrections factor for time basis functions
1192  for (int th=0 ; th<mp_ID->GetNbThreadsForImageComputation() ; th++)
1193  mp_corrBasisFunctions[th]=0;
1194 
1195  for (int rg=0 ; rg<mp_ID->GetNbRespGates() ; rg++)
1196  for (int cg=0 ; cg<mp_ID->GetNbCardGates() ; cg++)
1197  for (int rb=0 ; rb<m_nbRgateBF ; rb++)
1198  for (int cb=0 ; cb<m_nbCgateBF ; cb++)
1199  {
1200  int v;
1201  #pragma omp parallel for private(v) schedule(static, 1)
1202  for (v=0 ; v<mp_ID->GetNbVoxXYZ() ; v++)
1203  {
1204  // If a mask has been provided, check if the model applies in this voxel
1205  if(mp_maskModel != NULL && mp_maskModel[ v ] == 0)
1206  continue;
1207 
1208  int th = 0;
1209  #ifdef CASTOR_OMP
1210  th = omp_get_thread_num();
1211  #endif
1212 
1213  mp_corrBasisFunctions[th] += ap_ImageS->m6p_backwardImage[0][0][fr][rg][cg][v]
1214  * m2p_parametricImages[fb][v]
1215  * m2p_RGParametricImages[rb][v]
1216  * m2p_CGParametricImages[cb][v];
1217  }
1218  }
1219 
1220  // Reduce
1221  for (int th=1 ; th<mp_ID->GetNbThreadsForImageComputation() ; th++)
1222  mp_corrBasisFunctions[0] += mp_corrBasisFunctions[th];
1223 
1224  // Apply corrections and normalization to the temporal basis functions values
1225  if (mp_corrBasisFunctions[0] > 0.)
1226  m2p_nestedModelTimeBasisFunctions[fb][fr] *= mp_corrBasisFunctions[0]/parametric_image_norm;
1227  }
1228 
1229  }
1230  }
1231 
1232  if( m_nbRgateBF>1 )
1233  {
1234  if(m_verbose >=3) Cout("iLinearModel::NestedEM() -> Estimate Basis Functions - Respiratory gate basis functions estimation step" <<endl);
1235 
1236  for (int rb=0 ; rb<m_nbRgateBF ; rb++)
1237  {
1238  // Compute normalization related to the vowelwise time basis functions coefficients.
1239  parametric_image_norm = 0;
1240 
1241  for (int fb=0 ; fb<m_nbTimeBF ; fb++)
1242  for (int cb=0 ; cb<m_nbRgateBF ; cb++)
1243  for (int v=0 ; v<mp_ID->GetNbVoxXYZ() ; v++)
1244  if(mp_maskModel != NULL && mp_maskModel[ v ] == 1)
1245  parametric_image_norm += m2p_parametricImages[fb][v]
1246  * m2p_RGParametricImages[rb][v]
1247  * m2p_CGParametricImages[cb][v];
1248 
1249  parametric_image_norm *= mp_ID->GetNbTimeFrames()
1250  * mp_ID->GetNbCardGates();
1251 
1252  for (int rg=0 ; rg<mp_ID->GetNbRespGates() ; rg++)
1253  {
1254  // Initialization of corrections for time basis functions
1255  for (int th=0 ; th<mp_ID->GetNbThreadsForImageComputation() ; th++)
1256  mp_corrBasisFunctions[th]=0;
1257 
1258  for (int fr=0 ; fr<mp_ID->GetNbTimeFrames() ; fr++)
1259  for (int cg=0 ; cg<mp_ID->GetNbCardGates() ; cg++)
1260  for (int fb=0 ; fb<m_nbTimeBF ; fb++)
1261  for (int cb=0 ; cb<m_nbCgateBF ; cb++)
1262  {
1263  int v;
1264  #pragma omp parallel for private(v) schedule(static, 1)
1265  for (v=0 ; v<mp_ID->GetNbVoxXYZ() ; v++)
1266  {
1267  // If a mask has been provided, check if the model applies in this voxel
1268  if(mp_maskModel != NULL && mp_maskModel[ v ] == 0)
1269  continue;
1270 
1271  int th = 0;
1272  #ifdef CASTOR_OMP
1273  th = omp_get_thread_num();
1274  #endif
1275  mp_corrBasisFunctions[th] += ap_ImageS->m6p_backwardImage[0][0][fr][rg][cg][v]
1276  * m2p_parametricImages[fb][v]
1277  * m2p_RGParametricImages[rb][v]
1278  * m2p_CGParametricImages[cb][v];
1279  }
1280  }
1281 
1282  // Reduce
1283  for (int th=1 ; th<mp_ID->GetNbThreadsForImageComputation() ; th++)
1284  mp_corrBasisFunctions[0] += mp_corrBasisFunctions[th];
1285 
1286  // Apply corrections and normalization to the temporal basis functions values
1287  if (mp_corrBasisFunctions[0] > 0.)
1288  m2p_respBasisFunctions[rb][rg] *= mp_corrBasisFunctions[0]/parametric_image_norm;
1289 
1290  }
1291  }
1292  }
1293 
1294 
1295  if( m_nbCgateBF>1 )
1296  {
1297  if(m_verbose >=3) Cout("iLinearModel::NestedEM() -> Estimate Basis Functions - Cardiac gate basis functions estimation step" <<endl);
1298 
1299  for (int cb=0 ; cb<m_nbCgateBF ; cb++)
1300  {
1301  // Compute normalization related to the vowelwise time basis functions coefficients.
1302  parametric_image_norm = 0;
1303 
1304  for (int fb=0 ; fb<m_nbTimeBF ; fb++)
1305  for (int rb=0 ; rb<m_nbRgateBF ; rb++)
1306  for (int v=0 ; v<mp_ID->GetNbVoxXYZ() ; v++)
1307  if(mp_maskModel != NULL && mp_maskModel[ v ] == 1)
1308  parametric_image_norm += m2p_parametricImages[fb][v]
1309  * m2p_RGParametricImages[rb][v]
1310  * m2p_CGParametricImages[cb][v];
1311 
1312  parametric_image_norm *= mp_ID->GetNbTimeFrames()
1313  * mp_ID->GetNbRespGates();
1314 
1315  for (int cg=0 ; cg<mp_ID->GetNbCardGates() ; cg++)
1316  {
1317  // Initialization of corrections for time basis functions
1318  for (int th=0 ; th<mp_ID->GetNbThreadsForImageComputation() ; th++)
1319  mp_corrBasisFunctions[th]=0;
1320 
1321  for (int fr=0 ; fr<mp_ID->GetNbTimeFrames() ; fr++)
1322  for (int rg=0 ; rg<mp_ID->GetNbRespGates() ; rg++)
1323  for (int fb=0 ; fb<m_nbTimeBF ; fb++)
1324  for (int rb=0 ; rb<m_nbRgateBF ; rb++)
1325  {
1326  int v;
1327  #pragma omp parallel for private(v) schedule(static, 1)
1328  for (v=0 ; v<mp_ID->GetNbVoxXYZ() ; v++)
1329  {
1330  // If a mask has been provided, check if the model applies in this voxel
1331  if(mp_maskModel != NULL && mp_maskModel[ v ] == 0)
1332  continue;
1333 
1334  int th = 0;
1335  #ifdef CASTOR_OMP
1336  th = omp_get_thread_num();
1337  #endif
1338  mp_corrBasisFunctions[th] += ap_ImageS->m6p_backwardImage[0][0][fr][rg][cg][v]
1339  * m2p_parametricImages[fb][v]
1340  * m2p_RGParametricImages[rb][v]
1341  * m2p_CGParametricImages[cb][v];
1342  }
1343  }
1344 
1345  // Reduce
1346  for (int th=1 ; th<mp_ID->GetNbThreadsForImageComputation() ; th++)
1347  mp_corrBasisFunctions[0] += mp_corrBasisFunctions[th];
1348 
1349  // Apply corrections and normalization to the temporal basis functions values
1350  if (mp_corrBasisFunctions[0] > 0.)
1351  m2p_cardBasisFunctions[cb][cg] *= mp_corrBasisFunctions[0]/parametric_image_norm;
1352 
1353  }
1354  }
1355  }
1356 
1357  m_basisFunctionsUpdIdx++;
1358 
1359  // Some feedback :
1360  if(m_verbose >=3)
1361  {
1362  for (int fb=0 ; fb<m_nbTimeBF ; fb++)
1363  for (int fr=0 ; fr<mp_ID->GetNbTimeFrames() ; fr++)
1364  Cout( "iLinearModel::NestedEM() -> Basis function ["<<fb<<"] coefficients for frame ["<<fr<<"] :" << m2p_nestedModelTimeBasisFunctions[fb][fr] << endl);
1365 
1366  for (int rb=0 ; rb<m_nbRgateBF ; rb++)
1367  for (int rg=0 ; rg<mp_ID->GetNbRespGates() ; rg++)
1368  Cout( "iLinearModel::NestedEM() -> Basis function ["<<rb<<"] coefficients for resp gate ["<<rg<<"] :" << m2p_respBasisFunctions[rb][rg] << endl);
1369 
1370  for (int cb=0 ; cb<m_nbCgateBF ; cb++)
1371  for (int cg=0 ; cg<mp_ID->GetNbCardGates() ; cg++)
1372  Cout( "iLinearModel::NestedEM() -> Basis function ["<<cb<<"] coefficients for card gate ["<<cg<<"] :" << m2p_cardBasisFunctions[cb][cg] << endl);
1373  }
1374  } // end of if loop (basis functions update)
1375 
1376 
1377 
1378  // Step 2b : Image Coefficients estimation
1379  else
1380  {
1381  // Normalization general factor
1382  FLTNB basis_functions_norm = 0.;
1383 
1384  // Regularisation according to frame basis functions
1385  if( m_nbModelParam>1 )
1386  {
1387  if(m_verbose >=3) Cout("iLinearModel::NestedEM() -> Estimate Coefficients - Frame coeffs estimation step" <<endl);
1388 
1389  for (int fb=0 ; fb<m_nbTimeBF ; fb++)
1390  {
1391  // Initialization of voxelwise corrections coefficients for time basis functions coefficients
1392  for (int v=0 ; v<mp_ID->GetNbVoxXYZ() ; v++)
1393  mp_corrBasisCoeffs[v] = 0;
1394 
1395  basis_functions_norm = 0.;
1396 
1397 
1398  // Compute normalization related to the temporal basis functions.
1399  for (int fr=0 ; fr<mp_ID->GetNbTimeFrames() ; fr++)
1400  for (int rg=0 ; rg<mp_ID->GetNbRespGates() ; rg++)
1401  for (int cg=0 ; cg<mp_ID->GetNbCardGates() ; cg++)
1402  for (int rb=0 ; rb<m_nbRgateBF ; rb++)
1403  for (int cb=0 ; cb<m_nbCgateBF ; cb++)
1404  basis_functions_norm += m2p_nestedModelTimeBasisFunctions[fb][fr]
1405  * m2p_respBasisFunctions[rb][rg]
1406  * m2p_cardBasisFunctions[cb][cg];
1407 
1408  int v;
1409  #pragma omp parallel for private(v) schedule(static, 1)
1410  for (v=0 ; v<mp_ID->GetNbVoxXYZ() ; v++)
1411  {
1412  // If a mask has been provided, check if the model applies in this voxel
1413  if(mp_maskModel != NULL && mp_maskModel[ v ] == 0)
1414  continue;
1415 
1416  for (int fr=0 ; fr<mp_ID->GetNbTimeFrames() ; fr++)
1417  for (int rg=0 ; rg<mp_ID->GetNbRespGates() ; rg++)
1418  for (int cg=0 ; cg<mp_ID->GetNbCardGates() ; cg++)
1419  for (int rb=0 ; rb<m_nbRgateBF ; rb++)
1420  for (int cb=0 ; cb<m_nbCgateBF ; cb++)
1421  {
1422  mp_corrBasisCoeffs[v] += ap_ImageS->m6p_backwardImage[0][0][fr][rg][cg][v]
1423  * m2p_nestedModelTimeBasisFunctions[fb][fr]
1424  * m2p_respBasisFunctions[rb][rg]
1425  * m2p_cardBasisFunctions[cb][cg];
1426  }
1427 
1428  // Apply corrections and normalization to the frame model parametric image.
1429  if (mp_corrBasisCoeffs[v] >= 0.) // could be == 0
1430  m2p_parametricImages[fb][v] *= mp_corrBasisCoeffs[v]/basis_functions_norm;
1431  }
1432  }
1433  }
1434 
1435  // Regularisation according to respiratory basis functions
1436 
1437  if( m_nbRGModelParam>1 )
1438  {
1439  if(m_verbose >=3) Cout("iLinearModel::NestedEM() -> Estimate Coefficients - Respiratory Gate coeffs estimation step" <<endl);
1440 
1441  for (int rb=0 ; rb<m_nbRgateBF ; rb++)
1442  {
1443  // Initialization of voxelwise corrections for time basis functions coefficients
1444  for (int v=0 ; v<mp_ID->GetNbVoxXYZ() ; v++)
1445  mp_corrBasisCoeffs[v] = 0;
1446 
1447  basis_functions_norm = 0.;
1448 
1449  for (int fr=0 ; fr<mp_ID->GetNbTimeFrames() ; fr++)
1450  for (int rg=0 ; rg<mp_ID->GetNbRespGates() ; rg++)
1451  for (int cg=0 ; cg<mp_ID->GetNbCardGates() ; cg++)
1452  for (int fb=0 ; fb<m_nbTimeBF ; fb++)
1453  for (int cb=0 ; cb<m_nbCgateBF ; cb++)
1454  basis_functions_norm += m2p_nestedModelTimeBasisFunctions[fb][fr]
1455  * m2p_respBasisFunctions[rb][rg]
1456  * m2p_cardBasisFunctions[cb][cg];
1457 
1458  // Compute corrections for the time vowelwise basis functions coefficients
1459  int v;
1460  #pragma omp parallel for private(v) schedule(static, 1)
1461  for (v=0 ; v<mp_ID->GetNbVoxXYZ() ; v++)
1462  {
1463  // If a mask has been provided, check if the model applies in this voxel
1464  if(mp_maskModel != NULL && mp_maskModel[ v ] == 0)
1465  continue;
1466 
1467  for (int fr=0 ; fr<mp_ID->GetNbTimeFrames() ; fr++)
1468  for (int rg=0 ; rg<mp_ID->GetNbRespGates() ; rg++)
1469  for (int cg=0 ; cg<mp_ID->GetNbCardGates() ; cg++)
1470  for (int fb=0 ; fb<m_nbTimeBF ; fb++)
1471  for (int cb=0 ; cb<m_nbCgateBF ; cb++)
1472  {
1473  mp_corrBasisCoeffs[v] += ap_ImageS->m6p_backwardImage[0][0][fr][rg][cg][v]
1474  * m2p_nestedModelTimeBasisFunctions[fb][fr]
1475  * m2p_respBasisFunctions[rb][rg]
1476  * m2p_cardBasisFunctions[cb][cg];
1477  }
1478 
1479  // Apply corrections and normalization to the respiratory model parametric image.
1480  if (mp_corrBasisCoeffs[v] >= 0.)
1481  m2p_RGParametricImages[rb][v] *= mp_corrBasisCoeffs[v]/basis_functions_norm;
1482  }
1483  }
1484  }
1485 
1486 
1487  // Regularisation according to cardiac basis functions
1488  if( m_nbCGModelParam>1 )
1489  {
1490  if(m_verbose >=3) Cout("iLinearModel::NestedEM() -> Estimate Coefficients - Cardiac Gate coeffs estimation step" <<endl);
1491 
1492  for (int cb=0 ; cb<m_nbCgateBF ; cb++)
1493  {
1494  // Initialization of voxelwise corrections for time basis functions coefficients
1495  for (int v=0 ; v<mp_ID->GetNbVoxXYZ() ; v++)
1496  mp_corrBasisCoeffs[v] = 0;
1497 
1498  basis_functions_norm = 0.;
1499 
1500  for (int fr=0 ; fr<mp_ID->GetNbTimeFrames() ; fr++)
1501  for (int rg=0 ; rg<mp_ID->GetNbRespGates() ; rg++)
1502  for (int cg=0 ; cg<mp_ID->GetNbCardGates() ; cg++)
1503  for (int fb=0 ; fb<m_nbTimeBF ; fb++)
1504  for (int rb=0 ; rb<m_nbRgateBF ; rb++)
1505  basis_functions_norm += m2p_nestedModelTimeBasisFunctions[fb][fr]
1506  * m2p_respBasisFunctions[rb][rg]
1507  * m2p_cardBasisFunctions[cb][cg];
1508 
1509  // Compute corrections for the time vowelwise basis functions coefficients
1510  int v;
1511  #pragma omp parallel for private(v) schedule(static, 1)
1512  for (v=0 ; v<mp_ID->GetNbVoxXYZ() ; v++)
1513  {
1514  // If a mask has been provided, check if the model applies in this voxel
1515  if(mp_maskModel != NULL && mp_maskModel[ v ] == 0)
1516  continue;
1517 
1518  for (int fr=0 ; fr<mp_ID->GetNbTimeFrames() ; fr++)
1519  for (int rg=0 ; rg<mp_ID->GetNbRespGates() ; rg++)
1520  for (int cg=0 ; cg<mp_ID->GetNbCardGates() ; cg++)
1521  for (int fb=0 ; fb<m_nbTimeBF ; fb++)
1522  for (int rb=0 ; rb<m_nbRgateBF ; rb++)
1523  {
1524  mp_corrBasisCoeffs[v] += ap_ImageS->m6p_backwardImage[0][0][fr][rg][cg][v]
1525  * m2p_nestedModelTimeBasisFunctions[fb][fr]
1526  * m2p_respBasisFunctions[rb][rg]
1527  * m2p_cardBasisFunctions[cb][cg];
1528  }
1529 
1530  // Apply corrections and normalization to the cardiac model parametric image.
1531  if (mp_corrBasisCoeffs[v] >= 0.)
1532  m2p_CGParametricImages[cb][v] *= mp_corrBasisCoeffs[v]/basis_functions_norm;
1533  }
1534  }
1535  }
1536 
1537  m_basisFunctionsUpdIdx++;
1538 
1539  // Some feedback :
1540  if(m_verbose >=3)
1541  {
1542  if( m_nbModelParam>1 )
1543  for (int fb=0 ; fb<m_nbTimeBF ; fb++)
1544  {
1545  FLTNB avg = 0;
1546  for (int v=0 ; v<mp_ID->GetNbVoxXYZ() ; v++)
1547  if(mp_maskModel != NULL && mp_maskModel[ v ] == 1)
1548  avg += m2p_parametricImages[fb][v];
1549 
1550  avg /= m_nbVoxelsMask;
1551  Cout( "iLinearModel::NestedEM() -> Frame parametric image["<<fb<<"] avg value:" << avg << endl);
1552  }
1553 
1554  if( m_nbRGModelParam>1 )
1555  for (int rb=0 ; rb<m_nbRgateBF ; rb++)
1556  {
1557  FLTNB avg = 0;
1558 
1559  for (int v=0 ; v<mp_ID->GetNbVoxXYZ() ; v++)
1560  if(mp_maskModel != NULL && mp_maskModel[ v ] == 1)
1561  avg += m2p_RGParametricImages[rb][v];
1562 
1563  avg /= m_nbVoxelsMask;
1564  Cout( "iLinearModel::NestedEM() -> Resp gate parametric image["<<rb<<"] avg value:" << avg << endl);
1565  }
1566 
1567  if( m_nbCGModelParam>1 )
1568  for (int cb=0 ; cb<m_nbCgateBF ; cb++)
1569  {
1570  FLTNB avg = 0;
1571 
1572  for (int v=0 ; v<mp_ID->GetNbVoxXYZ() ; v++)
1573  if(mp_maskModel != NULL && mp_maskModel[ v ] == 1)
1574  avg += m2p_CGParametricImages[cb][v];
1575 
1576  avg /= m_nbVoxelsMask;
1577  Cout( "iLinearModel::NestedEM() -> Card gate parametric image["<<cb<<"] avg value:" << avg << endl);
1578  }
1579  }
1580 
1581  } // end of else loop (parametric images update)
1582 
1583  } // end of loop on cycles (linear model iterations)
1584 
1585  return 0;
1586 }
1587 
1588 
1589 
1590 // =====================================================================
1591 // ---------------------------------------------------------------------
1592 // ---------------------------------------------------------------------
1593 // =====================================================================
1594 /*
1595  \fn EstimateParametersWithNNLS
1596  \param ap_ImageS : pointer to the ImageSpace
1597  \param a_ite : index of the actual iteration (not used)
1598  \brief Estimate parametric images using the NNLS method
1599  \return 0 if success, other value otherwise.
1600 */
1601 int iLinearModel::EstimateParametersWithNNLS(oImageSpace* ap_ImageS, int a_ite)
1602 {
1603  if(m_verbose >=3) Cout("iLinearModel::EstimateParametersWithNNLS() ..." <<endl);
1604 
1605  // Regularisation according to frame basis functions
1606  if( m_nbModelParam>1 )
1607  {
1608  if(m_verbose >=3) Cout("iLinearModel::EstimateParametersWithNNLS() -> Estimate Coefficients - Frame coeffs estimation step" <<endl);
1609 
1610  for (int rg=0 ; rg<mp_ID->GetNbRespGates() ; rg++)
1611  for (int cg=0 ; cg<mp_ID->GetNbCardGates() ; cg++)
1612  {
1613  int v;
1614  #pragma omp parallel for private(v) schedule(guided)
1615  for (v=0 ; v<mp_ID->GetNbVoxXYZ() ; v++)
1616  {
1617  // If a mask has been provided, check if the model applies in this voxel
1618  if(mp_maskModel != NULL && mp_maskModel[ v ] == 0)
1619  continue;
1620 
1621  int th=0;
1622  #ifdef CASTOR_OMP
1623  th = omp_get_thread_num();
1624  #endif
1625 
1626  FLTNB** pp_nnls_A = m3p_nnlsA[ th ];
1627  FLTNB* p_nnls_B = m2p_nnlsB[ th ];
1628  FLTNB nnls_rnorm;
1629  FLTNB *p_nnls_zz=NULL;
1630 
1631  // Fill NNLS matrices
1632  for(int t=0; t<mp_ID->GetNbTimeFrames(); t++)
1633  {
1634  // Fill NNLS A matrix:
1635  for (int fb=0 ; fb<m_nbTimeBF ; fb++)
1636  pp_nnls_A[ fb ][ t ] = m2p_nestedModelTimeBasisFunctions[ fb ][ t ] * mp_w[ t ];
1637 
1638  // Fill NNLS B array: tissue
1639  p_nnls_B[ t ]=ap_ImageS->m4p_image[ t ][ rg ][ cg ][ v ] * mp_w[ t ];
1640  }
1641 
1642  // Launch NNLS
1643  if(NNLS(pp_nnls_A,
1644  mp_ID->GetNbTimeFrames(),
1645  m_nbTimeBF,
1646  p_nnls_B,
1647  m2p_nnlsX[ th ],
1648  &nnls_rnorm,
1649  m2p_nnlsWp[ th ],
1650  p_nnls_zz,
1651  m2p_nnlsIdx[ th ]) )
1652  continue;
1653 
1654  // Recover parameters
1655  for (int fb=0 ; fb<m_nbTimeBF ; fb++)
1656  m2p_parametricImages[fb][v] = m2p_nnlsX[ th ][ fb ];
1657 
1658  } // end of voxel loop
1659 
1660  }
1661  }
1662 
1663  // Regularisation according to respiratory basis functions
1664 
1665  if( m_nbRGModelParam>1 )
1666  {
1667  if(m_verbose >=3) Cout("iLinearModel::EstimateParametersWithNNLS() -> Estimate Coefficients - Respiratory Gate coeffs estimation step" <<endl);
1668 
1669  for (int fr=0 ; fr<mp_ID->GetNbTimeFrames() ; fr++)
1670  for (int cg=0 ; cg<mp_ID->GetNbCardGates() ; cg++)
1671  {
1672  int v;
1673  #pragma omp parallel for private(v) schedule(guided)
1674  for (v=0 ; v<mp_ID->GetNbVoxXYZ() ; v++)
1675  {
1676  // If a mask has been provided, check if the model applies in this voxel
1677  if(mp_maskModel != NULL && mp_maskModel[ v ] == 0)
1678  continue;
1679 
1680  int th=0;
1681  #ifdef CASTOR_OMP
1682  th = omp_get_thread_num();
1683  #endif
1684 
1685  FLTNB** pp_nnls_A = m3p_nnlsA[ th ];
1686  FLTNB* p_nnls_B = m2p_nnlsB[ th ];
1687  FLTNB nnls_rnorm;
1688  FLTNB *p_nnls_zz=NULL;
1689 
1690  // Fill NNLS matrices
1691  for (int rg=0 ; rg<mp_ID->GetNbRespGates() ; rg++)
1692  {
1693  // Fill NNLS A matrix:
1694  for (int rb=0 ; rb<m_nbRgateBF ; rb++)
1695  pp_nnls_A[ rb ][ rg ] = m2p_respBasisFunctions[ rb ][ rg ];
1696 
1697  // Fill NNLS B array: tissue
1698  p_nnls_B[ rg ]=ap_ImageS->m4p_image[ fr ][ rg ][ cg ][ v ];
1699  }
1700 
1701  // Launch NNLS
1702  if(NNLS(pp_nnls_A,
1703  mp_ID->GetNbRespGates(),
1704  m_nbRgateBF,
1705  p_nnls_B,
1706  m2p_nnlsX[ th ],
1707  &nnls_rnorm,
1708  m2p_nnlsWp[ th ],
1709  p_nnls_zz,
1710  m2p_nnlsIdx[ th ]) )
1711  continue;
1712 
1713  // Recover parameters
1714  for (int rb=0 ; rb<m_nbRgateBF ; rb++)
1715  m2p_RGParametricImages[rb][v] = m2p_nnlsX[ th ][ rb ];
1716 
1717  } // end of loop on voxels
1718 
1719  }
1720  }
1721 
1722 
1723  // Regularisation according to cardiac basis functions
1724  if( m_nbCGModelParam>1 )
1725  {
1726  if(m_verbose >=3) Cout("iLinearModel::EstimateParametersWithNNLS() -> Estimate Coefficients - Cardiac Gate coeffs estimation step" <<endl);
1727 
1728  for (int fr=0 ; fr<mp_ID->GetNbTimeFrames() ; fr++)
1729  for (int rg=0 ; rg<mp_ID->GetNbRespGates() ; rg++)
1730  {
1731  int v;
1732  #pragma omp parallel for private(v) schedule(guided)
1733  for (v=0 ; v<mp_ID->GetNbVoxXYZ() ; v++)
1734  {
1735  // If a mask has been provided, check if the model applies in this voxel
1736  if(mp_maskModel != NULL && mp_maskModel[ v ] == 0)
1737  continue;
1738 
1739  int th=0;
1740  #ifdef CASTOR_OMP
1741  th = omp_get_thread_num();
1742  #endif
1743 
1744  FLTNB** pp_nnls_A = m3p_nnlsA[ th ];
1745  FLTNB* p_nnls_B = m2p_nnlsB[ th ];
1746  FLTNB nnls_rnorm;
1747  FLTNB *p_nnls_zz=NULL;
1748 
1749  // Fill NNLS matrices
1750  for (int cg=0 ; cg<mp_ID->GetNbCardGates() ; cg++)
1751  {
1752  // Fill NNLS A matrix:
1753  for (int cb=0 ; cb<m_nbCgateBF ; cb++)
1754  pp_nnls_A[ cb ][ cg ] = m2p_cardBasisFunctions[ cb ][ cg ];
1755 
1756  // Fill NNLS B array: tissue
1757  p_nnls_B[ cg ]=ap_ImageS->m4p_image[ fr ][ rg ][ cg ][ v ];
1758  }
1759 
1760  // Launch NNLS
1761  if(NNLS(pp_nnls_A,
1762  mp_ID->GetNbCardGates(),
1763  m_nbCgateBF,
1764  p_nnls_B,
1765  m2p_nnlsX[ th ],
1766  &nnls_rnorm,
1767  m2p_nnlsWp[ th ],
1768  p_nnls_zz,
1769  m2p_nnlsIdx[ th ]) )
1770  continue;
1771 
1772  // Recover parameters
1773  for (int cb=0 ; cb<m_nbCgateBF ; cb++)
1774  m2p_CGParametricImages[cb][v] = m2p_nnlsX[ th ][ cb ];
1775 
1776  } // end of loop on voxels
1777 
1778  }
1779  }
1780 
1781 
1782  // Some feedback :
1783  if(m_verbose >=3)
1784  {
1785  if( m_nbModelParam>1 )
1786  for (int fb=0 ; fb<m_nbTimeBF ; fb++)
1787  {
1788  FLTNB avg = 0;
1789  for (int v=0 ; v<mp_ID->GetNbVoxXYZ() ; v++)
1790  if(mp_maskModel != NULL && mp_maskModel[ v ] == 1)
1791  avg += m2p_parametricImages[fb][v];
1792 
1793  avg /= m_nbVoxelsMask;
1794  Cout( "iLinearModel::EstimateParametersWithNNLS() -> Frame parametric image["<<fb<<"] avg value:" << avg << endl);
1795  }
1796 
1797  if( m_nbRGModelParam>1 )
1798  for (int rb=0 ; rb<m_nbRgateBF ; rb++)
1799  {
1800  FLTNB avg = 0;
1801 
1802  for (int v=0 ; v<mp_ID->GetNbVoxXYZ() ; v++)
1803  if(mp_maskModel != NULL && mp_maskModel[ v ] == 1)
1804  avg += m2p_RGParametricImages[rb][v];
1805 
1806  avg /= m_nbVoxelsMask;
1807  Cout( "iLinearModel::EstimateParametersWithNNLS() -> Resp gate parametric image["<<rb<<"] avg value:" << avg << endl);
1808  }
1809 
1810  if( m_nbCGModelParam>1 )
1811  for (int cb=0 ; cb<m_nbCgateBF ; cb++)
1812  {
1813  FLTNB avg = 0;
1814 
1815  for (int v=0 ; v<mp_ID->GetNbVoxXYZ() ; v++)
1816  if(mp_maskModel != NULL && mp_maskModel[ v ] == 1)
1817  avg += m2p_CGParametricImages[cb][v];
1818 
1819  avg /= m_nbVoxelsMask;
1820  Cout( "iLinearModel::EstimateParametersWithNNLS() -> Card gate parametric image["<<cb<<"] avg value:" << avg << endl);
1821  }
1822  }
1823 
1824  return 0;
1825 }
1826 
1827 
1828 // =====================================================================
1829 // ---------------------------------------------------------------------
1830 // ---------------------------------------------------------------------
1831 // =====================================================================
1832 /*
1833  \fn EstimateParametersWithLS (linear regresion)
1834  \param ap_ImageS : pointer to the ImageSpace
1835  \param a_ite : index of the actual iteration (not used)
1836  \brief Estimate parametric images using LS linear regression
1837  \return 0 if success, other value otherwise.
1838 */
1839 
1840 int iLinearModel::Patlak_LS(oImageSpace* ap_ImageS, int a_ite)
1841 {
1842  if(m_verbose >=3) Cout("iLinearPatlakModel::Patlak_LS ..." <<endl);
1843 
1844  int v;
1845  #pragma omp parallel for private(v) schedule(guided)
1846  for (v=0 ; v<mp_ID->GetNbVoxXYZ() ; v++)
1847  {
1848  // If a mask has been provided, check if the model applies in this voxel
1849  if(mp_maskModel != NULL && mp_maskModel[ v ] == 0)
1850  continue;
1851 
1852  FLTNB xmean = 0.,
1853  ymean = 0.,
1854  K = 0.,
1855  Vd = 0.,
1856  Cov= 0.,
1857  Var = 0.;
1858 
1859 
1860  // Compute means for this voxel
1861  for (int fr=0 ; fr<mp_ID->GetNbTimeFrames() ; fr++)
1862  {
1863  xmean += m2p_nestedModelTimeBasisFunctions[0][fr]/m2p_nestedModelTimeBasisFunctions[1][fr];
1864  ymean += ap_ImageS->m4p_image[fr][0][0][v]/m2p_nestedModelTimeBasisFunctions[1][fr];
1865  }
1866 
1867  xmean /= mp_ID->GetNbTimeFrames();
1868  ymean /= mp_ID->GetNbTimeFrames();
1869 
1870  FLTNB Yt,Xt ;
1871 
1872  // Compute covariance & variance
1873  for (int fr=0 ; fr<mp_ID->GetNbTimeFrames() ; fr++)
1874  {
1875  Yt = ap_ImageS->m4p_image[fr][0][0][v]/m2p_nestedModelTimeBasisFunctions[1][fr];
1876  Xt = m2p_nestedModelTimeBasisFunctions[0][fr]/m2p_nestedModelTimeBasisFunctions[1][fr];
1877  Cov += (Xt-xmean) * (Yt-ymean);
1878  Var += (Xt-xmean) * (Xt-xmean);
1879  }
1880 
1881  // Slope
1882  K = (Var != 0) ? Cov/Var : 0.;
1883 
1884  // Non-negativity constraint
1885  if(K<0) K = 0.;
1886 
1887  // Intercept
1888  Vd = ymean - K*xmean;
1889 
1890  // Non-negativity constraint
1891  if(Vd<0) Vd = 0.;
1892 
1893  m2p_parametricImages[0][v] = K;
1894  m2p_parametricImages[1][v] = Vd;
1895  }
1896 
1897  return 0;
1898 }
1899 
1900 
1901 // =====================================================================
1902 // ---------------------------------------------------------------------
1903 // ---------------------------------------------------------------------
1904 // =====================================================================
1905 /*
1906  \fn EstimateImageWithModel
1907  \param ap_ImageS : pointer to the ImageSpace
1908  \param a_ite : index of the actual iteration (not used)
1909  \param a_sset : index of the actual subset (not used)
1910  \brief Re-estimate image using the linear parametric images and basis functions
1911  \return 0 if success, other value otherwise.
1912 */
1913 int iLinearModel::EstimateImageWithModel(oImageSpace* ap_ImageS, int a_ite, int a_sset)
1914 {
1915  if(m_verbose >= 3) Cout("iLinearModel::EstimateImageWithModel ... " <<endl);
1916 
1917  #ifdef CASTOR_DEBUG
1918  if (!m_initialized)
1919  {
1920  Cerr("***** iLinearModel::EstimateImageWithModel() -> Called while not initialized !" << endl);
1921  Exit(EXIT_DEBUG);
1922  }
1923  #endif
1924 
1925  int v;
1926  #pragma omp parallel for private(v) schedule(static, 1)
1927  for (int fr=0 ; fr<mp_ID->GetNbTimeFrames() ; fr++)
1928  for (int rg=0 ; rg<mp_ID->GetNbRespGates() ; rg++)
1929  for (int cg=0 ; cg<mp_ID->GetNbCardGates() ; cg++)
1930  for (v=0 ; v<mp_ID->GetNbVoxXYZ() ; v++)
1931  {
1932  // If a mask has been provided, check if the model applies in this voxel
1933  if(mp_maskModel != NULL && mp_maskModel[ v ] == 0)
1934  continue;
1935 
1936  // Reset current estimated value
1937  ap_ImageS->m4p_image[fr][rg][cg][v] = 0.;
1938 
1939  for (int fb=0 ; fb<m_nbTimeBF ; fb++)
1940  for (int rb=0 ; rb<m_nbRgateBF ; rb++)
1941  for (int cb=0 ; cb<m_nbCgateBF ; cb++)
1942  {
1943  // Retrieve current estimation of image according to coeffs/basis functions
1944  ap_ImageS->m4p_image[fr][rg][cg][v] += m2p_parametricImages[fb][v]
1945  * m2p_nestedModelTimeBasisFunctions[fb][fr]
1946  * m2p_RGParametricImages[rb][v]
1947  * m2p_respBasisFunctions[rb][rg]
1948  * m2p_CGParametricImages[cb][v]
1949  * m2p_cardBasisFunctions[cb][cg] ;
1950  }
1951  }
1952 
1953  return 0;
1954 }
INTF_LERP_DISABLED
#define INTF_LERP_DISABLED
Definition: oInterfileIO.hh:96
vDynamicModel::NNLS
int NNLS(FLTNB **A, int m, int n, FLTNB *B, FLTNB *X, FLTNB *rnorm, FLTNB *wp, FLTNB *zzp, int *indexp)
Implementation of NNLS (non-negative least squares) algorithm Derived from Turku PET center libraries...
Definition: vDynamicModel.cc:836
iLinearModel::m_OptimisationMethod
int m_OptimisationMethod
Definition: iLinearModel.hh:195
FLTNB
#define FLTNB
Definition: gVariables.hh:81
vDynamicModel::m_listOptions
string m_listOptions
Definition: vDynamicModel.hh:407
vDynamicModel::m_nnlsN
uint16_t m_nnlsN
Definition: vDynamicModel.hh:435
iLinearModel.hh
Declaration of class iLinearModel.
iLinearModel::EstimateImageWithModel
int EstimateImageWithModel(oImageSpace *ap_Image, int a_ite, int a_sset)
Re-estimate image using the linear parametric images and basis functions.
Definition: iLinearModel.cc:1913
vDynamicModel::mp_ID
oImageDimensionsAndQuantification * mp_ID
Definition: vDynamicModel.hh:362
vDynamicModel
This is the mother class of dynamic model classes.
Definition: vDynamicModel.hh:62
vDynamicModel::m2p_nnlsB
FLTNB ** m2p_nnlsB
Definition: vDynamicModel.hh:429
iLinearModel::mp_corrBasisCoeffs
FLTNB * mp_corrBasisCoeffs
Definition: iLinearModel.hh:210
OPTIMISATION_METHOD_DR
#define OPTIMISATION_METHOD_DR
Definition: iLinearModel.hh:42
iLinearModel::InitializeSpecificToAllLinearModels
int InitializeSpecificToAllLinearModels()
This function is used to initialize the parametric images and basis functions for all Linear Models...
Definition: iLinearModel.cc:552
vDynamicModel::m2p_RGParametricImages
FLTNB ** m2p_RGParametricImages
Definition: vDynamicModel.hh:396
Exit
void Exit(int code)
Definition: sOutputManager.cc:46
vDynamicModel::m_fileOptions
string m_fileOptions
Definition: vDynamicModel.hh:406
iLinearModel::ReadAndCheckOptionsList
int ReadAndCheckOptionsList(string a_listOptions)
This function is used to read parameters from a string.
Definition: iLinearModel.cc:431
iLinearModel::iLinearModel
iLinearModel()
Constructor of iLinearModel. Simply set all data members to default values.
Definition: iLinearModel.cc:41
iLinearModel::ShowBasisFunctions
void ShowBasisFunctions()
This function is used to print the basis functions.
Definition: iLinearModel.cc:951
Cerr
#define Cerr(MESSAGE)
Definition: sOutputManager.hh:112
oImageSpace::m6p_backwardImage
FLTNB ****** m6p_backwardImage
Definition: oImageSpace.hh:94
oImageSpace::m4p_image
FLTNB **** m4p_image
Definition: oImageSpace.hh:80
iLinearModel::mp_corrBasisFunctions
FLTNB * mp_corrBasisFunctions
Definition: iLinearModel.hh:211
vDynamicModel::m2p_outputParImages
FLTNB ** m2p_outputParImages
Definition: vDynamicModel.hh:389
iLinearModel::CheckSpecificParametersForAllLinearModels
int CheckSpecificParametersForAllLinearModels()
This function is used to check parameters for all Linear Models. .
Definition: iLinearModel.cc:487
vDynamicModel::m2p_CGParametricImages
FLTNB ** m2p_CGParametricImages
Definition: vDynamicModel.hh:401
ReadDataASCIIFile
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:129
iLinearModel::ReadAndCheckConfigurationFileSpecific
virtual int ReadAndCheckConfigurationFileSpecific()
This function is used to read options from a configuration file.
Definition: iLinearModel.cc:252
iLinearModel::EstimateParametersWithNNLS
int EstimateParametersWithNNLS(oImageSpace *ap_ImageS, int a_ite)
Estimate parametric images using the NNLS method.
Definition: iLinearModel.cc:1601
vDynamicModel::mp_maskModel
FLTNB * mp_maskModel
Definition: vDynamicModel.hh:423
vDynamicModel::m2p_nnlsWp
FLTNB ** m2p_nnlsWp
Definition: vDynamicModel.hh:432
vDynamicModel::m_ModelSpecificBasisFunctionsRequired
bool m_ModelSpecificBasisFunctionsRequired
Definition: vDynamicModel.hh:414
vDynamicModel::m2p_nestedModelTimeBasisFunctions
FLTNB ** m2p_nestedModelTimeBasisFunctions
Definition: vDynamicModel.hh:383
KEYWORD_MANDATORY
#define KEYWORD_MANDATORY
Definition: gOptions.hh:47
iLinearModel::m_nbLinearModelCycles
uint32_t m_nbLinearModelCycles
Definition: iLinearModel.hh:212
KEYWORD_OPTIONAL
#define KEYWORD_OPTIONAL
Definition: gOptions.hh:49
OPTIMISATION_METHOD_LS
#define OPTIMISATION_METHOD_LS
Definition: iLinearModel.hh:48
iLinearModel::InitializeSpecific
int InitializeSpecific()
This function is used to initialize the parametric images and basis functions for all Linear Models...
Definition: iLinearModel.cc:983
iLinearModel::m2p_cardBasisFunctions
FLTNB ** m2p_cardBasisFunctions
Definition: iLinearModel.hh:205
iLinearModel::m_basisFunctionsUpdIdx
int m_basisFunctionsUpdIdx
Definition: iLinearModel.hh:217
oImageDimensionsAndQuantification::GetNbCardGates
int GetNbCardGates()
Get the number of cardiac gates.
Definition: oImageDimensionsAndQuantification.hh:699
vDynamicModel::m2p_parametricImages
FLTNB ** m2p_parametricImages
Definition: vDynamicModel.hh:374
oImageSpace
This class holds all the matrices in the image domain that can be used in the algorithm: image...
Definition: oImageSpace.hh:60
iLinearModel::CheckSpecificParameters
int CheckSpecificParameters()
This function is used to check whether all member variables have been correctly initialized or not...
Definition: iLinearModel.cc:460
oImageDimensionsAndQuantification::GetNbThreadsForImageComputation
int GetNbThreadsForImageComputation()
Get the number of threads used for image operations.
Definition: oImageDimensionsAndQuantification.hh:403
vDynamicModel::m2p_nnlsX
FLTNB ** m2p_nnlsX
Definition: vDynamicModel.hh:431
oImageDimensionsAndQuantification::GetNbTimeFrames
int GetNbTimeFrames()
Get the number of time frames.
Definition: oImageDimensionsAndQuantification.hh:477
oImageDimensionsAndQuantification::GetNbVoxXYZ
INTNB GetNbVoxXYZ()
Get the total number of voxels.
Definition: oImageDimensionsAndQuantification.hh:778
iLinearModel::~iLinearModel
~iLinearModel()
Destructor of iLinearModel.
Definition: iLinearModel.cc:76
iLinearModel::Patlak_LS
int Patlak_LS(oImageSpace *ap_ImageS, int a_ite)
Estimate parametric images using linear regression.
Definition: iLinearModel.cc:1840
EXIT_DEBUG
#define EXIT_DEBUG
Definition: gVariables.hh:97
OPTIMISATION_METHOD_NNLS
#define OPTIMISATION_METHOD_NNLS
Definition: iLinearModel.hh:46
oImageDimensionsAndQuantification::GetNbThreadsForProjection
int GetNbThreadsForProjection()
Get the number of threads used for projections.
Definition: oImageDimensionsAndQuantification.hh:396
OPTIMISATION_METHOD_NESTEM
#define OPTIMISATION_METHOD_NESTEM
Definition: iLinearModel.hh:44
iLinearModel::m_basisFunctionsUpdRatio
int m_basisFunctionsUpdRatio
Definition: iLinearModel.hh:215
vDynamicModel::m2p_nnlsMat
FLTNB ** m2p_nnlsMat
Definition: vDynamicModel.hh:430
vDynamicModel::m_noImageUpdateFlag
bool m_noImageUpdateFlag
Definition: vDynamicModel.hh:416
oImageDimensionsAndQuantification::GetNbRespGates
int GetNbRespGates()
Get the number of respiratory gates.
Definition: oImageDimensionsAndQuantification.hh:648
Cout
#define Cout(MESSAGE)
Definition: sOutputManager.hh:100
iLinearModel::m2p_respBasisFunctions
FLTNB ** m2p_respBasisFunctions
Definition: iLinearModel.hh:201
IntfReadImgDynCoeffFile
int IntfReadImgDynCoeffFile(const string &a_pathToHeaderFile, FLTNB **a2p_ImgMatrix, oImageDimensionsAndQuantification *ap_ID, int a_nbFbases, int vb, bool a_lerpFlag)
Function dedicated to Interfile image reading for dynamic coefficients images.
Definition: oInterfileIO.cc:733
iLinearModel::NestedEM
int NestedEM(oImageSpace *ap_ImageS, int a_ite)
Estimate parametric images and basis functions (if enabled) using the nested EM method.
Definition: iLinearModel.cc:1104
iLinearModel::m_basisFunctionsUpdStartIte
int m_basisFunctionsUpdStartIte
Definition: iLinearModel.hh:213
iLinearModel::ReadAndCheckConfigurationFileSpecificToAllLinearModels
int ReadAndCheckConfigurationFileSpecificToAllLinearModels()
This function is used to read parameters that are generic for all Linear Models. ...
Definition: iLinearModel.cc:276
vDynamicModel::mp_blackListedvoxelsImage
FLTNB * mp_blackListedvoxelsImage
Definition: vDynamicModel.hh:379
iLinearModel::EstimateModelParameters
int EstimateModelParameters(oImageSpace *ap_Image, int a_ite, int a_sset)
Estimate model parameters (parametric images and basis functions)
Definition: iLinearModel.cc:1022
vDynamicModel::m3p_nnlsA
FLTNB *** m3p_nnlsA
Definition: vDynamicModel.hh:428
iLinearModel::ShowHelp
void ShowHelp()
Print out specific help about the implementation of this model and its initialization.
Definition: iLinearModel.cc:139
Generated on Fri Nov 15 2019 15:00:59 for CASToR by   doxygen 1.8.13