# 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]