shifting the first detector index to the first crystal in the firts bucket

Author: danieldeidda

src/buildblock/GeometryBlocksOnCylindrical.cxx

@@ -169,12 +169,22 @@ build_crystal_maps()
 		see start_x*/
 
 		//calculate start_point to build the map.
+
+//    estimate the angle covered by half bucket, csi
+    float csi=_PI/num_transaxial_buckets;
+    float trans_blocks_gap=transaxial_block_spacing-num_transaxial_crystals_per_block*transaxial_crystal_spacing;
+    float csi_minus_csiGaps=csi-(csi/transaxial_block_spacing*2)*
+            (transaxial_crystal_spacing/2+trans_blocks_gap);
+//    distance between the center of the scannner and the first crystal in the bucket, r=Reffective/cos(csi)
+    float r=get_scanner_ptr()->get_effective_ring_radius()/cos(csi_minus_csiGaps);
+    
 	float start_z = 0;
 	float start_y = -1*get_scanner_ptr()->get_effective_ring_radius();
-	float start_x = -1*(
-								((num_transaxial_blocks_per_bucket-1)/2.)*transaxial_block_spacing
-							+ ((num_transaxial_crystals_per_block-1)/2.)*transaxial_crystal_spacing
-											); //the first crystal in the bucket
+	float start_x = -1*r*sin(csi_minus_csiGaps);//(
+//								((num_transaxial_blocks_per_bucket-1)/2.)*transaxial_block_spacing
+//							+ ((num_transaxial_crystals_per_block-1)/2.)*transaxial_crystal_spacing
+//											); //the first crystal in the bucket
+    
 	stir::CartesianCoordinate3D<float> start_point(start_z, start_y, start_x);
 
 	for (int ax_block_num=0; ax_block_num<num_axial_blocks; ++ax_block_num)
@@ -184,20 +194,12 @@ build_crystal_maps()
 					for (int trans_crys_num=0; trans_crys_num<num_transaxial_crystals_per_block; ++trans_crys_num)
 	{
 		// calculate detection position for a given detector
-		// note: in STIR convention, crystal(0,0,0) corresponds to card_coord(z=-l/2,y=-r,x=0)
+		// note: in STIR convention, crystal(0,0,0) corresponds to card_coord(z=0,y=-r,x=0)
 		int tangential_coord;
-		if (num_transaxial_blocks_per_bucket%2==0)
 			tangential_coord = trans_bucket_num*num_transaxial_blocks_per_bucket*num_transaxial_crystals_per_block
 														+ trans_block_num*num_transaxial_crystals_per_block
-														+ trans_crys_num
-														- num_transaxial_blocks_per_bucket/2*num_transaxial_crystals_per_block;
-		else
-			tangential_coord = trans_bucket_num*num_transaxial_blocks_per_bucket*num_transaxial_crystals_per_block
-														+ trans_block_num*num_transaxial_crystals_per_block
-														+ trans_crys_num
-														- num_transaxial_blocks_per_bucket/2*num_transaxial_crystals_per_block
-														- num_transaxial_crystals_per_block/2;
-
+														+ trans_crys_num;
+		
 		if (tangential_coord<0)
 					tangential_coord += num_detectors_per_ring;
 
@@ -210,7 +212,8 @@ build_crystal_maps()
 									ax_block_num*axial_block_spacing + ax_crys_num*axial_crystal_spacing,
 									0.,
 									trans_block_num*transaxial_block_spacing + trans_crys_num*transaxial_crystal_spacing);
-		float alpha = trans_bucket_num*(2*_PI)/num_transaxial_buckets;
+		float alpha = get_scanner_ptr()->get_intrinsic_azimuthal_tilt()+
+                trans_bucket_num*(2*_PI)/num_transaxial_buckets+csi_minus_csiGaps;
 
 		stir::Array<2, float> rotation_matrix = get_rotation_matrix(alpha);
 		// to match index range of CartesianCoordinate3D, which is 1 to 3

src/buildblock/ProjDataInfoGeneric.cxx

@@ -499,7 +499,7 @@ get_LOR(LORInAxialAndNoArcCorrSinogramCoordinates<float>& lor,
 	find_cartesian_coordinates_of_detection(_p1, _p2, bin);
 
 	LORAs2Points<float> lor_as_2_points(_p1, _p2);
-	const double R = get_ring_radius();
+	const double R = sqrt(max(square(_p1.x())+square(_p1.y()), square(_p2.x())+square(_p2.y())));
 
     lor_as_2_points.change_representation(lor, R);
 }

src/recon_test/test_blocks_on_cylindrical_projectors.cxx

