130 cout <<
"This projector is a line projector that uses linear interpolation between pixels." << endl;
131 cout <<
"It is implemented from the following published paper:" << endl;
132 cout <<
"P. M. Joseph, \"An improved algorithm for reprojecting rays through pixel images\", IEEE Trans. Med. Imaging, vol. 1, pp. 192-6, 1982." << endl;
156 if (
m_verbose>=2)
Cout(
"iProjectorJoseph::InitializeSpecific() -> Use Joseph projector" << endl);
167 INTNB nElts = ( mp_nbVox[ 0 ] + 2 ) * ( mp_nbVox[ 1 ] + 2 ) * ( mp_nbVox[ 2 ] + 2 );
169 ::memset(
mp_maskPad, 0,
sizeof( uint8_t ) * nElts );
170 for(
INTNB k = 1; k < mp_nbVox[ 2 ] + 1; ++k )
172 for(
INTNB j = 1; j < mp_nbVox[ 1 ] + 1; ++j )
174 for(
INTNB i = 1; i < mp_nbVox[ 0 ] + 1; ++i )
176 mp_maskPad[ i + j * ( ( mp_nbVox[ 0 ] + 2 ) ) + k * ( mp_nbVox[ 0 ] + 2 ) * ( mp_nbVox[ 1 ] + 2 ) ] = 1;
182 HPFLTNB tolerance_factor = 1.e-4;
208 max_nb_voxels_in_dimension *= 4;
210 return max_nb_voxels_in_dimension;
223 Cerr(
"***** iProjectorJoseph::ProjectWithoutTOF() -> Called while not initialized !" << endl);
228 #ifdef CASTOR_VERBOSE 231 string direction =
"";
232 if (a_direction==
FORWARD) direction =
"forward";
233 else direction =
"backward";
234 Cout(
"iProjectorJoseph::Project without TOF -> Project line '" << ap_ProjectionLine <<
"' in " << direction <<
" direction" << endl);
244 event2[ 0 ] - event1[ 0 ],
245 event2[ 1 ] - event1[ 1 ],
246 event2[ 2 ] - event1[ 2 ]
258 HPFLTNB alphaMin = 0.0, alphaMax = 1.0;
261 HPFLTNB pos[ 3 ] = { 0.0, 0.0, 0.0 };
262 HPFLTNB wfl[ 2 ] = { 0.0, 0.0 };
263 HPFLTNB wcl[ 2 ] = { 0.0, 0.0 };
264 HPFLTNB w[ 4 ] = { 0.0, 0.0, 0.0, 0.0 };
265 int16_t index[ 2 ] = { 0, 0 };
266 int32_t finalIdx = 0;
267 int8_t limitX1 = 1; int8_t limitX2 = 1;
268 int8_t limitY1 = 1; int8_t limitY2 = 1;
269 int8_t limitZ1 = 1; int8_t limitZ2 = 1;
272 if( ::fabs( r[ 0 ] ) > ::fabs( r[ 1 ] ) )
279 alphaMin = std::max( alphaMin, std::min( alphaX_0, alphaX_1 ) );
280 alphaMax = std::min( alphaMax, std::max( alphaX_0, alphaX_1 ) );
291 alphaMin = std::max( alphaMin, std::min( alphaY_0, alphaY_1 ) );
292 alphaMax = std::min( alphaMax, std::max( alphaY_0, alphaY_1 ) );
304 alphaMin = std::max( alphaMin, std::min( alphaZ_0, alphaZ_1 ) );
305 alphaMax = std::min( alphaMax, std::max( alphaZ_0, alphaZ_1 ) );
308 if( alphaMax <= alphaMin )
return 0;
310 if( r[ 1 ] == 0.0 )
if( event1[ 1 ] < -
mp_halfFOV[ 1 ] || event1[ 1 ] >
mp_halfFOV[ 1 ] )
return 0;
311 if( r[ 2 ] == 0.0 )
if( event1[ 2 ] < -
mp_halfFOV[ 2 ] || event1[ 2 ] >
mp_halfFOV[ 2 ] )
return 0;
314 HPFLTNB const factor( ::fabs( r[ 0 ] )
315 / ::sqrt( r2[ 0 ] + r2[ 1 ] + r2[ 2 ] ) );
326 int iMin = 0, iMax = 0;
332 else if( r[ 0 ] < 0.0 )
342 for(
int i = iMin; i < iMax; ++i )
346 pos[ 1 ] = event1[ 1 ] + step * ri[ 1 ];
347 pos[ 2 ] = event1[ 2 ] + step * ri[ 2 ];
354 finalIdx = i + ( index[ 0 ] - 1 ) *
mp_nbVox[ 0 ] + ( index[ 1 ] - 1 ) *
m_nbVoxXY;
357 wfl[ 0 ] = pos[ 1 ] - (
m_boundY + index[ 0 ] * mp_sizeVox[ 1 ] );
358 wfl[ 1 ] = pos[ 2 ] - (
m_boundZ + index[ 1 ] * mp_sizeVox[ 2 ] );
359 wfl[ 0 ] /= mp_sizeVox[ 1 ];
360 wfl[ 1 ] /= mp_sizeVox[ 2 ];
363 wcl[ 0 ] = 1.0 - wfl[ 0 ];
364 wcl[ 1 ] = 1.0 - wfl[ 1 ];
367 w[ 0 ] = wcl[ 0 ] * wcl[ 1 ];
368 w[ 1 ] = wcl[ 0 ] * wfl[ 1 ];
369 w[ 2 ] = wfl[ 0 ] * wfl[ 1 ];
370 w[ 3 ] = wfl[ 0 ] * wcl[ 1 ];
373 if( index[ 0 ] <= 0 ) limitY1 = 0;
374 if( index[ 0 ] >=
mp_nbVox[ 1 ] ) limitY2 = 0;
377 if( index[ 1 ] <= 0 ) limitZ1 = 0;
378 if( index[ 1 ] >=
mp_nbVox[ 2 ] ) limitZ2 = 0;
382 ap_ProjectionLine->
AddVoxel(a_direction, finalIdx * limitY1 * limitZ1, w[ 0 ] * weight * limitY1 * limitZ1);
383 ap_ProjectionLine->
AddVoxel(a_direction, ( finalIdx +
m_nbVoxXY ) * limitY1 * limitZ2, w[ 1 ] * weight * limitY1 * limitZ2);
384 ap_ProjectionLine->
AddVoxel(a_direction, ( finalIdx +
mp_nbVox[ 0 ] +
m_nbVoxXY ) * limitY2 * limitZ2, w[ 2 ] * weight * limitY2 * limitZ2);
385 ap_ProjectionLine->
AddVoxel(a_direction, ( finalIdx +
mp_nbVox[ 0 ] ) * limitY2 * limitZ1, w[ 3 ] * weight * limitY2 * limitZ1);
387 limitY1 = 1; limitY2 = 1;
388 limitZ1 = 1; limitZ2 = 1;
402 alphaMin = std::max( alphaMin, std::min( alphaX_0, alphaX_1 ) );
403 alphaMax = std::min( alphaMax, std::max( alphaX_0, alphaX_1 ) );
411 alphaMin = std::max( alphaMin, std::min( alphaY_0, alphaY_1 ) );
412 alphaMax = std::min( alphaMax, std::max( alphaY_0, alphaY_1 ) );
423 alphaMin = std::max( alphaMin, std::min( alphaZ_0, alphaZ_1 ) );
424 alphaMax = std::min( alphaMax, std::max( alphaZ_0, alphaZ_1 ) );
427 if( alphaMax <= alphaMin )
return 0;
429 if( r[ 0 ] == 0.0 )
if( event1[ 0 ] < -
mp_halfFOV[ 0 ] || event1[ 0 ] >
mp_halfFOV[ 0 ] )
return 0;
430 if( r[ 2 ] == 0.0 )
if( event1[ 2 ] < -
mp_halfFOV[ 2 ] || event1[ 2 ] >
mp_halfFOV[ 2 ] )
return 0;
433 HPFLTNB const factor( ::fabs( r[ 1 ] ) / ::sqrt( r2[ 0 ] + r2[ 1 ] + r2[ 2 ] ) );
444 int jMin = 0, jMax = 0;
450 else if( r[ 1 ] < 0.0 )
460 for(
int j = jMin; j < jMax; ++j )
464 pos[ 0 ] = event1[ 0 ] + step * ri[ 0 ];
465 pos[ 2 ] = event1[ 2 ] + step * ri[ 2 ];
472 finalIdx = ( index[ 0 ] - 1 ) + j *
mp_nbVox[ 0 ] + ( index[ 1 ] - 1 ) *
m_nbVoxXY;
475 wfl[ 0 ] = pos[ 0 ] - (
m_boundX + index[ 0 ] * mp_sizeVox[ 0 ] );
476 wfl[ 1 ] = pos[ 2 ] - (
m_boundZ + index[ 1 ] * mp_sizeVox[ 2 ] );
477 wfl[ 0 ] /= mp_sizeVox[ 0 ];
478 wfl[ 1 ] /= mp_sizeVox[ 2 ];
481 wcl[ 0 ] = 1.0 - wfl[ 0 ];
482 wcl[ 1 ] = 1.0 - wfl[ 1 ];
485 w[ 0 ] = wcl[ 0 ] * wcl[ 1 ];
486 w[ 1 ] = wcl[ 0 ] * wfl[ 1 ];
487 w[ 2 ] = wfl[ 0 ] * wfl[ 1 ];
488 w[ 3 ] = wfl[ 0 ] * wcl[ 1 ];
491 if( index[ 0 ] <= 0 ) limitX1 = 0;
492 if( index[ 0 ] >=
mp_nbVox[ 0 ] ) limitX2 = 0;
495 if( index[ 1 ] <= 0 ) limitZ1 = 0;
496 if( index[ 1 ] >=
mp_nbVox[ 2 ] ) limitZ2 = 0;
500 ap_ProjectionLine->
AddVoxel(a_direction, finalIdx * limitX1 * limitZ1, w[ 0 ] * weight * limitX1 * limitZ1);
501 ap_ProjectionLine->
AddVoxel(a_direction, ( finalIdx +
m_nbVoxXY ) * limitX1 * limitZ2, w[ 1 ] * weight * limitX1 * limitZ2);
502 ap_ProjectionLine->
AddVoxel(a_direction, ( finalIdx + 1 +
m_nbVoxXY ) * limitX2 * limitZ2, w[ 2 ] * weight * limitX2 * limitZ2);
503 ap_ProjectionLine->
AddVoxel(a_direction, ( finalIdx + 1 ) * limitX2 * limitZ1, w[ 3 ] * weight * limitX2 * limitZ1);
506 limitX1 = 1; limitX2 = 1;
507 limitZ1 = 1; limitZ2 = 1;
525 Cerr(
"***** iProjectorJoseph::ProjectTOFListmode() -> Called while not initialized !" << endl);
530 #ifdef CASTOR_VERBOSE 533 string direction =
"";
534 if (a_direction==
FORWARD) direction =
"forward";
535 else direction =
"backward";
536 Cout(
"iProjectorJoseph::Project with TOF measurement -> Project line '" << ap_ProjectionLine <<
"' in " << direction <<
" direction" << endl);
546 event2[ 0 ] - event1[ 0 ],
547 event2[ 1 ] - event1[ 1 ],
548 event2[ 2 ] - event1[ 2 ]
571 HPFLTNB tof_sigma_sqrt2 = sqrt(2.)*tof_sigma;
579 HPFLTNB lor_tof_center = lor_length * 0.5 + tof_delta;
582 HPFLTNB tof_edge_low[] = {0,0,0};
584 HPFLTNB tof_edge_high[] = {0,0,0};
589 for (
int ax=0;ax<3;ax++)
592 tof_center = event1[ax] + lor_tof_center * r[ax] / lor_length;
595 tof_edge_low[ax] = tof_center - tof_half_span * fabs(r[ax]) / lor_length;
598 if (tof_index>
mp_nbVox[ax]-1)
return 0;
602 tof_edge_high[ax] = tof_center + tof_half_span * fabs(r[ax]) / lor_length;
605 if (tof_index<0)
return 0;
611 HPFLTNB alphaMin = 0.0, alphaMax = 1.0;
614 HPFLTNB pos[ 3 ] = { 0.0, 0.0, 0.0 };
615 HPFLTNB wfl[ 2 ] = { 0.0, 0.0 };
616 HPFLTNB wcl[ 2 ] = { 0.0, 0.0 };
617 HPFLTNB w[ 4 ] = { 0.0, 0.0, 0.0, 0.0 };
618 int16_t index[ 2 ] = { 0, 0 };
619 int32_t finalIdx = 0;
620 int8_t limitX1 = 1; int8_t limitX2 = 1;
621 int8_t limitY1 = 1; int8_t limitY2 = 1;
622 int8_t limitZ1 = 1; int8_t limitZ2 = 1;
625 if( ::fabs( r[ 0 ] ) > ::fabs( r[ 1 ] ) )
630 HPFLTNB const alphaX_0 = ( tof_edge_low[ 0 ] - event1[ 0 ] ) / r[ 0 ];
631 HPFLTNB const alphaX_1 = ( tof_edge_high[ 0 ] - event1[ 0 ] ) / r[ 0 ];
632 alphaMin = std::max( alphaMin, std::min( alphaX_0, alphaX_1 ) );
633 alphaMax = std::min( alphaMax, std::max( alphaX_0, alphaX_1 ) );
640 HPFLTNB const alphaY_0 = ( tof_edge_low[ 1 ] -
mp_sizeVox[ 1 ] * 0.5 - event1[ 1 ] )
642 HPFLTNB const alphaY_1 = ( tof_edge_high[ 1 ] +
mp_sizeVox[ 1 ] * 0.5 - event1[ 1 ] )
644 alphaMin = std::max( alphaMin, std::min( alphaY_0, alphaY_1 ) );
645 alphaMax = std::min( alphaMax, std::max( alphaY_0, alphaY_1 ) );
653 HPFLTNB const alphaZ_0 = ( tof_edge_low[ 2 ] -
mp_sizeVox[ 2 ] * 0.5 - event1[ 2 ] )
655 HPFLTNB const alphaZ_1 = ( tof_edge_high[ 2 ] +
mp_sizeVox[ 2 ] * 0.5 - event1[ 2 ] )
657 alphaMin = std::max( alphaMin, std::min( alphaZ_0, alphaZ_1 ) );
658 alphaMax = std::min( alphaMax, std::max( alphaZ_0, alphaZ_1 ) );
661 if( alphaMax <= alphaMin )
return 0;
663 if( r[ 1 ] == 0.0 )
if( event1[ 1 ] < -
mp_halfFOV[ 1 ] || event1[ 1 ] >
mp_halfFOV[ 1 ] )
return 0;
664 if( r[ 2 ] == 0.0 )
if( event1[ 2 ] < -
mp_halfFOV[ 2 ] || event1[ 2 ] >
mp_halfFOV[ 2 ] )
return 0;
667 HPFLTNB const factor( ::fabs( r[ 0 ] ) / lor_length );
668 HPFLTNB const factor_for_tof( lor_length / r[ 0 ]) ;
679 int iMin = 0, iMax = 0;
685 else if( r[ 0 ] < 0.0 )
700 for(
int i = iMin; i < iMax; ++i )
704 pos[ 1 ] = event1[ 1 ] + step * ri[ 1 ];
705 pos[ 2 ] = event1[ 2 ] + step * ri[ 2 ];
712 wfl[ 0 ] = pos[ 1 ] - (
m_boundY + index[ 0 ] * mp_sizeVox[ 1 ] );
713 wfl[ 1 ] = pos[ 2 ] - (
m_boundZ + index[ 1 ] * mp_sizeVox[ 2 ] );
714 wfl[ 0 ] /= mp_sizeVox[ 1 ];
715 wfl[ 1 ] /= mp_sizeVox[ 2 ];
718 wcl[ 0 ] = 1.0 - wfl[ 0 ];
719 wcl[ 1 ] = 1.0 - wfl[ 1 ];
722 w[ 0 ] = wcl[ 0 ] * wcl[ 1 ];
723 w[ 1 ] = wcl[ 0 ] * wfl[ 1 ];
724 w[ 2 ] = wfl[ 0 ] * wfl[ 1 ];
725 w[ 3 ] = wfl[ 0 ] * wcl[ 1 ];
728 if( index[ 0 ] <= 0 ) limitY1 = 0;
729 if( index[ 0 ] >=
mp_nbVox[ 1 ] ) limitY2 = 0;
732 if( index[ 1 ] <= 0 ) limitZ1 = 0;
733 if( index[ 1 ] >=
mp_nbVox[ 2 ] ) limitZ2 = 0;
736 finalIdx = i + ( index[ 0 ] - 1 ) *
mp_nbVox[ 0 ] + ( index[ 1 ] - 1 ) *
m_nbVoxXY;
760 HPFLTNB temp_erf = std::min(tof_half_span, std::max(-tof_half_span, lor_tof_center -
m_TOFBinSizeInMm/2. - step * factor_for_tof));
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. - step * factor_for_tof));
763 HPFLTNB new_erf = erf(temp_erf/tof_sigma_sqrt2);
770 HPFLTNB temp = ( step * factor_for_tof - lor_tof_center ) / tof_sigma;
778 ap_ProjectionLine->
AddVoxel(a_direction, finalIdx * limitY1 * limitZ1, w[ 0 ] * weight * tof_weight * limitY1 * limitZ1);
779 ap_ProjectionLine->
AddVoxel(a_direction, ( finalIdx +
m_nbVoxXY ) * limitY1 * limitZ2, w[ 1 ] * weight * tof_weight * limitY1 * limitZ2);
780 ap_ProjectionLine->
AddVoxel(a_direction, ( finalIdx +
mp_nbVox[ 0 ] +
m_nbVoxXY ) * limitY2 * limitZ2, w[ 2 ] * weight * tof_weight * limitY2 * limitZ2);
781 ap_ProjectionLine->
AddVoxel(a_direction, ( finalIdx +
mp_nbVox[ 0 ] ) * limitY2 * limitZ1, w[ 3 ] * weight * tof_weight * limitY2 * limitZ1);
784 limitY1 = 1.0; limitY2 = 1.0;
785 limitZ1 = 1.0; limitZ2 = 1.0;
795 HPFLTNB const alphaX_0 = ( tof_edge_low[ 0 ] -
mp_sizeVox[ 0 ] * 0.5 - event1[ 0 ] )
797 HPFLTNB const alphaX_1 = ( tof_edge_high[ 0 ] +
mp_sizeVox[ 0 ] * 0.5 - event1[ 0 ] )
799 alphaMin = std::max( alphaMin, std::min( alphaX_0, alphaX_1 ) );
800 alphaMax = std::min( alphaMax, std::max( alphaX_0, alphaX_1 ) );
806 HPFLTNB const alphaY_0 = ( tof_edge_low[ 1 ] - event1[ 1 ] ) / r[ 1 ];
807 HPFLTNB const alphaY_1 = ( tof_edge_high[ 1 ] - event1[ 1 ] ) / r[ 1 ];
808 alphaMin = std::max( alphaMin, std::min( alphaY_0, alphaY_1 ) );
809 alphaMax = std::min( alphaMax, std::max( alphaY_0, alphaY_1 ) );
816 HPFLTNB const alphaZ_0 = ( tof_edge_low[ 2 ] -
mp_sizeVox[ 2 ] * 0.5 - event1[ 2 ] )
818 HPFLTNB const alphaZ_1 = ( tof_edge_high[ 2 ] +
mp_sizeVox[ 2 ] * 0.5 - event1[ 2 ] )
820 alphaMin = std::max( alphaMin, std::min( alphaZ_0, alphaZ_1 ) );
821 alphaMax = std::min( alphaMax, std::max( alphaZ_0, alphaZ_1 ) );
824 if( alphaMax <= alphaMin )
return 0;
826 if( r[ 0 ] == 0.0 )
if( event1[ 0 ] < -
mp_halfFOV[ 0 ] || event1[ 0 ] >
mp_halfFOV[ 0 ] )
return 0;
827 if( r[ 2 ] == 0.0 )
if( event1[ 2 ] < -
mp_halfFOV[ 2 ] || event1[ 2 ] >
mp_halfFOV[ 2 ] )
return 0;
830 HPFLTNB const factor( ::fabs( r[ 1 ] ) / lor_length );
831 HPFLTNB const factor_for_tof( lor_length / r[ 1 ] );
842 int jMin = 0, jMax = 0;
848 else if( r[ 1 ] < 0.0 )
863 for(
int j = jMin; j < jMax; ++j )
867 pos[ 0 ] = event1[ 0 ] + step * ri[ 0 ];
868 pos[ 2 ] = event1[ 2 ] + step * ri[ 2 ];
875 wfl[ 0 ] = pos[ 0 ] - (
m_boundX + index[ 0 ] * mp_sizeVox[ 0 ] );
876 wfl[ 1 ] = pos[ 2 ] - (
m_boundZ + index[ 1 ] * mp_sizeVox[ 2 ] );
877 wfl[ 0 ] /= mp_sizeVox[ 0 ];
878 wfl[ 1 ] /= mp_sizeVox[ 2 ];
881 wcl[ 0 ] = 1.0 - wfl[ 0 ];
882 wcl[ 1 ] = 1.0 - wfl[ 1 ];
885 w[ 0 ] = wcl[ 0 ] * wcl[ 1 ];
886 w[ 1 ] = wcl[ 0 ] * wfl[ 1 ];
887 w[ 2 ] = wfl[ 0 ] * wfl[ 1 ];
888 w[ 3 ] = wfl[ 0 ] * wcl[ 1 ];
891 if( index[ 0 ] <= 0 ) limitX1 = 0;
892 if( index[ 0 ] >=
mp_nbVox[ 0 ] ) limitX2 = 0;
895 if( index[ 1 ] <= 0 ) limitZ1 = 0;
896 if( index[ 1 ] >=
mp_nbVox[ 2 ] ) limitZ2 = 0;
899 finalIdx = ( index[ 0 ] - 1 ) + j *
mp_nbVox[ 0 ] + ( index[ 1 ] - 1 ) *
m_nbVoxXY;
923 HPFLTNB temp_erf = std::min(tof_half_span, std::max(-tof_half_span, lor_tof_center -
m_TOFBinSizeInMm/2. - step * factor_for_tof));
924 HPFLTNB prev_erf = erf(temp_erf/tof_sigma_sqrt2);
925 temp_erf = std::min(tof_half_span, std::max(-tof_half_span, lor_tof_center +
m_TOFBinSizeInMm/2. - step * factor_for_tof));
926 HPFLTNB new_erf = erf(temp_erf/tof_sigma_sqrt2);
933 HPFLTNB temp = ( step * factor_for_tof - lor_tof_center ) / tof_sigma;
941 ap_ProjectionLine->
AddVoxel(a_direction, finalIdx * limitX1 * limitZ1, w[ 0 ] * weight * tof_weight * limitX1 * limitZ1);
942 ap_ProjectionLine->
AddVoxel(a_direction, ( finalIdx +
m_nbVoxXY ) * limitX1 * limitZ2, w[ 1 ] * weight * tof_weight * limitX1 * limitZ2);
943 ap_ProjectionLine->
AddVoxel(a_direction, ( finalIdx + 1 +
m_nbVoxXY ) * limitX2 * limitZ2, w[ 2 ] * weight * tof_weight * limitX2 * limitZ2);
944 ap_ProjectionLine->
AddVoxel(a_direction, ( finalIdx + 1 ) * limitX2 * limitZ1, w[ 3 ] * weight * tof_weight * limitX2 * limitZ1);
947 limitX1 = 1; limitX2 = 1;
948 limitZ1 = 1; limitZ2 = 1;
964 Cerr(
"***** iProjectorJoseph::ProjectTOFHistogram() -> Called while not initialized !" << endl);
969 #ifdef CASTOR_VERBOSE 972 string direction =
"";
973 if (a_direction==
FORWARD) direction =
"forward";
974 else direction =
"backward";
975 Cout(
"iProjectorJoseph::Project with TOF bins -> Project line '" << ap_ProjectionLine <<
"' in " << direction <<
" direction" << endl);
985 event2[ 0 ] - event1[ 0 ],
986 event2[ 1 ] - event1[ 1 ],
987 event2[ 2 ] - event1[ 2 ]
992 HPFLTNB lor_length_half = lor_length * 0.5;
996 INTNB tof_half_nb_bins = tof_nb_bins/2;
1000 HPFLTNB tof_sigma_sqrt2 = sqrt(2.)*tof_sigma;
1006 HPFLTNB prev_erf = 0., new_erf = 0.;
1009 INTNB tof_bin_last = tof_half_nb_bins;
1010 INTNB tof_bin_first = -tof_half_nb_bins;
1018 INTNB tof_bin_first_for_voxel = 0, tof_bin_last_for_voxel = 0;
1031 HPFLTNB alphaMin = 0.0, alphaMax = 1.0;
1034 HPFLTNB pos[ 3 ] = { 0.0, 0.0, 0.0 };
1035 HPFLTNB wfl[ 2 ] = { 0.0, 0.0 };
1036 HPFLTNB wcl[ 2 ] = { 0.0, 0.0 };
1037 HPFLTNB w[ 4 ] = { 0.0, 0.0, 0.0, 0.0 };
1038 int16_t index[ 2 ] = { 0, 0 };
1039 int32_t finalIdx = 0;
1040 int8_t limitX1 = 1; int8_t limitX2 = 1;
1041 int8_t limitY1 = 1; int8_t limitY2 = 1;
1042 int8_t limitZ1 = 1; int8_t limitZ2 = 1;
1045 if( ::fabs( r[ 0 ] ) > ::fabs( r[ 1 ] ) )
1052 alphaMin = std::max( alphaMin, std::min( alphaX_0, alphaX_1 ) );
1053 alphaMax = std::min( alphaMax, std::max( alphaX_0, alphaX_1 ) );
1064 alphaMin = std::max( alphaMin, std::min( alphaY_0, alphaY_1 ) );
1065 alphaMax = std::min( alphaMax, std::max( alphaY_0, alphaY_1 ) );
1077 alphaMin = std::max( alphaMin, std::min( alphaZ_0, alphaZ_1 ) );
1078 alphaMax = std::min( alphaMax, std::max( alphaZ_0, alphaZ_1 ) );
1081 if( alphaMax <= alphaMin )
return 0;
1083 if( r[ 1 ] == 0.0 )
if( event1[ 1 ] < -
mp_halfFOV[ 1 ] || event1[ 1 ] >
mp_halfFOV[ 1 ] )
return 0;
1084 if( r[ 2 ] == 0.0 )
if( event1[ 2 ] < -
mp_halfFOV[ 2 ] || event1[ 2 ] >
mp_halfFOV[ 2 ] )
return 0;
1089 for (
INTNB tof_bin=0; tof_bin<tof_nb_bins; tof_bin++) tof_weights_temp[tof_bin] = 0.;
1092 HPFLTNB const factor( ::fabs( r[ 0 ] ) / lor_length );
1093 HPFLTNB const factor_for_tof( lor_length / r[ 0 ] );
1104 int iMin = 0, iMax = 0;
1110 else if( r[ 0 ] < 0.0 )
1120 for(
int i = iMin; i < iMax; ++i )
1124 pos[ 1 ] = event1[ 1 ] + step * ri[ 1 ];
1125 pos[ 2 ] = event1[ 2 ] + step * ri[ 2 ];
1131 finalIdx = i + ( index[ 0 ] - 1 ) *
mp_nbVox[ 0 ] + ( index[ 1 ] - 1 ) *
m_nbVoxXY;
1146 tof_bin_first_for_voxel = max(tof_bin_first , (
INTNB)(ceil( (step * factor_for_tof - tof_half_span - lor_length_half ) /
m_TOFBinSizeInMm )));
1147 tof_bin_last_for_voxel = min(tof_bin_last , (
INTNB)(floor( (step * factor_for_tof + tof_half_span - lor_length_half) /
m_TOFBinSizeInMm )));
1150 lor_tof_center = lor_length_half + tof_bin_first_for_voxel *
m_TOFBinSizeInMm;
1153 tof_bin_first_for_voxel += tof_half_nb_bins;
1154 tof_bin_last_for_voxel += tof_half_nb_bins;
1160 HPFLTNB temp = std::min(tof_half_span,std::max(-tof_half_span, lor_tof_center - m_TOFBinSizeInMm/2. - step * factor_for_tof));
1161 prev_erf = erf(temp/tof_sigma_sqrt2);
1166 for (
INTNB tof_bin=tof_bin_first_for_voxel; tof_bin<=tof_bin_last_for_voxel; tof_bin++)
1186 HPFLTNB temp = std::min(tof_half_span,std::max(-tof_half_span, lor_tof_center + m_TOFBinSizeInMm/2. - step * factor_for_tof));
1187 new_erf = erf(temp/tof_sigma_sqrt2);
1188 tof_weights_temp[tof_bin] = 0.5 * fabs(new_erf - prev_erf);
1198 HPFLTNB temp = ( step * factor_for_tof - lor_tof_center) / tof_sigma;
1210 wfl[ 0 ] = pos[ 1 ] - (
m_boundY + index[ 0 ] * mp_sizeVox[ 1 ] );
1211 wfl[ 1 ] = pos[ 2 ] - (
m_boundZ + index[ 1 ] * mp_sizeVox[ 2 ] );
1212 wfl[ 0 ] /= mp_sizeVox[ 1 ];
1213 wfl[ 1 ] /= mp_sizeVox[ 2 ];
1216 wcl[ 0 ] = 1.0 - wfl[ 0 ];
1217 wcl[ 1 ] = 1.0 - wfl[ 1 ];
1220 w[ 0 ] = wcl[ 0 ] * wcl[ 1 ];
1221 w[ 1 ] = wcl[ 0 ] * wfl[ 1 ];
1222 w[ 2 ] = wfl[ 0 ] * wfl[ 1 ];
1223 w[ 3 ] = wfl[ 0 ] * wcl[ 1 ];
1226 if( index[ 0 ] <= 0 ) limitY1 = 0;
1227 if( index[ 0 ] >=
mp_nbVox[ 1 ] ) limitY2 = 0;
1230 if( index[ 1 ] <= 0 ) limitZ1 = 0;
1231 if( index[ 1 ] >=
mp_nbVox[ 2 ] ) limitZ2 = 0;
1240 ap_ProjectionLine->
AddVoxelAllTOFBins(a_direction, finalIdx * limitY1 * limitZ1, w[ 0 ] * weight * limitY1 * limitZ1, tof_weights_temp, tof_bin_first_for_voxel, tof_bin_last_for_voxel);
1241 ap_ProjectionLine->
AddVoxelAllTOFBins(a_direction, ( finalIdx +
m_nbVoxXY ) * limitY1 * limitZ2, w[ 1 ] * weight * limitY1 * limitZ2 , tof_weights_temp, tof_bin_first_for_voxel, tof_bin_last_for_voxel);
1242 ap_ProjectionLine->
AddVoxelAllTOFBins(a_direction, ( finalIdx +
mp_nbVox[ 0 ] +
m_nbVoxXY ) * limitY2 * limitZ2, w[ 2 ] * weight * limitY2 * limitZ2 , tof_weights_temp, tof_bin_first_for_voxel, tof_bin_last_for_voxel);
1243 ap_ProjectionLine->
AddVoxelAllTOFBins(a_direction, ( finalIdx +
mp_nbVox[ 0 ] ) * limitY2 * limitZ1, w[ 3 ] * weight * limitY2 * limitZ1 , tof_weights_temp, tof_bin_first_for_voxel, tof_bin_last_for_voxel);
1246 limitY1 = 1.0; limitY2 = 1.0;
1247 limitZ1 = 1.0; limitZ2 = 1.0;
1249 delete[] tof_weights_temp;
1262 alphaMin = std::max( alphaMin, std::min( alphaX_0, alphaX_1 ) );
1263 alphaMax = std::min( alphaMax, std::max( alphaX_0, alphaX_1 ) );
1271 alphaMin = std::max( alphaMin, std::min( alphaY_0, alphaY_1 ) );
1272 alphaMax = std::min( alphaMax, std::max( alphaY_0, alphaY_1 ) );
1283 alphaMin = std::max( alphaMin, std::min( alphaZ_0, alphaZ_1 ) );
1284 alphaMax = std::min( alphaMax, std::max( alphaZ_0, alphaZ_1 ) );
1287 if( alphaMax <= alphaMin )
return 0;
1289 if( r[ 0 ] == 0.0 )
if( event1[ 0 ] < -
mp_halfFOV[ 0 ] || event1[ 0 ] >
mp_halfFOV[ 0 ] )
return 0;
1290 if( r[ 2 ] == 0.0 )
if( event1[ 2 ] < -
mp_halfFOV[ 2 ] || event1[ 2 ] >
mp_halfFOV[ 2 ] )
return 0;
1294 for (
INTNB tof_bin=0; tof_bin<tof_nb_bins; tof_bin++) tof_weights_temp[tof_bin] = 0.;
1297 HPFLTNB const factor( ::fabs( r[ 1 ] ) / lor_length );
1298 HPFLTNB const factor_for_tof( lor_length / r[ 1 ] );
1309 int jMin = 0, jMax = 0;
1315 else if( r[ 1 ] < 0.0 )
1325 for(
int j = jMin; j < jMax; ++j )
1329 pos[ 0 ] = event1[ 0 ] + step * ri[ 0 ];
1330 pos[ 2 ] = event1[ 2 ] + step * ri[ 2 ];
1336 finalIdx = ( index[ 0 ] - 1 ) + j *
mp_nbVox[ 0 ] + ( index[ 1 ] - 1 ) *
m_nbVoxXY;
1351 tof_bin_first_for_voxel = max(tof_bin_first , (
INTNB)(ceil( (step * factor_for_tof - tof_half_span - lor_length_half ) /
m_TOFBinSizeInMm )));
1352 tof_bin_last_for_voxel = min(tof_bin_last , (
INTNB)(floor( (step * factor_for_tof + tof_half_span - lor_length_half) /
m_TOFBinSizeInMm )));
1355 lor_tof_center = lor_length_half + tof_bin_first_for_voxel *
m_TOFBinSizeInMm;
1358 if(tof_bin_first_for_voxel>tof_bin_last || tof_bin_last_for_voxel<tof_bin_first)
continue;
1361 tof_bin_first_for_voxel += tof_half_nb_bins;
1362 tof_bin_last_for_voxel += tof_half_nb_bins;
1368 HPFLTNB temp = std::min(tof_half_span,std::max(-tof_half_span, lor_tof_center - m_TOFBinSizeInMm/2. - step * factor_for_tof));
1369 prev_erf = erf(temp/tof_sigma_sqrt2);
1374 for (
INTNB tof_bin=tof_bin_first_for_voxel; tof_bin<=tof_bin_last_for_voxel; tof_bin++)
1394 HPFLTNB temp = std::min(tof_half_span,std::max(-tof_half_span, lor_tof_center + m_TOFBinSizeInMm/2. - step * factor_for_tof));
1395 new_erf = erf(temp/tof_sigma_sqrt2);
1396 tof_weights_temp[tof_bin] = 0.5 * fabs(new_erf - prev_erf);
1406 HPFLTNB temp = ( step * factor_for_tof - lor_tof_center) / tof_sigma;
1418 wfl[ 0 ] = pos[ 0 ] - (
m_boundX + index[ 0 ] * mp_sizeVox[ 0 ] );
1419 wfl[ 1 ] = pos[ 2 ] - (
m_boundZ + index[ 1 ] * mp_sizeVox[ 2 ] );
1420 wfl[ 0 ] /= mp_sizeVox[ 0 ];
1421 wfl[ 1 ] /= mp_sizeVox[ 2 ];
1424 wcl[ 0 ] = 1.0 - wfl[ 0 ];
1425 wcl[ 1 ] = 1.0 - wfl[ 1 ];
1428 w[ 0 ] = wcl[ 0 ] * wcl[ 1 ];
1429 w[ 1 ] = wcl[ 0 ] * wfl[ 1 ];
1430 w[ 2 ] = wfl[ 0 ] * wfl[ 1 ];
1431 w[ 3 ] = wfl[ 0 ] * wcl[ 1 ];
1434 if( index[ 0 ] <= 0 ) limitX1 = 0;
1435 if( index[ 0 ] >=
mp_nbVox[ 0 ] ) limitX2 = 0;
1438 if( index[ 1 ] <= 0 ) limitZ1 = 0;
1439 if( index[ 1 ] >=
mp_nbVox[ 2 ] ) limitZ2 = 0;
1450 ap_ProjectionLine->
AddVoxelAllTOFBins(a_direction, finalIdx * limitX1 * limitZ1, w[ 0 ] * weight * limitX1 * limitZ1, tof_weights_temp, tof_bin_first_for_voxel, tof_bin_last_for_voxel);
1451 ap_ProjectionLine->
AddVoxelAllTOFBins(a_direction, ( finalIdx +
m_nbVoxXY ) * limitX1 * limitZ2, w[ 1 ] * weight * limitX1 * limitZ2, tof_weights_temp, tof_bin_first_for_voxel, tof_bin_last_for_voxel);
1452 ap_ProjectionLine->
AddVoxelAllTOFBins(a_direction, ( finalIdx + 1 +
m_nbVoxXY ) * limitX2 * limitZ2, w[ 2 ] * weight * limitX2 * limitZ2, tof_weights_temp, tof_bin_first_for_voxel, tof_bin_last_for_voxel);
1453 ap_ProjectionLine->
AddVoxelAllTOFBins(a_direction, ( finalIdx + 1 ) * limitX2 * limitZ1, w[ 3 ] * weight * limitX2 * limitZ1, tof_weights_temp, tof_bin_first_for_voxel, tof_bin_last_for_voxel);
1455 limitX1 = 1; limitX2 = 1;
1456 limitZ1 = 1; limitZ2 = 1;
1458 delete[] tof_weights_temp;
bool m_compatibleWithSPECTAttenuationCorrection
#define SPEED_OF_LIGHT_IN_MM_PER_PS
FLTNB m_TOFGaussianNormCoef
int InitializeSpecific()
This function is used to initialize specific stuff to the child projector.
void ShowHelpSpecific()
A function used to show help about the child module.
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
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.
int ReadConfigurationFile(const string &a_configurationFile)
A function used to read options from a configuration file.
bool m_TOFWeightingFcnPrecomputedFlag
bool m_TOFBinProperProcessingFlag
HPFLTNB * mp_TOFWeightingFcn
FLTNB GetLength()
This function is used to get the length of the line.
int ReadOptionsList(const string &a_optionsList)
A function used to read options from a list of options.
int CheckSpecificParameters()
A private function used to check the parameters settings specific to the child projector.
~iProjectorJoseph()
The destructor of iProjectorJoseph.
int ProjectTOFListmode(int a_direction, oProjectionLine *ap_ProjectionLine)
A function to project with TOF continuous information.
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...
FLTNB m_TOFResolutionInMm
INTNB EstimateMaxNumberOfVoxelsPerLine()
This function is used to compute and provide an estimate of the maximum number of voxels that could c...
Declaration of class iProjectorJoseph.
This class is designed to manage and store system matrix elements associated to a vEvent...
int ProjectWithoutTOF(int a_direction, oProjectionLine *ap_ProjectionLine)
A function to project without TOF.
FLTNB m_TOFPrecomputedSamplingFactor
Declaration of class sOutputManager.
int ProjectTOFHistogram(int a_direction, oProjectionLine *ap_ProjectionLine)
A function to project with TOF binned information.
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).
INTNB GetNbVoxX()
Get the number of voxels along the X axis.
INTNB GetNbVoxZ()
Get the number of voxels along the Z axis.
iProjectorJoseph()
The constructor of iProjectorJoseph.
bool m_compatibleWithCompression
INTNB GetNbVoxY()
Get the number of voxels along the Y axis.
INTNB m_TOFWeightingFcnNbSamples