151 if(
mp_Y[ th ] )
delete mp_Y[ th ] ;
152 if(
mp_X[ th ] )
delete mp_X[ th ] ;
203 cout <<
"-- This class implements a 2 compartments kinetic model " << endl;
204 cout <<
"-- or 1 Tissue Compartment Model (1TCM) for perfusion quantitation " << endl;
205 cout <<
"-- for radiotracers such as [15-O]H2O. " << endl;
207 cout <<
"-- *--------* K1 *--------* " << endl;
208 cout <<
"-- | | -----> | | " << endl;
209 cout <<
"-- | Cp | | Ct | " << endl;
210 cout <<
"-- | | <----- | | " << endl;
211 cout <<
"-- *--------* k2 *--------* " << endl;
213 cout <<
"-- This model contains 3 parameters, including 2 rate constants, " << endl;
214 cout <<
"-- K1 (v/s/v, where v is a volume unit), k2 (s-1) and the arterial volume fraction in tissue (Va), " << endl;
215 cout <<
"-- as described by the following equations : " << endl;
216 cout <<
"-- (1) Cpet ( t ) = Ct( t ) + Va*Cp( t ) " << endl;
217 cout <<
"-- (2) dCt( t ) / dt = K1*Cp( t ) - k2*Ct( t ) " << endl;
218 cout <<
"-- where Cpet( t ) is the image value in voxel for frame t, and" << endl;
219 cout <<
"-- Ct and Cp are activity concentration in tissue and plasma," << endl;
220 cout <<
"-- The model estimates K1, k2 and Va parametric images." << endl;
221 cout <<
"-- The input function (Cp) must be provided by the user." << endl;
223 cout <<
"--------------------------------------" << endl;
224 cout <<
" The model can be initialized using either a configuration text file, or a list of options :" << endl;
226 cout <<
" - CONFIGURATION FILE : (command-line options : _1TCM:path/to/conf/file.txt:)" << endl;
228 cout <<
" 'Input_function:' (mandatory) " << endl;
229 cout <<
" Enter the activity concentration values (kBq/v/s) of the plasma function (cp) for each time frame (tf)," << endl;
230 cout <<
" separated by ','. Time unit is seconds : " << endl;
231 cout <<
" -> Input functions: " << endl;
232 cout <<
" -> val_cp_tf1 , val_cp_tf2 , ..., val_cp_tfn" << endl;
234 cout <<
" 'Integral_input_function:' (optional) " << endl;
235 cout <<
" Enter the activity concentration values of the integral of the plasma function (Icp) for each time frame (tf)," << endl;
236 cout <<
" separated by ','. Time unit is seconds :" << endl;
237 cout <<
" -> Integral_input_function: " << endl;
238 cout <<
" -> val_Icp_tf1, val_Icp_tf2, ..., val_Icp_tfn" << endl;
239 cout <<
" NOTE: If not provided, the integral will be estimated from the input function" << endl;
240 cout <<
" Set the 'Integral method' flag below to select the method for TAC integration (default: WPO))" << endl;
242 cout <<
" 'Optimisation_method:' (optional) Define the method to use for parameters estimation. Only least-squares methods " << endl;
243 cout <<
" (default: NNLS) are currently available to estimate O = (K1, k2, Va)" << endl;
244 cout <<
" Ô = [ X'WX ]-1 X'Wy, with " << endl;
245 cout <<
" with y = data vector " << endl;
246 cout <<
" X = model matrix (Cp, Icp, Cpet) " << endl;
247 cout <<
" W = weights vector (frame duration) " << endl;
248 cout <<
" If the estimated parameter values are negative, the activity value for the related voxels will be kept to their original number" << endl;
249 cout <<
" 0 (default) = Non-negative least-square (NNLS)" << endl;
250 cout <<
" Derived from Turku PET center libraries, authors: Vesa Oikonen and Kaisa Sederholm" << endl;
251 cout <<
" (http://www.turkupetcentre.net/petanalysis/index.html)" << endl;
252 cout <<
" Routine based on the text and fortran code in C.L. Lawson and R.J. Hanson," << endl;
253 cout <<
" Solving Least Squares Problems, Prentice-Hall, Englewood Cliffs, New Jersey, 1974." << endl;
254 cout <<
" Note: K1 estimated values could still be negative as they are computed from a substraction of the estimated NNLS parameters." << endl;
255 cout <<
" Activity values will be kept to their original values for the voxels involved" << endl;
256 cout <<
" 1 = Standard Least Square (LS) optimization " << endl;
258 cout <<
" 'Integration_method:' (optional) Define the method to use for TAC integration over the time samples" << endl;
259 cout <<
" 0 (default) = Weighed parabola overlapping (WPO) (Z.Wang, D.Feng, Int. Sys. Sci 23 (1992), pp.1361-69)" << endl;
260 cout <<
" 1 = Trapezoidal " << endl;
261 cout <<
" 'Ridge_parameter:' (optional) Constant for Ridge Regression during Least-Square optimization (only available with Least-Square algorithm and not NNLS ) " << endl;
262 cout <<
" (default: 0) Bounds must be provided with the eponym options below in order to compute ridge weights and means for the new cost function: " << endl;
263 cout <<
" Ô = [ X'*W*X + t.Rw ]-1 [ X'*W*y ] " << endl;
264 cout <<
" + [ X'*W*X + t.Rw ]-1 [ t.Rw*Rm ] " << endl;
265 cout <<
" with y = data vector " << endl;
266 cout <<
" X = model matrix (Cp, Icp, Cpet) " << endl;
267 cout <<
" W = weights vector (frame duration) " << endl;
268 cout <<
" t = Ridge constant " << endl;
269 cout <<
" Rw= Ridge weights " << endl;
270 cout <<
" Rm= Ridge means " << endl;
271 cout <<
" 'Bounds:' (optional) Upper / Lower Bounds for each 3 parameters, to define ridge means mr0, and weights wr0," << endl;
272 cout <<
" such as mr0 = (Max+Min)/2 and wr0 = 1 / (Max-Min)^2 " << endl;
273 cout <<
" They must be provided as in the following syntax: " << endl;
274 cout <<
" Bounds: K1Max, K1Min, k2Max, k2Min, VaMax, VaMin" << endl;
275 cout <<
" Default: K1(Max,Min) = 0.1, 0. " << endl;
276 cout <<
" K2(Max,Min) = 0.1, 0. " << endl;
277 cout <<
" Va(Max,Min) = 1. , 0. " << endl;
278 cout <<
" 'VA_image:' (optional) Path to an interfile image containing the arterial volume fraction value in tissue for each voxel," << endl;
279 cout <<
" only K1 and k2 rate constants will be estimated (Default: All parameters are estimated) " << endl;
281 cout <<
" The following keywords are common to all dynamic models :" << endl;
282 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;
283 cout <<
" (Default x == 0) " << endl;
284 cout <<
" 'No image update: x' If set to 1, the reconstructed images for the next iteration/subset are not reestimated using the model" << endl;
285 cout <<
" (Default x == 0) (the code just performs standard independent reconstruction of each frames/gates) " << endl;
286 cout <<
" 'No parameters update: x' If set to 1, the parameters / functions of the model are not estimated with the image" << endl;
287 cout <<
" (Default x == 0) (this could be used to test The EstimateImageWithModel() function with specific user-provided parametric images) " << endl;
288 cout <<
" 'Save parametric images : x' Enable (1)/Disable(0) saving parametric images on disk" << endl;
289 cout <<
" (Default x == 1) " << endl;
290 cout <<
" 'Save blacklisted voxels images : x' Enable (1)/Disable(0) saving blacklisted voxels images on disk" << endl;
291 cout <<
" (Default x == 0) " << endl;
292 cout <<
" 'Mask:' Path to an interfile image containing a mask indicating if the model must be applied (1) or not (0) in each voxel" << endl;
293 cout <<
" (Default: model applies everywhere) " << endl;
296 cout <<
" - LINE OPTIONS : (command-line options : _1TCM,options:)" << endl;
298 cout <<
" Command-line options just require the samples of the plasma input curve, separated by commas. " << endl;
299 cout <<
" All other options will be set to default. The optimization algorithm will be NNLS, the integration method for TAC will be WPO " << endl;
300 cout <<
" -> 1cptModel,val_cp_tf1 , val_cp_tf2 , ..., val_cp_tfn" << endl;
318 if(
m_verbose >=3)
Cout(
"i1TCModel::ReadAndCheckConfigurationFileSpecific ..."<< endl);
331 string dtag =
"Integration_method";
338 Cerr(
"***** i1TCModel::ReadAndCheckConfigurationFileSpecific() -> Error while trying to read '"<< dtag <<
"' flag in "<<
m_fileOptions << endl);
343 dtag =
"Optimisation_method";
350 Cerr(
"***** i1TCModel::ReadAndCheckConfigurationFileSpecific() -> Error while trying to read '"<< dtag <<
"' flag in "<<
m_fileOptions << endl);
355 dtag =
"Save_blacklisted_voxels_images";
362 Cerr(
"***** vDynamicModel::ReadAndCheckConfigurationFile() -> Error while trying to read '"<< dtag <<
"' flag in " <<
m_fileOptions << endl);
367 const int nb_params = 3;
369 HPFLTNB p_bounds_params[ 6 ] = {0.1, 0., 0.1, 0., 1., 0.};
372 dtag =
"Ridge_Parameter";
379 Cerr(
"***** i1TCModel::ReadAndCheckConfigurationFileSpecific() -> Error while trying to read '"<< dtag <<
"' flag in "<<
m_fileOptions << endl);
394 Cerr(
"***** i1TCModel::ReadAndCheckConfigurationFileSpecific() -> Error while trying to read '"<< dtag <<
"' flag (upper/lower bounds) in "<<
m_fileOptions << endl);
402 for (
int p=0 ; p<nb_params ; p++ )
410 Cerr(
"***** i1TCModel::ReadAndCheckConfigurationFileSpecific() -> Error while trying to read configuration file at: " <<
m_fileOptions << endl);
432 if(
m_verbose >=2)
Cout(
"i1TCModel::ReadAndCheckOptionsList ..."<< endl);
457 if(
m_verbose >=2)
Cout(
"i1TCModel::CheckSpecificParameters ..."<< endl);
462 Cerr(
"***** i1TCModel::CheckSpecificParameters() -> ImageDimensions object has not been provided !" << endl);
469 Cerr(
"***** i1TCModel::CheckSpecificParameters() -> Wrong number of time frame basis functions !" << endl);
476 Cerr(
"***** i1TCModel::CheckSpecificParameters() -> Either a file or a list of options have to be selected to initialize the model, but not both ! " << endl);
483 Cerr(
"***** i1TCModel::CheckSpecificParameters -> Either a file or a list of options should have been provided at this point ! " << endl);
490 Cerr(
"***** i1TCModel::CheckSpecificParameters -> The implemented model is not compatible yet with gated reconstruction (parametric images will be the same for each gate)! " << endl);
498 Cerr(
"***** i1TCModel::CheckSpecificParameters() -> Error, ridge regression is not compatible with NNLS algorithm. Switch 'Optimisation method' flag to LS in order to use ridge-regression !" << endl);
520 if(
m_verbose >=2)
Cout(
"i1TCModel::InitializeSpecific() ..."<< endl);
525 Cerr(
"***** oDynamicModelManager::InitializeSpecific() -> Must call CheckParameters functions before Initialize() !" << endl);
566 string dtag =
"Input_function";
573 Cerr(
"***** i1TCModel::InitializeSpecific() -> Error while trying to read 1TCPM input functions !" << endl);
579 dtag =
"Integral_input_function";
586 Cerr(
"***** i1TCModel::InitializeSpecific() -> Error while trying to read 1TCPM input functions !" << endl);
593 string input_image =
"";
594 int return_value = 0;
596 dtag =
"Parametric_images_init";
603 if( return_value == 0)
613 Cerr(
"***** i1TCModel::InitializeSpecific() -> Error while trying to read the provided initialization parametric images : " << input_image << endl);
617 else if( return_value == 1)
619 Cerr(
"***** i1TCModel::InitializeSpecific() -> Error while trying to read input functions coefficients !" << endl);
632 string path_to_VAimage;
640 Cerr(
"***** i1TCModel::InitializeSpecific() -> Error while trying to read '"<< dtag <<
"' keyword in " <<
m_fileOptions << endl);
644 if(!path_to_VAimage.empty())
650 Cerr(
"***** i1TCModel::InitializeSpecific()-> Error reading Interfile : " << path_to_VAimage <<
" !" << endl);
660 Cerr(
"***** i1TCModel::InitializeSpecific() -> Voxels values in provided Va image ;"<< path_to_VAimage <<
" must be between [0,1]. One voxel value is "<<
m_fileOptions << endl);
673 Cerr(
"***** i1TCModel::InitializeSpecific() -> Error while trying to read configuration file at: " <<
m_fileOptions << endl);
687 "1cpt model configuration"))
689 Cerr(
"***** i1TCModel::InitializeSpecific() -> Failed to correctly read the list of options !" << endl);
764 Cerr(
"***** i1TCModel::InitializeSpecific() -> Error, plasma input curve has not been initialized !" << endl);
840 for (
int pl=0 ; pl<
m_nnlsN ; pl++)
841 for (
int pc=0 ; pc<
m_nnlsN ; pc++)
849 mp_RRnum[ th ]->SetMatriceElt( pl , 0 , 0. ) ;
864 Cout(
"i1TCModel::InitializeSpecific() -> Input Plasma TAC coefficients :" << endl);
869 Cout(
"i1TCModel::InitializeSpecific() -> Plasma TAC coefficients :" << endl);
876 Cout(
"i1TCModel::InitializeSpecific() -> Image update from estimated parametric images is disabled" << endl);
900 if(
m_verbose >=2)
Cout(
"i1TCModel::EstimateModelParameters ..." <<endl);
905 Cerr(
"***** i1TCModel::EstimateModelParameters() -> Called while not initialized !" << endl);
915 Cerr(
"***** i1TCModel::EstimateModelParameters() -> A problem occured when estimating 1cpt model parameters with LS method !" << endl);
923 Cerr(
"***** i1TCModel::EstimateModelParameters() -> A problem occured when estimating 1cpt model parameters with Basis functions method !" << endl);
931 Cerr(
"***** i1TCModel::EstimateModelParameters() -> A problem occured when estimating 1cpt model parameters with NNLS method !" << endl);
942 #pragma omp parallel for private(v) schedule(guided) 1004 if(
m_verbose >=2)
Cout(
"i1TCModel::EstimateModelParametersWithNNLS ..." <<endl);
1009 Cerr(
"***** i1TCModel::EstimateModelParametersWithNNLS() -> Called while not initialized !" << endl);
1019 #pragma omp parallel for private(v) schedule(guided) 1028 th = omp_get_thread_num();
1047 int *nnls_index =
new int[
m_nnlsN ];
1052 FLTNB *dptr=_p_nnlsMat;
1077 for(
size_t t=0; t<nb_tf; t++)
1078 ct[t] = ap_ImageS->
m4p_image[t][rg][cg][v];
1083 for(
size_t t=0; t<nb_tf; t++)
1086 if (ct[t]<0) ct[t]=0.;
1097 for(
size_t t=0; t<nb_tf; t++)
1099 if(cti[t]<0) cti[t] = 0.;
1102 _2p_nnlsA[0][t]=-cti[t]*
mp_w[t];
1104 _2p_nnlsA[1][t]=cpi[t]*mp_w[t];
1106 _2p_nnlsA[2][t]=cp[t]*mp_w[t];
1109 _p_nnlsB[t]=ct[t]*mp_w[t];
1112 if(
NNLS(_2p_nnlsA, nb_tf, m_nnlsN, _p_nnlsB, nnls_x, &nnls_rnorm,
1113 nnls_wp, nnls_zz, nnls_index) )
1135 for(
size_t t=0; t<nb_tf; t++)
1139 _2p_nnlsA[0][t] = -cti[t]*
mp_w[t];
1141 _2p_nnlsA[1][t] = cpi[t]*mp_w[t];
1144 _p_nnlsB[t] = ct[t]*mp_w[t];
1148 if(
NNLS(_2p_nnlsA, nb_tf, m_nnlsN, _p_nnlsB, nnls_x, &nnls_rnorm,
1149 nnls_wp, nnls_zz, nnls_index) )
1159 if(nnls_x)
delete[]nnls_x;
1160 if(nnls_wp)
delete[]nnls_wp;
1161 if(nnls_index)
delete[]nnls_index;
1183 if(
m_verbose >=2)
Cout(
"i1TCModel::EstimateModelParametersWithBF ..." <<endl);
1188 Cerr(
"***** i1TCModel::EstimateModelParametersWithBF() -> Called while not initialized !" << endl);
1193 Cerr(
"i1TCModel::EstimateModelParametersWithBF -> Not yet implemented !!!" <<endl);
1218 if(
m_verbose >=2)
Cout(
"i1TCModel::EstimateModelParametersWithLS ..." <<endl);
1223 Cerr(
"***** i1TCModel::EstimateModelParametersWithLS() -> Called while not initialized !" << endl);
1232 bool error_in_th_loop =
false;
1233 string error_msg =
"";
1238 #pragma omp parallel for private(v) schedule(guided) 1247 th = omp_get_thread_num();
1260 for(
size_t t=0; t<nb_tf; t++)
1261 ct[t] = ap_ImageS->
m4p_image[t][rg][cg][v];
1267 for(
size_t t=0; t<nb_tf; t++)
1270 if (ct[t]<0) ct[t]=0.;
1286 LS_Mod[0][t] = -cti[t];
1288 LS_Mod[1][t] = cpi[t];
1291 LS_Mod[2][t] = cp[t];
1305 error_in_th_loop =
true;
1306 error_msg =
"***** i1TCModel::EstimateModelParametersWithLS() -> An error occurred when minimizing WRSS function with ridge-regression!";
1319 error_in_th_loop =
true;
1320 error_msg =
"***** i1TCModel::EstimateModelParametersWithLS() -> An error occurred when minimizing WRSS function with standard LS!";
1335 if( LS_res )
delete[] LS_res;
1338 if ( error_in_th_loop ==
true )
1340 Cerr(error_msg << endl);
1389 Cerr(
"***** i1TCModel::RRLS() -> Called while not initialized !" << endl);
1398 Cerr(
"***** i1TCModel::RRLS() -> Function is called while the optimization method is not Least-square !" << endl);
1404 th = omp_get_thread_num();
1421 for( uint16_t t=0 ; t<a_nT ; t++ )
1425 for( uint16_t p=0 ; p<a_nP ; p++ )
1426 for( uint16_t t=0 ; t<a_nT ; t++ )
1433 cout <<
" ----------------------------------------------- " << endl;
1434 cout <<
" X->Describe() " << endl;
1435 cout <<
" ----------------------------------------------- " << endl;
1443 for( uint16_t p=0 ; p<a_nP ; p++ )
1444 for( uint16_t t=0 ; t<a_nT ; t++ )
1445 X->
SetMatriceElt( t , p , a2p_model[ p ][ t ] * ap_w[ t ] );
1451 cout <<
" ----------------------------------------------- " << endl;
1452 cout <<
" Xt->Describe() " << endl;
1454 cout <<
" ----------------------------------------------- " << endl;
1462 for (
int pl=0 ; pl<a_nP ; pl++)
1463 for (
int pc=0 ; pc<a_nP ; pc++)
1472 cout <<
" ----------------------------------------------- " << endl;
1473 cout <<
" X'WX->Describe() " << endl;
1475 cout <<
" ----------------------------------------------- " << endl;
1481 if ( XtX->
Inverse( LSden ) == 1 )
1483 for( uint16_t p=0 ; p<a_nP ; p++ )
1484 ap_result[ p ] = 0.;
1492 cout <<
" ----------------------------------------------- " << endl;
1493 cout <<
" RRLS denominator ->Describe() " << endl;
1495 cout <<
" ----------------------------------------------- " << endl;
1505 cout <<
" ----------------------------------------------- " << endl;
1506 cout <<
" LS numerator ->Describe() " << endl;
1508 cout <<
" ----------------------------------------------- " << endl;
1516 for (
int p=0 ; p<a_nP ; p++)
1523 cout <<
" ----------------------------------------------- " << endl;
1524 cout <<
" RRLS numerator ->Describe() " << endl;
1526 cout <<
" ----------------------------------------------- " << endl;
1539 cout <<
" ----------------------------------------------- " << endl;
1540 cout <<
" Theta->Describe() " << endl;
1546 for( uint16_t p=0 ; p<a_nP ; p++ )
1591 Cerr(
"***** i1TCModel::LS() -> Called while not initialized !" << endl);
1600 Cerr(
"***** i1TCModel::LS() -> Function is called while the optimization method is not Least-square !" << endl);
1606 th = omp_get_thread_num();
1619 for( uint16_t t=0 ; t<a_nT ; t++ )
1623 for( uint16_t p=0 ; p<a_nP ; p++ )
1624 for( uint16_t t=0 ; t<a_nT ; t++ )
1631 cout <<
" ----------------------------------------------- " << endl;
1632 cout <<
" X->Describe() " << endl;
1633 cout <<
" ----------------------------------------------- " << endl;
1641 for( uint16_t p=0 ; p<a_nP ; p++ )
1642 for( uint16_t t=0 ; t<a_nT ; t++ )
1643 X->
SetMatriceElt( t , p , a2p_model[ p ][ t ] * ap_w[ t ] );
1649 cout <<
" ----------------------------------------------- " << endl;
1650 cout <<
" Xt->Describe() " << endl;
1652 cout <<
" ----------------------------------------------- " << endl;
1662 cout <<
" ----------------------------------------------- " << endl;
1663 cout <<
" X'WX->Describe() " << endl;
1665 cout <<
" ----------------------------------------------- " << endl;
1675 cout <<
" ----------------------------------------------- " << endl;
1676 cout <<
" LS numerator ->Describe() " << endl;
1678 cout <<
" ----------------------------------------------- " << endl;
1684 if ( XtX->
Inverse( LSden ) == 1 )
1686 for( uint16_t p=0 ; p<a_nP ; p++ )
1687 ap_result[ p ] = 0.;
1695 cout <<
" ----------------------------------------------- " << endl;
1696 cout <<
" LS denominator ->Describe() " << endl;
1698 cout <<
" ----------------------------------------------- " << endl;
1708 cout <<
" ----------------------------------------------- " << endl;
1709 cout <<
" Theta->Describe() " << endl;
1715 for( uint16_t p=0 ; p<a_nP ; p++ )
1739 if(
m_verbose >= 2)
Cout(
"i1TCModel::EstimateImageWithModel ... " <<endl);
1744 Cerr(
"***** i1TCModel::EstimateImageWithModel() -> Called while not initialized !" << endl);
1755 DTavg[ 0 ] =
mp_DT2[ 0 ];
1763 #pragma omp parallel for private(v) schedule(guided) 1772 if( ( K1>0. || k2>0. || Va>0.)
1773 && ( K1>=0. && k2>=0. && Va>=0. ) )
1786 b = cti + DTavg[t]*b_ct ;
1789 / ( 1 + k2*DTavg[ t ] );
1801 / (1 + k2*DTavg[t]);
1804 cti +=
WPOinc( t, ct, b_ct, bb_ct, n_ct);
1807 cti +=
mp_DT2[t]*(ct + b_ct);
1811 ap_ImageS->
m4p_image[t][rg][cg][v] = (new_value > 0) ?
1849 FLTNB wpo_p=0., wpo_p_n=0, wpo_fd_n=0., wpo_bd=0.;
1850 uint32_t t = a_time;
1853 Itac =
mp_DT2[t]*tac*0.5;
1857 wpo_fd_n = 2*
mp_DT2[1] * ( 0.5*( tac + b_tac ) + wpo_p_n *
mp_wpoQ[t+1] );
1858 wpo_bd =
mp_DT2[t-1] * ( b_tac + tac );
1867 wpo_fd_n = 2*
mp_DT2[t] * ( 0.5*( tac + b_tac ) + wpo_p_n *
mp_wpoQ[t+1] );
1868 wpo_bd = 2*
mp_DT2[t-1] * ( 0.5*( b_tac + tac ) + wpo_p /
mp_wpoQ[t] );
1875 wpo_bd = 2*
mp_DT2[t-1] * ( 0.5*( b_tac + tac ) + wpo_p /
mp_wpoQ[t] );
1901 return WPO(ap_tac, ap_citac, a_th);
1903 return Trapz(ap_tac, ap_citac);
1933 for(uint32_t t=0 ; t<nb_tf ; t++)
1938 WPO_BD[t] = ap_tac[0] ;
1939 WPO_FD[t] = ap_tac[0] ;
1944 WPO_FD[t] = DT2[t-1] *( ( ap_tac[ t-1 ] ) );
1945 WPO_BD[t] = DT2[t-1] *( ( ap_tac[ t-1 ] + ap_tac[ t ] ));
1949 WPO_P[t] = ( ap_tac[ t-1 ] - ( DT2[t]*ap_tac[ t-2 ] + DT2[t-1]*ap_tac[ t ] ) / ( DT2[t] + DT2[t-1] ) ) / 6;
1950 WPO_FD[t] = 2*DT2[t-1] * ( 0.5 *( ap_tac[ t-1 ] + ap_tac[ t-2 ] ) + WPO_P[t] *
mp_wpoQ[t] );
1951 WPO_BD[t] = 2*DT2[t-1] * ( 0.5 *( ap_tac[ t-1 ] + ap_tac[ t ] ) + WPO_P[t] /
mp_wpoQ[t] );
1956 for(uint32_t t=0 ; t<nb_tf ; t++)
1960 itac = ap_tac[ t ] * DT2[ t ] ;
1961 else if( t == (nb_tf-1) )
1964 itac = (1-
mp_wpoA[ t ]) * WPO_BD[ t ] +
mp_wpoA[ t ]*WPO_FD[ t+1 ];
1968 ap_citac[ t ] = itac ;
1970 ap_citac[ t ] = ap_citac[ t-1] + itac;
1997 for(uint32_t t=0 ; t<nb_tf ; t++)
2004 itac = ap_tac[0] * DT2[0] *0.5;
2009 itac = ( ap_tac[t]+ap_tac[t-1] ) * DT2[t];
2010 ap_citac[t] += ap_citac[t-1] + itac;
#define INTF_LERP_DISABLED
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...
int Inverse(oMatrix *ap_MtxResult)
Inverse the matrix on which the function is called An error is returned if the matrix is not square...
int CheckSpecificParameters()
This function is used to check whether all member variables have been correctly initialized or not...
FLTNB * mp_parLowerBounds
i1TCModel()
Constructor of i1TCModel. Simply set all data members to default values.
int InitializeSpecific()
This function is used to initialize parametric images and basis functions.
bool m_ridgeRegressionFlag
oImageDimensionsAndQuantification * mp_ID
This is the mother class of dynamic model classes.
~i1TCModel()
Destructor of i1TCModel.
bool m_saveBlacklistedImageMaskFlag
int Trapz(FLTNB *ap_tac, FLTNB *ap_citac)
return integral for each time point t (WPO_S), and cumulative integral (cumWPO_S ) ...
int EstimateModelParametersWithBF(oImageSpace *ap_ImageS)
Estimate K1, k2, Va parametric images using basis functions approach.
int RRLS(uint16_t a_nP, uint16_t a_nT, FLTNB **a2p_model, FLTNB *ap_data, FLTNB *ap_w, FLTNB *ap_result)
Non-linear least square estimator with Ridge-Regression.
int WPO(FLTNB *ap_tac, FLTNB *ap_citac, int a_th)
return integral for each time point t (WPO_S), and cumulative integral (cumWPO_S ) ...
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...
FLTNB * mp_parUpperBounds
int ReadAndCheckOptionsList(string a_listOptions)
This function is used to read parameters from a string.
int EstimateModelParameters(oImageSpace *ap_Image, int a_ite, int a_sset)
Estimate K1, k2, Va parametric images.
int ReadAndCheckConfigurationFileSpecific()
This function is used to read options from a configuration file.
int EstimateImageWithModel(oImageSpace *ap_Image, int a_ite, int a_sset)
Estimate image using model parametric images and basis functions.
HPFLTNB GetMatriceElt(uint16_t l, uint16_t c)
void Describe()
Display the element of the matrix.
FLTNB ** m2p_nestedModelTimeBasisFunctions
#define KEYWORD_MANDATORY
void ShowHelp()
Print out specific help about the implementation of the model model and its initialization.
int EstimateModelParametersWithLS(oImageSpace *ap_ImageS)
Estimate K1, k2, Va parametric images using LS.
Declaration of class i1TCModel.
int GetNbCardGates()
Get the number of cardiac gates.
FLTNB ** m2p_parametricImages
This class holds all the matrices in the image domain that can be used in the algorithm: image...
int m_OptimisationMethodFlag
int GetNbThreadsForImageComputation()
Get the number of threads used for image operations.
int GetNbTimeFrames()
Get the number of time frames.
INTNB GetNbVoxXYZ()
Get the total number of voxels.
Structure designed for basic matrices operations.
int EstimateModelParametersWithNNLS(oImageSpace *ap_ImageS)
Estimate K1, k2, Va parametric images using NNLS.
int GetNbRespGates()
Get the number of respiratory gates.
int LS(uint16_t a_nP, uint16_t a_nT, FLTNB **a2p_model, FLTNB *ap_data, FLTNB *ap_w, FLTNB *ap_result)
Non-linear least square estimator.
int Multiplication(oMatrix *ap_Mtx, oMatrix *ap_MtxResult)
Multiply the member matrix with the matrix provided in 1st parameter Return the result in the matric ...
FLTNB WPOinc(uint32_t a_time, FLTNB tac, FLTNB b_tac, FLTNB bb_tac, FLTNB n_tac)
Estimate the next integration value for a specific time point of a tac using WPO. ...
FLTNB GetFrameDurationInSec(int a_bed, int a_frame)
Get the frame duration for the given bed, in seconds as a FLTNB.
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.
int SetMatriceElt(uint16_t l, uint16_t c, HPFLTNB a_val)
Set the matrix element corresponding to the argument indices with the provided value.
int ReadStringOption(const string &a_input, T *ap_return, int a_nbElts, const string &sep, const string &a_option)
Parse the 'a_input' string corresponding to the 'a_option' into 'a_nbElts' elements, using the 'sep' separator. The results are returned in the templated 'ap_return' dynamic templated array. Call "ConvertFromString()" to perform the correct conversion depending on the type of the data to convert.
int IntfReadImage(const string &a_pathToHeaderFile, FLTNB *ap_ImgMatrix, oImageDimensionsAndQuantification *ap_ID, int vb, bool a_lerpFlag)
Main function dedicated to Interfile 3D image loading.
FLTNB * mp_blackListedvoxelsImage
int IntegrateTAC(FLTNB *ap_tac, FLTNB *ap_citac, int a_th)
Call one of the TAC integration method.
int Transpose(oMatrix *a_MtxResult)
Transpose the elements of the matrix.