Raytracing2 #14

jjgomezcadenas · 2021-09-28T10:11:55Z

This PR completes the first version of raytracking functions, which for the time being include only the definition of a Ray, the definition of a Cylinder and the propagation of a ray to a cylinder. This functions are useful to correct the LORS in the NEMA5 test by the absorption in the phantom, and thus compute correctly the sensitivity.

The package Polynomials was added previously to compute roots of the ray propagator, but it turns out that hand calculation is more convenient here. Dependency not necesary now, thus removed.

Testing notebooks are added to the notebook respository. The testing notebook exercise some of the functions in the package but are not intended or guaranteed to be complete. All tests in the notebooks (and some times more) are included in the automatic testing.

andLaing

Some questions about naming and functionality. There seem to be raytracing packages already written but it's true that they're maybe a bit more complex than what we need.

andLaing · 2021-09-28T14:43:33Z

src/histos.jl

 """
 function hist1d(x::Vector{T}, bins::Vector{S}, norm=false) where {T<:Number,S<:Number}
-    h = fit(Histogram, x, bins)
+    h = StatsBase.fit(Histogram, x, bins)


This change doesn't seem necessary, especially since it's not used consistently even in histos.jl

This may be a fossile. At some stage there was a naming conflict.

andLaing · 2021-09-28T14:47:15Z

src/raytracing.jl

-    p::Vector{Float64}
+    e::Vector{Float64}
    d::Vector{Float64}
    u::Vector{Float64}


Why have you chosen to change to e rather than p? And without changing the docstring.

Also, A ray for me should be point and unit direction so it's not so clear what the additional "Direction vector" is for.

I decided to follow the nomenclature here:
https://www.cl.cam.ac.uk/teaching/1999/AGraphHCI/SMAG/node2.html

If one takes d as specifying a direction (not a unit direction) then one can compute the unit vector internally.

andLaing · 2021-09-28T14:48:12Z

src/raytracing.jl

-    function Ray(p , d)
-        new(p, d, unit_vector(d,p))
+    function Ray(e , d)
+        new(e, d, unit_vector(d,e))


The use of unit_vector here contradicts its meaning below.

Not sure I follow you here

The unit vector for a ray should be along the direction vector (so basically d / mag(d)). Your unit_vector function calculates a unit vector along the difference of the two inputs, this means that here you're effectively using the origin as if it were a direction and forming a triangle with the direction vector.

Right, seems there is also a mistake there.

andLaing · 2021-09-28T14:54:41Z

src/raytracing.jl

 function propagate_ray(r::Ray)
    function xray(t::Real)
-        return r.p + t * r.d
+        return r.e + t * r.d


For the propagation to be correct the function should use a unit vector so by the Ray definition above it should be r.e + t * r.u

Yes, I think you are right.

andLaing · 2021-09-28T15:03:56Z

src/raytracing.jl

 Return a unit vector in the direction defined by points p0 and p1  
 """
-unit_vector(p0::Vector{Real}, p1::Vector{Real}) = (p1 - p0) ./ mag(p1 - p0)
+unit_vector(p0::Vector{Float64}, p1::Vector{Float64}) = (p1 - p0) ./ norm(p1 - p0)


Broadcasting unnecessary.

Seems unusual that this type of function isn't in LinearAlgebra or similar

It may be, I didn't look.

andLaing · 2021-09-28T15:10:31Z

src/raytracing.jl



-length(c::Cylinder) = abs.(c.zmax - c.zmin)
+cylinder_length(c::Cylinder) = abs.(c.zmax - c.zmin)


c.zmax and c.zmin are Float64, there's no broadcasting to be done.

YES! good catch

andLaing · 2021-09-28T15:40:49Z

src/raytracing.jl

+    xray = ATools.propagate_ray(r)
+
+    # If no real roots return nothing 
+    arg = b^2 - 4*a*c


I think you should have a protection here after calculating the determinant. Otherwise a ray that doesn't intersect the barrel anywhere will error on the sqrt

