In [74]:
import tractography_latest as tract
reload(tract)
Out[74]:
In [2]:
test = tract.tiff_stack_to_array('/home/tractography_test/CTT/demo/data/')
In [41]:
from scipy import io
In [43]:
stackmat = io.loadmat('stack_image.mat')
In [51]:
np_stackmat = stackmat['stack']
In [52]:
print np_stackmat.shape
In [53]:
dogsigmaArr = [1]
gausigmaArr = [0.5]
img_data = np_stackmat
In [55]:
import os
import numpy as np
import math
import SimpleITK as sitk
from scipy import ndimage
import nibabel as nib
from PIL import Image
import scipy.misc
from scipy import signal
In [151]:
# Helper function similar to MATLAB
def cropvol(vol, sz):
print np.size(vol);
vol[0:sz, :, :, :] = 0;
vol[-1-sz+1:-1, :, :, :] = 0;
vol[:, 0:sz, :, :] = 0;
vol[:, -1-sz+1:-1, :, :] = 0;
vol[:, :, 0:sz, :] = 0;
vol[:, :, -1-sz+1:-1, :] = 0;
return vol;
for ii in range(len(dogsigmaArr)):
dogsigma = dogsigmaArr[ii];
print "Start DoG Sigma on " + str(dogsigma);
# Generate dog kernels
dogkercc = tract.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];
for jj in range(len(gausigmaArr)):
gausigma = gausigmaArr[jj];
print "Start Gauss Sigma with gausigma = " + str(gausigma);
print "Generating Gaussian kernel..."
gaussker = np.single(tract.gaussgen([gausigma, gausigma, gausigma]));
# Generate half size kernel
halfsz = (max(len(dogkercc), len(gaussker)) + 1) / 2;
# Compute gradients
grr = signal.convolve((img_data).astype('float64'), dogkerrr.astype('float64'), 'same');
grr.astype('float64');
#print grr[:, :, 0];
gcc = signal.convolve(img_data.astype('float64'), dogkercc.astype('float64'), 'same');
gcc.astype('float64');
#print gcc[:, :, 0];
gzz = signal.convolve(img_data.astype('float64'), dogkerzz.astype('float64'), 'same');
gzz.astype('float64');
#print gzz[:, :, 0];
# Compute gradient products
gprrrr = np.multiply(grr.astype('float64'), grr.astype('float64'));
#print gprrrr[:, :, 0];
gprrcc = np.multiply(grr.astype('float64'), gcc.astype('float64'));
#print gprrcc[:, :, 0];
gprrzz = np.multiply(grr.astype('float64'), gzz.astype('float64'));
#print gprrzz[:, :, 0]
gpcccc = np.multiply(gcc.astype('float64'), gcc.astype('float64'));
gpcczz = np.multiply(gcc.astype('float64'), gzz.astype('float64'));
gpzzzz = np.multiply(gzz.astype('float64'), gzz.astype('float64'));
# Compute gradient amplitudes
# print ga.dtype;
ga = np.sqrt(gprrrr.astype('float64') + gpcccc.astype('float64') + gpzzzz.astype('float64'));
#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.gz');
# Compute gradient vectors
gv = np.concatenate((grr[..., np.newaxis], gcc[..., np.newaxis], gzz[..., np.newaxis]), axis = 3);
#print gv[:, :, 0, 0];
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.gz');
#print gaussker[:, :, 0];
print "Blurring gradient products..."
gprrrrgauss = signal.convolve(gprrrr, gaussker, "same");
#print gprrrrgauss[:, :, 0];
gprrccgauss = signal.convolve(gprrcc, gaussker, "same");
#print gprrccgauss[:, :, 0];
gprrzzgauss = signal.convolve(gprrzz, gaussker, "same");
gpccccgauss = signal.convolve(gpcccc, gaussker, "same");
gpcczzgauss = signal.convolve(gpcczz, gaussker, "same");
gpzzzzgauss = signal.convolve(gpzzzz, gaussker, "same");
In [152]:
print gzz[0:9, 0:9, 0]
In [153]:
print gprrrr[0:9, 0:9, 0]
In [154]:
print grr[0:9, 0:9, 0]
In [155]:
print (grr[0, 0, 0]).astype('uint64') * (grr[0, 0, 0]).astype('uint64')
In [156]:
print ga[0:9, 0:9, 0]
In [157]:
reload(tract)
fsl, dtk = tract.generate_FSL_and_DTK_structure_tensor(np_stackmat, 'test_data', dogsigmaArr=[1], gausigmaArr=[0.5]);
In [9]:
mask = tract.tiff_stack_to_nii('/home/tractography_test/CTT/demo/mask-brain/', 'brainmask')
In [120]:
import nibabel as nib
MATLAB_output = nib.load("/home/tractography_test/CTT/demo/result/dog1gau0.5/dtk_tensor.nii.gz")
In [121]:
MATLAB_np_array = MATLAB_output.get_data()
In [122]:
print dtk.shape
In [158]:
dtk[:, :, :, 0]
Out[158]:
In [130]:
MATLAB_np_array[:, :, :, 0]
Out[130]:
In [169]:
import numpy as np
truth_boolean = np.isclose(dtk, MATLAB_np_array, rtol = 1e-4)
In [170]:
print MATLAB_np_array.shape
In [171]:
correct_number = np.sum(truth_boolean == True); # Total possible = 4,374,000
print correct_number
In [173]:
## Percentage correct (at the 1e-4 value)
print correct_number / (4374000.0)
In [ ]: