-
Notifications
You must be signed in to change notification settings - Fork 95
/
Copy pathGeometryBlocksOnCylindrical.cxx
148 lines (122 loc) · 6.96 KB
/
GeometryBlocksOnCylindrical.cxx
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
/*
Copyright 2017 ETH Zurich, Institute of Particle Physics and Astrophysics
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/
/*!
\file
\ingroup projdata
\brief Non-inline implementations of stir::GeometryBlocksOnCylindrical
\author Parisa Khateri
*/
#include "stir/DetectionPosition.h"
#include "stir/CartesianCoordinate3D.h"
#include "stir/Scanner.h"
#include "stir/shared_ptr.h"
#include "stir/GeometryBlocksOnCylindrical.h"
#include <string>
#include <cmath>
#include "stir/Array.h"
#include "stir/make_array.h"
#include "stir/numerics/MatrixFunction.h"
#include <map>
#include <iostream>
#include <fstream>
#include <iomanip>
#include <boost/format.hpp>
START_NAMESPACE_STIR
GeometryBlocksOnCylindrical::GeometryBlocksOnCylindrical(const Scanner& scanner)
{
if (scanner.check_consistency() == Succeeded::no)
error("Error in GeometryBlocksOnCylindrical: scanner configuration not accepted. Please check warnings.");
build_crystal_maps(scanner);
}
stir::Array<2, float>
GeometryBlocksOnCylindrical::get_rotation_matrix(float alpha) const
{
return stir::make_array(stir::make_1d_array(1.F, 0.F, 0.F),
stir::make_1d_array(0.F, std::cos(alpha), std::sin(alpha)),
stir::make_1d_array(0.F, -1 * std::sin(alpha), std::cos(alpha)));
}
void
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_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);
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;
}
set_detector_map(cartesian_coord_map_given_detection_position_keys);
}
END_NAMESPACE_STIR