In [2]:
## Following MATLAB code from http://capture-clarity.org/clarity-based-tractography/
In [3]:
import os
In [4]:
## 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 [5]:
pwd
Out[5]:
In [6]:
import numpy as np
In [7]:
import math
from scipy import ndimage
import nibabel as nib
In [8]:
## Change later on
file_path = "/Users/Tony/Documents/Git Folder/seelviz/Tony/ipynb/TIFF_stack"
directory = os.path.dirname(file_path)
In [9]:
cd TIFF_stack/
In [10]:
print len([name for name in os.listdir('.') if os.path.isfile(name)])
In [11]:
from PIL import Image
from scipy import signal
In [21]:
'''
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 [22]:
'''
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 [20]:
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('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.)
img_data = np.asarray(im); #data is hard coded to be np.ones
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
dogkerrr = np.transpose(dogkercc, (1, 0, 2));
dogkerzz = np.transpose(dogkercc, (0, 2, 1));
# Compute gradients
grr = signal.fftconvolve(img_data, dogkercc, 'same');
grr = np.transpose(grr, (1, 0, 2));
gcc = signal.fftconvolve(img_data, dogkerrr, 'same');
gcc = np.transpose(gcc, (1, 0, 2));
gzz = signal.fftconvolve(img_data, dogkerzz, 'same');
gzz = np.transpose(gzz, (1, 0, 2));
# Compute gradient products
gprrrr = np.multiply(grr, grr);
gprrcc = np.multiply(grr, gcc);
gprrzz = np.multiply(grr, gzz);
gpcccc = np.multiply(gcc, gcc);
gpcczz = np.multiply(gcc, gzz);
gpzzzz = np.multiply(gzz, gzz);
# Compute gradient amplitudes
ga = np.sqrt(gprrrr + gpcccc + gpzzzz);
print "Finished computing gradient amplitudes, saving gradient amplitude image..."
# 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 "Finished computing gradient vectors, saving gradient vector image..."
# 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 "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!'