Line data Source code
1 : #include "EllipCylinder.h"
2 :
3 : namespace elsa::phantoms
4 : {
5 :
6 : template <typename data_t>
7 : EllipCylinder<data_t>::EllipCylinder(Orientation o, data_t amplit, elsa::phantoms::Vec3i center,
8 : Vec2X<data_t> halfAxis, data_t length)
9 : : _orientation{o}, _amplit{amplit}, _center{center}, _halfAxis{halfAxis}, _length{length}
10 61 : {
11 :
12 61 : aSqr = _halfAxis[INDEX_A] * _halfAxis[INDEX_A];
13 61 : bSqr = _halfAxis[INDEX_B] * _halfAxis[INDEX_B];
14 61 : aSqrbSqr = aSqr * bSqr;
15 61 : };
16 :
17 : template <typename data_t>
18 : bool EllipCylinder<data_t>::isInEllipCylinder(const Vec3i& idx) const
19 1584 : {
20 : // Depending on the orientation only one part is compiled
21 1584 : if (_orientation == Orientation::X_AXIS) {
22 46 : return (double(idx[INDEX_Y] * idx[INDEX_Y]) * bSqr
23 46 : + double(idx[INDEX_Z] * idx[INDEX_Z]) * aSqr)
24 46 : <= aSqrbSqr;
25 1538 : } else if (_orientation == Orientation::Y_AXIS) {
26 32 : return (double(idx[INDEX_X] * idx[INDEX_X]) * bSqr
27 32 : + double(idx[INDEX_Z] * idx[INDEX_Z]) * aSqr)
28 32 : <= aSqrbSqr;
29 1506 : } else {
30 1506 : return (double(idx[INDEX_X] * idx[INDEX_X]) * bSqr
31 1506 : + double(idx[INDEX_Y] * idx[INDEX_Y]) * aSqr)
32 1506 : <= aSqrbSqr;
33 1506 : }
34 1584 : };
35 :
36 : /**
37 : * @brief returns the min and max value of the axis to allow length bigger than the dimension of
38 : * the datacontainer. If the length is inside the container, there is no clipping. Return value
39 : * are final voxel indices
40 : *
41 : */
42 : std::pair<index_t, index_t> checkBounds(index_t halfLength, index_t center,
43 : index_t maxDataContainer)
44 61 : {
45 61 : index_t min, max;
46 61 : if (halfLength + center < maxDataContainer) {
47 55 : max = halfLength + center;
48 55 : } else {
49 6 : max = maxDataContainer;
50 6 : }
51 :
52 61 : if (center - halfLength > 0) {
53 57 : min = center - halfLength;
54 57 : } else {
55 4 : min = 0;
56 4 : }
57 :
58 61 : return {min, max};
59 61 : }
60 :
61 : template <Blending b, typename data_t>
62 : void rasterize(EllipCylinder<data_t>& el, VolumeDescriptor& dd, DataContainer<data_t>& dc)
63 61 : {
64 : // check ellipse in origin, mirror 4 points, draw line in the orthogonal directions with
65 : // length
66 :
67 61 : Vec3i idx(3);
68 61 : idx << 0, 0, 0;
69 :
70 61 : auto strides = dd.getProductOfCoefficientsPerDimension();
71 :
72 : // global offests for a fast memory reuse in the for loops
73 61 : index_t aOffset = 0;
74 61 : index_t aOffsetNeg = 0;
75 :
76 61 : index_t bOffset = 0;
77 61 : index_t bOffsetNeg = 0;
78 :
79 61 : index_t cOffset = 0;
80 :
81 61 : Vec3i _center = el.getCenter();
82 61 : data_t _amplit = el.getAmplitude();
83 :
84 61 : int TEMP_A = INDEX_X;
85 61 : int TEMP_B = INDEX_Y;
86 61 : int TEMP_C = INDEX_Z;
87 :
88 61 : if (el.getOrientation() == Orientation::X_AXIS) {
89 4 : TEMP_A = INDEX_Y;
90 4 : TEMP_B = INDEX_Z;
91 4 : TEMP_C = INDEX_X;
92 :
93 57 : } else if (el.getOrientation() == Orientation::Y_AXIS) {
94 4 : TEMP_A = INDEX_X;
95 4 : TEMP_B = INDEX_Z;
96 4 : TEMP_C = INDEX_Y;
97 4 : }
98 :
99 61 : index_t halfLength = std::lround(el.getLength() / 2);
100 :
101 61 : auto [minC, maxC] = checkBounds(halfLength, _center[INDEX_C],
102 61 : dd.getNumberOfCoefficientsPerDimension()[INDEX_C] - 1);
103 :
104 : // check ellipse on AxB Plane, draw ellipse along the C axis from minC to maxC
105 264 : for (; el.isInEllipCylinder(idx); idx[TEMP_A]++) {
106 203 : aOffset = (idx[TEMP_A] + _center[TEMP_A]) * strides[TEMP_A];
107 203 : aOffsetNeg = (-idx[TEMP_A] + _center[TEMP_A]) * strides[TEMP_A];
108 :
109 1320 : for (; el.isInEllipCylinder(idx); idx[TEMP_B]++) {
110 1117 : bOffset = (idx[TEMP_B] + _center[TEMP_B]) * strides[TEMP_B];
111 1117 : bOffsetNeg = (-idx[TEMP_B] + _center[TEMP_B]) * strides[TEMP_B];
112 :
113 53674 : for (index_t line = minC; line <= maxC; line++) {
114 52557 : cOffset = line * strides[TEMP_C];
115 :
116 52557 : blend<b>(dc, aOffset + bOffset + cOffset, _amplit);
117 :
118 52557 : if (idx[TEMP_A] != 0) {
119 49246 : blend<b>(dc, aOffsetNeg + bOffset + cOffset, _amplit);
120 49246 : }
121 :
122 52557 : if (idx[TEMP_B] != 0) {
123 47556 : blend<b>(dc, aOffset + bOffsetNeg + cOffset, _amplit);
124 47556 : }
125 :
126 52557 : if (idx[TEMP_A] != 0 && idx[TEMP_B] != 0) {
127 44786 : blend<b>(dc, aOffsetNeg + bOffsetNeg + cOffset, _amplit);
128 44786 : }
129 52557 : };
130 1117 : idx[TEMP_C] = 0;
131 1117 : }
132 203 : idx[TEMP_B] = 0;
133 203 : }
134 61 : };
135 :
136 : // ------------------------------------------
137 : // explicit template instantiation
138 : template class EllipCylinder<float>;
139 : template class EllipCylinder<double>;
140 :
141 : template void rasterize<Blending::ADDITION, float>(EllipCylinder<float>& el,
142 : VolumeDescriptor& dd,
143 : DataContainer<float>& dc);
144 : template void rasterize<Blending::ADDITION, double>(EllipCylinder<double>& el,
145 : VolumeDescriptor& dd,
146 : DataContainer<double>& dc);
147 :
148 : template void rasterize<Blending::OVERWRITE, float>(EllipCylinder<float>& el,
149 : VolumeDescriptor& dd,
150 : DataContainer<float>& dc);
151 : template void rasterize<Blending::OVERWRITE, double>(EllipCylinder<double>& el,
152 : VolumeDescriptor& dd,
153 : DataContainer<double>& dc);
154 :
155 : } // namespace elsa::phantoms
|
