Skip to content
Snippets Groups Projects
Commit eea220a8 authored by Joshua Steer's avatar Joshua Steer
Browse files

Initial commit

parents
No related branches found
No related tags found
No related merge requests found
*.pyc
/*.egg-info
/*.egg
*.STL
AmpScan/*.pyc
\ No newline at end of file
import sys
import numpy as np
from core import AmpObject
from registration import regObject
from ampVis import qtVtkWindow
from pressSens import pressSense
from PyQt5.QtCore import QPoint, QSize, Qt, QTimer, QRect, pyqtSignal
from PyQt5.QtGui import (QColor, QFontMetrics, QImage, QPainter, QIcon,
QOpenGLVersionProfile)
from PyQt5.QtWidgets import (QAction, QApplication, QGridLayout,
QMainWindow, QMessageBox,
QOpenGLWidget, QFileDialog,
QSlider, QWidget)
class AmpScanGUI(QMainWindow):
def __init__(self, parent = None):
super(AmpScanGUI, self).__init__()
self.vtkWidget = qtVtkWindow()
self.mainWidget = QWidget()
self.AmpObj = None
self.CMap = np.array([[212.0, 221.0, 225.0],
[31.0, 73.0, 125.0]])/255.0
self.setCentralWidget(self.mainWidget)
self.createActions()
self.createMenus()
self.Layout = QGridLayout()
self.Layout.addWidget(self.vtkWidget, 0, 0)
self.mainWidget.setLayout(self.Layout)
self.setWindowTitle("AmpScan Visualiser")
self.resize(800, 800)
self.show()
def chooseOpenFile(self):
self.fname = QFileDialog.getOpenFileName(self, 'Open file',
filter="Meshes (*.stl)")
if self.AmpObj is not None:
self.vtkWidget.renderActors(self.AmpObj.actors, [])
self.AmpObj = AmpObject(self.fname[0], 'limb')
self.AmpObj.addActor(stype='limb')
self.AmpObj.lp_smooth(stype='limb')
self.vtkWidget.setnumViewports(1)
self.vtkWidget.setProjection()
self.vtkWidget.renderActors(self.AmpObj.actors, ['limb',])
def chooseSocket(self):
self.sockfname = QFileDialog.getOpenFileName(self, 'Open file',
filter="Meshes (*.stl)")
self.AmpObj.addData(self.sockfname[0], stype='socket')
self.AmpObj.addActor(stype='socket')
self.AmpObj.lp_smooth(stype='socket')
def align(self):
self.vtkWidget.setnumViewports(2)
self.vtkWidget.setView(view=[-1, 0, 0], viewport=1)
self.vtkWidget.setProjection(True, 0)
self.vtkWidget.setProjection(True, 1)
# self.vtkWidget.render(self.AmpObj.actors, dispActors=['limb',])
# self.vtkWidget.render(self.AmpObj.actors, dispActors=['socket',],
# viewport=1)
self.vtkWidget.renderActors(self.AmpObj.actors,
dispActors=['limb', 'socket'],
viewport=0)
self.vtkWidget.renderActors(self.AmpObj.actors,
dispActors=['limb', 'socket'],
viewport=1)
self.AmpObj.actors['limb'].setColor([1.0, 0.0, 0.0])
self.AmpObj.actors['limb'].setOpacity(0.5)
self.AmpObj.actors['socket'].setColor([0.0, 0.0, 1.0])
self.AmpObj.actors['socket'].setOpacity(0.5)
def register(self):
self.vtkWidget.setnumViewports(1)
self.vtkWidget.setProjection()
self.RegObj = regObject(self.AmpObj)
self.RegObj.registration(steps=5, baseline='socket', target='limb',
reg = 'reglimb', direct=True)
self.RegObj.addActor(stype='reglimb', CMap=self.CMap)
self.vtkWidget.renderActors(self.AmpObj.actors, ['reglimb',], shading=False)
self.vtkWidget.setScalarBar(self.AmpObj.actors['reglimb'])
def analyse(self):
self.RegObj.plot_slices()
def chooseFE(self):
FEname = QFileDialog.getOpenFileName(self, 'Open file',
filter="FE results (*.npy)")
self.vtkWidget.setnumViewports(1)
self.AmpObj.addFE([FEname[0],])
self.AmpObj.lp_smooth('FE', n=1)
self.AmpObj.addActor(stype='FE', CMap=self.CMap)
self.AmpObj.actors['FE'].setScalarRange(smin=0.0, smax=50)
self.vtkWidget.renderActors(self.AmpObj.actors, ['FE',])
self.vtkWidget.setScalarBar(self.AmpObj.actors['FE'])
def choosePress(self):
vName = QFileDialog.getOpenFileName(self, 'Open file',
filter="Sensor vertices (*.csv)")
pName = QFileDialog.getOpenFileName(self, 'Open file',
filter="Sensor pressures (*.csv)")
self.vtkWidget.setnumViewports(1)
self.pSense = pressSense()
self.pSense.calcFaces(d=5)
self.pSense.importVert(vName[0])
self.pSense.importPress(pName[0])
self.pSense.addActor(CMap=self.CMap)
self.AmpObj.actors['antS'] = self.pSense.actors['antS']
self.AmpObj.actors['socket'].setColor([1.0, 1.0, 1.0])
self.AmpObj.actors['socket'].setOpacity(1.0)
self.vtkWidget.renderActors(self.AmpObj.actors, ['socket', 'antS'])
self.vtkWidget.setScalarBar(self.AmpObj.actors['antS'])
def createActions(self):
self.openFile = QAction(QIcon('open.png'), 'Open', self,
shortcut='Ctrl+O',
triggered=self.chooseOpenFile)
self.openSocket = QAction(QIcon('open.png'), 'Open Socket', self,
triggered=self.chooseSocket)
self.openFE = QAction(QIcon('open.png'), 'Open FE', self,
triggered=self.chooseFE)
self.openPress = QAction(QIcon('open.png'), 'Open Press', self,
triggered=self.choosePress)
self.exitAct = QAction("E&xit", self, shortcut="Ctrl+Q",
triggered=self.close)
self.align = QAction(QIcon('open.png'), 'Align', self,
triggered=self.align)
self.rect = QAction(QIcon('open.png'), 'Register', self,
triggered=self.register)
self.analyse = QAction(QIcon('open.png'), 'Analyse', self,
triggered=self.analyse)
def createMenus(self):
self.fileMenu = self.menuBar().addMenu("&File")
self.fileMenu.addAction(self.openFile)
self.fileMenu.addAction(self.openSocket)
self.fileMenu.addSeparator()
self.fileMenu.addAction(self.exitAct)
self.alignMenu = self.menuBar().addMenu("&Align")
self.alignMenu.addAction(self.align)
self.regMenu = self.menuBar().addMenu("&Registration")
self.regMenu.addAction(self.rect)
self.feMenu = self.menuBar().addMenu("&FE Analysis")
self.feMenu.addAction(self.openFE)
self.analyseMenu = self.menuBar().addMenu("&Analyse")
self.analyseMenu.addAction(self.analyse)
self.kineticMenu = self.menuBar().addMenu("&Kinetic Measurements")
self.kineticMenu.addAction(self.openPress)
def runAmpScanGUI():
# if __name__ == "__main__":
app = QApplication(sys.argv)
mainWin = AmpScanGUI()
mainWin.show()
sys.exit(app.exec_())
\ No newline at end of file
# -*- coding: utf-8 -*-
"""
Created on Thu Dec 15 13:50:41 2016
@author: js22g12
"""
from .core import AmpObject
from .registration import regObject
from .AmpScanGUI import AmpScanGUI
from .socketDesignGUI import socketDesignGUI
# -*- coding: utf-8 -*-
"""
Created on Thu Sep 28 13:19:18 2017
@author: js22g12
Functions that deal with the visualisation of the limb and data
"""
import numpy as np
import vtk
from vtk.util import numpy_support
from vtk.qt.QVTKRenderWindowInteractor import QVTKRenderWindowInteractor
from PyQt5.QtWidgets import QWidget
class vtkRender(vtk.vtkRenderer):
"""
Minor modification to the vtkRenderer class to allow easy access to
the currently displayed actors within each vtkRenderer object
"""
def __init__(self):
super(vtkRender, self).__init__()
self.actors = []
class ampVTK(object):
"""
Methods for display of the AmpObj within a vtk window
"""
def __init__(self):
self.rens = []
self.cams = []
self.axes = []
self.scalar_bar = None
self.style = vtk.vtkInteractorStyleTrackballCamera()
self.cams.append(vtk.vtkCamera())
self.setView()
self.rens.append(vtkRender())
self.rens[0].SetBackground(0.1, 0.2, 0.4)
self.rens[0].SetActiveCamera(self.cams[0])
self.axes.append(vtk.vtkCubeAxesActor())
def renderActors(self, actors, dispActors=['limb'], viewport=0,
shading=True, zoom=1):
"""
Render the required AmpObj actor in the vtk viewport
"""
for actor in self.rens[viewport].actors:
self.rens[viewport].RemoveActor(actors[actor])
for actor in dispActors:
actors[actor].setShading(shading)
self.rens[viewport].AddActor(actors[actor])
self.rens[viewport].actors = dispActors
self.rens[viewport].ResetCamera()
self.cams[viewport].Zoom(zoom)
def setScalarBar(self, actor):
"""
Set scalar bar based upon lookup table
"""
if self.scalar_bar is not None:
self.rens[0].RemoveActor(self.scalar_bar)
self.scalar_bar = vtk.vtkScalarBarActor()
self.scalar_bar.SetLookupTable(actor.lut)
self.scalar_bar.SetOrientationToVertical()
self.scalar_bar.SetPosition(0.85, 0.7)
self.scalar_bar.SetPosition2(0.1, 0.3)
self.rens[0].AddActor(self.scalar_bar)
def setView(self, view = [0, -1, 0], viewport=0):
self.cams[viewport].SetPosition(view[0], view[1], view[2])
self.cams[viewport].SetViewUp(-0.0, 0.0, 1.0)
def setBackground(self, color=[0.1, 0.2, 0.4]):
for ren in self.rens:
ren.SetBackground(color)
def setProjection(self, perspective=False, viewport=0):
self.cams[viewport].SetParallelProjection(perspective)
def addAxes(self, actors, viewport=0, color = [1.0, 1.0, 1.0], font=None):
lim = []
for actor in actors:
lim.append(actors[actor].GetBounds())
lim = np.array(lim)
self.axes[viewport].SetBounds(tuple(lim.max(axis=0)))
self.axes[viewport].SetCamera(self.cams[viewport])
self.axes[viewport].SetFlyModeToClosestTriad()
for axes in xrange(3):
self.axes[viewport].GetTitleTextProperty(axes).SetColor(color)
self.axes[viewport].GetLabelTextProperty(axes).SetColor(color)
self.axes[viewport].GetTitleTextProperty(axes).SetFontFamilyToCourier()
self.axes[viewport].GetLabelTextProperty(axes).SetFontFamilyToCourier()
self.axes[viewport].GetXAxesLinesProperty().SetColor(color)
self.axes[viewport].GetYAxesLinesProperty().SetColor(color)
self.axes[viewport].GetZAxesLinesProperty().SetColor(color)
self.axes[viewport].SetGridLineLocation(self.axes[viewport].VTK_GRID_LINES_FURTHEST)
self.axes[viewport].XAxisMinorTickVisibilityOff()
self.axes[viewport].YAxisMinorTickVisibilityOff()
self.axes[viewport].ZAxisMinorTickVisibilityOff()
self.rens[viewport].AddActor(self.axes[viewport])
class qtVtkWindow(QVTKRenderWindowInteractor, ampVTK):
def __init__(self):
super(qtVtkWindow, self).__init__()
self.SetInteractorStyle(self.style)
self.GetRenderWindow().AddRenderer(self.rens[0])
self.iren = self.GetRenderWindow().GetInteractor()
self.iren.Initialize()
def setnumViewports(self, n):
"""
Function to set multiple viewports within the vtkWindow
Parameters
------------
n: int
number of viewports required
"""
dif = n - len(self.rens)
if dif == 0:
return
elif dif < 0:
for ren in self.rens[n:]:
self.GetRenderWindow().RemoveRenderer(ren)
self.rens = self.rens[:n]
elif dif > 0:
for i in xrange(dif):
self.rens.append(vtkRender())
self.axes.append(vtk.vtkCubeAxesActor())
self.GetRenderWindow().AddRenderer(self.rens[-1])
if len(self.cams) < len(self.rens):
self.cams.append(vtk.vtkCamera())
self.rens[-1].SetActiveCamera(self.cams[len(self.rens)-1])
for i, ren in enumerate(self.rens):
ren.SetViewport(float(i)/n, 0, float(i+1)/n, 1)
self.setBackground()
class vtkRenWin(vtk.vtkRenderWindow, ampVTK):
def __init__(self, winWidth=512, winHeight=512):
super(vtkRenWin, self).__init__()
self.winWidth = winWidth
self.winHeight = winHeight
self.SetSize(self.winWidth, self.winHeight)
self.OffScreenRenderingOn()
self.AddRenderer(self.rens[0])
self.Render()
def setnumViewports(self, n):
"""
Function to set multiple viewports within the vtkWindow
Parameters
------------
n: int
number of viewports required
"""
dif = n - len(self.rens)
if dif == 0:
return
elif dif < 0:
for ren in self.rens[n:]:
self.RemoveRenderer(ren)
self.rens = self.rens[:n]
elif dif > 0:
for i in xrange(dif):
self.rens.append(vtkRender())
self.axes.append(vtk.vtkCubeAxesActor())
self.AddRenderer(self.rens[-1])
if len(self.cams) < len(self.rens):
self.cams.append(vtk.vtkCamera())
self.rens[-1].SetActiveCamera(self.cams[len(self.rens)-1])
for i, ren in enumerate(self.rens):
ren.SetViewport(float(i)/n, 0, float(i+1)/n, 1)
self.setBackground()
def getImage(self):
self.vtkRGB = vtk.vtkUnsignedCharArray()
self.GetPixelData(0, 0, self.winWidth-1, self.winHeight-1,
1, self.vtkRGB)
self.vtkRGB.Squeeze()
self.im = np.flipud(np.resize(np.array(self.vtkRGB),
[self.winWidth, self.winHeight, 3])) / 255.0
def getScreenshot(self, fname, mag=10):
w2if = vtk.vtkWindowToImageFilter()
w2if.SetInput(self)
w2if.SetMagnification(mag)
w2if.Update()
writer = vtk.vtkTIFFWriter()
writer.SetFileName(fname)
writer.SetInputConnection(w2if.GetOutputPort())
writer.Write()
class visMixin(object):
def genIm(self, actor=['limb'], winWidth=512, winHeight=512,
views=[[0, -1, 0]], background=[1.0, 1.0, 1.0], projection=True,
shading=True, mag=10):
"""
"""
win = vtkRenWin(winWidth, winHeight)
win.setnumViewports(len(views))
win.setBackground(background)
win.setProjection(projection)
for i, view in enumerate(views):
win.addAxes(self.actors, color=[0.0, 0.0, 0.0], viewport=i)
win.setView(view, i)
win.setProjection(projection, viewport=i)
win.renderActors(self.actors, actor, viewport=i, shading=shading, zoom=1.3)
win.Render()
win.getScreenshot('test.tiff')
# return win.im
def addActor(self, stype=0, CMap=None):
"""
Function to insert a vtk actor into the actors dictionary within
the AmpObject
"""
if isinstance(stype, int):
stype = self.stype[stype]
data = getattr(self, stype)
self.actors[stype] = self.ampActor(data, CMap=CMap)
class ampActor(vtk.vtkActor):
"""
Class that inherits methods from vtk actor
Contains functions to set vertices, faces, scalars and color map
from numpy arrays
Add functions to add vert, add faces, cmap and make LUT
"""
def __init__(self, data, CMap=None, bands=None):
self.mesh = vtk.vtkPolyData()
self.points = vtk.vtkPoints()
self.polys = vtk.vtkCellArray()
self.setVert(data['vert'])
self.setFaces(data['faces'])
if CMap is not None:
self.setRect(data['values'])
self.setCMap(CMap)
self.Mapper = vtk.vtkPolyDataMapper()
self.Mapper.SetInputData(self.mesh)
if CMap is not None:
self.setScalarRange()
self.Mapper.SetLookupTable(self.lut)
self.SetMapper(self.Mapper)
def setVert(self, vert):
self.points.SetData(numpy_support.numpy_to_vtk(vert, deep=1))
self.mesh.SetPoints(self.points)
def setFaces(self, faces):
f = np.c_[np.tile(faces.shape[1], faces.shape[0]), faces].flatten()
self.polys.SetCells(len(faces),
numpy_support.numpy_to_vtkIdTypeArray(f, deep=1))
self.mesh.SetPolys(self.polys)
def setRect(self, rect):
self.scalars = numpy_support.numpy_to_vtk(rect, deep=1)
self.mesh.GetPointData().SetScalars(self.scalars)
def setOpacity(self, opacity=1.0):
self.GetProperty().SetOpacity(opacity)
def setColor(self, color=[1.0, 1.0, 1.0]):
self.GetProperty().SetColor(color)
def setScalarRange(self, smin=-8.0, smax=8.0):
self.Mapper.SetScalarRange(smin, smax)
def setCMap(self, CMap, bands=128):
self.ctf = vtk.vtkColorTransferFunction()
self.ctf.SetColorSpaceToDiverging()
for ind, point in zip(np.linspace(0, 1, len(CMap)), CMap):
self.ctf.AddRGBPoint(ind, point[0], point[1], point[2])
self.lut = vtk.vtkLookupTable()
self.lut.SetNumberOfTableValues(bands)
self.lut.Build()
for i in xrange(bands):
rgb = list(self.ctf.GetColor(float(i) / bands)) + [1]
self.lut.SetTableValue(i, rgb)
def setShading(self, shading=True):
if shading is True:
self.GetProperty().LightingOn()
if shading is False:
self.GetProperty().LightingOff()
\ No newline at end of file
# -*- coding: utf-8 -*-
"""
Created on Tue Jan 10 12:49:32 2017
@author: js22g12
"""
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from collections import defaultdict
class analyseMixin(object):
def plot_slices(self, axis='Z', SliceWidth=3, stype=0):
if isinstance(stype, int):
stype = self.stype[stype]
data = getattr(self, stype)
# Find the brim edges
ind = np.where(np.isnan(data['faceEdges'][:,1]))[0]
# Define max Z from lowest point on brim
maxZ = data['vert'][data['edges'][ind, :], 2].min()
fig = plt.figure()
fig.set_size_inches(12, 9)
ax1 = fig.add_subplot(221, projection='3d')
ax2 = fig.add_subplot(222)
#Z position of slices
slices = np.arange(data['vert'][:,2].min() + SliceWidth, maxZ, SliceWidth)
PolyArea = np.zeros([len(slices)])
for i, pos in enumerate(slices):
# Create a slice for each z position
poly = analyseMixin.create_slice(data, pos, axis)
ax1.plot(poly[:,0],
poly[:,1],
poly[:,2],
c='b')
#SlicePolys[i, :] = poly
# Compute area of slice
area = 0.5*np.abs(np.dot(poly[:,0], np.roll(poly[:,1], 1)) -
np.dot(poly[:,1], np.roll(poly[:,0], 1)))
PolyArea[i] = area
extents = np.array([getattr(ax1, 'get_{}lim'.format(dim))() for dim in 'xyz'])
sz = extents[:,1] - extents[:,0]
centers = np.mean(extents, axis=1)
maxsize = max(abs(sz))
r = maxsize/2
for ctr, dim in zip(centers, 'xyz'):
getattr(ax1, 'set_{}lim'.format(dim))(ctr - r, ctr + r)
ax1.set_axis_off()
ax2.plot(slices-slices[0], PolyArea)
# Rendering of the limb scan
ax3 = fig.add_subplot(2,2,3)
self.addActor(stype='limb')
Im = self.genIm(actor=['limb'])
ax3.imshow(Im, None)
ax3.set_axis_off()
# Rendering of the rectification map
ax4 = fig.add_subplot(2,2,4)
CMap = np.array([[212.0, 221.0, 225.0],
[31.0, 73.0, 125.0]])/255.0
self.addActor(stype='reglimb', CMap = CMap)
Im = self.genIm(actor=['reglimb'])
ax4.imshow(Im, None)
ax4.set_axis_off()
plt.tight_layout()
@staticmethod
def create_slice(data, plane, axis='Z'):
# Find all vertices below plane
ind = data['vert'][:,2][data['edges']] < plane
# Select edges with one vertex above and one below the slice plane
validEdgeInd = np.where(np.logical_xor(ind[:,0], ind[:,1]))[0]
validfE = data['faceEdges'][validEdgeInd, :].astype(int)
g = defaultdict(set)
faceOrder = np.zeros(len(validEdgeInd), dtype=int)
# Run eularian path algorithm to order faces
for v, w in validfE:
g[v].add(w)
g[w].add(v)
v = validfE[0,0]
i=0
while True:
try:
w = g[v].pop()
except KeyError:
break
g[w].remove(v)
faceOrder[i] = v
i+=1
v = w
# Get array of three edges attached to each face
validEdges = data['edgesFace'][faceOrder, :]
# Remove the edge that is not intersected by the plane
edges = validEdges[np.isin(validEdges, validEdgeInd)].reshape([-1,2])
# Remove the duplicate edge from order
e = edges.flatten()
odx = np.argsort(e)
inds = np.arange(1, len(e), 2)
row = np.unravel_index(odx, e.shape)[0]
mask = np.ones(len(e), dtype=bool)
mask[row[inds]] = False
sortE = e[mask]
# Add first edge to end of array
sortE = np.append(sortE, sortE[0])
polyEdge = data['edges'][sortE]
EdgePoints = np.c_[data['vert'][polyEdge[:,0], :],
data['vert'][polyEdge[:,1], :]]
# Create poly from
poly = analyseMixin.planeEdgeintersect(EdgePoints, plane, axis=axis)
return poly
def create_slice_fast(self, plane, axis='Z'):
# Find all vertices below plane
ind = self.vert[:,2][self.edges] < plane
# Select edges with one vertex above and one below the slice plane
validEdgeInd = np.where(np.logical_xor(ind[:,0], ind[:,1]))[0]
validfE = self.faceEdges[validEdgeInd, :].astype(int)
col0, col0ind = np.unique(validfE[:,0], return_index=True)
f0 = np.c_[col0, validfE[col0ind, 1]]
col0indrep = np.where(~np.isin(np.arange(len(validfE)), col0ind))[0]
f1 = validfE[col0indrep, :]
col1, col1ind = np.unique(validfE[:,1], return_index=True)
f1 = np.r_[np.c_[col1, validfE[col1ind, 0]], f1]
col1indrep = np.where(~np.isin(np.arange(len(validfE)), col1ind))[0]
f0 = np.r_[f0, validfE[col1indrep, :][:, [1,0]]]
f0_idx = np.argsort(f0[:,0])
f1_idx = np.argsort(f1[:,0])
h = dict(zip(f0[f0_idx, 0].tolist(),
np.c_[f0[f0_idx, 1], f1[f1_idx, 1]].tolist()))
faceOrder = np.zeros(len(validEdgeInd), dtype=int)
# Run eularian path algorithm to order faces
v = validfE[0,0]
for i in xrange(len(validEdgeInd)):
w = h[v].pop()
h[w].remove(v)
faceOrder[i] = v
v = w
# Get array of three edges attached to each face
validEdges = self.edgesFace[faceOrder, :]
# Remove the edge that is not intersected by the plane
edges = validEdges[np.isin(validEdges, validEdgeInd)].reshape([-1,2])
# Remove the duplicate edge from order
e = edges.flatten()
odx = np.argsort(e)
inds = np.arange(1, len(e), 2)
row = np.unravel_index(odx, e.shape)[0]
mask = np.ones(len(e), dtype=bool)
mask[row[inds]] = False
sortE = e[mask]
# Add first edge to end of array
sortE = np.append(sortE, sortE[0])
polyEdge = self.edges[sortE]
EdgePoints = np.c_[self.vert[polyEdge[:,0], :],
self.vert[polyEdge[:,1], :]]
# Create poly from
poly = analyseMixin.planeEdgeintersect(EdgePoints, plane, axis=axis)
return poly
@staticmethod
def planeEdgeintersect(edges, plane, axis='Z'):
axisInd = 0 if axis == 'X' else 1 if axis == 'Y' else 2
intersectPoints = np.zeros((edges.shape[0], 3))
intersectPoints[:, axisInd] = plane
axesInd = np.arange(0, 2, 1)[np.arange(0, 2, 1) != axisInd]
for i in axesInd:
intersectPoints[:, i] = (edges[:, i] +
(plane - edges[:, axisInd]) *
(edges[:, i+3] - edges[:, i]) /
(edges[:, axisInd+3] - edges[:, axisInd]))
return intersectPoints
# -*- coding: utf-8 -*-
"""
Created on Thu Sep 14 13:15:30 2017
@author: js22g12
"""
import numpy as np
from scipy import spatial
from scipy.optimize import minimize
class alignMixin(object):
def icp(self):
"""
Autmated alignment function between two
"""
tTree = spatial.cKDTree(self.baseline.vert)
rot = np.array([0,0,0], dtype=float)
res = minimize(self.calcDistError, rot, method='BFGS',
options={'gtol':1e-6, 'disp':True})
def calcDistError(self, rot, tTree):
Id = np.identity(3)
for i in xrange(3):
if rot[i] != 0:
ax = Id[i, :]
ang = np.deg2rad(rot[i])
dot = np.reshape(self.vert[:, 0] * ax[0] +
self.vert[:, 1] * ax[1] +
self.vert[:, 2] * ax[2], (-1, 1))
self.vert = (self.vert * np.cos(ang) +
np.cross(ax, self.vert) * np.sin(ang) +
np.reshape(ax, (1, -1)) * dot * (1-np.cos(ang)))
dist = tTree.query(self.vert, 10)[0]
dist = dist.min(axis=1)
return dist.sum()
\ No newline at end of file
# -*- coding: utf-8 -*-
"""
Created on Wed Sep 13 13:54:23 2017
@author: js22g12
Core functions for the AmpObject
Remove pd dependency and instead just use numpy arrays
Requires numpy 1.13
import os
path = (r'J:\\Shared Resources\\AmpScan IfLS Team\\'
'100 PYTHON\\STLReader')
path = r'\\filestore.soton.ac.uk\SharedResources\AmpScan IfLS Team\100 PYTHON\STLReader'
path = (r'C:\Users\js22g12\OneDrive - University of Southampton\Documents '
r'(OneDrive)\AmpScan\Code\2017_09\02_Code\AmpScan')
path = (r'C:\Users\Josh\OneDrive - University of Southampton\Documents '
r'(OneDrive)\AmpScan\Code\2017_09\02_Code\AmpScan')
filename = '01_PhantomShell_ICEM_3mm.stl'
filename2 = '01_PhantomShell_ICEM_3mm_write.stl'
os.chdir(path)
AmpObject
Read
Write
Centre
Translate
Rotate
lp_smooth
trimming
alignment
Residuum
Socket
SocketICP (overwrite)
SocketBrimLine
Registration
SoftTissueDepth
Bones
Liner
Finite Element Analysis
FE Mesh
FE Data
"""
import numpy as np
import struct
from autoAlign import alignMixin
from trim import trimMixin
from smooth import smoothMixin
from analyse import analyseMixin
from ampVis import visMixin
from fe import feMixin
from tsbSocketDesign import socketDesignMixin
class AmpObject(alignMixin, trimMixin, smoothMixin, analyseMixin,
visMixin, feMixin, socketDesignMixin):
def __init__(self, Data, stype):
self.stype = []
self.actors = {}
if stype in ['limb', 'socket', 'reglimb', 'regsocket', 'MRI']:
self.addData(Data, stype)
elif stype is 'AmpObj':
for d in Data.stype:
setattr(self, d, getattr(Data, d))
self.stype.append(d)
self.actors = Data.actors
elif stype is 'FE':
self.addFE([Data,])
else:
raise ValueError('dtype not supported, please choose from ' +
'limb, socket, reglimb, regsocket, MRI or AmpObj')
def addData(self, Data, stype):
if isinstance(Data, basestring):
self.stype.append(stype)
self.read_stl(Data, stype)
# Import stl as filename
elif isinstance(Data, dict):
self.stype.append(stype)
setattr(self, stype, Data)
def read_stl(self, filename, stype=0, unify=True, edges=True):
"""
Function to read .stl file from filename and import data into
the AmpObj
Parameters
-----------
filename: str
file path of the .stl file to read
unify: boolean, default True
unify the coincident vertices of each face
edges: boolean, default True
calculate the edges array automatically
"""
if isinstance(stype, int):
stype = self.stype[stype]
fh = open(filename, 'rb')
# Defined no of bytes for header and no of faces
HEADER_SIZE = 80
COUNT_SIZE = 4
# State the data type and length in bytes of the normals and vertices
data_type = np.dtype([('normals', np.float32, (3, )),
('vertices', np.float32, (9, )),
('atttr', '<i2', (1, ))])
# Read the header of the STL
fh.read(HEADER_SIZE).lower()
# Read the number of faces
NFaces, = struct.unpack('@i', fh.read(COUNT_SIZE))
# Read the remaining data and save as void, then close file
data = np.fromfile(fh, data_type)
fh.close()
# Write the data to a numpy arrays in AmpObj
vert = np.resize(np.array(data['vertices']), (NFaces*3, 3))
norm = np.array(data['normals'])
faces = np.reshape(range(NFaces*3), [-1,3])
setattr(self, stype, dict(zip(['vert', 'faces', 'norm'],
[vert, faces, norm])))
# Call function to unify vertices of the array
if unify is True:
self.unify_vertices(stype)
# Call function to calculate the edges array
if edges is True:
self.computeEdges(stype)
def unify_vertices(self, stype=0):
"""
Function to unify coincident vertices of the mesh to reduce
size of the vertices array enabling speed increases
"""
# Requires numpy 1.13
if isinstance(stype, int):
stype = self.stype[stype]
data = getattr(self, stype)
data['vert'], indC = np.unique(data['vert'], return_inverse=True, axis=0)
# Maps the new vertices index to the face array
data['faces'] = np.resize(indC[data['faces']], (len(data['norm']), 3))
def computeEdges(self, stype=0):
"""
Function to compute the edges array, the edges on each face,
and the faces on each edge
edges: numpy array N x 2 denotes the indicies of two vertices
on each edge
edgesFace: numpy array N x 3 denotes the indicies of the three edges
on each face
faceEdges: numpy array N x 2 denotes the indicies of the faces in each
edge, edges may have either 1 or 2 faces, if 1 then the second
index will be NaN
"""
if isinstance(stype, int):
stype = self.stype[stype]
data = getattr(self, stype)
# Get edges array
data['edges'] = np.reshape(data['faces'][:, [0, 1, 0, 2, 1, 2]], [-1, 2])
data['edges'] = np.sort(data['edges'], 1)
# Get edges on each face
data['edgesFace'] = np.reshape(range(len(data['faces'])*3), [-1,3])
# Unify the edges
data['edges'], indC = np.unique(data['edges'], return_inverse=True, axis=0)
#Remap the edgesFace array
data['edgesFace'] = indC[data['edgesFace'] ]
#Initiate the faceEdges array
data['faceEdges'] = np.empty([len(data['edges']), 2])
data['faceEdges'].fill(np.nan)
# Denote the face index for flattened edge array
fInd = np.repeat(np.array(range(len(data['faces']))), 3)
# Flatten edge array
eF = np.reshape(data['edgesFace'], [-1])
eFInd = np.unique(eF, return_index=True)[1]
logic = np.zeros([len(eF)], dtype=bool)
logic[eFInd] = True
data['faceEdges'][eF[logic], 0] = fInd[logic]
data['faceEdges'][eF[~logic], 1] = fInd[~logic]
def save(self, filename, stype=0):
"""
Function to save the AmpObj as a binary .stl file
Parameters
-----------
filename: str
file path of the .stl file to save to
"""
if isinstance(stype, int):
stype = self.stype[stype]
self.calc_norm(stype)
data = getattr(self, stype)
fv = data['vert'][np.reshape(data['faces'], len(data['faces'])*3)]
fh = open(filename, 'wb')
data_type = np.dtype([('normals', np.float32, (3, )),
('vertices', np.float32, (9, )),
('atttr', '<i2', (1, ))])
header = '%s' % (filename)
header = header[:80].ljust(80, ' ')
packed = struct.pack('@i', len(data['faces']))
fh.write(header)
fh.write(packed)
data_type = np.dtype([('normals', np.float32, (3, )),
('vertices', np.float32, (9, )),
('atttr', '<i2', (1, ))])
data_write = np.zeros(len(data['faces']), dtype=data_type)
data_write['normals'] = data['norm']
data_write['vertices'] = np.reshape(fv, (len(data['faces']), 9))
data_write.tofile(fh)
fh.close()
def calc_norm(self, stype=0):
"""
Calculate the normal of each face of the AmpObj
"""
if isinstance(stype, int):
stype = self.stype[stype]
data = getattr(self, stype)
norms = np.cross(data['vert'][data['faces'][:,1]] -
data['vert'][data['faces'][:,0]],
data['vert'][data['faces'][:,2]] -
data['vert'][data['faces'][:,0]])
mag = np.linalg.norm(norms, axis=1)
data['norm'] = np.divide(norms, mag[:,None])
def man_trans(self, trans, stype=0):
"""
Translate the AmpObj in 3D space
Parameters
-----------
trans: array-like
1x3 array of the tranlation in [x, y, z]
"""
if isinstance(stype, int):
stype = self.stype[stype]
data = getattr(self, stype)
data['vert'] += trans
def centre(self, stype=0):
"""
Centre the AmpObj based upon the mean of all the vertices
"""
if isinstance(stype, int):
stype = self.stype[stype]
data = getattr(self, stype)
self.man_trans(-data['vert'].mean(axis=0))
def man_rot(self, rot, stype=0):
"""
Rotate the AmpObj in 3D space and re-calculate the normals
Parameters
-----------
rot: array-like
1x3 array of the rotation around [x, y, z]
"""
if isinstance(stype, int):
stype = self.stype[stype]
data = getattr(self, stype)
Id = np.identity(3)
for i, r in enumerate(rot):
if r != 0:
ax = Id[i, :]
ang = np.deg2rad(rot[i])
dot = np.reshape(data['vert'][:, 0] * ax[0] +
data['vert'][:, 1] * ax[1] +
data['vert'][:, 2] * ax[2], (-1, 1))
data['vert'] = (data['vert'] * np.cos(ang) +
np.cross(ax, data['vert']) * np.sin(ang) +
np.reshape(ax, (1, -1)) * dot * (1-np.cos(ang)))
self.calc_norm()
# -*- coding: utf-8 -*-
"""
Created on Wed Nov 01 14:13:50 2017
@author: js22g12
"""
import numpy as np
class feMixin(object):
def addFE(self, files):
if len(files) == 1:
data = np.load(files[0]).item()
if len(files) == 3:
data = {}
names = ['vert', 'faces', 'values']
for n, f in zip(names, files):
data[n] = np.loadtxt(f)
self.FE = data
self.getSurf()
def getSurf(self):
valInd = self.FE['values'][:, 0].astype(int)
log = np.isin(self.FE['faces'], valInd)
f = self.FE['faces'][log].reshape([-1, 4])
log = np.zeros(len(self.FE['vert']), dtype=bool)
log[valInd] = True
fInd = np.cumsum(log) - 1
self.FE['vert'] = self.FE['vert'][log, :]
self.FE['faces'] = fInd[f].astype(np.int64)
self.FE['values'] = np.array(self.FE['values'][:, 1])
self.FE['edges'] = np.reshape(self.FE['faces'][:, [0, 1, 0, 2, 0, 3, 1, 2, 1, 3, 2, 3]], [-1, 2])
self.FE['edges'] = np.sort(self.FE['edges'], 1)
# Unify the edges
self.FE['edges'], indC = np.unique(self.FE['edges'], return_inverse=True, axis=0)
\ No newline at end of file
# -*- coding: utf-8 -*-
"""
Created on Wed Sep 27 12:21:37 2017
@author: js22g12
"""
import numpy as np
import vtk
from vtk.util import numpy_support
from vtk.qt.QVTKRenderWindowInteractor import QVTKRenderWindowInteractor
class pressSense(object):
def __init__(self):
self.actors = {}
self.antS = {}
def importVert(self, fname, pos='ant'):
self.antS['vert'] = np.loadtxt(fname, delimiter=',')
self.antS['vert'][:, 1] += -1.0
def importPress(self, fname):
sF = self.antS['sF']
val = np.loadtxt(fname, delimiter=',').flatten()
self.antS['values'] = np.zeros(len(self.antS['faces']))
self.antS['values'][sF] = val.flatten()[:, None]
def calcVert(self, mesh, cLine, sensePos='ant', limbSide='L', d=3):
"""
Function to compute vertex position of pressure sensors
Parameters
----------
vert: numpy array
vertices for mesh to project the pressure sensor onto
CLine: numpy array
vertices of the begining and end of the sensor centre line on the
socket
sensePos: str
Aspect of socket the sensor was positioned
limbSide:
"""
vert = mesh['vert']
faces = mesh['faces']
if sensePos == 'ant':
plane = np.array([0, -1, 0])
elif sensePos == 'lat' and limbSide == 'L':
plane = np.array([1, 0, 0])
elif sensePos == 'lat' and limbSide == 'R':
plane = np.array([-1, 0, 0])
rows, cols = [15*d + 1, 3*d + 1]
planeInd = [i for i, p in enumerate(plane) if p == 0]
lim = cLine[:, planeInd]
h = np.linspace(0, lim[1, 0]- lim[0, 0], 100)
v = np.linspace(0, lim[1, 1]- lim[0, 1], 100)
# Define co-ordinate transformation to translate to bottom
# of sensor and rotate around y axis
trans = np.zeros([4,4])
trans[planeInd, 3] = lim[0, :]
xTheta = -np.arctan((lim[0, 0] - lim[1, 0])/(lim[0, 1] - lim[1, 1]))
trans = np.identity(4) - trans
# rotation matrix
R = np.array([[np.cos(xTheta), 0, np.sin(xTheta),0],
[0, 1, 0, 0],
[-np.sin(xTheta), 0, np.cos(xTheta), 0],
[0, 0, 0, 1]])
M = np.dot(R, trans)
# Apply transformation
vertT = np.dot(M, np.transpose(np.c_[vert, np.ones(len(vert))]))
vertT = np.transpose(vertT[:-1, :])
# get faceCentroids )
# Equation of line
# Find nearest 10 centroids to line
# Test if line intersects faces using barycentric technique
# Intersection between line and face plane
def calcFaces(self, d=5):
"""
Function to compute face array for pressure sensor
FIX SO USE QUAD FACES
"""
# define no of rows and columns in sensors
rows, cols = [15*d, 3*d]
# Create face array
f = np.zeros([rows*cols*2, 3], dtype=np.int64)
# Create array to mapy from sensor data to face colour
sF = np.zeros([45, 2*d**2], dtype=np.int64)
# Define face
f0 = np.array([[0, cols+1, 1],
[1, cols+1, cols+2]])
for row in xrange(rows):
for col in xrange(cols):
ind = row*cols*2 + col*2
ind2 = row*(cols+1) + col
f[[ind, ind+1], :] = f0 + ind2
ind = np.arange(0, cols*2*d, cols*2)[:, None]
for row in xrange(15):
inds = np.array([np.arange(2*d) + d*row*cols*2]*d)
inds += np.arange(0, cols*2*d, cols*2)[:, None]
for col in xrange(3):
sF[row*3 + col, :] = inds.flatten()
inds += 2*d
self.antS['faces'] = f
self.antS['sF'] = sF
def calcFacesHex(self, d=5):
"""
Function to compute face array for pressure sensor
FIX SO USE QUAD FACES
"""
# define no of rows and columns in sensors
rows, cols = [15*d, 3*d]
# Create face array
f = np.zeros([rows*cols*2, 3], dtype=np.int64)
# Create array to mapy from sensor data to face colour
sF = np.zeros([45, 2*d**2], dtype=np.int64)
f0 = np.array([[0, cols+1, 1],
[1, cols+1, 2]])
for row in xrange(rows):
for col in xrange(cols):
ind = row*cols*2 + col*2
ind2 = row*(cols+1) + col
f[[ind, ind+1], :] = f0 + ind2
ind = np.arange(0, cols*2*d, cols*2)[:, None]
for row in xrange(15):
inds = np.array([np.arange(2*d) + d*row*cols*2]*d)
inds += np.arange(0, cols*2*d, cols*2)[:, None]
for col in xrange(3):
sF[row*3 + col, :] = inds.flatten()
inds += 2*d
self.antS['faces'] = f
self.antS['sF'] = sF
def intersectLineMesh(line, v, f):
"""
Function to calculate intersection between line and mesh
"""
def addActor(self, pos ='antS', CMap=None, connect=3):
"""
Function to insert a vtk actor into the actors dictionary within
the AmpObject
"""
self.actors[pos] = self.pressActor(self.antS, CMap=CMap, connect=connect)
class pressActor(vtk.vtkActor):
"""
Class that inherits methods from vtk actor
Contains functions to set vertices, faces, scalars and color map
from numpy arrays
Add functions to add vert, add faces, cmap and make LUT
"""
def __init__(self, data, CMap=None, bands=None, connect=3):
self.mesh = vtk.vtkPolyData()
self.points = vtk.vtkPoints()
self.polys = vtk.vtkCellArray()
self.setVert(data['vert'])
self.setFaces(data['faces'], connect)
if CMap is not None:
self.setPress(data['values'])
self.setCMap(CMap)
self.Mapper = vtk.vtkPolyDataMapper()
self.Mapper.SetInputData(self.mesh)
if CMap is not None:
self.setScalarRange()
self.Mapper.SetLookupTable(self.lut)
self.SetMapper(self.Mapper)
def setVert(self, vert):
self.points.SetData(numpy_support.numpy_to_vtk(vert, deep=1))
self.mesh.SetPoints(self.points)
def setFaces(self, faces, connect=3):
f = np.c_[np.tile(connect, len(faces)), faces].flatten()
self.polys.SetCells(len(faces),
numpy_support.numpy_to_vtkIdTypeArray(f, deep=1))
self.mesh.SetPolys(self.polys)
def setPress(self, press):
self.scalars = numpy_support.numpy_to_vtk(press, deep=1)
self.mesh.GetCellData().SetScalars(self.scalars)
def setScalarRange(self, smin=0.0, smax=100.0):
self.Mapper.SetScalarRange(smin, smax)
def setCMap(self, CMap, bands=128):
self.ctf = vtk.vtkColorTransferFunction()
self.ctf.SetColorSpaceToDiverging()
for ind, point in zip(np.linspace(0, 1, len(CMap)), CMap):
self.ctf.AddRGBPoint(ind, point[0], point[1], point[2])
self.lut = vtk.vtkLookupTable()
self.lut.SetNumberOfTableValues(bands)
self.lut.Build()
for i in xrange(bands):
rgb = list(self.ctf.GetColor(float(i) / bands)) + [1]
self.lut.SetTableValue(i, rgb)
def func1():
path = r'J:\Shared Resources\AmpScan IfLS Team\101 ImpAmp\Upgrade'
sockets = ['B_PTB.stl', 'B_TSB.stl', 'B_KBM.stl']
pressure = ['161108_Sub02_T001', '161108_Sub02_T002', '161108_Sub02_T006']
files = ['B_PTB_LOC', 'B_TSB_LOC', 'B_KBM_LOC']
LimbSide = 'L'
d = 5
i = 1
socket = path + '\\' + sockets[0]
Data = AmpObject(socket, 'limb')
vert = Data.limb['vert']
posData = np.genfromtxt(path + '\\' + files[0] + '.csv',
delimiter=',', skip_header=1)[:, 1::]
antData = np.loadtxt(path + '\\' + pressure[0] + '_ant.csv', delimiter = ',')
latData = np.loadtxt(path + '\\' + pressure[0] + '_lat.csv', delimiter = ',')
cLine = posData[[0,1], :]
def func2():
CMap = np.array([[212.0, 221.0, 225.0],
[31.0, 73.0, 125.0]])/255.0
vfname = (r'C:\Users\js22g12\OneDrive - University of Southampton'
r'\Documents (OneDrive)\AmpScan\Code\2017_09'
r'\02_Code\AmpScan\PTB_AntVert.csv')
pfname = (r'J:\Shared Resources\AmpScan IfLS Team\101 ImpAmp\Upgrade'
r'\161108_Sub02_T001_ant.csv')
pressData = pressSense()
pressData.calcFaces(d=5)
pressData.importVert(vfname)
pressData.importPress(pfname)
pressData.addActor(CMap=CMap)
# -*- coding: utf-8 -*-
"""
Created on Wed Sep 13 16:07:10 2017
@author: js22g12
"""
import numpy as np
import pandas as pd
from scipy import spatial
from core import AmpObject
"""
import os
path = ('J:\\Shared Resources\\AmpScan IfLS Team\\'
'100 PYTHON\\STLReader')
filename = '01_PhantomShell_ICEM_3mm.stl'
filename2 = '01_PhantomShell_ICEM_3mm_write.stl'
baseline = 'B_PTB.stl'
target = 'B_NormalLiner.stl'
os.chdir(path)
ampObj = AmpObject(filename)
regObj = regObject(ampObj)
ampObj.vert.loc[0,:]
regObj.vert.loc[0,:]
ampObj.icp()
regObj.icp()
regObj.registration()
ampObj.registration()
Data = AmpObject(target)
Data.getBaseline(baseline)
Reg = regObject(Data)
Reg.reg_fast()
Use mesh object as parent class with standard function
and then inherit for all other classes ie socket, residuum, registered
Either create from filename, object or fv
Standard functions:
Rotation
Translation
Read
Write
Normals
Slice analyse
Child classes:
Residuum
Socket
Registration (Target)
Soft tissue mesh (Bones, Liner)
FE mesh
"""
class regObject(AmpObject):
def __init__(self, Data=None, stype='AmpObj'):
super(regObject, self).__init__(Data, stype)
def registration(self, steps=1, baseline='limb',
target='socket', reg = 'reglimb', direct=True):
"""
Function to register the regObject to the baseline mesh
Parameters
----------
Steps: int, default 1
Number of iterations
"""
bData = getattr(self, baseline)
tData = getattr(self, target)
bV = bData['vert']
# Calculate the face centroids of the regObject
tData['fC'] = tData['vert'][tData['faces']].mean(axis=1)
# Construct knn tree
tTree = spatial.cKDTree(tData['fC'])
for step in np.arange(steps, 0, -1):
# Index of 10 centroids nearest to each baseline vertex
ind = tTree.query(bV, 10)[1]
D = np.zeros(bV.shape)
# Define normals for faces of nearest faces
norms = tData['norm'][ind]
# Get a point on each face
fPoints = tData['vert'][tData['faces'][ind, 0]]
# Calculate dot product between point on face and normals
d = np.einsum('ijk, ijk->ij', norms, fPoints)
t = d - np.einsum('ijk, ik->ij', norms, bV)
# Calculate new points
G = np.einsum('ijk, ij->ijk', norms, t)
GMag = np.sqrt(np.einsum('ijk, ijk->ij', G, G)).argmin(axis=1)
# Define vector from baseline point to intersect point
D = G[np.arange(len(G)), GMag, :]
bV = bV + D/step
regData = dict(bData)
regData['vert'] = bV
setattr(self, reg, regData)
self.calcError(baseline, reg, direct)
def calcError(self, baseline='limb', target='reglimb', direct=True):
# This is kinda slow
bData = getattr(self, baseline)
tData = getattr(self, target)
if direct is True:
values = np.linalg.norm(tData['vert'] - bData['vert'], axis=1)
# Calculate vertex normals on target from normal of surrounding faces
vNorm = np.zeros(tData['vert'].shape)
for face, norm in zip(tData['faces'], tData['norm']):
vNorm[face, :] += norm
vNorm = vNorm / np.linalg.norm(vNorm, axis=1)[:, None]
# Calculate the unit vector normal between corresponding vertices
# baseline and target
vector = (tData['vert'] - bData['vert'])/values[:, None]
# Calculate angle between the two unit vectors using normal of cross
# product between vNorm and vector and dot
normcrossP = np.linalg.norm(np.cross(vector, vNorm), axis=1)
dotP = np.einsum('ij,ij->i', vector, vNorm)
angle = np.arctan2(normcrossP, dotP)
polarity = np.ones(angle.shape)
polarity[angle < np.pi/2] =-1.0
tData['values'] = values * polarity
else:
tData['values'] = np.linalg.norm(tData['vert'] - bData['vert'],
axis=1)
# -*- coding: utf-8 -*-
"""
Created on Thu Sep 14 13:25:23 2017
@author: js22g12
"""
import numpy as np
import pandas as pd
class smoothMixin(object):
def lp_smooth(self, stype=0, n=1):
"""
Function to apply a simple laplacian smooth to the mesh
Parameters
---------
n: int, default 1
number of iterations of smoothing
"""
if isinstance(stype, int):
stype = self.stype[stype]
data = getattr(self, stype)
# Flatten the edges array to 1D
e = data['edges'].flatten()
# Get the indicies to sort edges
o_idx = e.argsort()
# Get indicies of sorted array where last of each vertex index
# occurs
ndx = np.searchsorted(e[o_idx], np.arange(len(data['vert'])),
side='right')
# Map indicies between flatted edges array and standard
row, col = np.unravel_index(o_idx, data['edges'].shape)
for i in np.arange(n):
# List all vertices
neighVerts = data['vert'][data['edges'][row, 1-col], :]
# Split into list with each array contating vertices
spl = np.split(neighVerts, ndx[0:-1])
# Get average of each array
vert = [vert.mean(axis=0) for vert in spl]
# Write to the AmpObj
data['vert'] = np.array(vert)
\ No newline at end of file
# -*- coding: utf-8 -*-
"""
Created on Tue Oct 10 16:21:12 2017
@author: js22g12
"""
import sys
import numpy as np
from core import AmpObject
from ampVis import qtVtkWindow
from tsbSocketDesign import dragSpline
from PyQt5.QtWidgets import (QAction, QApplication, QGridLayout,
QMainWindow, QFileDialog, QWidget)
from PyQt5.QtGui import QIcon
from PyQt5.QtCore import pyqtSignal
class GUI(QMainWindow):
mplChange = pyqtSignal()
def __init__(self, parent=None):
super(GUI, self).__init__()
self.AmpObj = None
self.CMap = np.array([[212.0, 221.0, 225.0],
[31.0, 73.0, 125.0]])/255.0
self.splineWidget = QWidget()
self.mainWidget = QWidget()
self.vtkWidget = qtVtkWindow()
self.vtkWidget.setBackground(color=[1,1,1])
self.setCentralWidget(self.mainWidget)
self.createActions()
self.createMenus()
self.Layout = QGridLayout()
self.splineWin = dragSpline(self.splineWidget, self.mplChange)
self.mplChange.connect(self.plotRect)
self.Layout.addWidget(self.splineWin, 0, 1, 1, 4)
self.Layout.addWidget(self.vtkWidget, 0, 0, 1, 1)
self.mainWidget.setLayout(self.Layout)
self.setWindowTitle("AmpScan Visualiser")
self.resize(1200, 800)
self.show()
def plotRect(self):
if self.AmpObj is None:
return
self.AmpObj.TSBSocket(self.splineWin.B, stype='reglimb')
self.AmpObj.actors['reglimb'].setRect(self.AmpObj.reglimb['values'])
self.vtkWidget.iren.Render()
def chooseOpenFile(self):
self.fname = QFileDialog.getOpenFileName(self, 'Open file',
filter="Meshes (*.stl)")
self.AmpObj = AmpObject(self.fname[0], 'reglimb')
self.AmpObj.centre()
self.AmpObj.TSBSocket(self.splineWin.B, stype='reglimb')
self.AmpObj.addActor(stype='reglimb', CMap=self.CMap)
self.AmpObj.actors['reglimb'].setScalarRange(0, 6)
self.vtkWidget.renderActors(self.AmpObj.actors, dispActors=['reglimb', ])
self.vtkWidget.setScalarBar(self.AmpObj.actors['reglimb'])
def createActions(self):
self.openFile = QAction(QIcon('open.png'), 'Open', self,
shortcut='Ctrl+O',
triggered=self.chooseOpenFile)
def createMenus(self):
self.fileMenu = self.menuBar().addMenu("&File")
self.fileMenu.addAction(self.openFile)
def socketDesignGUI():
# if __name__ == "__main__":
app = QApplication(sys.argv)
mainWin = GUI()
mainWin.show()
sys.exit(app.exec_())
\ No newline at end of file
# -*- coding: utf-8 -*-
"""
Created on Thu Sep 14 13:22:29 2017
@author: js22g12
"""
import numpy as np
import pandas as pd
class trimMixin(object):
def planarTrim(self):
print('This is where the planar trim function will go')
\ No newline at end of file
# -*- coding: utf-8 -*-
"""
Created on Wed Oct 11 09:05:19 2017
@author: js22g12
Backend for matplotlib to create a moveable spline figure
"""
from matplotlib.backends.backend_qt5agg import FigureCanvasQTAgg as FigureCanvas
from matplotlib.figure import Figure
from PyQt5.QtWidgets import QSizePolicy
import numpy as np
from scipy.special import binom
class socketDesignMixin(object):
def TSBSocket(self, B, stype=0):
if isinstance(stype, int):
stype = self.stype[stype]
data = getattr(self, stype)
rho = np.sqrt(data['vert'][:, 0]**2 + data['vert'][:,1]**2)
zRange = data['vert'][:, 2].max() - data['vert'][:, 2].min()
zB = (B[:, 0] * zRange) + data['vert'][:, 2].min()
perRed = np.interp(data['vert'][:, 2], zB, B[:, 1])
data['values'] = rho * (perRed * 0.01)
class mplCanvas(FigureCanvas):
"""Ultimately, this is a QWidget (as well as a FigureCanvasAgg, etc.)."""
def __init__(self, parent=None, width=4, height=5, dpi=100):
fig = Figure(figsize=(width, height), dpi=dpi)
self.axes = fig.add_subplot(111)
# fig.patch.set_facecolor((0.1, 0.2, 0.4))
# self.axes.set_facecolor((0.1, 0.2, 0.4))
self.compute_initial_figure()
super(mplCanvas, self).__init__(fig)
# FigureCanvas.__init__(self, fig)
self.setParent(parent)
FigureCanvas.setSizePolicy(self,
QSizePolicy.Expanding,
QSizePolicy.Expanding)
FigureCanvas.updateGeometry(self)
def compute_initial_figure(self):
pass
class dragSpline(mplCanvas):
lock = None #only one can be animated at a time
def __init__(self, parent, Signal=None):
super(dragSpline, self).__init__(parent)
self.Signal = Signal
self.points = np.array([[0.0, 8.0],
[0.5, 4.0],
[1.0, 0.0]])
self.weights = np.array([1.0, 5.0, 1.0])
self.axes.set_ylim([0, 1])
self.axes.set_xlim([0, 8])
self.axes.set_ylabel('Normalised distance along socket')
self.axes.set_xlabel('Percentage reduction in volume')
self.point = self.axes.plot(self.points[:, 1],
self.points[:, 0],
'ob', markersize=10)[0]
self.spline = self.axes.plot(0, 0)[0]
self.bezierCurve()
self.press = None
self.background = None
self.connect()
def connect(self):
'connect to all the events we need'
self.cidpress = self.point.figure.canvas.mpl_connect('button_press_event', self.on_press)
self.cidrelease = self.point.figure.canvas.mpl_connect('button_release_event', self.on_release)
self.cidmotion = self.point.figure.canvas.mpl_connect('motion_notify_event', self.on_motion)
def on_press(self, event):
if event.inaxes != self.point.axes: return
if dragSpline.lock is not None: return
contains, attrd = self.point.contains(event)
if not contains: return
self.press = self.point.get_xydata(), attrd['ind'][0], event.xdata, event.ydata
dragSpline.lock = self
# draw everything but the selected rectangle and store the pixel buffer
canvas = self.point.figure.canvas
axes = self.point.axes
self.point.set_animated(True)
self.spline.set_animated(True)
canvas.draw()
self.background = canvas.copy_from_bbox(self.point.axes.bbox)
# now redraw just the rectangle
axes.draw_artist(self.point)
axes.draw_artist(self.spline)
# and blit just the redrawn area
canvas.blit(axes.bbox)
def on_motion(self, event):
if dragSpline.lock is not self:
return
if event.inaxes != self.point.axes: return
# self.point.set_xdata, self.point._y, xpress, ypress = self.press
self.point.set_xdata(self.press[0][:, 0])
self.point.set_ydata(self.press[0][:, 1])
xpress = self.press[2]
ypress = self.press[3]
self.points[:, 0] = self.point.get_ydata()
self.points[self.press[1], 0] += event.ydata - ypress
if self.press[1] == 1:
self.points[:, 1] = self.point.get_xdata()
self.points[self.press[1], 1] += event.xdata - xpress
self.point.set_xdata(self.points[:, 1])
self.point.set_ydata(self.points[:, 0])
self.bezierCurve()
canvas = self.point.figure.canvas
axes = self.point.axes
# restore the background region
canvas.restore_region(self.background)
# redraw just the current rectangle
axes.draw_artist(self.point)
axes.draw_artist(self.spline)
# blit just the redrawn area
canvas.blit(axes.bbox)
if self.Signal is not None:
self.Signal.emit()
def on_release(self, event):
'on release we reset the press data'
if dragSpline.lock is not self:
return
self.press = None
dragSpline.lock = None
# turn off the rect animation property and reset the background
self.point.set_animated(False)
self.spline.set_animated(False)
self.background = None
# redraw the full figure
self.point.figure.canvas.draw()
def disconnect(self):
'disconnect all the stored connection ids'
self.point.figure.canvas.mpl_disconnect(self.cidpress)
self.point.figure.canvas.mpl_disconnect(self.cidrelease)
self.point.figure.canvas.mpl_disconnect(self.cidmotion)
def bezierCurve(self):
t = np.linspace(0, 1, 101)
num = np.zeros([101, 2])
dem = np.zeros([101, 2])
n = self.points.shape[0] - 1
for (i, point), weight in zip(enumerate(self.points), self.weights):
biCoeff = binom(n, i)
num = num + ((biCoeff*t**i) * ((1-t)**(n-i)))[:, None] * point * weight
dem = dem + ((biCoeff*t**i) * ((1-t)**(n-i)))[:, None] * weight
self.B = num/dem
self.spline.set_xdata(self.B[:, 1])
self.spline.set_ydata(self.B[:, 0])
include README.rst
\ No newline at end of file
AmpScan ReadME File
--------------------
Author: Joshua Steer
To install the package, open the cmd prompt and type the following script
J:
cd J:\Shared Resources\AmpScan IfLS Team\100 PYTHON\Code
pip install .
To use the package
import AmpScan
import os
read_filename = '01_PhantomShell_ICEM_3mm.stl'
write_filename = '01_PhantomShell_ICEM_3mm_write.stl'
Data = AmpObject(target, 'limb')
Data.addData(baseline, 'socket')
Data.lp_smooth()
Data.man_rot([5,5,5])
Reg = regObject(Data)
Reg.registration(steps=5, baseline='socket', target='limb', reg = 'reglimb')
Reg.save(saveStr, stype='reglimb')
Reg.plot_slices()
Structure:
AmpScan structure
core.py
class AmpObject(Data, stype)
func addData
func read_stl
func unify_vertices
func computeEdges
func save
func calc_norm
func man_trans
func man_rot
func centre
smooth.py
class smoothMixin
func lp_smooth
autoAlign.py
class alignMixin
func icp
func calcDistError
trim.py
class trimMixin
func planarTrim
analyse.py
class analyseMixin
func plot_slices
func create_slice
func planeEdgeIntersect
ampVis.py
class visMixin
func genIm
func addActor
class ampActor
func setVert
func setFaces
func setRect
func setColor
func setOpacity
func setCMap
class vtkRender
class vtkWindow
func render
func setScalarBar
func setViewports
func addAxes
registration.py
class regObject
func registration
TSBSocketDesign.py
class mplCanvas
class dragSpline
func connect
func on_press
func on_motion
func on_release
func disconnect
func bezierCurve
SocketDesignGUI.py
class GUI
func plotRect
func chooseOpenFile
func createActions
func createMenus
AmpScanGUI
class GUI
func chooseOpenFile
func chooseSocket
func align
func register
func analyse
func createActions
func createMenus
setup.py 0 → 100644
# -*- coding: utf-8 -*-
"""
Created on Thu Dec 15 13:43:43 2016
@author: js22g12
"""
from setuptools import setup
def readme():
with open('README.rst') as f:
return f.read()
setup(name='AmpScan',
version='2017.11',
description=('Package for analysis of '
'surface scan data of residual limbs'),
long_description=readme(),
author='Joshua Steer',
author_email='Joshua.Steer@soton.ac.uk',
license='MIT',
packages=['AmpScan',],
install_requires=['numpy', 'pandas', 'matplotlib', 'scipy'],
package_data={'example_stl':' 01_PhantomShell_ICEM_3mm.stl'},
include_package_data=True,
zip_safe=False)
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment