diff --git a/src/detectors/curved_detector.m b/src/detectors/curved_detector.m index f17c166..a496958 100644 --- a/src/detectors/curved_detector.m +++ b/src/detectors/curved_detector.m @@ -6,7 +6,7 @@ pixel_height (1, 1) double % The height of each pixel in z end methods (Access=private, Static) - function generator = get_ray_generator_static(ny_pixels, nz_pixels, pixel_angle, pixel_height, dist_to_detector, to_source_vec, source_pos) + function generator = get_ray_generator_static(ny_pixels, nz_pixels, pixel_angle, pixel_height, dist_to_detector, to_source_vec, source_pos, voxels) generator = @get_ray_attrs; % Create the function which returns the rays function [xray] = get_ray_attrs(y_pixel, z_pixel) assert(y_pixel <= ny_pixels && y_pixel > 0 && z_pixel <= nz_pixels && z_pixel > 0, ... @@ -15,7 +15,7 @@ final_length = sqrt(dist_to_detector.^2 + z_shift.^2); pixel_vec = (rotz(pixel_angle * (y_pixel - (ny_pixels+1)/2)) * to_source_vec.*dist_to_detector - ... [0;0;z_shift]) ./ final_length; - xray = ray(source_pos, -pixel_vec, final_length); + xray = ray(source_pos, -pixel_vec, final_length, voxels); end end end @@ -52,30 +52,32 @@ function reset(self) self.source_position = self.init_source_pos; end - function ray_generator = get_ray_generator(self, ray_per_pixel) + function ray_generator = get_ray_generator(self, voxels, ray_per_pixel) % Create a function which returns the rays which should be fired to hit each pixel. % Only 1 ray per pixel is supported at the moment, as anti-aliasing techniques are not yet implemented. arguments self curved_detector + voxels voxel_array ray_per_pixel int32 = 1 end - assert(nargin==1, "Only 1 ray per pixel is supported at the moment, as anti-aliasing techniques are not yet implemented.") + assert(nargin==2, "Only 1 ray per pixel is supported at the moment, as anti-aliasing techniques are not yet implemented.") % Create the function which returns the rays ray_generator = self.get_ray_generator_static(... self.ny_pixels, self.nz_pixels, self.pixel_angle, self.pixel_height, ... - self.dist_to_detector, self.to_source_vec, self.source_position... + self.dist_to_detector, self.to_source_vec, self.source_position, voxels... ); end - function pixel_generator = get_pixel_generator(self, angle_index, ray_per_pixel) + function pixel_generator = get_pixel_generator(self, angle_index, voxels, ray_per_pixel) % Create a function which returns the rays which should be fired to hit each pixel. % Only 1 ray per pixel is supported at the moment, as anti-aliasing techniques are not yet implemented. arguments self curved_detector angle_index double + voxels voxel_array ray_per_pixel int32 = 1 end - assert(nargin==2, "Only 1 ray per pixel is supported at the moment, as anti-aliasing techniques are not yet implemented.") + assert(nargin==3, "Only 1 ray per pixel is supported at the moment, as anti-aliasing techniques are not yet implemented.") if angle_index == 1; rot_mat = eye(3); else ; rot_mat = rotz(self.rot_angle * (angle_index - 1)); @@ -87,11 +89,10 @@ function reset(self) pixel_generator = @generator; static_ray_generator = curved_detector.get_ray_generator_static(... self.ny_pixels, self.nz_pixels, self.pixel_angle, self.pixel_height, ... - self.dist_to_detector, to_source_vec, source_pos... + self.dist_to_detector, to_source_vec, source_pos, voxels... ); - function pixel_value = generator(y_pixel, z_pixel, voxels) - xray = static_ray_generator(y_pixel, z_pixel); - pixel_value = xray.calculate_mu(voxels); + function pixel_value = generator(y_pixel, z_pixel) + pixel_value = static_ray_generator(y_pixel, z_pixel).calculate_mu(); end end diff --git a/src/detectors/detector.m b/src/detectors/detector.m index 615fda8..5919c29 100644 --- a/src/detectors/detector.m +++ b/src/detectors/detector.m @@ -24,8 +24,8 @@ rotate(self) reset(self) hit_pixel(self, ray) - get_ray_generator(self, ray_type, ray_per_pixel) - get_pixel_generator(self, angle_index, ray_type, ray_per_pixel) + get_ray_generator(self, voxels, ray_type, ray_per_pixel) + get_pixel_generator(self, angle_index, voxels, ray_type, ray_per_pixel) end methods @@ -87,12 +87,12 @@ function check_voxels(self, voxels) self.check_voxels(voxels); image = zeros(self.ny_pixels, self.nz_pixels, self.num_rotations); for i = 1:self.num_rotations - ray_generator = self.get_ray_generator(); + ray_generator = self.get_ray_generator(voxels); for j = 1:self.nz_pixels for k = 1:self.ny_pixels ray = ray_generator(k, j); - mu = ray.calculate_mu(voxels); + mu = ray.calculate_mu(); image(k, j, i) = mu; end @@ -111,10 +111,10 @@ function check_voxels(self, voxels) image = zeros(self.ny_pixels, self.nz_pixels, self.num_rotations); get_pixel_generator = @self.get_pixel_generator; for k = 1:self.num_rotations - pixel_calc = get_pixel_generator(k); + pixel_calc = get_pixel_generator(k, voxels); for j = 1:self.nz_pixels parfor i = 1:self.ny_pixels - image(i, j, k) = feval(pixel_calc, i, j, voxels); + image(i, j, k) = feval(pixel_calc, i, j); end end end @@ -131,7 +131,7 @@ function check_voxels(self, voxels) elseif self.scatter_type == 1 % Fast scatter scatter = self.conv_scatter(image); elseif self.scatter_type == 2 % Slow scatter - scatter = self.slow_scatter_p(voxels); + scatter = (self.slow_scatter_p(voxels)-image) / 2; end end @@ -141,12 +141,12 @@ function check_voxels(self, voxels) air = voxel_object(@(i,j,k) i==i, material_attenuation("air")); array = voxel_array(zeros(3, 1), zeros(3,1)+1e6, dtd/10, air); scan = zeros(self.ny_pixels, self.nz_pixels, self.num_rotations); - ray_generator = self.get_ray_generator(); + ray_generator = self.get_ray_generator(array); image_at_angle = zeros(self.ny_pixels, self.nz_pixels); for k = 1:self.nz_pixels for j = 1:self.ny_pixels ray = ray_generator(j, k); - mu = ray.calculate_mu(array); + mu = ray.calculate_mu(); image_at_angle(j, k) = mu; end end @@ -178,26 +178,25 @@ function check_voxels(self, voxels) % Do some Monte Carlo simulation of scatter ny = self.ny_pixels; nz = self.nz_pixels; scatter = zeros(ny, nz, self.num_rotations); - for sample = 1:self.scatter_factor - for k = 1:self.num_rotations - pixel_calc = self.get_pixel_generator(k, @scatter_ray); - scatter_idxs = zeros(ny, nz, 2); scatter_vals = zeros(ny, nz); - for j = 1:nz - parfor i = 1:ny - [pval, pixel, ~] = feval(pixel_calc, i, j, voxels); - if pval < inf % Collect the scatter values for adding later - scatter_idxs(i, j, :) = pixel; - scatter_vals(i, j) = pval; - end - end + for k = 1:self.num_rotations + pixel_calc = self.get_pixel_generator(k, voxels, @scatter_ray); + scatter_idxs = zeros(ny, nz, self.scatter_factor, 2); + scatter_vals = zeros(ny, nz, self.scatter_factor); + for j = 1:nz + parfor i = 1:ny + [pval, pixel, ~] = feval(pixel_calc, i, j); + scatter_idxs(i, j, :, :) = pixel; + scatter_vals(i, j, :) = pval; end + end - % Add the scatter to the image + % Add the scatter to the image + for sf = 1:self.scatter_factor for j = 1:nz for i = 1:ny - if scatter_vals(i, j) > 0 - scatter(scatter_idxs(i, j, 1), scatter_idxs(i, j, 2), k) = ... - scatter(scatter_idxs(i, j, 1), scatter_idxs(i, j, 2), k) + scatter_vals(i, j); + if ~isnan(scatter_vals(i, j, sf)) + scatter(scatter_idxs(i, j, sf, 1), scatter_idxs(i, j, sf, 2), k) = ... + scatter(scatter_idxs(i, j, sf, 1), scatter_idxs(i, j, sf, 2), k) + scatter_vals(i, j, sf); end end end @@ -211,14 +210,15 @@ function check_voxels(self, voxels) self.reset(); % Reset the detector to the initial position ny = self.ny_pixels; nz = self.nz_pixels; scatter = zeros(ny, nz, self.num_rotations); - for sample = 1:self.scatter_factor - for k = 1:self.num_rotations - pixel_calc = self.get_pixel_generator(k, @scatter_ray); - for j = 1:nz - for i = 1:ny - [pval, pixel, ~] = pixel_calc(i, j, voxels); - if pval < inf % Collect the attenuation values - scatter(pixel(1), pixel(2), k) = scatter(pixel(1), pixel(2), k) + pval; + for k = 1:self.num_rotations + pixel_calc = self.get_pixel_generator(k, voxels, @scatter_ray); + for j = 1:nz + for i = 1:ny + [pval, pixel, ~] = pixel_calc(i, j); + for sf = 1:self.scatter_factor + if ~isnan(pval(sf)) + scatter(pixel(sf, 1), pixel(sf, 2), k) = ... + scatter(pixel(sf, 1), pixel(sf, 2), k) + pval(sf); end end end diff --git a/src/detectors/parallel_detector.m b/src/detectors/parallel_detector.m index 5c13fb5..ba0a2b2 100644 --- a/src/detectors/parallel_detector.m +++ b/src/detectors/parallel_detector.m @@ -8,7 +8,7 @@ detector_vec (3, 1) double = [1;0;0] % Vector from left to right corner of detector end methods (Access=private, Static) - function generator = get_ray_generator_static(ray, ny_pixels, nz_pixels, centre, detector_vec, pixel_width, pixel_height, dist_to_detector, to_source_vec) + function generator = get_ray_generator_static(ray, ny_pixels, nz_pixels, centre, detector_vec, pixel_width, pixel_height, dist_to_detector, to_source_vec, voxels) generator = @get_ray_attrs; % Create the function which returns the rays function [xray] = get_ray_attrs(y_pixel, z_pixel) assert(y_pixel <= ny_pixels && y_pixel > 0 && z_pixel <= nz_pixels && z_pixel > 0, ... @@ -17,7 +17,7 @@ detector_vec .* (y_pixel - (ny_pixels+1)/2) .* pixel_width + ... [0;0;pixel_height] .* (z_pixel - (nz_pixels+1)/2); source_position = pixel_centre + to_source_vec .* dist_to_detector; - xray = ray(source_position, -to_source_vec, dist_to_detector); + xray = ray(source_position, -to_source_vec, dist_to_detector, voxels); end end end @@ -54,31 +54,34 @@ function reset(self) self.centre = self.init_centre; end - function ray_generator = get_ray_generator(self, ray_type, ray_per_pixel) + function ray_generator = get_ray_generator(self, voxels, ray_type, ray_per_pixel) % Create a function which returns the rays which should be fired to hit each pixel. % Only 1 ray per pixel is supported at the moment, as anti-aliasing techniques are not yet implemented. arguments self parallel_detector + voxels voxel_array ray_type = @ray ray_per_pixel int32 = 1 end - assert(nargin<3, "Only 1 ray per pixel is supported at the moment, as anti-aliasing techniques are not yet implemented.") + assert(nargin<4, "Only 1 ray per pixel is supported at the moment, as anti-aliasing techniques are not yet implemented.") % Create the function which returns the rays ray_generator = parallel_detector.get_ray_generator_static(... ray_type, self.ny_pixels, self.nz_pixels, self.centre, self.detector_vec, ... - self.pixel_dims(1), self.pixel_dims(2), self.dist_to_detector, self.to_source_vec... + self.pixel_dims(1), self.pixel_dims(2), self.dist_to_detector, ... + self.to_source_vec, voxels... ); end - function pixel_generator = get_pixel_generator(self, angle_index, ray_type, ray_per_pixel) + function pixel_generator = get_pixel_generator(self, angle_index, voxels, ray_type, ray_per_pixel) % Create a function which returns the rays which should be fired to hit each pixel. % Only 1 ray per pixel is supported at the moment, as anti-aliasing techniques are not yet implemented. arguments self parallel_detector angle_index double + voxels voxel_array ray_type = @ray ray_per_pixel int32 = 1 end - assert(nargin<4, "Only 1 ray per pixel is supported at the moment, as anti-aliasing techniques are not yet implemented.") + assert(nargin<5, "Only 1 ray per pixel is supported at the moment, as anti-aliasing techniques are not yet implemented.") if angle_index == 1; rot_mat = eye(3); else ; rot_mat = rotz(self.rot_angle * (angle_index - 1)); @@ -90,7 +93,8 @@ function reset(self) % Create the function which returns the rays static_ray_generator = parallel_detector.get_ray_generator_static(... ray_type, self.ny_pixels, self.nz_pixels, current_c, current_dv, ... - self.pixel_dims(1), self.pixel_dims(2), self.dist_to_detector, current_sv... + self.pixel_dims(1), self.pixel_dims(2), self.dist_to_detector, ... + current_sv, voxels... ); ray_type_str = func2str(ray_type); % Get the name of the ray type @@ -99,24 +103,32 @@ function reset(self) else; error('parallel_detector:InvalidRayType', "Must be either 'ray' or 'scatter_ray'."); end - function pixel_value = generator(y_pixel, z_pixel, voxels) + function pixel_value = generator(y_pixel, z_pixel) xray = static_ray_generator(y_pixel, z_pixel); - pixel_value = xray.calculate_mu(voxels); + pixel_value = xray.calculate_mu(); end - function [pixel_value, pixel, scattered] = scatter_generator(y_pixel, z_pixel, voxels) + function [pixel_values, pixels, scattered] = scatter_generator(y_pixel, z_pixel) xray = static_ray_generator(y_pixel, z_pixel); - xray = xray.calculate_mu(voxels); - - mu = xray.mu; scattered = xray.scatter_event > 0; - - hit = true; - if scattered; [pixel, hit] = self.hit_pixel(xray, current_dv); - else; pixel = [y_pixel, z_pixel]; - end + pixel_values = zeros(self.scatter_factor, 1); + pixels = zeros(self.scatter_factor, 2); + scattered = zeros(self.scatter_factor, 1); - if hit; pixel_value = mu; - else; pixel_value = inf; + for i = 1:self.scatter_factor + new_ray = xray.calculate_mu(); + + this_scatter = new_ray.scatter_event > 0; + + hit = true; + if this_scatter; [pixel, hit] = self.hit_pixel(new_ray, current_dv); + else; pixel = [y_pixel, z_pixel]; + end + pixels(i, :) = pixel; + scattered(i) = this_scatter; + if hit; pixel_values(i) = new_ray.mu; + else; pixel_values(i) = NaN; + end + % xray = xray.randomise_n_mfp(); end end end diff --git a/src/ray_tracing/ray.m b/src/ray_tracing/ray.m index 52ad9e6..ba848a2 100644 --- a/src/ray_tracing/ray.m +++ b/src/ray_tracing/ray.m @@ -3,24 +3,36 @@ start_point (3, 1) double % 3D point v1_to_v2 (3, 1) double % Vector from start_point to end_point energy double % energy of the ray in KeV + + voxels voxel_array % The voxel array the ray will be traced through + mu_dict % Dictionary of mu values for each material + lengths % Length of the ray in each voxel + indices % Indices of the voxels the ray intersects end properties (Access=private, Constant) - use_mex = ~~exist('ray_trace_mex', 'file'); % Use the MEX implementation of the photon_attenuation package if available + use_mex = ~~exist('ray_trace_mex', 'file'); % Use the MEX implementation of the ray tracing if available end methods - function self = ray(start_point, direction, dist_to_detector, energy) + function self = ray(start_point, direction, dist_to_detector, voxels, energy) arguments start_point (3, 1) double direction (3, 1) double dist_to_detector double + voxels voxel_array energy double = 30 %KeV end self.start_point = start_point; self.v1_to_v2 = direction .* dist_to_detector; self.energy = energy; + + % Using the assumption that the ray's energy is constant + % Create a dictionary of the values of mu for each material + self.mu_dict = voxels.get_mu_dict(energy); + [self.lengths, self.indices] = self.get_intersections(voxels); + self.voxels = voxels; end function [lengths, indices] = get_intersections(self, voxels) @@ -41,21 +53,15 @@ end end - function mu = calculate_mu (self, voxels) - arguments - self ray - voxels voxel_array - end - % Using the assumption that the ray's energy is constant - % Create a dictionary of the values of mu for each material - mu_dict = voxels.get_mu_dict(self.energy); - + function mu = calculate_mu (self) + % I don't see the point of this function if we precalculate everything else + % Maybe for testing and consistency with the scattered ray class + % Calculate the mu of the ray - [lengths, indices] = self.get_intersections(voxels); - if isempty(lengths) % No intersections + if isempty(self.lengths) % No intersections mu = 0; else - mu = sum(lengths .* voxels.get_saved_mu(indices, mu_dict)); + mu = sum(self.lengths .* self.voxels.get_saved_mu(self.indices, self.mu_dict)); end end end diff --git a/src/scatter/scatter_ray.m b/src/scatter/scatter_ray.m index ed40899..7ff27f7 100644 --- a/src/scatter/scatter_ray.m +++ b/src/scatter/scatter_ray.m @@ -5,65 +5,56 @@ scatter_event double = 0 % A number indicating the current scatter event direction (3, 1) double % 3D direction end_point (3, 1) double % 3D point + mfp_dict % Dictionary of mfp values for each material end - + methods - function self = scatter_ray(start_point, direction, dist_to_detector, energy) + function self = scatter_ray(start_point, direction, dist_to_detector, voxels, energy) arguments start_point double direction double dist_to_detector double + voxels voxel_array energy double = 30 % The energy of the ray in keV end - self@ray(start_point, direction, dist_to_detector, energy); + self@ray(start_point, direction, dist_to_detector, voxels, energy); self.direction = direction; self.end_point = start_point + self.v1_to_v2; + self.mfp_dict = voxels.get_mfp_dict(energy); self.n_mfp = -log(rand); % Control with rng(seed) for reproducibility end - function self = calculate_mu (self, voxels) % I have broken symmetry with the ray class -> maybe always return ray? - arguments - self ray - voxels voxel_array - end - % Get material mu and mean free path dictionaries - mu_dict = voxels.get_mu_dict(self.energy); - mfp_dict = voxels.get_mfp_dict(self.energy); - - % Calculate the mu of the ray - [ls, idxs] = self.get_intersections(voxels); - if isempty(ls); return; end % If there are no intersections, exit + function self = calculate_mu (self) % I have broken symmetry with the ray class -> maybe always return ray? - i = 1; % Initialize the index of the intersection - mfps = voxels.get_saved_mfp(idxs, mfp_dict); % Get the mean free path of the first intersection + if isempty(self.lengths); return; end % If there are no intersections, exit + ls = self.lengths; idxs = self.indices; % Get the lengths and indices of the intersections + mfps = self.voxels.get_saved_mfp(idxs, self.mfp_dict); % Get the mean free path of the first intersection for i = 1:length(ls) self.n_mfp = self.n_mfp - ls(i) / mfps(i); % Update the number of mean free paths until the next scatter event if self.n_mfp < 0 % If the number of mean free paths is less than 0, scatter the ray % Calculate the mu of the ray until the end of the current voxel - self.mu = self.mu + sum(ls(1:i) .* voxels.get_saved_mu(idxs(:, 1:i), mu_dict)); + self.mu = self.mu + sum(ls(1:i) .* self.voxels.get_saved_mu(idxs(:, 1:i), self.mu_dict)); % Remove the mu of the current voxel up to the scatter event - self.mu = self.mu + mfps(i) * self.n_mfp * voxels.get_saved_mu(idxs(i), mu_dict); + self.mu = self.mu + mfps(i) * self.n_mfp * self.voxels.get_saved_mu(idxs(i), self.mu_dict); % Get the new direction and energy of the ray, and update the start point - [direction, energy] = self.scatter(); + [ndirection, energy] = self.scatter(); start_point = self.start_point + (sum(ls(1:i)) + self.n_mfp * mfps(i)) .* self.direction; % Create a new ray with the new direction, energy, and start point - new_ray = scatter_ray(start_point, direction, sum(ls), energy); + new_ray = scatter_ray(start_point, ndirection, sum(ls), self.voxels, energy); new_ray.mu = self.mu; % Set the mu of the new ray to the mu of the old ray new_ray.scatter_event = self.scatter_event + 1; % Update the scatter event of the new ray - self = new_ray.calculate_mu(voxels); % Get the mu of the new ray + self = new_ray.calculate_mu(); % Get the mu of the new ray return % Return the mu of the new ray end - - i = i + 1; % Move to the next intersection end - self.mu = self.mu + sum(ls .* voxels.get_saved_mu(idxs, mu_dict)); % This case only occurs if the ray does not scatter + self.mu = self.mu + sum(ls .* self.voxels.get_saved_mu(idxs, self.mu_dict)); % This case only occurs if the ray does not scatter end function [direction, energy] = scatter(self) % Not sure if this needs to be separate? @@ -101,6 +92,10 @@ energy = ((constants.em_ee * E_0) / (constants.em_ee + E_0 * (1 - cos_theta))); end + + function self = randomise_n_mfp(self) + self.n_mfp = -log(rand); + end end end diff --git a/tests/detector_tests.m b/tests/detector_tests.m index 29dc1c9..1184563 100644 --- a/tests/detector_tests.m +++ b/tests/detector_tests.m @@ -81,14 +81,15 @@ function test_parallel_init(tc) function test_curved_ray_gen(tc) detector = curved_detector(9, [9*pi/60, 0.4], [60, 10], 4); - ray_generator = detector.get_ray_generator(); + empty_voxels = voxel_array([0;0;0], [1;1;1], 1); + ray_generator = detector.get_ray_generator(empty_voxels); rot_by_pixel = rotz(pi/60); unit_vector = rotz(pi/120) * [-1; 0; 0]; z_pos = @(i) (-2 + 0.2 + (0.4 .* (i-1)))/9; for i = 30:40 for j = 1:10 - my_ray = ray([0; 4.5; 0], rot_by_pixel^(i-1) * unit_vector, 9); + my_ray = ray([0; 4.5; 0], rot_by_pixel^(i-1) * unit_vector, 9, empty_voxels); gen_ray = ray_generator(i, j); tc.verifyEqual(gen_ray.start_point, my_ray.start_point); tc.verifyEqual(gen_ray.v1_to_v2, my_ray.v1_to_v2 + [0;0;z_pos(j)*9], 'RelTol', 1e-13); @@ -99,14 +100,15 @@ function test_curved_ray_gen(tc) function test_para_ray_gen(tc) unit_vector = [0; -1; 0]; detector = parallel_detector(2, [0.1, 0.35], [110, 20], 1); - ray_generator = detector.get_ray_generator(); + empty_voxels = voxel_array([0;0;0], [1;1;1], 1); + ray_generator = detector.get_ray_generator(empty_voxels); y_increment = [1; 0; 0] * 0.1; z_increment = [0; 0; 7] / 20; start = [-5.45; 1; -3.325]; for i = 50:60 for j = 5:15 - my_ray = ray(start + y_increment*(i-1) + z_increment*(j-1), unit_vector, 2); + my_ray = ray(start + y_increment*(i-1) + z_increment*(j-1), unit_vector, 2, empty_voxels); gen_ray = ray_generator(i, j); tc.verifyEqual(gen_ray.start_point, my_ray.start_point, 'AbsTol', 7e-15); tc.verifyEqual(gen_ray.v1_to_v2, my_ray.v1_to_v2); @@ -116,34 +118,34 @@ function test_para_ray_gen(tc) function test_para_pixel_gen_scatter(tc) detector = parallel_detector(2, [0.01, 0.01], [110, 20], 1); - pix_gen_scatter = detector.get_pixel_generator(1, @scatter_ray); - pix_gen = detector.get_pixel_generator(1); graphite = material_attenuation("graphite", 12, 1, 2.26); array = voxel_array(zeros(3, 1), [5; 5; 5], 1, voxel_object(@(i, j, k) i==i, graphite)); + pix_gen_scatter = detector.get_pixel_generator(1, array, @scatter_ray); + pix_gen = detector.get_pixel_generator(1, array); for i = 50:60 for j = 5:15 - [pixel_value, pixel_hit, scattered] = pix_gen_scatter(i, j, array); - [pixel_value_2] = pix_gen(i, j, array); + [scatter_pval, pixel_hit, scattered] = pix_gen_scatter(i, j); + [pval] = pix_gen(i, j); if scattered && any(pixel_hit) && all([i, j] == pixel_hit) % Unlikely case: Scattered, but small angle so still hits pixel - tc.verifyEqual(pixel_value, pixel_value_2, 'RelTol', 1e-5); + tc.verifyEqual(scatter_pval, pval, 'RelTol', 1e-5); elseif scattered && any(pixel_hit) % Scattered, but still hits pixel - tc.verifyNotEqual(pixel_value, pixel_hit); + tc.verifyNotEqual(scatter_pval, pixel_hit); elseif scattered % Scattered, but doesn't hit pixel - tc.verifyEqual(pixel_value, inf); + tc.verifyTrue(isnan(scatter_pval)) else % Doesn't scatter - tc.verifyEqual(pixel_value, pixel_value_2, 'RelTol', 1e-5); + tc.verifyEqual(scatter_pval, pval, 'RelTol', 1e-5); end end end func = @(a, b, c) "Not a ray"; - tc.verifyError(@() detector.get_pixel_generator(1, func), "parallel_detector:InvalidRayType"); + tc.verifyError(@() detector.get_pixel_generator(1, array, func), "parallel_detector:InvalidRayType"); end function test_get_scan_angles(tc) @@ -207,7 +209,8 @@ function test_curved_generate_image_p(tc) function test_hit_para_pixel(tc) detector = parallel_detector(2, [0.1, 0.35], [110, 20], 10); - ray_generator = detector.get_ray_generator(); + empty_voxels = voxel_array([0;0;0], [1;1;1], 1); + ray_generator = detector.get_ray_generator(empty_voxels); for i = 50:60 for j = 5:15 gen_ray = ray_generator(i, j); @@ -220,18 +223,18 @@ function test_hit_para_pixel(tc) pixel_centre = detector.centre + detector.detector_vec .* ... (-(111)/2) .* 0.1 + [0;0;0.35] .* (-(21)/2); source_position = pixel_centre + detector.to_source_vec .* 2; - xray = ray(source_position, -detector.to_source_vec, 2); % should miss + xray = ray(source_position, -detector.to_source_vec, 2, empty_voxels); % should miss [pixel, hit] = detector.hit_pixel(xray, detector.detector_vec); tc.verifyEqual(pixel, [0, 0]); tc.verifyFalse(hit); - xray = ray(source_position, rotz(pi/2) * -detector.to_source_vec, 2); % should not intersect + xray = ray(source_position, rotz(pi/2) * -detector.to_source_vec, 2, empty_voxels); % should not intersect [pixel, hit] = detector.hit_pixel(xray, detector.detector_vec); tc.verifyEqual(pixel, [0, 0]); tc.verifyFalse(hit); detector.rotate(); - ray_generator = detector.get_ray_generator(); + ray_generator = detector.get_ray_generator(empty_voxels); for i = 50:60 for j = 5:15 gen_ray = ray_generator(i, j); @@ -258,11 +261,11 @@ function test_scatter_image(tc) detector = parallel_detector(10, [1, 1], [5, 1], 4, "slow"); scatter_image = detector.generate_image_p(array); - tc.verifyEqual(image.*2, scatter_image, 'RelTol', 1e-15, 'AbsTol', 1e-15); + tc.verifyEqual(image, scatter_image, 'RelTol', 1e-15, 'AbsTol', 1e-15); detector = parallel_detector(10, [1, 1], [5, 1], 4, "slow", 2); scatter_image = detector.generate_image_p(array); - tc.verifyEqual(image.*2, scatter_image, 'RelTol', 1e-15, 'AbsTol', 1e-15); + tc.verifyEqual(image, scatter_image, 'RelTol', 1e-15, 'AbsTol', 1e-15); scatter = detector.slow_scatter(array); % Now scatter factor is 2 tc.verifyEqual(image, scatter, 'RelTol', 1e-15, 'AbsTol', 1e-15); diff --git a/tests/ray_tests.m b/tests/ray_tests.m index ca3b1dc..0ad43d1 100644 --- a/tests/ray_tests.m +++ b/tests/ray_tests.m @@ -1,4 +1,12 @@ classdef ray_tests < matlab.unittest.TestCase + properties + air + lead + water + water_array + lead_array + end + methods (TestClassSetup) % Shared setup for the entire test class @@ -6,48 +14,54 @@ methods (TestMethodSetup) % Setup for each test + function setup_test(tc) + tc.air = material_attenuation("air"); + tc.lead = material_attenuation("lead", 82, 1, 11.34); + tc.water = material_attenuation("water"); + abox = voxel_cube([0;0;0], [3;3;3], tc.water); + tc.water_array = voxel_array(zeros(3, 1), [5; 5; 5], 1, abox); + abox = voxel_cube([0;0;0], [3;3;3], tc.lead); + tc.lead_array = voxel_array(zeros(3, 1), [5; 5; 5], 1, abox); + end end methods function gen_get_intersections(tc, ray) - water = material_attenuation("water"); - my_box = voxel_cube([0;0;0], [3;3;3], water); - array = voxel_array(zeros(3, 1), [5; 5; 5], 1, my_box); threes = zeros(3, 5) + 3; - r = ray([-6;0;0], [1;0;0], 12); - [lengths, indices] = r.get_intersections(array); + r = ray([-6;0;0], [1;0;0], 12, tc.water_array); + [lengths, indices] = r.get_intersections(tc.water_array); exp = threes; exp(1, :) = [1, 2, 3, 4, 5]; tc.assertEqual(lengths, [1, 1, 1, 1, 1], 'AbsTol', 2e-15); tc.assertEqual(indices, exp); - r = ray([0;-6;0], [0;1;0], 12); - [lengths, indices] = r.get_intersections(array); + r = ray([0;-6;0], [0;1;0], 12, tc.water_array); + [lengths, indices] = r.get_intersections(tc.water_array); exp = threes; exp(2, :) = [1, 2, 3, 4, 5]; tc.assertEqual(lengths, [1, 1, 1, 1, 1], 'AbsTol', 2e-15); tc.assertEqual(indices, exp); - r = ray([0;0;-6], [0;0;1], 12); - [lengths, indices] = r.get_intersections(array); + r = ray([0;0;-6], [0;0;1], 12, tc.water_array); + [lengths, indices] = r.get_intersections(tc.water_array); exp = threes; exp(3, :) = [1, 2, 3, 4, 5]; tc.assertEqual(lengths, [1, 1, 1, 1, 1], 'AbsTol', 2e-15); tc.assertEqual(indices, exp); - r = ray([-6;-6;0], [1;1;0], 20); - [lengths, indices] = r.get_intersections(array); + r = ray([-6;-6;0], [1;1;0], 20, tc.water_array); + [lengths, indices] = r.get_intersections(tc.water_array); exp = threes; exp(2, :) = [1, 2, 3, 4, 5]; exp(1, :) = [1, 2, 3, 4, 5]; s2 = sqrt(2); tc.assertEqual(lengths, [s2, s2, s2, s2, s2], 'AbsTol', 2e-15); tc.assertEqual(indices, exp); - r = ray([0;0;-6], [0;0;1], 12); - [lengths, indices] = r.get_intersections(array); + r = ray([0;0;-6], [0;0;1], 12, tc.water_array); + lengths = r.lengths; indices = r.indices; exp = threes; exp(3, :) = [1, 2, 3, 4, 5]; tc.assertEqual(lengths, [1, 1, 1, 1, 1], 'AbsTol', 2e-15); tc.assertEqual(indices, exp); - r = ray([6;6;6], [-1;-1;-1], 22); % 3D diagonal backwards - [lengths, indices] = r.get_intersections(array); + r = ray([6;6;6], [-1;-1;-1], 22, tc.water_array); % 3D diagonal backwards + [lengths, indices] = r.get_intersections(tc.water_array); exp = [ 5, 4, 3, 2, 1; 5, 4, 3, 2, 1; @@ -57,8 +71,8 @@ function gen_get_intersections(tc, ray) tc.assertEqual(lengths, [s3, s3, s3, s3, s3], 'AbsTol', 5e-15); tc.assertEqual(indices, exp); - r = ray([6;6;6], [1;1;1], 22); % 3D diagonal away - [lengths, indices] = r.get_intersections(array); + r = ray([6;6;6], [1;1;1], 22, tc.water_array); % 3D diagonal away + [lengths, indices] = r.get_intersections(tc.water_array); tc.assertEqual(lengths, []); tc.assertEqual(indices, []); end @@ -67,12 +81,12 @@ function gen_get_intersections(tc, ray) methods (Test) % Test methods function test_ray_init(tc) - r = ray([0;0;0], [1;0;0], 10); + r = ray([0;0;0], [1;0;0], 10, tc.water_array); tc.assertEqual(r.start_point, [0;0;0]); tc.assertEqual(r.v1_to_v2, [10;0;0]); tc.assertEqual(r.energy, 30); % default energy - may change later - r2 = ray([0;-5;0], [0;1;0], 10, 100); + r2 = ray([0;-5;0], [0;1;0], 10, tc.water_array, 100); tc.assertEqual(r2.start_point, [0;-5;0]); tc.assertEqual(r2.v1_to_v2, [0;10;0]); tc.assertEqual(r2.energy, 100); @@ -83,29 +97,23 @@ function test_siddon_ray(tc) end function test_calculate_mu(tc) - water = material_attenuation("water"); - my_box = voxel_cube([0;0;0], [3;3;3], water); - array = voxel_array(zeros(3, 1), [5; 5; 5], 1, my_box); - water_attenuation = water.get_mu(30); % Default value of ray - % voxel_array.air - % array.air - % air= material_attenuation("air"); - % air_attenuation = air.get_mu(30); + water_attenuation = tc.water.get_mu(30); % Default value of ray + lead_attenuation = tc.lead.get_mu(30); % Default value of ray - r = ray([-6;0;0], [1;0;0], 12); - tc.assertEqual(r.calculate_mu(array), 3*water_attenuation, "RelTol", 3e-16); + r = ray([-6;0;0], [1;0;0], 12, tc.water_array); + tc.assertEqual(r.calculate_mu(), 3*water_attenuation, "RelTol", 3e-16); - r = ray([0;-6;0], [0;1;0], 12); - tc.assertEqual(r.calculate_mu(array), 3*water_attenuation, "RelTol", 3e-16); + r = ray([0;-6;0], [0;1;0], 12, tc.water_array); + tc.assertEqual(r.calculate_mu(), 3*water_attenuation, "RelTol", 3e-16); - r = ray([0;0;-6], [0;0;1], 12); - tc.assertEqual(r.calculate_mu(array), 3*water_attenuation, "RelTol", 3e-16); + r = ray([0;0;-6], [0;0;1], 12, tc.water_array); + tc.assertEqual(r.calculate_mu(), 3*water_attenuation, "RelTol", 3e-16); - r = ray([6;6;6], [-1;-1;-1], 22); % 3D diagonal backwards - tc.assertEqual(r.calculate_mu(array), 3*sqrt(3)*water_attenuation, "RelTol", 4e-16); + r = ray([6;6;6], [-1;-1;-1], 22, tc.lead_array); % 3D diagonal backwards + tc.assertEqual(r.calculate_mu(), 3*sqrt(3)*lead_attenuation, "RelTol", 4e-16); - r = ray([6;6;6], [1;1;1], 22); % 3D diagonal away - tc.assertEqual(r.calculate_mu(array), 0); + r = ray([6;6;6], [1;1;1], 22, tc.water_array); % 3D diagonal away + tc.assertEqual(r.calculate_mu(), 0); end function test_scatter_ray(tc) @@ -114,27 +122,23 @@ function test_scatter_ray(tc) start = [-6;0;0]; direction = [1;0;0]; dist_to_detector = 100; - r = ray(start, direction, dist_to_detector, energy); - - lead = material_attenuation("lead", 82, 1, 11.34); % Dense -> Likely to scatter - my_box = voxel_cube([0;0;0], [3;3;3], lead); - array = voxel_array(zeros(3, 1), [5; 5; 5], 1, my_box); + r = ray(start, direction, dist_to_detector, tc.lead_array, energy); - mu = r.calculate_mu(array); - [lengths, indices] = r.get_intersections(array); + mu = r.calculate_mu(); + [lengths, indices] = r.get_intersections(tc.lead_array); for i = 1:100 - rs = scatter_ray(start, direction, dist_to_detector, energy); % In loop so n_mfp is random - nrs = rs.calculate_mu(array); + rs = scatter_ray(start, direction, dist_to_detector, tc.lead_array, energy); % In loop so n_mfp is random + nrs = rs.calculate_mu(); total_mu = nrs.mu; if nrs.scatter_event == 1 % Check if it scattered once - mfp_dict = array.get_mfp_dict(energy); - mfp = array.get_saved_mfp(indices, mfp_dict); + mfp_dict = tc.lead_array.get_mfp_dict(energy); + mfp = tc.lead_array.get_saved_mfp(indices, mfp_dict); - mu_dict = array.get_mu_dict(energy); + mu_dict = tc.lead_array.get_mu_dict(energy); ray_n_mfp = lengths ./ mfp; - ray_mu = lengths .* array.get_saved_mu(indices, mu_dict); + ray_mu = lengths .* tc.lead_array.get_saved_mu(indices, mu_dict); tc.assertTrue(sum(ray_n_mfp) > rs.n_mfp); % Check it should have scattered tc.assertTrue(nrs.energy < energy); % energy is lost in scatter @@ -150,8 +154,11 @@ function test_scatter_ray(tc) end tc.assertEqual(nrs.start_point, start + scatter_length .* direction, 'RelTol', 1e-14, 'AbsTol', 1e-15); - scattered_ray = ray(nrs.start_point, nrs.direction, dist_to_detector, nrs.energy); - scatter_mu = scattered_ray.calculate_mu(array); + scattered_ray = ray(nrs.start_point, nrs.direction, dist_to_detector, tc.lead_array, nrs.energy); + scatter_mu = scattered_ray.calculate_mu(); + uptoscatter_mu + scatter_mu + total_mu tc.assertEqual(total_mu, scatter_mu + uptoscatter_mu, 'RelTol', 1e-14, 'AbsTol', 1e-14); tc.assertNotEqual(nrs.n_mfp, rs.n_mfp); % n_mfp should have changed @@ -164,8 +171,8 @@ function test_scatter_ray(tc) tc.assertEqual(nrs.end_point, start + dist_to_detector * direction); end end - r = ray([0;0;0], [1;0;0], 10); - rs = scatter_ray([0;0;0], [1;0;0], 10); + r = ray([0;0;0], [1;0;0], 10, tc.lead_array); + rs = scatter_ray([0;0;0], [1;0;0], 10, tc.lead_array); tc.assertEqual(r.energy, rs.energy); % Check the default energy is the same end diff --git a/tests/scatter_tests.m b/tests/scatter_tests.m index def460c..249a489 100644 --- a/tests/scatter_tests.m +++ b/tests/scatter_tests.m @@ -76,13 +76,13 @@ function test_cross_section(tc) function test_angle_distribution(tc) % Test the energy distribution of the scattered photons - energies = [10, 30, 60, 100, 300, 600, 1000]; + energies = [10, 30, 60, 100, 300]; vectors = [ [1 0 0] [0 1 0] [0 0 1] ]; - + empty_voxels = voxel_array([0;0;0], [1;1;1], 1); for ei = 1:length(energies)-2 for vi = 1:height(vectors) e1 = energies(ei); @@ -94,8 +94,8 @@ function test_angle_distribution(tc) energy_dist1 = zeros(1, num_rays); energy_dist2 = zeros(1, num_rays); for i = 1:num_rays - x_ray1 = scatter_ray([0; 0; 0], v, 1, e1); - x_ray2 = scatter_ray([0; 0; 0], v, 1, e2); + x_ray1 = scatter_ray([0; 0; 0], v, 1, empty_voxels, e1); + x_ray2 = scatter_ray([0; 0; 0], v, 1, empty_voxels, e2); [d1, e1s] = x_ray1.scatter(); [d2, e2s] = x_ray2.scatter();