In [1]:
## Following MATLAB code from http://capture-clarity.org/clarity-based-tractography/
In [1]:
import os
In [2]:
## 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 [3]:
pwd
Out[3]:
In [4]:
import numpy as np
In [5]:
import math
from scipy import ndimage
import nibabel as nib
In [6]:
## Change later on
file_path = "/Users/Kepler/Downloads/CAPTURE_v0.1/tractography copy"
directory = os.path.dirname(file_path)
In [7]:
cd TIFF_stack/
In [8]:
from PIL import Image
from numpy import matlib
from scipy import signal
In [9]:
'''
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 * np.max(sigma))
x = range(np.single(-halfsize), np.single(halfsize + 1)); # Python colon is not inclusive at end, while MATLAB is.
dim = len(sigma);
if dim == 1:
X = np.array(x); # Remember that, by default, numpy arrays are elementwise multiplicative
k = -X * np.exp(-X**2/(2 * sigma**2));
elif dim == 2:
[X, Y] = np.meshgrid(x, x);
k = -X * np.exp(-X**2/(2*sigma[0]^2) * np.exp(-Y**2))
elif dim == 3:
[X, Y, Z] = np.meshgrid(x, x, x);
X = X.transpose(0, 2, 1); # Obtained through vigorous testing (see below...)
Y = Y.transpose(2, 0, 1);
Z = Z.transpose(2, 1, 0);
X = X.astype(float);
Y = Y.astype(float);
Z = Z.astype(float);
k = -X * np.exp(np.divide(-np.power(X, 2), 2 * np.power(sigma[0], 2))) * np.exp(np.divide(-np.power(Y,2), 2 * np.power(sigma[1],2))) * np.exp(np.divide(-np.power(Z,2), 2 * np.power(sigma[2],2)))
else:
print 'Only supports up to 3 dimensions'
return np.divide(k, np.sum(np.abs(k[:])));
In [10]:
'''
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);
X = X.transpose(0, 2, 1); # Obtained through vigorous testing (see below...)
Y = Y.transpose(2, 0, 1);
Z = Z.transpose(2, 1, 0);
X = X.astype(float); # WHY PYTHON?
Y = Y.astype(float);
Z = Z.astype(float);
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 np.divide(k, np.sum(np.abs(k)));
In [40]:
# Set up results directory
if not os.path.exists(directory):
os.makedirs(directory)
# im = Image.open('page1.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.)
ones = np.single(100 * np.ones((10, 10, 10))); #data is hard coded to be np.ones
np.random.seed(1)
norm = np.single(20 + 2 * np.random.randn(10, 10, 10)) # image histogram is normally distributed with mean 20, std 2
uniform = np.single(100 * np.random.rand(10, 10, 10)) # image histogram is normally distributed between 0 and 100
# image has a line of intensity value 100 at col = 4, depth = 4
line = np.zeros((10, 10, 10))
line[:,4,4] = 100
line = np.single(line)
img_data = norm
print img_data[0][0][0];
for jj in range(len(dogsigmaArr)):
dogsigma = dogsigmaArr[jj];
print "Start DoG Sigma on " + str(dogsigma);
# Generate dog kernels
dogkercc = doggen([dogsigma, dogsigma, dogsigma]);
dogkercc = np.transpose(dogkercc, (0, 2, 1)); # annoying
#print dogkercc.shape;
#print dogkercc[:, :, 0];
dogkerrr = np.transpose(dogkercc, (1, 0, 2));
#print dogkerrr[:, :, 0];
dogkerzz = np.transpose(dogkercc, (0, 2, 1));
#print dogkerzz[:, :, 0];
# Compute gradients
grr = signal.fftconvolve(img_data, dogkercc, 'same');
grr = np.transpose(grr, (1, 0, 2));
#print grr[:, :, 0];
gcc = signal.fftconvolve(img_data, dogkerrr, 'same');
gcc = np.transpose(gcc, (1, 0, 2));
#print gcc[:, :, 0];
gzz = signal.fftconvolve(img_data, dogkerzz, 'same');
gzz = np.transpose(gzz, (1, 0, 2));
#print gzz[:, :, 0];
# Compute gradient products
gprrrr = np.multiply(grr, grr);
#print gprrrr[:, :, 0];
gprrcc = np.multiply(grr, gcc);
#print gprrcc[:, :, 0];
gprrzz = np.multiply(grr, gzz);
#print gprrzz[:, :, 0]
gpcccc = np.multiply(gcc, gcc);
gpcczz = np.multiply(gcc, gzz);
gpzzzz = np.multiply(gzz, gzz);
# Compute gradient amplitudes
# print ga.dtype;
ga = np.sqrt(gprrrr + gpcccc + gpzzzz);
#print ga[:, :, 0];
print "GA SHAPE:"
print 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);
gv = np.divide(gv, np.tile(ga[..., None], [1, 1, 1, 3]));
#print gv[:, :, 0, 1];
print "GV SHAPE:"
print gv.shape;
# 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 = np.single(gaussgen([gausigma, gausigma, gausigma]));
#print gaussker[:, :, 0];
print "Blurring gradient products..."
gprrrrgauss = signal.fftconvolve(gprrrr, gaussker, "same");
gprrccgauss = signal.fftconvolve(gprrcc, gaussker, "same");
gprrzzgauss = signal.fftconvolve(gprrzz, gaussker, "same");
gpccccgauss = signal.fftconvolve(gpcccc, gaussker, "same");
gpcczzgauss = signal.fftconvolve(gpcczz, gaussker, "same");
gpzzzzgauss = signal.fftconvolve(gpzzzz, gaussker, "same");
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, "dogsigma_" + str(jj) + "gausigma_" + str(kk) + 'tensorfsl.nii');
print 'Complete!'
In [41]:
## Compare tensor_fsl_data with imported values from MATLAB
print tensorfsl.shape # numpy
In [13]:
MATLAB_output = nib.load("/Users/Tony/Documents/Git Folder/seelviz/Tony/ipynb/test_MATLAB_10by10by10_outputs/dogsig1_gausig2.3/test_MATLAB_tensorfsl_dogsig1_gausig2.3.nii")
In [30]:
print MATLAB_output.shape
In [31]:
MATLAB_np_array = MATLAB_output.get_data()
print MATLAB_np_array.shape
In [32]:
print tensorfsl[:, :, :, 0];
In [33]:
print MATLAB_np_array[:, :, :, 0];
In [34]:
truth_boolean = np.isclose(tensorfsl, MATLAB_np_array)
print truth_boolean;
In [ ]: