In [1]:
import os
In [2]:
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from ndreg import *
import ndio.remote.neurodata as neurodata
import nibabel as nb
In [3]:
refToken = "ara_ccf2"
refImg = imgDownload(refToken)
In [4]:
refAnnoImg = imgDownload(refToken, channel="annotation")
In [5]:
randValues = np.random.rand(1000,3)
randValues = np.concatenate(([[0,0,0]],randValues))
randCmap = matplotlib.colors.ListedColormap (randValues)
In [6]:
inToken = "Aut1367"
nd = neurodata()
inImg = imgDownload(inToken, resolution=5)
In [7]:
rawestImg = sitk.GetArrayFromImage(inImg) ##convert to simpleITK image to normal numpy ndarray
In [8]:
print rawestImg.shape
In [9]:
rawestImg = sitk.GetArrayFromImage(inImg) ##convert to simpleITK image to normal numpy ndarray
In [ ]:
## Randomly sample 100k points
import random
x = rawestImg[:,0,0]
y = rawestImg[0,:,0]
z = rawestImg[0,0,:]
# mod by x to get x, mod by y to get y, mod by z to get z
xdimensions = len(x)
ydimensions = len(y)
zdimensions = len(z)
# index = random.sample(xrange(0,xdimensions*ydimensions*zdimensions),100000) #66473400 is multiplying xshape by yshape by zshape
#### random sampling of values > 250 ######
# X_val = []
# Y_val = []
# Z_val = []
# for i in index:
# z_val = int(i/(xdimensions*ydimensions))
# z_rem = i%(xdimensions*ydimensions)
# y_val = int(z_rem/xdimensions)
# x_val = int(z_rem/ydimensions)
# X_val.append(x_val)
# Y_val.append(y_val)
# Z_val.append(z_val)
# xlist = []
# ylist = []
# zlist = []
xyz = []
# for i in range(100000):
# value = 0
# while(value == 0):
# xval = random.sample(xrange(0,xdimensions), 1)[0]
# yval = random.sample(xrange(0,ydimensions), 1)[0]
# zval = random.sample(xrange(0,zdimensions), 1)[0]
# value = rawestImg[xval,yval,zval]
# if [xval, yval, zval] not in xyz and value > 250:
# xyz.append([xval, yval, zval])
# else:
# value = 0
# print xyz
#### obtaining 100000 brightest points #######
total = 100000
bright = np.max(rawestImg)
allpoints = []
brightpoints = []
for (x,y,z),value in np.ndenumerate(rawestImg):
# print (x,y,z)
# print (x,z,z)[0]
# print value
if value == bright:
brightpoints.append([(x,y,z)])
else:
allpoints.append([(x,y,z),value])
total -= len(brightpoints)
bright -= 1
if total < 0:
index = random.sample(xrange(0, len(brightpoints)), total)
for ind in index:
xyz.append(brightpoints[ind])
else:
xyz = xyz + brightpoints
# for item in brightpoints:
# outfile.write(",".join(item) + "\n")
while(total > 0):
del brightpoints[:]
for item in allpoints:
if item[3] == bright:
brightpoins.append(item)
allpoints.remove(item)
total -= len(brightpoints)
bright -= 1
if total < 0:
index = random.sample(xrange(0, len(brightpoints)), total)
for ind in index:
xyz.append(brightpoints[ind])
else:
xyz = xyz + brightpoints
# for item in brightpoints:
# outfile.write(",".join(item) + "\n")
print len(xyz)
In [ ]:
from plotly.offline import download_plotlyjs, init_notebook_mode, iplot
from plotly import tools
import plotly
plotly.offline.init_notebook_mode()
import plotly.graph_objs as go
# x = X_val
# y = Y_val
# z = Z_val
xlist = []
ylist = []
zlist = []
for i in xyz:
xlist.append(i[0])
ylist.append(i[1])
zlist.append(i[2])
x = xlist
y = ylist
z = zlist
trace1 = go.Scatter3d(
x = x,
y = y,
z = z,
mode='markers',
marker=dict(
size=1.2,
color='purple', # set color to an array/list of desired values
colorscale='Viridis', # choose a colorscale
opacity=0.15
)
)
data = [trace1]
layout = go.Layout(
margin=dict(
l=0,
r=0,
b=0,
t=0
)
)
fig = go.Figure(data=data, layout=layout)
print('Aut1367rawest' + "plotly")
iplot(fig)
In [ ]:
In [ ]:
plotly.offline.plot(fig, filename= 'AUT1367rawestluke200' + "plotly.html")
In [ ]:
## Clean out noise (Filter Image)
(values, bins) = np.histogram(sitk.GetArrayFromImage(inImg), bins=100, range=(0,500))
plt.plot(bins[:-1], values)
counts = np.bincount(values)
maximum = np.argmax(counts)
lowerThreshold = 100 #maximum
upperThreshold = sitk.GetArrayFromImage(inImg).max()+1
filteredImg = sitk.Threshold(inImg,lowerThreshold,upperThreshold,lowerThreshold) - lowerThreshold
In [ ]:
filterImg = sitk.GetArrayFromImage(filteredImg) ##convert to simpleITK image to normal numpy ndarray
In [ ]:
print filterImg[0][0]
In [ ]:
## Randomly sample 100k points after filtering
x = filterImg[:,0,0]
y = filterImg[0,:,0]
z = filterImg[0,0,:]
# mod by x to get x, mod by y to get y, mod by z to get z
xdimensions = len(x)
ydimensions = len(y)
zdimensions = len(z)
index = random.sample(xrange(0,xdimensions*ydimensions*zdimensions),100000) #66473400 is multiplying xshape by yshape by zshape
X_val = []
Y_val = []
Z_val = []
for i in index:
z_val = int(i/(xdimensions*ydimensions))
z_rem = i%(xdimensions*ydimensions)
y_val = int(z_rem/xdimensions)
x_val = int(z_rem/ydimensions)
X_val.append(x_val)
Y_val.append(y_val)
Z_val.append(z_val)
xlist = []
ylist = []
zlist = []
for i in range(100000):
xlist.append(random.sample(xrange(0,xdimensions), 1)[0])
ylist.append(random.sample(xrange(0,ydimensions), 1)[0])
zlist.append(random.sample(xrange(0,zdimensions), 1)[0])
print xlist
In [ ]:
print X_val
print Y_val
print Z_val
print len(X_val)
print len(Y_val)
print len(Z_val)
UNUSED:
spacingImg = inImg.GetSpacing()
spacing = tuple(i * 50 for i in spacingImg)
print spacingImg
print spacing
inImg.SetSpacing(spacingImg)
inImg_download = inImg # Aut1367 set to default spacing
inImg = imgResample(inImg, spacing=refImg.GetSpacing())
Img_reorient = imgReorient(inImg, "LPS", "RSA") #specific reoriented Aut1367
inImg_reorient = Img_reorient
refImg_ds = imgResample(refImg, spacing=spacing) # atlas downsample 50x
In [ ]:
from plotly.offline import download_plotlyjs, init_notebook_mode, iplot
from plotly import tools
import plotly
plotly.offline.init_notebook_mode()
import plotly.graph_objs as go
x = X_val
y = Y_val
z = Z_val
trace1 = go.Scatter3d(
x = x,
y = y,
z = z,
mode='markers',
marker=dict(
size=1.2,
color='purple', # set color to an array/list of desired values
colorscale='Viridis', # choose a colorscale
opacity=0.15
)
)
data = [trace1]
layout = go.Layout(
margin=dict(
l=0,
r=0,
b=0,
t=0
)
)
fig = go.Figure(data=data, layout=layout)
print('Aut1367filter' + "plotly")
plotly.offline.plot(fig, filename= 'AUT1367filter' + "plotly.html")
(values, bins) = np.histogram(filterImg, bins=100, range=(0,500)) plt.plot(bins[:-1], values)
counts = np.bincount(values) maximum = np.argmax(counts)
lowerThreshold = maximum upperThreshold = filterImg.max()+1
filterX2Img = sitk.Threshold(inImg,lowerThreshold,upperThreshold,lowerThreshold) - lowerThreshold
newImg = sitk.GetArrayFromImage(filterX2Img) ##convert to simpleITK image to normal numpy ndarray
In [ ]:
## Histogram Equalization
## Cut from generateHistogram from clarityviz
import cv2
im = filterImg
img = im[:,:,:]
shape = im.shape
#affine = im.get_affine()
x_value = shape[0]
y_value = shape[1]
z_value = shape[2]
#####################################################
imgflat = img.reshape(-1)
#img_grey = np.array(imgflat * 255, dtype = np.uint8)
#img_eq = exposure.equalize_hist(img_grey)#new_img = img_eq.reshape(x_value, y_value, z_value)
#globaleq = nib.Nifti1Image(new_img, np.eye(4))
#nb.save(globaleq, '/home/albert/Thumbo/AutAglobaleq.nii')
######################################################
#clahe = cv2.createCLAHE(clipLimit=2.0, tileGridSize=(8,8))
clahe = cv2.createCLAHE(clipLimit=2.0, tileGridSize=(8,8))
img_grey = np.array(imgflat * 255, dtype = np.uint8)
#threshed = cv2.adaptiveThreshold(img_grey, 255, cv2.ADAPTIVE_THRESH_MEAN_C, cv2.THRESH_BINARY, 3, 0)
cl1 = clahe.apply(img_grey)
#cv2.imwrite('clahe_2.jpg',cl1)
#cv2.startWindowThread()
#cv2.namedWindow("adaptive")
#cv2.imshow("adaptive", cl1)
#cv2.imshow("adaptive", threshed)
#plt.imshow(threshed)
localimgflat = cl1 #cl1.reshape(-1)
newer_img = localimgflat.reshape(x_value, y_value, z_value) ##this is the numpy.ndarray version
localeq = nb.Nifti1Image(newer_img, np.eye(4))
In [ ]:
print newer_img.shape
In [ ]:
x_histeq = newer_img[:,0,0]
y_histeq = newer_img[0,:,0]
z_histeq = newer_img[0,0,:]
## Randomly sample 100k points after filtering
xdimensions = len(x)
ydimensions = len(y)
zdimensions = len(z)
index = random.sample(xrange(0,xdimensions*ydimensions*zdimensions),100000) #66473400 is multiplying xshape by yshape by zshape
X_val = []
Y_val = []
Z_val = []
for i in index:
z_val = int(i/(xdimensions*ydimensions))
z_rem = i%(xdimensions*ydimensions)
y_val = int(z_rem/xdimensions)
x_val = int(z_rem/ydimensions)
X_val.append(x_val)
Y_val.append(y_val)
Z_val.append(z_val)
xlist = []
ylist = []
zlist = []
for i in range(100000):
xlist.append(random.sample(xrange(0,xdimensions), 1)[0])
ylist.append(random.sample(xrange(0,ydimensions), 1)[0])
zlist.append(random.sample(xrange(0,zdimensions), 1)[0])
for i in index:
z_val = int(i/(xdimensions*ydimensions))
z_rem = i%(xdimensions*ydimensions)
y_val = int(z_rem/xdimensions)
x_val = int(z_rem/ydimensions)
X_val.append(x_val)
Y_val.append(y_val)
Z_val.append(z_val)
trace1 = go.Scatter3d(
x = x,
y = y,
z = z,
mode='markers',
marker=dict(
size=1.2,
color='purple', # set color to an array/list of desired values
colorscale='Viridis', # choose a colorscale
opacity=0.15
)
)
data = [trace1]
layout = go.Layout(
margin=dict(
l=0,
r=0,
b=0,
t=0
)
)
fig = go.Figure(data=data, layout=layout)
print('Aut1367postfilterHistEq' + "plotly")
plotly.offline.plot(fig, filename= 'Aut1367postfilterHistEq' + "plotly.html")
In [ ]:
print "done"
In [ ]: