105 cout <<
"This projector is a line projector based on computations of overlap between a pair of detection elements and voxels." << endl;
106 cout <<
"It is implemented from the following published paper:" << endl;
107 cout <<
"B. De Man and S. Basu, \"Distance-driven projection and backprojection in three dimensions\", Phys. Med. Biol., vol. 49, pp. 2463-75, 2004." << endl;
131 if (
m_verbose>=2)
Cout(
"iProjectorDistanceDriven::InitializeSpecific() -> Use Distance Driven projector" << endl);
134 HPFLTNB tolerance_factor = 1.e-4;
164 Cerr(
"***** iProjectorDistanceDriven::ProjectWithoutTOF() -> Called while not initialized !" << endl);
169 #ifdef CASTOR_VERBOSE 172 string direction =
"";
173 if (a_direction==
FORWARD) direction =
"forward";
174 else direction =
"backward";
175 Cout(
"iProjectorDistanceDriven::Project without TOF -> Project line '" << ap_ProjectionLine <<
"' in " << direction <<
" direction" << endl);
233 &p1_x[1], &p1_y[1], &p1_z[1],
234 &p2_x[1], &p2_y[1], &p2_z[1]
237 Cerr(
"***** iProjectorDistanceDriven::ProjectWithoutTOF() -> A problem occurred while getting the edges' center positions !" << endl);
248 p1_x[ 0 ] = p1[ 0 ]; p1_y[ 0 ] = p1[ 1 ]; p1_z[ 0 ] = p1[ 2 ];
251 p2_x[ 0 ] = p2[ 0 ]; p2_y[ 0 ] = p2[ 1 ]; p2_z[ 0 ] = p2[ 2 ];
260 for(
int p = 0; p < 5; ++p )
262 r[ p ][ 0 ] = p2_x[ p ] - p1_x[ p ];
263 r[ p ][ 1 ] = p2_y[ p ] - p1_y[ p ];
264 r[ p ][ 2 ] = p2_z[ p ] - p1_z[ p ];
269 r[ 0 ][ 0 ] * r[ 0 ][ 0 ],
270 r[ 0 ][ 1 ] * r[ 0 ][ 1 ],
271 r[ 0 ][ 2 ] * r[ 0 ][ 2 ]
276 HPFLTNB alpha_min = 0.0, alpha_max = 1.0;
280 HPFLTNB alpha_x_min[ 4 ], alpha_x_max[ 4 ];
281 HPFLTNB alpha_y_min[ 4 ], alpha_y_max[ 4 ];
282 HPFLTNB alpha_z_min[ 4 ], alpha_z_max[ 4 ];
307 float const angle = std::acos( ( r[ 0 ][ 0 ] ) / ( std::sqrt( r2[ 0 ] + r2[ 1 ] ) ) ) * 180.0f / M_PI;
310 if( ( angle >= 0.0 && angle <= 45.0 ) || ( angle >= 135.0 && angle <= 180.0 ) )
314 alpha_x_min[ 0 ] = ( (-
mp_halfFOV[ 0 ]) - p1_x[ 0 ] ) / r[ 0 ][ 0 ];
315 alpha_x_min[ 1 ] = (
mp_halfFOV[ 0 ] - p1_x[ 0 ] ) / r[ 0 ][ 0 ];
316 alpha_x_max[ 0 ] = alpha_x_min[ 0 ];
317 alpha_x_max[ 1 ] = alpha_x_min[ 1 ];
318 alpha_min = std::max( alpha_min, std::min( alpha_x_min[ 0 ], alpha_x_min[ 1 ] ) );
319 alpha_max = std::min( alpha_max, std::max( alpha_x_max[ 0 ], alpha_x_max[ 1 ] ) );
323 if( r[ 1 ][ 1 ] != 0.0 )
325 alpha_y_min[ 0 ] = ( (-
mp_halfFOV[ 1 ]) - p1_y[ 1 ] ) / r[ 1 ][ 1 ];
326 alpha_y_min[ 1 ] = (
mp_halfFOV[ 1 ] - p1_y[ 1 ] ) / r[ 1 ][ 1 ];
327 alpha_y_max[ 0 ] = alpha_y_min[ 0 ];
328 alpha_y_max[ 1 ] = alpha_y_min[ 1 ];
332 alpha_y_min[ 0 ] = 0.0;
333 alpha_y_min[ 1 ] = 0.0;
334 alpha_y_max[ 0 ] = 1.0;
335 alpha_y_max[ 1 ] = 1.0;
338 if( r[ 2 ][ 1 ] != 0.0 )
340 alpha_y_min[ 2 ] = ( (-
mp_halfFOV[ 1 ]) - p1_y[ 2 ] ) / r[ 2 ][ 1 ];
341 alpha_y_min[ 3 ] = (
mp_halfFOV[ 1 ] - p1_y[ 2 ] ) / r[ 2 ][ 1 ];
342 alpha_y_max[ 2 ] = alpha_y_min[ 2 ];
343 alpha_y_max[ 3 ] = alpha_y_min[ 3 ];
347 alpha_y_min[ 2 ] = 0.0;
348 alpha_y_min[ 3 ] = 0.0;
349 alpha_y_max[ 2 ] = 1.0;
350 alpha_y_max[ 3 ] = 1.0;
354 alpha_min = std::max( alpha_min,
355 alpha_y_min[ std::minmax_element( alpha_y_min, alpha_y_min + 4 ).first - alpha_y_min ] );
357 alpha_max = std::min( alpha_max,
358 alpha_y_max[ std::minmax_element( alpha_y_max, alpha_y_max + 4 ).second - alpha_y_max ] );
361 if( r[ 3 ][ 2 ] != 0.0 )
363 alpha_z_min[ 0 ] = ( (-
mp_halfFOV[ 2 ]) - p1_z[ 3 ] ) / r[ 3 ][ 2 ];
364 alpha_z_min[ 1 ] = (
mp_halfFOV[ 2 ] - p1_z[ 3 ] ) / r[ 3 ][ 2 ];
365 alpha_z_max[ 0 ] = alpha_z_min[ 0 ];
366 alpha_z_max[ 1 ] = alpha_z_min[ 1 ];
370 alpha_z_min[ 0 ] = 0.0;
371 alpha_z_min[ 1 ] = 0.0;
372 alpha_z_max[ 0 ] = 1.0;
373 alpha_z_max[ 1 ] = 1.0;
376 if( r[ 4 ][ 2 ] != 0.0 )
378 alpha_z_min[ 2 ] = ( (-
mp_halfFOV[ 2 ]) - p1_z[ 4 ] ) / r[ 4 ][ 2 ];
379 alpha_z_min[ 3 ] = (
mp_halfFOV[ 2 ] - p1_z[ 4 ] ) / r[ 4 ][ 2 ];
380 alpha_z_max[ 2 ] = alpha_z_min[ 2 ];
381 alpha_z_max[ 3 ] = alpha_z_min[ 3 ];
385 alpha_z_min[ 2 ] = 0.0;
386 alpha_z_min[ 3 ] = 0.0;
387 alpha_z_max[ 2 ] = 1.0;
388 alpha_z_max[ 3 ] = 1.0;
392 alpha_min = std::max( alpha_min,
393 alpha_z_min[ std::minmax_element( alpha_z_min, alpha_z_min + 4 ).first
396 alpha_max = std::min( alpha_max,
397 alpha_z_max[ std::minmax_element( alpha_z_max, alpha_z_max + 4 ).second
400 if( alpha_max <= alpha_min )
return 0;
403 int16_t i_min = 0, i_max = 0;
404 if( r[ 0 ][ 0 ] > 0.0 )
409 else if( r[ 0 ][ 0 ] < 0.0 )
417 HPFLTNB const factor( ::fabs( r[ 0 ][ 0 ] ) / ::sqrt( r2[ 0 ] + r2[ 1 ] + r2[ 2 ] ) );
422 for( uint8_t p = 0; p < 4; ++p )
425 ri[ p ][ 1 ] = r[ p + 1 ][ 1 ] / r[ p + 1 ][ 0 ];
426 ri[ p ][ 2 ] = r[ p + 1 ][ 2 ] / r[ p + 1 ][ 0 ];
430 int8_t incr_orient_trans = 0;
431 int8_t idx_orient_trans = 0;
434 int8_t
const incr_orient_axial = p2_z[ 3 ] < p2_z[ 4 ] ? 1 : -1;
438 int8_t
const idx_orient_axial = p2_z[ 3 ] < p2_z[ 4 ] ? 1 : 0;
444 for( int16_t i = i_min; i < i_max; ++i )
449 pos[ 0 ] = p1_y[ 1 ] + step * ri[ 0 ][ 1 ];
450 pos[ 1 ] = p1_y[ 2 ] + step * ri[ 1 ][ 1 ];
452 pos[ 2 ] = p1_z[ 3 ] + step * ri[ 2 ][ 2 ];
453 pos[ 3 ] = p1_z[ 4 ] + step * ri[ 3 ][ 2 ];
456 w1 = ::fabs( pos[ 0 ] - pos[ 1 ] );
457 w2 = ::fabs( pos[ 2 ] - pos[ 3 ] );
459 if( ( w1 * w2 ) != 0 ) final_weight = weight / ( w1 * w2 );
472 bool test_index = (index[ 0 ] >= 0 && index[ 2 ] >= 0 && index[ 1 ] >= 0 && index[ 3 ] >= 0
476 if( index[ 0 ] < index[ 1 ] )
478 incr_orient_trans = 1;
479 idx_orient_trans = 1;
481 else if( index[ 0 ] > index[ 1 ] )
483 incr_orient_trans = -1;
484 idx_orient_trans = 0;
488 int16_t
const n_distance_y = ::abs( index[ 0 ] - index[ 1 ] ) + 2;
490 distance_y[ 0 ] = pos[ 0 ];
491 distance_y[ n_distance_y - 1 ] = pos[ 1 ];
493 if( n_distance_y > 2 )
495 for( int16_t d = 0; d < n_distance_y - 2; ++d )
496 distance_y[ d + 1 ] = (-
mp_halfFOV[ 1 ]) + ( index[ 0 ] + idx_orient_trans + ( incr_orient_trans * d ) ) * mp_sizeVox[ 1 ];
500 for( int16_t d = 0; d < n_distance_y - 1; ++d )
501 distance_y[ d ] = ::fabs( distance_y[ d + 1 ] - distance_y[ d ] );
504 int16_t
const n_distance_z = ::abs( index[ 2 ] - index[ 3 ] ) + 2;
506 distance_z[ 0 ] = pos[ 2 ];
507 distance_z[ n_distance_z - 1 ] = pos[ 3 ];
509 if( n_distance_z > 2 )
511 for( int16_t d = 0; d < n_distance_z - 2; ++d )
512 distance_z[ d + 1 ] = (-
mp_halfFOV[ 2 ]) + ( index[ 2 ] + idx_orient_axial + ( incr_orient_axial * d ) ) * mp_sizeVox[ 2 ];
516 for( int16_t d = 0; d < n_distance_z - 1; ++d )
517 distance_z[ d ] = ::fabs( distance_z[ d + 1 ] - distance_z[ d ] );
520 int16_t index_y, index_z;
521 for( int16_t jj = 0; jj < n_distance_z - 1; ++jj )
523 index_z = index[ 2 ] + jj * incr_orient_axial;
524 if( index_z < 0 || index_z >
mp_nbVox[ 2 ] - 1 )
continue;
526 for( int16_t ii = 0; ii < n_distance_y - 1; ++ii )
528 index_y = index[ 0 ] + ii * incr_orient_trans;
529 if( index_y < 0 || index_y >
mp_nbVox[ 1 ] - 1 )
continue;
531 ap_ProjectionLine->
AddVoxel(a_direction, i + index_y *
mp_nbVox[ 0 ] + index_z * m_nbVoxXY, final_weight * distance_z[ jj ] * distance_y[ ii ]);
540 if( r[ 1 ][ 0 ] != 0.0 )
542 alpha_x_min[ 0 ] = ( (-
mp_halfFOV[ 0 ]) - p1_x[ 1 ] ) / r[ 1 ][ 0 ];
543 alpha_x_min[ 1 ] = (
mp_halfFOV[ 0 ] - p1_x[ 1 ] ) / r[ 1 ][ 0 ];
544 alpha_x_max[ 0 ] = alpha_x_min[ 0 ];
545 alpha_x_max[ 1 ] = alpha_x_min[ 1 ];
549 alpha_x_min[ 0 ] = 0.0;
550 alpha_x_min[ 1 ] = 0.0;
551 alpha_x_max[ 0 ] = 1.0;
552 alpha_x_max[ 1 ] = 1.0;
555 if( r[ 2 ][ 0 ] != 0.0 )
557 alpha_x_min[ 2 ] = ( (-
mp_halfFOV[ 0 ]) - p1_x[ 2 ] ) / r[ 2 ][ 0 ];
558 alpha_x_min[ 3 ] = (
mp_halfFOV[ 0 ] - p1_x[ 2 ] ) / r[ 2 ][ 0 ];
559 alpha_x_max[ 2 ] = alpha_x_min[ 2 ];
560 alpha_x_max[ 3 ] = alpha_x_min[ 3 ];
564 alpha_x_min[ 2 ] = 0.0;
565 alpha_x_min[ 3 ] = 0.0;
566 alpha_x_max[ 2 ] = 1.0;
567 alpha_x_max[ 3 ] = 1.0;
571 alpha_min = std::max( alpha_min,
572 alpha_x_min[ std::minmax_element( alpha_x_min, alpha_x_min + 4 ).first - alpha_x_min ] );
574 alpha_max = std::min( alpha_max,
575 alpha_x_max[ std::minmax_element( alpha_x_max, alpha_x_max + 4 ).second - alpha_x_max ] );
579 alpha_y_min[ 0 ] = ( (-
mp_halfFOV[ 1 ]) - p1_y[ 0 ] ) / r[ 0 ][ 1 ];
580 alpha_y_min[ 1 ] = (
mp_halfFOV[ 1 ] - p1_y[ 0 ] ) / r[ 0 ][ 1 ];
581 alpha_y_max[ 0 ] = alpha_y_min[ 0 ];
582 alpha_y_max[ 1 ] = alpha_y_min[ 1 ];
583 alpha_min = std::max( alpha_min,
584 std::min( alpha_y_min[ 0 ], alpha_y_min[ 1 ] ) );
585 alpha_max = std::min( alpha_max,
586 std::max( alpha_y_max[ 0 ], alpha_y_max[ 1 ] ) );
589 if( r[ 3 ][ 2 ] != 0.0 )
591 alpha_z_min[ 0 ] = ( (-
mp_halfFOV[ 2 ]) - p1_z[ 3 ] ) / r[ 3 ][ 2 ];
592 alpha_z_min[ 1 ] = (
mp_halfFOV[ 2 ] - p1_z[ 3 ] ) / r[ 3 ][ 2 ];
593 alpha_z_max[ 0 ] = alpha_z_min[ 0 ];
594 alpha_z_max[ 1 ] = alpha_z_min[ 1 ];
598 alpha_z_min[ 0 ] = 0.0;
599 alpha_z_min[ 1 ] = 0.0;
600 alpha_z_max[ 0 ] = 1.0;
601 alpha_z_max[ 1 ] = 1.0;
604 if( r[ 4 ][ 2 ] != 0.0 )
606 alpha_z_min[ 2 ] = ( (-
mp_halfFOV[ 2 ]) - p1_z[ 4 ] ) / r[ 4 ][ 2 ];
607 alpha_z_min[ 3 ] = (
mp_halfFOV[ 2 ] - p1_z[ 4 ] ) / r[ 4 ][ 2 ];
608 alpha_z_max[ 2 ] = alpha_z_min[ 2 ];
609 alpha_z_max[ 3 ] = alpha_z_min[ 3 ];
613 alpha_z_min[ 2 ] = 0.0;
614 alpha_z_min[ 3 ] = 0.0;
615 alpha_z_max[ 2 ] = 1.0;
616 alpha_z_max[ 3 ] = 1.0;
620 alpha_min = std::max( alpha_min, alpha_z_min[ std::minmax_element( alpha_z_min, alpha_z_min + 4 ).first - alpha_z_min ] );
622 alpha_max = std::min( alpha_max, alpha_z_max[ std::minmax_element( alpha_z_max, alpha_z_max + 4 ).second - alpha_z_max ] );
624 if( alpha_max <= alpha_min )
return 0;
627 int16_t j_min = 0, j_max = 0;
628 if( r[ 0 ][ 1 ] > 0.0 )
633 else if( r[ 0 ][ 1 ] < 0.0 )
641 HPFLTNB const factor( ::fabs( r[ 0 ][ 1 ] ) / ::sqrt( r2[ 0 ] + r2[ 1 ] + r2[ 2 ] ) );
646 for( uint8_t p = 0; p < 4; ++p )
648 ri[ p ][ 0 ] = r[ p + 1 ][ 0 ] / r[ p + 1 ][ 1 ];
650 ri[ p ][ 2 ] = r[ p + 1 ][ 2 ] / r[ p + 1 ][ 1 ];
654 int8_t incr_orient_trans = 0;
655 int8_t idx_orient_trans = 0;
658 int8_t
const incr_orient_axial = p2_z[ 3 ] < p2_z[ 4 ] ? 1 : -1;
660 int8_t
const idx_orient_axial = p2_z[ 3 ] < p2_z[ 4 ] ? 1 : 0;
666 for( int16_t j = j_min; j < j_max; ++j )
671 pos[ 0 ] = p1_x[ 1 ] + step * ri[ 0 ][ 0 ];
672 pos[ 1 ] = p1_x[ 2 ] + step * ri[ 1 ][ 0 ];
674 pos[ 2 ] = p1_z[ 3 ] + step * ri[ 2 ][ 2 ];
675 pos[ 3 ] = p1_z[ 4 ] + step * ri[ 3 ][ 2 ];
678 w1 = ::fabs( pos[ 0 ] - pos[ 1 ] );
679 w2 = ::fabs( pos[ 2 ] - pos[ 3 ] );
681 if( ( w1 * w2 ) != 0 ) final_weight = weight / ( w1 * w2 );
694 bool test_index = (index[ 0 ] >= 0 && index[ 2 ] >= 0 && index[ 1 ] >= 0 && index[ 3 ] >= 0
698 if( index[ 0 ] < index[ 1 ] )
700 incr_orient_trans = 1;
701 idx_orient_trans = 1;
703 else if( index[ 0 ] > index[ 1 ] )
705 incr_orient_trans = -1;
706 idx_orient_trans = 0;
710 int16_t
const n_distance_x = ::abs( index[ 0 ] - index[ 1 ] ) + 2;
712 distance_x[ 0 ] = pos[ 0 ];
713 distance_x[ n_distance_x - 1 ] = pos[ 1 ];
715 if( n_distance_x > 2 )
717 for( int16_t d = 0; d < n_distance_x - 2; ++d )
718 distance_x[ d + 1 ] = (-
mp_halfFOV[ 0 ]) + ( index[ 0 ] + idx_orient_trans + ( incr_orient_trans * d ) ) * mp_sizeVox[ 0 ];
722 for( int16_t d = 0; d < n_distance_x - 1; ++d )
723 distance_x[ d ] = ::fabs( distance_x[ d + 1 ] - distance_x[ d ] );
726 int16_t
const n_distance_z = ::abs( index[ 2 ] - index[ 3 ] ) + 2;
728 distance_z[ 0 ] = pos[ 2 ];
729 distance_z[ n_distance_z - 1 ] = pos[ 3 ];
731 if( n_distance_z > 2 )
733 for( int16_t d = 0; d < n_distance_z - 2; ++d )
734 distance_z[ d + 1 ] = (-
mp_halfFOV[ 2 ]) + ( index[ 2 ] + idx_orient_axial + ( incr_orient_axial * d ) ) * mp_sizeVox[ 2 ];
738 for( int16_t d = 0; d < n_distance_z - 1; ++d )
739 distance_z[ d ] = ::fabs( distance_z[ d + 1 ] - distance_z[ d ] );
742 int16_t index_x, index_z;
743 for( int16_t jj = 0; jj < n_distance_z - 1; ++jj )
745 index_z = index[ 2 ] + jj * incr_orient_axial;
746 if( index_z < 0 || index_z >
mp_nbVox[ 2 ] - 1 )
continue;
748 for( int16_t ii = 0; ii < n_distance_x - 1; ++ii )
750 index_x = index[ 0 ] + ii * incr_orient_trans;
751 if( index_x < 0 || index_x >
mp_nbVox[ 0 ] - 1 )
continue;
753 ap_ProjectionLine->
AddVoxel(a_direction, index_x + j *
mp_nbVox[ 0 ] + index_z * m_nbVoxXY, final_weight * distance_z[ jj ] * distance_x[ ii ]);
772 Cerr(
"***** iProjectorDistanceDriven::ProjectTOFListmode() -> Called while not initialized !" << endl);
777 #ifdef CASTOR_VERBOSE 780 string direction =
"";
781 if (a_direction==
FORWARD) direction =
"forward";
782 else direction =
"backward";
783 Cout(
"iProjectorDistanceDriven::Project with TOF measurement -> Project line '" << ap_ProjectionLine <<
"' in " << direction <<
" direction" << endl);
840 &p1_x[1], &p1_y[1], &p1_z[1],
841 &p2_x[1], &p2_y[1], &p2_z[1]
844 Cerr(
"***** iProjectorDistanceDriven::ProjectWithoutTOF() -> A problem occurred while getting the edges' center positions !" << endl);
855 p1_x[ 0 ] = p1[ 0 ]; p1_y[ 0 ] = p1[ 1 ]; p1_z[ 0 ] = p1[ 2 ];
858 p2_x[ 0 ] = p2[ 0 ]; p2_y[ 0 ] = p2[ 1 ]; p2_z[ 0 ] = p2[ 2 ];
867 for(
int p = 0; p < 5; ++p )
869 r[ p ][ 0 ] = p2_x[ p ] - p1_x[ p ];
870 r[ p ][ 1 ] = p2_y[ p ] - p1_y[ p ];
871 r[ p ][ 2 ] = p2_z[ p ] - p1_z[ p ];
876 r[ 0 ][ 0 ] * r[ 0 ][ 0 ],
877 r[ 0 ][ 1 ] * r[ 0 ][ 1 ],
878 r[ 0 ][ 2 ] * r[ 0 ][ 2 ]
883 HPFLTNB alpha_min = 0.0, alpha_max = 1.0;
887 HPFLTNB alpha_x_min[ 4 ], alpha_x_max[ 4 ];
888 HPFLTNB alpha_y_min[ 4 ], alpha_y_max[ 4 ];
889 HPFLTNB alpha_z_min[ 4 ], alpha_z_max[ 4 ];
923 HPFLTNB tof_sigma_sqrt2 = sqrt(2.)*tof_sigma;
931 HPFLTNB lor_tof_center = length_LOR * 0.5 + tof_delta;
934 HPFLTNB tof_edge_low[] = {0,0,0};
936 HPFLTNB tof_edge_high[] = {0,0,0};
937 HPFLTNB tof_center[] = {0,0,0};
941 for (
int ax=0;ax<3;ax++)
944 tof_center[ax] = p1[ax] + lor_tof_center * r[0][ax] / length_LOR;
947 tof_edge_low[ax] = tof_center[ax] - tof_half_span * fabs(r[0][ax]) / length_LOR;
950 if (tof_index>
mp_nbVox[ax]-1)
return 0;
954 tof_edge_high[ax] = tof_center[ax] + tof_half_span * fabs(r[0][ax]) / length_LOR;
957 if (tof_index<0)
return 0;
962 float const angle = std::acos( ( r[ 0 ][ 0 ] ) / ( std::sqrt( r2[ 0 ] + r2[ 1 ] ) ) ) * 180.0f / M_PI;
965 if( ( angle >= 0.0 && angle <= 45.0 ) || ( angle >= 135.0 && angle <= 180.0 ) )
970 alpha_x_min[ 0 ] = ( tof_edge_low[ 0 ] - p1_x[ 0 ] ) / r[ 0 ][ 0 ];
971 alpha_x_min[ 1 ] = ( tof_edge_high[ 0 ] - p1_x[ 0 ] ) / r[ 0 ][ 0 ];
972 alpha_x_max[ 0 ] = alpha_x_min[ 0 ];
973 alpha_x_max[ 1 ] = alpha_x_min[ 1 ];
974 alpha_min = std::max( alpha_min, std::min( alpha_x_min[ 0 ], alpha_x_min[ 1 ] ) );
975 alpha_max = std::min( alpha_max, std::max( alpha_x_max[ 0 ], alpha_x_max[ 1 ] ) );
979 if( r[ 1 ][ 1 ] != 0 )
981 alpha_y_min[ 0 ] = ( tof_edge_low[ 1 ] - p1_y[ 1 ] ) / r[ 1 ][ 1 ];
982 alpha_y_min[ 1 ] = ( tof_edge_high[ 1 ] - p1_y[ 1 ] ) / r[ 1 ][ 1 ];
983 alpha_y_max[ 0 ] = alpha_y_min[ 0 ];
984 alpha_y_max[ 1 ] = alpha_y_min[ 1 ];
988 alpha_y_min[ 0 ] = 0.0;
989 alpha_y_min[ 1 ] = 0.0;
990 alpha_y_max[ 0 ] = 1.0;
991 alpha_y_max[ 1 ] = 1.0;
994 if( r[ 2 ][ 1 ] != 0 )
996 alpha_y_min[ 2 ] = ( tof_edge_low[ 1 ] - p1_y[ 2 ] ) / r[ 2 ][ 1 ];
997 alpha_y_min[ 3 ] = ( tof_edge_high[ 1 ] - p1_y[ 2 ] ) / r[ 2 ][ 1 ];
998 alpha_y_max[ 2 ] = alpha_y_min[ 2 ];
999 alpha_y_max[ 3 ] = alpha_y_min[ 3 ];
1003 alpha_y_min[ 2 ] = 0.0;
1004 alpha_y_min[ 3 ] = 0.0;
1005 alpha_y_max[ 2 ] = 1.0;
1006 alpha_y_max[ 3 ] = 1.0;
1010 alpha_min = std::max( alpha_min,
1011 alpha_y_min[ std::minmax_element( alpha_y_min, alpha_y_min + 4 ).first - alpha_y_min ] );
1013 alpha_max = std::min( alpha_max,
1014 alpha_y_max[ std::minmax_element( alpha_y_max, alpha_y_max + 4 ).second - alpha_y_max ] );
1017 if( r[ 3 ][ 2 ] != 0 )
1019 alpha_z_min[ 0 ] = ( tof_edge_low[ 2 ] - p1_z[ 3 ] ) / r[ 3 ][ 2 ];
1020 alpha_z_min[ 1 ] = ( tof_edge_high[ 2 ] - p1_z[ 3 ] ) / r[ 3 ][ 2 ];
1021 alpha_z_max[ 0 ] = alpha_z_min[ 0 ];
1022 alpha_z_max[ 1 ] = alpha_z_min[ 1 ];
1026 alpha_z_min[ 0 ] = 0.0;
1027 alpha_z_min[ 1 ] = 0.0;
1028 alpha_z_max[ 0 ] = 1.0;
1029 alpha_z_max[ 1 ] = 1.0;
1032 if( r[ 4 ][ 2 ] != 0 )
1034 alpha_z_min[ 2 ] = ( tof_edge_low[ 2 ] - p1_z[ 4 ] ) / r[ 4 ][ 2 ];
1035 alpha_z_min[ 3 ] = ( tof_edge_high[ 2 ] - p1_z[ 4 ] ) / r[ 4 ][ 2 ];
1036 alpha_z_max[ 2 ] = alpha_z_min[ 2 ];
1037 alpha_z_max[ 3 ] = alpha_z_min[ 3 ];
1041 alpha_z_min[ 2 ] = 0.0;
1042 alpha_z_min[ 3 ] = 0.0;
1043 alpha_z_max[ 2 ] = 1.0;
1044 alpha_z_max[ 3 ] = 1.0;
1048 alpha_min = std::max( alpha_min,
1049 alpha_z_min[ std::minmax_element( alpha_z_min, alpha_z_min + 4 ).first
1052 alpha_max = std::min( alpha_max,
1053 alpha_z_max[ std::minmax_element( alpha_z_max, alpha_z_max + 4 ).second
1056 if( alpha_max <= alpha_min )
return 0;
1059 int16_t i_min = 0, i_max = 0;
1060 if( r[ 0 ][ 0 ] > 0.0 )
1065 else if( r[ 0 ][ 0 ] < 0.0 )
1073 HPFLTNB const factor( ::fabs( r[ 0 ][ 0 ] ) / ::sqrt( r2[ 0 ] + r2[ 1 ] + r2[ 2 ] ) );
1075 HPFLTNB const factor_for_tof( length_LOR / r[ 0 ][ 0 ] );
1079 for( uint8_t p = 0; p < 4; ++p )
1082 ri[ p ][ 1 ] = r[ p + 1 ][ 1 ] / r[ p + 1 ][ 0 ];
1083 ri[ p ][ 2 ] = r[ p + 1 ][ 2 ] / r[ p + 1 ][ 0 ];
1087 int8_t incr_orient_trans = 0;
1088 int8_t idx_orient_trans = 0;
1091 int8_t
const incr_orient_axial = p2_z[ 3 ] < p2_z[ 4 ] ? 1 : -1;
1095 int8_t
const idx_orient_axial = p2_z[ 3 ] < p2_z[ 4 ] ? 1 : 0;
1101 for( int16_t i = i_min; i < i_max; ++i )
1106 pos[ 0 ] = p1_y[ 1 ] + step * ri[ 0 ][ 1 ];
1107 pos[ 1 ] = p1_y[ 2 ] + step * ri[ 1 ][ 1 ];
1109 pos[ 2 ] = p1_z[ 3 ] + step * ri[ 2 ][ 2 ];
1110 pos[ 3 ] = p1_z[ 4 ] + step * ri[ 3 ][ 2 ];
1113 w1 = ::fabs( pos[ 0 ] - pos[ 1 ] );
1114 w2 = ::fabs( pos[ 2 ] - pos[ 3 ] );
1116 if( ( w1 * w2 ) != 0 ) final_weight = weight / ( w1 * w2 );
1129 bool test_index = (index[ 0 ] >= 0 && index[ 2 ] >= 0 && index[ 1 ] >= 0 && index[ 3 ] >= 0
1133 if( index[ 0 ] < index[ 1 ] )
1135 incr_orient_trans = 1;
1136 idx_orient_trans = 1;
1138 else if( index[ 0 ] > index[ 1 ] )
1140 incr_orient_trans = -1;
1141 idx_orient_trans = 0;
1145 int16_t
const n_distance_y = ::abs( index[ 0 ] - index[ 1 ] ) + 2;
1147 distance_y[ 0 ] = pos[ 0 ];
1148 distance_y[ n_distance_y - 1 ] = pos[ 1 ];
1150 if( n_distance_y > 2 )
1152 for( int16_t d = 0; d < n_distance_y - 2; ++d )
1153 distance_y[ d + 1 ] = (-
mp_halfFOV[ 1 ]) + ( index[ 0 ] + idx_orient_trans + ( incr_orient_trans * d ) ) * mp_sizeVox[ 1 ];
1157 for( int16_t d = 0; d < n_distance_y - 1; ++d )
1158 distance_y[ d ] = ::fabs( distance_y[ d + 1 ] - distance_y[ d ] );
1161 int16_t
const n_distance_z = ::abs( index[ 2 ] - index[ 3 ] ) + 2;
1163 distance_z[ 0 ] = pos[ 2 ];
1164 distance_z[ n_distance_z - 1 ] = pos[ 3 ];
1166 if( n_distance_z > 2 )
1168 for( int16_t d = 0; d < n_distance_z - 2; ++d )
1169 distance_z[ d + 1 ] = (-
mp_halfFOV[ 2 ]) + ( index[ 2 ] + idx_orient_axial + ( incr_orient_axial * d ) ) * mp_sizeVox[ 2 ];
1173 for( int16_t d = 0; d < n_distance_z - 1; ++d )
1174 distance_z[ d ] = ::fabs( distance_z[ d + 1 ] - distance_z[ d ] );
1191 HPFLTNB temp_erf = std::min(tof_half_span, std::max(-tof_half_span, lor_tof_center -
m_TOFBinSizeInMm/2. - step * factor_for_tof));
1192 HPFLTNB prev_erf = erf(temp_erf/tof_sigma_sqrt2);
1193 temp_erf = std::min(tof_half_span, std::max(-tof_half_span, lor_tof_center +
m_TOFBinSizeInMm/2. - step * factor_for_tof));
1194 HPFLTNB new_erf = erf(temp_erf/tof_sigma_sqrt2);
1201 HPFLTNB temp = (step * factor_for_tof - lor_tof_center) / tof_sigma;
1207 int16_t index_y, index_z;
1208 for( int16_t jj = 0; jj < n_distance_z - 1; ++jj )
1210 index_z = index[ 2 ] + jj * incr_orient_axial;
1211 if( index_z < 0 || index_z >
mp_nbVox[ 2 ] - 1 )
continue;
1213 for( int16_t ii = 0; ii < n_distance_y - 1; ++ii )
1215 index_y = index[ 0 ] + ii * incr_orient_trans;
1216 if( index_y < 0 || index_y >
mp_nbVox[ 1 ] - 1 )
continue;
1218 ap_ProjectionLine->
AddVoxel(a_direction, i + index_y *
mp_nbVox[ 0 ] + index_z * m_nbVoxXY, final_weight * distance_z[ jj ] * distance_y[ ii ] * tof_weight);
1228 if( r[ 1 ][ 0 ] != 0.0 )
1230 alpha_x_min[ 0 ] = ( tof_edge_low[ 0 ] - p1_x[ 1 ] ) / r[ 1 ][ 0 ];
1231 alpha_x_min[ 1 ] = ( tof_edge_high[ 0 ] - p1_x[ 1 ] ) / r[ 1 ][ 0 ];
1232 alpha_x_max[ 0 ] = alpha_x_min[ 0 ];
1233 alpha_x_max[ 1 ] = alpha_x_min[ 1 ];
1237 alpha_x_min[ 0 ] = 0.0;
1238 alpha_x_min[ 1 ] = 0.0;
1239 alpha_x_max[ 0 ] = 1.0;
1240 alpha_x_max[ 1 ] = 1.0;
1243 if( r[ 2 ][ 0 ] != 0 )
1245 alpha_x_min[ 2 ] = ( tof_edge_low[ 0 ] - p1_x[ 2 ] ) / r[ 2 ][ 0 ];
1246 alpha_x_min[ 3 ] = ( tof_edge_high[ 0 ] - p1_x[ 2 ] ) / r[ 2 ][ 0 ];
1247 alpha_x_max[ 2 ] = alpha_x_min[ 2 ];
1248 alpha_x_max[ 3 ] = alpha_x_min[ 3 ];
1252 alpha_x_min[ 2 ] = 0.0;
1253 alpha_x_min[ 3 ] = 0.0;
1254 alpha_x_max[ 2 ] = 1.0;
1255 alpha_x_max[ 3 ] = 1.0;
1259 alpha_min = std::max( alpha_min,
1260 alpha_x_min[ std::minmax_element( alpha_x_min, alpha_x_min + 4 ).first - alpha_x_min ] );
1262 alpha_max = std::min( alpha_max,
1263 alpha_x_max[ std::minmax_element( alpha_x_max, alpha_x_max + 4 ).second - alpha_x_max ] );
1267 alpha_y_min[ 0 ] = ( tof_edge_low[1] - p1_y[ 0 ] ) / r[ 0 ][ 1 ];
1268 alpha_y_min[ 1 ] = ( tof_edge_high[1] - p1_y[ 0 ] ) / r[ 0 ][ 1 ];
1269 alpha_y_max[ 0 ] = alpha_y_min[ 0 ];
1270 alpha_y_max[ 1 ] = alpha_y_min[ 1 ];
1271 alpha_min = std::max( alpha_min,
1272 std::min( alpha_y_min[ 0 ], alpha_y_min[ 1 ] ) );
1273 alpha_max = std::min( alpha_max,
1274 std::max( alpha_y_max[ 0 ], alpha_y_max[ 1 ] ) );
1277 if( r[ 3 ][ 2 ] != 0 )
1279 alpha_z_min[ 0 ] = ( tof_edge_low[ 2 ] - p1_z[ 3 ] ) / r[ 3 ][ 2 ];
1280 alpha_z_min[ 1 ] = ( tof_edge_high[ 2 ]- p1_z[ 3 ] ) / r[ 3 ][ 2 ];
1281 alpha_z_max[ 0 ] = alpha_z_min[ 0 ];
1282 alpha_z_max[ 1 ] = alpha_z_min[ 1 ];
1286 alpha_z_min[ 0 ] = 0.0;
1287 alpha_z_min[ 1 ] = 0.0;
1288 alpha_z_max[ 0 ] = 1.0;
1289 alpha_z_max[ 1 ] = 1.0;
1292 if( r[ 4 ][ 2 ] != 0 )
1294 alpha_z_min[ 2 ] = ( tof_edge_low[ 2 ] - p1_z[ 4 ] ) / r[ 4 ][ 2 ];
1295 alpha_z_min[ 3 ] = ( tof_edge_high[ 2 ] - p1_z[ 4 ] ) / r[ 4 ][ 2 ];
1296 alpha_z_max[ 2 ] = alpha_z_min[ 2 ];
1297 alpha_z_max[ 3 ] = alpha_z_min[ 3 ];
1301 alpha_z_min[ 2 ] = 0.0;
1302 alpha_z_min[ 3 ] = 0.0;
1303 alpha_z_max[ 2 ] = 1.0;
1304 alpha_z_max[ 3 ] = 1.0;
1308 alpha_min = std::max( alpha_min, alpha_z_min[ std::minmax_element( alpha_z_min, alpha_z_min + 4 ).first - alpha_z_min ] );
1310 alpha_max = std::min( alpha_max, alpha_z_max[ std::minmax_element( alpha_z_max, alpha_z_max + 4 ).second - alpha_z_max ] );
1312 if( alpha_max <= alpha_min )
return 0;
1315 int16_t j_min = 0, j_max = 0;
1316 if( r[ 0 ][ 1 ] > 0.0 )
1321 else if( r[ 0 ][ 1 ] < 0.0 )
1329 HPFLTNB const factor( ::fabs( r[ 0 ][ 1 ] ) / ::sqrt( r2[ 0 ] + r2[ 1 ] + r2[ 2 ] ) );
1331 HPFLTNB const factor_for_tof( length_LOR / r[ 0 ][ 1 ] );
1335 for( uint8_t p = 0; p < 4; ++p )
1337 ri[ p ][ 0 ] = r[ p + 1 ][ 0 ] / r[ p + 1 ][ 1 ];
1339 ri[ p ][ 2 ] = r[ p + 1 ][ 2 ] / r[ p + 1 ][ 1 ];
1343 int8_t incr_orient_trans = 0;
1344 int8_t idx_orient_trans = 0;
1347 int8_t
const incr_orient_axial = p2_z[ 3 ] < p2_z[ 4 ] ? 1 : -1;
1349 int8_t
const idx_orient_axial = p2_z[ 3 ] < p2_z[ 4 ] ? 1 : 0;
1355 for( int16_t j = j_min; j < j_max; ++j )
1360 pos[ 0 ] = p1_x[ 1 ] + step * ri[ 0 ][ 0 ];
1361 pos[ 1 ] = p1_x[ 2 ] + step * ri[ 1 ][ 0 ];
1363 pos[ 2 ] = p1_z[ 3 ] + step * ri[ 2 ][ 2 ];
1364 pos[ 3 ] = p1_z[ 4 ] + step * ri[ 3 ][ 2 ];
1367 w1 = ::fabs( pos[ 0 ] - pos[ 1 ] );
1368 w2 = ::fabs( pos[ 2 ] - pos[ 3 ] );
1370 if( ( w1 * w2 ) != 0 ) final_weight = weight / ( w1 * w2 );
1383 bool test_index = (index[ 0 ] >= 0 && index[ 2 ] >= 0 && index[ 1 ] >= 0 && index[ 3 ] >= 0
1387 if( index[ 0 ] < index[ 1 ] )
1389 incr_orient_trans = 1;
1390 idx_orient_trans = 1;
1392 else if( index[ 0 ] > index[ 1 ] )
1394 incr_orient_trans = -1;
1395 idx_orient_trans = 0;
1399 int16_t
const n_distance_x = ::abs( index[ 0 ] - index[ 1 ] ) + 2;
1401 distance_x[ 0 ] = pos[ 0 ];
1402 distance_x[ n_distance_x - 1 ] = pos[ 1 ];
1404 if( n_distance_x > 2 )
1406 for( int16_t d = 0; d < n_distance_x - 2; ++d )
1407 distance_x[ d + 1 ] = (-
mp_halfFOV[ 0 ]) + ( index[ 0 ] + idx_orient_trans + ( incr_orient_trans * d ) ) * mp_sizeVox[ 0 ];
1411 for( int16_t d = 0; d < n_distance_x - 1; ++d )
1412 distance_x[ d ] = ::fabs( distance_x[ d + 1 ] - distance_x[ d ] );
1415 int16_t
const n_distance_z = ::abs( index[ 2 ] - index[ 3 ] ) + 2;
1417 distance_z[ 0 ] = pos[ 2 ];
1418 distance_z[ n_distance_z - 1 ] = pos[ 3 ];
1420 if( n_distance_z > 2 )
1422 for( int16_t d = 0; d < n_distance_z - 2; ++d )
1423 distance_z[ d + 1 ] = (-
mp_halfFOV[ 2 ]) + ( index[ 2 ] + idx_orient_axial + ( incr_orient_axial * d ) ) * mp_sizeVox[ 2 ];
1427 for( int16_t d = 0; d < n_distance_z - 1; ++d )
1428 distance_z[ d ] = ::fabs( distance_z[ d + 1 ] - distance_z[ d ] );
1445 HPFLTNB temp_erf = std::min(tof_half_span, std::max(-tof_half_span, lor_tof_center -
m_TOFBinSizeInMm/2. - step * factor_for_tof));
1446 HPFLTNB prev_erf = erf(temp_erf/tof_sigma_sqrt2);
1447 temp_erf = std::min(tof_half_span, std::max(-tof_half_span, lor_tof_center +
m_TOFBinSizeInMm/2. - step * factor_for_tof));
1448 HPFLTNB new_erf = erf(temp_erf/tof_sigma_sqrt2);
1455 HPFLTNB temp = (step * factor_for_tof - lor_tof_center) / tof_sigma;
1461 int16_t index_x, index_z;
1462 for( int16_t jj = 0; jj < n_distance_z - 1; ++jj )
1464 index_z = index[ 2 ] + jj * incr_orient_axial;
1465 if( index_z < 0 || index_z >
mp_nbVox[ 2 ] - 1 )
continue;
1467 for( int16_t ii = 0; ii < n_distance_x - 1; ++ii )
1469 index_x = index[ 0 ] + ii * incr_orient_trans;
1470 if( index_x < 0 || index_x >
mp_nbVox[ 0 ] - 1 )
continue;
1472 ap_ProjectionLine->
AddVoxel(a_direction, index_x + j *
mp_nbVox[ 0 ] + index_z * m_nbVoxXY, final_weight * distance_z[ jj ] * distance_x[ ii ] * tof_weight);
1491 Cerr(
"***** iProjectorDistanceDriven::ProjectTOFHistogram() -> Called while not initialized !" << endl);
1496 #ifdef CASTOR_VERBOSE 1499 string direction =
"";
1500 if (a_direction==
FORWARD) direction =
"forward";
1501 else direction =
"backward";
1502 Cout(
"iProjectorDistanceDriven::Project with TOF bins -> Project line '" << ap_ProjectionLine <<
"' in " << direction <<
" direction" << endl);
1548 HPFLTNB length_LOR_half = length_LOR * 0.5;
1552 INTNB tof_half_nb_bins = tof_nb_bins/2;
1556 HPFLTNB tof_sigma_sqrt2 = sqrt(2.)*tof_sigma;
1562 HPFLTNB prev_erf = 0., new_erf = 0.;
1565 INTNB tof_bin_last = tof_half_nb_bins;
1566 INTNB tof_bin_first = -tof_half_nb_bins;
1574 INTNB tof_bin_first_for_voxel = 0, tof_bin_last_for_voxel = 0;
1590 &p1_x[1], &p1_y[1], &p1_z[1],
1591 &p2_x[1], &p2_y[1], &p2_z[1]
1594 Cerr(
"***** iProjectorDistanceDriven::ProjectWithoutTOF() -> A problem occurred while getting the edges' center positions !" << endl);
1605 p1_x[ 0 ] = p1[ 0 ]; p1_y[ 0 ] = p1[ 1 ]; p1_z[ 0 ] = p1[ 2 ];
1608 p2_x[ 0 ] = p2[ 0 ]; p2_y[ 0 ] = p2[ 1 ]; p2_z[ 0 ] = p2[ 2 ];
1617 for(
int p = 0; p < 5; ++p )
1619 r[ p ][ 0 ] = p2_x[ p ] - p1_x[ p ];
1620 r[ p ][ 1 ] = p2_y[ p ] - p1_y[ p ];
1621 r[ p ][ 2 ] = p2_z[ p ] - p1_z[ p ];
1626 r[ 0 ][ 0 ] * r[ 0 ][ 0 ],
1627 r[ 0 ][ 1 ] * r[ 0 ][ 1 ],
1628 r[ 0 ][ 2 ] * r[ 0 ][ 2 ]
1633 HPFLTNB alpha_min = 0.0, alpha_max = 1.0;
1637 HPFLTNB alpha_x_min[ 4 ], alpha_x_max[ 4 ];
1638 HPFLTNB alpha_y_min[ 4 ], alpha_y_max[ 4 ];
1639 HPFLTNB alpha_z_min[ 4 ], alpha_z_max[ 4 ];
1664 float const angle = std::acos( ( r[ 0 ][ 0 ] ) / ( std::sqrt( r2[ 0 ] + r2[ 1 ] ) ) ) * 180.0f / M_PI;
1667 if( ( angle >= 0.0 && angle <= 45.0 ) || ( angle >= 135.0 && angle <= 180.0 ) )
1671 alpha_x_min[ 0 ] = ( (-
mp_halfFOV[ 0 ]) - p1_x[ 0 ] ) / r[ 0 ][ 0 ];
1672 alpha_x_min[ 1 ] = (
mp_halfFOV[ 0 ] - p1_x[ 0 ] ) / r[ 0 ][ 0 ];
1673 alpha_x_max[ 0 ] = alpha_x_min[ 0 ];
1674 alpha_x_max[ 1 ] = alpha_x_min[ 1 ];
1675 alpha_min = std::max( alpha_min, std::min( alpha_x_min[ 0 ], alpha_x_min[ 1 ] ) );
1676 alpha_max = std::min( alpha_max, std::max( alpha_x_max[ 0 ], alpha_x_max[ 1 ] ) );
1680 if( r[ 1 ][ 1 ] != 0.0 )
1682 alpha_y_min[ 0 ] = ( (-
mp_halfFOV[ 1 ]) - p1_y[ 1 ] ) / r[ 1 ][ 1 ];
1683 alpha_y_min[ 1 ] = (
mp_halfFOV[ 1 ] - p1_y[ 1 ] ) / r[ 1 ][ 1 ];
1684 alpha_y_max[ 0 ] = alpha_y_min[ 0 ];
1685 alpha_y_max[ 1 ] = alpha_y_min[ 1 ];
1689 alpha_y_min[ 0 ] = 0.0;
1690 alpha_y_min[ 1 ] = 0.0;
1691 alpha_y_max[ 0 ] = 1.0;
1692 alpha_y_max[ 1 ] = 1.0;
1695 if( r[ 2 ][ 1 ] != 0.0 )
1697 alpha_y_min[ 2 ] = ( (-
mp_halfFOV[ 1 ]) - p1_y[ 2 ] ) / r[ 2 ][ 1 ];
1698 alpha_y_min[ 3 ] = (
mp_halfFOV[ 1 ] - p1_y[ 2 ] ) / r[ 2 ][ 1 ];
1699 alpha_y_max[ 2 ] = alpha_y_min[ 2 ];
1700 alpha_y_max[ 3 ] = alpha_y_min[ 3 ];
1704 alpha_y_min[ 2 ] = 0.0;
1705 alpha_y_min[ 3 ] = 0.0;
1706 alpha_y_max[ 2 ] = 1.0;
1707 alpha_y_max[ 3 ] = 1.0;
1711 alpha_min = std::max( alpha_min,
1712 alpha_y_min[ std::minmax_element( alpha_y_min, alpha_y_min + 4 ).first - alpha_y_min ] );
1714 alpha_max = std::min( alpha_max,
1715 alpha_y_max[ std::minmax_element( alpha_y_max, alpha_y_max + 4 ).second - alpha_y_max ] );
1718 if( r[ 3 ][ 2 ] != 0.0 )
1720 alpha_z_min[ 0 ] = ( (-
mp_halfFOV[ 2 ]) - p1_z[ 3 ] ) / r[ 3 ][ 2 ];
1721 alpha_z_min[ 1 ] = (
mp_halfFOV[ 2 ] - p1_z[ 3 ] ) / r[ 3 ][ 2 ];
1722 alpha_z_max[ 0 ] = alpha_z_min[ 0 ];
1723 alpha_z_max[ 1 ] = alpha_z_min[ 1 ];
1727 alpha_z_min[ 0 ] = 0.0;
1728 alpha_z_min[ 1 ] = 0.0;
1729 alpha_z_max[ 0 ] = 1.0;
1730 alpha_z_max[ 1 ] = 1.0;
1733 if( r[ 4 ][ 2 ] != 0.0 )
1735 alpha_z_min[ 2 ] = ( (-
mp_halfFOV[ 2 ]) - p1_z[ 4 ] ) / r[ 4 ][ 2 ];
1736 alpha_z_min[ 3 ] = (
mp_halfFOV[ 2 ] - p1_z[ 4 ] ) / r[ 4 ][ 2 ];
1737 alpha_z_max[ 2 ] = alpha_z_min[ 2 ];
1738 alpha_z_max[ 3 ] = alpha_z_min[ 3 ];
1742 alpha_z_min[ 2 ] = 0.0;
1743 alpha_z_min[ 3 ] = 0.0;
1744 alpha_z_max[ 2 ] = 1.0;
1745 alpha_z_max[ 3 ] = 1.0;
1749 alpha_min = std::max( alpha_min,
1750 alpha_z_min[ std::minmax_element( alpha_z_min, alpha_z_min + 4 ).first
1753 alpha_max = std::min( alpha_max,
1754 alpha_z_max[ std::minmax_element( alpha_z_max, alpha_z_max + 4 ).second
1757 if( alpha_max <= alpha_min )
return 0;
1762 for (
INTNB tof_bin=0; tof_bin<tof_nb_bins; tof_bin++) tof_weights_temp[tof_bin] = 0.;
1765 int16_t i_min = 0, i_max = 0;
1766 if( r[ 0 ][ 0 ] > 0.0 )
1771 else if( r[ 0 ][ 0 ] < 0.0 )
1779 HPFLTNB const factor( ::fabs( r[ 0 ][ 0 ] ) / ::sqrt( r2[ 0 ] + r2[ 1 ] + r2[ 2 ] ) );
1780 HPFLTNB const factor_for_tof( length_LOR / r[ 0 ][ 0 ] );
1785 for( uint8_t p = 0; p < 4; ++p )
1788 ri[ p ][ 1 ] = r[ p + 1 ][ 1 ] / r[ p + 1 ][ 0 ];
1789 ri[ p ][ 2 ] = r[ p + 1 ][ 2 ] / r[ p + 1 ][ 0 ];
1793 int8_t incr_orient_trans = 0;
1794 int8_t idx_orient_trans = 0;
1797 int8_t
const incr_orient_axial = p2_z[ 3 ] < p2_z[ 4 ] ? 1 : -1;
1801 int8_t
const idx_orient_axial = p2_z[ 3 ] < p2_z[ 4 ] ? 1 : 0;
1807 for( int16_t i = i_min; i < i_max; ++i )
1812 pos[ 0 ] = p1_y[ 1 ] + step * ri[ 0 ][ 1 ];
1813 pos[ 1 ] = p1_y[ 2 ] + step * ri[ 1 ][ 1 ];
1815 pos[ 2 ] = p1_z[ 3 ] + step * ri[ 2 ][ 2 ];
1816 pos[ 3 ] = p1_z[ 4 ] + step * ri[ 3 ][ 2 ];
1819 w1 = ::fabs( pos[ 0 ] - pos[ 1 ] );
1820 w2 = ::fabs( pos[ 2 ] - pos[ 3 ] );
1822 if( ( w1 * w2 ) != 0 ) final_weight = weight / ( w1 * w2 );
1835 bool test_index = (index[ 0 ] >= 0 && index[ 2 ] >= 0 && index[ 1 ] >= 0 && index[ 3 ] >= 0
1839 if( index[ 0 ] < index[ 1 ] )
1841 incr_orient_trans = 1;
1842 idx_orient_trans = 1;
1844 else if( index[ 0 ] > index[ 1 ] )
1846 incr_orient_trans = -1;
1847 idx_orient_trans = 0;
1851 int16_t
const n_distance_y = ::abs( index[ 0 ] - index[ 1 ] ) + 2;
1853 distance_y[ 0 ] = pos[ 0 ];
1854 distance_y[ n_distance_y - 1 ] = pos[ 1 ];
1856 if( n_distance_y > 2 )
1858 for( int16_t d = 0; d < n_distance_y - 2; ++d )
1859 distance_y[ d + 1 ] = (-
mp_halfFOV[ 1 ]) + ( index[ 0 ] + idx_orient_trans + ( incr_orient_trans * d ) ) * mp_sizeVox[ 1 ];
1863 for( int16_t d = 0; d < n_distance_y - 1; ++d )
1864 distance_y[ d ] = ::fabs( distance_y[ d + 1 ] - distance_y[ d ] );
1867 int16_t
const n_distance_z = ::abs( index[ 2 ] - index[ 3 ] ) + 2;
1869 distance_z[ 0 ] = pos[ 2 ];
1870 distance_z[ n_distance_z - 1 ] = pos[ 3 ];
1872 if( n_distance_z > 2 )
1874 for( int16_t d = 0; d < n_distance_z - 2; ++d )
1875 distance_z[ d + 1 ] = (-
mp_halfFOV[ 2 ]) + ( index[ 2 ] + idx_orient_axial + ( incr_orient_axial * d ) ) * mp_sizeVox[ 2 ];
1879 for( int16_t d = 0; d < n_distance_z - 1; ++d )
1880 distance_z[ d ] = ::fabs( distance_z[ d + 1 ] - distance_z[ d ] );
1892 tof_bin_first_for_voxel = max(tof_bin_first , (
INTNB)(ceil( (step * factor_for_tof - tof_half_span - length_LOR_half ) /
m_TOFBinSizeInMm )));
1893 tof_bin_last_for_voxel = min(tof_bin_last , (
INTNB)(floor( (step * factor_for_tof + tof_half_span - length_LOR_half) /
m_TOFBinSizeInMm )));
1896 lor_tof_center = length_LOR_half + tof_bin_first_for_voxel *
m_TOFBinSizeInMm;
1899 tof_bin_first_for_voxel += tof_half_nb_bins;
1900 tof_bin_last_for_voxel += tof_half_nb_bins;
1906 HPFLTNB temp = std::min(tof_half_span,std::max(-tof_half_span, lor_tof_center - m_TOFBinSizeInMm/2. - step * factor_for_tof));
1907 prev_erf = erf(temp/tof_sigma_sqrt2);
1912 for (
INTNB tof_bin=tof_bin_first_for_voxel; tof_bin<=tof_bin_last_for_voxel; tof_bin++)
1932 HPFLTNB temp = std::min(tof_half_span,std::max(-tof_half_span, lor_tof_center + m_TOFBinSizeInMm/2. - step * factor_for_tof));
1933 new_erf = erf(temp/tof_sigma_sqrt2);
1934 tof_weights_temp[tof_bin] = 0.5 * fabs(new_erf - prev_erf);
1944 HPFLTNB temp = (step * factor_for_tof - lor_tof_center) / tof_sigma;
1963 int16_t index_y, index_z;
1964 for( int16_t jj = 0; jj < n_distance_z - 1; ++jj )
1966 index_z = index[ 2 ] + jj * incr_orient_axial;
1967 if( index_z < 0 || index_z >
mp_nbVox[ 2 ] - 1 )
continue;
1969 for( int16_t ii = 0; ii < n_distance_y - 1; ++ii )
1971 index_y = index[ 0 ] + ii * incr_orient_trans;
1972 if( index_y < 0 || index_y >
mp_nbVox[ 1 ] - 1 )
continue;
1974 ap_ProjectionLine->
AddVoxelAllTOFBins(a_direction, i + index_y *
mp_nbVox[ 0 ] + index_z * m_nbVoxXY, final_weight * distance_z[ jj ] * distance_y[ ii ], tof_weights_temp, tof_bin_first_for_voxel, tof_bin_last_for_voxel);
1978 delete[] tof_weights_temp;
1984 if( r[ 1 ][ 0 ] != 0.0 )
1986 alpha_x_min[ 0 ] = ( (-
mp_halfFOV[ 0 ]) - p1_x[ 1 ] ) / r[ 1 ][ 0 ];
1987 alpha_x_min[ 1 ] = (
mp_halfFOV[ 0 ] - p1_x[ 1 ] ) / r[ 1 ][ 0 ];
1988 alpha_x_max[ 0 ] = alpha_x_min[ 0 ];
1989 alpha_x_max[ 1 ] = alpha_x_min[ 1 ];
1993 alpha_x_min[ 0 ] = 0.0;
1994 alpha_x_min[ 1 ] = 0.0;
1995 alpha_x_max[ 0 ] = 1.0;
1996 alpha_x_max[ 1 ] = 1.0;
1999 if( r[ 2 ][ 0 ] != 0.0 )
2001 alpha_x_min[ 2 ] = ( (-
mp_halfFOV[ 0 ]) - p1_x[ 2 ] ) / r[ 2 ][ 0 ];
2002 alpha_x_min[ 3 ] = (
mp_halfFOV[ 0 ] - p1_x[ 2 ] ) / r[ 2 ][ 0 ];
2003 alpha_x_max[ 2 ] = alpha_x_min[ 2 ];
2004 alpha_x_max[ 3 ] = alpha_x_min[ 3 ];
2008 alpha_x_min[ 2 ] = 0.0;
2009 alpha_x_min[ 3 ] = 0.0;
2010 alpha_x_max[ 2 ] = 1.0;
2011 alpha_x_max[ 3 ] = 1.0;
2015 alpha_min = std::max( alpha_min,
2016 alpha_x_min[ std::minmax_element( alpha_x_min, alpha_x_min + 4 ).first - alpha_x_min ] );
2018 alpha_max = std::min( alpha_max,
2019 alpha_x_max[ std::minmax_element( alpha_x_max, alpha_x_max + 4 ).second - alpha_x_max ] );
2023 alpha_y_min[ 0 ] = ( (-
mp_halfFOV[ 1 ]) - p1_y[ 0 ] ) / r[ 0 ][ 1 ];
2024 alpha_y_min[ 1 ] = (
mp_halfFOV[ 1 ] - p1_y[ 0 ] ) / r[ 0 ][ 1 ];
2025 alpha_y_max[ 0 ] = alpha_y_min[ 0 ];
2026 alpha_y_max[ 1 ] = alpha_y_min[ 1 ];
2027 alpha_min = std::max( alpha_min,
2028 std::min( alpha_y_min[ 0 ], alpha_y_min[ 1 ] ) );
2029 alpha_max = std::min( alpha_max,
2030 std::max( alpha_y_max[ 0 ], alpha_y_max[ 1 ] ) );
2033 if( r[ 3 ][ 2 ] != 0.0 )
2035 alpha_z_min[ 0 ] = ( (-
mp_halfFOV[ 2 ]) - p1_z[ 3 ] ) / r[ 3 ][ 2 ];
2036 alpha_z_min[ 1 ] = (
mp_halfFOV[ 2 ] - p1_z[ 3 ] ) / r[ 3 ][ 2 ];
2037 alpha_z_max[ 0 ] = alpha_z_min[ 0 ];
2038 alpha_z_max[ 1 ] = alpha_z_min[ 1 ];
2042 alpha_z_min[ 0 ] = 0.0;
2043 alpha_z_min[ 1 ] = 0.0;
2044 alpha_z_max[ 0 ] = 1.0;
2045 alpha_z_max[ 1 ] = 1.0;
2048 if( r[ 4 ][ 2 ] != 0.0 )
2050 alpha_z_min[ 2 ] = ( (-
mp_halfFOV[ 2 ]) - p1_z[ 4 ] ) / r[ 4 ][ 2 ];
2051 alpha_z_min[ 3 ] = (
mp_halfFOV[ 2 ] - p1_z[ 4 ] ) / r[ 4 ][ 2 ];
2052 alpha_z_max[ 2 ] = alpha_z_min[ 2 ];
2053 alpha_z_max[ 3 ] = alpha_z_min[ 3 ];
2057 alpha_z_min[ 2 ] = 0.0;
2058 alpha_z_min[ 3 ] = 0.0;
2059 alpha_z_max[ 2 ] = 1.0;
2060 alpha_z_max[ 3 ] = 1.0;
2064 alpha_min = std::max( alpha_min, alpha_z_min[ std::minmax_element( alpha_z_min, alpha_z_min + 4 ).first - alpha_z_min ] );
2066 alpha_max = std::min( alpha_max, alpha_z_max[ std::minmax_element( alpha_z_max, alpha_z_max + 4 ).second - alpha_z_max ] );
2068 if( alpha_max <= alpha_min )
return 0;
2073 for (
INTNB tof_bin=0; tof_bin<tof_nb_bins; tof_bin++) tof_weights_temp[tof_bin] = 0.;
2076 int16_t j_min = 0, j_max = 0;
2077 if( r[ 0 ][ 1 ] > 0.0 )
2082 else if( r[ 0 ][ 1 ] < 0.0 )
2090 HPFLTNB const factor( ::fabs( r[ 0 ][ 1 ] ) / ::sqrt( r2[ 0 ] + r2[ 1 ] + r2[ 2 ] ) );
2091 HPFLTNB const factor_for_tof( length_LOR / r[ 0 ][ 1 ] );
2096 for( uint8_t p = 0; p < 4; ++p )
2098 ri[ p ][ 0 ] = r[ p + 1 ][ 0 ] / r[ p + 1 ][ 1 ];
2100 ri[ p ][ 2 ] = r[ p + 1 ][ 2 ] / r[ p + 1 ][ 1 ];
2104 int8_t incr_orient_trans = 0;
2105 int8_t idx_orient_trans = 0;
2108 int8_t
const incr_orient_axial = p2_z[ 3 ] < p2_z[ 4 ] ? 1 : -1;
2110 int8_t
const idx_orient_axial = p2_z[ 3 ] < p2_z[ 4 ] ? 1 : 0;
2116 for( int16_t j = j_min; j < j_max; ++j )
2121 pos[ 0 ] = p1_x[ 1 ] + step * ri[ 0 ][ 0 ];
2122 pos[ 1 ] = p1_x[ 2 ] + step * ri[ 1 ][ 0 ];
2124 pos[ 2 ] = p1_z[ 3 ] + step * ri[ 2 ][ 2 ];
2125 pos[ 3 ] = p1_z[ 4 ] + step * ri[ 3 ][ 2 ];
2128 w1 = ::fabs( pos[ 0 ] - pos[ 1 ] );
2129 w2 = ::fabs( pos[ 2 ] - pos[ 3 ] );
2131 if( ( w1 * w2 ) != 0 ) final_weight = weight / ( w1 * w2 );
2143 bool test_index = (index[ 0 ] >= 0 && index[ 2 ] >= 0 && index[ 1 ] >= 0 && index[ 3 ] >= 0
2147 if( index[ 0 ] < index[ 1 ] )
2149 incr_orient_trans = 1;
2150 idx_orient_trans = 1;
2152 else if( index[ 0 ] > index[ 1 ] )
2154 incr_orient_trans = -1;
2155 idx_orient_trans = 0;
2159 int16_t
const n_distance_x = ::abs( index[ 0 ] - index[ 1 ] ) + 2;
2161 distance_x[ 0 ] = pos[ 0 ];
2162 distance_x[ n_distance_x - 1 ] = pos[ 1 ];
2164 if( n_distance_x > 2 )
2166 for( int16_t d = 0; d < n_distance_x - 2; ++d )
2167 distance_x[ d + 1 ] = (-
mp_halfFOV[ 0 ]) + ( index[ 0 ] + idx_orient_trans + ( incr_orient_trans * d ) ) * mp_sizeVox[ 0 ];
2171 for( int16_t d = 0; d < n_distance_x - 1; ++d )
2172 distance_x[ d ] = ::fabs( distance_x[ d + 1 ] - distance_x[ d ] );
2175 int16_t
const n_distance_z = ::abs( index[ 2 ] - index[ 3 ] ) + 2;
2177 distance_z[ 0 ] = pos[ 2 ];
2178 distance_z[ n_distance_z - 1 ] = pos[ 3 ];
2180 if( n_distance_z > 2 )
2182 for( int16_t d = 0; d < n_distance_z - 2; ++d )
2183 distance_z[ d + 1 ] = (-
mp_halfFOV[ 2 ]) + ( index[ 2 ] + idx_orient_axial + ( incr_orient_axial * d ) ) * mp_sizeVox[ 2 ];
2187 for( int16_t d = 0; d < n_distance_z - 1; ++d )
2188 distance_z[ d ] = ::fabs( distance_z[ d + 1 ] - distance_z[ d ] );
2200 tof_bin_first_for_voxel = max(tof_bin_first , (
INTNB)(ceil( (step * factor_for_tof - tof_half_span - length_LOR_half ) /
m_TOFBinSizeInMm )));
2201 tof_bin_last_for_voxel = min(tof_bin_last , (
INTNB)(floor( (step * factor_for_tof + tof_half_span - length_LOR_half) /
m_TOFBinSizeInMm )));
2204 lor_tof_center = length_LOR_half + tof_bin_first_for_voxel *
m_TOFBinSizeInMm;
2207 tof_bin_first_for_voxel += tof_half_nb_bins;
2208 tof_bin_last_for_voxel += tof_half_nb_bins;
2214 HPFLTNB temp = std::min(tof_half_span,std::max(-tof_half_span, lor_tof_center - m_TOFBinSizeInMm/2. - step * factor_for_tof));
2215 prev_erf = erf(temp/tof_sigma_sqrt2);
2220 for (
INTNB tof_bin=tof_bin_first_for_voxel; tof_bin<=tof_bin_last_for_voxel; tof_bin++)
2240 HPFLTNB temp = std::min(tof_half_span,std::max(-tof_half_span, lor_tof_center + m_TOFBinSizeInMm/2. - step * factor_for_tof));
2241 new_erf = erf(temp/tof_sigma_sqrt2);
2242 tof_weights_temp[tof_bin] = 0.5 * fabs(new_erf - prev_erf);
2252 HPFLTNB temp = (step * factor_for_tof - lor_tof_center) / tof_sigma;
2271 int16_t index_x, index_z;
2272 for( int16_t jj = 0; jj < n_distance_z - 1; ++jj )
2274 index_z = index[ 2 ] + jj * incr_orient_axial;
2275 if( index_z < 0 || index_z >
mp_nbVox[ 2 ] - 1 )
continue;
2277 for( int16_t ii = 0; ii < n_distance_x - 1; ++ii )
2279 index_x = index[ 0 ] + ii * incr_orient_trans;
2280 if( index_x < 0 || index_x >
mp_nbVox[ 0 ] - 1 )
continue;
2282 ap_ProjectionLine->
AddVoxelAllTOFBins(a_direction, index_x + j *
mp_nbVox[ 0 ] + index_z * m_nbVoxXY, final_weight * distance_z[ jj ] * distance_x[ ii ], tof_weights_temp, tof_bin_first_for_voxel, tof_bin_last_for_voxel);
2287 delete[] tof_weights_temp;
bool m_compatibleWithSPECTAttenuationCorrection
#define SPEED_OF_LIGHT_IN_MM_PER_PS
FLTNB m_TOFGaussianNormCoef
iProjectorDistanceDriven()
The constructor of iProjectorDistanceDriven.
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.
INTNB EstimateMaxNumberOfVoxelsPerLine()
This function is used to compute and provide an estimate of the maximum number of voxels that could c...
#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 ReadOptionsList(const string &a_optionsList)
A function used to read options from a list of options.
int ProjectWithoutTOF(int a_direction, oProjectionLine *ap_ProjectionLine)
A function to project without TOF.
bool m_TOFWeightingFcnPrecomputedFlag
bool m_TOFBinProperProcessingFlag
HPFLTNB * mp_TOFWeightingFcn
FLTNB GetLength()
This function is used to get the length of the line.
int CheckSpecificParameters()
A private function used to check the parameters settings specific to the child projector.
int GetIndex2()
This function is used to get the index associated to point 2.
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...
~iProjectorDistanceDriven()
The destructor of iProjectorDistanceDriven.
int ReadConfigurationFile(const string &a_configurationFile)
A function used to read options from a configuration file.
FLTNB m_TOFResolutionInMm
int InitializeSpecific()
This function is used to initialize specific stuff to the child projector.
INTNB GetNbVoxXY()
Get the number of voxels in a slice.
Declaration of class iProjectorDistanceDriven.
This class is designed to manage and store system matrix elements associated to a vEvent...
FLTNB m_TOFPrecomputedSamplingFactor
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).
bool m_compatibleWithCompression
int GetIndex1()
This function is used to get the index associated to point 1.
int ProjectTOFListmode(int a_direction, oProjectionLine *ap_ProjectionLine)
A function to project with TOF continuous information.
void ShowHelpSpecific()
A function used to show help about the child projector.
INTNB m_TOFWeightingFcnNbSamples
virtual int GetEdgesCenterPositions(int a_index1, int a_index2, FLTNB ap_pos_line_point1[3], FLTNB ap_pos_line_point2[3], FLTNB ap_pos_point1_x[4], FLTNB ap_pos_point1_y[4], FLTNB ap_pos_point1_z[4], FLTNB ap_pos_point2_x[4], FLTNB ap_pos_point2_y[4], FLTNB ap_pos_point2_z[4])=0
This is a pure virtual method that must be implemented by children. Get the cartesian coordinaters ...