-
Notifications
You must be signed in to change notification settings - Fork 2
/
calc_t1map_cpp.py
85 lines (63 loc) · 2.19 KB
/
calc_t1map_cpp.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
# -*- coding: utf-8 -*-
"""
voxelwise T1 mapping in C++
contributors: Yoon-Chul Kim, Khu Rai Kim, Hyelee Lee
"""
import numpy as np
import matplotlib.pyplot as plt
import os
import sys
import warnings
import copy
import pickle
import skimage.morphology as sm
import time
import pydicom as dicom
sys.path.append(os.path.join(os.path.dirname(__file__), '.', 'x64/Release'))
import t1_mapping as t1 # t1_mapping.pyd
from skimage import data, io, filters
from numpy import arange, sin, pi
from scipy.optimize import curve_fit
from math import floor, sqrt
def func_orig(x, a, b, c):
return a*(1-np.exp(-b*x)) + c
def calculate_T1map_cpp_rd(ir_img, inversiontime, multicore_flag, zoom=2, zoom_lenz=16):
'''
implementation of Barral's method
'''
if inversiontime[-1] == 0:
inversiontime = inversiontime[0:-1]
if ir_img.shape[2] > inversiontime.shape[0]:
ir_img = ir_img[:,:,0:ir_img.shape[2]-1]
zoom_lenz = 16
shape = ir_img.shape
flat_t1map = t1.fit_t1_barral(ir_img.flatten().astype(np.float64),
np.array(inversiontime, dtype=np.float64),
shape[0], shape[1], shape[2],
zoom, zoom_lenz, multicore_flag)
t1_map = np.reshape(flat_t1map, [shape[0], shape[1], 3])
return t1_map
def calculate_T1map_cpp_lm(ir_img, inversiontime, multicore_flag):
nx, ny, nti = ir_img.shape
y = np.zeros(nti)
prtno = 1 # post
if prtno == 0: # pre
p0_initial = [350, 0.001, -150]
else: # post
p0_initial = [350, 0.005, -150]
if inversiontime[-1] == 0:
nTI = 7
inversiontime = inversiontime[0:-1]
y = y[0:-1]
else:
nTI = 8
err_tol = 1e-7
ir_img = ir_img[:,:, :nTI]
shape = ir_img.shape
flat_t1map = t1.fit_t1(ir_img.flatten(),
np.array(inversiontime),
np.array(p0_initial),
shape[0], shape[1], shape[2],
err_tol, multicore_flag)
t1_map = np.reshape(flat_t1map, [nx, ny, 3])
return t1_map