I will check, thanks

I will check

andLaing · 2021-09-28T15:51:55Z

src/raytracing.jl

+
+
+"""Computes intersection between a ray and a cylinder"""
+function propagate(r::Ray, cy::Cylinder)


The name is misleading again. The function seems to look for intersections between the ray and cylinder. i think that might be a better thing to base the name on. intersection or find_intersections?

Yes, right again

andLaing · 2021-09-28T16:09:06Z

src/raytracing.jl

+
+    # belong-to-cylinder condition 
+    zc1 = cy.zmin < z1 < cy.zmax
+    zc2 = cy.zmin < z2 < cy.zmax


It's not clear to me exactly what the function is meant to do. You keep both t1 and t2 which made me think that you wanted intersections for back-to-back gammas but then the function only returns one position.

If you assume that the origin is inside the cylinder (should be true in our application), positive t is the only relevant intersection if we only want to follow the ray. If both are possible but only one returned the meaning of the function is ambiguous. If you only want the forward intersections then a lot of the if else below can be removed in favour of taking the positive t. If you want a function that gives you the two/one/or no intersection(s) of the back-to-back gammas then I'd return both with nothing or missing in the case of no valid intersection.

https://www.cl.cam.ac.uk/teaching/1999/AGraphHCI/SMAG/node2.html

I understand the calculations for where the intersections should be. My question is what you want to return. If you only want the forward going ray then you can check if one of the t is positive, if you only want to consider origins within the cylinder.

In the file you link they already talk about t > 0 being the only really relevant root for a ray my question about what you want is inspired by the back-to-back nature of the specific problem the code is written for.

Let me see if I follow you. The argument is that the solution with t<0 also works, given that there are two back-to-back gammas?

If you assume back-to-back gammas both t > 0 and t < 0 can be valid (if you are inside the cylinder you can only have one positive and 1 negative or neither defined). In the case that you want both I'd check each intersect for validity (intersect within the limits of interest) and return them both with something like a missing or nothing for any of them that aren't valid. If you're only interested in the forward ray (t > 0) you can already ignore one of them from earlier in the function and save yourself some of the if else

andLaing · 2021-09-28T16:09:57Z

src/raytracing.jl

+
+
+"""Returns True if point in end-caps of cyinder"""
+in_endcaps(c:: Cylinder, p::Vector{Float64}) = isapprox(p[3], c.zmin) || isapprox(p[3], c.zmax)


This alone doesn't test if the point is in an endcap. Only that the one-dimensional position is level with the endcaps

That is what was intended, name probably not very well chosen

jjgomezcadenas added 4 commits September 28, 2021 09:44

Adding Polynomials for raytracing.

0ccf5e0

Functions for propagation of a ray to a cylinder with tests.

00522ba

Remove dependence from Polynomials.

f48029e

The package Polynomials was added previously to compute roots of the ray propagator, but it turns out that hand calculation is more convenient here. Dependency not necesary now, thus removed.

Added two notebooks.

715917e

Testing notebooks are added to the notebook respository. The testing notebook exercise some of the functions in the package but are not intended or guaranteed to be complete. All tests in the notebooks (and some times more) are included in the automatic testing.

andLaing reviewed Sep 28, 2021

View reviewed changes



		length(c::Cylinder) = abs.(c.zmax - c.zmin)
		cylinder_length(c::Cylinder) = abs.(c.zmax - c.zmin)



		"""Computes intersection between a ray and a cylinder"""
		function propagate(r::Ray, cy::Cylinder)



		"""Returns True if point in end-caps of cyinder"""
		in_endcaps(c:: Cylinder, p::Vector{Float64}) = isapprox(p[3], c.zmin) \|\| isapprox(p[3], c.zmax)

Raytracing2 #14

Are you sure you want to change the base?

Raytracing2 #14

Uh oh!

Conversation

jjgomezcadenas commented Sep 28, 2021

Uh oh!

andLaing left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants