Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

BlocksOnCylindrical refactor and add axial bucket gaps #1400

Draft
wants to merge 2 commits into
base: master
Choose a base branch
from
Draft

BlocksOnCylindrical refactor and add axial bucket gaps #1400

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
5c3c804
Refactor scanner methods
robbietuk Mar 19, 2024
84b641c
Refactor build_crystal_maps and fix some issues with initial x position
robbietuk Mar 19, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
192 changes: 122 additions & 70 deletions src/buildblock/GeometryBlocksOnCylindrical.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 <map>
#include <iostream>
#include <fstream>
Expand Down Expand Up @@ -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<float> 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<float> 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<float> 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<float> transformed_coord = start_point + transformation_matrix;
stir::CartesianCoordinate3D<float> 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<float> translation_matrix(
axial_translation,
0.,
get_crystal_in_bucket_transaxial_translation(scanner, trans_block_num, trans_crys_num));

stir::CartesianCoordinate3D<float> 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<float>
GeometryBlocksOnCylindrical::calculate_crystal_rotation(const CartesianCoordinate3D<float>& 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
Loading
Loading