In [1]:
## Following MATLAB code from http://capture-clarity.org/clarity-based-tractography/
In [2]:
import os
In [3]:
## Parameters (the script loops through all parameters and saves each result automatically)
dogsigmaArr = [1]; # Sigma values for derivative of gaussian filter, recommended value: 0.6 - 1.3 (based on actual data)
gausigmaArr = [2.3]; # Sigma values for gaussian filter, recommended value: 1.3 - 2.3 (based on actual data)
angleArr = [25]; # Angle thresholds for fiber tracking, recommended value: 20 - 30
In [4]:
pwd
Out[4]:
In [5]:
import numpy as np
In [6]:
import math
from scipy import ndimage
import nibabel as nib
In [7]:
## Change later on
file_path = "bigtiff"
directory = os.path.dirname(file_path)
In [8]:
directory = '/root/bigtiff'
In [9]:
cd /root/bigtiff/
In [10]:
print len([name for name in os.listdir('.') if os.path.isfile(name)])
In [11]:
from PIL import Image
In [12]:
im = Image.open('sample3.tiff') # Needs to be changed to dynamically go down list of fnDataArr (currently just loads same test image)
# Omitted: channel data (red/green - our CLARITY data was single channel, so no channel data loaded.)
img_data = np.asarray(im);
print img_data.shape
In [16]:
print "Start Computing Structure Tensor..."
fnDataArr = len([name for name in os.listdir('.') if os.path.isfile(name)])
for ii in range(fnDataArr):
print "Iteration: " + str(ii);
# Set up results directory
if not os.path.exists(directory):
os.makedirs(directory)
im = Image.open('sample3.tiff') # Needs to be changed to dynamically go down list of fnDataArr (currently just loads same test image)
# Omitted: channel data (red/green - our CLARITY data was single channel, so no channel data loaded.)
img_data = np.ones((100, 100, 100));
for jj in range(len(dogsigmaArr)):
dogsigma = dogsigmaArr[jj];
print "Start DoG Sigma on " + str(dogsigma);
# Generate dog kernels
dogkercc = doggen([dogsigma, dogsigma, dogsigma]);
dogkerrr = np.transpose(dogkercc, (1, 0, 2));
dogkerzz = np.transpose(dogkercc, (0, 2, 1));
# Compute gradients
grr = ndimage.convolve(img_data, dogkercc);
gcc = ndimage.convolve(img_data, dogkerrr);
gzz = ndimage.convolve(img_data, dogkerzz);
# Compute gradient products
gprrrr = grr * grr;
gprrcc = grr * gcc;
gprrzz = grr * gzz;
gpcccc = gcc * gcc;
gpcczz = gcc * gzz;
gpzzzz = gzz * gzz;
# Compute gradient amplitudes
ga = (gprrrr + gpcccc + gpzzzz)**(1/2);
print "ga shape: " + str(ga.shape);
# Convert numpy ndarray object to Nifti data type
gradient_amplitudes_data = nib.Nifti1Image(ga, affine=np.eye(4));
# Save gradient amplitudes image
nib.save(gradient_amplitudes_data, 'gradient_amplitudes.nii');
# Compute gradient vectors
gv = np.concatenate((grr[..., np.newaxis], gcc[..., np.newaxis], gzz[..., np.newaxis]), axis = 3);
# Convert numpy ndarray object to Nifti data type
gradient_vectors_data = nib.Nifti1Image(gv, affine=np.eye(4));
# Save gradient vectors
nib.save(gradient_vectors_data, 'gradient_vectors.nii');
# Compute structure tensor
for kk in range(len(gausigmaArr)):
gausigma = gausigmaArr[kk];
print "Start Gauss Sigma with gausigma = " + str(gausigma);
print "Generating Gaussian kernel..."
gaussker = gaussgen([gausigma, gausigma, gausigma]);
print "Blurring gradient products..."
gprrrrgauss = ndimage.convolve(gprrrr, gaussker);
gprrccgauss = ndimage.convolve(gprrcc, gaussker);
gprrzzgauss = ndimage.convolve(gprrzz, gaussker);
gpccccgauss = ndimage.convolve(gpcccc, gaussker);
gpcczzgauss = ndimage.convolve(gpcczz, gaussker);
gpzzzzgauss = ndimage.convolve(gpzzzz, gaussker);
print "Saving a copy for this Gaussian sigma..."
tensorfsl = np.concatenate((gprrrrgauss[..., np.newaxis], gprrccgauss[..., np.newaxis], gprrzzgauss[..., np.newaxis], gpccccgauss[..., np.newaxis], gpcczzgauss[..., np.newaxis], gpzzzzgauss[..., np.newaxis]), axis = 3);
# Convert numpy ndarray object to Nifti data type
tensor_fsl_data = nib.Nifti1Image(tensorfsl, affine=np.eye(4));
nib.save(tensor_fsl_data, "slice_" + str(ii) + "dogsigma_" + str(jj) + "gausigma_" + str(kk) + 'tensorfsl.nii');
print 'Complete!'
In [14]:
'''
Function to generate derivatives of Gaussian kernels, in either 1D, 2D, or 3D.
Source code in MATLAB obtained from Qiyuan Tian, Stanford University, September 2015
Edited to work in Python by Tony
'''
def doggen(sigma):
halfsize = np.ceil(3 * max(sigma))
x = range(np.single(-halfsize), np.single(halfsize + 1)); # Python colon is not inclusive at end, while MATLAB is.
dim = len(sigma);
print sigma
if dim == 1:
X = np.array(x); # Remember that, by default, numpy arrays are elementwise multiplicative
return -X * np.exp(-X**2/(2 * sigma**2));
elif dim == 2:
[X, Y] = np.meshgrid(x, x);
return -X * np.exp(-X**2/(2*sigma[0]^2) * np.exp(-Y**2))
elif dim == 3:
[X, Y, Z] = np.meshgrid(x, x, x);
return X * np.exp(np.divide(-X**2, 2 * sigma[0]**2)) * np.exp(np.divide(-Y**2, 2 * sigma[1]**2)) * np.exp(np.divide(-Z**2, 2 * sigma[2]**2))
else:
print 'Only supports up to 3 dimensions'
In [15]:
'''
Function to generate Gaussian kernels, in 1D, 2D and 3D.
Source code in MATLAB obtained from Qiyuan Tian, Stanford University, September 2015
Edited to work in Python by Tony.
'''
def gaussgen(sigma):
halfsize = np.ceil(3 * max(sigma));
x = range(np.single(-halfsize), np.single(halfsize + 1));
dim = len(sigma);
if dim == 1:
k = np.exp(-x**2 / (2 * sigma^2));
elif dim == 2:
[X, Y] = np.meshgrid(x, x);
k = np.exp(-X**2 / (2 * sigma[0]**2)) * np.exp(-Y**2 / (2 * sigma[1]**2));
elif dim == 3:
[X, Y, Z] = np.meshgrid(x, x, x);
k = np.exp(-X**2 / (2 * sigma[0]**2)) * np.exp(-Y**2 / (2 * sigma[1]**2)) * np.exp(-Z**2 / (2 * sigma[2]**2));
else:
print 'Only supports up to dimension 3'
return k / sum(abs(k));
In [ ]:
In [ ]:
In [ ]: