Tissue classification using MAP-MRF

dipy/segment/energy_mrf.py

@@ -0,0 +1,57 @@
+from __future__ import division, print_function, absolute_import
+import numpy as np
+
+
+def total_energy(masked_image, masked_segmentation,
+                 mu, var, index, label, beta):
+
+    energytotal = log_likelihood(masked_image, mu, var, index, label)
+    energytotal += gibbs_energy(masked_segmentation, index, label, beta)
+
+    return energytotal
+
+
+def log_likelihood(img, mu, var, index, label):
+
+    loglike = ((img[index] - mu[label]) ** 2) / (2 * var[label])
+
+    loglike += np.log(np.sqrt(var[label]))
+
+    return loglike
+
+
+def gibbs_energy(seg, index, label, beta):
+
+    energy = 0
+
+    if label == seg[index[0] - 1, index[1], index[2]]:
+        energy = energy - beta
+    else:
+        energy = energy + beta
+
+    if label == seg[index[0] + 1, index[1], index[2]]:
+        energy = energy - beta
+    else:
+        energy = energy + beta
+
+    if label == seg[index[0], index[1] - 1, index[2]]:
+        energy = energy - beta
+    else:
+        energy = energy + beta
+
+    if label == seg[index[0], index[1] + 1, index[2]]:
+        energy = energy - beta
+    else:
+        energy = energy + beta
+
+    if label == seg[index[0], index[1], index[2] - 1]:
+        energy = energy - beta
+    else:
+        energy = energy + beta
+
+    if label == seg[index[0], index[1], index[2] + 1]:
+        energy = energy - beta
+    else:
+        energy = energy + beta
+
+    return energy

dipy/segment/icm_map.py

@@ -0,0 +1,59 @@
+
+
+# Use ICM to segment T1 image with HMRF
+
+import numpy as np
+import nibabel as nib
+
+img = nib.load('/Users/jvillalo/Documents/GSoC_2015/Code/Data/3587_BL_T1_to_MNI_Linear_6p.nii.gz')
+dataimg = img.get_data()
+
+mask_image = nib.load('/Users/jvillalo/Documents/GSoC_2015/Code/Data/3587_mask.nii.gz')
+datamask = mask_image.get_data()
+
+from dipy.segment.mask import applymask
+masked_img = applymask(dataimg,datamask)
+print('masked_img.shape (%d, %d, %d)' % masked_img.shape)
+shape=masked_img.shape[:3]
+
+# Still must do Otsu for 3 classes. 
+#from init_param import otsu_param
+#mu_tissue1, sig_tissue1, mu_tissue2, sig_tissue2 = otsu_param(masked_img, 4, 4)
+
+seg = nib.load('/Users/jvillalo/Documents/GSoC_2015/Code/Data/FAST/3587_BL_T1_to_MNI_Linear_6p_seg.nii.gz')
+seg_initial = seg.get_data()
+
+nclass = 3
+nh = 6   #neighborhood
+niter = 1
+totalE = np.zeros((shape[0],shape[1],shape[2],nclass))
+
+from dipy.segment.rois_stats import seg_stats
+seg_stats(masked_img, seg_initial, 3)
+
+from dipy.core.ndindex import ndindex
+
+for idx in np.ndindex(shape):
+    if not masked_img[idx]:
+        continue
+
+    while True:
+
+        mu, std = seg_stats(masked_img, seg_initial, 3)
+
+        for l in range(0,nclass-3):
+            totalE = 1/(2*vars[l])*(masked_img-mus[l])^2 + np.log(sigs[l]) - beta*Himg[:,:,:,l]
+        np.amin(totalE[-1])
+
+    N = niter+1        
+    if N == 1:
+        break
+
+
+
+
+
+
+
+
+

dipy/segment/init_param.py

@@ -0,0 +1,19 @@
+
+import numpy as np
+from dipy.segment.mask import multi_median
+from dipy.segment.threshold import otsu
+
+def otsu_param(input_volume, median_radius=4, numpass=4):
+
+    input_filtered = multi_median(input_volume, median_radius, numpass)
+    thresh = otsu(input_filtered)
+    tissue1 = input_filtered > thresh
+    tissue2 = input_filtered < thresh
+
+    mu_tissue1 = np.mean(tissue1)
+    sig_tissue1 = np.std(tissue1)
+    mu_tissue2 = np.mean(tissue2)
+    sig_tissue2 = np.std(tissue2)
+
+    return mu_tissue1, sig_tissue1, mu_tissue2, sig_tissue2
+

dipy/segment/neighborh.py

@@ -0,0 +1,24 @@
+
+
+# Compute the neighborhood. For now only 6. La
+
+import numpy as np
+
+def neighbor3D(input_volume,nhood):
+
+    volume_shape = input_volume.shape[:3]
+    Nhood = np.ones((volume_shape[0],volume_shape[1],volume_shape[2],nhood))
+
+
+
+    a = np.arange(shape[0])
+    b = np.arange(shape[1])
+    c = np.arange(shape[2])
+
+    A = a[1:-1]
+    B = b[1:-1]
+    C = c[1:-1]    
+
+    Nhood1 = input_volume[]
+
+

