# T1 Mapping in Python?

I'm trying to write a simple program to perform T1 mapping of a few MRI Images. My understanding is that all I need is to plot the pixel intensity values by inversion time then fit a curve. The output of fitting the curve is what gives me the T1 values.

This was my attempt:

# Loop through every pixel in image
for row in range(rows):
for col in range(cols):
x = []
y = []
# Add all the inversion time and pixel intesities to list
for image in dicom_images:
x.append(image.InversionTime)
y.append(image.pixel_array[row][col])
# Fit curve
t1 = np.polyfit(x, y, 2)


Assuming what I did above is right so far I am now stuck because I don't know how to take the three coefficients produced by np.polyfit and convert them into a single pixel value. I also don't know if there are additional steps I need to take.

I found some MATLAB examples online but they were really big and hard to follow. This is my first time doing T1 mapping so I'm trying to keep it simple.

Here are some plots of the data points:

• So it has multiple images over time, and you're taking an individual pixel value from multiple images and plotting the curve of it (which is like a logarithmic decay?) and fitting a polynomial to it? If it is logarithmic decay, then why fit a polynomial? Can you pick one of the pixels and plt.plot(x, y) it (and show the image here) to at least show that it's decaying as expected? – endolith Mar 31 at 15:00
• I don't think it has logarithmic decay. It should appear as an exponential curve that is increasing and fitting the curve is how you estimate the T1 value. I've added a link to a paper which explains the math I am trying to replicate. – nikebol906 Mar 31 at 17:45
• Can you plot one (or more) of the pixels over time? We need to clarify that first – endolith Mar 31 at 19:35
• I've plotted three different variation. In all of them the x-axis is the same because the image is taken at the same time but there are different pixel intensity values which are plotted on the y-axis. – nikebol906 Mar 31 at 21:08
• Well, those don't look like anything in particular. :) – endolith Mar 31 at 21:56

## 1 Answer

1. You are fitting a polynomial to a logarithmic decay curve, which is probably wrong. You probably want to fit a log curve to a log curve. You can use the naïve method of np.polyfit(log(x), y, 1) which will probably work fine, or you can weight the points to make it more accurate, or you can use scipy.optimize.curve_fit. All are described here: https://stackoverflow.com/q/3433486/125507
2. The coefficients returned by the polyfit are the coefficients of the polynomial. So with deg=2, you're fitting a parabola to the decay curve, which again is probably wrong: p(x) = p * x**2 + p * x + p With the np.polyfit(log(x), y, 1) method you'll be fitting a straight line to a straight line, so there will only be a slope and an offset: p(x) = p * x + p The slope would then be proportional to the decay rate of the log curve (somehow, which you can figure out for the type of value you need).
3. After getting this working, you can probably vectorize it for speed, fitting all pixels simultaneously instead of using for loops. np.polyfit can handle this, for instance: "Several data sets of sample points sharing the same x-coordinates can be fitted at once by passing in a 2D-array that contains one dataset per column."
• I didn't fully understand what np.polyfit did and as a result I thought 2 was the correct input. I see now that using scipy.optimize.curve_fit is the better option. Although I'm still confused as to how I'd use the outputted coefficients to produce an image. I know you talked about it on but I'm still a bit confused. Could you maybe show me some example code? – nikebol906 Mar 31 at 17:51
• @nikebol906 So you accepted this answer, but did you figure out the issue? – endolith Apr 7 at 15:22
• I just used curve_fit from scipy instead of np.polyfit. I also just changed the code a lot because I realized my attempt to was too minimalistic. – nikebol906 Apr 7 at 16:43