diff --git a/src/buildblock/GeometryBlocksOnCylindrical.cxx b/src/buildblock/GeometryBlocksOnCylindrical.cxx index eb13fb2344..2c66384fac 100644 --- a/src/buildblock/GeometryBlocksOnCylindrical.cxx +++ b/src/buildblock/GeometryBlocksOnCylindrical.cxx @@ -35,6 +35,7 @@ limitations under the License. #include "stir/Array.h" #include "stir/make_array.h" #include "stir/numerics/MatrixFunction.h" +#include "stir/warning.h" #include #include #include @@ -63,86 +64,137 @@ GeometryBlocksOnCylindrical::build_crystal_maps(const Scanner& scanner) { // local variables to describe scanner int num_axial_crystals_per_block = scanner.get_num_axial_crystals_per_block(); + int num_axial_buckets = scanner.get_num_axial_buckets(); + int num_axial_blocks_per_bucket = scanner.get_num_axial_blocks_per_bucket(); int num_transaxial_crystals_per_block = scanner.get_num_transaxial_crystals_per_block(); int num_transaxial_blocks_per_bucket = scanner.get_num_transaxial_blocks_per_bucket(); - int num_axial_blocks_per_bucket = scanner.get_num_axial_blocks_per_bucket(); int num_transaxial_buckets = scanner.get_num_transaxial_buckets(); - int num_axial_buckets = scanner.get_num_axial_buckets(); - int num_detectors_per_ring = scanner.get_num_detectors_per_ring(); - float axial_block_spacing = scanner.get_axial_block_spacing(); - float transaxial_block_spacing = scanner.get_transaxial_block_spacing(); - float axial_crystal_spacing = scanner.get_axial_crystal_spacing(); - float transaxial_crystal_spacing = scanner.get_transaxial_crystal_spacing(); det_pos_to_coord_type cartesian_coord_map_given_detection_position_keys; - /*Building starts from a bucket perpendicular to y axis, from its first crystal. - 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 ax_blocks_gap = axial_block_spacing - num_axial_crystals_per_block * axial_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 = scanner.get_effective_ring_radius() / cos(csi_minus_csiGaps); - - float start_z = -(axial_block_spacing * (num_axial_blocks_per_bucket)*num_axial_buckets - axial_crystal_spacing - - ax_blocks_gap * (num_axial_blocks_per_bucket * num_axial_buckets - 1)) - / 2; - float start_y = -1 * scanner.get_effective_ring_radius(); - 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 start_point(start_z, start_y, start_x); + float csi_minus_csiGaps = get_csi_minus_csi_gaps(scanner); + + // calculate first_crystal_offset to build the map + stir::CartesianCoordinate3D first_crystal_offset(get_initial_axial_z_offset(scanner), + -scanner.get_effective_ring_radius(), + get_initial_axial_x_offset_for_each_bucket(scanner)); + int radial_coord = 0; + + // Lopp over axial geometry for (int ax_bucket_num = 0; ax_bucket_num < num_axial_buckets; ++ax_bucket_num) for (int ax_block_num = 0; ax_block_num < num_axial_blocks_per_bucket; ++ax_block_num) for (int ax_crys_num = 0; ax_crys_num < num_axial_crystals_per_block; ++ax_crys_num) - for (int trans_bucket_num = 0; trans_bucket_num < num_transaxial_buckets; ++trans_bucket_num) - for (int trans_block_num = 0; trans_block_num < num_transaxial_blocks_per_bucket; ++trans_block_num) - 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=0,y=-r,x=0) - int tangential_coord; - 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; - - if (tangential_coord < 0) - tangential_coord += num_detectors_per_ring; - - int axial_coord = ax_bucket_num * num_axial_blocks_per_bucket * num_axial_crystals_per_block - + ax_block_num * num_axial_crystals_per_block + ax_crys_num; - int radial_coord = 0; - stir::DetectionPosition<> det_pos(tangential_coord, axial_coord, radial_coord); - - // calculate cartesian coordinate for a given detector - stir::CartesianCoordinate3D transformation_matrix( - 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 = scanner.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 - rotation_matrix.set_min_index(1); - rotation_matrix[1].set_min_index(1); - rotation_matrix[2].set_min_index(1); - rotation_matrix[3].set_min_index(1); - - stir::CartesianCoordinate3D transformed_coord = start_point + transformation_matrix; - stir::CartesianCoordinate3D cart_coord = stir::matrix_multiply(rotation_matrix, transformed_coord); - - cartesian_coord_map_given_detection_position_keys[det_pos] = cart_coord; - } + { + int axial_coord = get_axial_coord(scanner, ax_bucket_num, ax_block_num, ax_crys_num); + float axial_translation = get_axial_translation(scanner, ax_bucket_num, ax_block_num, ax_crys_num); + + // Loop over transaxial geometry + for (int trans_bucket_num = 0; trans_bucket_num < num_transaxial_buckets; ++trans_bucket_num) + for (int trans_block_num = 0; trans_block_num < num_transaxial_blocks_per_bucket; ++trans_block_num) + 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=0,y=-r,x=0) + int transaxial_coord = get_transaxial_coord(scanner, trans_bucket_num, trans_block_num, trans_crys_num); + stir::DetectionPosition<> det_pos(transaxial_coord, axial_coord, radial_coord); + + // The translation matrix from the first crystal in the block + stir::CartesianCoordinate3D translation_matrix( + axial_translation, + 0., + get_crystal_in_bucket_transaxial_translation(scanner, trans_block_num, trans_crys_num)); + + stir::CartesianCoordinate3D transformed_coord = first_crystal_offset + translation_matrix; + + // Calculate the rotation of the crystal + float alpha = scanner.get_intrinsic_azimuthal_tilt() + trans_bucket_num * (2 * _PI) / num_transaxial_buckets + + csi_minus_csiGaps; + cartesian_coord_map_given_detection_position_keys[det_pos] + = calculate_crystal_rotation(transformed_coord, alpha); + } + } set_detector_map(cartesian_coord_map_given_detection_position_keys); } +CartesianCoordinate3D +GeometryBlocksOnCylindrical::calculate_crystal_rotation(const CartesianCoordinate3D& crystal_position, + const float alpha) const +{ + stir::Array<2, float> rotation_matrix = get_rotation_matrix(alpha); + // to match index range of CartesianCoordinate3D, which is 1 to 3 + rotation_matrix.set_min_index(1); + rotation_matrix[1].set_min_index(1); + rotation_matrix[2].set_min_index(1); + rotation_matrix[3].set_min_index(1); + return stir::matrix_multiply(rotation_matrix, crystal_position); +} + +int +GeometryBlocksOnCylindrical::get_transaxial_coord(const Scanner& scanner, + int transaxial_bucket_num, + int transaxial_block_num, + int transaxial_crystal_num) +{ + return transaxial_bucket_num * scanner.get_num_transaxial_blocks_per_bucket() * scanner.get_num_transaxial_crystals_per_block() + + transaxial_block_num * scanner.get_num_transaxial_crystals_per_block() + transaxial_crystal_num; +} + +int +GeometryBlocksOnCylindrical::get_axial_coord(const Scanner& scanner, + int axial_bucket_num, + int axial_block_num, + int axial_crystal_num) +{ + return axial_bucket_num * scanner.get_num_axial_blocks_per_bucket() * scanner.get_num_axial_crystals_per_block() + + axial_block_num * scanner.get_num_axial_crystals_per_block() + axial_crystal_num; +} + +float +GeometryBlocksOnCylindrical::get_crystal_in_bucket_transaxial_translation(const Scanner& scanner, + int transaxial_block_num, + int transaxial_crystal_num) +{ + // Currently, only supports 1 transaxial bucket per angle + return transaxial_block_num * scanner.get_transaxial_block_spacing() + + transaxial_crystal_num * scanner.get_transaxial_crystal_spacing(); +} + +float +GeometryBlocksOnCylindrical::get_axial_translation(const Scanner& scanner, + int axial_bucket_num, + int axial_block_num, + int axial_crystal_num) +{ + return // axial_bucket_num * scanner.get_axial_bucket_spacing() + + axial_block_num * scanner.get_axial_block_spacing() + axial_crystal_num * scanner.get_axial_crystal_spacing(); +} + +float +GeometryBlocksOnCylindrical::get_initial_axial_z_offset(const Scanner& scanner) +{ + // Crystals in a block are centered, blocks in a bucket are centered, and buckets are centered in the z axis. + // This centers the scanner in z + float crystals_in_block_offset = (scanner.get_num_axial_crystals_per_block() - 1) * scanner.get_axial_crystal_spacing(); + float blocks_in_bucket_offset = (scanner.get_num_axial_blocks_per_bucket() - 1) * scanner.get_axial_block_spacing(); + // float bucket_offset = (scanner.get_num_axial_buckets() - 1) * scanner.get_axial_bucket_spacing(); + float bucket_offset = 0; + + // Negative because the scanner is centered at z=0 and increases axial coordinates increase + // 1/2 because it is half the distance from the center to the edge of the scanner + return -(1.0 / 2) * (crystals_in_block_offset + blocks_in_bucket_offset + bucket_offset); +} + +float +GeometryBlocksOnCylindrical::get_initial_axial_x_offset_for_each_bucket(const Scanner& scanner) +{ + // This is the old method... This is probably wrong + // float csi_minus_csiGaps = get_csi_minus_csi_gaps(scanner); + // float r = scanner.get_effective_ring_radius() / cos(csi_minus_csiGaps); + // return -1 * r * sin(csi_minus_csiGaps); + + auto first_crystal_coord = get_crystal_in_bucket_transaxial_translation(scanner, 0, 0); + auto last_crystal_coord = get_crystal_in_bucket_transaxial_translation( + scanner, scanner.get_num_transaxial_blocks_per_bucket() - 1, scanner.get_num_transaxial_crystals_per_block() - 1); + return -(1.0 / 2) * (first_crystal_coord + last_crystal_coord); +} END_NAMESPACE_STIR diff --git a/src/buildblock/Scanner.cxx b/src/buildblock/Scanner.cxx index b66dee692e..68394cc37d 100644 --- a/src/buildblock/Scanner.cxx +++ b/src/buildblock/Scanner.cxx @@ -29,6 +29,7 @@ \author Ottavia Bertolli \author PARAPET project \author Parisa Khateri + \author Robert Twyman Skelly */ #include "stir/Scanner.h" @@ -1490,6 +1491,7 @@ Scanner::set_params(Type type_v, void Scanner::set_scanner_geometry(const std::string& new_scanner_geometry) { + info(boost::format("Scanner::set_scanner_geometry: setting scanner_geometry to %1%") % new_scanner_geometry); scanner_geometry = new_scanner_geometry; _already_setup = false; } @@ -1667,12 +1669,12 @@ Scanner::check_consistency() const // exclusion of generic as 'get_num_axial_crystals_per_block()' is sometimes false for asymmetric detectors and not // important for generic - if (dets_axial != (get_num_rings() + get_num_virtual_axial_crystals_per_block()) && scanner_geometry != "Generic") + if (dets_axial != get_num_axial_crystals() && scanner_geometry != "Generic") { warning("Scanner %s: inconsistent axial block info: %d vs %d", this->get_name().c_str(), dets_axial, - get_num_rings() + get_num_virtual_axial_crystals_per_block()); + get_num_axial_crystals()); return Succeeded::no; } } @@ -1682,17 +1684,6 @@ Scanner::check_consistency() const warning("Scanner %s: axial bucket info is not set (probably irrelevant unless you use dead-time correction that needs this " "info)", this->get_name().c_str()); - else - { - const int blocks_axial = get_num_axial_buckets() * get_num_axial_blocks_per_bucket(); - // exclusion of generic as 'get_num_axial_blocks_per_bucket()' is sometimes false for asymmetric detectors and not - // important for generic - if (blocks_axial != get_num_axial_blocks() && scanner_geometry != "Generic") - { - warning("Scanner %s: inconsistent axial block/bucket info", this->get_name().c_str()); - return Succeeded::no; - } - } } // checks on singles units { @@ -1795,13 +1786,14 @@ Scanner::check_consistency() const if (round(get_transaxial_block_spacing() * get_num_transaxial_blocks_per_bucket() * 1000.0) / 1000.0 < round(2 * inner_ring_radius * tan(_PI / 2 / get_num_transaxial_buckets()) * 1000.0) / 1000.0) { - warning("Scanner %s: inconsistent transaxial spacing:\n" - "\ttransaxial_block_spacing %f muliplied by num_transaxial_blocks_per_bucket %d should fit into a polygon that " - "encircles a cylinder with inner_ring_radius %f", - this->get_name().c_str(), - get_transaxial_block_spacing(), - get_num_transaxial_blocks_per_bucket(), - get_inner_ring_radius()); + warning( + "Scanner %s: inconsistent transaxial spacing:\n" + "\ttransaxial_block_spacing %f multiplied by num_transaxial_blocks_per_bucket %d should fit into a polygon that " + "encircles a cylinder with inner_ring_radius %f", + this->get_name().c_str(), + get_transaxial_block_spacing(), + get_num_transaxial_blocks_per_bucket(), + get_inner_ring_radius()); return Succeeded::no; } } @@ -1962,7 +1954,8 @@ Scanner::parameter_info() const } // block/bucket description - s << " Number of blocks per bucket in transaxial direction := " << get_num_transaxial_blocks_per_bucket() << '\n' + s << " Number of buckets in axial direction := " << get_num_axial_buckets() << '\n' + << " Number of blocks per bucket in transaxial direction := " << get_num_transaxial_blocks_per_bucket() << '\n' << " Number of blocks per bucket in axial direction := " << get_num_axial_blocks_per_bucket() << '\n' << " Number of crystals per block in axial direction := " << get_num_axial_crystals_per_block() << '\n' << " Number of crystals per block in transaxial direction := " << get_num_transaxial_crystals_per_block() << '\n' @@ -1978,14 +1971,14 @@ Scanner::parameter_info() const { s << " Scanner geometry (BlocksOnCylindrical/Cylindrical/Generic) := " << get_scanner_geometry() << '\n'; } - if (get_axial_crystal_spacing() >= 0) - s << " Distance between crystals in axial direction (cm) := " << get_axial_crystal_spacing() / 10 << '\n'; - if (get_transaxial_crystal_spacing() >= 0) - s << " Distance between crystals in transaxial direction (cm) := " << get_transaxial_crystal_spacing() / 10 << '\n'; - if (get_axial_block_spacing() >= 0) - s << " Distance between blocks in axial direction (cm) := " << get_axial_block_spacing() / 10 << '\n'; - if (get_transaxial_block_spacing() >= 0) - s << " Distance between blocks in transaxial direction (cm) := " << get_transaxial_block_spacing() / 10 << '\n'; + if (get_axial_crystal_spacing() >= 0 || get_num_axial_crystals_per_block() > 1) + s << " Crystal spacing/pitch in axial direction (cm) := " << get_axial_crystal_spacing() / 10 << '\n'; + if (get_transaxial_crystal_spacing() >= 0 || get_num_transaxial_crystals_per_block() > 1) + s << " Crystal spacing/pitch in transaxial direction (cm) := " << get_transaxial_crystal_spacing() / 10 << '\n'; + if (get_axial_block_spacing() >= 0 || get_num_axial_blocks_per_bucket() > 1) + s << " Block spacing/pitch in axial direction (cm) := " << get_axial_block_spacing() / 10 << '\n'; + if (get_transaxial_block_spacing() >= 0 || get_num_transaxial_blocks_per_bucket() > 1) + s << " Block spacing/pitch in transaxial direction (cm) := " << get_transaxial_block_spacing() / 10 << '\n'; s << "End scanner parameters:=\n"; @@ -2070,13 +2063,21 @@ Scanner::ask_parameters() float BinSize = ask_num("Enter default (tangential) bin size after arc-correction (in mm):", 0.F, 60.F, 3.75F); float intrTilt = ask_num("Enter intrinsic_tilt (in degrees):", -180.F, 360.F, 0.F); - int TransBlocksPerBucket = ask_num("Enter number of transaxial blocks per bucket: ", 0, 10, 2); - int AxialBlocksPerBucket = ask_num("Enter number of axial blocks per bucket: ", 0, 10, 6); - int AxialCrystalsPerBlock = ask_num("Enter number of axial crystals per block: ", 0, 16, 8); - int TransaxialCrystalsPerBlock = ask_num("Enter number of transaxial crystals per block: ", 0, 16, 8); - int AxialCrstalsPerSinglesUnit = ask_num("Enter number of axial crystals per singles unit: ", 0, NoRings, 1); + + // Helper function: if x / y = 3.33, divideAndRoundUp rounds this up to 4. If x / y = 3.0, it would leave it as 3. + auto divideAndRoundUp = [](int x, int y) -> int { return (x + y - 1) / y; }; + int AxialCrystalsPerBlock = ask_num("Enter number of axial crystals per block: ", 1, NoRings, 8); + int AxialBlocksPerBucket + = ask_num("Enter number of axial blocks per bucket: ", 1, divideAndRoundUp(NoRings, AxialCrystalsPerBlock), 6); + int AxialCrystalsPerSinglesUnit = ask_num("Enter number of axial crystals per singles unit: ", 1, NoRings, 1); + + int TransaxialCrystalsPerBlock = ask_num("Enter number of transaxial crystals per block: ", 1, num_detectors_per_ring, 8); + int TransBlocksPerBucket = ask_num("Enter number of transaxial blocks per bucket: ", + 1, + divideAndRoundUp(num_detectors_per_ring, TransaxialCrystalsPerBlock), + 2); int TransaxialCrystalsPerSinglesUnit - = ask_num("Enter number of transaxial crystals per singles unit: ", 0, num_detectors_per_ring, 1); + = ask_num("Enter number of transaxial crystals per singles unit: ", 1, num_detectors_per_ring, 1); short int Num_TOF_bins = ask_num("Number of TOF time bins :", 0, 800, 0); float Size_TOF_bin = ask_num("Size of timing bin (ps) :", 0.0f, 100.0f, 0.0f); @@ -2095,10 +2096,11 @@ Scanner::ask_parameters() const string ScannerGeometry = ask_string("Enter the scanner geometry ( BlocksOnCylindrical / Cylindrical / Generic ) :", "Cylindrical"); - float AxialCrystalSpacing = ask_num("Enter crystal spacing in axial direction (in mm): ", 0.F, 30.F, 6.75F); - float TransaxialCrystalSpacing = ask_num("Enter crystal spacing in transaxial direction (in mm): ", 0.F, 30.F, 6.75F); - float AxialBlockSpacing = ask_num("Enter block spacing in axial direction (in mm): ", 0.F, 360.F, 54.F); - float TransaxialBlockSpacing = ask_num("Enter block spacing in transaxial direction (in mm): ", 0.F, 360.F, 54.F); + float AxialCrystalSpacing = ask_num("Enter crystal spacing/pitch in axial direction (in mm): ", 0.F, 30.F, 6.75F); + float AxialBlockSpacing = ask_num("Enter block spacing/pitch in axial direction (in mm): ", 0.F, 360.F, 54.F); + float AxialBucketSpacing = ask_num("Enter bucket spacing/pitch in axial direction (in mm): ", 0.F, 360.F, 54.F); + float TransaxialCrystalSpacing = ask_num("Enter crystal spacing/pitch in transaxial direction (in mm): ", 0.F, 30.F, 6.75F); + float TransaxialBlockSpacing = ask_num("Enter block spacing/pitch in transaxial direction (in mm): ", 0.F, 360.F, 54.F); string crystal_map_file_name = ""; if (ScannerGeometry == "Generic") @@ -2123,7 +2125,7 @@ Scanner::ask_parameters() TransBlocksPerBucket, AxialCrystalsPerBlock, TransaxialCrystalsPerBlock, - AxialCrstalsPerSinglesUnit, + AxialCrystalsPerSinglesUnit, TransaxialCrystalsPerSinglesUnit, num_detector_layers, EnergyResolution, diff --git a/src/include/stir/GeometryBlocksOnCylindrical.h b/src/include/stir/GeometryBlocksOnCylindrical.h index 84677de4ff..6f076ca4de 100644 --- a/src/include/stir/GeometryBlocksOnCylindrical.h +++ b/src/include/stir/GeometryBlocksOnCylindrical.h @@ -52,6 +52,45 @@ class GeometryBlocksOnCylindrical : public DetectorCoordinateMap public: GeometryBlocksOnCylindrical(const Scanner& scanner); + //! Calculates the transaxial coordinate of a crystal given a scanner and the crystal's indices + static int + get_transaxial_coord(const Scanner& scanner, int transaxial_bucket_num, int transaxial_block_num, int transaxial_crystal_num); + + //! Calculates the axial coordinate of a crystal given a scanner and the crystal's indices + static int get_axial_coord(const Scanner& scanner, int axial_bucket_num, int axial_block_num, int axial_crystal_num); + + //! Calculates the transaxial translation of a crystal given a scanner and the crystal's indices + static float + get_crystal_in_bucket_transaxial_translation(const Scanner& scanner, int transaxial_block_num, int transaxial_crystal_num); + + //! Calculates the axial translation of a crystal given a scanner and the crystal's indices + static float get_axial_translation(const Scanner& scanner, int axial_bucket_num, int axial_block_num, int axial_crystal_num); + + //! Calculate the initial axial z offset to center the scanner on 0,0,0 + static float get_initial_axial_z_offset(const Scanner& scanner); + + //! Calculate the initial transaxial x offset to center the scanner on 0,0,0 + static float get_initial_axial_x_offset_for_each_bucket(const Scanner& scanner); + + static float get_csi_minus_csi_gaps(const Scanner& scanner) + { + //! Calculate the CSI, angle covered by half a bucket + // 2 * PI / num_transaxial_buckets / 2 (simplified) + float csi = _PI / scanner.get_num_transaxial_buckets(); // TODO, this assumes 1 transaxial bucket per angle + + // The difference between the transaxial block spacing and the sum of all transaxial crystal spacing's in the block + float trans_blocks_gap = scanner.get_transaxial_block_spacing() + - scanner.get_num_transaxial_crystals_per_block() * scanner.get_transaxial_crystal_spacing(); + // Calculate the angle covered by the gaps between the blocks + float csi_gaps + = 2 * csi * (scanner.get_transaxial_crystal_spacing() / 2 + trans_blocks_gap) / scanner.get_transaxial_block_spacing(); + return csi - csi_gaps; + }; + + //! + CartesianCoordinate3D calculate_crystal_rotation(const CartesianCoordinate3D& crystal_position, + const float alpha) const; + private: //! Get rotation matrix for a given angle around z axis stir::Array<2, float> get_rotation_matrix(float alpha) const; diff --git a/src/include/stir/Scanner.h b/src/include/stir/Scanner.h index 79c8cbde82..7f60066066 100644 --- a/src/include/stir/Scanner.h +++ b/src/include/stir/Scanner.h @@ -5,7 +5,8 @@ Copyright (C) 2016, University of Hull Copyright (C) 2016, 2019, 2021, 2023 UCL Copyright 2017 ETH Zurich, Institute of Particle Physics and Astrophysics - Copyright (C 2017-2018, University of Leeds + Copyright (C) 2017-2018, University of Leeds + Copyright (C) 2024, Prescient Imaging This file is part of STIR. SPDX-License-Identifier: Apache-2.0 AND License-ref-PARAPET-license @@ -28,6 +29,7 @@ \author Palak Wadhwa \author PARAPET project \author Parisa Khateri + \author Robert Twyman Skelly */ #ifndef __stir_buildblock_SCANNER_H__ @@ -332,7 +334,7 @@ class Scanner algorithm. It uses the same coordinate system as ProjDataInfo::get_phi(). */ inline float get_intrinsic_azimuthal_tilt() const; - //! \name Info on crystals per block etc. + //! \name Info on crystals per block/bucket etc. //@{ //! get number of transaxial blocks per bucket inline int get_num_transaxial_blocks_per_bucket() const; @@ -348,6 +350,10 @@ class Scanner inline int get_num_axial_crystals_per_bucket() const; //! get number of crystal layers (for DOI) inline int get_num_detector_layers() const; + /*! With virtual crystals between blocks, it will miss the end virtual crystals. This method takes this into account.*/ + inline int get_num_axial_crystals() const; + //! get the number of transaxial crystals + inline int get_num_transaxial_crystals() const; //! get number of axial blocks inline int get_num_axial_blocks() const; //! get number of axial blocks @@ -416,13 +422,13 @@ class Scanner //! get scanner geometry /*! \see set_scanner_geometry */ inline std::string get_scanner_geometry() const; - //! get crystal spacing in axial direction + //! get crystal spacing/pitch in axial direction in mm inline float get_axial_crystal_spacing() const; - //! get crystal spacing in transaxial direction + //! get crystal spacing/pitch in transaxial direction in mm inline float get_transaxial_crystal_spacing() const; - //! get block spacing in axial direction + //! get block spacing/pitch in axial direction in mm inline float get_axial_block_spacing() const; - //! get block spacing in transaxial direction + //! get block spacing/pitch in transaxial direction in mm inline float get_transaxial_block_spacing() const; //@} (end of get block geometry info) @@ -478,7 +484,7 @@ class Scanner inline void set_default_bin_size(const float& new_size); //! in degrees inline void set_intrinsic_azimuthal_tilt(const float new_tilt); - //! \name Info on crystals per block etc. + //! \name Info on crystals per block/bucket etc. //@{ //! set number of transaxial blocks per bucket inline void set_num_transaxial_blocks_per_bucket(const int& new_num); diff --git a/src/include/stir/Scanner.inl b/src/include/stir/Scanner.inl index 874f97e45c..f9f142a66f 100644 --- a/src/include/stir/Scanner.inl +++ b/src/include/stir/Scanner.inl @@ -22,6 +22,7 @@ \author Long Zhang (set*() functions) \author PARAPET project \author Parisa Khateri + \author Robert Twyman Skelly */ @@ -143,13 +144,13 @@ Scanner::get_num_transaxial_crystals_per_block() const int Scanner::get_num_axial_crystals_per_bucket() const { - return get_num_axial_blocks_per_bucket() * get_num_axial_crystals_per_block(); + return get_num_axial_crystals_per_block() * get_num_axial_blocks_per_bucket(); } int Scanner::get_num_transaxial_crystals_per_bucket() const { - return get_num_transaxial_blocks_per_bucket() * get_num_transaxial_crystals_per_block(); + return get_num_transaxial_crystals_per_block() * get_num_transaxial_blocks_per_bucket(); } int @@ -159,29 +160,41 @@ Scanner::get_num_detector_layers() const } int -Scanner::get_num_axial_blocks() const +Scanner::get_num_axial_crystals() const { - // when using virtual crystals between blocks, there won't be one at the end, so we + // when using virtual crystals between blocks, there won't be any at the end, so we // need to take this into account. - return (num_rings + get_num_virtual_axial_crystals_per_block()) / num_axial_crystals_per_block; + return num_rings + get_num_virtual_axial_crystals_per_block(); +} + +int +Scanner::get_num_transaxial_crystals() const +{ + return num_detectors_per_ring; +} + +int +Scanner::get_num_axial_blocks() const +{ + return get_num_axial_crystals() / get_num_axial_crystals_per_block(); } int Scanner::get_num_transaxial_blocks() const { - return num_detectors_per_ring / num_transaxial_crystals_per_block; + return get_num_transaxial_crystals() / get_num_transaxial_crystals_per_block(); } int Scanner::get_num_axial_buckets() const { - return get_num_axial_blocks() / num_axial_blocks_per_bucket; + return get_num_axial_blocks() / get_num_axial_blocks_per_bucket(); } int Scanner::get_num_transaxial_buckets() const { - return get_num_transaxial_blocks() / num_transaxial_blocks_per_bucket; + return get_num_transaxial_blocks() / get_num_transaxial_blocks_per_bucket(); } int @@ -205,7 +218,7 @@ Scanner::get_num_axial_singles_units() const } else { - return (num_rings + get_num_virtual_axial_crystals_per_block()) / num_axial_crystals_per_singles_unit; + return get_num_axial_crystals() / num_axial_crystals_per_singles_unit; } } @@ -218,7 +231,8 @@ Scanner::get_num_transaxial_singles_units() const } else { - return num_detectors_per_ring / num_transaxial_crystals_per_singles_unit; + // TODO: RTS believe this should account for virtual crystals, check! + return get_num_detectors_per_ring() / num_transaxial_crystals_per_singles_unit; } } diff --git a/src/recon_test/test_geometry_blocks_on_cylindrical.cxx b/src/recon_test/test_geometry_blocks_on_cylindrical.cxx index e88c4e08d1..81ac8f500e 100644 --- a/src/recon_test/test_geometry_blocks_on_cylindrical.cxx +++ b/src/recon_test/test_geometry_blocks_on_cylindrical.cxx @@ -27,6 +27,7 @@ #include "stir/GeometryBlocksOnCylindrical.h" #include "stir/IO/write_to_file.h" #include +// #include START_NAMESPACE_STIR @@ -40,12 +41,26 @@ class GeometryBlocksOnCylindricalTests : public RunTests void run_tests() override; private: + //! Loop through all transaxial and axial indices to ensure the coordinates are monotonic, indicating a spiralling + // generation of the coordinates of the coordinates around the scanner geometry + void run_monotonic_coordinates_generation_test(); /*! \brief Tests multiple axial blocks/bucket configurations to ensure the detector map's axial indices and coordinates - * are monotonic - */ + * are monotonic */ void run_monotonic_axial_coordinates_in_detector_map_test(); //! Tests the axial indices and coordinates are monotonic in the detector map static Succeeded monotonic_axial_coordinates_in_detector_map_test(const shared_ptr& scanner_sptr); + + //! Sets up the scanner for the test + static Succeeded setup_scanner_for_test(shared_ptr scanner_sptr); + + void run_assert_scanner_centred_on_origin_test(); + + /*! This is a test of the start_z position. The code was refactored and this test ensures that the new calculation is + * equivalent to the old calculation */ + void validate_start_z_with_old_calculation(); + + void validate_first_bucket_is_centred_on_x_axis(); + ; }; void @@ -53,34 +68,44 @@ GeometryBlocksOnCylindricalTests::run_monotonic_axial_coordinates_in_detector_ma { auto scanner_sptr = std::make_shared(Scanner::SAFIRDualRingPrototype); scanner_sptr->set_scanner_geometry("BlocksOnCylindrical"); - scanner_sptr->set_transaxial_block_spacing(scanner_sptr->get_transaxial_crystal_spacing() - * scanner_sptr->get_num_transaxial_crystals_per_block()); - int num_axial_buckets = 1; // TODO add for loop when support is added - for (int num_axial_crystals_per_blocks = 1; num_axial_crystals_per_blocks < 3; ++num_axial_crystals_per_blocks) - for (int num_axial_blocks_per_bucket = 1; num_axial_blocks_per_bucket < 3; ++num_axial_blocks_per_bucket) + // for (int num_axial_crystals_per_blocks = 1; num_axial_crystals_per_blocks < 3; ++num_axial_crystals_per_blocks) + // for (int num_axial_blocks_per_bucket = 1; num_axial_blocks_per_bucket < 3; ++num_axial_blocks_per_bucket) + // for (int num_axial_buckets = 1; num_axial_buckets < 3; ++num_axial_buckets) + + // TESTING CONFIG: + int num_axial_buckets = 2; + int num_axial_blocks_per_bucket = 2; + int num_axial_crystals_per_blocks = 2; + { + int num_rings = num_axial_crystals_per_blocks * num_axial_blocks_per_bucket * num_axial_buckets; + scanner_sptr->set_num_axial_crystals_per_block(num_axial_crystals_per_blocks); + scanner_sptr->set_num_axial_blocks_per_bucket(num_axial_blocks_per_bucket); + scanner_sptr->set_num_rings(num_rings); + + scanner_sptr->set_axial_block_spacing(scanner_sptr->get_axial_crystal_spacing() + * (scanner_sptr->get_num_axial_crystals_per_block() + 0.5)); + + // if (num_axial_buckets > 1) + // scanner_sptr->set_axial_bucket_spacing(scanner_sptr->get_axial_block_spacing() * 1.5); + // else + // scanner_sptr->set_axial_bucket_spacing(-1); + + if (monotonic_axial_coordinates_in_detector_map_test(scanner_sptr) == Succeeded::no) { - scanner_sptr->set_num_axial_crystals_per_block(num_axial_crystals_per_blocks); - scanner_sptr->set_num_axial_blocks_per_bucket(num_axial_blocks_per_bucket); - scanner_sptr->set_num_rings(scanner_sptr->get_num_axial_crystals_per_bucket() * num_axial_buckets); - scanner_sptr->set_axial_block_spacing(scanner_sptr->get_axial_crystal_spacing() - * (scanner_sptr->get_num_axial_crystals_per_block() + 0.5)); - - if (monotonic_axial_coordinates_in_detector_map_test(scanner_sptr) == Succeeded::no) - { - warning(boost::format("Monothonic axial coordinates test failed for:\n" - "\taxial_crystal_per_block =\t%1%\n" - "\taxial_blocks_per_bucket =\t%2%\n" - "\tnum_axial_buckets =\t\t\t%3%") - % num_axial_crystals_per_blocks % num_axial_blocks_per_bucket % num_axial_buckets); - everything_ok = false; - return; - } + warning(boost::format("Monothonic axial coordinates test failed for:\n" + "\taxial_crystal_per_block =\t%1%\n" + "\taxial_blocks_per_bucket =\t%2%\n" + "\tnum_axial_buckets =\t\t\t%3%") + % num_axial_crystals_per_blocks % num_axial_blocks_per_bucket % scanner_sptr->get_num_axial_buckets()); + everything_ok = false; + return; } + } } Succeeded -GeometryBlocksOnCylindricalTests::monotonic_axial_coordinates_in_detector_map_test(const shared_ptr& scanner_sptr) +GeometryBlocksOnCylindricalTests::setup_scanner_for_test(shared_ptr scanner_sptr) { if (scanner_sptr->get_scanner_geometry() != "BlocksOnCylindrical") { @@ -88,10 +113,9 @@ GeometryBlocksOnCylindricalTests::monotonic_axial_coordinates_in_detector_map_te return Succeeded::no; } - shared_ptr detector_map_sptr; try { - detector_map_sptr.reset(new GeometryBlocksOnCylindrical(*scanner_sptr)); + scanner_sptr->set_up(); } catch (const std::runtime_error& e) { @@ -100,41 +124,255 @@ GeometryBlocksOnCylindricalTests::monotonic_axial_coordinates_in_detector_map_te % e.what()); return Succeeded::no; } + return Succeeded::yes; +} + +Succeeded +GeometryBlocksOnCylindricalTests::monotonic_axial_coordinates_in_detector_map_test(const shared_ptr& scanner_sptr) +{ + if (setup_scanner_for_test(scanner_sptr) == Succeeded::no) + return Succeeded::no; unsigned min_axial_pos = 0; float prev_min_axial_coord = -std::numeric_limits::max(); + shared_ptr detector_map_sptr = scanner_sptr->get_detector_map_sptr(); for (unsigned axial_idx = 0; axial_idx < detector_map_sptr->get_num_axial_coords(); ++axial_idx) - for (unsigned tangential_idx = 0; tangential_idx < detector_map_sptr->get_num_tangential_coords(); ++tangential_idx) - for (unsigned radial_idx = 0; radial_idx < detector_map_sptr->get_num_radial_coords(); ++radial_idx) + { + const DetectionPosition<> det_pos = DetectionPosition<>(0, axial_idx, 0); + CartesianCoordinate3D coord = detector_map_sptr->get_coordinate_for_det_pos(det_pos); + // std::cerr << "coord.z() = " << coord.z() << "\tprev_min_axial_coord = " << prev_min_axial_coord + // << "\tdelta = " << coord.z() - prev_min_axial_coord << std::endl; + if (coord.z() > prev_min_axial_coord) + { + min_axial_pos = axial_idx; + prev_min_axial_coord = coord.z(); + } + else if (coord.z() < prev_min_axial_coord) { - const DetectionPosition<> det_pos = DetectionPosition<>(tangential_idx, axial_idx, radial_idx); - CartesianCoordinate3D coord = detector_map_sptr->get_coordinate_for_det_pos(det_pos); - if (coord.z() > prev_min_axial_coord) - { - min_axial_pos = axial_idx; - prev_min_axial_coord = coord.z(); - } - else if (coord.z() < prev_min_axial_coord) - { - float delta = coord.z() - prev_min_axial_coord; - warning(boost::format("Axial Coordinates are not monotonic.\n" - "Next axial index =\t\t%1%, Next axial coord (mm) =\t\t%2% (%3%)\n" - "Previous axial index =\t%4%, Previous axial coord (mm) =\t%5%") - % axial_idx % coord.z() % delta % min_axial_pos % prev_min_axial_coord); - return Succeeded::no; - } + float delta = coord.z() - prev_min_axial_coord; + warning(boost::format("Axial Coordinates are not monotonic.\n" + "Next axial index =\t\t%1%, Next axial coord (mm) =\t\t%2% (%3%)\n" + "Previous axial index =\t%4%, Previous axial coord (mm) =\t%5%") + % axial_idx % coord.z() % delta % min_axial_pos % prev_min_axial_coord); + return Succeeded::no; } + } return Succeeded::yes; } +void +GeometryBlocksOnCylindricalTests::run_assert_scanner_centred_on_origin_test() +{ + // auto scanner_sptr = std::make_shared(Scanner::SAFIRDualRingPrototype); + // scanner_sptr->set_scanner_geometry("BlocksOnCylindrical"); + // if (setup_scanner_for_test(scanner_sptr) == Succeeded::no) + // { + // warning("GeometryBlocksOnCylindricalTests::run_assert_scanner_centred_on_origin_test: " + // "Scanner not set up correctly for test"); + // everything_ok = false; + // } + // + // auto detector_map_sptr = scanner_sptr->get_detector_map_sptr(); + // for (int transaxial_idx = 0; transaxial_idx < scanner_sptr->get_num_detectors_per_ring(); ++transaxial_idx) + // { + // const DetectionPosition<> det_pos = DetectionPosition<>(transaxial_idx, 0, 0); + // auto cart_coord = detector_map_sptr->get_coordinate_for_det_pos(det_pos); + // if (cart_coord.z() != 0.0) + // { + // warning(boost::format("Scanner Z component of the first index is > 0 \n" + // "Transaxial index =\t%1%\n" + // "Axial index =\t%2%\n" + // "Cartesian Coordinate =\t%2%") + // % transaxial_idx % cart_coord); + // everything_ok = false; + // } + // } +} + +void +GeometryBlocksOnCylindricalTests::run_monotonic_coordinates_generation_test() +{ + auto scanner_sptr = std::make_shared(Scanner::SAFIRDualRingPrototype); + scanner_sptr->set_scanner_geometry("BlocksOnCylindrical"); + + int prev_max_axial_coord = -1; + int prev_max_transaxial_coord = -1; + + for (int ax_bucket_num = 0; ax_bucket_num < scanner_sptr->get_num_axial_buckets(); ++ax_bucket_num) + for (int ax_block_num = 0; ax_block_num < scanner_sptr->get_num_axial_blocks_per_bucket(); ++ax_block_num) + for (int ax_crystal_num = 0; ax_crystal_num < scanner_sptr->get_num_axial_crystals_per_block(); ++ax_crystal_num) + for (int trans_bucket_num = 0; trans_bucket_num < scanner_sptr->get_num_transaxial_buckets(); ++trans_bucket_num) + for (int trans_block_num = 0; trans_block_num < scanner_sptr->get_num_transaxial_blocks_per_bucket(); ++trans_block_num) + for (int trans_crys_num = 0; trans_crys_num < scanner_sptr->get_num_transaxial_crystals_per_block(); ++trans_crys_num) + { + int axial_coord + = GeometryBlocksOnCylindrical::get_axial_coord(*scanner_sptr, ax_bucket_num, ax_block_num, ax_crystal_num); + int transaxial_coord; + transaxial_coord = GeometryBlocksOnCylindrical::get_transaxial_coord( + *scanner_sptr, trans_bucket_num, trans_block_num, trans_crys_num); + + // Test that the axial coordinates are monotonic + if (prev_max_axial_coord > axial_coord) + { + warning(boost::format("Axial Coordinates are not monotonic.\n" + "Next axial index =\t\t%1%\n" + "Previous axial index =\t%2%\n" + "Next Transaxial index =\t%3%\n" + "Previous Transaxial index =\t%4%") + % axial_coord % prev_max_axial_coord % transaxial_coord % prev_max_transaxial_coord); + everything_ok = false; + return; + } + + // Test that the transaxial coordinates are monotonic, it will reset to 0 when the axial index increases + if (prev_max_transaxial_coord > transaxial_coord && axial_coord <= prev_max_axial_coord) + { + warning(boost::format("Transaxial Coordinates are not monotonic.\n" + "Next axial index =\t\t%1%\n" + "Previous axial index =\t%2%\n" + "Next Transaxial index =\t%3%\n" + "Previous Transaxial index =\t%4%") + % axial_coord % prev_max_axial_coord % transaxial_coord % prev_max_transaxial_coord); + everything_ok = false; + return; + } + } +} + +void +GeometryBlocksOnCylindricalTests::validate_start_z_with_old_calculation() +{ + auto scanner_sptr = std::make_shared(Scanner::SAFIRDualRingPrototype); + scanner_sptr->set_scanner_geometry("BlocksOnCylindrical"); + + // calculate the start_z with the old equation. + // The code was moved here and slightly refactored to get values from the scanner object + float axial_blocks_gap = scanner_sptr->get_axial_block_spacing() + - scanner_sptr->get_num_axial_crystals_per_block() * scanner_sptr->get_axial_crystal_spacing(); + float old_code_calculation + = -(scanner_sptr->get_axial_block_spacing() * (scanner_sptr->get_num_axial_blocks_per_bucket()) + * scanner_sptr->get_num_axial_buckets() + - scanner_sptr->get_axial_crystal_spacing() + - axial_blocks_gap * (scanner_sptr->get_num_axial_blocks_per_bucket() * scanner_sptr->get_num_axial_buckets() - 1)) + / 2; + + // new method to calculate start_z + float start_z_method = GeometryBlocksOnCylindrical::get_initial_axial_z_offset(*scanner_sptr); + + if (std::abs(old_code_calculation - start_z_method) > 1e-6) + { + warning(boost::format("Old code calculation =\t%1%\n" + "New code calculation =\t%2%\n" + "Difference =\t%3%") + % old_code_calculation % start_z_method % (old_code_calculation - start_z_method)); + everything_ok = false; + } +} + +void +GeometryBlocksOnCylindricalTests::validate_first_bucket_is_centred_on_x_axis() +{ + auto scanner_sptr = std::make_shared(Scanner::SAFIRDualRingPrototype); + scanner_sptr->set_scanner_geometry("BlocksOnCylindrical"); + scanner_sptr->set_num_transaxial_blocks_per_bucket(1); + scanner_sptr->set_intrinsic_azimuthal_tilt(-0.235307813); // Required by test calculations + + GeometryBlocksOnCylindrical detector_map_builder(*scanner_sptr); + scanner_sptr->get_num_transaxial_blocks_per_bucket(); + int first_transaxial_index_in_first_bucket = detector_map_builder.get_transaxial_coord(*scanner_sptr, 0, 0, 0); + int last_transaxial_index_in_first_bucket + = detector_map_builder.get_transaxial_coord(*scanner_sptr, + 0, + scanner_sptr->get_num_transaxial_blocks_per_bucket() - 1, + scanner_sptr->get_num_transaxial_crystals_per_block()); + auto first_transaxial_coord_in_first_bucket + = detector_map_builder.get_coordinate_for_det_pos(DetectionPosition<>(first_transaxial_index_in_first_bucket, 0, 0)); + auto last_transaxial_coord_in_first_bucket + = detector_map_builder.get_coordinate_for_det_pos(DetectionPosition<>(last_transaxial_index_in_first_bucket, 0, 0)); + + { // Testing the coordinates of the first axial row of the first bucket. + // For this scanner, with the set intrinsic_azimuthal_tilt, the first bucket's tangential coordinates should + // be centered at x=0 with a constant z,y position + + // Get the coordinates of the first block/bucket + std::vector> crystal_coords(scanner_sptr->get_num_transaxial_crystals_per_bucket()); + for (int i = 0; i < crystal_coords.size(); i++) + { + crystal_coords[i] = detector_map_builder.get_coordinate_for_det_pos(DetectionPosition<>(i, 0, 0)); + } + + { + float previous_x = -std::numeric_limits::max(); + for (int i = 1; i < crystal_coords.size(); i++) + { + // Check if the crystal coords are all at the same z,y position + check_if_equal(crystal_coords[0].z(), crystal_coords[i].z()); + check_if_equal(crystal_coords[0].y(), crystal_coords[i].y()); + // Check the x position is monotonic + check_if_less(previous_x, + crystal_coords[i].x(), + boost::str(boost::format("Crystal %1% x position is more than " + "the previous crystal x position") + % i)); + } + } + + // Check if x position is centered on 0.0 + for (int i = 0; i < int(crystal_coords.size() / 2); i++) + { + float left_crystal_x = crystal_coords[i].x(); + float right_crystal_x = crystal_coords[crystal_coords.size() - 1 - i].x(); + float offset_x = abs(left_crystal_x + right_crystal_x) / 2; + float acceptable_error = (abs(left_crystal_x) + abs(right_crystal_x)) * 1e-6; + check_if_less(offset_x, + acceptable_error, + boost::str(boost::format("Left and right crystal x positions are not evenly " + "distributed over x axis.\tLeft crystal x = %1%\tRight crystal x = " + "%2%\tOffset = %3%") + % left_crystal_x % right_crystal_x % offset_x)); + } + + // Check the central crystal x position is centered on the x axis if there is an odd number of crystals + if (crystal_coords.size() % 2 == 1) + { + float central_x = crystal_coords[crystal_coords.size() / 2].x(); + float acceptable_error = abs(crystal_coords[0].x()) * 1e-6; + check_if_less(abs(central_x), + acceptable_error, + boost::str(boost::format("With an odd number of crystals, the central crystal x position is not " + "centered on the x axis.\tCentral crystal x = %1%\t Acceptable error = %2%") + % central_x % acceptable_error)); + } + } + + // { // Loop over all transaxial buckets, blocks and crystals and save the coordinates to a csv file with the format + // // tang,axial,radial,z,y,x + // std::ofstream file; + // + // const char* filename = "crystal_coords.csv"; + // file.open(filename); + // for (int i = 0; i < scanner_sptr->get_num_detectors_per_ring(); i++) + // { + // auto coord = detector_map_builder.get_coordinate_for_det_pos(DetectionPosition<>(i, 0, 0)); + // file << i << ",0,0," << coord.z() << "," << coord.y() << "," << coord.x() << std::endl; + // } + // file.close(); + // std::cerr << "Crystal coordinates written to " << std::filesystem::current_path() / filename << std::endl; + // } +} + void GeometryBlocksOnCylindricalTests::run_tests() { HighResWallClockTimer timer; timer.start(); + run_monotonic_coordinates_generation_test(); run_monotonic_axial_coordinates_in_detector_map_test(); + run_assert_scanner_centred_on_origin_test(); + validate_start_z_with_old_calculation(); + validate_first_bucket_is_centred_on_x_axis(); timer.stop(); } END_NAMESPACE_STIR