@@ -106,13 +106,22 @@ BlocksTests::run_plane_symmetry_test(){
 
 //    rotate by 30 degrees
     phi2=30*_PI/180;
-    VoxelsOnCartesianGrid<float>  image2=image;
+    VoxelsOnCartesianGrid<float>  image2=*image.get_empty_copy();
     const Array<2,float> direction2=make_array(make_1d_array(1.F,0.F,0.F),
     make_1d_array(0.F,cos(float(_PI)-phi2),sin(float(_PI)-phi2)),
     make_1d_array(0.F,-sin(float(_PI)-phi2),cos(float(_PI)-phi2)));
-    plane.set_direction_vectors(direction2);
+    
+    Ellipsoid
+      plane2(CartesianCoordinate3D<float>(/*edge_z*/25*grid_spacing.z(),
+                                             /*edge_y*/91*grid_spacing.y(),
+                                             /*edge_x*/5*grid_spacing.x()),
+		        /*centre*/CartesianCoordinate3D<float>((image.get_min_index()+image.get_max_index())/2*grid_spacing.z(),
+                                                       0*grid_spacing.y(),
+                                                       0),
+                direction2);
+//    plane.set_direction_vectors(direction2);
 
-    plane.construct_volume(image2, make_coordinate(3,3,3));
+    plane2.construct_volume(image2, make_coordinate(3,3,3));
 
 //    create projadata info
 
@@ -196,20 +205,21 @@ BlocksTests::run_plane_symmetry_test(){
     forw_projector2_sptr->forward_project(*projdata2, *image2_sptr);
 
     int view1_num = 0, view2_num = 0;
-    LORInAxialAndNoArcCorrSinogramCoordinates<float> lorB;
+    LORInAxialAndNoArcCorrSinogramCoordinates<float> lorB1;
     for(int i=0;i<projdata->get_max_view_num();i++){
         Bin bin(0,i,0,0);
-        proj_data_info_blocks_ptr->get_LOR(lorB,bin);
-        if(abs(lorB.phi()-phi1)/phi1<=1E-2){
+        proj_data_info_blocks_ptr->get_LOR(lorB1,bin);
+        if(abs(lorB1.phi()-phi1)/phi1<=1E-2){
             view1_num=i;
             break;
         }
     }
 
+    LORInAxialAndNoArcCorrSinogramCoordinates<float> lorB2;
     for(int i=0;i<projdata2->get_max_view_num();i++){
         Bin bin(0,i,0,0);
-        proj_data_info_blocks_ptr->get_LOR(lorB,bin);
-        if(abs(lorB.phi()-phi2)/phi2<=1E-2){
+        proj_data_info_blocks_ptr->get_LOR(lorB2,bin);
+        if(abs(lorB2.phi()-phi2)/phi2<=1E-2){
             view2_num=i;
             break;
         }
@@ -231,16 +241,16 @@ BlocksTests::run_plane_symmetry_test(){
         for(int tang=projdata2->get_min_tangential_pos_num();tang<projdata2->get_max_tangential_pos_num();tang++){
 
             if((max2-projdata2->get_sinogram(0,0).at(view2_num).at(tang))/max2<1E-3) {
-                tang1_num=tang;
+                tang2_num=tang;
                 break;
             }
         }
 
         float bin1=projdata->get_sinogram(0,0).at(view1_num).at(tang1_num);
         float bin2=projdata2->get_sinogram(0,0).at(view2_num).at(tang2_num);
-    set_tolerance(10E-3);
-    check_if_equal(bin1, max1,"the value seen in the block at 30 degrees should be the same as the max value of the sinogram");
-    check_if_equal(bin2, max2,"the value seen in the block at 60 degrees should be the same as the max value of the sinogram");
+    set_tolerance(10E-2);
+    check_if_equal(bin1, max1,"the value seen in the block at 60 degrees should be the same as the max value of the sinogram");
+    check_if_equal(bin2, max2,"the value seen in the block at 30 degrees should be the same as the max value of the sinogram");
 }
 
 /*! The following is a test for symmetries: a simulated image is created with spherical source in front of each detector block,

src/test/test_proj_data_info.cxx

@@ -489,7 +489,7 @@ ProjDataInfoTests::run_Blocks_DOI_test()
     int Bring1, Bring2, Bdet1,Bdet2, BDring1, BDring2, BDdet1, BDdet2;
     CartesianCoordinate3D< float> b1,b2,bd1,bd2;
     float doi=scannerBlocksDOI_ptr->get_average_depth_of_interaction();
-//    timer.reset(); timer.start();
+    timer.reset(); timer.start();
 
     for (int seg =proj_data_info_blocks_doi0_ptr->get_min_segment_num(); seg <=proj_data_info_blocks_doi0_ptr->get_max_segment_num();  ++ seg)
       for (int ax =proj_data_info_blocks_doi0_ptr->get_min_axial_pos_num(seg); ax <=proj_data_info_blocks_doi0_ptr->get_max_axial_pos_num(seg); ++ax)
@@ -545,15 +545,30 @@ ProjDataInfoTests::run_coordinate_test()
     scannerBlocks_ptr->set_num_axial_crystals_per_block(1);
     scannerBlocks_ptr->set_axial_block_spacing(scannerBlocks_ptr->get_axial_crystal_spacing()*
                                                scannerBlocks_ptr->get_num_axial_crystals_per_block());
-    scannerBlocks_ptr->set_transaxial_block_spacing(scannerBlocks_ptr->get_transaxial_crystal_spacing()*
-                                                    scannerBlocks_ptr->get_num_transaxial_crystals_per_block());
     scannerBlocks_ptr->set_num_transaxial_crystals_per_block(1);
     scannerBlocks_ptr->set_num_axial_blocks_per_bucket(1);
     scannerBlocks_ptr->set_num_transaxial_blocks_per_bucket(1);
+        scannerBlocks_ptr->set_transaxial_block_spacing(scannerBlocks_ptr->get_transaxial_crystal_spacing()*
+                                                    scannerBlocks_ptr->get_num_transaxial_crystals_per_block());
     scannerBlocks_ptr->set_num_rings(1);
 
     scannerBlocks_ptr->set_scanner_geometry("BlocksOnCylindrical");
 
+    int num_transaxial_blocks_per_bucket = scannerBlocks_ptr->get_num_transaxial_blocks_per_bucket();
+    int num_transaxial_crystals_per_block = scannerBlocks_ptr->get_num_transaxial_crystals_per_block();
+	int num_transaxial_buckets = scannerBlocks_ptr->get_num_transaxial_blocks()/num_transaxial_blocks_per_bucket;
+    float transaxial_block_spacing = scannerBlocks_ptr->get_transaxial_block_spacing();
+    float transaxial_crystal_spacing = scannerBlocks_ptr->get_transaxial_crystal_spacing();
+//    estimate the angle covered by a bucket, alpha
+    
+    float csi=_PI/num_transaxial_buckets;
+    float trans_blocks_gap=transaxial_block_spacing-num_transaxial_crystals_per_block*transaxial_crystal_spacing;
+    float csi_minus_csiGaps=csi-(csi/transaxial_block_spacing/2)*
+            (transaxial_crystal_spacing/2+trans_blocks_gap);
+    
+    float dx=scannerBlocks_ptr->get_effective_ring_radius()*sin(csi_minus_csiGaps);
+    float dy=scannerBlocks_ptr->get_effective_ring_radius()-scannerBlocks_ptr->get_effective_ring_radius()*cos(csi_minus_csiGaps);
+
     shared_ptr<Scanner> scannerCyl_ptr;
     scannerCyl_ptr.reset(new Scanner (Scanner::SAFIRDualRingPrototype));
     scannerCyl_ptr->set_num_axial_crystals_per_block(1);
@@ -604,7 +619,7 @@ ProjDataInfoTests::run_coordinate_test()
     int Bring1, Bring2, Bdet1,Bdet2, Cring1, Cring2, Cdet1, Cdet2;
     int RTring1, RTring2, RTdet1,RTdet2;
     CartesianCoordinate3D< float> b1,b2,c1,c2,roundt1, roundt2;
-//    timer.reset(); timer.start();
+    timer.reset(); timer.start();
     LORInAxialAndNoArcCorrSinogramCoordinates<float> lorB;
     LORInAxialAndNoArcCorrSinogramCoordinates<float> lorC;
 
@@ -796,7 +811,14 @@ ProjDataInfoTests::run_coordinate_test_for_realistic_scanner()
 
     int Bring1, Bring2, Bdet1,Bdet2, Cring1, Cring2, Cdet1, Cdet2;
     CartesianCoordinate3D< float> b1,b2,c1,c2;
-//    timer.reset(); timer.start();
+    
+    //    estimate the angle covered by half bucket, csi
+        float csi=_PI/scannerBlocks_ptr->get_num_transaxial_buckets();
+    //    distance between the center of the scannner and the first crystal in the bucket, r=Reffective/cos(csi)
+        float r=scannerBlocks_ptr->get_effective_ring_radius()/cos(csi);
+        float max_tolerance=abs(scannerBlocks_ptr->get_effective_ring_radius()-r);
+        
+    timer.reset(); timer.start();
     LORInAxialAndNoArcCorrSinogramCoordinates<float> lorB;
     LORInAxialAndNoArcCorrSinogramCoordinates<float> lorC;
 
@@ -842,14 +864,15 @@ ProjDataInfoTests::run_coordinate_test_for_realistic_scanner()
                 proj_data_info_blocks_ptr->find_cartesian_coordinates_of_detection(b1,b2,bin);
 
                 // we expect to be differences of the order of the mm in x and y due to the difference in geometry
-                set_tolerance(10E-1);
 
-                check_if_equal(b1.z(),c1.z(), " ");
-                check_if_equal(b2.z(),c2.z(), " ");
-                check_if_equal(b1.y(),c1.y(), " ");
-                check_if_equal(b2.y(),c2.y(), " ");
-                check_if_equal(b1.x(),c1.x(), " ");
-                check_if_equal(b2.x(),c2.x(), " ");
+                set_tolerance(max_tolerance);
+                
+                check_if_equal(b1.z(),c1.z(), " checking cartesian coordinate z1");
+                check_if_equal(b2.z(),c2.z(), " checking cartesian coordinate z2");
+                check_if_equal(b1.y(),c1.y(), " checking cartesian coordinate y1");
+                check_if_equal(b2.y(),c2.y(), " checking cartesian coordinate y2");
+                check_if_equal(b1.x(),c1.x(), " checking cartesian coordinate x1");
+                check_if_equal(b2.x(),c2.x(), " checking cartesian coordinate x2");
 
               }
     timer.stop(); std::cerr<< "-- CPU Time " << timer.value() << '\n';
@@ -948,31 +971,30 @@ run_lor_get_s_test(){
 //    num_det_per_ring:2PI=det_id_diff:PI/2
     int det_id_diff=scannerBlocks_ptr->get_num_detectors_per_ring()*(_PI/2)/(2*_PI);
     int Ctb=scannerCyl_ptr->get_num_transaxial_crystals_per_block();
-    float crystal_trans_spacing=scannerBlocks_ptr->get_transaxial_crystal_spacing();
+    float transaxial_crystal_spacing=scannerBlocks_ptr->get_transaxial_crystal_spacing();
     float block_trans_spacing=scannerBlocks_ptr->get_transaxial_block_spacing();
     float prev_s=0;
     float prev_phi=0;
     for (int i=0; i<scannerCyl_ptr->get_num_transaxial_crystals_per_block();i++){
 
-        Cring1=0;Cdet1=i+2*Ctb-Ctb/2;
-        Cring2=0;Cdet2=2*Ctb-Ctb/2+det_id_diff+Ctb-1 -i;
+        Cring1=0;Cdet1=i+2*Ctb;
+        Cring2=0;Cdet2=2*Ctb+det_id_diff+Ctb-1 -i;
 
         proj_data_info_blocks_ptr->get_bin_for_det_pair(bin, Cdet1,Cring1,Cdet2,Cring1);
         proj_data_info_blocks_ptr->get_LOR(lorB,bin);
         float R=block_trans_spacing*(sin(_PI*5/12)+sin(_PI/4)+sin(_PI/12));
         float s=R*cos(_PI/3)+
-                crystal_trans_spacing/2*sin(_PI/4)+
-                (i)*crystal_trans_spacing*sin(_PI/4);
+                transaxial_crystal_spacing/2*sin(_PI/4)+
+                (i)*transaxial_crystal_spacing*sin(_PI/4);
 
-        float s_step=crystal_trans_spacing*sin(_PI/4);
+        float s_step=transaxial_crystal_spacing*sin(_PI/4);
 
 //        the following fails at the moment
 //        check_if_equal(s, lorB.s(),std::to_string(i)+ " Blocks get_s()  is different");
 //        the first value we expect to be different
+        set_tolerance(10E-3);
         if (i>0){
-        check_if_equal(s_step, lorB.s()-prev_s,std::to_string(i)+ " Blocks get_s() the step is different");//+
-//                       std::to_string(crystal_trans_spacing*sin(_PI/4)) + " lorB.s() the step is "+std::to_string(lorB.s()-prev_s)+
-//                       + "phi step is "+std::to_string(lorB.phi()-prev_phi));
+        check_if_equal(s_step, lorB.s()-prev_s,std::to_string(i)+ " Blocks get_s() the step is different");
         check_if_equal(0.F, lorB.phi()-prev_phi, " Blocks get_phi() should be always the same as we are considering parallel LORs");
         }
         prev_s=lorB.s();