![]() |
CASToR
1.0
Tomographic Reconstruction (PET/SPECT)
|
00001 00002 /* 00003 Implementation of class iProjectorClassicSiddon 00004 00005 - separators: done 00006 - doxygen: done 00007 - default initialization: done 00008 - CASTOR_DEBUG: 00009 - CASTOR_VERBOSE: 00010 */ 00011 00018 #include <algorithm> 00019 #include <cstring> 00020 00021 #include "iProjectorClassicSiddon.hh" 00022 #include "sOutputManager.hh" 00023 00024 // ===================================================================== 00025 // --------------------------------------------------------------------- 00026 // --------------------------------------------------------------------- 00027 // ===================================================================== 00028 00029 iProjectorClassicSiddon::iProjectorClassicSiddon() : vProjector() 00030 { 00031 // This projector is compatible with SPECT attenuation correction because 00032 // all voxels contributing to a given line are ordered from point 1 (focal) 00033 // to point 2 (scanner) 00034 m_compatibleWithSPECTAttenuationCorrection = true; 00035 } 00036 00037 // ===================================================================== 00038 // --------------------------------------------------------------------- 00039 // --------------------------------------------------------------------- 00040 // ===================================================================== 00041 00042 iProjectorClassicSiddon::~iProjectorClassicSiddon() 00043 { 00044 } 00045 00046 // ===================================================================== 00047 // --------------------------------------------------------------------- 00048 // --------------------------------------------------------------------- 00049 // ===================================================================== 00050 00051 int iProjectorClassicSiddon::ReadConfigurationFile(const string& a_configurationFile) 00052 { 00053 // No options for classic siddon 00054 ; 00055 // Normal end 00056 return 0; 00057 } 00058 00059 // ===================================================================== 00060 // --------------------------------------------------------------------- 00061 // --------------------------------------------------------------------- 00062 // ===================================================================== 00063 00064 int iProjectorClassicSiddon::ReadOptionsList(const string& a_optionsList) 00065 { 00066 // No options for classic siddon 00067 ; 00068 // Normal end 00069 return 0; 00070 } 00071 00072 // ===================================================================== 00073 // --------------------------------------------------------------------- 00074 // --------------------------------------------------------------------- 00075 // ===================================================================== 00076 00077 void iProjectorClassicSiddon::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 implemented from the following published paper:" << endl; 00081 cout << "R. L. Siddon, \"Fast calculation of the exact radiological path for a three-dimensional CT array\", Med. Phys., vol. 12, pp. 252-5, 1985." << endl; 00082 cout << "All voxels are correctly ordered in the line, so this projector can be used with SPECT attenuation correction." << endl; 00083 } 00084 00085 // ===================================================================== 00086 // --------------------------------------------------------------------- 00087 // --------------------------------------------------------------------- 00088 // ===================================================================== 00089 00090 int iProjectorClassicSiddon::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 iProjectorClassicSiddon::InitializeSpecific() 00104 { 00105 // Verbose 00106 if (m_verbose>=2) Cout("iProjectorClassiclSiddon::InitializeSpecific() -> Use classic Siddon projector" << endl); 00107 // Normal end 00108 return 0; 00109 } 00110 00111 // ===================================================================== 00112 // --------------------------------------------------------------------- 00113 // --------------------------------------------------------------------- 00114 // ===================================================================== 00115 00116 INTNB iProjectorClassicSiddon::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 iProjectorClassicSiddon::ProjectWithoutTOF(int a_direction, oProjectionLine* ap_ProjectionLine ) 00135 { 00136 #ifdef CASTOR_DEBUG 00137 if (!m_initialized) 00138 { 00139 Cerr("***** iProjectorClassicSiddon::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("iProjectorClassicSiddon::Project without TOF -> Project line '" << ap_ProjectionLine << "' in " << direction << " direction" << endl); 00151 } 00152 #endif 00153 00154 // Simpler way now, hopefully it works 00155 FLTNB* event1_position = ap_ProjectionLine->GetPosition1(); 00156 FLTNB* event2_position = ap_ProjectionLine->GetPosition2(); 00157 /*ORIG CODE 00158 long double event1[3] = { event1_position[0], event1_position[1], event1_position[2] }; 00159 long double event2[3] = { event2_position[0], event2_position[1], event2_position[2] }; 00160 00161 // ************************************** 00162 // STEP 1: LOR length calculation 00163 // ************************************** 00164 long double length_LOR = ap_ProjectionLine->GetLength(); 00165 00166 // ************************************** 00167 // STEP 2: Compute entrance voxel indexes 00168 // ************************************** 00169 long double alphaFirst[] = { 0.0, 0.0, 0.0 }; 00170 long double alphaLast[] = { 0.0, 0.0, 0.0 }; 00171 00172 long double alphaMin = 0.0, alphaMax = 1.0; 00173 long double delta_pos[] = { 00174 event2[ 0 ] - event1[ 0 ], 00175 event2[ 1 ] - event1[ 1 ], 00176 event2[ 2 ] - event1[ 2 ] 00177 }; 00178 * */ 00179 00180 FLTNB event1[3] = { event1_position[0], event1_position[1], event1_position[2] }; 00181 FLTNB event2[3] = { event2_position[0], event2_position[1], event2_position[2] }; 00182 00183 // ************************************** 00184 // STEP 1: LOR length calculation 00185 // ************************************** 00186 FLTNB length_LOR = ap_ProjectionLine->GetLength(); 00187 00188 // ************************************** 00189 // STEP 2: Compute entrance voxel indexes 00190 // ************************************** 00191 FLTNB alphaFirst[] = { 0.0, 0.0, 0.0 }; 00192 FLTNB alphaLast[] = { 0.0, 0.0, 0.0 }; 00193 00194 FLTNB alphaMin = 0.0, alphaMax = 1.0; 00195 FLTNB delta_pos[] = { 00196 event2[ 0 ] - event1[ 0 ], 00197 event2[ 1 ] - event1[ 1 ], 00198 event2[ 2 ] - event1[ 2 ] 00199 }; 00200 00201 // Computation of alphaMin et alphaMax values (entrance and exit point of the ray) 00202 for( int i = 0; i < 3; ++i ) 00203 { 00204 if( delta_pos[ i ] != 0.0 ) 00205 { 00206 alphaFirst[i] = (-mp_halfFOV[i] - event1[i]) / (delta_pos[ i ]); 00207 alphaLast[i] = ( mp_halfFOV[i] - event1[i]) / (delta_pos[ i ]); 00208 alphaMin = (std::max)(alphaMin,(std::min)(alphaFirst[i],alphaLast[i])); 00209 alphaMax = (std::min)(alphaMax,(std::max)(alphaFirst[i],alphaLast[i])); 00210 } 00211 } 00212 00213 // if alphaMax is less than or equal to alphaMin no intersection 00214 // and return an empty buffer 00215 if( alphaMax <= alphaMin ) return 0; 00216 00217 // Now we have to find the indices of the particular plane 00218 // (iMin,iMax), (jMin,jMax), (kMin,kMax) 00219 int iMin = 0, iMax = 0; 00220 int jMin = 0, jMax = 0; 00221 int kMin = 0, kMax = 0; 00222 00223 // For the X-axis 00224 if( delta_pos[ 0 ] > 0.0 ) 00225 { 00226 iMin = ::ceil( ( mp_nbVox[ 0 ] + 1 ) - ( mp_halfFOV[ 0 ] - alphaMin * delta_pos[ 0 ] - event1[ 0 ] ) / mp_sizeVox[0] ); 00227 iMax = ::floor( 1 + ( event1[ 0 ] + alphaMax * delta_pos[ 0 ] - (-mp_halfFOV[ 0 ]) ) / mp_sizeVox[0] ); 00228 } 00229 else if( delta_pos[ 0 ] < 0.0 ) 00230 { 00231 iMin = ::ceil( ( mp_nbVox[ 0 ] + 1 ) - ( mp_halfFOV[ 0 ] - alphaMax * delta_pos[ 0 ] - event1[ 0 ] ) / mp_sizeVox[0] ); 00232 iMax = ::floor( 1 + ( event1[ 0 ] + alphaMin * delta_pos[ 0 ] - (-mp_halfFOV[ 0 ]) ) / mp_sizeVox[0] ); 00233 } 00234 if( delta_pos[ 0 ] == 0. ) 00235 { 00236 iMin = 1, iMax = 0; 00237 } 00238 00239 // For the Y-axis 00240 if( delta_pos[ 1 ] > 0. ) 00241 { 00242 jMin = ::ceil( ( mp_nbVox[ 1 ] + 1 ) - ( mp_halfFOV[ 1 ] - alphaMin * delta_pos[ 1 ] - event1[ 1 ] ) / mp_sizeVox[1] ); 00243 jMax = ::floor( 1 + ( event1[ 1 ] + alphaMax * delta_pos[ 1 ] - (-mp_halfFOV[ 1 ]) ) / mp_sizeVox[1] ); 00244 } 00245 else if( delta_pos[ 1 ] < 0. ) 00246 { 00247 jMin = ::ceil( ( mp_nbVox[ 1 ] + 1 ) - ( mp_halfFOV[ 1 ] - alphaMax * delta_pos[ 1 ] - event1[ 1 ] ) / mp_sizeVox[1] ); 00248 jMax = ::floor( 1 + ( event1[ 1 ] + alphaMin * delta_pos[ 1 ] - (-mp_halfFOV[ 1 ]) ) / mp_sizeVox[1] ); 00249 } 00250 else if( delta_pos[ 1 ] == 0. ) 00251 { 00252 jMin = 1, jMax = 0; 00253 } 00254 00255 // For the Z-axis 00256 if( delta_pos[ 2 ] > 0. ) 00257 { 00258 kMin = ::ceil( ( mp_nbVox[ 2 ] + 1 ) - ( mp_halfFOV[ 2 ] - alphaMin * delta_pos[ 2 ] - event1[ 2 ] ) / mp_sizeVox[2] ); 00259 kMax = ::floor( 1 + ( event1[ 2 ] + alphaMax * delta_pos[ 2 ] - (-mp_halfFOV[ 2 ]) ) / mp_sizeVox[2] ); 00260 } 00261 else if( delta_pos[ 2 ] < 0. ) 00262 { 00263 kMin = ::ceil( ( mp_nbVox[ 2 ] + 1 ) - ( mp_halfFOV[ 2 ] - alphaMax * delta_pos[ 2 ] - event1[ 2 ] ) / mp_sizeVox[2] ); 00264 kMax = ::floor( 1 + ( event1[ 2 ] + alphaMin * delta_pos[ 2 ] - (-mp_halfFOV[ 2 ]) ) / mp_sizeVox[2] ); 00265 } 00266 else if( delta_pos[ 2 ] == 0. ) 00267 { 00268 kMin = 1, kMax = 0; 00269 } 00270 00271 // Computing the last term n number of intersection 00272 int n = ( iMax - iMin + 1 ) + ( jMax - jMin + 1 ) 00273 + ( kMax - kMin + 1 ); 00274 00275 // We create a buffer storing the merging data 00276 // We merge alphaMin, alphaMax, alphaX, alphaY and alphaZ 00277 FLTNB *alpha = new FLTNB[ mp_nbVox[ 0 ] + mp_nbVox[ 1 ] + mp_nbVox[ 2 ] + 3 ]; 00278 FLTNB *tmpAlpha = new FLTNB[ mp_nbVox[ 0 ] + mp_nbVox[ 1 ] + mp_nbVox[ 2 ] + 3 ]; 00279 FLTNB *alphaX = new FLTNB[ ( mp_nbVox[ 0 ] + 1 ) ]; 00280 FLTNB *alphaY = new FLTNB[ ( mp_nbVox[ 1 ] + 1 ) ]; 00281 FLTNB *alphaZ = new FLTNB[ ( mp_nbVox[ 2 ] + 1 ) ]; 00282 00283 INTNB iElement = iMax - iMin + 1; 00284 if( iElement > 0 ) 00285 { 00286 FLTNB *idx = alphaX; 00287 if( delta_pos[ 0 ] > 0. ) 00288 { 00289 for( int i = iMin; i <= iMax; ++i ) 00290 *idx++ = ( ( (-mp_halfFOV[ 0 ]) + ( i - 1 ) * mp_sizeVox[0] ) - event1[ 0 ] ) / delta_pos[ 0 ]; 00291 } 00292 else if( delta_pos[ 0 ] < 0. ) 00293 { 00294 for( int i = iMax; i >= iMin; --i ) 00295 *idx++ = ( ( (-mp_halfFOV[ 0 ]) + ( i - 1 ) * mp_sizeVox[0] ) - event1[ 0 ] ) / delta_pos[ 0 ]; 00296 } 00297 } 00298 00299 // For alphaY 00300 INTNB jElement = jMax - jMin + 1; 00301 if( jElement > 0 ) 00302 { 00303 FLTNB *idx = alphaY; 00304 if( delta_pos[ 1 ] > 0. ) 00305 { 00306 for( int j = jMin; j <= jMax; ++j ) 00307 *idx++ = ( ( (-mp_halfFOV[ 1 ]) + ( j - 1 ) * mp_sizeVox[1] ) - event1[ 1 ] ) / delta_pos[ 1 ]; 00308 } 00309 else if( delta_pos[ 1 ] < 0. ) 00310 { 00311 for( int j = jMax; j >= jMin; --j ) 00312 *idx++ = ( ( (-mp_halfFOV[ 1 ]) + ( j - 1 ) * mp_sizeVox[1] ) - event1[ 1 ] ) / delta_pos[ 1 ]; 00313 } 00314 } 00315 00316 // For alphaZ 00317 INTNB kElement = kMax - kMin + 1; 00318 if( kElement > 0 ) 00319 { 00320 FLTNB *idx = alphaZ; 00321 if( delta_pos[ 2 ] > 0. ) 00322 { 00323 for( int k = kMin; k <= kMax; ++k ) 00324 *idx++ = ( ( (-mp_halfFOV[ 2 ]) + ( k - 1 ) * mp_sizeVox[2] ) - event1[ 2 ] ) / delta_pos[ 2 ]; 00325 } 00326 else if( delta_pos[ 2 ] < 0. ) 00327 { 00328 for( int k = kMax; k >= kMin; --k ) 00329 *idx++ = ( ( (-mp_halfFOV[ 2 ]) + ( k - 1 ) * mp_sizeVox[2] ) - event1[ 2 ] ) / delta_pos[ 2 ]; 00330 } 00331 } 00332 00333 std::merge( 00334 alphaX, alphaX + iElement, 00335 tmpAlpha, tmpAlpha, 00336 alpha ); 00337 00338 ::memcpy( tmpAlpha, alpha, ( iElement ) * sizeof( FLTNB ) ); 00339 00340 std::merge( 00341 alphaY, alphaY + jElement, 00342 tmpAlpha, tmpAlpha + iElement, 00343 alpha ); 00344 00345 ::memcpy( tmpAlpha, alpha, ( iElement + jElement ) * sizeof( FLTNB ) ); 00346 00347 std::merge( 00348 alphaZ, alphaZ + kElement, 00349 tmpAlpha, tmpAlpha + iElement + jElement, 00350 alpha ); 00351 00352 // Computing the index of the voxels 00353 FLTNB alphaMid = 0.0; 00354 INTNB i = 0, j = 0, k = 0; // indices of the voxel 00355 FLTNB *pAlpha = &alpha[ 1 ]; 00356 FLTNB *pAlphaPrevious = &alpha[ 0 ]; 00357 FLTNB coeff = 0.0; 00358 INTNB numVox = 0; 00359 // Loop over the number of crossed planes 00360 for( int nP = 0; nP < n - 1; ++nP, ++pAlpha, ++pAlphaPrevious ) 00361 { 00362 alphaMid = ( *pAlpha + *pAlphaPrevious ) * 0.5; 00363 00364 i = 1 + ( event1[ 0 ] + alphaMid * delta_pos[ 0 ] - (-mp_halfFOV[ 0 ]) ) / mp_sizeVox[0]; 00365 if( i < 1 || i > mp_nbVox[ 0 ] ) continue; 00366 00367 j = 1 + ( event1[ 1 ] + alphaMid * delta_pos[ 1 ] - (-mp_halfFOV[ 1 ]) ) / mp_sizeVox[1]; 00368 if( j < 1 || j > mp_nbVox[ 1 ] ) continue; 00369 00370 k = 1 + ( event1[ 2 ] + alphaMid * delta_pos[ 2 ] - (-mp_halfFOV[ 2 ]) ) / mp_sizeVox[2]; 00371 if( k < 1 || k > mp_nbVox[ 2 ] ) continue; 00372 00373 numVox = ( i - 1 ) + ( j - 1 ) * mp_nbVox[0] + ( k - 1 ) * mp_nbVox[0] * mp_nbVox[1]; 00374 coeff = length_LOR * ( *pAlpha - *pAlphaPrevious ); 00375 00376 ap_ProjectionLine->AddVoxel(a_direction, numVox, coeff); 00377 } 00378 00379 delete[] alpha; 00380 delete[] tmpAlpha; 00381 delete[] alphaX; 00382 delete[] alphaY; 00383 delete[] alphaZ; 00384 00385 #ifdef CASTOR_VERBOSE 00386 if (m_verbose>=10) 00387 { 00388 Cout("iProjectorClassicSiddon::Project without TOF -> Exit function" << endl); 00389 } 00390 #endif 00391 00392 return 0; 00393 } 00394 00395 // ===================================================================== 00396 // --------------------------------------------------------------------- 00397 // --------------------------------------------------------------------- 00398 // ===================================================================== 00399 00400 int iProjectorClassicSiddon::ProjectWithTOFPos(int a_Projector, oProjectionLine* ap_ProjectionLine) 00401 { 00402 Cerr("***** iProjectorClassicSiddon::ProjectWithTOFPos() -> Not yet implemented !" << endl); 00403 return 1; 00404 } 00405 00406 // ===================================================================== 00407 // --------------------------------------------------------------------- 00408 // --------------------------------------------------------------------- 00409 // ===================================================================== 00410 00411 int iProjectorClassicSiddon::ProjectWithTOFBin(int a_Projector, oProjectionLine* ap_ProjectionLine) 00412 { 00413 Cerr("***** iProjectorClassicSiddon::ProjectWithTOFBin() -> Not yet implemented !" << endl); 00414 return 1; 00415 } 00416 00417 // ===================================================================== 00418 // --------------------------------------------------------------------- 00419 // --------------------------------------------------------------------- 00420 // =====================================================================