![]() |
CASToR
1.0
Tomographic Reconstruction (PET/SPECT)
|
00001 00002 /* 00003 Implementation of class iProjectorIncrementalSiddonMulti 00004 00005 - separators: done 00006 - doxygen: done 00007 - default initialization: done 00008 - CASTOR_DEBUG: 00009 - CASTOR_VERBOSE: 00010 */ 00011 00018 #include "iProjectorIncrementalSiddonMulti.hh" 00019 #include "sOutputManager.hh" 00020 00021 // ===================================================================== 00022 // --------------------------------------------------------------------- 00023 // --------------------------------------------------------------------- 00024 // ===================================================================== 00025 00026 iProjectorIncrementalSiddonMulti::iProjectorIncrementalSiddonMulti() : vProjector() 00027 { 00028 m_nbLines = -1; 00029 // This projector is not compatible with SPECT attenuation correction because 00030 // the voxels contributing to the line are not strictly ordered with respect to 00031 // their distance to point 2 (due to the use of multiple lines that are 00032 // stack one after the other) 00033 m_compatibleWithSPECTAttenuationCorrection = false; 00034 } 00035 00036 // ===================================================================== 00037 // --------------------------------------------------------------------- 00038 // --------------------------------------------------------------------- 00039 // ===================================================================== 00040 00041 iProjectorIncrementalSiddonMulti::~iProjectorIncrementalSiddonMulti() 00042 { 00043 } 00044 00045 // ===================================================================== 00046 // --------------------------------------------------------------------- 00047 // --------------------------------------------------------------------- 00048 // ===================================================================== 00049 00050 int iProjectorIncrementalSiddonMulti::ReadConfigurationFile(const string& a_configurationFile) 00051 { 00052 // Read the transaxial FWHM option 00053 string key_word = "number of lines"; 00054 if (ReadDataASCIIFile(a_configurationFile, key_word, &m_nbLines, 1, KEYWORD_MANDATORY)) 00055 { 00056 Cerr("***** iProjectorIncrementalSiddonMulti::ReadConfigurationFile() -> Failed to get the '" << key_word << "' keyword !" << endl); 00057 return 1; 00058 } 00059 // Normal end 00060 return 0; 00061 } 00062 00063 // ===================================================================== 00064 // --------------------------------------------------------------------- 00065 // --------------------------------------------------------------------- 00066 // ===================================================================== 00067 00068 int iProjectorIncrementalSiddonMulti::ReadOptionsList(const string& a_optionsList) 00069 { 00070 // Read them 00071 if (ReadStringOption(a_optionsList, &m_nbLines, 1, ",", "Multi-Siddon configuration")) 00072 { 00073 Cerr("***** iProjectorIncrementalSiddonMulti::ReadConfigurationFile() -> Failed to correctly read the list of options !" << endl); 00074 return 1; 00075 } 00076 // Normal end 00077 return 0; 00078 } 00079 00080 // ===================================================================== 00081 // --------------------------------------------------------------------- 00082 // --------------------------------------------------------------------- 00083 // ===================================================================== 00084 00085 void iProjectorIncrementalSiddonMulti::ShowHelpSpecific() 00086 { 00087 cout << "This projector uses multiple ray-tracing for a single event in order to estimate the solid angle contribution." << endl; 00088 cout << "For each line of an event, the end points of the line are randomly chosen inside the detector element." << endl; 00089 cout << "The ray-tracing is performed with the incremental Siddon algorithm (see incrementalSiddon projector)." << endl; 00090 cout << "The only parameter of this projector is the number of lines to use per event:" << endl; 00091 cout << " number of lines: the number of lines used per event" << endl; 00092 } 00093 00094 // ===================================================================== 00095 // --------------------------------------------------------------------- 00096 // --------------------------------------------------------------------- 00097 // ===================================================================== 00098 00099 int iProjectorIncrementalSiddonMulti::CheckSpecificParameters() 00100 { 00101 // Send an error if less than 1 line 00102 if (m_nbLines<1) 00103 { 00104 Cerr("***** iProjectorIncrementalSiddonMulti::CheckSpecificParameters() -> The provided number of lines is less than 1 !" << endl); 00105 return 1; 00106 } 00107 // Normal end 00108 return 0; 00109 } 00110 00111 // ===================================================================== 00112 // --------------------------------------------------------------------- 00113 // --------------------------------------------------------------------- 00114 // ===================================================================== 00115 00116 int iProjectorIncrementalSiddonMulti::InitializeSpecific() 00117 { 00118 // Verbose 00119 if (m_verbose>=2) Cout("iProjectorIncrementalSiddonMulti::Initialize() -> Use incremental Siddon projector with " << m_nbLines << " lines per event" << endl); 00120 // Normal end 00121 return 0; 00122 } 00123 00124 // ===================================================================== 00125 // --------------------------------------------------------------------- 00126 // --------------------------------------------------------------------- 00127 // ===================================================================== 00128 00129 INTNB iProjectorIncrementalSiddonMulti::EstimateMaxNumberOfVoxelsPerLine() 00130 { 00131 // Find the maximum number of voxels along a given dimension 00132 INTNB max_nb_voxels_in_dimension = mp_ImageDimensionsAndQuantification->GetNbVoxX(); 00133 if (mp_ImageDimensionsAndQuantification->GetNbVoxY()>max_nb_voxels_in_dimension) max_nb_voxels_in_dimension = mp_ImageDimensionsAndQuantification->GetNbVoxY(); 00134 if (mp_ImageDimensionsAndQuantification->GetNbVoxZ()>max_nb_voxels_in_dimension) max_nb_voxels_in_dimension = mp_ImageDimensionsAndQuantification->GetNbVoxZ(); 00135 // We should have at most 4 voxels in a given plane, so multiply by 4 00136 // (note: this is not true however it ensures no overflow and is already quite optimized for RAM usage !) 00137 max_nb_voxels_in_dimension *= 4; 00138 // Finally multiply by the number of lines 00139 max_nb_voxels_in_dimension *= ((INTNB)m_nbLines); 00140 // Return the value 00141 return max_nb_voxels_in_dimension; 00142 } 00143 00144 // ===================================================================== 00145 // --------------------------------------------------------------------- 00146 // --------------------------------------------------------------------- 00147 // ===================================================================== 00148 00149 int iProjectorIncrementalSiddonMulti::ProjectWithoutTOF(int a_direction, oProjectionLine* ap_ProjectionLine ) 00150 { 00151 #ifdef CASTOR_DEBUG 00152 if (!m_initialized) 00153 { 00154 Cerr("***** iProjectorIncrementalSiddonMulti::ProjectWithoutTOF() -> Called while not initialized !" << endl); 00155 Exit(EXIT_DEBUG); 00156 } 00157 #endif 00158 00159 #ifdef CASTOR_VERBOSE 00160 if (m_verbose>=10) 00161 { 00162 string direction = ""; 00163 if (a_direction==FORWARD) direction = "forward"; 00164 else direction = "backward"; 00165 Cout("iProjectorIncrementalSiddonMulti::Project without TOF -> Project line '" << ap_ProjectionLine << "' in " << direction << " direction" << endl); 00166 } 00167 #endif 00168 00169 // Loop on the number of lines 00170 for (int line=0; line<m_nbLines; line++) 00171 { 00172 // Get random positions from the scanner (mean depth of interaction intrinsicaly taken into account), taking POI into account if any 00173 if (mp_Scanner->GetRdmPositionsAndOrientations( ap_ProjectionLine->GetIndex1(), ap_ProjectionLine->GetIndex2(), 00174 ap_ProjectionLine->GetPosition1(), ap_ProjectionLine->GetPosition2(), 00175 ap_ProjectionLine->GetOrientation1(), ap_ProjectionLine->GetOrientation2() )) 00176 { 00177 Cerr("***** vProjector::Project() -> A problem occured while getting positions and orientations from scanner !" << endl); 00178 return 1; 00179 } 00180 00181 // Get end points position 00182 FLTNB* event1Float = ap_ProjectionLine->GetPosition1(); 00183 FLTNB* event2Float = ap_ProjectionLine->GetPosition2(); 00184 double event1[3] = { event1Float[0], event1Float[1], event1Float[2] }; 00185 double event2[3] = { event2Float[0], event2Float[1], event2Float[2] }; 00186 00187 // ************************************** 00188 // STEP 1: LOR length calculation 00189 // ************************************** 00190 double length_LOR = ap_ProjectionLine->GetLength(); 00191 00192 // ************************************** 00193 // STEP 2: Compute entrance voxel indexes 00194 // ************************************** 00195 double alphaFirst[] = { 0.0, 0.0, 0.0 }; 00196 double alphaLast[] = { 0.0, 0.0, 0.0 }; 00197 00198 double alphaMin = 0.0f, alphaMax = 1.0f; 00199 double delta_pos[] = { 00200 event2[ 0 ] - event1[ 0 ], 00201 event2[ 1 ] - event1[ 1 ], 00202 event2[ 2 ] - event1[ 2 ] 00203 }; 00204 00205 // Computation of alphaMin et alphaMax values (entrance and exit point of the ray) 00206 for( int i = 0; i < 3; ++i ) 00207 { 00208 if( delta_pos[ i ] != 0.0 ) 00209 { 00210 alphaFirst[i] = (-mp_halfFOV[i] - event1[i]) / (delta_pos[ i ]); 00211 alphaLast[i] = ( mp_halfFOV[i] - event1[i]) / (delta_pos[ i ]); 00212 alphaMin = (std::max)(alphaMin,(std::min)(alphaFirst[i],alphaLast[i])); 00213 alphaMax = (std::min)(alphaMax,(std::max)(alphaFirst[i],alphaLast[i])); 00214 } 00215 } 00216 00217 // if alphaMax is less than or equal to alphaMin no intersection 00218 // and return an empty buffer 00219 if( alphaMax <= alphaMin ) return 0; 00220 00221 // Now we have to find the indices of the particular plane 00222 // (iMin,iMax), (jMin,jMax), (kMin,kMax) 00223 int iMin = 0, iMax = 0; 00224 int jMin = 0, jMax = 0; 00225 int kMin = 0, kMax = 0; 00226 00227 // For the X-axis 00228 if( delta_pos[ 0 ] > 0.0f) 00229 { 00230 iMin = ::ceil( ( mp_nbVox[ 0 ] + 1 ) - ( mp_halfFOV[ 0 ] - alphaMin * delta_pos[ 0 ] - event1[ 0 ] ) / mp_sizeVox[0] ); 00231 iMax = ::floor( 1 + ( event1[ 0 ] + alphaMax * delta_pos[ 0 ] - (-mp_halfFOV[ 0 ]) ) / mp_sizeVox[0] ); 00232 } 00233 else if( delta_pos[ 0 ] < 0.0 ) 00234 { 00235 iMin = ::ceil( ( mp_nbVox[ 0 ] + 1 ) - ( mp_halfFOV[ 0 ] - alphaMax * delta_pos[ 0 ] - event1[ 0 ] ) / mp_sizeVox[0] ); 00236 iMax = ::floor( 1 + ( event1[ 0 ] + alphaMin * delta_pos[ 0 ] - (-mp_halfFOV[ 0 ]) ) / mp_sizeVox[0] ); 00237 } 00238 if( delta_pos[ 0 ] == 0 ) 00239 { 00240 iMin = 1, iMax = 0; 00241 } 00242 00243 // For the Y-axis 00244 if( delta_pos[ 1 ] > 0 ) 00245 { 00246 jMin = ::ceil( ( mp_nbVox[ 1 ] + 1 ) - ( mp_halfFOV[ 1 ] - alphaMin * delta_pos[ 1 ] - event1[ 1 ] ) / mp_sizeVox[1] ); 00247 jMax = ::floor( 1 + ( event1[ 1 ] + alphaMax * delta_pos[ 1 ] - (-mp_halfFOV[ 1 ]) ) / mp_sizeVox[1] ); 00248 } 00249 else if( delta_pos[ 1 ] < 0 ) 00250 { 00251 jMin = ::ceil( ( mp_nbVox[ 1 ] + 1 ) - ( mp_halfFOV[ 1 ] - alphaMax * delta_pos[ 1 ] - event1[ 1 ] ) / mp_sizeVox[1] ); 00252 jMax = ::floor( 1 + ( event1[ 1 ] + alphaMin * delta_pos[ 1 ] - (-mp_halfFOV[ 1 ]) ) / mp_sizeVox[1] ); 00253 } 00254 else if( delta_pos[ 1 ] == 0 ) 00255 { 00256 jMin = 1, jMax = 0; 00257 } 00258 00259 // For the Z-axis 00260 if( delta_pos[ 2 ] > 0 ) 00261 { 00262 kMin = ::ceil( ( mp_nbVox[ 2 ] + 1 ) - ( mp_halfFOV[ 2 ] - alphaMin * delta_pos[ 2 ] - event1[ 2 ] ) / mp_sizeVox[2] ); 00263 kMax = ::floor( 1 + ( event1[ 2 ] + alphaMax * delta_pos[ 2 ] - (-mp_halfFOV[ 2 ]) ) / mp_sizeVox[2] ); 00264 } 00265 else if( delta_pos[ 2 ] < 0 ) 00266 { 00267 kMin = ::ceil( ( mp_nbVox[ 2 ] + 1 ) - ( mp_halfFOV[ 2 ] - alphaMax * delta_pos[ 2 ] - event1[ 2 ] ) / mp_sizeVox[2] ); 00268 kMax = ::floor( 1 + ( event1[ 2 ] + alphaMin * delta_pos[ 2 ] - (-mp_halfFOV[ 2 ]) ) / mp_sizeVox[2] ); 00269 } 00270 else if( delta_pos[ 2 ] == 0 ) 00271 { 00272 kMin = 1, kMax = 0; 00273 } 00274 00275 // Computing the last term n number of intersection 00276 INTNB n = ( iMax - iMin + 1 ) + ( jMax - jMin + 1 ) 00277 + ( kMax - kMin + 1 ); 00278 00279 // Array storing the first alpha in X, Y and Z 00280 double alpha_XYZ[ 3 ] = { 1.0, 1.0, 1.0 }; 00281 00282 // Computing the first alpha in X 00283 if( delta_pos[ 0 ] > 0 ) 00284 alpha_XYZ[ 0 ] = ( ( (-mp_halfFOV[ 0 ]) + ( iMin - 1 ) * mp_sizeVox[0] ) - event1[ 0 ] ) / delta_pos[ 0 ]; 00285 else if( delta_pos[ 0 ] < 0 ) 00286 alpha_XYZ[ 0 ] = ( ( (-mp_halfFOV[ 0 ]) + ( iMax - 1 ) * mp_sizeVox[0] ) - event1[ 0 ] ) / delta_pos[ 0 ]; 00287 00288 // Computing the first alpha in Y 00289 if( delta_pos[ 1 ] > 0 ) 00290 alpha_XYZ[ 1 ] = ( ( (-mp_halfFOV[ 1 ]) + ( jMin - 1 ) * mp_sizeVox[1] ) - event1[ 1 ] ) / delta_pos[ 1 ]; 00291 else if( delta_pos[ 1 ] < 0 ) 00292 alpha_XYZ[ 1 ] = ( ( (-mp_halfFOV[ 1 ]) + ( jMax - 1 ) * mp_sizeVox[1] ) - event1[ 1 ] ) / delta_pos[ 1 ]; 00293 00294 // Computing the first alpha in Z 00295 if( delta_pos[ 2 ] > 0 ) 00296 alpha_XYZ[ 2 ] = ( ( (-mp_halfFOV[ 2 ]) + ( kMin - 1 ) * mp_sizeVox[2] ) - event1[ 2 ] ) / delta_pos[ 2 ]; 00297 else if( delta_pos[ 2 ] < 0 ) 00298 alpha_XYZ[ 2 ] = ( ( (-mp_halfFOV[ 2 ]) + ( kMax - 1 ) * mp_sizeVox[2] ) - event1[ 2 ] ) / delta_pos[ 2 ]; 00299 00300 // Computing the alpha updating 00301 double alpha_u[ 3 ] = { 00302 mp_sizeVox[0] / std::fabs( delta_pos[ 0 ] ), 00303 mp_sizeVox[1] / std::fabs( delta_pos[ 1 ] ), 00304 mp_sizeVox[2] / std::fabs( delta_pos[ 2 ] ) 00305 }; 00306 00307 // Computing the index updating 00308 INTNB index_u[ 3 ] = { 00309 (delta_pos[ 0 ] < 0) ? -1 : 1, 00310 (delta_pos[ 1 ] < 0) ? -1 : 1, 00311 (delta_pos[ 2 ] < 0) ? -1 : 1 00312 }; 00313 00314 // Check which alpha is the min/max and increment 00315 if( alpha_XYZ[ 0 ] == alphaMin ) alpha_XYZ[ 0 ] += alpha_u[ 0 ]; 00316 if( alpha_XYZ[ 1 ] == alphaMin ) alpha_XYZ[ 1 ] += alpha_u[ 1 ]; 00317 if( alpha_XYZ[ 2 ] == alphaMin ) alpha_XYZ[ 2 ] += alpha_u[ 2 ]; 00318 00319 // Computing the minimum value in the alpha_XYZ buffer 00320 double const min_alpha_XYZ = (std::min)( alpha_XYZ[ 0 ], 00321 (std::min)( alpha_XYZ[ 1 ], alpha_XYZ[ 2 ] ) ); 00322 00323 // Computing the first index of intersection 00324 double const alphaMid = ( min_alpha_XYZ + alphaMin ) * 0.5; 00325 INTNB index_ijk[ 3 ] = { 00326 1 + (int)( ( event1[ 0 ] + alphaMid * delta_pos[ 0 ] - (-mp_halfFOV[ 0 ]) ) / mp_sizeVox[0] ), 00327 1 + (int)( ( event1[ 1 ] + alphaMid * delta_pos[ 1 ] - (-mp_halfFOV[ 1 ]) ) / mp_sizeVox[1] ), 00328 1 + (int)( ( event1[ 2 ] + alphaMid * delta_pos[ 2 ] - (-mp_halfFOV[ 2 ]) ) / mp_sizeVox[2] ) 00329 }; 00330 00331 INTNB const w = mp_nbVox[ 0 ]; 00332 INTNB const wh = w * mp_nbVox[ 1 ]; 00333 00334 // Loop over the number of plane to cross 00335 double alpha_c = alphaMin; 00336 FLTNB coeff = 0.0; 00337 INTNB numVox = 0; 00338 for( int nP = 0; nP < n - 1; ++nP ) 00339 { 00340 if( ( alpha_XYZ[ 0 ] <= alpha_XYZ[ 1 ] ) 00341 && ( alpha_XYZ[ 0 ] <= alpha_XYZ[ 2 ] ) ) 00342 { 00343 // Storing values 00344 if( ( alpha_XYZ[ 0 ] >= alphaMin ) 00345 && ( index_ijk[ 0 ] - 1 >= 0 ) 00346 && ( index_ijk[ 0 ] - 1 <= mp_nbVox[ 0 ] - 1 ) 00347 && ( index_ijk[ 1 ] - 1 >= 0 ) 00348 && ( index_ijk[ 1 ] - 1 <= mp_nbVox[ 1 ] - 1 ) 00349 && ( index_ijk[ 2 ] - 1 >= 0 ) 00350 && ( index_ijk[ 2 ] - 1 <= mp_nbVox[ 2 ] - 1 ) ) 00351 { 00352 coeff = ( alpha_XYZ[ 0 ] - alpha_c ) * length_LOR; 00353 numVox = ( index_ijk[ 0 ] - 1 ) + ( ( index_ijk[ 1 ] - 1 ) ) * w + ( ( index_ijk[ 2 ] - 1 ) ) * wh; 00354 ap_ProjectionLine->AddVoxel(a_direction, numVox, coeff/((FLTNB)m_nbLines)); 00355 } 00356 00357 // Increment values 00358 alpha_c = alpha_XYZ[ 0 ]; 00359 alpha_XYZ[ 0 ] += alpha_u[ 0 ]; 00360 index_ijk[ 0 ] += index_u[ 0 ]; 00361 } 00362 else if( ( alpha_XYZ[ 1 ] < alpha_XYZ[ 0 ] ) 00363 && ( alpha_XYZ[ 1 ] <= alpha_XYZ[ 2 ] ) ) 00364 { 00365 // Storing values 00366 if( ( alpha_XYZ[ 1 ] >= alphaMin ) 00367 && ( index_ijk[ 0 ] - 1 >= 0 ) 00368 && ( index_ijk[ 0 ] - 1 <= mp_nbVox[ 0 ] - 1 ) 00369 && ( index_ijk[ 1 ] - 1 >= 0 ) 00370 && ( index_ijk[ 1 ] - 1 <= mp_nbVox[ 1 ] - 1 ) 00371 && ( index_ijk[ 2 ] - 1 >= 0 ) 00372 && ( index_ijk[ 2 ] - 1 <= mp_nbVox[ 2 ] - 1 ) ) 00373 { 00374 coeff = ( alpha_XYZ[ 1 ] - alpha_c ) * length_LOR; 00375 numVox = ( index_ijk[ 0 ] - 1 ) + ( ( index_ijk[ 1 ] - 1 ) ) * w + ( ( index_ijk[ 2 ] - 1 ) ) * wh; 00376 ap_ProjectionLine->AddVoxel(a_direction, numVox, coeff/((FLTNB)m_nbLines)); 00377 } 00378 00379 // Increment values 00380 alpha_c = alpha_XYZ[ 1 ]; 00381 alpha_XYZ[ 1 ] += alpha_u[ 1 ]; 00382 index_ijk[ 1 ] += index_u[ 1 ]; 00383 } 00384 else if( ( alpha_XYZ[ 2 ] < alpha_XYZ[ 0 ] ) 00385 && ( alpha_XYZ[ 2 ] < alpha_XYZ[ 1 ] ) ) 00386 { 00387 // Storing values 00388 if( ( alpha_XYZ[ 2 ] >= alphaMin ) 00389 && ( index_ijk[ 0 ] - 1 >= 0 ) 00390 && ( index_ijk[ 0 ] - 1 <= mp_nbVox[ 0 ] - 1 ) 00391 && ( index_ijk[ 1 ] - 1 >= 0 ) 00392 && ( index_ijk[ 1 ] - 1 <= mp_nbVox[ 1 ] - 1 ) 00393 && ( index_ijk[ 2 ] - 1 >= 0 ) 00394 && ( index_ijk[ 2 ] - 1 <= mp_nbVox[ 2 ] - 1 ) ) 00395 { 00396 coeff = ( alpha_XYZ[ 2 ] - alpha_c ) * length_LOR; 00397 numVox = ( index_ijk[ 0 ] - 1 ) + ( ( index_ijk[ 1 ] - 1 ) ) * w + ( ( index_ijk[ 2 ] - 1 ) ) * wh; 00398 ap_ProjectionLine->AddVoxel(a_direction, numVox, coeff/((FLTNB)m_nbLines)); 00399 } 00400 00401 // Increment values 00402 alpha_c = alpha_XYZ[ 2 ]; 00403 alpha_XYZ[ 2 ] += alpha_u[ 2 ]; 00404 index_ijk[ 2 ] += index_u[ 2 ]; 00405 } 00406 } 00407 } 00408 00409 return 0; 00410 } 00411 00412 // ===================================================================== 00413 // --------------------------------------------------------------------- 00414 // --------------------------------------------------------------------- 00415 // ===================================================================== 00416 00417 int iProjectorIncrementalSiddonMulti::ProjectWithTOFPos(int a_Projector, oProjectionLine* ap_ProjectionLine) 00418 { 00419 Cerr("***** iProjectorIncrementalSiddonMulti::ProjectWithTOFPos() -> Not yet implemented !" << endl); 00420 return 1; 00421 } 00422 00423 // ===================================================================== 00424 // --------------------------------------------------------------------- 00425 // --------------------------------------------------------------------- 00426 // ===================================================================== 00427 00428 int iProjectorIncrementalSiddonMulti::ProjectWithTOFBin(int a_Projector, oProjectionLine* ap_ProjectionLine) 00429 { 00430 Cerr("***** iProjectorIncrementalSiddonMulti::ProjectWithTOFBin() -> Not yet implemented !" << endl); 00431 return 1; 00432 } 00433 00434 // ===================================================================== 00435 // --------------------------------------------------------------------- 00436 // --------------------------------------------------------------------- 00437 // =====================================================================