MRI

The following problem was posted on the Julia Discourse:

So I have three magnetic resonance (MR) images: $X, Y, Z$. They come from three different pulse sequences which interact with NMR relaxation parameters ) in known ways. That is, we have equations which given the underlying NMR parameters (specifically, T1, T2 and PD) and the scanner parameters $p^X, p^Y, p^Z$ corresponding to each image $X,Y,Z$, we can model the output intensity value observed in the acquired image at some voxel [$(i, j, k)$, e.g., $X_{i,j,k}$, since we are in 3D].

However, we don’t know the $T_1$, $T_2$ and $PD$ parameters but we do have $X,Y,Z$ and we can separately solve for each $p^{(\cdot)}$. If you haven’t followed any of the previous mediocre explanation, then the problem simply reduces to the following system of equations, for which I am trying to solve $T_1$, $T_2$, and $P_D$ where everything else is given.

$$ \begin{align} \log(X_{i,j,k}) &= \log(PD) + p_0^X + p_1^X T_1 + p_2^X T_1^2 \\ \log(Y_{i,j,k}) &= \log(PD) + p_0^Y + p_1^Y T_1 + p_2^Y / T_2 \\ \log(Z_{i,j,k}) &= \log(PD) + p_0^Z + p_1^Z T_1 + p_2^Z / T_2 \\ \end{align} $$

In previous iterations of this project, the person who built the code used non-linear least squares (in MATLAB, FWIW) to solve for the values of $T_1$, $T_2$ and $PD$. However, the person who built the code left our lab, and with that has gone the knowledge of why that was chosen.

Previously I was using python’s scipy implementation of least squares  to solve the problem at every voxel; however, even with parallelization this takes hours.

Since these are polynomial equations we can use homotopy continuation to solve the problem. Let's setup the system first:

In [13]:
using HomotopyContinuation, DynamicPolynomials
@polyvar x y z PD p[1:3,1:3] T[1:2]

f1 = PD + p[1,1] + p[1,2]*T[1] + p[1,3] * T[1]^2 - x
# note that we eliminate the denominator in equation 2 and 3
f2 = (PD + p[2,1] + p[2,2]*T[1])*T[2] + p[2,3] - y * T[2]
f3 = (PD + p[3,1] + p[3,2]*T[1])*T[2] + p[3,3] - z * T[2]

f = [f1, f2, f3]
Out[13]:
3-element Array{Polynomial{true,Int64},1}:
 p₁₋₃T₁² + p₁₋₂T₁ - x + PD + p₁₋₁     
 p₂₋₂T₁T₂ - yT₂ + PDT₂ + p₂₋₁T₂ + p₂₋₃
 p₃₋₂T₁T₂ - zT₂ + PDT₂ + p₃₋₁T₂ + p₃₋₃

Since we need to solve the same instance many many times it makes sense to first compute the generic number of solutions of this polynomial system and then to use a parameter homotopy.

We need some generic points to generate a generic instance

In [50]:
params = [vec(p); x; y; z]
generic_params = randn(ComplexF64, length(params)) # Note we should sample them randomly
f_generic = subs.(f, Ref(params => generic_params))
Out[50]:
3-element Array{Polynomial{true,Complex{Float64}},1}:
 (0.06560018192766423 + 1.226509609243317im)T₁² + PD + (-0.40076878384840564 + 0.3936976384006835im)T₁ + (-0.4034430867072931 + 2.130445251499418im)   
 PDT₂ + (0.17548682027400092 + 1.0868781469710116im)T₁T₂ + (1.0068779513453254 - 0.9349141017423815im)T₂ + (0.6579055458504889 + 0.11880590069556508im)
 PDT₂ + (-0.6911103187839271 - 1.1844653755707422im)T₁T₂ + (0.216353296242008 - 0.3087205812565972im)T₂ + (-1.3312172733617187 + 1.9140562766862028im) 
In [53]:
generic_result = solve(f_generic)
Out[53]:
AffineResult with 8 tracked paths
==================================
• 2 non-singular finite solutions (0 real)
• 0 singular finite solutions (0 real)
• 6 solutions at infinity
• 0 failed paths
• random seed: 808747

We see that the system has only 2 solutions (instead the maximal possible 8 solutions).

In [54]:
generic_solutions = solutions(generic_result)
Out[54]:
2-element Array{Array{Complex{Float64},1},1}:
 [0.100221+0.505434im, -0.256837+1.29393im, 0.809632-0.787296im]
 [-2.12345+1.73512im, 0.438372-2.01187im, -0.398561+0.217356im] 

Now let's prepare to solve our actual problems:

In [61]:
# Construct a path tracker to track the two solutions to our specific instances
# p₁ and 
tracker = pathtracker(f, generic_solutions;
            parameters=params,
            # we just use som dummy parameters for now
            targetparameters=generic_params,
            startparameters=generic_params,
            affine=true)
Out[61]:
PathTracker()

Out of lack of real world data let's sample some random parameters again

In [56]:
specific_params = randn(length(params));

Note that the parameters order is

In [57]:
params
Out[57]:
12-element Array{PolyVar{true},1}:
 p₁₋₁
 p₂₋₁
 p₃₋₁
 p₁₋₂
 p₂₋₂
 p₃₋₂
 p₁₋₃
 p₂₋₃
 p₃₋₃
 x   
 y   
 z   
In [78]:
# Setup the new parameters
set_parameters!(tracker.homotopy, generic_params, specific_params);

# track to the two solutions
r1 = track(tracker, generic_solutions[1])
r2 = track(tracker, generic_solutions[2])
Out[78]:
PathTrackerResult{Array{Complex{Float64},1}}:
 • returncode → success
 • x → Complex{Float64}[0.454103+9.56595e-24im, 1.3872-3.60672e-23im, 0.132065+1.31272e-24im]
 • t → 0.0 + 0.0im
 • accuracy → 2.7287366859420435e-15
 • accepted_steps → 17
 • rejected_steps → 2
In [88]:
# Lets check that the tracking was successfull for both
r1.returncode == PathTrackerStatus.success && r2.returncode == PathTrackerStatus.success
Out[88]:
true

And we get our two specific solutions:

In [86]:
r1.x
Out[86]:
3-element Array{Complex{Float64},1}:
  -0.4432752095469262 - 1.2672493700951796e-19im
 -0.18955251455468036 + 9.362841922611108e-19im 
  0.20328610096481764 - 7.293918355685165e-20im 
In [87]:
r2.x
Out[87]:
3-element Array{Complex{Float64},1}:
 0.45410326818290414 + 9.565951274599459e-24im 
  1.3871995772087389 - 3.60671834508408e-23im  
  0.1320653938006548 + 1.3127235321293564e-24im

We also see that these are complex vectors, but we are interested in the real solutions.

In [93]:
if maximum(abs  imag, r1.x) < 1e-8
    println("First solution: ", real.(r1.x))
end
if maximum(abs  imag, r2.x) < 1e-8
    println("Second solution: ", real.(r2.x))
end
First solution: [-0.443275, -0.189553, 0.203286]
Second solution: [0.454103, 1.3872, 0.132065]