T1 Mapping in Python?

Question

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:

Answer 1

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[0] * x**2 + p[1] * x + p[2] 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[0] * x + p[1] 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."