100 cout <<
"This projector is a simple line projector that computes the exact path length of a line through the voxels." << endl;
101 cout <<
"It is basically an incremental version of the original algorithm proposed by R. L. Siddon (see the classicSiddon projector)." << endl;
102 cout <<
"It is implemented from the following published paper:" << endl;
103 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;
127 if (
m_verbose>=2)
Cout(
"iProjectorIncrementalSiddon::InitializeSpecific() -> Use incremental Siddon projector" << endl);
145 max_nb_voxels_in_dimension *= 4;
147 return max_nb_voxels_in_dimension;
160 Cerr(
"***** iProjectorIncrementalSiddon::ProjectWithoutTOF() -> Called while not initialized !" << endl);
165 #ifdef CASTOR_VERBOSE 168 string direction =
"";
169 if (a_direction==
FORWARD) direction =
"forward";
170 else direction =
"backward";
171 Cout(
"iProjectorIncrementalSiddon::Project without TOF -> Project line '" << ap_ProjectionLine <<
"' in " << direction <<
" direction" << endl);
178 HPFLTNB event1[3] = { event1Float[0], event1Float[1], event1Float[2] };
179 HPFLTNB event2[3] = { event2Float[0], event2Float[1], event2Float[2] };
189 HPFLTNB alphaFirst[] = { 0.0, 0.0, 0.0 };
190 HPFLTNB alphaLast[] = { 0.0, 0.0, 0.0 };
192 HPFLTNB alphaMin = 0.0f, alphaMax = 1.0f;
194 event2[ 0 ] - event1[ 0 ],
195 event2[ 1 ] - event1[ 1 ],
196 event2[ 2 ] - event1[ 2 ]
200 for(
int i = 0; i < 3; ++i )
202 if( delta_pos[ i ] != 0.0 )
204 alphaFirst[i] = (-
mp_halfFOV[i] - event1[i]) / (delta_pos[ i ]);
205 alphaLast[i] = (
mp_halfFOV[i] - event1[i]) / (delta_pos[ i ]);
206 alphaMin = (std::max)(alphaMin,(std::min)(alphaFirst[i],alphaLast[i]));
207 alphaMax = (std::min)(alphaMax,(std::max)(alphaFirst[i],alphaLast[i]));
213 if( alphaMax <= alphaMin )
return 0;
217 int iMin = 0, iMax = 0;
218 int jMin = 0, jMax = 0;
219 int kMin = 0, kMax = 0;
222 if( delta_pos[ 0 ] > 0.0f )
225 iMax = ::floor( 1 + ( event1[ 0 ] + alphaMax * delta_pos[ 0 ] - (-
mp_halfFOV[ 0 ]) ) /
mp_sizeVox[0] );
227 else if( delta_pos[ 0 ] < 0.0f )
230 iMax = ::floor( 1 + ( event1[ 0 ] + alphaMin * delta_pos[ 0 ] - (-
mp_halfFOV[ 0 ]) ) /
mp_sizeVox[0] );
232 if( delta_pos[ 0 ] == 0 )
238 if( delta_pos[ 1 ] > 0 )
241 jMax = ::floor( 1 + ( event1[ 1 ] + alphaMax * delta_pos[ 1 ] - (-
mp_halfFOV[ 1 ]) ) /
mp_sizeVox[1] );
243 else if( delta_pos[ 1 ] < 0 )
246 jMax = ::floor( 1 + ( event1[ 1 ] + alphaMin * delta_pos[ 1 ] - (-
mp_halfFOV[ 1 ]) ) /
mp_sizeVox[1] );
248 else if( delta_pos[ 1 ] == 0 )
254 if( delta_pos[ 2 ] > 0 )
257 kMax = ::floor( 1 + ( event1[ 2 ] + alphaMax * delta_pos[ 2 ] - (-
mp_halfFOV[ 2 ]) ) /
mp_sizeVox[2] );
259 else if( delta_pos[ 2 ] < 0 )
262 kMax = ::floor( 1 + ( event1[ 2 ] + alphaMin * delta_pos[ 2 ] - (-
mp_halfFOV[ 2 ]) ) /
mp_sizeVox[2] );
264 else if( delta_pos[ 2 ] == 0 )
270 INTNB n = ( iMax - iMin + 1 ) + ( jMax - jMin + 1 )
271 + ( kMax - kMin + 1 );
274 HPFLTNB alpha_XYZ[ 3 ] = { 1.0f, 1.0f, 1.0f };
277 if( delta_pos[ 0 ] > 0 )
278 alpha_XYZ[ 0 ] = ( ( (-
mp_halfFOV[ 0 ]) + ( iMin - 1 ) *
mp_sizeVox[0] ) - event1[ 0 ] ) / delta_pos[ 0 ];
279 else if( delta_pos[ 0 ] < 0 )
280 alpha_XYZ[ 0 ] = ( ( (-
mp_halfFOV[ 0 ]) + ( iMax - 1 ) *
mp_sizeVox[0] ) - event1[ 0 ] ) / delta_pos[ 0 ];
283 if( delta_pos[ 1 ] > 0 )
284 alpha_XYZ[ 1 ] = ( ( (-
mp_halfFOV[ 1 ]) + ( jMin - 1 ) *
mp_sizeVox[1] ) - event1[ 1 ] ) / delta_pos[ 1 ];
285 else if( delta_pos[ 1 ] < 0 )
286 alpha_XYZ[ 1 ] = ( ( (-
mp_halfFOV[ 1 ]) + ( jMax - 1 ) *
mp_sizeVox[1] ) - event1[ 1 ] ) / delta_pos[ 1 ];
289 if( delta_pos[ 2 ] > 0 )
290 alpha_XYZ[ 2 ] = ( ( (-
mp_halfFOV[ 2 ]) + ( kMin - 1 ) *
mp_sizeVox[2] ) - event1[ 2 ] ) / delta_pos[ 2 ];
291 else if( delta_pos[ 2 ] < 0 )
292 alpha_XYZ[ 2 ] = ( ( (-
mp_halfFOV[ 2 ]) + ( kMax - 1 ) *
mp_sizeVox[2] ) - event1[ 2 ] ) / delta_pos[ 2 ];
302 INTNB index_u[ 3 ] = {
303 (delta_pos[ 0 ] < 0) ? -1 : 1,
304 (delta_pos[ 1 ] < 0) ? -1 : 1,
305 (delta_pos[ 2 ] < 0) ? -1 : 1
309 if( alpha_XYZ[ 0 ] == alphaMin ) alpha_XYZ[ 0 ] += alpha_u[ 0 ];
310 if( alpha_XYZ[ 1 ] == alphaMin ) alpha_XYZ[ 1 ] += alpha_u[ 1 ];
311 if( alpha_XYZ[ 2 ] == alphaMin ) alpha_XYZ[ 2 ] += alpha_u[ 2 ];
314 HPFLTNB const min_alpha_XYZ = (std::min)( alpha_XYZ[ 0 ],
315 (std::min)( alpha_XYZ[ 1 ], alpha_XYZ[ 2 ] ) );
318 HPFLTNB alphaMid = ( min_alpha_XYZ + alphaMin ) / 2.0f;
319 INTNB index_ijk[ 3 ] = {
333 for(
int nP = 0; nP < n - 1; ++nP )
335 if( ( alpha_XYZ[ 0 ] <= alpha_XYZ[ 1 ] )
336 && ( alpha_XYZ[ 0 ] <= alpha_XYZ[ 2 ] ) )
339 if( ( alpha_XYZ[ 0 ] >= alphaMin )
340 && ( index_ijk[ 0 ] - 1 >= 0 )
341 && ( index_ijk[ 0 ] - 1 <= mp_nbVox[ 0 ] - 1 )
342 && ( index_ijk[ 1 ] - 1 >= 0 )
343 && ( index_ijk[ 1 ] - 1 <= mp_nbVox[ 1 ] - 1 )
344 && ( index_ijk[ 2 ] - 1 >= 0 )
345 && ( index_ijk[ 2 ] - 1 <= mp_nbVox[ 2 ] - 1 ) )
347 coeff = ( alpha_XYZ[ 0 ] - alpha_c ) * length_LOR;
348 numVox = ( index_ijk[ 0 ] - 1 ) + ( ( index_ijk[ 1 ] - 1 ) ) * w + ( ( index_ijk[ 2 ] - 1 ) ) * wh;
353 alpha_c = alpha_XYZ[ 0 ];
354 alpha_XYZ[ 0 ] += alpha_u[ 0 ];
355 index_ijk[ 0 ] += index_u[ 0 ];
357 else if( ( alpha_XYZ[ 1 ] < alpha_XYZ[ 0 ] )
358 && ( alpha_XYZ[ 1 ] <= alpha_XYZ[ 2 ] ) )
361 if( ( alpha_XYZ[ 1 ] >= alphaMin )
362 && ( index_ijk[ 0 ] - 1 >= 0 )
363 && ( index_ijk[ 0 ] - 1 <= mp_nbVox[ 0 ] - 1 )
364 && ( index_ijk[ 1 ] - 1 >= 0 )
365 && ( index_ijk[ 1 ] - 1 <= mp_nbVox[ 1 ] - 1 )
366 && ( index_ijk[ 2 ] - 1 >= 0 )
367 && ( index_ijk[ 2 ] - 1 <= mp_nbVox[ 2 ] - 1 ) )
369 coeff = ( alpha_XYZ[ 1 ] - alpha_c ) * length_LOR;
370 numVox = ( index_ijk[ 0 ] - 1 ) + ( ( index_ijk[ 1 ] - 1 ) ) * w + ( ( index_ijk[ 2 ] - 1 ) ) * wh;
375 alpha_c = alpha_XYZ[ 1 ];
376 alpha_XYZ[ 1 ] += alpha_u[ 1 ];
377 index_ijk[ 1 ] += index_u[ 1 ];
379 else if( ( alpha_XYZ[ 2 ] < alpha_XYZ[ 0 ] )
380 && ( alpha_XYZ[ 2 ] < alpha_XYZ[ 1 ] ) )
383 if( ( alpha_XYZ[ 2 ] >= alphaMin )
384 && ( index_ijk[ 0 ] - 1 >= 0 )
385 && ( index_ijk[ 0 ] - 1 <= mp_nbVox[ 0 ] - 1 )
386 && ( index_ijk[ 1 ] - 1 >= 0 )
387 && ( index_ijk[ 1 ] - 1 <= mp_nbVox[ 1 ] - 1 )
388 && ( index_ijk[ 2 ] - 1 >= 0 )
389 && ( index_ijk[ 2 ] - 1 <= mp_nbVox[ 2 ] - 1 ) )
391 coeff = ( alpha_XYZ[ 2 ] - alpha_c ) * length_LOR;
392 numVox = ( index_ijk[ 0 ] - 1 ) + ( ( index_ijk[ 1 ] - 1 ) ) * w + ( ( index_ijk[ 2 ] - 1 ) ) * wh;
397 alpha_c = alpha_XYZ[ 2 ];
398 alpha_XYZ[ 2 ] += alpha_u[ 2 ];
399 index_ijk[ 2 ] += index_u[ 2 ];
416 Cerr(
"***** iProjectorIncrementalSiddon::ProjectTOFListmode() -> Called while not initialized !" << endl);
421 #ifdef CASTOR_VERBOSE 424 string direction =
"";
425 if (a_direction==
FORWARD) direction =
"forward";
426 else direction =
"backward";
427 Cout(
"iProjectorIncrementalSiddon::Project with TOF position -> Project line '" << ap_ProjectionLine <<
"' in " << direction <<
" direction" << endl);
434 HPFLTNB event1[3] = { event1Float[0], event1Float[1], event1Float[2] };
435 HPFLTNB event2[3] = { event2Float[0], event2Float[1], event2Float[2] };
445 HPFLTNB alphaFirst[] = { 0.0, 0.0, 0.0 };
446 HPFLTNB alphaLast[] = { 0.0, 0.0, 0.0 };
448 HPFLTNB alphaMin = 0.0f, alphaMax = 1.0f;
450 event2[ 0 ] - event1[ 0 ],
451 event2[ 1 ] - event1[ 1 ],
452 event2[ 2 ] - event1[ 2 ]
463 HPFLTNB tof_sigma_sqrt2 = sqrt(2.)*tof_sigma;
471 HPFLTNB lor_tof_center = length_LOR * 0.5 + tof_delta;
474 HPFLTNB tof_edge_low[] = {0,0,0};
476 HPFLTNB tof_edge_high[] = {0,0,0};
478 INTNB tof_index_low[] = {0,0,0};
479 INTNB tof_index_high[] = {0,0,0};
482 for (
int ax=0;ax<3;ax++)
485 tof_center = event1[ax] + lor_tof_center * delta_pos[ax] / length_LOR;
488 tof_edge_low[ax] = tof_center - tof_half_span * fabs(delta_pos[ax]) / length_LOR;
491 if (tof_index_low[ax]>
mp_nbVox[ax]-1)
return 0;
495 tof_edge_high[ax] = tof_center + tof_half_span * fabs(delta_pos[ax]) / length_LOR;
498 if (tof_index_high[ax]<0)
return 0;
505 for(
int i = 0; i < 3; ++i )
507 if( delta_pos[ i ] != 0.0 )
509 alphaFirst[i] = (-
mp_halfFOV[i] - event1[i]) / (delta_pos[ i ]);
510 alphaLast[i] = (
mp_halfFOV[i] - event1[i]) / (delta_pos[ i ]);
511 alphaMin = (std::max)(alphaMin,(std::min)(alphaFirst[i],alphaLast[i]));
512 alphaMax = (std::min)(alphaMax,(std::max)(alphaFirst[i],alphaLast[i]));
518 if( alphaMax <= alphaMin )
return 0;
522 int iMin = 0, iMax = 0;
523 int jMin = 0, jMax = 0;
524 int kMin = 0, kMax = 0;
527 if( delta_pos[ 0 ] > 0.0f )
530 iMax = ::floor( 1 + ( event1[ 0 ] + alphaMax * delta_pos[ 0 ] - (-
mp_halfFOV[ 0 ]) ) /
mp_sizeVox[0] );
532 else if( delta_pos[ 0 ] < 0.0f )
535 iMax = ::floor( 1 + ( event1[ 0 ] + alphaMin * delta_pos[ 0 ] - (-
mp_halfFOV[ 0 ]) ) /
mp_sizeVox[0] );
537 if( delta_pos[ 0 ] == 0 )
543 if( delta_pos[ 1 ] > 0 )
546 jMax = ::floor( 1 + ( event1[ 1 ] + alphaMax * delta_pos[ 1 ] - (-
mp_halfFOV[ 1 ]) ) /
mp_sizeVox[1] );
548 else if( delta_pos[ 1 ] < 0 )
551 jMax = ::floor( 1 + ( event1[ 1 ] + alphaMin * delta_pos[ 1 ] - (-
mp_halfFOV[ 1 ]) ) /
mp_sizeVox[1] );
553 else if( delta_pos[ 1 ] == 0 )
559 if( delta_pos[ 2 ] > 0 )
562 kMax = ::floor( 1 + ( event1[ 2 ] + alphaMax * delta_pos[ 2 ] - (-
mp_halfFOV[ 2 ]) ) /
mp_sizeVox[2] );
564 else if( delta_pos[ 2 ] < 0 )
567 kMax = ::floor( 1 + ( event1[ 2 ] + alphaMin * delta_pos[ 2 ] - (-
mp_halfFOV[ 2 ]) ) /
mp_sizeVox[2] );
569 else if( delta_pos[ 2 ] == 0 )
575 INTNB n = ( iMax - iMin + 1 ) + ( jMax - jMin + 1 )
576 + ( kMax - kMin + 1 );
579 HPFLTNB alpha_XYZ[ 3 ] = { 1.0f, 1.0f, 1.0f };
582 if( delta_pos[ 0 ] > 0 )
583 alpha_XYZ[ 0 ] = ( ( (-
mp_halfFOV[ 0 ]) + ( iMin - 1 ) *
mp_sizeVox[0] ) - event1[ 0 ] ) / delta_pos[ 0 ];
584 else if( delta_pos[ 0 ] < 0 )
585 alpha_XYZ[ 0 ] = ( ( (-
mp_halfFOV[ 0 ]) + ( iMax - 1 ) *
mp_sizeVox[0] ) - event1[ 0 ] ) / delta_pos[ 0 ];
588 if( delta_pos[ 1 ] > 0 )
589 alpha_XYZ[ 1 ] = ( ( (-
mp_halfFOV[ 1 ]) + ( jMin - 1 ) *
mp_sizeVox[1] ) - event1[ 1 ] ) / delta_pos[ 1 ];
590 else if( delta_pos[ 1 ] < 0 )
591 alpha_XYZ[ 1 ] = ( ( (-
mp_halfFOV[ 1 ]) + ( jMax - 1 ) *
mp_sizeVox[1] ) - event1[ 1 ] ) / delta_pos[ 1 ];
594 if( delta_pos[ 2 ] > 0 )
595 alpha_XYZ[ 2 ] = ( ( (-
mp_halfFOV[ 2 ]) + ( kMin - 1 ) *
mp_sizeVox[2] ) - event1[ 2 ] ) / delta_pos[ 2 ];
596 else if( delta_pos[ 2 ] < 0 )
597 alpha_XYZ[ 2 ] = ( ( (-
mp_halfFOV[ 2 ]) + ( kMax - 1 ) *
mp_sizeVox[2] ) - event1[ 2 ] ) / delta_pos[ 2 ];
607 INTNB index_u[ 3 ] = {
608 (delta_pos[ 0 ] < 0) ? -1 : 1,
609 (delta_pos[ 1 ] < 0) ? -1 : 1,
610 (delta_pos[ 2 ] < 0) ? -1 : 1
614 if( alpha_XYZ[ 0 ] == alphaMin ) alpha_XYZ[ 0 ] += alpha_u[ 0 ];
615 if( alpha_XYZ[ 1 ] == alphaMin ) alpha_XYZ[ 1 ] += alpha_u[ 1 ];
616 if( alpha_XYZ[ 2 ] == alphaMin ) alpha_XYZ[ 2 ] += alpha_u[ 2 ];
619 HPFLTNB const min_alpha_XYZ = (std::min)( alpha_XYZ[ 0 ], (std::min)( alpha_XYZ[ 1 ], alpha_XYZ[ 2 ] ) );
622 HPFLTNB alphaMid = ( min_alpha_XYZ + alphaMin ) / 2.0f;
623 INTNB index_ijk[ 3 ] = {
641 for(
int nP = 0; nP < n - 1; ++nP )
643 if( ( alpha_XYZ[ 0 ] <= alpha_XYZ[ 1 ] )
644 && ( alpha_XYZ[ 0 ] <= alpha_XYZ[ 2 ] ) )
647 if( ( alpha_XYZ[ 0 ] >= alphaMin )
648 && ( index_ijk[ 0 ] - 1 >= tof_index_low[0] )
649 && ( index_ijk[ 0 ] - 1 <= tof_index_high[0] )
650 && ( index_ijk[ 1 ] - 1 >= tof_index_low[1] )
651 && ( index_ijk[ 1 ] - 1 <= tof_index_high[1] )
652 && ( index_ijk[ 2 ] - 1 >= tof_index_low[2] )
653 && ( index_ijk[ 2 ] - 1 <= tof_index_high[2] ) )
662 numVox = ( index_ijk[ 0 ] - 1 ) + ( ( index_ijk[ 1 ] - 1 ) ) * w + ( ( index_ijk[ 2 ] - 1 ) ) * wh;
665 alphaMid = ( alpha_XYZ[ 0 ] + alpha_c ) * 0.5;
668 coeff = ( alpha_XYZ[ 0 ] - alpha_c ) * length_LOR;
686 HPFLTNB temp_erf = std::min(tof_half_span, std::max(-tof_half_span, lor_tof_center -
m_TOFBinSizeInMm/2. - alphaMid * length_LOR));
687 HPFLTNB prev_erf = erf(temp_erf/tof_sigma_sqrt2);
688 temp_erf = std::min(tof_half_span, std::max(-tof_half_span, lor_tof_center +
m_TOFBinSizeInMm/2. - alphaMid * length_LOR));
689 HPFLTNB new_erf = erf(temp_erf/tof_sigma_sqrt2);
704 HPFLTNB temp = (alphaMid * length_LOR - lor_tof_center) / tof_sigma;
714 alpha_c = alpha_XYZ[ 0 ];
715 alpha_XYZ[ 0 ] += alpha_u[ 0 ];
716 index_ijk[ 0 ] += index_u[ 0 ];
718 else if( ( alpha_XYZ[ 1 ] < alpha_XYZ[ 0 ] )
719 && ( alpha_XYZ[ 1 ] <= alpha_XYZ[ 2 ] ) )
722 if( ( alpha_XYZ[ 1 ] >= alphaMin )
723 && ( index_ijk[ 0 ] - 1 >= tof_index_low[0] )
724 && ( index_ijk[ 0 ] - 1 <= tof_index_high[0] )
725 && ( index_ijk[ 1 ] - 1 >= tof_index_low[1] )
726 && ( index_ijk[ 1 ] - 1 <= tof_index_high[1] )
727 && ( index_ijk[ 2 ] - 1 >= tof_index_low[2] )
728 && ( index_ijk[ 2 ] - 1 <= tof_index_high[2] ) )
737 numVox = ( index_ijk[ 0 ] - 1 ) + ( ( index_ijk[ 1 ] - 1 ) ) * w + ( ( index_ijk[ 2 ] - 1 ) ) * wh;
739 alphaMid = ( alpha_XYZ[ 1 ] + alpha_c ) / 2.;
742 coeff = ( alpha_XYZ[ 1 ] - alpha_c ) * length_LOR;
760 HPFLTNB temp_erf = std::min(tof_half_span, std::max(-tof_half_span, lor_tof_center -
m_TOFBinSizeInMm/2. - alphaMid * length_LOR));
761 HPFLTNB prev_erf = erf(temp_erf/tof_sigma_sqrt2);
762 temp_erf = std::min(tof_half_span, std::max(-tof_half_span, lor_tof_center +
m_TOFBinSizeInMm/2. - alphaMid * length_LOR));
763 HPFLTNB new_erf = erf(temp_erf/tof_sigma_sqrt2);
778 HPFLTNB temp = (alphaMid * length_LOR - lor_tof_center) / tof_sigma;
789 alpha_c = alpha_XYZ[ 1 ];
790 alpha_XYZ[ 1 ] += alpha_u[ 1 ];
791 index_ijk[ 1 ] += index_u[ 1 ];
793 else if( ( alpha_XYZ[ 2 ] < alpha_XYZ[ 0 ] )
794 && ( alpha_XYZ[ 2 ] < alpha_XYZ[ 1 ] ) )
797 if( ( alpha_XYZ[ 2 ] >= alphaMin )
798 && ( index_ijk[ 0 ] - 1 >= tof_index_low[0] )
799 && ( index_ijk[ 0 ] - 1 <= tof_index_high[0] )
800 && ( index_ijk[ 1 ] - 1 >= tof_index_low[1] )
801 && ( index_ijk[ 1 ] - 1 <= tof_index_high[1] )
802 && ( index_ijk[ 2 ] - 1 >= tof_index_low[2] )
803 && ( index_ijk[ 2 ] - 1 <= tof_index_high[2] ) )
812 numVox = ( index_ijk[ 0 ] - 1 ) + ( ( index_ijk[ 1 ] - 1 ) ) * w + ( ( index_ijk[ 2 ] - 1 ) ) * wh;
814 alphaMid = ( alpha_XYZ[ 2 ] + alpha_c ) / 2.;
817 coeff = ( alpha_XYZ[ 2 ] - alpha_c ) * length_LOR;
835 HPFLTNB temp_erf = std::min(tof_half_span, std::max(-tof_half_span, lor_tof_center -
m_TOFBinSizeInMm/2. - alphaMid * length_LOR));
836 HPFLTNB prev_erf = erf(temp_erf/tof_sigma_sqrt2);
837 temp_erf = std::min(tof_half_span, std::max(-tof_half_span, lor_tof_center +
m_TOFBinSizeInMm/2. - alphaMid * length_LOR));
838 HPFLTNB new_erf = erf(temp_erf/tof_sigma_sqrt2);
853 HPFLTNB temp = (alphaMid * length_LOR - lor_tof_center) / tof_sigma;
864 alpha_c = alpha_XYZ[ 2 ];
865 alpha_XYZ[ 2 ] += alpha_u[ 2 ];
866 index_ijk[ 2 ] += index_u[ 2 ];
883 Cerr(
"***** iProjectorIncrementalSiddon::ProjectTOFHistogram() -> Called while not initialized !" << endl);
888 #ifdef CASTOR_VERBOSE 891 string direction =
"";
892 if (a_direction==
FORWARD) direction =
"forward";
893 else direction =
"backward";
894 Cout(
"iProjectorIncrementalSiddon::Project with TOF bins -> Project line '" << ap_ProjectionLine <<
"' in " << direction <<
" direction" << endl);
901 HPFLTNB event1[3] = { event1Float[0], event1Float[1], event1Float[2] };
902 HPFLTNB event2[3] = { event2Float[0], event2Float[1], event2Float[2] };
908 HPFLTNB length_LOR_half = length_LOR * 0.5;
913 HPFLTNB alphaFirst[] = { 0.0, 0.0, 0.0 };
914 HPFLTNB alphaLast[] = { 0.0, 0.0, 0.0 };
916 HPFLTNB alphaMin = 0.0f, alphaMax = 1.0f;
918 event2[ 0 ] - event1[ 0 ],
919 event2[ 1 ] - event1[ 1 ],
920 event2[ 2 ] - event1[ 2 ]
925 INTNB tof_half_nb_bins = tof_nb_bins/2;
929 HPFLTNB tof_sigma_sqrt2 = sqrt(2.)*tof_sigma;
935 HPFLTNB prev_erf = 0., new_erf = 0.;
938 INTNB tof_bin_last = tof_half_nb_bins;
939 INTNB tof_bin_first = -tof_half_nb_bins;
947 INTNB tof_bin_first_for_voxel = 0, tof_bin_last_for_voxel = 0;
950 for(
int i = 0; i < 3; ++i )
952 if( delta_pos[ i ] != 0.0 )
954 alphaFirst[i] = (-
mp_halfFOV[i] - event1[i]) / (delta_pos[ i ]);
955 alphaLast[i] = (
mp_halfFOV[i] - event1[i]) / (delta_pos[ i ]);
956 alphaMin = (std::max)(alphaMin,(std::min)(alphaFirst[i],alphaLast[i]));
957 alphaMax = (std::min)(alphaMax,(std::max)(alphaFirst[i],alphaLast[i]));
963 if( alphaMax <= alphaMin )
return 0;
968 for (
INTNB tof_bin=0; tof_bin<tof_nb_bins; tof_bin++) tof_weights_temp[tof_bin] = 0.;
972 int iMin = 0, iMax = 0;
973 int jMin = 0, jMax = 0;
974 int kMin = 0, kMax = 0;
977 if( delta_pos[ 0 ] > 0.0f )
980 iMax = ::floor( 1 + ( event1[ 0 ] + alphaMax * delta_pos[ 0 ] - (-
mp_halfFOV[ 0 ]) ) /
mp_sizeVox[0] );
982 else if( delta_pos[ 0 ] < 0.0f )
985 iMax = ::floor( 1 + ( event1[ 0 ] + alphaMin * delta_pos[ 0 ] - (-
mp_halfFOV[ 0 ]) ) /
mp_sizeVox[0] );
987 if( delta_pos[ 0 ] == 0 )
993 if( delta_pos[ 1 ] > 0 )
996 jMax = ::floor( 1 + ( event1[ 1 ] + alphaMax * delta_pos[ 1 ] - (-
mp_halfFOV[ 1 ]) ) /
mp_sizeVox[1] );
998 else if( delta_pos[ 1 ] < 0 )
1001 jMax = ::floor( 1 + ( event1[ 1 ] + alphaMin * delta_pos[ 1 ] - (-
mp_halfFOV[ 1 ]) ) /
mp_sizeVox[1] );
1003 else if( delta_pos[ 1 ] == 0 )
1009 if( delta_pos[ 2 ] > 0 )
1012 kMax = ::floor( 1 + ( event1[ 2 ] + alphaMax * delta_pos[ 2 ] - (-
mp_halfFOV[ 2 ]) ) /
mp_sizeVox[2] );
1014 else if( delta_pos[ 2 ] < 0 )
1017 kMax = ::floor( 1 + ( event1[ 2 ] + alphaMin * delta_pos[ 2 ] - (-
mp_halfFOV[ 2 ]) ) /
mp_sizeVox[2] );
1019 else if( delta_pos[ 2 ] == 0 )
1025 INTNB n = ( iMax - iMin + 1 ) + ( jMax - jMin + 1 )
1026 + ( kMax - kMin + 1 );
1029 HPFLTNB alpha_XYZ[ 3 ] = { 1.0f, 1.0f, 1.0f };
1032 if( delta_pos[ 0 ] > 0 )
1033 alpha_XYZ[ 0 ] = ( ( (-
mp_halfFOV[ 0 ]) + ( iMin - 1 ) *
mp_sizeVox[0] ) - event1[ 0 ] ) / delta_pos[ 0 ];
1034 else if( delta_pos[ 0 ] < 0 )
1035 alpha_XYZ[ 0 ] = ( ( (-
mp_halfFOV[ 0 ]) + ( iMax - 1 ) *
mp_sizeVox[0] ) - event1[ 0 ] ) / delta_pos[ 0 ];
1038 if( delta_pos[ 1 ] > 0 )
1039 alpha_XYZ[ 1 ] = ( ( (-
mp_halfFOV[ 1 ]) + ( jMin - 1 ) *
mp_sizeVox[1] ) - event1[ 1 ] ) / delta_pos[ 1 ];
1040 else if( delta_pos[ 1 ] < 0 )
1041 alpha_XYZ[ 1 ] = ( ( (-
mp_halfFOV[ 1 ]) + ( jMax - 1 ) *
mp_sizeVox[1] ) - event1[ 1 ] ) / delta_pos[ 1 ];
1044 if( delta_pos[ 2 ] > 0 )
1045 alpha_XYZ[ 2 ] = ( ( (-
mp_halfFOV[ 2 ]) + ( kMin - 1 ) *
mp_sizeVox[2] ) - event1[ 2 ] ) / delta_pos[ 2 ];
1046 else if( delta_pos[ 2 ] < 0 )
1047 alpha_XYZ[ 2 ] = ( ( (-
mp_halfFOV[ 2 ]) + ( kMax - 1 ) *
mp_sizeVox[2] ) - event1[ 2 ] ) / delta_pos[ 2 ];
1057 INTNB index_u[ 3 ] = {
1058 (delta_pos[ 0 ] < 0) ? -1 : 1,
1059 (delta_pos[ 1 ] < 0) ? -1 : 1,
1060 (delta_pos[ 2 ] < 0) ? -1 : 1
1064 if( alpha_XYZ[ 0 ] == alphaMin ) alpha_XYZ[ 0 ] += alpha_u[ 0 ];
1065 if( alpha_XYZ[ 1 ] == alphaMin ) alpha_XYZ[ 1 ] += alpha_u[ 1 ];
1066 if( alpha_XYZ[ 2 ] == alphaMin ) alpha_XYZ[ 2 ] += alpha_u[ 2 ];
1069 HPFLTNB const min_alpha_XYZ = (std::min)( alpha_XYZ[ 0 ], (std::min)( alpha_XYZ[ 1 ], alpha_XYZ[ 2 ] ) );
1072 HPFLTNB alphaMid = ( min_alpha_XYZ + alphaMin ) / 2.0f;
1073 INTNB index_ijk[ 3 ] = {
1087 for(
int nP = 0; nP < n - 1; ++nP )
1089 if( ( alpha_XYZ[ 0 ] <= alpha_XYZ[ 1 ] )
1090 && ( alpha_XYZ[ 0 ] <= alpha_XYZ[ 2 ] ) )
1093 if( ( alpha_XYZ[ 0 ] >= alphaMin )
1094 && ( index_ijk[ 0 ] - 1 >= 0 )
1095 && ( index_ijk[ 0 ] - 1 <= mp_nbVox[ 0 ] - 1 )
1096 && ( index_ijk[ 1 ] - 1 >= 0 )
1097 && ( index_ijk[ 1 ] - 1 <= mp_nbVox[ 1 ] - 1 )
1098 && ( index_ijk[ 2 ] - 1 >= 0 )
1099 && ( index_ijk[ 2 ] - 1 <= mp_nbVox[ 2 ] - 1 ) )
1101 numVox = ( index_ijk[ 0 ] - 1 ) + ( ( index_ijk[ 1 ] - 1 ) ) * w + ( ( index_ijk[ 2 ] - 1 ) ) * wh;
1106 coeff = ( alpha_XYZ[ 0 ] - alpha_c ) * length_LOR;
1108 alphaMid = ( alpha_XYZ[ 0 ] + alpha_c ) * 0.5;
1120 tof_bin_first_for_voxel = max(tof_bin_first , (
INTNB)(ceil( (alphaMid * length_LOR - tof_half_span - length_LOR_half ) /
m_TOFBinSizeInMm )));
1121 tof_bin_last_for_voxel = min(tof_bin_last , (
INTNB)(floor( (alphaMid * length_LOR + tof_half_span - length_LOR_half) /
m_TOFBinSizeInMm )));
1124 lor_tof_center = length_LOR_half + tof_bin_first_for_voxel *
m_TOFBinSizeInMm;
1127 tof_bin_first_for_voxel += tof_half_nb_bins;
1128 tof_bin_last_for_voxel += tof_half_nb_bins;
1134 HPFLTNB temp = std::min(tof_half_span,std::max(-tof_half_span, lor_tof_center - m_TOFBinSizeInMm/2. - alphaMid * length_LOR));
1135 prev_erf = erf(temp/tof_sigma_sqrt2);
1140 for (
INTNB tof_bin=tof_bin_first_for_voxel; tof_bin<=tof_bin_last_for_voxel; tof_bin++)
1160 HPFLTNB temp = std::min(tof_half_span, std::max(-tof_half_span, lor_tof_center + m_TOFBinSizeInMm/2. - alphaMid * length_LOR));
1161 new_erf = erf(temp/tof_sigma_sqrt2);
1162 tof_weights_temp[tof_bin] = 0.5 * fabs(new_erf - prev_erf);
1172 HPFLTNB temp = ( alphaMid * length_LOR - lor_tof_center) / tof_sigma;
1190 ap_ProjectionLine->
AddVoxelAllTOFBins(a_direction, numVox, coeff, tof_weights_temp, tof_bin_first_for_voxel, tof_bin_last_for_voxel);
1195 alpha_c = alpha_XYZ[ 0 ];
1196 alpha_XYZ[ 0 ] += alpha_u[ 0 ];
1197 index_ijk[ 0 ] += index_u[ 0 ];
1199 else if( ( alpha_XYZ[ 1 ] < alpha_XYZ[ 0 ] )
1200 && ( alpha_XYZ[ 1 ] <= alpha_XYZ[ 2 ] ) )
1203 if( ( alpha_XYZ[ 1 ] >= alphaMin )
1204 && ( index_ijk[ 0 ] - 1 >= 0 )
1205 && ( index_ijk[ 0 ] - 1 <= mp_nbVox[ 0 ] - 1 )
1206 && ( index_ijk[ 1 ] - 1 >= 0 )
1207 && ( index_ijk[ 1 ] - 1 <= mp_nbVox[ 1 ] - 1 )
1208 && ( index_ijk[ 2 ] - 1 >= 0 )
1209 && ( index_ijk[ 2 ] - 1 <= mp_nbVox[ 2 ] - 1 ) )
1211 numVox = ( index_ijk[ 0 ] - 1 ) + ( ( index_ijk[ 1 ] - 1 ) ) * w + ( ( index_ijk[ 2 ] - 1 ) ) * wh;
1216 coeff = ( alpha_XYZ[ 1 ] - alpha_c ) * length_LOR;
1218 alphaMid = ( alpha_XYZ[ 1 ] + alpha_c ) * 0.5;
1230 tof_bin_first_for_voxel = max(tof_bin_first , (
INTNB)(ceil( (alphaMid * length_LOR - tof_half_span - length_LOR_half ) /
m_TOFBinSizeInMm )));
1231 tof_bin_last_for_voxel = min(tof_bin_last , (
INTNB)(floor( (alphaMid * length_LOR + tof_half_span - length_LOR_half) /
m_TOFBinSizeInMm )));
1234 lor_tof_center = length_LOR_half + tof_bin_first_for_voxel *
m_TOFBinSizeInMm;
1237 tof_bin_first_for_voxel += tof_half_nb_bins;
1238 tof_bin_last_for_voxel += tof_half_nb_bins;
1244 HPFLTNB temp = std::min(tof_half_span,std::max(-tof_half_span, lor_tof_center - m_TOFBinSizeInMm/2. - alphaMid * length_LOR));
1245 prev_erf = erf(temp/tof_sigma_sqrt2);
1250 for (
INTNB tof_bin=tof_bin_first_for_voxel; tof_bin<=tof_bin_last_for_voxel; tof_bin++)
1270 HPFLTNB temp = std::min(tof_half_span, std::max(-tof_half_span, lor_tof_center + m_TOFBinSizeInMm/2. - alphaMid * length_LOR));
1271 new_erf = erf(temp/tof_sigma_sqrt2);
1272 tof_weights_temp[tof_bin] = 0.5 * fabs(new_erf - prev_erf);
1282 HPFLTNB temp = ( alphaMid * length_LOR - lor_tof_center) / tof_sigma;
1300 ap_ProjectionLine->
AddVoxelAllTOFBins(a_direction, numVox, coeff, tof_weights_temp, tof_bin_first_for_voxel, tof_bin_last_for_voxel);
1305 alpha_c = alpha_XYZ[ 1 ];
1306 alpha_XYZ[ 1 ] += alpha_u[ 1 ];
1307 index_ijk[ 1 ] += index_u[ 1 ];
1309 else if( ( alpha_XYZ[ 2 ] < alpha_XYZ[ 0 ] )
1310 && ( alpha_XYZ[ 2 ] < alpha_XYZ[ 1 ] ) )
1313 if( ( alpha_XYZ[ 2 ] >= alphaMin )
1314 && ( index_ijk[ 0 ] - 1 >= 0 )
1315 && ( index_ijk[ 0 ] - 1 <= mp_nbVox[ 0 ] - 1 )
1316 && ( index_ijk[ 1 ] - 1 >= 0 )
1317 && ( index_ijk[ 1 ] - 1 <= mp_nbVox[ 1 ] - 1 )
1318 && ( index_ijk[ 2 ] - 1 >= 0 )
1319 && ( index_ijk[ 2 ] - 1 <= mp_nbVox[ 2 ] - 1 ) )
1321 numVox = ( index_ijk[ 0 ] - 1 ) + ( ( index_ijk[ 1 ] - 1 ) ) * w + ( ( index_ijk[ 2 ] - 1 ) ) * wh;
1326 coeff = ( alpha_XYZ[ 2 ] - alpha_c ) * length_LOR;
1328 alphaMid = ( alpha_XYZ[ 2 ] + alpha_c ) * 0.5;
1340 tof_bin_first_for_voxel = max(tof_bin_first , (
INTNB)(ceil( (alphaMid * length_LOR - tof_half_span - length_LOR_half ) /
m_TOFBinSizeInMm )));
1341 tof_bin_last_for_voxel = min(tof_bin_last , (
INTNB)(floor( (alphaMid * length_LOR + tof_half_span - length_LOR_half) /
m_TOFBinSizeInMm )));
1344 lor_tof_center = length_LOR_half + tof_bin_first_for_voxel *
m_TOFBinSizeInMm;
1347 tof_bin_first_for_voxel += tof_half_nb_bins;
1348 tof_bin_last_for_voxel += tof_half_nb_bins;
1355 HPFLTNB temp = std::min(tof_half_span,std::max(-tof_half_span, lor_tof_center - m_TOFBinSizeInMm/2. - alphaMid * length_LOR));
1356 prev_erf = erf(temp/tof_sigma_sqrt2);
1361 for (
INTNB tof_bin=tof_bin_first_for_voxel; tof_bin<=tof_bin_last_for_voxel; tof_bin++)
1381 HPFLTNB temp = std::min(tof_half_span, std::max(-tof_half_span, lor_tof_center + m_TOFBinSizeInMm/2. - alphaMid * length_LOR));
1382 new_erf = erf(temp/tof_sigma_sqrt2);
1383 tof_weights_temp[tof_bin] = 0.5 * fabs(new_erf - prev_erf);
1393 HPFLTNB temp = ( alphaMid * length_LOR - lor_tof_center) / tof_sigma;
1410 ap_ProjectionLine->
AddVoxelAllTOFBins(a_direction, numVox, coeff, tof_weights_temp, tof_bin_first_for_voxel, tof_bin_last_for_voxel);
1415 alpha_c = alpha_XYZ[ 2 ];
1416 alpha_XYZ[ 2 ] += alpha_u[ 2 ];
1417 index_ijk[ 2 ] += index_u[ 2 ];
1421 delete[] tof_weights_temp;
bool m_compatibleWithSPECTAttenuationCorrection
int ProjectTOFListmode(int a_direction, oProjectionLine *ap_ProjectionLine)
A function to project with TOF continuous information.
~iProjectorIncrementalSiddon()
The destructor of iProjectorIncrementalSiddon.
#define SPEED_OF_LIGHT_IN_MM_PER_PS
FLTNB m_TOFGaussianNormCoef
void AddVoxelAllTOFBins(int a_direction, INTNB a_voxelIndex, FLTNB a_voxelWeight, HPFLTNB *a_tofWeights, INTNB a_tofBinFirst, INTNB a_tofBinLast)
Add a voxel contribution to the line for all relevant TOF bins.
#define TWO_SQRT_TWO_LN_2
oImageDimensionsAndQuantification * mp_ImageDimensionsAndQuantification
void ShowHelpSpecific()
A function used to show help about the child module.
FLTNB GetTOFMeasurementInPs()
This function is used to get the TOF measurement in ps.
This class is designed to generically described any on-the-fly projector.
bool m_TOFWeightingFcnPrecomputedFlag
bool m_TOFBinProperProcessingFlag
HPFLTNB * mp_TOFWeightingFcn
Declaration of class iProjectorIncrementalSiddon.
FLTNB GetLength()
This function is used to get the length of the line.
FLTNB * GetPosition1()
This function is used to get the pointer to the mp_position1 (3-values tab).
void AddVoxel(int a_direction, INTNB a_voxelIndice, FLTNB a_voxelWeight)
This function is used to add a voxel contribution to the line, assuming TOF bin 0 (i...
int ReadConfigurationFile(const string &a_configurationFile)
A function used to read options from a configuration file.
FLTNB m_TOFResolutionInMm
iProjectorIncrementalSiddon()
The constructor of iProjectorIncrementalSiddon.
INTNB EstimateMaxNumberOfVoxelsPerLine()
This function is used to compute and provide an estimate of the maximum number of voxels that could c...
int InitializeSpecific()
This function is used to initialize specific stuff to the child projector.
This class is designed to manage and store system matrix elements associated to a vEvent...
FLTNB m_TOFPrecomputedSamplingFactor
int ProjectWithoutTOF(int a_direction, oProjectionLine *ap_ProjectionLine)
A function to project without TOF.
int ProjectTOFHistogram(int a_direction, oProjectionLine *ap_ProjectionLine)
A function to project with TOF binned information.
Declaration of class sOutputManager.
int GetNbTOFBins()
This function is used to get the number of TOF bins in use.
FLTNB * GetPosition2()
This function is used to get the pointer to the mp_position2 (3-values tab).
int ReadOptionsList(const string &a_optionsList)
A function used to read options from a list of options.
INTNB GetNbVoxX()
Get the number of voxels along the X axis.
INTNB GetNbVoxZ()
Get the number of voxels along the Z axis.
bool m_compatibleWithCompression
int CheckSpecificParameters()
A private function used to check the parameters settings specific to the child projector.
INTNB GetNbVoxY()
Get the number of voxels along the Y axis.
INTNB m_TOFWeightingFcnNbSamples