dipy/segment/rois_stats.py

@@ -0,0 +1,37 @@
+from __future__ import division, print_function, absolute_import
+import numpy as np
+
+
+def seg_stats(input_image, seg_image, nclass):
+    r""" Mean and standard variation for 3 tissue classes
+
+    1 is CSF
+    2 is gray matter
+    3 is white matter
+
+    Parameters
+    ----------
+    input_image : ndarray
+        blah blah
+
+
+
+    Returns
+    -------
+    mu, std : float
+        Mean and standard deviation for every class
+
+
+    """
+    mu = np.zeros(3)
+    std = np.zeros(3)
+
+    for i in range(1, nclass + 1):
+
+        H = input_image[seg_image == i]
+
+        mu[i - 1] = np.mean(H, -1)
+        std[i - 1] = np.std(H, -1)
+
+
+    return mu, std

dipy/segment/tests/test_icm_map.py

@@ -0,0 +1,70 @@
+# Use ICM to segment T1 image with MRF
+
+import numpy as np
+import nibabel as nib
+
+from dipy.segment.mask import applymask
+from dipy.core.ndindex import ndindex
+from dipy.segment.rois_stats import seg_stats
+from dipy.segment.energy_mrf import total_energy
+
+# dname = '/Users/jvillalo/Documents/GSoC_2015/Code/Data/'
+dname = '/home/eleftherios/Dropbox/DIPY_GSoC_2015/'
+
+img = nib.load(dname + '3587_BL_T1_to_MNI_Linear_6p.nii.gz')
+dataimg = img.get_data()
+
+mask_image = nib.load(dname + '3587_mask.nii.gz')
+datamask = mask_image.get_data()
+
+masked_img = applymask(dataimg, datamask)
+print('masked_img.shape (%d, %d, %d)' % masked_img.shape)
+shape = masked_img.shape[:3]
+
+seg = nib.load(dname + '3587_BL_T1_to_MNI_Linear_6p_seg.nii.gz')
+seg_init = seg.get_data()
+seg_init_masked = applymask(seg_init, datamask)
+
+print("computing the statistics of the ROIs (CSF, GM, WM)")
+mu, std = seg_stats(masked_img, seg_init_masked, 3)
+
+print(mu)
+print(std)
+
+nclass = 3
+L = range(1, nclass + 1)
+niter = 1
+beta = 0.05
+totalE = np.zeros(nclass)
+N = 0
+
+segmented = np.zeros(dataimg.shape)
+
+while True:
+
+    mu, std = seg_stats(masked_img, seg_init, 3)
+    var = std ** 2
+
+    for idx in ndindex(shape):
+        # print(idx)
+        if not masked_img[idx]:
+            continue
+        for l in range(0, nclass):
+
+            totalE[l] = total_energy(masked_img, seg_init_masked,
+                                     mu, var, idx, l, beta)
+
+        segmented[idx] = L[np.argmin(totalE)]
+
+    N = N + 1
+    if N == niter:
+        break
+
+
+# print('Show results')
+figure()
+imshow(seg_init_masked[:, :, 95])
+figure()
+imshow(segmented[:, :, 95])
+
+

dipy/segment/threshold.py

@@ -37,3 +37,46 @@ def otsu(image, nbins=256):
     idx = np.argmax(variance12)
     threshold = bin_centers[:-1][idx]
     return threshold
+
+    def otsu_3(image, nbins=256):
+
+            """
+    Return two threshold values based on Otsu's method.
+    It is intended to work on a masked T1 brain image.
+    It is intended to give threshold values for CSF/WM and WM/GM. 
+
+    Copied from scikit-image to remove dependency.
+
+    Parameters
+    ----------
+    image : array
+        Input image.
+    nbins : int
+        Number of bins used to calculate histogram. This value is ignored for
+        integer arrays.
+
+    Returns
+    -------
+    threshold : float
+        Threshold values.
+    """
+
+    hist, bin_centers = np.histogram(image, nbins)
+    hist = hist.astype(np.float)
+
+    # class probabilities for all possible thresholds
+    weight1 = np.cumsum(hist)
+    weight2 = np.cumsum(hist[::-1])[::-1]
+
+    # class means for all possible thresholds
+    mean1 = np.cumsum(hist * bin_centers[1:]) / weight1
+    mean2 = (np.cumsum((hist * bin_centers[1:])[::-1]) / weight2[::-1])[::-1]
+
+    # Clip ends to align class 1 and class 2 variables:
+    # The last value of `weight1`/`mean1` should pair with zero values in
+    # `weight2`/`mean2`, which do not exist.
+    variance12 = weight1[:-1] * weight2[1:] * (mean1[:-1] - mean2[1:])**2
+
+    idx = np.argmax(variance12)
+    threshold = bin_centers[:-1][idx]
+    return threshold