import matplotlib
import SimpleITK as sitk
import matplotlib.pyplot as plt
from matplotlib.widgets import Slider
import numpy as np
from stl import mesh
from mpl_toolkits import mplot3d
import cv2
import scipy.ndimage as ndimage
import pandas
import timeit
from sympy.geometry import *
from math import atan2
from sympy import plot
# my_mesh = mesh.Mesh.from_file('/home/morrile2/Documents/MULTIS_test/MRI_Segmentation/Muscle_verycoarse.stl')
reader = sitk.ImageFileReader()
reader.SetFileName('/home/morrile2/Documents/MULTIS Data/MRI_Segmentation/Femur_T1.nii')
img_original = reader.Execute()
img_T1 = sitk.Cast(sitk.RescaleIntensity(img_original), sitk.sitkUInt8)
img_o = sitk.GetArrayFromImage(img_T1)
reader2 = sitk.ImageFileReader()
reader2.SetFileName('/home/morrile2/Documents/MULTIS Data/MRI_Segmentation/Femur_T1-label.nii')
img_label = reader2.Execute()
img_l = sitk.GetArrayFromImage(img_label)
img_combined = sitk.LabelOverlay(img_T1, img_label, opacity = 0.8)
img_c = sitk.GetArrayFromImage(img_combined)
img_origin = img_original.GetOrigin()
img_spacing = img_original.GetSpacing()
img_dim = img_original.GetSize()
# array = np.empty()
item = []
f = open('/home/morrile2/Documents/MULTIS Data/MRI_Segmentation/Muscle_verycoarse.dat')
d = f.readlines()
item = ([x.split() for x in d])
f.close()
def column(matrix, i):
return [row[i] for row in matrix]
array = np.array(item)
for p, x in enumerate(array):
for i, n in enumerate(x):
array[p][i] = float(n)
# Create a new plot
figure = plt.figure()
ax = figure.add_subplot(111,projection='3d')
array_nodes = array[1:int(array[0][0])+1]
x = column(array_nodes,1)
y = column(array_nodes,2)
z = column(array_nodes,3)
# print(len(x))
array_elements = array[int(array[0][0]+1):int(array[0][0]+array[0][1]+1)]
def get_voxels(e_num, IMG):
# print(e_num)
# voxels = []
img_elem = []
vox_all = []
# Getting voxels that are within a specific element
n1 = int(array_elements[e_num-1][2])
n2 = int(array_elements[e_num-1][3])
n3 = int(array_elements[e_num-1][4])
n4 = int(array_elements[e_num-1][5])
# ax.scatter([x[n1 - 1], x[n2 - 1], x[n3 - 1], x[n4 - 1]],
# [y[n1 - 1], y[n2 - 1], y[n3 - 1], y[n4 - 1]],
# [z[n1 - 1], z[n2 - 1], z[n3 - 1], z[n4 - 1]])
#
# ax.set_ylim([-120, -25.3])
# ax.set_xlim([-3, 124])
# ax.set_zlim([-251, -52])
# ax.set_title('3D view of element '+ str(e_num))
# ax.hold(True)
nodes = [n1, n2, n3, n4]
# print(nodes)
x_range = []
y_range = []
z_range = []
for n in nodes:
x_range.append(array_nodes[int(n-1)][1])
y_range.append(array_nodes[int(n-1)][2])
z_range.append(array_nodes[int(n-1)][3])
# Grid setup for node points (boundaries of image voxels)
X = np.linspace(-img_origin[0]+img_spacing[0]/2, -img_origin[0]-img_spacing[0]/2+(-(img_dim[0]-1)*img_spacing[0]), num = img_dim[0]+1)
Y = np.linspace(-img_origin[1]-img_spacing[1]/2+(-(img_dim[1]-1)*img_spacing[1]), -img_origin[1]+img_spacing[1]/2, num = img_dim[1]+1)
Z = np.linspace(img_origin[2]-img_spacing[2]/2, img_origin[2]+(img_dim[2]-1)*img_spacing[2]+img_spacing[2]/2, num=img_dim[2]+1)
# Grid setup for center points (follow image voxels)
XC = np.linspace(-img_origin[0], -img_origin[0]+(-(img_dim[0]-1)*img_spacing[0]), num = img_dim[0])
YC = np.linspace(-img_origin[1]+(-(img_dim[1]-1)*img_spacing[1]), -img_origin[1], num = img_dim[1])
ZC = np.linspace(img_origin[2], img_origin[2]+(img_dim[2]-1)*img_spacing[2], num=img_dim[2])
# print(img_dim, img_spacing)
z_planes = (np.where(np.logical_and(Z>=min(z_range), Z<=max(z_range))))
xx, yy = np.meshgrid(X, Y)
for p in z_planes[0]:
# for p in [41]: #Testing only
n_above = []
n_below = []
for i, n in enumerate(z_range):
if n>Z[p]:
n_above.append(i)
elif n<Z[p]:
n_below.append(i)
else:
print("Error")
# Parametric definition of the lines above and below the z plane to find where the tetrahedral intersects with the z planes
x_intersect = []
y_intersect = []
z_intersect = []
voxels = []
for i1 in n_above:
for i2 in n_below:
vector = [x_range[i1]-x_range[i2], y_range[i1]-y_range[i2], z_range[i1]-z_range[i2]]
t = (Z[p]-z_range[i1])/vector[2]
x_intersect.append(x_range[i1]+vector[0]*t)
y_intersect.append(y_range[i1]+vector[1]*t)
z_intersect.append(Z[p])
# Check for intersections in polygon intersection with z-plane
# Remove line segments going through the middle of the polygon (not a boundary of the polygon)
if len(x_intersect)>3:
poly_lines = []
for n1 in range(len(x_intersect)):
for n2 in range(n1 + 1, len(x_intersect), 1):
poly_lines.append(
Segment(Point(x_intersect[n1], y_intersect[n1]), Point(x_intersect[n2], y_intersect[n2])))
poly_sides = []
for l1 in range(len(poly_lines)):
for l2 in range(l1+1, len(poly_lines)):
intersect = poly_lines[l1].intersect(poly_lines[l2])
if len(intersect) == 0:
poly_sides.append(poly_lines[l1])
poly_sides.append(poly_lines[l2])
# print(poly_lines[l1], poly_lines[l2])
# Keep all line segments for triangle intersection with z-plane
elif len(x_intersect)==3:
poly_sides = []
for n1 in range(len(x_intersect)):
for n2 in range(n1 + 1, len(x_intersect), 1):
poly_sides.append(Segment(Point(x_intersect[n1], y_intersect[n1]), Point(x_intersect[n2], y_intersect[n2])))
# print(len(poly_sides))
y_lines = (np.where(np.logical_and(Y>=min(y_intersect), Y<=max(y_intersect))))
# x_i, y_i = np.meshgrid([min(x_intersect), max(x_intersect)], [min(y_intersect), max(y_intersect)])
# ax.plot_surface(x_i, y_i, 0*x_i+0*y_i+z_intersect[0], alpha = 0.4)
#
# plt.figure(figsize=(15,11))
# plt.scatter(x_intersect, y_intersect)
boundaries = []
for YY in y_lines[0]:
x_horiz = []
y_horiz = []
for a in range(len(x_intersect)):
for b in range(a+1,len(x_intersect),1):
x_line_intersect = x_intersect[a]+(Y[YY]-y_intersect[a])/(y_intersect[a]-y_intersect[b])*(x_intersect[a]-x_intersect[b])
if x_line_intersect > min([x_intersect[a], x_intersect[b]]) and x_line_intersect < max([x_intersect[a], x_intersect[b]]):
x_horiz.append(x_line_intersect)
y_horiz.append(Y[YY])
x_points = np.where(np.logical_and(X > min(x_horiz), X < max(x_horiz)))
bound_info = [[Y[YY], min(x_horiz), max(x_horiz)]]
boundaries = boundaries+bound_info
# plt.scatter(X[x_points[0]], Y[x_points[0]], color = 'g')
# print(y_lines, x_points)
for XX in x_points[0]:
# voxels.append(IMG[p, YY, XX])
if (any([p,YY,XX]==x for x in voxels)) == False:
voxels.append([p,YY,XX])
if (any([p, YY-1, XX] == x for x in voxels)) == False:
voxels.append([p, YY-1, XX])
if (any([p, YY, XX-1] == x for x in voxels)) == False:
voxels.append([p, YY, XX-1])
if (any([p, YY-1, XX-1] == x for x in voxels)) == False:
voxels.append([p, YY-1, XX-1])
img_elem.append([XX, YY])
h = plt.scatter(X[XX],Y[YY], color='g')
# plt.scatter(XC[XX], YC[YY], color='y')
# plt.scatter(XC[XX-1], YC[YY], color='y')
# plt.scatter(XC[XX], YC[YY-1], color='y')
# plt.scatter(XC[XX-1], YC[YY-1], color='y')
# plt.scatter([min(x_horiz), max(x_horiz)], [Y[YY], Y[YY]], color='r')
# # plt.legend(['Polygon', 'Voxel centers', 'X-Boundaries'], loc='upper left', bbox_to_anchor=(1.05,1))
# plt.subplots_adjust(right=0.8)
# plt.title('2D-Slice at Z = '+ str('%0.2f' %Z[p]) + ' mm (Slice ' + str(p) + ')')
# print(img_elem)#node points of image within boundary of slice
count_full = 0
count_partial = 0
boundaries = np.array(boundaries).transpose()
# print(voxels)
for pi, ii in enumerate(voxels):
if (any([ii[2], ii[1]]==x for x in img_elem)) and (any([ii[2]+1, ii[1]]==x for x in img_elem)) and (any([ii[2], ii[1]+1]==x for x in img_elem)) and (any([ii[2]+1, ii[1]+1]==x for x in img_elem)):
area_ratio = 1.0
count_full +=1
else:
P1 = Point(XC[ii[2]]-.5*img_spacing[0], YC[ii[1]]-0.5*img_spacing[0])
P2 = Point(XC[ii[2]]+.5*img_spacing[0], YC[ii[1]]-0.5*img_spacing[0])
P3 = Point(XC[ii[2]]+.5*img_spacing[0], YC[ii[1]]+0.5*img_spacing[0])
P4 = Point(XC[ii[2]]-.5*img_spacing[0], YC[ii[1]]+0.5*img_spacing[0])
# P1.astype(float)
# print(P1)
pp = [P1, P2, P3, P4]
# square = Polygon(P1, P2, P3, P4)
# print('******************')
square = Polygon(*pp)
for num, s in enumerate(poly_sides):
intersections = square.intersection(s)
if len(intersections) != 0:
for pts in intersections:
pp.append(Point(pts[0], pts[1]))
# Add any points from the polygon that are within the voxel (2D)
for p in range(len(x_intersect)):
if (square.encloses(Point(x_intersect[p], y_intersect[p]))) == True:
pp.append(Point(x_intersect[p], y_intersect[p]))
final_pp = []
# Remove any node points that are not in the z slice polygon
original_nodes = []
for i, n in enumerate(pp):
if i<4:
try:
ind_y = np.where(np.logical_and(boundaries[0] > n[1]-0.00001, boundaries[0] < n[1]+0.00001))
if n[0] > boundaries[1][ind_y[0]] and n[0] < boundaries[2][ind_y[0]]:
original_nodes.append(n)
final_pp.append(n)
except:
pass
else:
final_pp.append(n)
Cent = GetCenterPoint(final_pp)
# plotCenterPoints(Cent)
angles = []
for P in final_pp:
angles.append(atan2(P[1]-Cent[1], P[0]-Cent[0])*180/3.14)
sorted_info = [i[1] for i in sorted(enumerate(zip(angles, final_pp)), key=lambda x: x[1])]
angles_sort, final_pp_sort = zip(*sorted_info)
partial_vox = Polygon(*final_pp_sort)
area_ratio = float(partial_vox.area/square.area)
count_partial += 1
voxels[pi] = [img_l[ii[0], ii[1], ii[2]], area_ratio]
vox_all = vox_all + voxels
print(count_full+count_partial, vox_all)
# plt.show()
def plotCenterPoints(points):
# print(points)
plt.scatter(points[0], points[1], color='k')
def GetCenterPoint(points):
pointsSum = [0, 0]
for p in points:
pointsSum[0] += p[0]
pointsSum[1] += p[1]
pointsSum[0]/=len(points)
pointsSum[1]/=len(points)
return pointsSum
s_time = timeit.default_timer()
# array_elements = [array_elements[16686]]
for ii in [array_elements[16686]]:
if len(ii) == 6:
get_voxels(int(ii[0]), img_o)
elapsed = timeit.default_timer()-s_time
print(elapsed)
# sitk.WriteImage(img_ROI, '/home/morrile2/Documents/MULTIS_test/MRI_Segmentation/Femur_T1-label_small.nii')
# def show2Dslices(img):
# fig, ax = plt.subplots()
# plt.subplots_adjust(left=0.25, bottom=0.25)
#
# # print(-img_origin[1]+img_spacing[1]+(-img_dim[1]*img_spacing[1]), -img_origin[1]-img_spacing[1])
# l = plt.imshow(img[0,:,:], cmap='gray', extent=[-img_origin[0]-img_spacing[0]/2, -img_origin[0]+img_spacing[0]/2+(-img_dim[0]*img_spacing[0]), -img_origin[1]+img_spacing[1]+(-img_dim[1]*img_spacing[1]), -img_origin[1]-img_spacing[1]])
#
# # ax.autoscale(False)
#
# axcolor = 'lightgoldenrodyellow'
# axImgIdx = plt.axes([0.25, 0.1, 0.65, 0.03], axisbg=axcolor)
#
# # print(img.shape)
# sImgIdx = Slider(axImgIdx, 'Img Index', 0, img.shape[0]-1, valinit=int(0), valfmt='%d')
#
# zPos = ax.annotate('z = '+ str(Z[int(sImgIdx.val)]), (0.01, .01), (0.01, .01), xycoords='figure fraction')
#
# def update(val):
# value = int(sImgIdx.val)
# l.set_data(img[value,:,:])
# zPos.set_text('z = '+ str(Z[int(sImgIdx.val)]))
# fig.canvas.draw_idle()
#
# sImgIdx.on_changed(update)
# plt.show()
# show2Dslices(img_c)