Mc scatter #50
Merged
Mc scatter #50
Conversation
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
I like Pamela, Angela, Sandra, and RITA... ...for sampling Rayleigh scatter angles. The RITA that I am implementing is a generalized RITA that is not yet tailored to the specific PDF and CDF of Rayleigh scatter. I will have to write those functions and pass them into my generalized RITA code.
To apply it to Rayleigh scatter, I will need to define the true PDF and CDF as functions to the RITA algorithm
Turned the RITA parameter arrays into a class, rita.py:RITA. This allows for a dictionary, Dict[str, RITA] that maps the material (e.g. "bone", "soft tissue") to the RITA sampler. Based on my reading of the PENELOPE-2006 specifications, the distribution we are sampling, \pi(x^2), is distributed like [F(x,Z)]^2. However, for a given photon with dimensionless energy \kappa, we only want to sample values on the interval (0, (20.6074 * 2\kappa)^2). So, it appears that a RITA object is a function of both the material AND the energy of the incoming photon, since a RITA object is specified by the PDF to approximate as well as the interval to approximate over. This dependency on material AND energy is bad, since we don't want to recalculate a RITA object for every single photon energy level within a material. The way around this need for recalculation is to rely on resampling. We initialize a sinlge RITA object for each material, based on the maximum possible energy of a photon from the spectrum that we are using. Then, when sampling for photons with lower energy levels, we reject any RITA samples that fall outside the desired interval _for_that_photon_.
Next thing to do is a major refactoring. I want to have the Rayleigh sampling function and the Compton sampling function to be in separate files, for ease of maintenance.
Labelled variables with their units, e.g [eV], [g/cm^3]. The only iffy thing is the units of the density data in variable 'volume' in project_kernel.cu. I'm -pretty- sure that [g/cm^3] is correct, but I am awaiting confirmation on that.
In __init__.py, __all__ is supposed to be a list of strings (https://docs.python.org/3/tutorial/modules.html). Initial (and partial) data dump of Mean Free Path data for scatter simulation. Mainly a proof of concept. Will be updating shortly to have all of the data (which is A LOT). Issues with the MFP data: I do not currently have data corresponding to the following materials (see material_coefficients.py): - soft tissue - iron - teflon - anything called "external"
All of the MFP data is now there. Next up is RITA parameters
Loaded and formatted the RITA paramter information from MC-GPU. Also created a dictionary that maps material names to a RITA sampler for that material.
This should be the last MC-GPU "data dump." Next step is to refactor/rewrite the photon-tracking and scatter-sampling functions.
My Rayleigh and Compton sampling function are set. The last few things to do are: - Figure out the geometry for efficiently generating random [initial_dir]s - calculate the initial position ***in the volume*** from the the X-Ray source position and the initial direction - Handle interpolation to get from photon_energy to MFP values - Handle crossing of segment boundaries - Handle "leaving the volume" detection - Figure out the geometry for determining which [photons that left the volume] hit the detector, and where Once I handle all of that, the next step is to do Variance Reduction-type stuff as laid out in Sisniega, and -then- translate to CUDA
"Area density" refers to the product of [linear distance of the ray through material 'm'] and [density of the material], and it's the quantity used in mass attenuation calculations.
Figured out a clean geometric solution to representing the 'plane vector' of the detector plane, complete with thorough comments. Is very much better than the "analytically find the global minimum of a function of two variables" method that I was trying to use for the "distance from the origin" part of the detector plane vector. Still some TODOs before implementing variance reduction techniques, namely: 1. Math for determining which pixel got hit, in function track_single_photon_no_vr(...) 2. Converting the photon energies into image intensities in function simulate_scatter_no_vr(...)
Wrapper for all the data and functionality needed for calculations relating to: - Does a photon's path hit a plane, and where - Does the photon hit a specific region on the plane
1. Moved PlaneSurface class to a separate file 2. Modified functions to use the PlaneSurface class instead of an amalgamation of np.ndarrays 3. Fixed the Woodcock tracking logic and added comments that refer to useful resources for understanding Woodcock tracking
Was having nagging doubts about the purported units of various quantities. What was formerly "intensity" (purported units: [energy] / [area * time]) in project_kernel(...) is now properly called "deposited_energy" (purported and actual units: [energy], namely electron-volts). Additionally, tightened up the mass attenuation GPU code to eliminate unnecessary division operations (boooo!!!!) as well as reduce the number of [multiplications by -1] (minor-er booo). The scatter estimation code now is such that it can be added to the X-ray primary (from the ray-casting) prior to the neglog transform. TODOs: - Is the (if self.collected_energy) check mislabeled? It seems like it would instead be a conversion from [deposited energy] to [intensity] - Run the current code to see if it compiles, and if it does compile, if it produces a scatter-looking image
ScatterNet is no longer a thing
Confirmed voxel-centered coordinate system from the projector.py:SingleProjector.project(...) arguments list
Plus a slight improvement to the get_scattered_dir(...) function -- now, only normalizes the direction unit vector when the square of the magnitude gets too out of whack (by 1.0e-14)
I think that the psurface_check_ray_intersection(...) function might end up being the only PlaneSurface function that I need to translate from Python to CUDA. I can ignore the other PlaneSurface functions if I end up using the index_from_ijk matrix to obtain the pixel coord.s in the kernel code.
Also fixed some stuff with the interaction-type sampling in Python. Currently having trouble combining the mass-attenuated primary image with the scatter image; the background color of the final image doesn't seem the shade that it should be. I think the reason that combining the image and noise is faulty is because the mass-attenuation code and the scatter code are returning different quantities. In project_kernel(...), deposited_energy[...] has units of eV, BUT SCALED by 1) the attenuation and 2) the pdf for each energy's spectral bin. To achieve the ACTUAL value of the deposited energy, I think the currently-returned value needs to be multiplied by the photon count used in the projection. Once I have taken that step, I think the primary and scatter should mesh if the photon count used for both is equal. That will be something to investigate.
A few more TODOs: - get_mat_mfp_data(...) and get_wc_mfp_data(...). These are basically wrappers for binary search of a sorted array - Rewriting a calculation in sample_Compton(...) to remove divisions in favor of multiplications. This will primarily be a pencil-and-paper task before finally altering the code. - Implementing the RANECU seed initialization
scatter_correction_kernel.cu will implement the scatter correction described in Sisniega et al. (2015), "High-fidelity artifact correction for cone-beam CT imaging of the brain"
The image difference between the "images" and "images with noise" suggests that the noise for the wire-box-around-cube phantom affects the regions surrounding edges in the volume.
I had compared the wrong images :(
I realized an issue with an assumption that I was making. I had assumed that (rt_kinv) @ (u,v,1)^T gives a vector pointing all the way from the source to the pixel [u,v]. However, this readout clearly demonstrates otherwise, since the source-to-detector distance is around 1200 mm >> 1.04mm:
vector in voxel-space along ray from source-point to pixel at [0,0] on the detector plane:
(rt_kinv) @ (0.500000, 0.500000, 1) = (-0.123871, 0.160038, 1.000000)^T
vector in voxel-space along ray from source-point to pixel at [1,0] on the detector plane:
(rt_kinv) @ (1.500000, 0.500000, 1) = (-0.123871, 0.159779, 1.000000)^T
vector in voxel-space along ray from source-point to pixel at [0,1] on the detector plane:
(rt_kinv) @ (0.500000, 1.500000, 1) = (-0.123613, 0.160038, 1.000000)^T
vector in voxel-space along ray from source-point to pixel at [1,1] on the detector plane:
(rt_kinv) @ (1.500000, 1.500000, 1) = (-0.123613, 0.159779, 1.000000)^T
The way to recover my previous logic is to make myself right: wherever I used a quantity [(rt_kinv) @ (u,v,1)^T], I need to scale that vector up by a certain factor.
After some investigating, I determined that the vector [(rt_kinv) @ (W/2,H/2,1)^T], which points along the ray used to measure the source-to-detector distance, always has magnitude 1 (even when printed as "%1.10f"!). Accordingly, my previous code and logic should work once I scale up vectors by a factor to source_to_detector_distance.
After a run with photon_count=100,000, I obtained the following readouts:
NOISE:
pixel hit data: [pixel_x, pixel_y], initial_energy -> energy
[174, 668], 39000.0 -> 32409.08329237009
[1230, 304], 21500.0 -> 21375.622025088112
[262, 441], 57500.0 -> 44863.22686201331
[807, 240], 51000.0 -> 50901.5307396224
[201, 786], 62500.0 -> 62235.16826503576
[492, 600], 69000.0 -> 68345.39473533675
[1063, 546], 55500.0 -> 53375.739028593176
[299, 909], 32500.0 -> 27064.73843459087
[422, 269], 71500.0 -> 70988.29566130895
[106, 922], 56000.0 -> 55241.5741016763
[962, 88], 50500.0 -> 47400.34269330423
[1209, 901], 26500.0 -> 21581.528628207667
[207, 819], 44000.0 -> 43266.048531480315
[856, 861], 33500.0 -> 32480.67659711316
[807, 288], 54500.0 -> 54202.66406735478
[940, 410], 46000.0 -> 45592.6781458572
[352, 849], 59500.0 -> 58324.3978956708
[928, 692], 53500.0 -> 53047.301695281865
[810, 115], 25000.0 -> 24859.62207585512
[898, 777], 35000.0 -> 34818.37665189842
[564, 424], 71500.0 -> 71044.01041902752
[1003, 164], 48000.0 -> 47184.715205093824
[755, 557], 37000.0 -> 36357.37124382156
[1020, 484], 49500.0 -> 48205.98102014361
[478, 164], 58500.0 -> 58297.87841406196
[575, 580], 22500.0 -> 19020.49839474008
[340, 190], 36000.0 -> 25506.753995187417
[1213, 618], 59500.0 -> 59414.525009229765
X-RAY PRIMARY:
X-ray primary data: [pixel_x, pixel_y], deposited_energy
[174, 668], 4924489.974975586
[1230, 304], 4924489.974975586
[262, 441], 4924489.974975586
[807, 240], 4924489.974975586
[201, 786], 4924489.974975586
[492, 600], 1023486.8049621582
[1063, 546], 4924489.974975586
[299, 909], 4924489.974975586
[422, 269], 2564818.3822631836
[106, 922], 4924489.974975586
[962, 88], 4924489.974975586
[1209, 901], 4924489.974975586
[207, 819], 4924489.974975586
[856, 861], 2955480.0033569336
[807, 288], 4924489.974975586
[940, 410], 4924489.974975586
[352, 849], 4924489.974975586
[928, 692], 4924489.974975586
[810, 115], 4924489.974975586
[898, 777], 4924489.974975586
[564, 424], 1024391.7465209961
[1003, 164], 4924489.974975586
[755, 557], 4924489.974975586
[1020, 484], 4784326.553344727
[478, 164], 4924489.974975586
[575, 580], 1024166.8701171875
[340, 190], 4924489.974975586
[1213, 618], 4924489.974975586
My hypothesis for way these values differ by orders of magnitude (which then cause the noise to not show up on the purported "primary+noise" image) is that I scaled the results of project_kernel(...) incorrectly.
I think I was correct that project_kernel(...) computes the deposited energy per unit photon. However, I incorrectly tried to remedy this by multiplying the result by photon_count. Instead, I should multiply the result [deposited energy per unit photon] by the quantity [photon_count] x [solid angle covered by the pixel]. Under that system, the photon_count parameter for the Projector class refers to the number of photons emitted by the X-Ray source irrespective of direction, of which a small fraction eventually hit the volume. Indeed, over [0.5 x photon_count] of the emitted photons travel away from the volume being imaged.
The next step, which I am about to do, is to modify project_kernel(...) so that it takes as additional arguments:
- source_to_detector distance
- photon_count
The source_to_detector_distance is necessary for calculating the solid angle covered by the pixel. There may be other additional information needed to calculate the solid angle that I did not mention above, but I will add those to the project_kernel(...) parameter list when needed.
Really nice results for 10^6 photons. It seems like the noise gets drowned out by the X-ray primary at the 5x10^6 photon level. Thinking about whether that makes sense is a TODO
…inear-combination
index_from_ijk version:
index_from_ijk is a 3x4 transform. I care about the upper-left 2x3 submatrix.
hit_x = (hit.data[0] - source_ijk.data[0]) / source_to_detector_distance
hit_y = (hit.data[1] - source_ijk.data[1]) / source_to_detector_distance
hit_z = (hit.data[2] - source_ijk.data[2]) / source_to_detector_distance
pixel_x = index_from_ijk[0,0] * hit_x + index_from_ijk[0,1] * hit_y + index_from_ijk[0,2] * hit_z
pixel_y = index_from_ijk[1,0] * hit_x + index_from_ijk[1,1] * hit_y + index_from_ijk[1,2] * hit_z
Linear-combination version: already implemented in the PlaneSurface class
I want to scale by [photon count] x [proportion of sphere covered by pixel] == [photon count] x [solid angle subtended by pixel] / [solid angle subtended by sphere] == [photon_count] x [\Omega] / [4 \pi]
Previously, there were a few theoretical cases where wacky input coordinate systems might possibly cause it to mess up. (The phantom is not one of those wacky cases.)
I forgot to do the arctan(...) step. Also removed an erroneous factor of [photon_count] in the Python wrapper
A bit more work to be done to have PyCuda properly call the scatter kernel. Additionally, I still need to verify the math for combining the scatter signal with the primary signal. I think I should return the 'intensity' output to being [deposited energy per photon]. That output would not be distorted by the differences in solid angle between the pixels, and thus is useful for the intensity_upper_bound clipping.
I implemented the collected-energy conversion logic. Something that would be quite nice would be to intelligibly return the scatter signal as some sort of "signal per photon" value. Would it be "signal per simulated photon"? Would it be "signal per detected scattered photon" (I think so)? Thinking about this is a TODO
I have implemented the PyCuda wrapping of the scatter kernel TODO: can I reuse the intensity_gpu allocation for the scatter kernel's output? I'd like to reduce the amount of memory I'm using on the GPU
- corrected the MEMSIZE field for cuda_scatter_structs.py. My predicted struct sizes were different than those returned by sizeof(...). - corrected _convert_to_millimeters(...) in mcgpu_mfp_data.py. - In PlaneSurface, changed the assertion from "equals zero" to "is close to zero" to account for floating-point numerical errors Some substantive bug fixes: - my linear interpolation used subtraction instead of addition for scatter_kernel.cu:get_mfp_data(...) (yikes!) - edited get_voxel_1D so that once a photon is moved to the volume, it actually registers as being inside Also confirmed with printf(...) that the physical constants (MFP, Compton, etc.) are properly copied to the GPU. TODO: - check what sort of choices for thread and block dimensions are fastest - if I use multiple blocks, my seed_input needs to be different for each. Handling that _intelligently_ is definitely a TODO
Also, note that on the given phantom, the results seem quite strange, with slices of scatter. Additionally, scatter seems more intense on the left of the image.
Trying to compare to MCGPU
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
Meant to merge into dev