![]() |
CASToR
1.0
Tomographic Reconstruction (PET/SPECT)
|
00001 00002 /* 00003 Implementation of class iProjectorIncrementalSiddon 00004 00005 - separators: done 00006 - doxygen: done 00007 - default initialization: done 00008 - CASTOR_DEBUG: 00009 - CASTOR_VERBOSE: 00010 */ 00011 00018 #include "iProjectorIncrementalSiddon.hh" 00019 #include "sOutputManager.hh" 00020 00021 // ===================================================================== 00022 // --------------------------------------------------------------------- 00023 // --------------------------------------------------------------------- 00024 // ===================================================================== 00025 00026 iProjectorIncrementalSiddon::iProjectorIncrementalSiddon() : vProjector() 00027 { 00028 // This projector is compatible with SPECT attenuation correction because 00029 // all voxels contributing to a given line are ordered from point 1 (focal) 00030 // to point 2 (scanner). Note that there can be some aliasing problem with 00031 // this incremental method, but which will depend on the voxel size and number. 00032 // However, these problems are quite negligible and pretty uncommon. So we still 00033 // say that this projector is compatible. 00034 m_compatibleWithSPECTAttenuationCorrection = true; 00035 } 00036 00037 // ===================================================================== 00038 // --------------------------------------------------------------------- 00039 // --------------------------------------------------------------------- 00040 // ===================================================================== 00041 00042 iProjectorIncrementalSiddon::~iProjectorIncrementalSiddon() 00043 { 00044 } 00045 00046 // ===================================================================== 00047 // --------------------------------------------------------------------- 00048 // --------------------------------------------------------------------- 00049 // ===================================================================== 00050 00051 int iProjectorIncrementalSiddon::ReadConfigurationFile(const string& a_configurationFile) 00052 { 00053 // No options for incremental siddon 00054 ; 00055 // Normal end 00056 return 0; 00057 } 00058 00059 // ===================================================================== 00060 // --------------------------------------------------------------------- 00061 // --------------------------------------------------------------------- 00062 // ===================================================================== 00063 00064 int iProjectorIncrementalSiddon::ReadOptionsList(const string& a_optionsList) 00065 { 00066 // No options for incremental siddon 00067 ; 00068 // Normal end 00069 return 0; 00070 } 00071 00072 // ===================================================================== 00073 // --------------------------------------------------------------------- 00074 // --------------------------------------------------------------------- 00075 // ===================================================================== 00076 00077 void iProjectorIncrementalSiddon::ShowHelpSpecific() 00078 { 00079 cout << "This projector is a simple line projector that computes the exact path length of a line through the voxels." << endl; 00080 cout << "It is basically an incremental version of the original algorithm proposed by R. L. Siddon (see the classicSiddon projector)." << endl; 00081 cout << "It is implemented from the following published paper:" << endl; 00082 cout << "F. Jacobs et al, \"A fast algorithm to calculate the exact radiological path through a pixel or voxel space\", J. Comput. Inf. Technol., vol. 6, pp. 89-94, 1998." << endl; 00083 } 00084 00085 // ===================================================================== 00086 // --------------------------------------------------------------------- 00087 // --------------------------------------------------------------------- 00088 // ===================================================================== 00089 00090 int iProjectorIncrementalSiddon::CheckSpecificParameters() 00091 { 00092 // Nothing to check for this projector 00093 ; 00094 // Normal end 00095 return 0; 00096 } 00097 00098 // ===================================================================== 00099 // --------------------------------------------------------------------- 00100 // --------------------------------------------------------------------- 00101 // ===================================================================== 00102 00103 int iProjectorIncrementalSiddon::InitializeSpecific() 00104 { 00105 // Verbose 00106 if (m_verbose>=2) Cout("iProjectorIncrementalSiddon::InitializeSpecific() -> Use incremental Siddon projector" << endl); 00107 // Normal end 00108 return 0; 00109 } 00110 00111 // ===================================================================== 00112 // --------------------------------------------------------------------- 00113 // --------------------------------------------------------------------- 00114 // ===================================================================== 00115 00116 INTNB iProjectorIncrementalSiddon::EstimateMaxNumberOfVoxelsPerLine() 00117 { 00118 // Find the maximum number of voxels along a given dimension 00119 INTNB max_nb_voxels_in_dimension = mp_ImageDimensionsAndQuantification->GetNbVoxX(); 00120 if (mp_ImageDimensionsAndQuantification->GetNbVoxY()>max_nb_voxels_in_dimension) max_nb_voxels_in_dimension = mp_ImageDimensionsAndQuantification->GetNbVoxY(); 00121 if (mp_ImageDimensionsAndQuantification->GetNbVoxZ()>max_nb_voxels_in_dimension) max_nb_voxels_in_dimension = mp_ImageDimensionsAndQuantification->GetNbVoxZ(); 00122 // We should have at most 4 voxels in a given plane, so multiply by 4 00123 // (note: this is not true however it ensures no overflow and is already quite optimized for RAM usage !) 00124 max_nb_voxels_in_dimension *= 4; 00125 // Return the value 00126 return max_nb_voxels_in_dimension; 00127 } 00128 00129 // ===================================================================== 00130 // --------------------------------------------------------------------- 00131 // --------------------------------------------------------------------- 00132 // ===================================================================== 00133 00134 int iProjectorIncrementalSiddon::ProjectWithoutTOF(int a_direction, oProjectionLine* ap_ProjectionLine ) 00135 { 00136 #ifdef CASTOR_DEBUG 00137 if (!m_initialized) 00138 { 00139 Cerr("***** iProjectorIncrementalSiddon::ProjectWithoutTOF() -> Called while not initialized !" << endl); 00140 Exit(EXIT_DEBUG); 00141 } 00142 #endif 00143 00144 #ifdef CASTOR_VERBOSE 00145 if (m_verbose>=10) 00146 { 00147 string direction = ""; 00148 if (a_direction==FORWARD) direction = "forward"; 00149 else direction = "backward"; 00150 Cout("iProjectorIncrementalSiddon::Project without TOF -> Project line '" << ap_ProjectionLine << "' in " << direction << " direction" << endl); 00151 } 00152 #endif 00153 00154 // Simpler way now, hopefully it works 00155 FLTNB* event1Float = ap_ProjectionLine->GetPosition1(); 00156 FLTNB* event2Float = ap_ProjectionLine->GetPosition2(); 00157 double event1[3] = { event1Float[0], event1Float[1], event1Float[2] }; 00158 double event2[3] = { event2Float[0], event2Float[1], event2Float[2] }; 00159 00160 // ************************************** 00161 // STEP 1: LOR length calculation 00162 // ************************************** 00163 double length_LOR = ap_ProjectionLine->GetLength(); 00164 00165 // ************************************** 00166 // STEP 2: Compute entrance voxel indexes 00167 // ************************************** 00168 double alphaFirst[] = { 0.0, 0.0, 0.0 }; 00169 double alphaLast[] = { 0.0, 0.0, 0.0 }; 00170 00171 double alphaMin = 0.0f, alphaMax = 1.0f; 00172 double delta_pos[] = { 00173 event2[ 0 ] - event1[ 0 ], 00174 event2[ 1 ] - event1[ 1 ], 00175 event2[ 2 ] - event1[ 2 ] 00176 }; 00177 00178 // Computation of alphaMin et alphaMax values (entrance and exit point of the ray) 00179 for( int i = 0; i < 3; ++i ) 00180 { 00181 if( delta_pos[ i ] != 0.0 ) 00182 { 00183 alphaFirst[i] = (-mp_halfFOV[i] - event1[i]) / (delta_pos[ i ]); 00184 alphaLast[i] = ( mp_halfFOV[i] - event1[i]) / (delta_pos[ i ]); 00185 alphaMin = (std::max)(alphaMin,(std::min)(alphaFirst[i],alphaLast[i])); 00186 alphaMax = (std::min)(alphaMax,(std::max)(alphaFirst[i],alphaLast[i])); 00187 } 00188 } 00189 00190 // if alphaMax is less than or equal to alphaMin no intersection 00191 // and return an empty buffer 00192 if( alphaMax <= alphaMin ) return 0; 00193 00194 // Now we have to find the indices of the particular plane 00195 // (iMin,iMax), (jMin,jMax), (kMin,kMax) 00196 int iMin = 0, iMax = 0; 00197 int jMin = 0, jMax = 0; 00198 int kMin = 0, kMax = 0; 00199 00200 // For the X-axis 00201 if( delta_pos[ 0 ] > 0.0f ) 00202 { 00203 iMin = ::ceil( ( mp_nbVox[ 0 ] + 1 ) - ( mp_halfFOV[ 0 ] - alphaMin * delta_pos[ 0 ] - event1[ 0 ] ) / mp_sizeVox[0] ); 00204 iMax = ::floor( 1 + ( event1[ 0 ] + alphaMax * delta_pos[ 0 ] - (-mp_halfFOV[ 0 ]) ) / mp_sizeVox[0] ); 00205 } 00206 else if( delta_pos[ 0 ] < 0.0f ) 00207 { 00208 iMin = ::ceil( ( mp_nbVox[ 0 ] + 1 ) - ( mp_halfFOV[ 0 ] - alphaMax * delta_pos[ 0 ] - event1[ 0 ] ) / mp_sizeVox[0] ); 00209 iMax = ::floor( 1 + ( event1[ 0 ] + alphaMin * delta_pos[ 0 ] - (-mp_halfFOV[ 0 ]) ) / mp_sizeVox[0] ); 00210 } 00211 if( delta_pos[ 0 ] == 0 ) 00212 { 00213 iMin = 1, iMax = 0; 00214 } 00215 00216 // For the Y-axis 00217 if( delta_pos[ 1 ] > 0 ) 00218 { 00219 jMin = ::ceil( ( mp_nbVox[ 1 ] + 1 ) - ( mp_halfFOV[ 1 ] - alphaMin * delta_pos[ 1 ] - event1[ 1 ] ) / mp_sizeVox[1] ); 00220 jMax = ::floor( 1 + ( event1[ 1 ] + alphaMax * delta_pos[ 1 ] - (-mp_halfFOV[ 1 ]) ) / mp_sizeVox[1] ); 00221 } 00222 else if( delta_pos[ 1 ] < 0 ) 00223 { 00224 jMin = ::ceil( ( mp_nbVox[ 1 ] + 1 ) - ( mp_halfFOV[ 1 ] - alphaMax * delta_pos[ 1 ] - event1[ 1 ] ) / mp_sizeVox[1] ); 00225 jMax = ::floor( 1 + ( event1[ 1 ] + alphaMin * delta_pos[ 1 ] - (-mp_halfFOV[ 1 ]) ) / mp_sizeVox[1] ); 00226 } 00227 else if( delta_pos[ 1 ] == 0 ) 00228 { 00229 jMin = 1, jMax = 0; 00230 } 00231 00232 // For the Z-axis 00233 if( delta_pos[ 2 ] > 0 ) 00234 { 00235 kMin = ::ceil( ( mp_nbVox[ 2 ] + 1 ) - ( mp_halfFOV[ 2 ] - alphaMin * delta_pos[ 2 ] - event1[ 2 ] ) / mp_sizeVox[2] ); 00236 kMax = ::floor( 1 + ( event1[ 2 ] + alphaMax * delta_pos[ 2 ] - (-mp_halfFOV[ 2 ]) ) / mp_sizeVox[2] ); 00237 } 00238 else if( delta_pos[ 2 ] < 0 ) 00239 { 00240 kMin = ::ceil( ( mp_nbVox[ 2 ] + 1 ) - ( mp_halfFOV[ 2 ] - alphaMax * delta_pos[ 2 ] - event1[ 2 ] ) / mp_sizeVox[2] ); 00241 kMax = ::floor( 1 + ( event1[ 2 ] + alphaMin * delta_pos[ 2 ] - (-mp_halfFOV[ 2 ]) ) / mp_sizeVox[2] ); 00242 } 00243 else if( delta_pos[ 2 ] == 0 ) 00244 { 00245 kMin = 1, kMax = 0; 00246 } 00247 00248 // Computing the last term n number of intersection 00249 INTNB n = ( iMax - iMin + 1 ) + ( jMax - jMin + 1 ) 00250 + ( kMax - kMin + 1 ); 00251 00252 // Array storing the first alpha in X, Y and Z 00253 double alpha_XYZ[ 3 ] = { 1.0f, 1.0f, 1.0f }; 00254 00255 // Computing the first alpha in X 00256 if( delta_pos[ 0 ] > 0 ) 00257 alpha_XYZ[ 0 ] = ( ( (-mp_halfFOV[ 0 ]) + ( iMin - 1 ) * mp_sizeVox[0] ) - event1[ 0 ] ) / delta_pos[ 0 ]; 00258 else if( delta_pos[ 0 ] < 0 ) 00259 alpha_XYZ[ 0 ] = ( ( (-mp_halfFOV[ 0 ]) + ( iMax - 1 ) * mp_sizeVox[0] ) - event1[ 0 ] ) / delta_pos[ 0 ]; 00260 00261 // Computing the first alpha in Y 00262 if( delta_pos[ 1 ] > 0 ) 00263 alpha_XYZ[ 1 ] = ( ( (-mp_halfFOV[ 1 ]) + ( jMin - 1 ) * mp_sizeVox[1] ) - event1[ 1 ] ) / delta_pos[ 1 ]; 00264 else if( delta_pos[ 1 ] < 0 ) 00265 alpha_XYZ[ 1 ] = ( ( (-mp_halfFOV[ 1 ]) + ( jMax - 1 ) * mp_sizeVox[1] ) - event1[ 1 ] ) / delta_pos[ 1 ]; 00266 00267 // Computing the first alpha in Z 00268 if( delta_pos[ 2 ] > 0 ) 00269 alpha_XYZ[ 2 ] = ( ( (-mp_halfFOV[ 2 ]) + ( kMin - 1 ) * mp_sizeVox[2] ) - event1[ 2 ] ) / delta_pos[ 2 ]; 00270 else if( delta_pos[ 2 ] < 0 ) 00271 alpha_XYZ[ 2 ] = ( ( (-mp_halfFOV[ 2 ]) + ( kMax - 1 ) * mp_sizeVox[2] ) - event1[ 2 ] ) / delta_pos[ 2 ]; 00272 00273 // Computing the alpha updating 00274 double alpha_u[ 3 ] = { 00275 mp_sizeVox[0] / std::fabs( delta_pos[ 0 ] ), 00276 mp_sizeVox[1] / std::fabs( delta_pos[ 1 ] ), 00277 mp_sizeVox[2] / std::fabs( delta_pos[ 2 ] ) 00278 }; 00279 00280 // Computing the index updating 00281 INTNB index_u[ 3 ] = { 00282 (delta_pos[ 0 ] < 0) ? -1 : 1, 00283 (delta_pos[ 1 ] < 0) ? -1 : 1, 00284 (delta_pos[ 2 ] < 0) ? -1 : 1 00285 }; 00286 00287 // Check which alpha is the min/max and increment 00288 if( alpha_XYZ[ 0 ] == alphaMin ) alpha_XYZ[ 0 ] += alpha_u[ 0 ]; 00289 if( alpha_XYZ[ 1 ] == alphaMin ) alpha_XYZ[ 1 ] += alpha_u[ 1 ]; 00290 if( alpha_XYZ[ 2 ] == alphaMin ) alpha_XYZ[ 2 ] += alpha_u[ 2 ]; 00291 00292 // Computing the minimum value in the alpha_XYZ buffer 00293 double const min_alpha_XYZ = (std::min)( alpha_XYZ[ 0 ], 00294 (std::min)( alpha_XYZ[ 1 ], alpha_XYZ[ 2 ] ) ); 00295 00296 // Computing the first index of intersection 00297 double const alphaMid = ( min_alpha_XYZ + alphaMin ) / 2.0f; 00298 INTNB index_ijk[ 3 ] = { 00299 1 + (int)( ( event1[ 0 ] + alphaMid * delta_pos[ 0 ] - (-mp_halfFOV[ 0 ]) ) / mp_sizeVox[0] ), 00300 1 + (int)( ( event1[ 1 ] + alphaMid * delta_pos[ 1 ] - (-mp_halfFOV[ 1 ]) ) / mp_sizeVox[1] ), 00301 1 + (int)( ( event1[ 2 ] + alphaMid * delta_pos[ 2 ] - (-mp_halfFOV[ 2 ]) ) / mp_sizeVox[2] ) 00302 }; 00303 00304 INTNB const w = mp_nbVox[ 0 ]; 00305 INTNB const wh = w * mp_nbVox[ 1 ]; 00306 00307 // Loop over the number of plane to cross 00308 double alpha_c = alphaMin; 00309 FLTNB coeff = 0.0f; 00310 INTNB numVox = 0; 00311 for( int nP = 0; nP < n - 1; ++nP ) 00312 { 00313 if( ( alpha_XYZ[ 0 ] <= alpha_XYZ[ 1 ] ) 00314 && ( alpha_XYZ[ 0 ] <= alpha_XYZ[ 2 ] ) ) 00315 { 00316 // Storing values 00317 if( ( alpha_XYZ[ 0 ] >= alphaMin ) 00318 && ( index_ijk[ 0 ] - 1 >= 0 ) 00319 && ( index_ijk[ 0 ] - 1 <= mp_nbVox[ 0 ] - 1 ) 00320 && ( index_ijk[ 1 ] - 1 >= 0 ) 00321 && ( index_ijk[ 1 ] - 1 <= mp_nbVox[ 1 ] - 1 ) 00322 && ( index_ijk[ 2 ] - 1 >= 0 ) 00323 && ( index_ijk[ 2 ] - 1 <= mp_nbVox[ 2 ] - 1 ) ) 00324 { 00325 coeff = ( alpha_XYZ[ 0 ] - alpha_c ) * length_LOR; 00326 numVox = ( index_ijk[ 0 ] - 1 ) + ( ( index_ijk[ 1 ] - 1 ) ) * w + ( ( index_ijk[ 2 ] - 1 ) ) * wh; 00327 ap_ProjectionLine->AddVoxel(a_direction, numVox, coeff); 00328 } 00329 00330 // Increment values 00331 alpha_c = alpha_XYZ[ 0 ]; 00332 alpha_XYZ[ 0 ] += alpha_u[ 0 ]; 00333 index_ijk[ 0 ] += index_u[ 0 ]; 00334 } 00335 else if( ( alpha_XYZ[ 1 ] < alpha_XYZ[ 0 ] ) 00336 && ( alpha_XYZ[ 1 ] <= alpha_XYZ[ 2 ] ) ) 00337 { 00338 // Storing values 00339 if( ( alpha_XYZ[ 1 ] >= alphaMin ) 00340 && ( index_ijk[ 0 ] - 1 >= 0 ) 00341 && ( index_ijk[ 0 ] - 1 <= mp_nbVox[ 0 ] - 1 ) 00342 && ( index_ijk[ 1 ] - 1 >= 0 ) 00343 && ( index_ijk[ 1 ] - 1 <= mp_nbVox[ 1 ] - 1 ) 00344 && ( index_ijk[ 2 ] - 1 >= 0 ) 00345 && ( index_ijk[ 2 ] - 1 <= mp_nbVox[ 2 ] - 1 ) ) 00346 { 00347 coeff = ( alpha_XYZ[ 1 ] - alpha_c ) * length_LOR; 00348 numVox = ( index_ijk[ 0 ] - 1 ) + ( ( index_ijk[ 1 ] - 1 ) ) * w + ( ( index_ijk[ 2 ] - 1 ) ) * wh; 00349 ap_ProjectionLine->AddVoxel(a_direction, numVox, coeff); 00350 } 00351 00352 // Increment values 00353 alpha_c = alpha_XYZ[ 1 ]; 00354 alpha_XYZ[ 1 ] += alpha_u[ 1 ]; 00355 index_ijk[ 1 ] += index_u[ 1 ]; 00356 } 00357 else if( ( alpha_XYZ[ 2 ] < alpha_XYZ[ 0 ] ) 00358 && ( alpha_XYZ[ 2 ] < alpha_XYZ[ 1 ] ) ) 00359 { 00360 // Storing values 00361 if( ( alpha_XYZ[ 2 ] >= alphaMin ) 00362 && ( index_ijk[ 0 ] - 1 >= 0 ) 00363 && ( index_ijk[ 0 ] - 1 <= mp_nbVox[ 0 ] - 1 ) 00364 && ( index_ijk[ 1 ] - 1 >= 0 ) 00365 && ( index_ijk[ 1 ] - 1 <= mp_nbVox[ 1 ] - 1 ) 00366 && ( index_ijk[ 2 ] - 1 >= 0 ) 00367 && ( index_ijk[ 2 ] - 1 <= mp_nbVox[ 2 ] - 1 ) ) 00368 { 00369 coeff = ( alpha_XYZ[ 2 ] - alpha_c ) * length_LOR; 00370 numVox = ( index_ijk[ 0 ] - 1 ) + ( ( index_ijk[ 1 ] - 1 ) ) * w + ( ( index_ijk[ 2 ] - 1 ) ) * wh; 00371 ap_ProjectionLine->AddVoxel(a_direction, numVox, coeff); 00372 } 00373 00374 // Increment values 00375 alpha_c = alpha_XYZ[ 2 ]; 00376 alpha_XYZ[ 2 ] += alpha_u[ 2 ]; 00377 index_ijk[ 2 ] += index_u[ 2 ]; 00378 } 00379 } 00380 00381 return 0; 00382 } 00383 00384 // ===================================================================== 00385 // --------------------------------------------------------------------- 00386 // --------------------------------------------------------------------- 00387 // ===================================================================== 00388 00389 int iProjectorIncrementalSiddon::ProjectWithTOFPos(int a_Projector, oProjectionLine* ap_ProjectionLine) 00390 { 00391 Cerr("***** iProjectorIncrementalSiddon::ProjectWithTOFPos() -> Not yet implemented !" << endl); 00392 return 1; 00393 } 00394 00395 // ===================================================================== 00396 // --------------------------------------------------------------------- 00397 // --------------------------------------------------------------------- 00398 // ===================================================================== 00399 00400 int iProjectorIncrementalSiddon::ProjectWithTOFBin(int a_Projector, oProjectionLine* ap_ProjectionLine) 00401 { 00402 Cerr("***** iProjectorIncrementalSiddon::ProjectWithTOFBin() -> Not yet implemented !" << endl); 00403 return 1; 00404 } 00405 00406 // ===================================================================== 00407 // --------------------------------------------------------------------- 00408 // --------------------------------------------------------------------- 00409 // =====================================================================