diff --git a/mrBOLD/UI/guimenu/windowMenu.m b/mrBOLD/UI/guimenu/windowMenu.m index 1d73923d..195a8421 100755 --- a/mrBOLD/UI/guimenu/windowMenu.m +++ b/mrBOLD/UI/guimenu/windowMenu.m @@ -72,4 +72,13 @@ uimenu(winmenu, 'Label', 'Open 3D Window', 'Separator', 'on',... 'CallBack', 'open3DWindow;'); +% 3D mrMeshPy submenu. (ADG - Jul2018). +% +%% avoiding this for now so we can get debug info in a terminal window +%%uimenu(winmenu, 'Label', 'Open mrMeshPy 3D Viewer', 'Separator', 'on',... +%% 'CallBack', 'launchMeshPyViewer;'); +uimenu(winmenu, 'Label', 'Open mrMeshPy Control Window', ... + 'CallBack', 'gui_3dWindow_MeshPy;'); + + return \ No newline at end of file diff --git a/mrMesh/meshviewer/meshColorOverlay.m b/mrMesh/meshviewer/meshColorOverlay.m index e8cbdcdb..f9d29089 100755 --- a/mrMesh/meshviewer/meshColorOverlay.m +++ b/mrMesh/meshviewer/meshColorOverlay.m @@ -358,4 +358,4 @@ msh = mrmSet(msh,'colors',newColors'); end -return; +return; \ No newline at end of file diff --git a/mrMeshPy/.gitignore b/mrMeshPy/.gitignore new file mode 100644 index 00000000..66e848bc --- /dev/null +++ b/mrMeshPy/.gitignore @@ -0,0 +1,9 @@ +#python compile +*.pyc + +#matlab in use +*.m~ + +#mac specific +*.DS_store* + diff --git a/mrMeshPy/QVTKRenderWindowInteractor.py b/mrMeshPy/QVTKRenderWindowInteractor.py new file mode 100644 index 00000000..45bdf24b --- /dev/null +++ b/mrMeshPy/QVTKRenderWindowInteractor.py @@ -0,0 +1,635 @@ +#!/usr/bin/python + + +""" +** adapted for mrMeshPy ** AG/MH 2017 + +A simple VTK widget for PyQt or PySide. +See http://www.trolltech.com for Qt documentation, +http://www.riverbankcomputing.co.uk for PyQt, and +http://pyside.github.io for PySide. + +This class is based on the vtkGenericRenderWindowInteractor and is +therefore fairly powerful. It should also play nicely with the +vtk3DWidget code. + +Created by Prabhu Ramachandran, May 2002 +Based on David Gobbi's QVTKRenderWidget.py + +Changes by Gerard Vermeulen Feb. 2003 + Win32 support. + +Changes by Gerard Vermeulen, May 2003 + Bug fixes and better integration with the Qt framework. + +Changes by Phil Thompson, Nov. 2006 + Ported to PyQt v4. + Added support for wheel events. + +Changes by Phil Thompson, Oct. 2007 + Bug fixes. + +Changes by Phil Thompson, Mar. 2008 + Added cursor support. + +Changes by Rodrigo Mologni, Sep. 2013 (Credit to Daniele Esposti) + Bug fix to PySide: Converts PyCObject to void pointer. + +Changes by Greg Schussman, Aug. 2014 + The keyPressEvent function now passes keysym instead of None. + +Changes by Alex Tsui, Apr. 2015 + Port from PyQt4 to PyQt5. + +Changes by Fabian Wenzel, Jan. 2016 + Support for Python3 + + +Changes by Mark Hymers / Andre Gouws October 2017 + Iren initialisation change to stop seg fault +""" + +PyQtImpl = None + +# Check whether a specific PyQt implementation was chosen +try: + import vtk.qt + PyQtImpl = vtk.qt.PyQtImpl +except ImportError: + pass + +# Check whether a specific QVTKRenderWindowInteractor base +# class was chosen, can be set to "QGLWidget" +QVTKRWIBase = "QGLWidget" +try: + import vtk.qt + QVTKRWIBase = vtk.qt.QVTKRWIBase +except ImportError: + pass + +if PyQtImpl is None: + # Autodetect the PyQt implementation to use + try: + import PyQt5 + PyQtImpl = "PyQt5" + except ImportError: + try: + import PyQt4 + PyQtImpl = "PyQt4" + except ImportError: + try: + import PySide + PyQtImpl = "PySide" + except ImportError: + raise ImportError("Cannot load either PyQt or PySide") + +if PyQtImpl == "PyQt5": + if QVTKRWIBase == "QGLWidget": + from PyQt5.QtOpenGL import QGLWidget + from PyQt5.QtWidgets import QWidget + from PyQt5.QtWidgets import QSizePolicy + from PyQt5.QtWidgets import QApplication + from PyQt5.QtCore import Qt + from PyQt5.QtCore import QTimer + from PyQt5.QtCore import QObject + from PyQt5.QtCore import QSize + from PyQt5.QtCore import QEvent +elif PyQtImpl == "PyQt4": + if QVTKRWIBase == "QGLWidget": + from PyQt4.QtOpenGL import QGLWidget + from PyQt4.QtGui import QWidget + from PyQt4.QtGui import QSizePolicy + from PyQt4.QtGui import QApplication + from PyQt4.QtCore import Qt + from PyQt4.QtCore import QTimer + from PyQt4.QtCore import QObject + from PyQt4.QtCore import QSize + from PyQt4.QtCore import QEvent +elif PyQtImpl == "PySide": + if QVTKRWIBase == "QGLWidget": + from PySide.QtOpenGL import QGLWidget + from PySide.QtGui import QWidget + from PySide.QtGui import QSizePolicy + from PySide.QtGui import QApplication + from PySide.QtCore import Qt + from PySide.QtCore import QTimer + from PySide.QtCore import QObject + from PySide.QtCore import QSize + from PySide.QtCore import QEvent +else: + raise ImportError("Unknown PyQt implementation " + repr(PyQtImpl)) + +# Define types for base class, based on string +if QVTKRWIBase == "QWidget": + QVTKRWIBaseClass = QWidget +elif QVTKRWIBase == "QGLWidget": + QVTKRWIBaseClass = QGLWidget +else: + raise ImportError("Unknown base class for QVTKRenderWindowInteractor " + QVTKRWIBase) + +import vtk + +class QVTKRenderWindowInteractor(QVTKRWIBaseClass): + + """ A QVTKRenderWindowInteractor for Python and Qt. Uses a + vtkGenericRenderWindowInteractor to handle the interactions. Use + GetRenderWindow() to get the vtkRenderWindow. Create with the + keyword stereo=1 in order to generate a stereo-capable window. + + The user interface is summarized in vtkInteractorStyle.h: + + - Keypress j / Keypress t: toggle between joystick (position + sensitive) and trackball (motion sensitive) styles. In joystick + style, motion occurs continuously as long as a mouse button is + pressed. In trackball style, motion occurs when the mouse button + is pressed and the mouse pointer moves. + + - Keypress c / Keypress o: toggle between camera and object + (actor) modes. In camera mode, mouse events affect the camera + position and focal point. In object mode, mouse events affect + the actor that is under the mouse pointer. + + - Button 1: rotate the camera around its focal point (if camera + mode) or rotate the actor around its origin (if actor mode). The + rotation is in the direction defined from the center of the + renderer's viewport towards the mouse position. In joystick mode, + the magnitude of the rotation is determined by the distance the + mouse is from the center of the render window. + + - Button 2: pan the camera (if camera mode) or translate the actor + (if object mode). In joystick mode, the direction of pan or + translation is from the center of the viewport towards the mouse + position. In trackball mode, the direction of motion is the + direction the mouse moves. (Note: with 2-button mice, pan is + defined as -Button 1.) + + - Button 3: zoom the camera (if camera mode) or scale the actor + (if object mode). Zoom in/increase scale if the mouse position is + in the top half of the viewport; zoom out/decrease scale if the + mouse position is in the bottom half. In joystick mode, the amount + of zoom is controlled by the distance of the mouse pointer from + the horizontal centerline of the window. + + - Keypress 3: toggle the render window into and out of stereo + mode. By default, red-blue stereo pairs are created. Some systems + support Crystal Eyes LCD stereo glasses; you have to invoke + SetStereoTypeToCrystalEyes() on the rendering window. Note: to + use stereo you also need to pass a stereo=1 keyword argument to + the constructor. + + - Keypress e: exit the application. + + - Keypress f: fly to the picked point + + - Keypress p: perform a pick operation. The render window interactor + has an internal instance of vtkCellPicker that it uses to pick. + + - Keypress r: reset the camera view along the current view + direction. Centers the actors and moves the camera so that all actors + are visible. + + - Keypress s: modify the representation of all actors so that they + are surfaces. + + - Keypress u: invoke the user-defined function. Typically, this + keypress will bring up an interactor that you can type commands in. + + - Keypress w: modify the representation of all actors so that they + are wireframe. + """ + + # Map between VTK and Qt cursors. + _CURSOR_MAP = { + 0: Qt.ArrowCursor, # VTK_CURSOR_DEFAULT + 1: Qt.ArrowCursor, # VTK_CURSOR_ARROW + 2: Qt.SizeBDiagCursor, # VTK_CURSOR_SIZENE + 3: Qt.SizeFDiagCursor, # VTK_CURSOR_SIZENWSE + 4: Qt.SizeBDiagCursor, # VTK_CURSOR_SIZESW + 5: Qt.SizeFDiagCursor, # VTK_CURSOR_SIZESE + 6: Qt.SizeVerCursor, # VTK_CURSOR_SIZENS + 7: Qt.SizeHorCursor, # VTK_CURSOR_SIZEWE + 8: Qt.SizeAllCursor, # VTK_CURSOR_SIZEALL + 9: Qt.PointingHandCursor, # VTK_CURSOR_HAND + 10: Qt.CrossCursor, # VTK_CURSOR_CROSSHAIR + } + + def __init__(self, parent=None, **kw): + # the current button + self._ActiveButton = Qt.NoButton + + # private attributes + self.__saveX = 0 + self.__saveY = 0 + self.__saveModifiers = Qt.NoModifier + self.__saveButtons = Qt.NoButton + self.__wheelDelta = 0 + + # do special handling of some keywords: + # stereo, rw + + try: + stereo = bool(kw['stereo']) + except KeyError: + stereo = False + + try: + rw = kw['rw'] + except KeyError: + rw = None + + # create base qt-level widget + if QVTKRWIBase == "QWidget": + if "wflags" in kw: + wflags = kw['wflags'] + else: + wflags = Qt.WindowFlags() + QWidget.__init__(self, parent, wflags | Qt.MSWindowsOwnDC) + elif QVTKRWIBase == "QGLWidget": + QGLWidget.__init__(self, parent) + + if rw: # user-supplied render window + self._RenderWindow = rw + else: + self._RenderWindow = vtk.vtkRenderWindow() + + WId = self.winId() + + # Python2 + if type(WId).__name__ == 'PyCObject': + from ctypes import pythonapi, c_void_p, py_object + + pythonapi.PyCObject_AsVoidPtr.restype = c_void_p + pythonapi.PyCObject_AsVoidPtr.argtypes = [py_object] + + WId = pythonapi.PyCObject_AsVoidPtr(WId) + + # Python3 + elif type(WId).__name__ == 'PyCapsule': + from ctypes import pythonapi, c_void_p, py_object, c_char_p + + pythonapi.PyCapsule_GetName.restype = c_char_p + pythonapi.PyCapsule_GetName.argtypes = [py_object] + + name = pythonapi.PyCapsule_GetName(WId) + + pythonapi.PyCapsule_GetPointer.restype = c_void_p + pythonapi.PyCapsule_GetPointer.argtypes = [py_object, c_char_p] + + WId = pythonapi.PyCapsule_GetPointer(WId, name) + + self._RenderWindow.SetWindowInfo(str(int(WId))) + + if stereo: # stereo mode + self._RenderWindow.StereoCapableWindowOn() + self._RenderWindow.SetStereoTypeToCrystalEyes() + + try: + self._Iren = kw['iren'] + except KeyError: + self._Iren = vtk.vtkGenericRenderWindowInteractor() + self._Iren.SetRenderWindow(self._RenderWindow) + + # do all the necessary qt setup + self.setAttribute(Qt.WA_OpaquePaintEvent) + self.setAttribute(Qt.WA_PaintOnScreen) + self.setMouseTracking(True) # get all mouse events + self.setFocusPolicy(Qt.WheelFocus) + self.setSizePolicy(QSizePolicy(QSizePolicy.Expanding, QSizePolicy.Expanding)) + + self._Timer = QTimer(self) + self._Timer.timeout.connect(self.TimerEvent) + + self._Iren.AddObserver('CreateTimerEvent', self.CreateTimer) + self._Iren.AddObserver('DestroyTimerEvent', self.DestroyTimer) + self._Iren.GetRenderWindow().AddObserver('CursorChangedEvent', + self.CursorChangedEvent) + + #Create a hidden child widget and connect its destroyed signal to its + #parent ``Finalize`` slot. The hidden children will be destroyed before + #its parent thus allowing cleanup of VTK elements. + self._hidden = QWidget(self) + self._hidden.hide() + self._hidden.destroyed.connect(self.Finalize) + + def __getattr__(self, attr): + """Makes the object behave like a vtkGenericRenderWindowInteractor""" + if attr == '__vtk__': + return lambda t=self._Iren: t + elif hasattr(self._Iren, attr): + return getattr(self._Iren, attr) + else: + raise AttributeError(self.__class__.__name__ + + " has no attribute named " + attr) + + def Finalize(self): + ''' + Call internal cleanup method on VTK objects + ''' + self._RenderWindow.Finalize() + + def CreateTimer(self, obj, evt): + self._Timer.start(10) + + def DestroyTimer(self, obj, evt): + self._Timer.stop() + return 1 + + def TimerEvent(self): + self._Iren.TimerEvent() + + def CursorChangedEvent(self, obj, evt): + """Called when the CursorChangedEvent fires on the render window.""" + # This indirection is needed since when the event fires, the current + # cursor is not yet set so we defer this by which time the current + # cursor should have been set. + QTimer.singleShot(0, self.ShowCursor) + + def HideCursor(self): + """Hides the cursor.""" + self.setCursor(Qt.BlankCursor) + + def ShowCursor(self): + """Shows the cursor.""" + vtk_cursor = self._Iren.GetRenderWindow().GetCurrentCursor() + qt_cursor = self._CURSOR_MAP.get(vtk_cursor, Qt.ArrowCursor) + self.setCursor(qt_cursor) + + def closeEvent(self, evt): + self.Finalize() + + def sizeHint(self): + return QSize(400, 400) + + def paintEngine(self): + return None + + def paintEvent(self, ev): + self._Iren.Render() + + def resizeEvent(self, ev): + w = self.width() + h = self.height() + vtk.vtkRenderWindow.SetSize(self._RenderWindow, w, h) + self._Iren.SetSize(w, h) + self._Iren.ConfigureEvent() + self.update() + + def _GetCtrlShift(self, ev): + ctrl = shift = False + + if hasattr(ev, 'modifiers'): + if ev.modifiers() & Qt.ShiftModifier: + shift = True + if ev.modifiers() & Qt.ControlModifier: + ctrl = True + else: + if self.__saveModifiers & Qt.ShiftModifier: + shift = True + if self.__saveModifiers & Qt.ControlModifier: + ctrl = True + + return ctrl, shift + + def enterEvent(self, ev): + ctrl, shift = self._GetCtrlShift(ev) + self._Iren.SetEventInformationFlipY(self.__saveX, self.__saveY, + ctrl, shift, chr(0), 0, None) + self._Iren.EnterEvent() + + def leaveEvent(self, ev): + ctrl, shift = self._GetCtrlShift(ev) + self._Iren.SetEventInformationFlipY(self.__saveX, self.__saveY, + ctrl, shift, chr(0), 0, None) + self._Iren.LeaveEvent() + + def mousePressEvent(self, ev): + ctrl, shift = self._GetCtrlShift(ev) + repeat = 0 + if ev.type() == QEvent.MouseButtonDblClick: + repeat = 1 + self._Iren.SetEventInformationFlipY(ev.x(), ev.y(), + ctrl, shift, chr(0), repeat, None) + + self._ActiveButton = ev.button() + + if self._ActiveButton == Qt.LeftButton: + self._Iren.LeftButtonPressEvent() + elif self._ActiveButton == Qt.RightButton: + self._Iren.RightButtonPressEvent() + elif self._ActiveButton == Qt.MidButton: + self._Iren.MiddleButtonPressEvent() + + def mouseReleaseEvent(self, ev): + ctrl, shift = self._GetCtrlShift(ev) + self._Iren.SetEventInformationFlipY(ev.x(), ev.y(), + ctrl, shift, chr(0), 0, None) + + if self._ActiveButton == Qt.LeftButton: + self._Iren.LeftButtonReleaseEvent() + elif self._ActiveButton == Qt.RightButton: + self._Iren.RightButtonReleaseEvent() + elif self._ActiveButton == Qt.MidButton: + self._Iren.MiddleButtonReleaseEvent() + + def mouseMoveEvent(self, ev): + self.__saveModifiers = ev.modifiers() + self.__saveButtons = ev.buttons() + self.__saveX = ev.x() + self.__saveY = ev.y() + + ctrl, shift = self._GetCtrlShift(ev) + self._Iren.SetEventInformationFlipY(ev.x(), ev.y(), + ctrl, shift, chr(0), 0, None) + self._Iren.MouseMoveEvent() + + def keyPressEvent(self, ev): + ctrl, shift = self._GetCtrlShift(ev) + if ev.key() < 256: + key = str(ev.text()) + else: + key = chr(0) + + keySym = _qt_key_to_key_sym(ev.key()) + if shift and len(keySym) == 1 and keySym.isalpha(): + keySym = keySym.upper() + + self._Iren.SetEventInformationFlipY(self.__saveX, self.__saveY, + ctrl, shift, key, 0, keySym) + self._Iren.KeyPressEvent() + self._Iren.CharEvent() + + def keyReleaseEvent(self, ev): + ctrl, shift = self._GetCtrlShift(ev) + if ev.key() < 256: + key = chr(ev.key()) + else: + key = chr(0) + + self._Iren.SetEventInformationFlipY(self.__saveX, self.__saveY, + ctrl, shift, key, 0, None) + self._Iren.KeyReleaseEvent() + + def wheelEvent(self, ev): + if hasattr(ev, 'delta'): + self.__wheelDelta += ev.delta() + else: + self.__wheelDelta += ev.angleDelta().y() + + if self.__wheelDelta >= 120: + self._Iren.MouseWheelForwardEvent() + self.__wheelDelta = 0 + elif self.__wheelDelta <= -120: + self._Iren.MouseWheelBackwardEvent() + self.__wheelDelta = 0 + + def GetRenderWindow(self): + return self._RenderWindow + + def Render(self): + self.update() + + +def QVTKRenderWidgetConeExample(): + """A simple example that uses the QVTKRenderWindowInteractor class.""" + + # every QT app needs an app + app = QApplication(['QVTKRenderWindowInteractor']) + + # create the widget + widget = QVTKRenderWindowInteractor() + widget.Initialize() + widget.Start() + # if you dont want the 'q' key to exit comment this. + widget.AddObserver("ExitEvent", lambda o, e, a=app: a.quit()) + + ren = vtk.vtkRenderer() + widget.GetRenderWindow().AddRenderer(ren) + + cone = vtk.vtkConeSource() + cone.SetResolution(8) + + coneMapper = vtk.vtkPolyDataMapper() + coneMapper.SetInputConnection(cone.GetOutputPort()) + + coneActor = vtk.vtkActor() + coneActor.SetMapper(coneMapper) + + ren.AddActor(coneActor) + + # show the widget + widget.show() + # start event processing + app.exec_() + + +_keysyms = { + Qt.Key_Backspace: 'BackSpace', + Qt.Key_Tab: 'Tab', + Qt.Key_Backtab: 'Tab', + # Qt.Key_Clear : 'Clear', + Qt.Key_Return: 'Return', + Qt.Key_Enter: 'Return', + Qt.Key_Shift: 'Shift_L', + Qt.Key_Control: 'Control_L', + Qt.Key_Alt: 'Alt_L', + Qt.Key_Pause: 'Pause', + Qt.Key_CapsLock: 'Caps_Lock', + Qt.Key_Escape: 'Escape', + Qt.Key_Space: 'space', + # Qt.Key_Prior : 'Prior', + # Qt.Key_Next : 'Next', + Qt.Key_End: 'End', + Qt.Key_Home: 'Home', + Qt.Key_Left: 'Left', + Qt.Key_Up: 'Up', + Qt.Key_Right: 'Right', + Qt.Key_Down: 'Down', + Qt.Key_SysReq: 'Snapshot', + Qt.Key_Insert: 'Insert', + Qt.Key_Delete: 'Delete', + Qt.Key_Help: 'Help', + Qt.Key_0: '0', + Qt.Key_1: '1', + Qt.Key_2: '2', + Qt.Key_3: '3', + Qt.Key_4: '4', + Qt.Key_5: '5', + Qt.Key_6: '6', + Qt.Key_7: '7', + Qt.Key_8: '8', + Qt.Key_9: '9', + Qt.Key_A: 'a', + Qt.Key_B: 'b', + Qt.Key_C: 'c', + Qt.Key_D: 'd', + Qt.Key_E: 'e', + Qt.Key_F: 'f', + Qt.Key_G: 'g', + Qt.Key_H: 'h', + Qt.Key_I: 'i', + Qt.Key_J: 'j', + Qt.Key_K: 'k', + Qt.Key_L: 'l', + Qt.Key_M: 'm', + Qt.Key_N: 'n', + Qt.Key_O: 'o', + Qt.Key_P: 'p', + Qt.Key_Q: 'q', + Qt.Key_R: 'r', + Qt.Key_S: 's', + Qt.Key_T: 't', + Qt.Key_U: 'u', + Qt.Key_V: 'v', + Qt.Key_W: 'w', + Qt.Key_X: 'x', + Qt.Key_Y: 'y', + Qt.Key_Z: 'z', + Qt.Key_Asterisk: 'asterisk', + Qt.Key_Plus: 'plus', + Qt.Key_Minus: 'minus', + Qt.Key_Period: 'period', + Qt.Key_Slash: 'slash', + Qt.Key_F1: 'F1', + Qt.Key_F2: 'F2', + Qt.Key_F3: 'F3', + Qt.Key_F4: 'F4', + Qt.Key_F5: 'F5', + Qt.Key_F6: 'F6', + Qt.Key_F7: 'F7', + Qt.Key_F8: 'F8', + Qt.Key_F9: 'F9', + Qt.Key_F10: 'F10', + Qt.Key_F11: 'F11', + Qt.Key_F12: 'F12', + Qt.Key_F13: 'F13', + Qt.Key_F14: 'F14', + Qt.Key_F15: 'F15', + Qt.Key_F16: 'F16', + Qt.Key_F17: 'F17', + Qt.Key_F18: 'F18', + Qt.Key_F19: 'F19', + Qt.Key_F20: 'F20', + Qt.Key_F21: 'F21', + Qt.Key_F22: 'F22', + Qt.Key_F23: 'F23', + Qt.Key_F24: 'F24', + Qt.Key_NumLock: 'Num_Lock', + Qt.Key_ScrollLock: 'Scroll_Lock', + } + +def _qt_key_to_key_sym(key): + """ Convert a Qt key into a vtk keysym. + + This is essentially copied from the c++ implementation in + GUISupport/Qt/QVTKInteractorAdapter.cxx. + """ + + if key not in _keysyms: + return None + + return _keysyms[key] + + +if __name__ == "__main__": + print(PyQtImpl) + QVTKRenderWidgetConeExample() diff --git a/mrMeshPy/README.md b/mrMeshPy/README.md new file mode 100644 index 00000000..ed0be735 --- /dev/null +++ b/mrMeshPy/README.md @@ -0,0 +1,25 @@ +# mrMeshPy +A Python, VTK6 and Qt5 dependent viewer to replace mrMesh (Stanford). + +Instruction video can be see at https://youtu.be/NH_wEimGxVQ + +This Readme should be read in conjunction with the License.md file in this directory. + +This program is reliant on the following packages: + - Python VTK 6 + - Python Qt 5 + - numpy + - scipy + +HEALTH WARNING -- Current caveats for users familiar with mrMesh + +As you can imagine, this program is very much a work in progress. You should be aware of the following current limitations of the program as it stands (and how it differs from mrMesh): + +- - Support for VOLUME{1} only – mrMeshPy does not currently support loading of multiple VOLUME views and will always try to reference VOLUME{1} in the matlab workspace. + +- - Crash = restart – if for some reason either matlab or mrMeshPy crash, please restart both matlab and mrMeshPy or things might get horribly out of synch – we are working on this :) + +- - The menu bar (top) in the viewer is incomplete, but the ROIs submenu works + +- - Other functions are slowly being added – see/edit TODO.txt for info. + diff --git a/mrMeshPy/README.txt b/mrMeshPy/README.txt new file mode 100644 index 00000000..17a9a5e5 --- /dev/null +++ b/mrMeshPy/README.txt @@ -0,0 +1,2 @@ +Test upload to get persmissions set correctly +ADG diff --git a/mrMeshPy/TODO.txt b/mrMeshPy/TODO.txt new file mode 100644 index 00000000..c82fa2c5 --- /dev/null +++ b/mrMeshPy/TODO.txt @@ -0,0 +1,23 @@ +# TODO.txt +mrMeshPy TODO / wishlist + +Priority(0 crucial, 1 high, 2 low, 3 wishlist) -- Action -- Assigned to + +(1) -- enable different TCP adresses - currently hard-coded to 9999 -- ADG +(1) -- improve error handling -- ADG / MH +(1) -- improve TP comms (error state, etc) -- ADG / MH +(1) -- review meshBuild in matlab (i.e. not vtk dependent) -- ADG / AW +(1) -- review meshBuild in mrMeshPy (i.e. IS vtk dependent, but will be built into viewer) -- ADG / AW + +(2) -- ROI transforms to vista support currently only support "ALL LAYERS" -- AG +(2) -- add exisiting function: SAVE CAMERA -- AG +(2) -- add exisiting function: SET CAMERA -- AG +(2) -- improve tracking of mesh instances between matlab and mrMeshPy in case of crash -- AG +(2) -- add support for additional VOLUME instances in matlab -- AG +(2) -- fix parent UI callback global in menu code -- AG +(2) -- review timeouts in TCP server -- AG / MH + +(3) -- add “current session save” -- AG +(3) -- incorporate code for handle removal - AG +(3) -- incorporate code for removing “spikes” (have this in DV3D) - AG +(3) -- incorporate code for surface clipping (have this in DV3D) - AG diff --git a/mrMeshPy/matlabRoutines/createMrMeshPySendTCPCommand.m b/mrMeshPy/matlabRoutines/createMrMeshPySendTCPCommand.m new file mode 100755 index 00000000..fde2f73b --- /dev/null +++ b/mrMeshPy/matlabRoutines/createMrMeshPySendTCPCommand.m @@ -0,0 +1,50 @@ +function cmdStringTCP = createMrMeshPySendTCPCommand(cmdStruct) +% function returnOutput = createMrMeshPySendCommand(cmdStruct) +% This function will generate a TCP compatible command string for use +% with the mrMeshPy interface +% +%It requires a structure as input, with the following attributes: +% cmdStruct.CommandToSend should be on of hte possible commands listed in +% mrMeshPyCommandsList.m +% cmdStruct.TargetMeshSession should be an integer pointing to the mesh +% instance the command is referenced to +% cmdStruct.Args is a cell array with a number of string in it: each +% string is a comma-separated set of values or descriptors +% and describe the data / actions mrMeshpy should expect in +% the subsequent transactions +% Andre' Gouws 2017 + +% the code below incrementally builds a TCP command string from different +% command components - it may appear overkill to do it this way but I just +% wanted to be transparent about how the command is constructed for future +% developers + +cmdString = 'cmd'; % all commands start with a 'cmd' here +cmdString = [cmdString, ';']; % seperator + +cmdString = [cmdString, 'None']; % a placeholder - we almost certainly forgot something obvious so lets leave a gap ;) +cmdString = [cmdString, ';']; % seperator + +cmdString = [cmdString, 'd']; % a dataType placeholder - for the placeholder above ;) +cmdString = [cmdString, ';']; % seperator + +cmdString = [cmdString, 'None']; % a data counter for the place holder above +cmdString = [cmdString, ';']; % seperator + +% Place holders done, lets actually add our data to the string +cmdString = [cmdString, cmdStruct.CommandToSend]; % the command to send - must be in mrMeshPyCommandsList +cmdString = [cmdString, ';']; % seperator + +% major change here - no longer an instance number but a unique string +% identifier +cmdString = [cmdString, num2str(cmdStruct.TargetMeshSession)]; % unique identifier of mrMesh sub-window +cmdString = [cmdString, ';']; % seperator + +% Process the extra arg strings by looping through our Arg cell array + +for i = 1:length(cmdStruct.Args) + cmdString = [cmdString, cmdStruct.Args{i}]; % instance number of mrMesh sub-window + cmdString = [cmdString, ';']; % seperator +end + +cmdStringTCP = cmdString; diff --git a/mrMeshPy/matlabRoutines/gui_3dWindow_MeshPy.fig b/mrMeshPy/matlabRoutines/gui_3dWindow_MeshPy.fig new file mode 100755 index 00000000..5352062f Binary files /dev/null and b/mrMeshPy/matlabRoutines/gui_3dWindow_MeshPy.fig differ diff --git a/mrMeshPy/matlabRoutines/gui_3dWindow_MeshPy.m b/mrMeshPy/matlabRoutines/gui_3dWindow_MeshPy.m new file mode 100755 index 00000000..4e95f8b4 --- /dev/null +++ b/mrMeshPy/matlabRoutines/gui_3dWindow_MeshPy.m @@ -0,0 +1,405 @@ +function varargout = gui_3dWindow_MeshPy(varargin) +% GUI_3DWINDOW_MESHPY MATLAB code for gui_3dWindow_MeshPy.fig +% GUI_3DWINDOW_MESHPY, by itself, creates a new GUI_3DWINDOW_MESHPY or raises the existing +% singleton*. +% +% H = GUI_3DWINDOW_MESHPY returns the handle to a new GUI_3DWINDOW_MESHPY or the handle to +% the existing singleton*. +% +% GUI_3DWINDOW_MESHPY('CALLBACK',hObject,eventData,handles,...) calls the local +% function named CALLBACK in GUI_3DWINDOW_MESHPY.M with the given input arguments. +% +% GUI_3DWINDOW_MESHPY('Property','Value',...) creates a new GUI_3DWINDOW_MESHPY or raises the +% existing singleton*. Starting from the left, property value pairs are +% applied to the GUI before gui_3dWindow_MeshPy_OpeningFcn gets called. An +% unrecognized property name or invalid value makes property application +% stop. All inputs are passed to gui_3dWindow_MeshPy_OpeningFcn via varargin. +% +% *See GUI Options on GUIDE's Tools menu. Choose "GUI allows only one +% instance to run (singleton)". +% +% See also: GUIDE, GUIDATA, GUIHANDLES + +% Edit the above text to modify the response to help gui_3dWindow_MeshPy + +% Last Modified by GUIDE v2.5 26-Sep-2017 13:56:57 +% Andre' Gouws 2017 +debug = 0; + +if debug + pwd + fileparts(mfilename('fullpath')) +end + +% If this function is called, we are going to assume that the user wants to +% use mrMeshPy and not mrMesh - we we will add some altered routines to the +% top of the search path + +mrMeshMFileDir = fileparts(mfilename('fullpath')); %get the directory that holds THIS script +addpath(genpath(mrMeshMFileDir)); % add it to the top of the search path. +% TODO - there must be a better way of doing this - maybe just drop mrMesh! + + +% Begin initialization code - DO NOT EDIT +gui_Singleton = 1; +gui_State = struct('gui_Name', mfilename, ... + 'gui_Singleton', gui_Singleton, ... + 'gui_OpeningFcn', @gui_3dWindow_MeshPy_OpeningFcn, ... + 'gui_OutputFcn', @gui_3dWindow_MeshPy_OutputFcn, ... + 'gui_LayoutFcn', [] , ... + 'gui_Callback', []); +if nargin && ischar(varargin{1}) + gui_State.gui_Callback = str2func(varargin{1}); +end + +if nargout + [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:}); +else + gui_mainfcn(gui_State, varargin{:}); +end +% End initialization code - DO NOT EDIT + +try + % Hack - TODO fix - this puts the VOLUME in the scope of this gui + VOLUME = evalin('base','VOLUME','VOLUME'); %% HACK TODO fix +catch + disp('no VOLUME loaded yet') +end + +% --- Executes just before gui_3dWindow_MeshPy is made visible. +function gui_3dWindow_MeshPy_OpeningFcn(hObject, eventdata, handles, varargin) +% This function has no output args, see OutputFcn. +% hObject handle to figure +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) +% varargin command line arguments to gui_3dWindow_MeshPy (see VARARGIN) + +% Choose default command line output for gui_3dWindow_MeshPy +handles.output = hObject; + +% Update handles structure +guidata(hObject, handles); + +% UIWAIT makes gui_3dWindow_MeshPy wait for user response (see UIRESUME) +% uiwait(handles.figure1); + + +% --- Outputs from this function are returned to the command line. +function varargout = gui_3dWindow_MeshPy_OutputFcn(hObject, eventdata, handles) +% varargout cell array for returning output args (see VARARGOUT); +% hObject handle to figure +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) + +% Get default command line output from handles structure +varargout{1} = handles.output; + + + +%% disabled for now - TODO +% --- Executes on button press in launch_button. +function launch_button_Callback(hObject, eventdata, handles) + +[myfile,mydir]= uigetfile({'*.mat','MAT-files (*.mat)'}); +meshFilePath = [mydir,myfile]; +myPid = num2str(feature('getpid')); + +if ~isfield(VOLUME{1},'meshNum3d') %% TODO - this willneed to reflect the current volume and x-ref the correct mesh + meshInstance = '1'; +else + meshInstance = num2str(VOLUME{1}.meshNum3d + 1); +end + +evalstr = ['/home/andre/mrMeshPy/launchMeshPy.sh /home/andre/mrMeshPy/meshPy_v03.py ',meshFilePath,' ',myPid,' ',meshInstance,' &']; +disp(['ran command: ', evalstr]); + +system(evalstr); +%% + + + +% --- Executes on button press in pushbutton_update. +function pushbutton_update_Callback(hObject, eventdata, handles) + +mrGlobals; +set( findall(handles.uipanel1, '-property', 'Enable'), 'Enable', 'off') + +try + currMesh = VOLUME{1}.meshNum3d; + + [VOLUME{1},~,~,~,VOLUME{1}.mesh{currMesh}] = meshColorOverlay(VOLUME{1},0); + mrMeshPySend('updateMeshData',VOLUME{1}); +catch + disp 'error in update mesh routine'; +end + +set( findall(handles.uipanel1, '-property', 'Enable'), 'Enable', 'On') + + + +% --- Executes on button press in pushbutton_LoadMesh. +function pushbutton_LoadMesh_Callback(hObject, eventdata, handles) +% hObject handle to pushbutton_LoadMesh (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) + +mrGlobals; + +% keep tack of whether the VOLUME actually changes - i.e. make sure a new +% mesh has been loaded and the user hasn't hit cancel +try + currMeshCount = length(VOLUME{1}.mesh); +catch + currMeshCount = 0; +end + +% load the mesh to the VOLUME struct +VOLUME{1} = meshLoad(VOLUME{1},'./'); %start in the current directory for now %TODO later give options? + +try length(VOLUME{1}.mesh) + + if length(VOLUME{1}.mesh) > currMeshCount %a new mesh has been added + + % create a unique ID for the mesh based on a timestamp (clock) + VOLUME{1}.mesh{VOLUME{1}.meshNum3d}.mrMeshPyID = makeUniqueID; + + % send the newly loaded mesh to the viewer via the VOLUME + mrMeshPySend('sendNewMeshData',VOLUME{1}); + + handles = guidata(hObject); % Update! + currString = get(handles.popupmenu_Meshes,'string') + + if strcmp(currString,'None') + newstring = char(['mesh-',VOLUME{1}.mesh{VOLUME{1}.meshNum3d}.mrMeshPyID]); + else + newstring = char(currString,['mesh-',VOLUME{1}.mesh{VOLUME{1}.meshNum3d}.mrMeshPyID]); + end + %disp 'here1' + %VOLUME{1}.meshNum3d + set(handles.popupmenu_Meshes,'value',VOLUME{1}.meshNum3d) ; + set(handles.popupmenu_Meshes,'string',newstring); + %disp 'here2' + else % no new mesh added + disp('User cancelled mesh load or there was an error loading ...'); + end + + +catch % no new mesh added + disp('User cancelled mesh load or there was an error loading ...'); +end + + + +% --- Executes on selection change in popupmenu_Meshes. +function popupmenu_Meshes_Callback(hObject, eventdata, handles) +% hObject handle to popupmenu_Meshes (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) + +% Hints: contents = cellstr(get(hObject,'String')) returns popupmenu_Meshes contents as cell array +% contents{get(hObject,'Value')} returns selected item from popupmenu_Meshes + +mrGlobals; + +% % %assignin('base','hObj',hObject) +% % %contents = cellstr(get(hObject,'String')) +% % %meshNum = contents{get(hObject,'Value')} + +meshNum = hObject.Value; %should be the index + +%%%meshNum = meshNum(5:end); %TODO improve +VOLUME{1}.meshNum3d = meshNum; + + +% --- Executes during object creation, after setting all properties. +function popupmenu_Meshes_CreateFcn(hObject, eventdata, handles) +% hObject handle to popupmenu_Meshes (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles empty - handles not created until after all CreateFcns called + +% Hint: popupmenu controls usually have a white background on Windows. +% See ISPC and COMPUTER. +if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) + set(hObject,'BackgroundColor','white'); +end + + +% --- Executes on button press in pushbutton_getROI. +function pushbutton_getROI_Callback(hObject, eventdata, handles) +% hObject handle to pushbutton_getROI (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) + +mrGlobals; + +mrMeshPySend('checkMeshROI',VOLUME{1}); + + +% --- Executes on button press in pushbutton_smooth. +function pushbutton_smooth_Callback(hObject, eventdata, handles) +% hObject handle to pushbutton_smooth (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) + +mrGlobals; + +handles = guidata(hObject); % Update! +relax = get(handles.edit_relaxationFactor,'String') +iterations = get(handles.edit_iterations,'String') + +relax = str2num(relax); +iterations = str2num(iterations); + +%assignin('base','iterations',iterations); +%assignin('base','relax',relax); + +currMeshID = VOLUME{1}.mesh{VOLUME{1}.meshNum3d}.mrMeshPyID; + +%disp('here1') + +% send (with VOLUME also) +mrMeshPySend('smoothMesh',{currMeshID,iterations,relax,VOLUME{1}}); + + +function edit_relaxationFactor_Callback(hObject, eventdata, handles) +% hObject handle to edit_relaxationFactor (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) + +% Hints: get(hObject,'String') returns contents of edit_relaxationFactor as text +% str2double(get(hObject,'String')) returns contents of edit_relaxationFactor as a double + + +% --- Executes during object creation, after setting all properties. +function edit_relaxationFactor_CreateFcn(hObject, eventdata, handles) +% hObject handle to edit_relaxationFactor (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles empty - handles not created until after all CreateFcns called + +% Hint: edit controls usually have a white background on Windows. +% See ISPC and COMPUTER. +if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) + set(hObject,'BackgroundColor','white'); +end + + + +function edit4_Callback(hObject, eventdata, handles) +% hObject handle to edit4 (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) + +% Hints: get(hObject,'String') returns contents of edit4 as text +% str2double(get(hObject,'String')) returns contents of edit4 as a double + + +% --- Executes during object creation, after setting all properties. +function edit4_CreateFcn(hObject, eventdata, handles) +% hObject handle to edit4 (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles empty - handles not created until after all CreateFcns called + +% Hint: edit controls usually have a white background on Windows. +% See ISPC and COMPUTER. +if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) + set(hObject,'BackgroundColor','white'); +end + + + +function edit_iterations_Callback(hObject, eventdata, handles) +% hObject handle to edit_iterations (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) + +% Hints: get(hObject,'String') returns contents of edit_iterations as text +% str2double(get(hObject,'String')) returns contents of edit_iterations as a double + + +% --- Executes during object creation, after setting all properties. +function edit_iterations_CreateFcn(hObject, eventdata, handles) +% hObject handle to edit_iterations (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles empty - handles not created until after all CreateFcns called + +% Hint: edit controls usually have a white background on Windows. +% See ISPC and COMPUTER. +if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) + set(hObject,'BackgroundColor','white'); +end + + +% --- Executes on button press in pushbutton_buildMesh. +function pushbutton_buildMesh_Callback(hObject, eventdata, handles) +% hObject handle to pushbutton_buildMesh (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) + +mrGlobals; +handles = guidata(hObject) +set( findall(handles.uipanel1, '-property', 'Enable'), 'Enable', 'off') + +%try + % get user option or left or right + [s,v] = listdlg('PromptString','Select hemisphere:',... + 'SelectionMode','single',... + 'ListString',{'left','right'}); + + if s == 1 + hemi = 'left'; + elseif s == 2 + hemi = 'right'; + else + disp('error selecting mesh'); + return + end + + % build the mesh + VOLUME{1} = meshBuild_mrMeshPy(VOLUME{1}, hemi); + + % ask if we would like to load to mrMeshPy? + + % Include the desired Default answer + options.Interpreter = 'tex'; + options.Default = 'Yes'; + % Use the TeX interpreter in the question + qstring = 'Would you like to load the mesh to mrMeshPy?'; + loadNow = questdlg(qstring,'Mesh ready .. load?',... + 'Yes','No',options) + + if strcmp(loadNow,'No') + return + else + %assume yes + + % create a unique ID for the mesh based on a timestamp (clock) + VOLUME{1}.mesh{VOLUME{1}.meshNum3d}.mrMeshPyID = makeUniqueID; + + % send the newly loaded mesh to the viewer + mrMeshPySend('sendNewMeshData',VOLUME{1}); + + handles = guidata(hObject); % Update! + currString = get(handles.popupmenu_Meshes,'string') + + if strcmp(currString,'None') + newstring = char(['mesh-',VOLUME{1}.mesh{VOLUME{1}.meshNum3d}.mrMeshPyID]); + else + newstring = char(currString,['mesh-',VOLUME{1}.mesh{VOLUME{1}.meshNum3d}.mrMeshPyID]); + end + + set(handles.popupmenu_Meshes,'value',VOLUME{1}.meshNum3d) ; + set(handles.popupmenu_Meshes,'string',newstring); + end + +%catch +% disp 'error in build mesh routine'; +%end + +set( findall(handles.uipanel1, '-property', 'Enable'), 'Enable', 'On') + + + + + + + diff --git a/mrMeshPy/matlabRoutines/launchMeshBuild.sh b/mrMeshPy/matlabRoutines/launchMeshBuild.sh new file mode 100755 index 00000000..1cc758b3 --- /dev/null +++ b/mrMeshPy/matlabRoutines/launchMeshBuild.sh @@ -0,0 +1,10 @@ +#!/bin/bash + +''' +# wrapper script to launch mrMeshPy from MatLab TODO implement this later +Mark Hymers, Andre Gouws 2017 +''' + +unset LD_LIBRARY_PATH +exec "$@" + diff --git a/mrMeshPy/matlabRoutines/launchMeshPy.m b/mrMeshPy/matlabRoutines/launchMeshPy.m new file mode 100644 index 00000000..66a646ef --- /dev/null +++ b/mrMeshPy/matlabRoutines/launchMeshPy.m @@ -0,0 +1,10 @@ +#!/bin/bash + +''' +# wrapper script to launch mrMeshPy from MatLab TODO implement this later +Mark Hymers, Andre Gouws 2017 +''' + +unset LD_LIBRARY_PATH +exec /usr/bin/python "$@" & + diff --git a/mrMeshPy/matlabRoutines/launchMeshPy.sh b/mrMeshPy/matlabRoutines/launchMeshPy.sh new file mode 100755 index 00000000..66a646ef --- /dev/null +++ b/mrMeshPy/matlabRoutines/launchMeshPy.sh @@ -0,0 +1,10 @@ +#!/bin/bash + +''' +# wrapper script to launch mrMeshPy from MatLab TODO implement this later +Mark Hymers, Andre Gouws 2017 +''' + +unset LD_LIBRARY_PATH +exec /usr/bin/python "$@" & + diff --git a/mrMeshPy/matlabRoutines/launchMeshPyViewer.m b/mrMeshPy/matlabRoutines/launchMeshPyViewer.m new file mode 100644 index 00000000..628e77ff --- /dev/null +++ b/mrMeshPy/matlabRoutines/launchMeshPyViewer.m @@ -0,0 +1,9 @@ +% a launcher for mrMeshPy viewer from within matlab + +%% new code +% get directory of matlab routine - same place as pyMeshBuild to call later +meshBuildPath = which(mfilename); +[meshBuildDir,~,~] = fileparts(meshBuildPath); + +cmdString = [meshBuildDir,'/launchMeshPy.sh ',meshBuildDir,'/../mrMeshPy.py'] %TODO set python path? +system(cmdString); diff --git a/mrMeshPy/matlabRoutines/makeUniqueID.m b/mrMeshPy/matlabRoutines/makeUniqueID.m new file mode 100755 index 00000000..d7331a65 --- /dev/null +++ b/mrMeshPy/matlabRoutines/makeUniqueID.m @@ -0,0 +1,11 @@ +function uID = makeUniqueID; +% function uID = makeUniqueID(clockString) +% +% takes the system clock and creates a unique (hopefully) string to label +% the mesh and allow us to keep track between processes + +a = num2str(fix(clock)); % get clock string +a = strrep(a,' ','_'); % tidy it up 1 +a = strrep(a,'__','_'); % ang again + +uID = a; diff --git a/mrMeshPy/matlabRoutines/meshBuildFromClass_mrMeshPy.m b/mrMeshPy/matlabRoutines/meshBuildFromClass_mrMeshPy.m new file mode 100755 index 00000000..0b876354 --- /dev/null +++ b/mrMeshPy/matlabRoutines/meshBuildFromClass_mrMeshPy.m @@ -0,0 +1,144 @@ +function [msh,class] = meshBuildFromClass_mrMeshPy(voxels,mmPerVox,hemisphere) +% Build a VTK mesh from a class file +% +% [msh,class] = meshBuildFromClass_mrMeshPy(voxels,[mmPerVox],[hemisphere='left']) +% +% voxels: Either +% - the file name of a white matter class file, or +% - the voxel classification data returned from readClassFile (see code) +% If the file name (class file or nifti file) is entered, then the +% classification data in the file can be returned. +% +% mmPerVox: defaults to [1 1 1] (mm) +% hemisphere: 'right' or 'left' or 'both' +% +% White matter values in the class or NIFTI file are the voxels with the +% value 16. +% +% See also: meshBuild, meshVisualize, +% +% Examples: +% fName=fullfile(mrvDataRootPath,'anatomy','anatomyV','left','left.Class'); +% msh = meshBuildFromClass(fName); +% msh = meshSmooth(msh); +% msh = meshColor(msh); +% meshVisualize(msh); +% +% fName ='X:\anatomy\nakadomari\right\20050901_fixV1\right.Class'; +% mmPerVox = [1 1 1]; +% msh = meshBuildFromClass(fName, mmPerVox, 'right'); +% msh = meshSmooth(msh); +% msh = meshColor(msh); +% meshVisualize(msh); +% +% Guillaume Bertello (c) Stanford VISTA Team + +% PROGRAMMING TODO: +% Perhaps we should replace this function with the Matlab isosurface +% routine. See mrmBuildMeshMatlab. In general, We would like to get rid +% of the VTK dependent mex-files for smooth, build, and color. These can +% be replaced with matlab mesh functions. We want to reduce the mex file +% stuff to only the mrMesh related calls to the window. +% +% We would like the mesh vertices to coregister with the vAnatomy or NIFTI +% T1 data. We need to understand this better. +% + +if ieNotDefined('mmPerVox'), mmPerVox = [1 1 1]; end +if ieNotDefined('hemisphere'), hemisphere = 'left'; end +if isempty(voxels) || ischar(voxels) + switch hemisphere + case 'both' + fprintf('[%s]: Loading %s hemisphere white matter voxels...\n', mfilename),'right'; + headerOnly = 0; voiOnly = 0; + class = readClassFile(voxels,headerOnly,voiOnly,'right'); + voxelsR = uint8(class.data == class.type.white); + + fprintf('[%s]: Loading %s hemisphere white matter voxels...\n', mfilename,'left'); + headerOnly = 0; voiOnly = 0; + class = readClassFile(voxels,headerOnly,voiOnly,'left'); + voxelsL = uint8(class.data == class.type.white); + + voxels = voxelsL + voxelsR; + case {'left','right'} + fprintf('[%s]: Loading %s hemisphere white matter voxels...\n', mfilename,hemisphere); + headerOnly = 0; voiOnly = 0; + class = readClassFile(voxels,headerOnly,voiOnly,hemisphere); + voxels = uint8(class.data == class.type.white); + otherwise + error('[%s]: Unknown hemisphere label', mfilename) + end +elseif(isstruct(voxels)) + voxels = uint8(voxels.data == voxels.type.white); +end + +%% AG EDIT -TODO +%assignin('base','voxels',voxels); + +% build_mesh is a dll in VISTASRC. It converts classification data into +% vertices and triangles. It could be replaced by the Matlab isosurface +% routine. +fprintf('[%s]: Building a %s hemisphere mesh ...', mfilename, hemisphere) +%%% NEXT LINE REMOVED - now use mrMeshPy's pyMeshBuild instead +%%% msh = build_mesh(voxels,mmPerVox); % Vertices (class) are in mm space +%%% + +%% new code +% get directory of matlab routine - same place as pyMeshBuild to call later +meshBuildPath = which(mfilename); +[meshBuildDir,~,~] = fileparts(meshBuildPath); + +%reshape the voxel array +voxels = permute(voxels,[3,2,1]); + +% create a tmp file to write the data to +voxFileForMrMeshPy = [tempname,'.mat']; +mshFileFromMrMeshPy = [tempname,'.mat']; + + +%run the pyMeshBuild.py program to generate the meshes +if ismac + % Mac - same as linux? #TODO + % save the voxel data to a tmp file + eval(['save ',voxFileForMrMeshPy,' voxels mmPerVox;']); + %run the pyMeshBuild app + cmdString = [meshBuildDir,'/launchMeshBuild.sh ',meshBuildDir,'/pyMeshBuild_mac.py ',voxFileForMrMeshPy,' ',mshFileFromMrMeshPy] %TODO set python path? + system(cmdString); + +elseif isunix + % Linux + % save the voxel data to a tmp file + eval(['save ',voxFileForMrMeshPy,' voxels mmPerVox;']); + %run the pyMeshBuild app + cmdString = [meshBuildDir,'/launchMeshBuild.sh ',meshBuildDir,'/pyMeshBuild.py ',voxFileForMrMeshPy,' ',mshFileFromMrMeshPy] + system(cmdString); + +elseif ispc + % Windows + disp('Platform not supported') + barf for now +else + disp('Platform not supported') + return +end + +msh = load(mshFileFromMrMeshPy); + +% AG EDIT -TODO +assignin('base','msh1',msh); + +msh = meshFormat(msh); % Converts old format to new. + +% AG EDIT -TODO +assignin('base','msh2',msh); + +% Set the mesh origin, by default, to the center of the object. +vertices = meshGet(msh,'vertices'); +msh = meshSet(msh,'origin',-mean(vertices,2)'); +msh = meshSet(msh,'mmPerVox',mmPerVox); +fprintf('[%s]: done. \n', mfilename); + + +assignin('base','msh3',msh); + +return; diff --git a/mrMeshPy/matlabRoutines/meshBuild_mrMeshPy.m b/mrMeshPy/matlabRoutines/meshBuild_mrMeshPy.m new file mode 100755 index 00000000..17579940 --- /dev/null +++ b/mrMeshPy/matlabRoutines/meshBuild_mrMeshPy.m @@ -0,0 +1,246 @@ +function [vw,newMeshNum] = meshBuild_mrMeshPy(vw,hemisphere) +% COPY OF ORIGINAL meshBuild adapted for mrMeshPy - AG 2017 +% Build a 3D mesh for visualization and analysis +% +% [vw,newMeshNum,meshBuildPath] = meshBuild(vw,[hemisphere],meshBuildPath); +% +% Using mrVista data, build a mesh, save it in a file in the anatomy +% directory, and add the mesh to the 3D Control Window pull down +% options. +% +% vw: A VISTASOFT view structure +% hemisphere: left, right or both. Default: left +% +% A mrMeshPy window is opened as well, showing the computed mesh. +% +% Example: +% [VOLUME{1},newMeshNum,meshBuildPath] = meshBuild(VOLUME{1},'left','/home/user/mrMeshPy/matlabRoutines'); +% VOLUME{1} = viewSet(VOLUME{1},'currentmeshn',newMeshNum); +% +% See also: meshBuildFromClass, meshBuildFromNiftiClass, meshSmooth, +% meshColor +% +% 11/05 ras: also saves mesh path. +% +% (c) Stanford VISTA Team 2008 +% +% +% Programming TODO. Check this! +% We have (or had?) trouble for building 'both' meshes. We need a new +% procedure. +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% COPY OF ORIGINAL meshBuild adapted for mrMeshPy - AG 2017 - in trying to +% keep as many things consistent with the previous mrMesh, I have stripped +% out this part of the stream and changed the bits we need for mrMeshPy +% This cascade calls adapted versions of mrmBuild and meshBuildFromClass +% Andre Gouws 2017 + + + +disp 'meshBuild called' + +% Be sure anatomy is loaded (we need it for the mmPerVox field) +if isempty(vw.anat), vw = loadAnat(vw); end +if ieNotDefined('hemisphere'), hemisphere = 'left'; end + +newMeshNum = viewGet(vw,'nmesh') + 1; + +% Parameters we establish in this routine +[meshName,numGrayLayers,hemiNum] = readParams(newMeshNum,hemisphere); +if isempty(meshName), newMeshNum = newMeshNum - 1; return; end % User pressed cancel. + +% mmPerVox = viewGet(vw,'mmPerVoxel'); + +wbar = mrvWaitbar(0.1, ... + sprintf('meshBuild: Combining white and gray matter...')); + +% Load left, right, or both hemispheres. +if (hemiNum==1) + [voxels,vw] = meshGetWhite(vw, 'left', numGrayLayers); +elseif (hemiNum==2) + [voxels,vw] = meshGetWhite(vw, 'right', numGrayLayers); +elseif (hemiNum == 0) + [voxels,vw] = meshGetWhite(vw, 'left', numGrayLayers); + [voxels,vw] = meshGetWhite(vw, 'right', numGrayLayers,voxels); +end + +% host = 'localhost'; +% windowID = -1; + +% We build a smoothed (mesh) and an unsmoothed mesh (tenseMesh) with these calls +mrvWaitbar(0.35,wbar,sprintf('Building mesh')); +[newMesh, tenseMesh] = mrmBuild_mrMeshPy(voxels,viewGet(vw,'mmPerVox'),1); + +% Must have a name +newMesh = meshSet(newMesh,'name',meshName); +tenseMesh = meshSet(tenseMesh,'name',sprintf('%s-tense',meshName)); + +% mrvWaitbar(0.65,wbar,sprintf('meshBuild: Unsmoothed mesh vertex to gray mapping')); +initVertices = meshGet(tenseMesh,'vertices'); +newMesh = meshSet(newMesh,'initialvertices',initVertices); +vertexGrayMap = mrmMapVerticesToGray(... + initVertices, ... + viewGet(vw,'nodes'), ... + viewGet(vw,'mmPerVox'),... + viewGet(vw,'edges')); + +newMesh = meshSet(newMesh,'vertexGrayMap',vertexGrayMap); +newMesh = meshSet(newMesh,'name',meshName); +newMesh = meshSet(newMesh,'nGrayLayers',numGrayLayers); + +mrvWaitbar(0.9,wbar,sprintf('meshBuild: Saving mesh file %s',meshGet(newMesh,'name'))); + +% Save mesh file +[newMesh newMesh.path] = mrmWriteMeshFile(newMesh); + +mrvWaitbar(1,wbar,sprintf('meshBuild: Done')); +pause(0.5); +close(wbar); + +% Now refresh the UI +vw = viewSet(vw,'add and select mesh',newMesh); + +return; + +%--------------------------------------- +function classFile = verifyClassFile(vw,hemisphere) + +classFile = viewGet(vw,'classFileName',hemisphere); +str = sprintf('Class %s',classFile); + +r=questdlg(str); +if ~strcmp(r,'Yes') + switch hemisphere + case 'left' + vw = viewSet(vw,'leftClassFileName',[]); + case 'right' + vw = viewSet(vw,'rightClassFileName',[]); + end + classFile = viewGet(vw,'classFileName',hemisphere); +end + +return; + +%--------------------------------------- +function voxels = classExtractWhite(voxels,data,voi,whiteValue) +% + +% ras 05/07: the indexing of data seems off to me -- is this correct? +voxels(voi(1):voi(2), voi(3):voi(4), voi(5):voi(6)) = ... + voxels(voi(1):voi(2), voi(3):voi(4), voi(5):voi(6)) ... + | (data(voi(1):voi(2), voi(3):voi(4), voi(5):voi(6)) == whiteValue); + +return; + + +%---------------------------------------- +function [meshName,numGrayLayers,hemiNum,alpha,restrictVOI,relaxIterations] = ... + readParams(newMeshNum,hemisphere) +% +% readParams +% +% Internal routine to read the parameters for meshBuild +% +meshName = sprintf('%sSmooth',hemisphere); +numGrayLayers = 0; +switch hemisphere + case 'left' + hemiNum = 1; + case 'right' + hemiNum = 2; + case 'both' + hemiNum = 0; +end + +% transparency level (transparency is off by default, but if it gets turned +% on, this alpha parameter will have an effect). +alpha = 200; +restrictVOI = 1; +relaxIterations = 0.2; + +prompt = {'Mesh Name:',... + 'Number of Gray Layers (0-4):',... + 'Hemisphere (0=both, 1=left, 2=right):',... + % 'Default alpha (0-255):',... + % 'Inflation (0=none, 1=lots):',... + % 'Restrict to class VOI (0|1):'}; + }; +defAns = {meshName,... + num2str(numGrayLayers),... + num2str(hemiNum),... + % num2str(alpha),... + % num2str(relaxIterations),... + % num2str(restrictVOI)}; + }; + +resp = inputdlg(prompt, 'meshBuild Parameters', 1, defAns); + +if(~isempty(resp)) + meshName = resp{1}; + numGrayLayers = str2num(resp{2}); + hemiNum = str2num(resp{3}); + % alpha = str2num(resp{4}); + % relaxIterations = round(str2num(resp{5})*160); % Arbitrary choice, scales iters [0,160] + % restrictVOI = str2num(resp{6}); +else + meshName = []; + numGrayLayers = []; + hemiNum = []; + % alpha = []; + % relaxIterations = []; % Arbitrary choice, scales iters [0,160] + % restrictVOI = []; +end + +return; + + +%--------------------------------- +function [voxels,vw] = meshGetWhite(vw, hemiName, numGrayLayers, voxels) +% +% +% + +if ieNotDefined('vw'), error('You must send in a volume vw'); end +if ieNotDefined('hemiName'), error('You must define right,left or both'); end +if ieNotDefined('numGrayLayers'), numGrayLayers = 0; end + +classFile = verifyClassFile(vw,hemiName); +if isempty(classFile), + close(wbar); newMeshNum = -1; + voxels = []; + return; +end +classFileParam = [hemiName,'ClassFile']; +vw = viewSet(vw,classFileParam,classFile); + +classData = viewGet(vw,'classdata',hemiName); +if ieNotDefined('voxels'), + voxelsOld = uint8(zeros(classData.header.xsize, ... + classData.header.ysize, ... + classData.header.zsize)); +else + voxelsOld = voxels; +end +voxels = zeros(classData.header.xsize, ... + classData.header.ysize, ... + classData.header.zsize); + +% Restrict the white matter volume to a size equal to the ROI in which it +% was selected +voxels = classExtractWhite(voxels,... + classData.data,classData.header.voi,classData.type.white); +% msh = meshColor(meshSmooth(meshBuildFromClass(voxels,[1 1 1]))); +% meshVisualize(msh); + +% Add the gray matter +if(numGrayLayers>0) + + [nodes,edges,classData] = mrgGrowGray(classData,numGrayLayers); + voxels = ... + uint8( (classData.data == classData.type.white) | ... + (classData.data == classData.type.gray)); +end + +voxels = uint8(voxels | voxelsOld); + +return; diff --git a/mrMeshPy/matlabRoutines/meshColorOverlay.m b/mrMeshPy/matlabRoutines/meshColorOverlay.m new file mode 100755 index 00000000..bc84d11e --- /dev/null +++ b/mrMeshPy/matlabRoutines/meshColorOverlay.m @@ -0,0 +1,371 @@ +function [vw, dataMask, dataMaskIndices, newColors, msh, roiColors] = meshColorOverlay(vw,showData,dataOverlayScale,dataThreshold) +% Compute and possibly show color overlays in the mrMesh Window. +% +% [vw, dataMask,dataMaskIndices,newColors] = ... +% meshColorOverlay(vw,showData,dataOverlayScale,dataThreshold); +% +% Compute and show color overlays (phase, co, etc) of the mesh in the +% current 3D view (3dWindow). +% +% vw is a VOLUME{} structure containing the mesh +% showData (default = 1) indicates whether to show the new colors into +% the mrMesh window. +% roiPerimThickness (mm) +% dataOverlayScale +% dataThreshold +% +%Example: +% vw = VOLUME{1}; +% meshColorOverlay(vw); % Just renders the mesh data +% +% Returns value (and renders). Can be used to build up images. +% [junk, dm, dmIndices, newColors] = meshColorOverlay(VOLUME{1},0); +% +% Author: Ress (June, 2003) +% +% 09/2005 SOD: change to handle datasets with NaNs +% 10/2005 SOD: put in switches to handle phase data +% 2007.08.23 RFD: fixed bug that was causing the connection matrix to be +% recomputed several times in here. Should be substantially faster now. +% +% 2017.09.20 ADG: added tracking of newly computed colors as an additional +% attribute of the mesh structure for access without calling mrMesh +% +% 2010, Stanford VISTA + +if notDefined('vw'), vw = getSelectedGray; end +if notDefined('showData'), showData = 1; end +if(~exist('dataOverlayScale', 'var') ||isempty(dataOverlayScale)), dataOverlayScale = 1; end +if(~exist('dataThreshold', 'var')), dataThreshold = []; end + +msh = viewGet(vw,'currentmesh'); + +% NOTE: sometimes a view may be set to map, amp, co, or ph mode, but +% a parameter map or corAnal hasn't been loaded. (Esp. because of loading +% prefs). In this case, auto-set the view to anatomy mode, so these +% functions don't go looking for data that aren't there. (ras, 11/05) +switch vw.ui.displayMode + case 'anat' % do nothing + case 'map' + if isempty(vw.map) || isempty(vw.map{viewGet(vw, 'curScan')}) + vw.ui.displayMode = 'anat'; + end + case {'co', 'amp', 'ph'} + if isempty(vw.co) || isempty(vw.co{viewGet(vw, 'curScan')}) + vw.ui.displayMode = 'anat'; + end +end + +% Still experimenting with these settings. +prefs = mrmPreferences; +overlayModDepth = prefs.overlayModulationDepth; +dataSmoothIterations = prefs.dataSmoothIterations; + +clusterThreshold = prefs.clusterThreshold; + +% The vertex gray map can now be N x M, where N is the the number of data +% 'layers' that map to each of the M vertices. +vertexGrayMap = meshGet(msh,'vertexGrayMap'); + +if isempty(vertexGrayMap) || (strcmp(prefs.layerMapMode,'all')) + mapToAllLayers = true; + if size(vertexGrayMap,1)==1 + vertexGrayMap = mrmMapVerticesToGray(msh.initVertices, viewGet(vw,'nodes'), ... + viewGet(vw,'mmPerVox'), viewGet(vw,'edges')); + msh.vertexGrayMap = vertexGrayMap; + vw = viewSet(vw,'currentmesh',msh); + updateGlobal(vw); + end +else + mapToAllLayers = false; +end + +% Initialize the new color overlay to middle-grey. This reduces +% edge artifacts when we smooth. +sz = size(meshGet(msh,'colors'),2); + +% Center new colors around 127. Don't forget alpha channel. These are the +% vertices that are painted in the 3D view. +newColors = ones(4,sz) + 127; + +% VertInds will be a logical map that tells us which vertices have at least +% one data layer that is not null. Some vertices (find(vertInds==0)) have +% no data associated with them. Note that this mask is different from the +% dataMask, which tells us which data values are below the thresholds that +% the user has set. +vertInds = (vertexGrayMap(1,:) > 0); +if(mapToAllLayers) + for ii=2:size(vertexGrayMap,1) + vertInds = vertInds | (vertexGrayMap(ii,:) > 0); + end +end +vertInds = logical(vertInds); + +% Create the Color Overlay and Data Mask. The dataMaskIndices are the +% indices into the data (e.g., co{1}(dataMaskIndices) are the entries that +% have valid data values given the thresholds. The color overlay describes +% the RGB values for ALL of the data, not just the valid points described +% in dataMaskIndices. +%[dataMaskIndices,colorOverlay] = meshCODM(vw); +% 09/2005 SOD: added a flag (phaseFlag) which indicates whether the +% data (data) is phase (circular) or not. This will have an impact +% on the kind of smoothing that is done later on. +[dataMaskIndices, data, cmap, dataRange, phaseFlag] = meshCODM(vw,clusterThreshold); + +%%%%% Grab data for the overlay from the appropiate gray nodes %%%%% +% The new colors are assigned from the color overlay. The color +% overlay was found from the VOLUME data. The mapping from volume +% locations to vertex indices is managed by vertexGrayMap. Each vertex is +% mapped to a gray node- except the vertices that have no nearest gray +% node- they are zeroed-out in vertexGrayMap and here are flagged by +% vertInds. Remeber- the indices to vertexGrayMap are in vertex space- +% each value in vertexGrayMap corresponds to an entry in +% mesh.vertices and specifies which gray matter node is closest to that +% vertex. +% +% ras 05/06 comment: the logic here is pretty painful. Why not consolidate +% the 2 prefs dealing with mapping into a single, expandable preference, +% and make this a switch statement? We start with possible mappings +% {'layer1', 'mean', 'max'} but should allow for other ones: 'layer2', +% 'min', 'max'. Also, needs to be more modular b/c it's clear other pieces +% of code aren't doing things consistently with this (see ROI drawing +% code), and it's very hard to figure out what goes wrong when things +% break. +if(mapToAllLayers && size(vertexGrayMap,1)>1 && ~isempty(data)) + + % Initialize dataOverlay with the mean data value to avoid bias when + % smoothing. + %dataOverlay = zeros(size(vertexGrayMap(vertInds))); + allV = vertexGrayMap > 0; + + % ignore non-finite data such as NaN's + tmpdata = data(vertexGrayMap(allV(:))); + dataOverlay = repmat(mean(tmpdata(isfinite(tmpdata))),1,sz); + % dataOverlay = repmat(mean(data(vertexGrayMap(allV(:)))),1,sz); + + switch prefs.overlayLayerMapMode + case 'mean' + % take mean data across layers + + if ~phaseFlag + n = sum(allV,1); + for ii=1:size(vertexGrayMap,1) + % 03/2006 SOD: If we average over the vertexGrayMap layers we + % should start from layer 1 and not with the already existing + % dataOverlay. The existing dataOverlay is a mean of all data + % and will bias the layer-based average. + if ii==1 + dataOverlay(allV(ii,:)) = data(vertexGrayMap(ii,allV(ii,:))); + else + dataOverlay(allV(ii,:)) = dataOverlay(allV(ii,:)) + data(vertexGrayMap(ii,allV(ii,:))); + end + end + dataOverlay(vertInds) = dataOverlay(vertInds)./n(vertInds); + + else % for phase data go complex average and go back + % recompute complex dataOverlay, taking only finite data + tmpdata = -exp(1i*data(vertexGrayMap(allV(:)))); + dataOverlay = repmat(mean(tmpdata(isfinite(tmpdata))),1,sz); + %dataOverlay = repmat(mean(-exp(i*data(vertexGrayMap(allV(:))))),1,sz); + n = sum(allV,1); + for ii=1:size(vertexGrayMap,1) + if ii==1 + dataOverlay(allV(ii,:)) = -exp(1i*data(vertexGrayMap(ii,allV(ii,:)))); + else + dataOverlay(allV(ii,:)) = dataOverlay(allV(ii,:)) + -exp(1i*data(vertexGrayMap(ii,allV(ii,:)))); + end + end + dataOverlay(vertInds) = angle(dataOverlay(vertInds)./n(vertInds))+pi; + end + + case 'max' + % take max value across layers + curData = zeros(size(dataOverlay)); + for ii=1:size(vertexGrayMap,1) + curData(:) = -9999; + curData(allV(ii,:)) = data(vertexGrayMap(ii,allV(ii,:))); + biggerInds = curData>dataOverlay; + dataOverlay(biggerInds) = curData(biggerInds); + end + case 'min' + % take min value across layers + curData = zeros(size(dataOverlay)); + for ii=1:size(vertexGrayMap,1) + curData(:) = 9999; + curData(allV(ii,:)) = data(vertexGrayMap(ii,allV(ii,:))); + smallerInds = curDatadataOverlay; + dataOverlay(biggerInds) = curData(biggerInds).*signs(biggerInds); + end + end +else + % layer 1 mapping only + dataOverlay = zeros(size(vertInds)); + if ~isempty(data) + dataOverlay(vertInds) = data(vertexGrayMap(1,vertInds)); + end +end + +% The dataMaskIndices are indices into the gray matter nodes. vertexGrayMap +% maps each vertex to some gray matter nodes and vertInds just tells us which +% entries in vertexGrayMap are non-zero (zero entries mean that vertex has +% no gray node close to it). So, we use ismember to create a new data mask +% that is in vertex space rather than node space. (ie. for each non-zero +% gray-node in vertexGrayMap, is it a memeber of dataMaskIndices?) +tmp = ismember(vertexGrayMap, dataMaskIndices); +dataMask = any(tmp,1); +dataMask(~vertInds) = 0; + +if(~isempty(dataOverlayScale)) + dataOverlay = dataOverlay.*dataOverlayScale; +end +if(~isempty(dataThreshold)) + if(numel(dataThreshold)==1) + dataMask = dataMask & dataOverlay>dataThreshold; + else + dataMask = dataMask & dataOverlay>dataThreshold(1) & dataOverlay0) + for t=1:dataSmoothIterations + % Smooth and re-threshold the data mask + dataMask = connectionBasedSmooth(conMat, double(dataMask)); + dataMask = dataMask>=0.5; + % New: smooth the data, not the colors (RFD 2005.08.15) + % 09/2005 SOD: smoothing the data can only be done if the + % data is non-circular (ie not valid for phase maps). So + % two ways depending on the data: + if ~phaseFlag + dataOverlay = connectionBasedSmooth(conMat,double(dataOverlay)); + else + % phase data, so go complex, smooth and go back + dataOverlay = -exp(1i*dataOverlay); + dataOverlay = connectionBasedSmooth(conMat,double(dataOverlay)); + dataOverlay = angle(dataOverlay)+pi; + end + end +end + +% compute newColors now after data smoothing +if ~isequal(vw.ui.displayMode, 'anat') && sum(data) ~= 0 + % 09/2005 SOD: remove NaNs from vertInds. This is otherwise done by + % meshData2Colors, changing the size of the output colors, which + % will consequently result in an error due to subscripted + % assignment dimension mismatch. + vertInds(isnan(dataOverlay))=0; + newColors(1:3,vertInds) = meshData2Colors(dataOverlay(vertInds), cmap, dataRange); +end + +% Assign the anatomy colors to the locations where there are no data +% values. +oldColors = meshGet(msh,'colors'); + +% ras, 08/2007: added a flag to allow you to set the transparency +% according to the 'co' field. This way, weaker activations appear fainter. +% This seems to work nicely in Freesurfer-generated figures from papers. +% (I put this before the ROI drawing, b/c we don't want to fade the ROIs.) +if prefs.coTransparency==1 && ~isequal(vw.ui.displayMode, 'co') + co = viewGet(vw, 'scanco'); % coherence for this scan + + if ~isempty(co) + % map the co values to each mesh vertex (need a mrVista2 + % function here) + co = meshMapData(msh, co, 0, 'layer1'); + + % convert the coherence levels to a scan weight: the cothresh should + % be very transparent, but reasonably high co should be opaque. + % Let's try getting this range from the clip mode of the coherence + % display mode (i.e., vw.ui.coMode.clipMode.) Usually this is [0 1], + % but you can set it to max out at a lower value (and won't need to + % make the mrmPreferences specification more complex). + if checkfields(vw, 'ui', 'coMode', 'clipMode') + clim = vw.ui.coMode.clipMode ./ 2; + if isempty(clim) || isequal(clim, 'auto') + clim = [min(co(:)) max(co(:))]; + end + else + clim = [viewGet(vw, 'cothresh') max(co(:))/2]; + end + + % clip and rescale to [0, 1] + co(co < clim(1)) = clim(1); + co(co > clim(2)) = clim(2); + w = normalize(co, 0, 1); + + %% manually alpha the underlay colors through the overlay + + % we don't want to modify new colors outside the data mask: + w( ~dataMask ) = 1; + + % replicate w to affect the [R G B] values (but not the [alpha]) + % channel, which is not working, and which is why I need to + % manually compute the alpha layering): + w = [repmat(w(:)', [3 1]); ones(1, length(w))]; + + % main alpha layering (don't worry about dataMask, taken care of + % above) + newColors = w .* newColors + (1-w) .* double(oldColors); + end +end + + +%%%%% Modulate overlay colors to allow underlying curvature to show %%%%% +% This is useful when your mesh is completely painted, but +% you need some clue about the surface curvature (which is presumable what +% the original mesh colors represent). -Bob +if (overlayModDepth>0) +% newColors(:,dataMask) = (1-overlayModDepth)*newColors(:,dataMask) ... +% + overlayModDepth*double(oldColors(:,dataMask)); +% newColors(newColors>255) = 255; +% newColors(newColors<0) = 0; + + newColors = (1-overlayModDepth)*newColors ... + + overlayModDepth*double(oldColors); + newColors(newColors>255) = 255; + newColors(newColors<0) = 0; +end + +%%%% Draw ROIs: now in separate function %%%%% +[newColors, roiVertInds, dataMask, vw, roiColors] = meshDrawROIs(vw, newColors, vertexGrayMap, [], dataMask, prefs); + +% Apply the data mask +newColors(1:3,~dataMask) = oldColors(1:3,~dataMask); + +% Place the new colors in the rendering +newColors = uint8(round(newColors)); + +%%%% For mrMeshPy functionality, save the newly computed ("current") colors +%%%% to the mesh structure - currently removing alpha channel support TODO +% ADG YNiC 2017 +msh.currentColors = newColors; +msh.currentColors(4,:) = 255; +disp('Assigned new attribute to msh instance - currentColors - for mrMeshPy'); + +% Sometimes we just want the values. Usually, we want to show the data. +if showData + msh = mrmSet(msh,'colors',newColors'); +end + +return; diff --git a/mrMeshPy/matlabRoutines/mrMeshPyBuildMesh.m b/mrMeshPy/matlabRoutines/mrMeshPyBuildMesh.m new file mode 100755 index 00000000..01107c4f --- /dev/null +++ b/mrMeshPy/matlabRoutines/mrMeshPyBuildMesh.m @@ -0,0 +1 @@ +%% coming later \ No newline at end of file diff --git a/mrMeshPy/matlabRoutines/mrMeshPySend.m b/mrMeshPy/matlabRoutines/mrMeshPySend.m new file mode 100755 index 00000000..5c8aa86a --- /dev/null +++ b/mrMeshPy/matlabRoutines/mrMeshPySend.m @@ -0,0 +1,443 @@ +function myOutput = mrMeshPySend(the_command, the_data) +% function myOutput = mrMeshPySend(the_command, the_data) +% +% +% This module interprets the commands sent from mrVista output to +% mrMeshPy (python) through the mrMeshPyServer .. +% +% Matlab sends a string in one transaction giving a command to the +% visualisation module. This command either performs an explicit +% function in the viewer (e.g. rotate the camera 90 degrees) or the +% commnd describes the configuration/content of a large data chunk to +% will be sent in the subsequent transaction so that we know how to +% unpack the data, and how to process it (e.g. 70,000 floating point +% numbers which are scalar values to show as an amplitude map. +% +% N.B. - currently command strings have a maximum length of 1024 bytes. +% +% Commands are specifically ordered, semi-colon seperated strings which are +% unpacked to describe what the user is trying to do / send from matlab. +% Commands have a MINIMUM LENGTH of 6 arguments and have the following +% structure and item order (zero-indexed) +% +% 0 - "cmd" -- always this, identifies it as a cmd :) +% 1-3 -- place holders for later faetures +% 4 - commandName - should match a command in mp_Commands file +% 5 - theMeshInstance - integer pointing to the the mesh window that we +% want to operate on +% 6 onwards - commandArgs - a list of comma-separated pairs of arguments +% to characterise the processing of the incoming data +% blob or apply some settings to the viewport - +% CAN BE EMPTY but must be set to [] +% +% +% Andre' Gouws 2017 + +%% new mesh +if strcmp(the_command, 'sendNewMeshData') == 1 + + % get the target mesh from the VOLUME + VOLUME{1} = the_data(1); + msh = VOLUME{1}.mesh{VOLUME{1}.meshNum3d}; % expects a mesh (msh) structure + + % send the new vertices as a vector ([x1 y1 z1 x2 y2 z2 ..... xn yn zn]) + verticesVector = reshape(msh.vertices,1,size(msh.vertices,1)*size(msh.vertices,2)); + verticesLength = size(verticesVector,2); + + % send the new vertices as a vector ([t1a t1b t1c t2a t2b t2c ..... tna tnb tnc]) + trianglesVector = reshape(msh.triangles,1,size(msh.triangles,1)*size(msh.triangles,2)); + trianglesLength = size(trianglesVector,2); + + + currColors = uint8(msh.colors); + + % ---- set up the TCP command -- working example ---------------------- + + % now that we have all the data we need lets generate a TCP command + % that mrMeshPy can interpret - the server will be expecting the + % command in a specific format that a wrapper script can generate for + % us if we pass it a struct of commands/args as below + + theCommandStruct = struct('cmdIdentifier','cmd'); %always start a command with this + theCommandStruct.CommandToSend = 'loadNewMesh'; % a command listed in mrMeshPyCommandsList.m + theCommandStruct.TargetMeshSession = msh.mrMeshPyID; % function to create unique ID based on system clock + theCommandStruct.Args = {}; % we will use a cell array to pass all the extra arguments we want + + % set up the first extra arg: mrMeshPy will be expecting a data blob + % called "vertices" of double type ("d") of length "verticesLength" + % NB a command argument set is single string with commas to separate + % different components of the argument + theCommandStruct.Args{1} = ['vertices,d,',num2str(verticesLength)]; + + % 3rd set of arguments- same logic as Arg1 above + theCommandStruct.Args{end+1} = ['triangles,d,',num2str(trianglesLength)]; + + % mrMeshPy will be expecting a data blob of R values for the colour + % lookup table + theCommandStruct.Args{end+1} = ['r_rgb,B,',num2str(length(currColors(1,:)))]; + + % mrMeshPy will be expecting a data blob of G values for the colour + % lookup table + theCommandStruct.Args{end+1} = ['g_rgb,B,',num2str(length(currColors(2,:)))]; + + % mrMeshPy will be expecting a data blob of B values for the colour + % lookup table + theCommandStruct.Args{end+1} = ['b_rgb,B,',num2str(length(currColors(3,:)))]; + + % mrMeshPy will be expecting a data blob of B values for the colour + % lookup table + theCommandStruct.Args{end+1} = ['a_rgb,B,',num2str(length(currColors(4,:)))]; + + % get our wrapper to generate our command for us + initialTCPCommand = createMrMeshPySendTCPCommand(theCommandStruct) + + %open the port + t = tcpclient('localhost', 9999, 'Timeout', 20); + + % write our initial setup / staging command that explains what's coming + % in the following transactions + write(t,uint8([initialTCPCommand])); + + % write the data that will be processed based on the previous command + write(t,verticesVector); + write(t,trianglesVector); + write(t,currColors(1,:)); + write(t,currColors(2,:)); + write(t,currColors(3,:)); + write(t,currColors(4,:)); + + mesh_reply = read(t); + while isempty(mesh_reply) + mesh_reply = read(t); + disp 'waiting' + pause(0.1); + end + + disp(char(mesh_reply)); + + +%% smooth existing mesh +elseif strcmp(the_command, 'smoothMesh') == 1 + %simple smoothing routine + + %the_data + + targetMesh = the_data{1}; + iterations = the_data{2}; + relaxationfactor = the_data{3}; + + VOLUME{1} = the_data{4}; %place in scope in case we need a re-synch + + theCommandStruct = struct('cmdIdentifier','cmd'); %always start a command with this + theCommandStruct.CommandToSend = 'smoothMesh'; % a command listed in mrMeshPyCommandsList.m + theCommandStruct.TargetMeshSession = targetMesh; % BIG change here - string identifier now + theCommandStruct.Args = {}; % we will use a cell array to pass all the extra arguments we want + + theCommandStruct.Args{1} = ['iterations,',num2str(iterations),',relaxationfactor,',num2str(relaxationfactor)]; + + % get our wrapper to generate our command for us + initialTCPCommand = createMrMeshPySendTCPCommand(theCommandStruct) + + %open the port + t = tcpclient('localhost', 9999, 'Timeout', 20); + + % write our initial setup / staging command that explains what's coming + % in the following transactions + write(t,uint8([initialTCPCommand])); + + % done! - no extra data blob to send here! + + mesh_reply = read(t); + while isempty(mesh_reply) + mesh_reply = read(t); + disp 'waiting' + pause(0.1); + end + + disp(char(mesh_reply)); + + actionToTake = mrMeshReplyHandler(mesh_reply, targetMesh); + + if actionToTake == 101 %reload required + mrMeshPySend('sendNewMeshData',VOLUME{1}); + end + + %% TODO - auto-remap the colour data? + % we do however need to make sure we re-map the correct colors to their + % corresponding vertices - this gets screwed up in the smoothing + %mrMeshPySend('updateMeshData',currView); + + +%% update the currently displayed data on the mesh to reflect the VOLUME +% view +elseif strcmp(the_command, 'updateMeshData') == 1 + % TODO - description + + currView = the_data; % expects a View structure + + try + newColors = uint8(currView.mesh{currView.meshNum3d}.currentColors); + catch + disp 'here' + pause; + newColors = uint8(currView.mesh{currView.meshNum3d}.colors); + end + + currScanNum = currView.curScan; % TODO needed? + currOverlayType = currView.ui.displayMode; % TODO needed? + + theCommandStruct = struct('cmdIdentifier','cmd'); %always start a command with this + theCommandStruct.CommandToSend = 'updateMeshData'; % a command listed in mrMeshPyCommandsList.m + currMesh = currView.meshNum3d; % keep track of current / apt / target mesh + theCommandStruct.TargetMeshSession = currView.mesh{currMesh}.mrMeshPyID; % unique string id for this mesh + theCommandStruct.Args = {}; % we will use a cell array to pass all the extra arguments we want + + % mrMeshPy will be expecting a data blob of R values for the colour + % lookup table + theCommandStruct.Args{end+1} = ['r_rgb,B,',num2str(length(newColors(1,:)))]; + + % mrMeshPy will be expecting a data blob of G values for the colour + % lookup table + theCommandStruct.Args{end+1} = ['g_rgb,B,',num2str(length(newColors(2,:)))]; + + % mrMeshPy will be expecting a data blob of B values for the colour + % lookup table + theCommandStruct.Args{end+1} = ['b_rgb,B,',num2str(length(newColors(3,:)))]; + + % mrMeshPy will be expecting a data blob of B values for the colour + % lookup table + theCommandStruct.Args{end+1} = ['a_rgb,B,',num2str(length(newColors(4,:)))]; + + + % get our wrapper to generate our command for us + initialTCPCommand = createMrMeshPySendTCPCommand(theCommandStruct) + + %open the port + t = tcpclient('localhost', 9999, 'Timeout', 20); + + % write our initial setup / staging command that explains what's coming + % in the following transactions + write(t,uint8([initialTCPCommand])); + + % write the data that will be processed based on the previous command + write(t,newColors(1,:)); + write(t,newColors(2,:)); + write(t,newColors(3,:)); + write(t,newColors(4,:)); + + mesh_reply = read(t); + while isempty(mesh_reply) + mesh_reply = read(t); + disp 'waiting' + pause(0.1); + end + + actionToTake = mrMeshReplyHandler(mesh_reply, currView.mesh{currMesh}.mrMeshPyID); + + if actionToTake == 101 %reload required + mrMeshPySend('sendNewMeshData',currView); + end + + +%% get an roi drawn on the mesh back into matlab +elseif strcmp(the_command, 'checkMeshROI') == 1 + % extract a set of ROI vertices from mrMeshPy as long as one is ready + + currView = the_data % expects a View structure + currMesh = currView.meshNum3d; % keep track of current / apt / target mesh + + theCommandStruct = struct('cmdIdentifier','cmd'); %always start a command with this + theCommandStruct.CommandToSend = 'checkMeshROI'; % a command listed in mrMeshPyCommandsList.m + theCommandStruct.TargetMeshSession = currView.mesh{currMesh}.mrMeshPyID; % unique string pointing to vtk session + theCommandStruct.Args = {}; % we will use a cell array to pass all the extra arguments we want + + % get our wrapper to generate our command for us + initialTCPCommand = createMrMeshPySendTCPCommand(theCommandStruct); + + %open the port + t = tcpclient('localhost', 9999, 'Timeout', 10); + + % write our initial setup / staging command that explains what's coming + % in the following transactions + write(t,uint8([initialTCPCommand])); + + % wait til the server returns something before moving on + mesh_reply = read(t); + while isempty(mesh_reply) + mesh_reply = read(t); + disp 'waiting' + pause(0.1); + end + + disp(char(mesh_reply)); + + + % TODO add timeout? + + disp(['Got response - ', mesh_reply]); + assignin('base','reply',mesh_reply); + + % unpack the + if strcmp(char(mesh_reply(1:8)),'RoiReady') + replyStr = char(mesh_reply); + replyArgs = strsplit(replyStr,',') + else + disp('Error: No ROI data received? ...') + return + end + + % get ready to receive data + incomingDataType = replyArgs{2}; + incomingDataCount = str2num(replyArgs{3}); + + if strcmp(incomingDataType, 'double') + expectedBytes = incomingDataCount*8 % 8 bytes per double + else + disp('Error: No ROI data received? ...') + return + end + pause(1.0); + + theCommandStruct = struct('cmdIdentifier','cmd'); %always start a command with this + theCommandStruct.CommandToSend = 'sendROIVertices'; % a command listed in mrMeshPyCommandsList.m + theCommandStruct.TargetMeshSession = currView.mesh{currMesh}.mrMeshPyID; % unique string pointing to vtk session + theCommandStruct.Args = {}; % we will use a cell array to pass all the extra arguments we want% tell meshPy we're ready for the vertices + + % get our wrapper to generate our command for us + TCPCommand = createMrMeshPySendTCPCommand(theCommandStruct); + + %open the port + t = tcpclient('localhost', 9999, 'Timeout', 100); + + % write our initial setup / staging command that explains what's coming + % in the following transactions + write(t,uint8([TCPCommand])); + + + reply = read(t, expectedBytes) + while isempty(reply) + reply = read(t, expectedBytes); + disp 'waiting' + pause(0.1); + end + + %debug + assignin('base','reply',reply); + + vertices = typecast(uint8(reply), 'double'); + vertices = uint32(vertices) + 1; + + + + %% TODO layer 1 only code? +% % % verts = adjustPerimeter(vertices, [], currView)'; +% % % msh = currView.mesh{currMesh} +% % % +% % % vert = meshGet(msh,'initialVertices'); +% % % coords = vert([2 1 3],vertices); +% % % for ii=1:3, +% % % coords(ii,:) = coords(ii,:)./msh.mmPerVox(ii); +% % % end +% % % coords = round(coords); + + + %% currently support "all layers" mode only + verts = adjustPerimeter(vertices, [], currView)'; + + msh = currView.mesh{currMesh} + grayInds = msh.vertexGrayMap(1,verts); + curLayer = unique(grayInds(grayInds>0)); + allLayers = curLayer; + nodes = viewGet(currView, 'nodes'); + edges = viewGet(currView, 'edges'); + curLayerNum = 1; + while(~isempty(curLayer)) + nextLayer = []; + curLayerNum = curLayerNum+1; + for ii=1:length(curLayer) + offset = nodes(5,curLayer(ii)); + if offset>length(edges), continue; end + numConnected = nodes(4,curLayer(ii)); + neighbors = edges(offset:offset+numConnected-1); + nextLayer = [nextLayer, neighbors(nodes(6,neighbors)==curLayerNum)]; + end + nextLayer = unique(nextLayer); + allLayers = [allLayers, int32(nextLayer)]; + curLayer = nextLayer; + end + coords = currView.coords(:,allLayers); + + % user dialog to name the incoming ROI + + prompt={'Enter a name for the ROI:'}; + name='Incoming ROI name'; + numlines=[1,100]; + defaultanswer={'ROI_new'}; + options.Resize='on'; + options.WindowStyle='normal'; + options.Interpreter='tex'; + + roiName = inputdlg(prompt,name,numlines,defaultanswer); + if isempty(roiName) + roiName = 'ROI_new'; + else + roiName = roiName{1}; + end + + currView = newROI(currView,roiName,1,[],coords); + [currView,~,~,~,currView.mesh{currView.meshNum3d}] = meshColorOverlay(currView,0); + mrMeshPySend('updateMeshData',currView); + + % return + myOutput = currView; + + +%% a simple animation test +elseif strcmp(the_command, 'rotateMeshAnimation') == 1 + + + %data - rotate 90 steps, 4 degrees at a time + rotations = ones(1,90.0)*4; + + % ---- set up the TCP command + + theCommandStruct = struct('cmdIdentifier','cmd'); %always start a command with this + theCommandStruct.CommandToSend = 'rotateMeshAnimation'; % a command listed in mrMeshPyCommandsList.m + theCommandStruct.TargetMeshSession = the_data - 1; % integer pointing to target mrMeshypy sub-window, minus one for zero-index in python + theCommandStruct.Args = {}; % we will use a cell array to pass all the extra arguments we want + + % set up the first extra arg: mrMeshPy will be expecting a data blob + % called "vertices" of double type ("d") of length "verticesLength" + % NB a command argument set is single string with commas to separate + % different components of the argument + theCommandStruct.Args{1} = ['rotate,d,',num2str(length(rotations))]; + + % get our wrapper to generate our command for us + initialTCPCommand = createMrMeshPySendTCPCommand(theCommandStruct); + + %open the port + t = tcpclient('localhost', 9999, 'Timeout', 20); + + % write our initial setup / staging command that explains what's coming + % in the following transactions + write(t,uint8([initialTCPCommand])); + + % write the data that will be processed based on the previous command + write(t,rotations); + + %% wait til the server returns something before moving on - TODO - get some useful info back from server not just junk + y = read(t); + while isempty(y) + y = read(t); + disp 'waiting' + pause(0.1); + end + y + +%% or else the command is not recognised +else + disp 'ERROR: Matlab generated a command that is not recognised - not sending!!' + +end + + diff --git a/mrMeshPy/matlabRoutines/mrMeshReplyHandler.m b/mrMeshPy/matlabRoutines/mrMeshReplyHandler.m new file mode 100755 index 00000000..c7a7eba4 --- /dev/null +++ b/mrMeshPy/matlabRoutines/mrMeshReplyHandler.m @@ -0,0 +1,35 @@ +function replyAction = mrMeshReplyHandler(mesh_reply, mrMeshPyID) +%function replyAction = mrMeshReplyHandler(mesh_reply, mrMeshPyID) +% +% This function processes replies from the mrMeshPy TCP server and takes +% action according to the response received + +disp('running reply handler'); + +if strcmp(char(mesh_reply), 'Mesh smooth failed') + choice = questdlg(['No mesh with ID ',mrMeshPyID,'? - shall I try reloading it?'], ... + 'Missing Mesh ', ... + 'Yes - try reload','No - Cancel','No - Cancel'); + % Handle response + switch choice + case 'Yes - try reload' + replyAction = 101; %TODO 101 indicates a reload?? + case 'No - Cancel' + replyAction = 0; + end + +elseif strcmp(char(mesh_reply), 'Mesh update failed') + choice = questdlg(['No mesh with ID ',mrMeshPyID,'? - shall I try reloading it?'], ... + 'Missing Mesh ', ... + 'Yes - try reload','No - Cancel','No - Cancel'); + % Handle response + switch choice + case 'Yes - try reload' + replyAction = 101; %TODO 101 indicates a reload?? + case 'No - Cancel' + replyAction = 0; + end +else + disp('No handler for this reply yet.. continuing.'); + replyAction = 0; +end \ No newline at end of file diff --git a/mrMeshPy/matlabRoutines/mrmBuild_mrMeshPy.m b/mrMeshPy/matlabRoutines/mrmBuild_mrMeshPy.m new file mode 100755 index 00000000..dd7e3c23 --- /dev/null +++ b/mrMeshPy/matlabRoutines/mrmBuild_mrMeshPy.m @@ -0,0 +1,46 @@ +function [msh, unsmoothedMsh] = mrmBuild_mrMeshPy(wm,mmPerVox,smoothFlag,visualizeFlag,varargin) +%Build a smoothed and unsmooth mesh from a white matter voxels +% +% [msh, lights, unsmoothedMsh] = mrmBuild(wm,mmPerVox,smoothFlag,visualizeFlag,varargin) +% +% If you do not want the smoothed mesh, set smoothFlag = 0 +% +% This routine will replace mrmBuildMesh shortly. +% + +if ieNotDefined('wm'), error('White matter voxels required.'); end +if ieNotDefined('mmPerVox'), mmPerVox = [1 1 1]; end +if ieNotDefined('smoothFlag'), smoothFlag = 1; end +%if ieNotDefined('visualizeFlag'), visualizeFlag = 1; end %AG edit - TODO +if ieNotDefined('visualizeFlag'), visualizeFlag = 0; end + +lights = []; + +nVar = length(varargin); +if mod(nVar,2), error('Parameters must be (name,value)'); end + +% Build the one that is only marching cubes +unsmoothedMsh = meshBuildFromClass_mrMeshPy(wm,mmPerVox); + +% We should think about moving the proper format of msh into buildMesh (the +% mex file). This requires recompilation, we think, so for now we are +% handling this issue outside of the mex files, in Matlab. +unsmoothedMsh = meshFormat(unsmoothedMsh); +% meshVisualize(unsmoothedMsh); + +msh = unsmoothedMsh; + +% % % % % % Let VTK smooth the mesh +% % % % % if smoothFlag, +% % % % % for ii=1:2:nVar +% % % % % msh = meshSet(msh,varargin{ii},str2num(varargin{ii+1})); +% % % % % end +% % % % % msh = meshSmooth(unsmoothedMsh); +% % % % % msh = meshColor(msh); +% % % % % end + +% Return the user the lights +if visualizeFlag, [msh,lights] = meshVisualize(msh); end + +return; + diff --git a/mrMeshPy/matlabRoutines/prepareMrMeshPyColourMap.m b/mrMeshPy/matlabRoutines/prepareMrMeshPyColourMap.m new file mode 100755 index 00000000..03d70307 --- /dev/null +++ b/mrMeshPy/matlabRoutines/prepareMrMeshPyColourMap.m @@ -0,0 +1,26 @@ +function rgb_matrix = prepareMrMeshPyColourMap(cmap) +% function rgb_matrix = prepareMrMeshPyColourMap(cmap) +% so mrVista uses the convetion of color mapping of scalars to lookup +% tables with 256 entries - the first 128 are used for grayscale mapping of +% the anatomy - light gray for gyri and dark gray for sulci; the 2nd half +% of the table (entries 129-256) are r/g/b/ moldulated channels that the +% scalar data are mapped to relative to a min/max range and some clipping / +% windowing. We'll rescale these 128 colour values to 1:1024 range and then +% send to mrMeshPy + +% extract just the colour section of the map +cols = cmap(129:end,:); + +% upsample the 128 samples to 1024 samples for each colour vector +r_vec = interp(cols(:,1),8); +g_vec = interp(cols(:,2),8); +b_vec = interp(cols(:,3),8); + +%put it all togtether and clean up +rgb_matrix = [r_vec,g_vec,b_vec]; + +% interpolation causes a slight under/over-shoot so clip to 0-1 range. +rgb_matrix(rgb_matrix<0) = 0; +rgb_matrix(rgb_matrix>1) = 1; + + diff --git a/mrMeshPy/matlabRoutines/pyMeshBuild.py b/mrMeshPy/matlabRoutines/pyMeshBuild.py new file mode 100755 index 00000000..0b93a2c5 --- /dev/null +++ b/mrMeshPy/matlabRoutines/pyMeshBuild.py @@ -0,0 +1,340 @@ +#!/usr/bin/python + +# TODO - python path - assuming condo here + +## HEALTH WARNING - BETA CODE IN DEVELOPMENT ## + +''' +This standalone application will build a mesh from a nifti classification file. +To keep the procedure as similar as possible to the way mrMesh used to do this, +we will keep this as a standalon application. Matlab reads in the segmented +nifti file using vistasofts own nifti class handler meshBuild>mrmBuild> +meshBuildFromClass - we just dont use the old build_mesh mex file - we do that +bit and any smoothing in this application and send a mesh struture back to +matlab. + +AG 2017 +''' + +import os,sys +import scipy +import vtk + +from numpy import * +from scipy.io import loadmat, savemat +from vtk.util import numpy_support + +debug = False + +#TODO error handling +fileToLoad = sys.argv[1] +fileToSave = sys.argv[2] + +# load the voxel data that has been dumped to disk +voxels = scipy.io.loadmat(fileToLoad) +mmPerVox = voxels['mmPerVox'][0] +if debug: print mmPerVox + +voxels = voxels['voxels'] #unpack + + +if debug: print voxels + +if debug: print shape(voxels) + +extent = shape(voxels) +if debug: print extent +if debug: print extent[0] +if debug: print extent[1] +if debug: print extent[2] + + +''' +### ------------------------------------------------------------------------------ +### this is faster but for now exactly replicate the way mrMesh sets up the volume array +# import voxels to vtk +dataImporter = vtk.vtkImageImport() +data_string = voxels.tostring() +dataImporter.CopyImportVoidPointer(data_string, len(data_string)) +dataImporter.SetDataScalarTypeToUnsignedChar() +dataImporter.SetDataExtent(0, extent[2]-1, 0, extent[1]-1, 0, extent[0]-1) # TODO have to work this out +dataImporter.SetWholeExtent(0, extent[2]-1, 0, extent[1]-1, 0, extent[0]-1) # TODO have to work this out +dataImporter.SetDataSpacing(mmPerVox[0],mmPerVox[1],mmPerVox[2]) # TODO have to work this out +dataImporter.Update() +if debug: print dataImporter.GetOutput() + +pClassData = dataImporter.GetOutput() +pClassData.Modified() + +### ------------------------------------------------------------------------------ +''' + + +### ------- the way mrMesh did it in mesh_build -------------------------------- +pArray = map(ord,voxels.tostring()) #unpack +pDims = shape(voxels) +if debug: print pDims + +# need a reshuffle of indices ?? -- doesn't matter if cube (same xyz dims) +pDims = [pDims[2],pDims[1],pDims[0]] + +scale = mmPerVox +iSizes = [pDims[0]+2, pDims[1]+2, pDims[2]+2] + +nTotalValues = iSizes[0] * iSizes[1] * iSizes[2] + +pClassValues = vtk.vtkUnsignedCharArray() +pClassValues.SetNumberOfValues(nTotalValues) + +pClassData = vtk.vtkStructuredPoints() +pClassData.SetDimensions(iSizes[0], iSizes[1], iSizes[2]) +pClassData.SetOrigin(-scale[0], -scale[1], -scale[2]) #??? +pClassData.SetOrigin(-1, -1, -1) #??? +pClassData.SetSpacing(scale[0], scale[1], scale[2]) + + +for iSrcZ in range(pDims[2]): + + for iSrcY in range(pDims[1]): + + iSrcIndex = iSrcZ * pDims[1] * pDims[0] + iSrcY * pDims[0] + iDstIndex = (iSrcZ+1) * iSizes[1] * iSizes[0] + (iSrcY+1) * iSizes[0] + 1 + + for iSrcX in range(pDims[0]): + + fTemp = int(pArray[iSrcIndex]) + #print fTemp, 'iSrcIndex', iSrcIndex, 'iDstIndex', iDstIndex + if fTemp>0: + pClassValues.SetValue(iDstIndex, 0) + else: + pClassValues.SetValue(iDstIndex, 1) + + + iSrcIndex+=1 + iDstIndex+=1 + + +pClassData.GetPointData().SetScalars(pClassValues) + +pClassData.Modified() + + + + + +if debug: + spw = vtk.vtkStructuredPointsWriter() + spw.SetFileTypeToASCII() + spw.SetInputData(pClassData) + spw.SetFileName("/tmp/test-mrMeshPy-structuredPoints.vtk") + spw.Write() + spw.Update() + + + + +### ------ Data volume is loaded and constructed - extract some surgfaces ------------- + +mc = vtk.vtkMarchingCubes() +# mc = vtk.vtkContourFilter() #- could use a contour filter instead? + +# mc.SetInputConnection(dataImporter.GetOutputPort()) # later - for use with direct imagedata import + +mc.SetInputData(pClassData) +mc.SetValue(0,0.5) #extract 0-th surface at 0.5? +mc.ComputeGradientsOff() +mc.ComputeNormalsOff() +mc.ComputeScalarsOff() +mc.Update() + + +if debug: + print mc.GetOutput() + write = vtk.vtkPolyDataWriter() + write.SetFileName('/htmp/test-mrMeshPy-marchingCubesOutput.txt') + write.SetFileTypeToASCII() + write.SetInputData(mc.GetOutput()) + write.Write() + write.Update() + + +# ---- extract center surface - edges are normally extracted to (the cube around the edge of the volume)-------- + +confilter = vtk.vtkPolyDataConnectivityFilter() +confilter.SetInputConnection(mc.GetOutputPort()) +confilter.SetExtractionModeToClosestPointRegion() +confilter.SetClosestPoint(extent[0]/2.0,extent[1]/2.0,extent[2]/2.0) # center of volume +confilter.Update() + + + +# ---- Normals --------------------- + +# normals already computed by mc algorithm so this code is obsolete +normals = vtk.vtkPolyDataNormals() +normals.ComputePointNormalsOn() +normals.SplittingOff() +normals.SetInputConnection(confilter.GetOutputPort()) +#normals.SetInputData(discrete.GetOutput()) +normals.Update() +print normals.GetOutput() + +norm = normals.GetOutput().GetPointData().GetNormals() +output_normals = array(numpy_support.vtk_to_numpy(norm).transpose(),'d') + +####if debug: print output_normals + + + +# ---- Initial vertices - unsmoothed --------------------- +init_verts = normals.GetOutput().GetPoints().GetData() +output_init_verts = array(numpy_support.vtk_to_numpy(init_verts).transpose(),'d') + +if debug: print output_init_verts + + +# ---- Polys (triangles) --------------------- +triangles = normals.GetOutput().GetPolys().GetData() +tmp_triangles = numpy_support.vtk_to_numpy(triangles) + +# N.B. the polygon data returned here have 4 values for poly - the first is the number +# of vertices that describe the polygo (ironically always 3) and the next 3 are the +# indices of the vertices that make up the polygon + +# so first we need to reshape data from a vector +tmp_triangles = reshape(tmp_triangles,(len(tmp_triangles)/4,4)) + +# and then we drop the first column (all 3's) +output_triangles = array((tmp_triangles[:,1:4]).transpose(),'d') #remember zero index here, add one for matlab + +if debug: print output_triangles + + + + + +# -------- smoothed version of mesh ---------------- + +smooth = vtk.vtkSmoothPolyDataFilter() +smooth.SetNumberOfIterations(32) #standard value sused in old mrMesh +smooth.SetRelaxationFactor(0.5) #standard value sused in old mrMesh +smooth.FeatureEdgeSmoothingOff() +smooth.SetFeatureAngle(45) +smooth.SetEdgeAngle(15) +smooth.SetBoundarySmoothing(1) +smooth.SetInputConnection(normals.GetOutputPort()) +smooth.Update() + +# different smoothing option? +''' +smooth = vtk.vtkWindowedSincPolyDataFilter() +smooth.SetInputConnection(mc.GetOutputPort()) +smooth.SetNumberOfIterations(30) +smooth.SetPassBand(0.5) +smooth.SetFeatureAngle(45) +smooth.SetEdgeAngle(15) +smooth.SetBoundarySmoothing(1) +smooth.SetFeatureEdgeSmoothing(0) +smooth.Update() +''' + + +# ---- Vertices - smoothed --------------------- +smooth_verts = smooth.GetOutput().GetPoints().GetData() +output_smooth_verts = array(numpy_support.vtk_to_numpy(smooth_verts).transpose(),'d') + +if debug: print output_smooth_verts + + +# ---- Curvature --------------------- +curvature = vtk.vtkCurvatures() +curvature.SetInputConnection(smooth.GetOutputPort()) +curvature.SetCurvatureTypeToMean() +curvature.Update() + +curv = curvature.GetOutput().GetPointData().GetScalars() +output_curvature = array(numpy_support.vtk_to_numpy(curv).transpose(),'d') + +if debug: print min(output_curvature) +if debug: print max(output_curvature) +if debug: print output_curvature + + + +# -------- colours based on curvature ------------ + +# turn curvature into color +tmp_colors = output_curvature.copy() + +#min_curv = min(tmp_colors) +#max_curv = max(tmp_colors) +#tmp_colors = (tmp_colors -min_curv) / (max_curv-min_curv) *255 + +tmp_colors[tmp_colors>=0] = 85 #standard value sused in old mrMesh +tmp_colors[tmp_colors<0] = 160 #standard value sused in old mrMesh + +output_colors = vstack((tmp_colors, tmp_colors, tmp_colors, ones((1,len(tmp_colors)))*255)) +output_colors = array(output_colors,'d') + +if debug: print output_colors + +# OK we have all the data we need now, lets write it out to file + +data = {} #empty dictionary +data['initVertices'] = output_init_verts +data['initialvertices'] = output_init_verts +data['vertices'] = output_smooth_verts +data['colors'] = output_colors +data['normals'] = output_normals +data['triangles'] = output_triangles +data['curvature'] = output_curvature + +# save it out +savemat(fileToSave,data) + + + +# data have been sent, but let's view them here + +pdm = vtk.vtkPolyDataMapper() +pdm.SetInputConnection(confilter.GetOutputPort()) + +actor = vtk.vtkActor() +actor.SetMapper(pdm) + +ren = vtk.vtkRenderer() +renWin = vtk.vtkRenderWindow() +renWin.AddRenderer(ren) + +iren = vtk.vtkRenderWindowInteractor() +iren.SetRenderWindow(renWin) + +ren.AddActor(actor) +ren.SetBackground(1,1,1) +renWin.SetSize(500,500) + +iren.Initialize() +iren.Start() + + + +pdm = vtk.vtkPolyDataMapper() +pdm.SetInputConnection(curvature.GetOutputPort()) + +actor = vtk.vtkActor() +actor.SetMapper(pdm) + +ren = vtk.vtkRenderer() +renWin = vtk.vtkRenderWindow() +renWin.AddRenderer(ren) + +iren = vtk.vtkRenderWindowInteractor() +iren.SetRenderWindow(renWin) + +ren.AddActor(actor) +ren.SetBackground(1,1,1) +renWin.SetSize(500,500) + +iren.Initialize() +iren.Start() + diff --git a/mrMeshPy/matlabRoutines/pyMeshBuild_ImageData.py b/mrMeshPy/matlabRoutines/pyMeshBuild_ImageData.py new file mode 100755 index 00000000..1133f020 --- /dev/null +++ b/mrMeshPy/matlabRoutines/pyMeshBuild_ImageData.py @@ -0,0 +1,326 @@ +#!/usr/bin/python + +# TODO - python path - assuming condo here + +## HEALTH WARNING - BETA CODE IN DEVELOPMENT ## + +''' +This standalone application will build a mesh from a nifti classification file. +To keep the procedure as similar as possible to the way mrMesh used to do this, +we will keep this as a standalon application. Matlab reads in the segmented +nifti file using vistasofts own nifti class handler meshBuild>mrmBuild> +meshBuildFromClass - we just dont use the old build_mesh mex file - we do that +bit and any smoothing in this application and send a mesh struture back to +matlab. + +AG 2017 +''' + +import os,sys +import scipy +import vtk + +from numpy import * +from scipy.io import loadmat, savemat +from vtk.util import numpy_support + +debug = False + +#TODO error handling +fileToLoad = sys.argv[1] +fileToSave = sys.argv[2] + +# load the voxel data that has been dumped to disk +voxels = scipy.io.loadmat(fileToLoad) +mmPerVox = voxels['mmPerVox'][0] +if debug: print mmPerVox + +voxels = voxels['voxels'] #unpack + + +if debug: print voxels + +if debug: print shape(voxels) + +extent = shape(voxels) +if debug: print extent +if debug: print extent[0] +if debug: print extent[1] +if debug: print extent[2] + + + +### ------------------------------------------------------------------------------ +### this is faster but for now exactly replicate the way mrMesh sets up the volume array +# import voxels to vtk +dataImporter = vtk.vtkImageImport() +data_string = voxels.tostring() +dataImporter.CopyImportVoidPointer(data_string, len(data_string)) +dataImporter.SetDataScalarTypeToUnsignedChar() +dataImporter.SetDataExtent(0, extent[2]-1, 0, extent[1]-1, 0, extent[0]-1) # TODO have to work this out +dataImporter.SetWholeExtent(0, extent[2]-1, 0, extent[1]-1, 0, extent[0]-1) # TODO have to work this out +dataImporter.SetDataSpacing(mmPerVox[0],mmPerVox[1],mmPerVox[2]) # TODO have to work this out +dataImporter.Update() +if debug: print dataImporter.GetOutput() + +### ------------------------------------------------------------------------------ +''' + + +### ------- the way mrMesh did it in mesh_build -------------------------------- +pArray = map(ord,voxels.tostring()) #unpack + +pDims = shape(voxels) +scale = mmPerVox +iSizes = [pDims[0]+2, pDims[1]+2, pDims[2]+2] + +nTotalValues = iSizes[0] * iSizes[1] * iSizes[2] + +pClassValues = vtk.vtkUnsignedCharArray() +pClassValues.SetNumberOfValues(nTotalValues) + +pClassData = vtk.vtkStructuredPoints() +pClassData.SetDimensions(iSizes[0], iSizes[1], iSizes[2]) +pClassData.SetOrigin(-scale[0], -scale[1], -scale[2]) #??? +pClassData.SetOrigin(-1, -1, -1) #??? +pClassData.SetSpacing(scale[0], scale[1], scale[2]) + + +for iSrcZ in range(pDims[2]): + + for iSrcY in range(pDims[1]): + + iSrcIndex = iSrcZ * pDims[1] * pDims[0] + iSrcY * pDims[0] + iDstIndex = (iSrcZ+1) * iSizes[1] * iSizes[0] + (iSrcY+1) * iSizes[0] + 1 + + for iSrcX in range(pDims[0]): + + fTemp = int(pArray[iSrcIndex]) + #if debug: print fTemp, 'iSrcIndex', iSrcIndex, 'iDstIndex', iDstIndex + if fTemp>0: + pClassValues.SetValue(iDstIndex, 0) + else: + pClassValues.SetValue(iDstIndex, 1) + + + iSrcIndex+=1 + iDstIndex+=1 + + +pClassData.GetPointData().SetScalars(pClassValues) + +pClassData.Modified() + +if debug: + spw = vtk.vtkStructuredPointsWriter() + spw.SetFileTypeToASCII() + spw.SetInputData(pClassData) + spw.SetFileName("/tmp/test-mrMeshPy-structuredPoints.vtk") + spw.Write() + spw.Update() + + +''' + +### ------ Data volume is loaded and constructed - extract some surgfaces ------------- + +mc = vtk.vtkMarchingCubes() # mc = vtk.vtkContourFilter() #- could use a contour filter instead? + +mc.SetInputConnection(dataImporter.GetOutputPort()) # later - for use with direct imagedata import #mc.SetInputData(pClassData) +mc.SetValue(0,0.5) #extract 0-th surface at 0.5? +mc.ComputeGradientsOff() +mc.ComputeNormalsOff() +mc.ComputeScalarsOff() +mc.Update() + + +if debug: + print mc.GetOutput() + write = vtk.vtkPolyDataWriter() + write.SetFileName('/htmp/test-mrMeshPy-marchingCubesOutput.txt') + write.SetFileTypeToASCII() + write.SetInputData(mc.GetOutput()) + write.Write() + write.Update() + + +# ---- extract center surface - edges are normally extracted to (the cube around the edge of the volume)-------- + +confilter = vtk.vtkPolyDataConnectivityFilter() +confilter.SetInputConnection(mc.GetOutputPort()) +confilter.SetExtractionModeToClosestPointRegion() +confilter.SetClosestPoint(extent[0]/2.0,extent[1]/2.0,extent[2]/2.0) # center of volume +confilter.Update() + + + +# ---- Normals --------------------- + +# normals already computed by mc algorithm so this code is obsolete +normals = vtk.vtkPolyDataNormals() +normals.ComputePointNormalsOn() +normals.SplittingOff() +normals.SetInputConnection(confilter.GetOutputPort()) +#normals.SetInputData(discrete.GetOutput()) +normals.Update() +print normals.GetOutput() + +norm = normals.GetOutput().GetPointData().GetNormals() +output_normals = array(numpy_support.vtk_to_numpy(norm).transpose(),'d') + +####if debug: print output_normals + + + +# ---- Initial vertices - unsmoothed --------------------- +init_verts = normals.GetOutput().GetPoints().GetData() +output_init_verts = array(numpy_support.vtk_to_numpy(init_verts).transpose(),'d') + +if debug: print output_init_verts + + +# ---- Polys (triangles) --------------------- +triangles = normals.GetOutput().GetPolys().GetData() +tmp_triangles = numpy_support.vtk_to_numpy(triangles) + +# N.B. the polygon data returned here have 4 values for poly - the first is the number +# of vertices that describe the polygo (ironically always 3) and the next 3 are the +# indices of the vertices that make up the polygon + +# so first we need to reshape data from a vector +tmp_triangles = reshape(tmp_triangles,(len(tmp_triangles)/4,4)) + +# and then we drop the first column (all 3's) +output_triangles = array((tmp_triangles[:,1:4]).transpose(),'d') #remember zero index here, add one for matlab + +if debug: print output_triangles + + + + + +# -------- smoothed version of mesh ---------------- + +smooth = vtk.vtkSmoothPolyDataFilter() +smooth.SetNumberOfIterations(32) #standard value sused in old mrMesh +smooth.SetRelaxationFactor(0.5) #standard value sused in old mrMesh +smooth.FeatureEdgeSmoothingOff() +smooth.SetFeatureAngle(45) +smooth.SetEdgeAngle(15) +smooth.SetBoundarySmoothing(1) +smooth.SetInputConnection(normals.GetOutputPort()) +smooth.Update() + +# different smoothing option? +''' +smooth = vtk.vtkWindowedSincPolyDataFilter() +smooth.SetInputConnection(mc.GetOutputPort()) +smooth.SetNumberOfIterations(30) +smooth.SetPassBand(0.5) +smooth.SetFeatureAngle(45) +smooth.SetEdgeAngle(15) +smooth.SetBoundarySmoothing(1) +smooth.SetFeatureEdgeSmoothing(0) +smooth.Update() +''' + + +# ---- Vertices - smoothed --------------------- +smooth_verts = smooth.GetOutput().GetPoints().GetData() +output_smooth_verts = array(numpy_support.vtk_to_numpy(smooth_verts).transpose(),'d') + +if debug: print output_smooth_verts + + +# ---- Curvature --------------------- +curvature = vtk.vtkCurvatures() +curvature.SetInputConnection(smooth.GetOutputPort()) +curvature.SetCurvatureTypeToMean() +curvature.Update() + +curv = curvature.GetOutput().GetPointData().GetScalars() +output_curvature = array(numpy_support.vtk_to_numpy(curv).transpose(),'d') + +if debug: print min(output_curvature) +if debug: print max(output_curvature) +if debug: print output_curvature + + + +# -------- colours based on curvature ------------ + +# turn curvature into color +tmp_colors = output_curvature.copy() + +#min_curv = min(tmp_colors) +#max_curv = max(tmp_colors) +#tmp_colors = (tmp_colors -min_curv) / (max_curv-min_curv) *255 + +tmp_colors[tmp_colors>=0] = 160 #standard value sused in old mrMesh +tmp_colors[tmp_colors<0] = 85 #standard value sused in old mrMesh + +output_colors = vstack((tmp_colors, tmp_colors, tmp_colors, ones((1,len(tmp_colors)))*255)) +output_colors = array(output_colors,'d') + +if debug: print output_colors + +# OK we have all the data we need now, lets write it out to file + +data = {} #empty dictionary +data['initVertices'] = output_init_verts +data['initialvertices'] = output_init_verts +data['vertices'] = output_smooth_verts +data['colors'] = output_colors +data['normals'] = output_normals +data['triangles'] = output_triangles +data['curvature'] = output_curvature + +# save it out +savemat(fileToSave,data) + + + +# data have been sent, but let's view them here + +pdm = vtk.vtkPolyDataMapper() +pdm.SetInputConnection(confilter.GetOutputPort()) + +actor = vtk.vtkActor() +actor.SetMapper(pdm) + +ren = vtk.vtkRenderer() +renWin = vtk.vtkRenderWindow() +renWin.AddRenderer(ren) + +iren = vtk.vtkRenderWindowInteractor() +iren.SetRenderWindow(renWin) + +ren.AddActor(actor) +ren.SetBackground(1,1,1) +renWin.SetSize(500,500) + +iren.Initialize() +iren.Start() + + + +pdm = vtk.vtkPolyDataMapper() +pdm.SetInputConnection(curvature.GetOutputPort()) + +actor = vtk.vtkActor() +actor.SetMapper(pdm) + +ren = vtk.vtkRenderer() +renWin = vtk.vtkRenderWindow() +renWin.AddRenderer(ren) + +iren = vtk.vtkRenderWindowInteractor() +iren.SetRenderWindow(renWin) + +ren.AddActor(actor) +ren.SetBackground(1,1,1) +renWin.SetSize(500,500) + +iren.Initialize() +iren.Start() + diff --git a/mrMeshPy/matlabRoutines/pyMeshBuild_ImageData_mac.py b/mrMeshPy/matlabRoutines/pyMeshBuild_ImageData_mac.py new file mode 100755 index 00000000..9c09985b --- /dev/null +++ b/mrMeshPy/matlabRoutines/pyMeshBuild_ImageData_mac.py @@ -0,0 +1,326 @@ +#!/Users/lcladm/anaconda2/bin/python + +# TODO - python path - assuming condo here + +## HEALTH WARNING - BETA CODE IN DEVELOPMENT ## + +''' +This standalone application will build a mesh from a nifti classification file. +To keep the procedure as similar as possible to the way mrMesh used to do this, +we will keep this as a standalon application. Matlab reads in the segmented +nifti file using vistasofts own nifti class handler meshBuild>mrmBuild> +meshBuildFromClass - we just dont use the old build_mesh mex file - we do that +bit and any smoothing in this application and send a mesh struture back to +matlab. + +AG 2017 +''' + +import os,sys +import scipy +import vtk + +from numpy import * +from scipy.io import loadmat, savemat +from vtk.util import numpy_support + +debug = False + +#TODO error handling +fileToLoad = sys.argv[1] +fileToSave = sys.argv[2] + +# load the voxel data that has been dumped to disk +voxels = scipy.io.loadmat(fileToLoad) +mmPerVox = voxels['mmPerVox'][0] +if debug: print mmPerVox + +voxels = voxels['voxels'] #unpack + + +if debug: print voxels + +if debug: print shape(voxels) + +extent = shape(voxels) +if debug: print extent +if debug: print extent[0] +if debug: print extent[1] +if debug: print extent[2] + + + +### ------------------------------------------------------------------------------ +### this is faster but for now exactly replicate the way mrMesh sets up the volume array +# import voxels to vtk +dataImporter = vtk.vtkImageImport() +data_string = voxels.tostring() +dataImporter.CopyImportVoidPointer(data_string, len(data_string)) +dataImporter.SetDataScalarTypeToUnsignedChar() +dataImporter.SetDataExtent(0, extent[2]-1, 0, extent[1]-1, 0, extent[0]-1) # TODO have to work this out +dataImporter.SetWholeExtent(0, extent[2]-1, 0, extent[1]-1, 0, extent[0]-1) # TODO have to work this out +dataImporter.SetDataSpacing(mmPerVox[0],mmPerVox[1],mmPerVox[2]) # TODO have to work this out +dataImporter.Update() +if debug: print dataImporter.GetOutput() + +### ------------------------------------------------------------------------------ +''' + + +### ------- the way mrMesh did it in mesh_build -------------------------------- +pArray = map(ord,voxels.tostring()) #unpack + +pDims = shape(voxels) +scale = mmPerVox +iSizes = [pDims[0]+2, pDims[1]+2, pDims[2]+2] + +nTotalValues = iSizes[0] * iSizes[1] * iSizes[2] + +pClassValues = vtk.vtkUnsignedCharArray() +pClassValues.SetNumberOfValues(nTotalValues) + +pClassData = vtk.vtkStructuredPoints() +pClassData.SetDimensions(iSizes[0], iSizes[1], iSizes[2]) +pClassData.SetOrigin(-scale[0], -scale[1], -scale[2]) #??? +pClassData.SetOrigin(-1, -1, -1) #??? +pClassData.SetSpacing(scale[0], scale[1], scale[2]) + + +for iSrcZ in range(pDims[2]): + + for iSrcY in range(pDims[1]): + + iSrcIndex = iSrcZ * pDims[1] * pDims[0] + iSrcY * pDims[0] + iDstIndex = (iSrcZ+1) * iSizes[1] * iSizes[0] + (iSrcY+1) * iSizes[0] + 1 + + for iSrcX in range(pDims[0]): + + fTemp = int(pArray[iSrcIndex]) + #if debug: print fTemp, 'iSrcIndex', iSrcIndex, 'iDstIndex', iDstIndex + if fTemp>0: + pClassValues.SetValue(iDstIndex, 0) + else: + pClassValues.SetValue(iDstIndex, 1) + + + iSrcIndex+=1 + iDstIndex+=1 + + +pClassData.GetPointData().SetScalars(pClassValues) + +pClassData.Modified() + +if debug: + spw = vtk.vtkStructuredPointsWriter() + spw.SetFileTypeToASCII() + spw.SetInputData(pClassData) + spw.SetFileName("/tmp/test-mrMeshPy-structuredPoints.vtk") + spw.Write() + spw.Update() + + +''' + +### ------ Data volume is loaded and constructed - extract some surgfaces ------------- + +mc = vtk.vtkMarchingCubes() # mc = vtk.vtkContourFilter() #- could use a contour filter instead? + +mc.SetInputConnection(dataImporter.GetOutputPort()) # later - for use with direct imagedata import #mc.SetInputData(pClassData) +mc.SetValue(0,0.5) #extract 0-th surface at 0.5? +mc.ComputeGradientsOff() +mc.ComputeNormalsOff() +mc.ComputeScalarsOff() +mc.Update() + + +if debug: + print mc.GetOutput() + write = vtk.vtkPolyDataWriter() + write.SetFileName('/htmp/test-mrMeshPy-marchingCubesOutput.txt') + write.SetFileTypeToASCII() + write.SetInputData(mc.GetOutput()) + write.Write() + write.Update() + + +# ---- extract center surface - edges are normally extracted to (the cube around the edge of the volume)-------- + +confilter = vtk.vtkPolyDataConnectivityFilter() +confilter.SetInputConnection(mc.GetOutputPort()) +confilter.SetExtractionModeToClosestPointRegion() +confilter.SetClosestPoint(extent[0]/2.0,extent[1]/2.0,extent[2]/2.0) # center of volume +confilter.Update() + + + +# ---- Normals --------------------- + +# normals already computed by mc algorithm so this code is obsolete +normals = vtk.vtkPolyDataNormals() +normals.ComputePointNormalsOn() +normals.SplittingOff() +normals.SetInputConnection(confilter.GetOutputPort()) +#normals.SetInputData(discrete.GetOutput()) +normals.Update() +print normals.GetOutput() + +norm = normals.GetOutput().GetPointData().GetNormals() +output_normals = array(numpy_support.vtk_to_numpy(norm).transpose(),'d') + +####if debug: print output_normals + + + +# ---- Initial vertices - unsmoothed --------------------- +init_verts = normals.GetOutput().GetPoints().GetData() +output_init_verts = array(numpy_support.vtk_to_numpy(init_verts).transpose(),'d') + +if debug: print output_init_verts + + +# ---- Polys (triangles) --------------------- +triangles = normals.GetOutput().GetPolys().GetData() +tmp_triangles = numpy_support.vtk_to_numpy(triangles) + +# N.B. the polygon data returned here have 4 values for poly - the first is the number +# of vertices that describe the polygo (ironically always 3) and the next 3 are the +# indices of the vertices that make up the polygon + +# so first we need to reshape data from a vector +tmp_triangles = reshape(tmp_triangles,(len(tmp_triangles)/4,4)) + +# and then we drop the first column (all 3's) +output_triangles = array((tmp_triangles[:,1:4]).transpose(),'d') #remember zero index here, add one for matlab + +if debug: print output_triangles + + + + + +# -------- smoothed version of mesh ---------------- + +smooth = vtk.vtkSmoothPolyDataFilter() +smooth.SetNumberOfIterations(32) #standard value sused in old mrMesh +smooth.SetRelaxationFactor(0.5) #standard value sused in old mrMesh +smooth.FeatureEdgeSmoothingOff() +smooth.SetFeatureAngle(45) +smooth.SetEdgeAngle(15) +smooth.SetBoundarySmoothing(1) +smooth.SetInputConnection(normals.GetOutputPort()) +smooth.Update() + +# different smoothing option? +''' +smooth = vtk.vtkWindowedSincPolyDataFilter() +smooth.SetInputConnection(mc.GetOutputPort()) +smooth.SetNumberOfIterations(30) +smooth.SetPassBand(0.5) +smooth.SetFeatureAngle(45) +smooth.SetEdgeAngle(15) +smooth.SetBoundarySmoothing(1) +smooth.SetFeatureEdgeSmoothing(0) +smooth.Update() +''' + + +# ---- Vertices - smoothed --------------------- +smooth_verts = smooth.GetOutput().GetPoints().GetData() +output_smooth_verts = array(numpy_support.vtk_to_numpy(smooth_verts).transpose(),'d') + +if debug: print output_smooth_verts + + +# ---- Curvature --------------------- +curvature = vtk.vtkCurvatures() +curvature.SetInputConnection(smooth.GetOutputPort()) +curvature.SetCurvatureTypeToMean() +curvature.Update() + +curv = curvature.GetOutput().GetPointData().GetScalars() +output_curvature = array(numpy_support.vtk_to_numpy(curv).transpose(),'d') + +if debug: print min(output_curvature) +if debug: print max(output_curvature) +if debug: print output_curvature + + + +# -------- colours based on curvature ------------ + +# turn curvature into color +tmp_colors = output_curvature.copy() + +#min_curv = min(tmp_colors) +#max_curv = max(tmp_colors) +#tmp_colors = (tmp_colors -min_curv) / (max_curv-min_curv) *255 + +tmp_colors[tmp_colors>=0] = 160 #standard value sused in old mrMesh +tmp_colors[tmp_colors<0] = 85 #standard value sused in old mrMesh + +output_colors = vstack((tmp_colors, tmp_colors, tmp_colors, ones((1,len(tmp_colors)))*255)) +output_colors = array(output_colors,'d') + +if debug: print output_colors + +# OK we have all the data we need now, lets write it out to file + +data = {} #empty dictionary +data['initVertices'] = output_init_verts +data['initialvertices'] = output_init_verts +data['vertices'] = output_smooth_verts +data['colors'] = output_colors +data['normals'] = output_normals +data['triangles'] = output_triangles +data['curvature'] = output_curvature + +# save it out +savemat(fileToSave,data) + + + +# data have been sent, but let's view them here + +pdm = vtk.vtkPolyDataMapper() +pdm.SetInputConnection(confilter.GetOutputPort()) + +actor = vtk.vtkActor() +actor.SetMapper(pdm) + +ren = vtk.vtkRenderer() +renWin = vtk.vtkRenderWindow() +renWin.AddRenderer(ren) + +iren = vtk.vtkRenderWindowInteractor() +iren.SetRenderWindow(renWin) + +ren.AddActor(actor) +ren.SetBackground(1,1,1) +renWin.SetSize(500,500) + +iren.Initialize() +iren.Start() + + + +pdm = vtk.vtkPolyDataMapper() +pdm.SetInputConnection(curvature.GetOutputPort()) + +actor = vtk.vtkActor() +actor.SetMapper(pdm) + +ren = vtk.vtkRenderer() +renWin = vtk.vtkRenderWindow() +renWin.AddRenderer(ren) + +iren = vtk.vtkRenderWindowInteractor() +iren.SetRenderWindow(renWin) + +ren.AddActor(actor) +ren.SetBackground(1,1,1) +renWin.SetSize(500,500) + +iren.Initialize() +iren.Start() + diff --git a/mrMeshPy/matlabRoutines/pyMeshBuild_mac.py b/mrMeshPy/matlabRoutines/pyMeshBuild_mac.py new file mode 100755 index 00000000..9ab81482 --- /dev/null +++ b/mrMeshPy/matlabRoutines/pyMeshBuild_mac.py @@ -0,0 +1,329 @@ +#!/Users/lcladm/anaconda2/bin/python + +# TODO - python path - assuming condo here + +## HEALTH WARNING - BETA CODE IN DEVELOPMENT ## + +''' +This standalone application will build a mesh from a nifti classification file. +To keep the procedure as similar as possible to the way mrMesh used to do this, +we will keep this as a standalon application. Matlab reads in the segmented +nifti file using vistasofts own nifti class handler meshBuild>mrmBuild> +meshBuildFromClass - we just dont use the old build_mesh mex file - we do that +bit and any smoothing in this application and send a mesh struture back to +matlab. + +AG 2017 +''' + +import os,sys +import scipy +import vtk + +from numpy import * +from scipy.io import loadmat, savemat +from vtk.util import numpy_support + +debug = False + +#TODO error handling +fileToLoad = sys.argv[1] +fileToSave = sys.argv[2] + +# load the voxel data that has been dumped to disk +voxels = scipy.io.loadmat(fileToLoad) +mmPerVox = voxels['mmPerVox'][0] +if debug: print mmPerVox + +voxels = voxels['voxels'] #unpack + + +if debug: print voxels + +if debug: print shape(voxels) + +extent = shape(voxels) +if debug: print extent +if debug: print extent[0] +if debug: print extent[1] +if debug: print extent[2] + + +''' +### ------------------------------------------------------------------------------ +### this is faster but for now exactly replicate the way mrMesh sets up the volume array +# import voxels to vtk +dataImporter = vtk.vtkImageImport() +data_string = voxels.tostring() +dataImporter.CopyImportVoidPointer(data_string, len(data_string)) +dataImporter.SetDataScalarTypeToUnsignedChar() +dataImporter.SetDataExtent(0, extent[2]-1, 0, extent[1]-1, 0, extent[0]-1) # TODO have to work this out +dataImporter.SetWholeExtent(0, extent[2]-1, 0, extent[1]-1, 0, extent[0]-1) # TODO have to work this out +dataImporter.SetDataSpacing(mmPerVox[0],mmPerVox[1],mmPerVox[2]) # TODO have to work this out +dataImporter.Update() +if debug: print dataImporter.GetOutput() + +### ------------------------------------------------------------------------------ +''' + + +### ------- the way mrMesh did it in mesh_build -------------------------------- +pArray = map(ord,voxels.tostring()) #unpack + +pDims = shape(voxels) +scale = mmPerVox +iSizes = [pDims[0]+2, pDims[1]+2, pDims[2]+2] + +nTotalValues = iSizes[0] * iSizes[1] * iSizes[2] + +pClassValues = vtk.vtkUnsignedCharArray() +pClassValues.SetNumberOfValues(nTotalValues) + +pClassData = vtk.vtkStructuredPoints() +pClassData.SetDimensions(iSizes[0], iSizes[1], iSizes[2]) +pClassData.SetOrigin(-scale[0], -scale[1], -scale[2]) #??? +pClassData.SetOrigin(-1, -1, -1) #??? +pClassData.SetSpacing(scale[0], scale[1], scale[2]) + + +for iSrcZ in range(pDims[2]): + + for iSrcY in range(pDims[1]): + + iSrcIndex = iSrcZ * pDims[1] * pDims[0] + iSrcY * pDims[0] + iDstIndex = (iSrcZ+1) * iSizes[1] * iSizes[0] + (iSrcY+1) * iSizes[0] + 1 + + for iSrcX in range(pDims[0]): + + fTemp = int(pArray[iSrcIndex]) + #if debug: print fTemp, 'iSrcIndex', iSrcIndex, 'iDstIndex', iDstIndex + if fTemp>0: + pClassValues.SetValue(iDstIndex, 0) + else: + pClassValues.SetValue(iDstIndex, 1) + + + iSrcIndex+=1 + iDstIndex+=1 + + +pClassData.GetPointData().SetScalars(pClassValues) + +pClassData.Modified() + +if debug: + spw = vtk.vtkStructuredPointsWriter() + spw.SetFileTypeToASCII() + spw.SetInputData(pClassData) + spw.SetFileName("/tmp/test-mrMeshPy-structuredPoints.vtk") + spw.Write() + spw.Update() + + + + +### ------ Data volume is loaded and constructed - extract some surgfaces ------------- + +mc = vtk.vtkMarchingCubes() +# mc = vtk.vtkContourFilter() #- could use a contour filter instead? + +# mc.SetInputConnection(dataImporter.GetOutputPort()) # later - for use with direct imagedata import + +mc.SetInputData(pClassData) +mc.SetValue(0,0.5) #extract 0-th surface at 0.5? +mc.ComputeGradientsOff() +mc.ComputeNormalsOff() +mc.ComputeScalarsOff() +mc.Update() + + +if debug: + print mc.GetOutput() + write = vtk.vtkPolyDataWriter() + write.SetFileName('/htmp/test-mrMeshPy-marchingCubesOutput.txt') + write.SetFileTypeToASCII() + write.SetInputData(mc.GetOutput()) + write.Write() + write.Update() + + +# ---- extract center surface - edges are normally extracted to (the cube around the edge of the volume)-------- + +confilter = vtk.vtkPolyDataConnectivityFilter() +confilter.SetInputConnection(mc.GetOutputPort()) +confilter.SetExtractionModeToClosestPointRegion() +confilter.SetClosestPoint(extent[0]/2.0,extent[1]/2.0,extent[2]/2.0) # center of volume +confilter.Update() + + + +# ---- Normals --------------------- + +# normals already computed by mc algorithm so this code is obsolete +normals = vtk.vtkPolyDataNormals() +normals.ComputePointNormalsOn() +normals.SplittingOff() +normals.SetInputConnection(confilter.GetOutputPort()) +#normals.SetInputData(discrete.GetOutput()) +normals.Update() +print normals.GetOutput() + +norm = normals.GetOutput().GetPointData().GetNormals() +output_normals = array(numpy_support.vtk_to_numpy(norm).transpose(),'d') + +####if debug: print output_normals + + + +# ---- Initial vertices - unsmoothed --------------------- +init_verts = normals.GetOutput().GetPoints().GetData() +output_init_verts = array(numpy_support.vtk_to_numpy(init_verts).transpose(),'d') + +if debug: print output_init_verts + + +# ---- Polys (triangles) --------------------- +triangles = normals.GetOutput().GetPolys().GetData() +tmp_triangles = numpy_support.vtk_to_numpy(triangles) + +# N.B. the polygon data returned here have 4 values for poly - the first is the number +# of vertices that describe the polygo (ironically always 3) and the next 3 are the +# indices of the vertices that make up the polygon + +# so first we need to reshape data from a vector +tmp_triangles = reshape(tmp_triangles,(len(tmp_triangles)/4,4)) + +# and then we drop the first column (all 3's) +output_triangles = array((tmp_triangles[:,1:4]).transpose(),'d') #remember zero index here, add one for matlab + +if debug: print output_triangles + + + + + +# -------- smoothed version of mesh ---------------- + +smooth = vtk.vtkSmoothPolyDataFilter() +smooth.SetNumberOfIterations(32) #standard value sused in old mrMesh +smooth.SetRelaxationFactor(0.5) #standard value sused in old mrMesh +smooth.FeatureEdgeSmoothingOff() +smooth.SetFeatureAngle(45) +smooth.SetEdgeAngle(15) +smooth.SetBoundarySmoothing(1) +smooth.SetInputConnection(normals.GetOutputPort()) +smooth.Update() + +# different smoothing option? +''' +smooth = vtk.vtkWindowedSincPolyDataFilter() +smooth.SetInputConnection(mc.GetOutputPort()) +smooth.SetNumberOfIterations(30) +smooth.SetPassBand(0.5) +smooth.SetFeatureAngle(45) +smooth.SetEdgeAngle(15) +smooth.SetBoundarySmoothing(1) +smooth.SetFeatureEdgeSmoothing(0) +smooth.Update() +''' + + +# ---- Vertices - smoothed --------------------- +smooth_verts = smooth.GetOutput().GetPoints().GetData() +output_smooth_verts = array(numpy_support.vtk_to_numpy(smooth_verts).transpose(),'d') + +if debug: print output_smooth_verts + + +# ---- Curvature --------------------- +curvature = vtk.vtkCurvatures() +curvature.SetInputConnection(smooth.GetOutputPort()) +curvature.SetCurvatureTypeToMean() +curvature.Update() + +curv = curvature.GetOutput().GetPointData().GetScalars() +output_curvature = array(numpy_support.vtk_to_numpy(curv).transpose(),'d') + +if debug: print min(output_curvature) +if debug: print max(output_curvature) +if debug: print output_curvature + + + +# -------- colours based on curvature ------------ + +# turn curvature into color +tmp_colors = output_curvature.copy() + +#min_curv = min(tmp_colors) +#max_curv = max(tmp_colors) +#tmp_colors = (tmp_colors -min_curv) / (max_curv-min_curv) *255 + +tmp_colors[tmp_colors>=0] = 85 #standard value sused in old mrMesh +tmp_colors[tmp_colors<0] = 160 #standard value sused in old mrMesh + +output_colors = vstack((tmp_colors, tmp_colors, tmp_colors, ones((1,len(tmp_colors)))*255)) +output_colors = array(output_colors,'d') + +if debug: print output_colors + +# OK we have all the data we need now, lets write it out to file + +data = {} #empty dictionary +data['initVertices'] = output_init_verts +data['initialvertices'] = output_init_verts +data['vertices'] = output_smooth_verts +data['colors'] = output_colors +data['normals'] = output_normals +data['triangles'] = output_triangles +data['curvature'] = output_curvature + +# save it out +savemat(fileToSave,data) + + + +# data have been sent, but let's view them here + +pdm = vtk.vtkPolyDataMapper() +pdm.SetInputConnection(confilter.GetOutputPort()) + +actor = vtk.vtkActor() +actor.SetMapper(pdm) + +ren = vtk.vtkRenderer() +renWin = vtk.vtkRenderWindow() +renWin.AddRenderer(ren) + +iren = vtk.vtkRenderWindowInteractor() +iren.SetRenderWindow(renWin) + +ren.AddActor(actor) +ren.SetBackground(1,1,1) +renWin.SetSize(500,500) + +iren.Initialize() +iren.Start() + + + +pdm = vtk.vtkPolyDataMapper() +pdm.SetInputConnection(curvature.GetOutputPort()) + +actor = vtk.vtkActor() +actor.SetMapper(pdm) + +ren = vtk.vtkRenderer() +renWin = vtk.vtkRenderWindow() +renWin.AddRenderer(ren) + +iren = vtk.vtkRenderWindowInteractor() +iren.SetRenderWindow(renWin) + +ren.AddActor(actor) +ren.SetBackground(1,1,1) +renWin.SetSize(500,500) + +iren.Initialize() +iren.Start() + diff --git a/mrMeshPy/matlabRoutines/pyMeshBuild_test.py b/mrMeshPy/matlabRoutines/pyMeshBuild_test.py new file mode 100755 index 00000000..7bdaf975 --- /dev/null +++ b/mrMeshPy/matlabRoutines/pyMeshBuild_test.py @@ -0,0 +1,372 @@ +#!/usr/bin/python + +# TODO - python path - assuming condo here + +## HEALTH WARNING - BETA CODE IN DEVELOPMENT ## + +''' +This standalone application will build a mesh from a nifti classification file. +To keep the procedure as similar as possible to the way mrMesh used to do this, +we will keep this as a standalon application. Matlab reads in the segmented +nifti file using vistasofts own nifti class handler meshBuild>mrmBuild> +meshBuildFromClass - we just dont use the old build_mesh mex file - we do that +bit and any smoothing in this application and send a mesh struture back to +matlab. + +AG 2017 +''' + +import os,sys +import scipy +import vtk + +from numpy import * +from scipy.io import loadmat, savemat +from vtk.util import numpy_support + +debug = True + +#TODO error handling +fileToLoad = sys.argv[1] +fileToSave = sys.argv[2] + +# load the voxel data that has been dumped to disk +voxels = scipy.io.loadmat(fileToLoad) +mmPerVox = voxels['mmPerVox'][0] +if debug: print mmPerVox + +voxels = voxels['voxels'] #unpack + + +if debug: print voxels + +if debug: print shape(voxels) + +extent = shape(voxels) +if debug: print extent +if debug: print extent[0] +if debug: print extent[1] +if debug: print extent[2] + + +''' +### ------------------------------------------------------------------------------ +### this is faster but for now exactly replicate the way mrMesh sets up the volume array +# import voxels to vtk +dataImporter = vtk.vtkImageImport() +data_string = voxels.tostring() +dataImporter.CopyImportVoidPointer(data_string, len(data_string)) +dataImporter.SetDataScalarTypeToUnsignedChar() +dataImporter.SetDataExtent(0, extent[2]-1, 0, extent[1]-1, 0, extent[0]-1) # TODO have to work this out +dataImporter.SetWholeExtent(0, extent[2]-1, 0, extent[1]-1, 0, extent[0]-1) # TODO have to work this out +dataImporter.SetDataSpacing(mmPerVox[0],mmPerVox[1],mmPerVox[2]) # TODO have to work this out +dataImporter.Update() +if debug: print dataImporter.GetOutput() + +### ------------------------------------------------------------------------------ +''' + + +### ------- the way mrMesh did it in mesh_build -------------------------------- +pArray = map(ord,voxels.tostring()) #unpack + +pDims = shape(voxels) +scale = mmPerVox +iSizes = [pDims[0]+2, pDims[1]+2, pDims[2]+2] + +nTotalValues = iSizes[0] * iSizes[1] * iSizes[2] + +pClassValues = vtk.vtkUnsignedCharArray() +pClassValues.SetNumberOfValues(nTotalValues) + +pClassData = vtk.vtkStructuredPoints() +pClassData.SetDimensions(iSizes[0], iSizes[1], iSizes[2]) +pClassData.SetOrigin(-scale[0], -scale[1], -scale[2]) #??? +pClassData.SetOrigin(-1, -1, -1) #??? +pClassData.SetSpacing(scale[0], scale[1], scale[2]) + + +for iSrcZ in range(pDims[2]): + + for iSrcY in range(pDims[1]): + + iSrcIndex = iSrcZ * pDims[1] * pDims[0] + iSrcY * pDims[0] + iDstIndex = (iSrcZ+1) * iSizes[1] * iSizes[0] + (iSrcY+1) * iSizes[0] + 1 + + for iSrcX in range(pDims[0]): + + fTemp = int(pArray[iSrcIndex]) + #if debug: print fTemp, 'iSrcIndex', iSrcIndex, 'iDstIndex', iDstIndex + if fTemp>0: + pClassValues.SetValue(iDstIndex, 0) + else: + pClassValues.SetValue(iDstIndex, 1) + + + iSrcIndex+=1 + iDstIndex+=1 + + +pClassData.GetPointData().SetScalars(pClassValues) + +pClassData.Modified() + +if debug: + spw = vtk.vtkStructuredPointsWriter() + spw.SetFileTypeToASCII() + spw.SetInputData(pClassData) + spw.SetFileName("/tmp/test-mrMeshPy-structuredPoints.vtk") + spw.Write() + spw.Update() + + + + +### ------ Data volume is loaded and constructed - extract some surfaces ------------- + +mc = vtk.vtkMarchingCubes() +# mc = vtk.vtkContourFilter() #- could use a contour filter instead? + +# mc.SetInputConnection(dataImporter.GetOutputPort()) # later - for use with direct imagedata import + +mc.SetInputData(pClassData) +mc.SetValue(0,0.5) #extract 0-th surface at 0.5? +mc.ComputeGradientsOff() +mc.ComputeNormalsOff() +mc.ComputeScalarsOff() +mc.Update() + + +if debug: + print mc.GetOutput() + write = vtk.vtkPolyDataWriter() + write.SetFileName('/htmp/test-mrMeshPy-marchingCubesOutput.txt') + write.SetFileTypeToASCII() + write.SetInputData(mc.GetOutput()) + write.Write() + write.Update() + + +# ---- To exactly replicate the mrMesh routines, we also need to reverse the triangles (??) +# its unclear why they did this but we will do it because they did +# its a little trcky in vtk so we pull the vertex indices out of the vtk cell array +# and then rebuild it after shuffling the order + +polyDat = mc.GetOutput() +theCellArray = polyDat.GetPolys() + +numpyCellValues = numpy_support.vtk_to_numpy(theCellArray.GetData()) #convert to numpy +if debug: print(numpyCellValues) + +nar = numpyCellValues.reshape(theCellArray.GetNumberOfCells(),4) # reshape to n x 4 +nar2 = zeros(nar.shape) # create a temporary holder + +nar2[:,0]=nar[:,0] # keep col 1 in place (always the number 3) +nar2[:,1]=nar[:,3] # swap vertex 3 with 1 +nar2[:,2]=nar[:,2] # keep vertex 2 in place +nar2[:,3]=nar[:,1] # swap vertex 1 with 3 + +#triangles now reversed + +nar2r = nar2.reshape(1,theCellArray.GetNumberOfCells()*4) # reshape to a vector + +# overwite the original values in the cell array +for i in range(theCellArray.GetNumberOfCells()*4): + theCellArray.ReplaceCell(i-1,1,[int(nar2r[0][i])]) #TODO why off by one? + +na = numpy_support.vtk_to_numpy(theCellArray.GetData()) +if debug: print(na) + +polyDat.Modified() + + + +# ---- extract just "center surface" - edges are normally extracted too (the cube around the edge of the volume) -------- + +confilter = vtk.vtkPolyDataConnectivityFilter() +confilter.SetInputData(polyDat) +confilter.SetExtractionModeToClosestPointRegion() +confilter.SetClosestPoint(extent[0]/2.0,extent[1]/2.0,extent[2]/2.0) # center of volume +confilter.Update() + + + +# ---- Normals --------------------- + +# normals already computed by mc algorithm so this code is obsolete +normals = vtk.vtkPolyDataNormals() +normals.ComputePointNormalsOn() +normals.SplittingOff() +normals.SetInputConnection(confilter.GetOutputPort()) +######normals.SetInputData(polyDat) +normals.Update() +print normals.GetOutput() + +norm = normals.GetOutput().GetPointData().GetNormals() +output_normals = array(numpy_support.vtk_to_numpy(norm).transpose(),'d') + +####if debug: print output_normals + + + +# ---- Initial vertices - unsmoothed --------------------- +init_verts = normals.GetOutput().GetPoints().GetData() +output_init_verts = array(numpy_support.vtk_to_numpy(init_verts).transpose(),'d') + +if debug: print output_init_verts + + +# ---- Polys (triangles) --------------------- +triangles = normals.GetOutput().GetPolys().GetData() +tmp_triangles = numpy_support.vtk_to_numpy(triangles) + +# N.B. the polygon data returned here have 4 values for poly - the first is the number +# of vertices that describe the polygo (ironically always 3) and the next 3 are the +# indices of the vertices that make up the polygon + +# so first we need to reshape data from a vector +tmp_triangles = reshape(tmp_triangles,(len(tmp_triangles)/4,4)) + +# and then we drop the first column (all 3's) +output_triangles = array((tmp_triangles[:,1:4]).transpose(),'d') #remember zero index here, add one for matlab + +if debug: print output_triangles + + + + + +# -------- smoothed version of mesh ---------------- + +smooth = vtk.vtkSmoothPolyDataFilter() +smooth.SetNumberOfIterations(32) #standard value sused in old mrMesh +smooth.SetRelaxationFactor(0.5) #standard value sused in old mrMesh +smooth.FeatureEdgeSmoothingOff() +smooth.SetFeatureAngle(45) +smooth.SetEdgeAngle(15) +smooth.SetBoundarySmoothing(1) +smooth.SetInputConnection(normals.GetOutputPort()) +smooth.Update() + +# different smoothing option? +''' +smooth = vtk.vtkWindowedSincPolyDataFilter() +smooth.SetInputConnection(mc.GetOutputPort()) +smooth.SetNumberOfIterations(30) +smooth.SetPassBand(0.5) + +# different smoothing option? +smooth.SetFeatureAngle(45) +smooth.SetEdgeAngle(15) +smooth.SetBoundarySmoothing(1) +smooth.SetFeatureEdgeSmoothing(0) +smooth.Update() +''' + + +#### NEW NORMALS !! +norm = smooth.GetOutput().GetPointData().GetNormals() +output_normals = array(numpy_support.vtk_to_numpy(norm).transpose(),'d') + + + +# ---- Vertices - smoothed --------------------- +smooth_verts = smooth.GetOutput().GetPoints().GetData() +output_smooth_verts = array(numpy_support.vtk_to_numpy(smooth_verts).transpose(),'d') + +if debug: print output_smooth_verts + + +# ---- Curvature --------------------- +curvature = vtk.vtkCurvatures() +curvature.SetInputConnection(smooth.GetOutputPort()) +curvature.SetCurvatureTypeToMean() +curvature.Update() + +curv = curvature.GetOutput().GetPointData().GetScalars() +output_curvature = array(numpy_support.vtk_to_numpy(curv).transpose(),'d') + +if debug: print min(output_curvature) +if debug: print max(output_curvature) +if debug: print output_curvature + + + +# -------- colours based on curvature ------------ + +# turn curvature into color +tmp_colors = output_curvature.copy() + +#min_curv = min(tmp_colors) +#max_curv = max(tmp_colors) +#tmp_colors = (tmp_colors -min_curv) / (max_curv-min_curv) *255 + +tmp_colors[tmp_colors>=0] = 85 #standard value sused in old mrMesh +tmp_colors[tmp_colors<0] = 160 #standard value sused in old mrMesh + +output_colors = vstack((tmp_colors, tmp_colors, tmp_colors, ones((1,len(tmp_colors)))*255)) +output_colors = array(output_colors,'d') + +if debug: print output_colors + +# OK we have all the data we need now, lets write it out to file + +data = {} #empty dictionary +data['initVertices'] = output_init_verts +data['initialvertices'] = output_init_verts +data['vertices'] = output_smooth_verts +data['colors'] = output_colors +data['normals'] = output_normals +data['triangles'] = output_triangles +data['curvature'] = output_curvature + +# save it out +savemat(fileToSave,data) + + + +# data have been sent, but let's view them here + +pdm = vtk.vtkPolyDataMapper() +pdm.SetInputConnection(confilter.GetOutputPort()) +######pdm.SetInputData(normals.GetOutput()) + +actor = vtk.vtkActor() +actor.SetMapper(pdm) + +ren = vtk.vtkRenderer() +renWin = vtk.vtkRenderWindow() +renWin.AddRenderer(ren) + +iren = vtk.vtkRenderWindowInteractor() +iren.SetRenderWindow(renWin) + +ren.AddActor(actor) +ren.SetBackground(1,1,1) +renWin.SetSize(500,500) + +iren.Initialize() +iren.Start() + + + +pdm = vtk.vtkPolyDataMapper() +pdm.SetInputConnection(curvature.GetOutputPort()) + +actor = vtk.vtkActor() +actor.SetMapper(pdm) + +ren = vtk.vtkRenderer() +renWin = vtk.vtkRenderWindow() +renWin.AddRenderer(ren) + +iren = vtk.vtkRenderWindowInteractor() +iren.SetRenderWindow(renWin) + +ren.AddActor(actor) +ren.SetBackground(1,1,1) +renWin.SetSize(500,500) + +iren.Initialize() +iren.Start() + diff --git a/mrMeshPy/matlabRoutines/scratchpad.m b/mrMeshPy/matlabRoutines/scratchpad.m new file mode 100755 index 00000000..f5c135eb --- /dev/null +++ b/mrMeshPy/matlabRoutines/scratchpad.m @@ -0,0 +1,35 @@ +% a scratchpad of test commands used during mrMeshPy development + +% % % Callback: MSH = viewGet(VOLUME{1}, 'Mesh'); +% % % vertexGrayMap = mrmMapVerticesToGray( meshGet(MSH, 'vertices'), viewGet(VOLUME{1}, 'nodes'), viewGet(VOLUME{1}, 'mmPerVox'), viewGet(VOLUME{1}, 'edges') ); +% % % MSH = meshSet(MSH, 'vertexgraymap', vertexGrayMap); +% % % VOLUME{1} = viewSet(VOLUME{1}, 'Mesh', MSH); +% % % clear MSH vertexGrayMap +% % % +% % % +% % % +% % % msh = mrmBuildMeshMatlab(cFile,1.0); +% % % msh.colors(4,:) = 255; +% % % msh.vertexGrayMap = mrmMapVerticesToGray( meshGet(msh, 'vertices'), viewGet(VOLUME{1}, 'nodes'), viewGet(VOLUME{1}, 'mmPerVox'), viewGet(VOLUME{1}, 'edges') ) +% % % VOLUME{1} = rmfield(VOLUME{1},'mesh') + + +filename = '/groups/Projects/P1252/Data/Anatomy/R3517/key_files/t1_class_edit.nii.gz'; +msh = meshBuildFromNiftiClass_mrMeshPy(filename,'right'); %puts voxels in workspace +voxels = permute(voxels,[3,2,1]); +save /tmp/voxels.mat voxels; +system('/groups/examples/mrMeshPy/mrMeshPy/matlabRoutines/launchMeshBuild.sh /groups/examples/mrMeshPy/testCode/testMeshBuild.py') +msh = load('/tmp/temp.mat'); +msh = meshFormat(msh); % Converts old format to new. +msh.initialvertices = msh.initVertices; +vertices = meshGet(msh,'vertices'); +msh = meshSet(msh,'origin',-mean(vertices,2)'); +msh = meshSet(msh,'mmPerVox',mmPerVox); +vertexGrayMap = mrmMapVerticesToGray(... + meshGet(msh, 'initialvertices'),... + viewGet(vw, 'nodes'),... + viewGet(vw, 'mmPerVox'),... + viewGet(vw, 'edges')); +msh = meshSet(msh, 'vertexgraymap', vertexGrayMap); +save ./test3.mat msh + diff --git a/mrMeshPy/mp_Commands.py b/mrMeshPy/mp_Commands.py new file mode 100644 index 00000000..a3f7badd --- /dev/null +++ b/mrMeshPy/mp_Commands.py @@ -0,0 +1,146 @@ +#!/usr/bin/python + +''' +Command module for mrMeshPy viewer + +Commands and data are passed here from matab via the mrMeshPyServer. + +Matlab sends a string in one transaction giving a command to the +visualisation module. This command either performs an explicit +function in the viewer (e.g. rotate the camera 90 degrees) or the +commnd describes the configuration/content of a large data chunk to +will be sent in the subsequent transaction so that we know how to +unpack the data, and how to process it (e.g. 70,000 floating point +numbers which are scalar values to show as an amplitude map. + +Each command string is interpreted by the mp_commandInterpret module. + +N.B. - currently command strings have a maximum length of 1024 bytes. + +Commands are specifically ordered, semi-colon seperated strings which are +unpacked to describe what the user is trying to do / send from matlab. +Commands have a MINIMUM LENGTH of 6 arguments and have the following +structure and item order (zero-indexed) + +0 - "cmd" -- always this, identifies it as a cmd :) +1-3 - place holders +4 - commandName - should match a command in mp_Commands file +5 - theMeshInstance - integer pointing to the the mesh window that we + want to operate on +6 onwards - commandArgs - a list of comma-separated pairs of arguments + to characterise the processing of the incoming data + blob or apply some settings to the viewport - + CAN BE EMPTY but must be set to [] + + + +Andre' Gouws 2017 +''' + +import vtk + + +import scipy.io #so we can read .mat files +import vtk +import vtk.util.numpy_support +from numpy import * + + +#local modules +from mp_setupVTKWindow import mrMeshVTKWindow +from mp_VTKRoutines import * +from mp_SendFunctions import * + + +debug = True + + +# master command handler + +def run_mp_command(commandName, commandArgs, theMeshInstance, mainWindowUI, the_TCPserver): + + if commandName == 'loadNewMesh': + mainWindowUI.statusbar.showMessage(' ... attempting to load new mesh ...') + + # TODO - index will now be a new entry at the end of the exisitng .ui.vtkInstances list + newIndex = len(mainWindowUI.vtkInstances) #could be zero + + # create an entry in the vtkDict to link the unique mesh ID to where it is stored + # in the vtkInstances list + mainWindowUI.vtkDict[theMeshInstance] = newIndex + + # add a new tab with a new wVTK window + mrMeshVTKWindow(mainWindowUI, theMeshInstance, 'None') + + mainWindowUI.tabWidget.setCurrentIndex(newIndex) #zero indexed + mainWindowUI.tabWidget.update() + + #load data and generate the mesh + loadNewMesh(theMeshInstance, commandArgs, mainWindowUI, the_TCPserver) + mainWindowUI.statusbar.showMessage(' ... New mesh Loaded ...') + #the_TCPserver.socket.write(str('send useful message back here TODO')) + the_TCPserver.socket.write(str('1001')) + if debug: print mainWindowUI.vtkDict + + + + elif commandName == 'smoothMesh': + mainWindowUI.statusbar.showMessage(' ... attempting to smooth mesh with id %s ...' %(theMeshInstance)) + #load data and generate the mesh + err = smoothMesh(theMeshInstance, commandArgs, mainWindowUI, the_TCPserver) + if err == 0: + mainWindowUI.statusbar.showMessage(' ... Finished smoothing mesh with id %s ...' %(theMeshInstance)) + the_TCPserver.socket.write(str('Mesh smooth complete')) + else: + mainWindowUI.statusbar.showMessage(' ... Error trying to smooth mesh with id %s ...' %(theMeshInstance)) + the_TCPserver.socket.write(str('Mesh smooth failed')) + + + elif commandName == 'updateMeshData': + mainWindowUI.statusbar.showMessage(' ... updating mesh with id %s with current View settings ...' %(theMeshInstance)) + #load data and send to the mesh + err = updateMeshData(theMeshInstance, commandArgs, mainWindowUI, the_TCPserver) + if err == 0: + mainWindowUI.statusbar.showMessage(' ... Finished: updated data for mesh id %s ...' %(theMeshInstance)) + the_TCPserver.socket.write(str('Mesh update complete')) + else: + mainWindowUI.statusbar.showMessage(' ... Error trying to update mesh with id %s ...' %(theMeshInstance)) + the_TCPserver.socket.write(str('Mesh update failed')) + + + elif commandName == 'checkMeshROI': + mainWindowUI.statusbar.showMessage(' ... MATLAB requested an ROI from mesh id %s ...' %(theMeshInstance)) + #get roi data (if exists) and send to matlab + error = sendROIInfo(theMeshInstance, commandArgs, mainWindowUI, the_TCPserver) #returns 1 or 0 + if error == 0: + mainWindowUI.statusbar.showMessage(' ... ROI ready to send to MATLAB from mesh id %s...' %(theMeshInstance)) + else: + mainWindowUI.statusbar.showMessage(' ... No ROI to send to MATLAB from mesh id %s...' %(theMeshInstance)) + the_TCPserver.socket.write(str('send useful message back here TODO')) + + + elif commandName == 'sendROIVertices': + mainWindowUI.statusbar.showMessage(' ... MATLAB requested an ROI from mesh id %s ...' %(theMeshInstance)) + #get roi data (if exists) and send to matlab + error = sendROIVertices(theMeshInstance, commandArgs, mainWindowUI, the_TCPserver) #returns 1 or 0 + if error == 0: + mainWindowUI.statusbar.showMessage(' ... ROI ready to send to MATLAB from mesh id %s...' %(theMeshInstance)) + else: + mainWindowUI.statusbar.showMessage(' ... No ROI to send to MATLAB from mesh id %s...' %(theMeshInstance)) + the_TCPserver.socket.write(str('send useful message back here TODO')) + + + elif commandName == 'rotateMeshAnimation': + mainWindowUI.statusbar.showMessage(' ... doing rotation animation ...') + #really just for testing initially + rotateMeshAnimation(theMeshInstance, commandArgs, mainWindowUI, the_TCPserver) + mainWindowUI.statusbar.showMessage(' ... Finished rotation animation ...') + the_TCPserver.socket.write(str('send useful message back here TODO')) + + + else: + print('mrMeshPy received a command it did not recognise') + the_TCPserver.socket.write(str('send cmd error message back here TODO')) + + + diff --git a/mrMeshPy/mp_MainWindow.py b/mrMeshPy/mp_MainWindow.py new file mode 100644 index 00000000..ad44b91d --- /dev/null +++ b/mrMeshPy/mp_MainWindow.py @@ -0,0 +1,102 @@ +# -*- coding: utf-8 -*- + +# Form implementation generated from reading ui file 'mp_MainWindow.ui' +# +# Created by: PyQt5 UI code generator 5.7 +# +# WARNING! All changes made in this file will be lost! + +from PyQt5 import QtCore, QtGui, QtWidgets + +class Ui_MrMeshMainWindow(object): + def setupUi(self, MrMeshMainWindow): + MrMeshMainWindow.setObjectName("MrMeshMainWindow") + MrMeshMainWindow.resize(708, 550) + self.centralwidget = QtWidgets.QWidget(MrMeshMainWindow) + sizePolicy = QtWidgets.QSizePolicy(QtWidgets.QSizePolicy.Expanding, QtWidgets.QSizePolicy.Expanding) + sizePolicy.setHorizontalStretch(0) + sizePolicy.setVerticalStretch(0) + sizePolicy.setHeightForWidth(self.centralwidget.sizePolicy().hasHeightForWidth()) + self.centralwidget.setSizePolicy(sizePolicy) + self.centralwidget.setMinimumSize(QtCore.QSize(662, 0)) + self.centralwidget.setObjectName("centralwidget") + self.horizontalLayout = QtWidgets.QHBoxLayout(self.centralwidget) + self.horizontalLayout.setObjectName("horizontalLayout") + self.tabWidget = QtWidgets.QTabWidget(self.centralwidget) + sizePolicy = QtWidgets.QSizePolicy(QtWidgets.QSizePolicy.MinimumExpanding, QtWidgets.QSizePolicy.MinimumExpanding) + sizePolicy.setHorizontalStretch(0) + sizePolicy.setVerticalStretch(0) + sizePolicy.setHeightForWidth(self.tabWidget.sizePolicy().hasHeightForWidth()) + self.tabWidget.setSizePolicy(sizePolicy) + self.tabWidget.setAutoFillBackground(False) + self.tabWidget.setTabsClosable(False) + self.tabWidget.setMovable(False) + self.tabWidget.setObjectName("tabWidget") + self.horizontalLayout.addWidget(self.tabWidget) + MrMeshMainWindow.setCentralWidget(self.centralwidget) + self.statusbar = QtWidgets.QStatusBar(MrMeshMainWindow) + self.statusbar.setObjectName("statusbar") + MrMeshMainWindow.setStatusBar(self.statusbar) + self.menubar = QtWidgets.QMenuBar(MrMeshMainWindow) + self.menubar.setGeometry(QtCore.QRect(0, 0, 708, 24)) + self.menubar.setNativeMenuBar(False) + self.menubar.setObjectName("menubar") + self.menuMenu = QtWidgets.QMenu(self.menubar) + self.menuMenu.setObjectName("menuMenu") + self.menuROIs = QtWidgets.QMenu(self.menubar) + self.menuROIs.setObjectName("menuROIs") + self.menuTesting = QtWidgets.QMenu(self.menubar) + self.menuTesting.setObjectName("menuTesting") + self.menuExport = QtWidgets.QMenu(self.menubar) + self.menuExport.setObjectName("menuExport") + MrMeshMainWindow.setMenuBar(self.menubar) + self.actionSoon = QtWidgets.QAction(MrMeshMainWindow) + self.actionSoon.setObjectName("actionSoon") + self.actionEnableDraw = QtWidgets.QAction(MrMeshMainWindow) + self.actionEnableDraw.setObjectName("actionEnableDraw") + self.actionDisableDraw = QtWidgets.QAction(MrMeshMainWindow) + self.actionDisableDraw.setObjectName("actionDisableDraw") + self.actionCloseAndFill = QtWidgets.QAction(MrMeshMainWindow) + self.actionCloseAndFill.setObjectName("actionCloseAndFill") + self.actionSend_10_bytes_TCP = QtWidgets.QAction(MrMeshMainWindow) + self.actionSend_10_bytes_TCP.setObjectName("actionSend_10_bytes_TCP") + self.actionStart_New_ROI = QtWidgets.QAction(MrMeshMainWindow) + self.actionStart_New_ROI.setObjectName("actionStart_New_ROI") + self.actionExport_Surface_to_stl = QtWidgets.QAction(MrMeshMainWindow) + self.actionExport_Surface_to_stl.setObjectName("actionExport_Surface_to_stl") + self.actionUnder_development = QtWidgets.QAction(MrMeshMainWindow) + self.actionUnder_development.setObjectName("actionUnder_development") + self.menuMenu.addAction(self.actionUnder_development) + self.menuROIs.addAction(self.actionEnableDraw) + self.menuROIs.addAction(self.actionDisableDraw) + self.menuROIs.addAction(self.actionCloseAndFill) + self.menuROIs.addSeparator() + self.menuROIs.addSeparator() + self.menuROIs.addAction(self.actionStart_New_ROI) + self.menuTesting.addAction(self.actionSend_10_bytes_TCP) + self.menuExport.addAction(self.actionExport_Surface_to_stl) + self.menubar.addAction(self.menuMenu.menuAction()) + self.menubar.addAction(self.menuROIs.menuAction()) + self.menubar.addAction(self.menuTesting.menuAction()) + self.menubar.addAction(self.menuExport.menuAction()) + + self.retranslateUi(MrMeshMainWindow) + self.tabWidget.setCurrentIndex(-1) + QtCore.QMetaObject.connectSlotsByName(MrMeshMainWindow) + + def retranslateUi(self, MrMeshMainWindow): + _translate = QtCore.QCoreApplication.translate + MrMeshMainWindow.setWindowTitle(_translate("MrMeshMainWindow", "mrMeshPy v0.1")) + self.menuMenu.setTitle(_translate("MrMeshMainWindow", "Menu")) + self.menuROIs.setTitle(_translate("MrMeshMainWindow", "ROIs")) + self.menuTesting.setTitle(_translate("MrMeshMainWindow", "Testing")) + self.menuExport.setTitle(_translate("MrMeshMainWindow", "Export")) + self.actionSoon.setText(_translate("MrMeshMainWindow", "Soon")) + self.actionEnableDraw.setText(_translate("MrMeshMainWindow", "Enable drawing")) + self.actionDisableDraw.setText(_translate("MrMeshMainWindow", "Disable drawing")) + self.actionCloseAndFill.setText(_translate("MrMeshMainWindow", "Close and Fill")) + self.actionSend_10_bytes_TCP.setText(_translate("MrMeshMainWindow", "Send 10 bytes TCP")) + self.actionStart_New_ROI.setText(_translate("MrMeshMainWindow", "Start New ROI")) + self.actionExport_Surface_to_stl.setText(_translate("MrMeshMainWindow", "Export Surface to .stl")) + self.actionUnder_development.setText(_translate("MrMeshMainWindow", "Under development")) + diff --git a/mrMeshPy/mp_MainWindow.ui b/mrMeshPy/mp_MainWindow.ui new file mode 100644 index 00000000..3cbf9eb6 --- /dev/null +++ b/mrMeshPy/mp_MainWindow.ui @@ -0,0 +1,144 @@ + + + MrMeshMainWindow + + + + 0 + 0 + 708 + 550 + + + + mrMeshPy v0.1 + + + + + 0 + 0 + + + + + 662 + 0 + + + + + + + + 0 + 0 + + + + false + + + -1 + + + false + + + false + + + + + + + + + + 0 + 0 + 708 + 24 + + + + false + + + + Menu + + + + + + ROIs + + + + + + + + + + + Testing + + + + + + Export + + + + + + + + + + + Soon + + + + + Enable drawing + + + + + Disable drawing + + + + + Close and Fill + + + + + Send 10 bytes TCP + + + + + Start New ROI + + + + + Export Surface to .stl + + + + + Under development + + + + + + diff --git a/mrMeshPy/mp_MenuCallbacks.py b/mrMeshPy/mp_MenuCallbacks.py new file mode 100644 index 00000000..6ba99419 --- /dev/null +++ b/mrMeshPy/mp_MenuCallbacks.py @@ -0,0 +1,130 @@ +#!/usr/bin/python + +''' +Drop menu callback handlers for mrMeshPy + +Andre' Gouws 2017 + +#TODO - workaround that doesn't require the parent UI to be global - needed for drawing routines +''' + + +import vtk + +from mp_VTKDrawing import * + +debug = True + + +global the_parent_UI + + +def Ui_setupMenuBarCallbacks(parent_UI): + + global the_parent_UI + + the_parent_UI = parent_UI; + + if debug: print('setting up menu bar callbacks') + + parent_UI.actionEnableDraw.triggered.connect(cb_MenuEnableDraw) + parent_UI.actionDisableDraw.triggered.connect(cb_MenuDisableDraw) + parent_UI.actionCloseAndFill.triggered.connect(cb_MenuCloseAndFill) + parent_UI.actionStart_New_ROI.triggered.connect(startNewROI) + parent_UI.actionSend_10_bytes_TCP.triggered.connect(send10Bytes) + parent_UI.actionExport_Surface_to_stl.triggered.connect(cb_ExportToSTL) + + + +def testMenuMessage(): + print('test message from mrMeshPy menu! - wait for it!') + global the_parent_UI + + +def cb_MenuEnableDraw(): + global the_parent_UI + #enable draw mode + currentIndex = the_parent_UI.tabWidget.currentIndex() + + the_parent_UI.vtkInstances[currentIndex]._Iren.inDrawMode = 1 #TODO + the_parent_UI.vtkInstances[currentIndex].inDrawMode = the_parent_UI.vtkInstances[currentIndex]._Iren.inDrawMode + style = vtk.vtkInteractorStyleUser() + the_parent_UI.vtkInstances[currentIndex].SetInteractorStyle(style) + the_parent_UI.statusbar.showMessage("In DRAW ROI Mode! - rotation/zoom disabled - (Menu-ROIs-Disable Draw to continue) ..") + + # it has become apparent that users may start drawing a new ROI before an old one is removed. + # this can cause problems if the old temporary overlay surface is still present + # what we'll do is detect if a closed roi surface is present, and - if so - remove it. + + try: + #remove the overlaid roi actor + the_parent_UI.vtkInstances[currentIndex]._Iren.ren.RemoveActor(the_parent_UI.vtkInstances[currentIndex]._Iren.roiActor) + print "Removed an exiting temporary ROI surface to start a new one." + except: + pass + + +def cb_ExportToSTL(): + global the_parent_UI + # export the current surface to stl file format + # get the current viewed instance + currentIndex = the_parent_UI.tabWidget.currentIndex() + + # extract the polydata from the surface as it is currently rendered for export + # get the actor + act = the_parent_UI.vtkInstances[currentIndex].curr_actor + polydata = act.GetMapper().GetInput() + polydata.Modified() + + # open an stl exporter + stl = vtk.vtkSTLWriter() + stl.SetInputData(polydata) + stl.SetFileName("/tmp/brain.stl") + stl.Write() + +def cb_MenuDisableDraw(): + global the_parent_UI + #disable draw mode + currentIndex = the_parent_UI.tabWidget.currentIndex() + + the_parent_UI.vtkInstances[currentIndex]._Iren.inDrawMode = 0 #TODO + the_parent_UI.vtkInstances[currentIndex].inDrawMode = the_parent_UI.vtkInstances[currentIndex]._Iren.inDrawMode + style = vtk.vtkInteractorStyleTrackballCamera() + the_parent_UI.vtkInstances[currentIndex].SetInteractorStyle(style) + the_parent_UI.statusbar.showMessage("Draw mode disabled - normal rotation/zoom resumed ..") + +def cb_MenuCloseAndFill(): + global the_parent_UI + currentIndex = the_parent_UI.tabWidget.currentIndex() + obj = the_parent_UI.vtkInstances[currentIndex]._Iren + obj.inDrawMode = 1 #tmp reset + drawingMakeROI(obj, None) + style = vtk.vtkInteractorStyleTrackballCamera() + the_parent_UI.vtkInstances[currentIndex].SetInteractorStyle(style) + the_parent_UI.statusbar.showMessage("Closed and filled ROI ..") + + +def startNewROI(): + global the_parent_UI + currentIndex = the_parent_UI.tabWidget.currentIndex() + try: + #remove the overlaid roi actor + the_parent_UI.vtkInstances[currentIndex]._Iren.ren.RemoveActor(the_parent_UI.vtkInstances[currentIndex]._Iren.roiActor) + + # reset all pointers to ROI data + the_parent_UI.vtkInstances[currentIndex]._Iren.inDrawMode = 0 #TODO + the_parent_UI.vtkInstances[currentIndex]._Iren.pickedPointIds = [] #place holder for picked vtk point IDs so we can track + the_parent_UI.vtkInstances[currentIndex]._Iren.pickedPointOrigValues = [] #place holder for picked vtk point IDs so we can track + the_parent_UI.vtkInstances[currentIndex]._Iren.pickedPoints = vtk.vtkPoints() #place holder for picked vtk point IDs so we can track + style = vtk.vtkInteractorStyleTrackballCamera() + the_parent_UI.vtkInstances[currentIndex].SetInteractorStyle(style) + the_parent_UI.statusbar.showMessage("Old ROI drawing removed. New ROI started - enable drawing mode to continue ...") + except: + print "Error while trying to remove an exisiting ROI actor - does one actually exist yet?" + + +def send10Bytes(): + global the_parent_UI + the_parent_UI.TCP_server.socket.write('Ready,uint8,%i' %129860*1000) + + diff --git a/mrMeshPy/mp_SendFunctions.py b/mrMeshPy/mp_SendFunctions.py new file mode 100644 index 00000000..048ffee0 --- /dev/null +++ b/mrMeshPy/mp_SendFunctions.py @@ -0,0 +1,50 @@ +#!/usr/bin/python + +''' +Send module for mrMeshPy viewer + +The mrMeshPyServer uses the routines here to send data back to mrVista over the TCP socket. + +Andre' Gouws 2017 +''' + +from numpy import * +import time + +debug = True + + + +def sendROIInfo(theMeshInstance, commandArgs, mainWindowUI, the_TCPserver): + # when a request for ROI info comes in, this sends an initial confirmation that an ROI is available + + try: + targetVTKWindow = mainWindowUI.vtkInstances[mainWindowUI.vtkDict[theMeshInstance]] + data_to_send = targetVTKWindow._Iren.filledROIPoints + + num2send = len(data_to_send) + if debug: print(num2send) + + the_TCPserver.socket.write(str('RoiReady,double,%i,' %num2send)) + return 0 + except: + the_TCPserver.socket.write(str('NoROIData')) + return 1 + + + +def sendROIVertices(theMeshInstance, commandArgs, mainWindowUI, the_TCPserver): + # when a request for ROI info comes we send the data after the confirmation (sendROIInfo above) + targetVTKWindow = mainWindowUI.vtkInstances[mainWindowUI.vtkDict[theMeshInstance]] + data_to_send = targetVTKWindow._Iren.filledROIPoints + + the_TCPserver.socket.waitForReadyRead(10) + + formatData = array(data_to_send,'d') + formatString = formatData.tostring() + if debug: print(formatData) + if debug: print(formatString) + if debug: print(len(formatString)) + the_TCPserver.socket.write(formatString) + + diff --git a/mrMeshPy/mp_VTKDrawing.py b/mrMeshPy/mp_VTKDrawing.py new file mode 100644 index 00000000..99b58c92 --- /dev/null +++ b/mrMeshPy/mp_VTKDrawing.py @@ -0,0 +1,254 @@ +#!/usr/bin/python + +''' +Commands that allow surface Drawing for meshes + +Currently the user enters "draw mode" which stops any further +zooming / rotation in the vtk window. Left clicks on the surface +make black dots at selected vertices. Right clicking closes the +loop and extracts a surface enclosed by the boundary. + +The user can toggle the draw mode on and off to re-enable +rotation or zooming. + +Currently the current mode / reset (to start a new ROI) is done +via the drop-down menu - #TODO we could enable key bindings. + +Andre' Gouws 2017 +''' + + +import vtk +from numpy import * + +debug = True + + + +def drawingPickPoint(obj, ev): + + if obj.inDrawMode == 1: + + try: #TODO will this capture non-pick events? + #get the pick position + if debug: print "started pick" + obj.GetPicker().Pick(obj.GetEventPosition()[0], obj.GetEventPosition()[1], 0, obj.GetRenderWindow().GetRenderers().GetFirstRenderer()) + if debug: print "ended pick" + if debug: print(obj.GetPicker().GetPointId()) + currPickPointID = obj.GetPicker().GetPointId() + obj.pickedPointIds.append(currPickPointID) # append to place holder for picked vtk point IDs so we can track + + ##draw on the selected vertices -just change each selected vertex for now -- we'll join the dots later + + # get the required surface + + currPolyData = obj.GetPicker().GetActor().GetMapper().GetInput() + if debug: print(currPolyData.GetPointData().GetScalars().GetValue(currPickPointID)) + + # append orig value to place holder for picked vtk point scalar values so we can revert + obj.pickedPointOrigValues.append(currPolyData.GetPointData().GetScalars().GetTuple(currPickPointID)) + + if debug: print(currPolyData.GetPointData().GetScalars().GetTuple(currPickPointID)) + + # now change the colour at the picked vertex + currPolyData.GetPointData().GetScalars().SetTuple(currPickPointID,(0,0,0,255)) #black for char array + if debug: print(currPolyData.GetPointData().GetScalars().GetTuple(currPickPointID)) + + # notify the stream that there has been a change + currPolyData.Modified() + obj.GetPicker().GetActor().GetMapper().Modified() + obj.Render() + + # handle non-pick events or append pick point to list for later drawing + if obj.GetPicker().GetPointId() != -1: #catch non-picks + obj.pickedPoints.InsertNextPoint(obj.GetPicker().GetPickPosition()[0], obj.GetPicker().GetPickPosition()[1], obj.GetPicker().GetPickPosition()[2]) + else: + print('no point here! .. skipping') + except: + print('ignored left mouse click - outside the object? - no point EXACTLY here? - Try again') + + else: + if debug: print('ignored left mouse click - not in draw mode') + pass + + + +def drawingMakeROI(obj, ev): + + if obj.inDrawMode == 1: + print('closing and filling ROI') + + # get the required surface + currPolyData = obj.GetPicker().GetActor().GetMapper().GetInput() + + # ... and points + currPoints = obj.pickedPoints + + # we're not actually going to change the original: we're going to clip out the region we want and + # show it as a separate surface "patch" - this allows easy control of transparency/colour etc + + # a vtk selection tool that loops around our points + if debug: print "before select polydata" + selecta = vtk.vtkSelectPolyData() + selecta.SetInputData(currPolyData) + selecta.SetLoop(currPoints) + selecta.GenerateSelectionScalarsOn() + ## selecta.SetSelectionModeToSmallestRegion() + selecta.SetSelectionModeToLargestRegion() + selecta.Update() + if debug: print "after select polydata" + + + ## TODO - consider a different clipping function + ## the one below appears to clip front AND back surfaces? + #loop = vtk.vtkImplicitSelectionLoop() + #loop.SetLoop(currPoints) + #clipROI = vtk.vtkExtractPolyDataGeometry() + #clipROI.SetInputData(currPolyData) + #clipROI.SetImplicitFunction(loop) + #clipROI.ExtractBoundaryCellsOn() + #clipROI.PassPointsOff() + #clipROI.Update() + + + # the tool that actually clips out the selected region + clipROI = vtk.vtkClipPolyData() + clipROI.SetInputConnection(selecta.GetOutputPort()) + clipROI.Update() + + + if debug: print(clipROI.GetOutput()) + + # --- try to get an ROI + if clipROI.GetOutput().GetNumberOfPoints() == 0: + #loop could not clip data for some reason - we need to clean up + obj.pickedPoints = vtk.vtkPoints() #new + #stop drawing mode and reset the interactor to normal behaviour + obj.inDrawMode = 0 + + print 'Error while trying to close the ROI loop - resetting colours - please try again.' + print obj.ScalarsCopyForRevert + + tmpScalarsCopyForRevert = vtk.vtkUnsignedCharArray() + tmpScalarsCopyForRevert.DeepCopy(obj.ScalarsCopyForRevert) + + obj.curr_polydata.GetPointData().SetScalars(obj.ScalarsCopyForRevert) + obj.curr_polydata.Modified() + + obj.ScalarsCopyForRevert = tmpScalarsCopyForRevert; #AND LOCK IN CURRENT STATE + + obj.curr_smoother.Update() + obj.curr_mapper.SetColorModeToDefault() + obj.curr_mapper.Modified() + + style = vtk.vtkInteractorStyleTrackballCamera() + obj.SetInteractorStyle(style) + + obj.ren.Render() + obj.Render() + + obj.parent_ui.statusbar.showMessage("Failed to close and fill ROI! Perhaps the surface is too spikey? - try smoothing.") + + else: + + ## -- track the point IDs + # unfortunately the clip tool does not keep a reference of the original surface + # point IDs (that we need to send bakc to other programs later + # -- we need a pretty nasty bit or code to track the xyz positions extracted + # so that we can go back to the original surface and label these + + + # first, get the xyz position of all points in original surface + origDat = currPolyData.GetPoints().GetData() + origArr = [] + + for i in range(currPolyData.GetNumberOfPoints()): + origArr.append(origDat.GetTuple(i)) + + + # next get the xyz position of select points + selectedDat = clipROI.GetOutput().GetPoints().GetData() + selectedArr = [] #place holder for loop + + if debug: print(selectedDat) + + for i in range(clipROI.GetOutput().GetNumberOfPoints()): + selectedArr.append(selectedDat.GetTuple(i)) #xyz + + ## now find matching indices in original using numpy manipulation + numpyOrigArr = array(origArr) + numpySelectedArr = array(selectedArr) + + #if debug: + # numpyOrigArr.tofile('orig.txt',',') + # numpySelectedArr.tofile('new.txt',',') + + # TODO - we only check the first 6 characters of the coordinate data - this may cause precision issues later - check + origAll = char.array(numpyOrigArr[:,0],'S6') + '-' + char.array(numpyOrigArr[:,1],'S6') + '-' + char.array(numpyOrigArr[:,2],'S6') + selectedAll = char.array(numpySelectedArr[:,0],'S6') + '-' + char.array(numpySelectedArr[:,1],'S6') + '-' + char.array(numpySelectedArr[:,2],'S6') + + # and pipe it back to the original + obj.filledROIPoints = where(in1d(origAll, selectedAll))[0] + obj.ROI_ready = 1 + + if debug: print(obj.filledROIPoints) + if debug: print(size(obj.filledROIPoints)) + + + ## -- back to drawing + # draw the new surface + roiMapper = vtk.vtkPolyDataMapper() + roiMapper.SetInputConnection(clipROI.GetOutputPort()) + roiMapper.ScalarVisibilityOff() + roiMapper.Update() + + roiActor = vtk.vtkActor() + roiActor.SetVisibility(1) + roiActor.SetMapper(roiMapper) + roiActor.GetProperty().SetOpacity(0.5) + roiActor.GetProperty().SetColor(1.0,1.0,1.0) + + obj.roiActor = roiActor; + obj.ren.AddActor(obj.roiActor) + obj.Render() + + #stop drawing mode and reset hte interactor to normal behaviour + obj.inDrawMode = 0 + + #remove our black dots + tmpScalarsCopyForRevert = vtk.vtkUnsignedCharArray() + tmpScalarsCopyForRevert.DeepCopy(obj.ScalarsCopyForRevert) + + obj.curr_polydata.GetPointData().SetScalars(obj.ScalarsCopyForRevert) + obj.curr_polydata.Modified() + + obj.ScalarsCopyForRevert = tmpScalarsCopyForRevert; #AND LOCK IN CURRENT STATE + + obj.curr_smoother.Update() + obj.curr_mapper.SetColorModeToDefault() + obj.curr_mapper.Modified() + + style = vtk.vtkInteractorStyleTrackballCamera() + obj.SetInteractorStyle(style) + + # and clean up for the next loop? TODO? - or button to reset? + obj.pickedPoints = vtk.vtkPoints() #new + + else: + if debug: print('ignored right mouse click - not in draw mode') + pass + + + + + + + + + + + + + + + diff --git a/mrMeshPy/mp_VTKProcessing.py b/mrMeshPy/mp_VTKProcessing.py new file mode 100644 index 00000000..9d3f1a6d --- /dev/null +++ b/mrMeshPy/mp_VTKProcessing.py @@ -0,0 +1,149 @@ +## TODO header + +import vtk +from numpy import pi + +debug = True + + +def VTK_smoothing(the_smoother, the_mapper, iterations, relaxation_factor): + + # Standard way to perform mesh smoothing via relaxation in VTK. + # I'm pretty sure this is a direct replication of what mrMesh was + # doing before. + + if debug: print 'starting smoothing' + the_smoother.SetNumberOfIterations(iterations) + the_smoother.SetRelaxationFactor(relaxation_factor) + the_smoother.Modified() + the_smoother.Update() + + if debug: print 'finished smoothing -setting up mapper' + the_mapper.SetInputConnection(the_smoother.GetOutputPort()) + the_mapper.SetScalarModeToUsePointData() + the_mapper.SetColorModeToDefault() + the_mapper.Modified() + + if debug: print 'finished mapper -setting up actor' + newActor = vtk.vtkActor() + newActor.SetMapper(the_mapper) + if debug: print 'finished actor' + + return newActor + + +def VTK_updateMesh(currVTKInstance, colorData, mainWindowUI): + + # user has send some new scalar values to be rendered on the mesh BUT COLORS ARE + # ALREADY CALCULATED IN MATLAB + + if debug: print(colorData) + + currVTKInstance.curr_polydata.GetPointData().SetScalars(colorData) + currVTKInstance.curr_polydata.Modified() + + currVTKInstance.curr_smoother.Update() + currVTKInstance.curr_mapper.SetColorModeToDefault() + currVTKInstance.curr_mapper.Modified() + + # in case of error when drawing ROIs we can revert the color map + # turns out that later processes access the inherited renderwindowinteractor (?) + # so lets put all the above in the scope of that + currVTKInstance._Iren.ScalarsCopyForRevert = vtk.vtkUnsignedCharArray() + currVTKInstance._Iren.ScalarsCopyForRevert.DeepCopy(colorData) + + newActor = vtk.vtkActor() + newActor.SetMapper(currVTKInstance.curr_mapper) + + return newActor + + + if debug: print('colorData processed in VTK_updateMeshDirect') + + + +''' OBSOLETE if we use the direct color import - i.e. all colour handling done in matlab +def VTK_updateMesh(currVTKInstance, newLUT, phaseData, cohData, cohThr, mainWindowUI): + + # user has send some new scalar values to be rendered on the mesh + + ## TODO - other data types - assuming phase for now + # rescale phase data into range: 0-1024 to match lookup table + rescaledPhaseData = phaseData/(2.0*pi)*1024.0 + + curvature = currVTKInstance.curr_curvature + + scalars = vtk.vtkFloatArray() + + for i in range(len(rescaledPhaseData)): + if cohData[i] >= float(cohThr): # apply threshold + scalars.InsertNextValue(rescaledPhaseData[i]) + else: + if curvature[i] < 0: + scalars.InsertNextValue(1024+50) + else: + scalars.InsertNextValue(1024+150) + + currVTKInstance.curr_polydata.GetPointData().SetScalars(scalars) + currVTKInstance.curr_polydata.Modified() + + currVTKInstance.curr_smoother.Update() + + currVTKInstance.curr_mapper.SetLookupTable(newLUT) + currVTKInstance.curr_mapper.SetScalarRange(0,1224) + currVTKInstance.curr_mapper.Modified() + + newActor = vtk.vtkActor() + newActor.SetMapper(currVTKInstance.curr_mapper) + + return newActor + + + if debug: print('here') + + + +def VTK_buildLookupTable(r_vec, g_vec, b_vec): + + # we need to build custom colour lookup table that can show colour + # data for vertices above threshold and grayscale anatomy at + # vertices that do nor reach threshold. + # I build a lookup table that is made up of two parts: the lower + # end of the table is a RGB colour map (entries 1-1024)- when we get + # the scalar values from vista for the overlay data (e.g. phase) we + # rescale the incoming data into the range 1-1024 so that the scalar + # values map explicitly onto a colour table value. + # The 'upper' part of the table (1025-1224) has grayscale values. In + # theory we could just have 2 extra entries in the table for light + # or dark gray, but interpolation in vtk can cause some odd effects + # and color blending so we pad and extra 200 values onto the table + + + #create an arbitrary LUT with 1224 entries, we'll overwrite this + cLUT = vtk.vtkLookupTable() + cLUT.SetHueRange(0,1) + cLUT.SetValueRange(1,1) + cLUT.SetSaturationRange(1,1) + cLUT.SetNumberOfColors(1224) + cLUT.Build() #build it + + # now overwrite it with the incoming RGB data for entries 1-1024 + for i in range(1024): + cLUT.SetTableValue(i,(r_vec[i],g_vec[i],b_vec[i],1)) + + # and overwrite 1025-1124 with dark gray + for i in range(1024,1124): + cLUT.SetTableValue(i,(0.3,0.3,0.3,1)) + + # and the overwrite 1125-1224 with light gray + for i in range(1124,1224): + cLUT.SetTableValue(i,(0.7,0.7,0.7,1)) + + ## TODO - allow modulation of light/dark gray levels? + + # rebuild the table with these new values + cLUT.Modified() + + + return cLUT +''' diff --git a/mrMeshPy/mp_VTKRoutines.py b/mrMeshPy/mp_VTKRoutines.py new file mode 100644 index 00000000..5989ab50 --- /dev/null +++ b/mrMeshPy/mp_VTKRoutines.py @@ -0,0 +1,430 @@ +#!/usr/bin/python + +''' +VTK engine room for mrMeshPy viewer + +The main vtk processing is done by functions here - although some hardcore +processing is handled in subroutines of other imported modules. + +A core concept here is the tracking (kepping in scope) or the "targetVTKWindow" + - this is a vtkRenderWindowInteractor instance in the main program UI (user +interface) - by creatoing multiple instances of vtk windows we can load +multiple meshes. Some functions reference this specifically with a reference +index passed from mrVista --- mainWindowUI.vtkInstances[int(theMeshInstance)] +while others just referene the most recently added instance (e.g. when adding +a new mesh) --- mainWindowUI.vtkInstances[-1] + +Note that it is the mainWindowUI that is passed to all functions so that all +funcitons have the content of the main window in scope. + +Andre' Gouws 2017 +''' + + +import vtk +from numpy import * +import time + +from vtk.util import numpy_support + +debug = True + +# local modules +from mp_unpackIncomingData import unpackData +from mp_VTKProcessing import * +from mp_VTKDrawing import * + + + +def loadNewMesh(currVTKInstance, commandArgs, mainWindowUI, the_TCPserver): + + #first get all the data we are expecting from the server + ## NB this assumes that the order of sending by the server is + # 1) vertices + # 2) triangles + # 3) color data r (rgba) for each vertex + # 4) color data g (rgba) for each vertex + # 5) color data b (rgba) for each vertex + # 6) color data a (rgba) for each vertex + + + if debug: + print('received request for new mesh with Args:') + print(commandArgs) + + # sanity check + if ('vertices' in commandArgs[0]) and ('triangles' in commandArgs[1]): + pass + else: + return "error - expecting vertices, then triangles!" + + + # load the surfaces data + verticesArgs = commandArgs[0].strip().split(',') + vertices = unpackData(verticesArgs[1], int(verticesArgs[2]), the_TCPserver) + vertices = array(vertices,'f') + vertices = vertices.reshape((len(vertices)/3,3)) + + trianglesArgs = commandArgs[1].strip().split(',') + triangles = unpackData(trianglesArgs[1], int(trianglesArgs[2]), the_TCPserver) + triangles = array(triangles,'f') + if debug: print(triangles) + triangles = triangles.reshape((len(triangles)/3,3)) + if debug: print(triangles) + + # load the surface colour data + rVecArgs = commandArgs[2].strip().split(',') + r_vec = unpackData(rVecArgs[1], int(rVecArgs[2]), the_TCPserver) + r_vec = array(r_vec,'uint8') + if debug: print(r_vec) + + gVecArgs = commandArgs[3].strip().split(',') + g_vec = unpackData(gVecArgs[1], int(gVecArgs[2]), the_TCPserver) + g_vec = array(g_vec,'uint8') + + bVecArgs = commandArgs[4].strip().split(',') + b_vec = unpackData(bVecArgs[1], int(bVecArgs[2]), the_TCPserver) + b_vec = array(b_vec,'uint8') + + aVecArgs = commandArgs[5].strip().split(',') + a_vec = unpackData(aVecArgs[1], int(aVecArgs[2]), the_TCPserver) + a_vec = array(a_vec,'uint8') + + if debug: + print(len(r_vec)) + print(len(g_vec)) + print(len(b_vec)) + print(len(a_vec)) + + #combine into numpy array + colorDat = squeeze(array(squeeze([r_vec,g_vec,b_vec,a_vec]),'B',order='F').transpose()) + + # convert this to a VTK unsigned char array + scalars = numpy_support.numpy_to_vtk(colorDat,0) + + curr_scalars = vtk.vtkUnsignedCharArray() + curr_scalars.DeepCopy(scalars) + + ## ---- ok, we hav the data, lets turn it into vtk stuff + + # Process vertices + points = vtk.vtkPoints() + for i in range(vertices.shape[0]): + points.InsertPoint(i,vertices[i][0],vertices[i][1],vertices[i][2]) + + # Process faces (triangles) + polys = vtk.vtkCellArray() + + + nTriangles = triangles.shape[0] + for i in range(nTriangles): + polys.InsertNextCell(3) + for j in range(3): + polys.InsertCellPoint(int(triangles[i][j])) + + # check + if debug: print(points) + if debug: print(polys) + if debug: print(scalars) + if debug: print(currVTKInstance) + + # Assemble as PolyData + polyData = vtk.vtkPolyData() + polyData.SetPoints(points) + polyData.SetPolys(polys) + polyData.GetPointData().SetScalars(scalars) + + ## TODO ? smoothing on first load? + + smooth = vtk.vtkSmoothPolyDataFilter() + smooth = vtk.vtkSmoothPolyDataFilter() + smooth.SetNumberOfIterations(0) + smooth.SetRelaxationFactor(0.0) + smooth.FeatureEdgeSmoothingOff() + smooth.SetInputData(polyData) + + pdm = vtk.vtkPolyDataMapper() + pdm.SetScalarModeToUsePointData() + pdm.SetInputConnection(smooth.GetOutputPort()) + + actor = vtk.vtkActor() + actor.SetMapper(pdm) + + iren = mainWindowUI.vtkInstances[-1] + + ## ---- engine room for drawing on the surface + + # add a picker that allows is top pick points on the surface + picker = vtk.vtkCellPicker() + picker.SetTolerance(0.0001) + mainWindowUI.vtkInstances[-1].SetPicker(picker) + mainWindowUI.vtkInstances[-1]._Iren.pickedPointIds = [] #place holder for picked vtk point IDs so we can track + mainWindowUI.vtkInstances[-1].pickedPointIds = mainWindowUI.vtkInstances[-1]._Iren.pickedPointIds + mainWindowUI.vtkInstances[-1]._Iren.pickedPointOrigValues = [] #place holder for picked vtk point IDs so we can track + mainWindowUI.vtkInstances[-1].pickedPointOrigValues = mainWindowUI.vtkInstances[-1]._Iren.pickedPointOrigValues + mainWindowUI.vtkInstances[-1]._Iren.pickedPoints = vtk.vtkPoints() #place holder for picked vtk point IDs so we can track + mainWindowUI.vtkInstances[-1].pickedPoints = mainWindowUI.vtkInstances[-1]._Iren.pickedPoints + mainWindowUI.vtkInstances[-1]._Iren.inDrawMode = 0 #TODO + mainWindowUI.vtkInstances[-1].inDrawMode = mainWindowUI.vtkInstances[-1]._Iren.inDrawMode + + # drawing functions imported from mp_VTKDrawing + mainWindowUI.vtkInstances[-1].AddObserver('LeftButtonPressEvent', drawingPickPoint, 1.0) + mainWindowUI.vtkInstances[-1].AddObserver('RightButtonPressEvent', drawingMakeROI, 1.0) + + ren = mainWindowUI.vtkInstances[-1].ren + mainWindowUI.vtkInstances[-1]._Iren.ren = ren + + # ADD A LIGHT SOURCE TODO: MAKE THIS OPTIONAL/DEFAULT? + lightKit = vtk.vtkLightKit() + lightKit.SetKeyLightIntensity(0.5) + # TODO: SOME OPTIONS TO EXPLORE + #lightKit.MaintainLuminanceOn() + #lightKit.SetKeyLightIntensity(1.0) + ## warmth of the lights + #lightKit.SetKeyLightWarmth(0.65) + #lightKit.SetFillLightWarmth(0.6) + #lightKit.SetHeadLightWarmth(0.45) + ## intensity ratios + ## back lights will be very dimm + lightKit.SetKeyToFillRatio(1.) + lightKit.SetKeyToHeadRatio(2.) + lightKit.SetKeyToBackRatio(1.) + lightKit.AddLightsToRenderer(ren) + + ren.AddActor(actor) + ren.SetBackground(1,1,1) + ren.ResetCamera() + ren.Render() + mainWindowUI.vtkInstances[-1].Render() + + # lets put some of the data objects in the scope of the + # main window so that they can be manipulated later. + mainWindowUI.vtkInstances[-1].curr_actor = actor + mainWindowUI.vtkInstances[-1].curr_smoother = smooth + mainWindowUI.vtkInstances[-1].curr_polydata = polyData + mainWindowUI.vtkInstances[-1].curr_mapper = pdm + mainWindowUI.vtkInstances[-1].curr_camera = ren.GetActiveCamera() + # and the raw mesh coordinate data.. why not + mainWindowUI.vtkInstances[-1].curr_points = points + mainWindowUI.vtkInstances[-1].curr_polys = polys + mainWindowUI.vtkInstances[-1].curr_scalars = curr_scalars #Deep copied + + + # turns out that later processes access the inherited renderwindowinteractor (?) + # so lets put all the above in the scope of that too + mainWindowUI.vtkInstances[-1]._Iren.curr_actor = actor + mainWindowUI.vtkInstances[-1]._Iren.curr_smoother = smooth + mainWindowUI.vtkInstances[-1]._Iren.curr_polydata = polyData + mainWindowUI.vtkInstances[-1]._Iren.curr_mapper = pdm + mainWindowUI.vtkInstances[-1]._Iren.curr_camera = ren.GetActiveCamera() + mainWindowUI.vtkInstances[-1]._Iren.curr_points = points + mainWindowUI.vtkInstances[-1]._Iren.curr_polys = polys + mainWindowUI.vtkInstances[-1]._Iren.curr_scalars = curr_scalars #Deep copied + + # and so we can access ui controls (e.g. statusbar) from the inherited window + mainWindowUI.vtkInstances[-1]._Iren.parent_ui = mainWindowUI + + + def KeyPress(obj, evt): + key = obj.GetKeySym() + if key == 'l': + currVTKinstance = len(mainWindowUI.vtkInstances) + print(key) + print(mainWindowUI.vtkInstances[currVTKinstance-1]) + + + #let's also track key presses per instance esp for the draw routine :) + mainWindowUI.vtkInstances[-1].AddObserver("KeyPressEvent",KeyPress) + + mainWindowUI.tabWidget.setCurrentIndex(len(mainWindowUI.vtkInstances)-1) #zero index + + + +def smoothMesh(theMeshInstance, commandArgs, mainWindowUI, the_TCPserver): + + #lets try to get the apt window + try: + targetVTKWindow = mainWindowUI.vtkInstances[mainWindowUI.vtkDict[theMeshInstance]] + except: + print ('No mesh instance with id:%s currently available - may need a re-synch' %theMeshInstance) + #return error + return 1 + + # lets show the correct tab + mainWindowUI.tabWidget.setCurrentIndex(int(mainWindowUI.vtkDict[theMeshInstance])) + #mainWindowUI.tabWidget.repaint() + mainWindowUI.tabWidget.update() + + + #lets get the original data + the_smoother = targetVTKWindow.curr_smoother + the_mapper = targetVTKWindow.curr_mapper + + + if debug: print(targetVTKWindow.curr_actor.GetMapper().GetInput().GetPointData().GetScalars()) + if debug: print(targetVTKWindow.curr_actor.GetMapper().GetInput().GetPointData().GetScalars().GetTuple(1000)) + + #expecting a string that reads something like 'iterations,200,relaxationfactor,1.2' + # sanity check + if ('iterations' in commandArgs[0]) and ('relaxationfactor' in commandArgs[0]): + smoothingArgs = commandArgs[0].strip().split(',') + iterations = int(smoothingArgs[1]) + relaxationfactor = float(smoothingArgs[3]) + else: + return "error - expecting vertices, then curvature, then triangles!" + + if debug: print 'starting smoothing callback' + newActor = VTK_smoothing(the_smoother, the_mapper, iterations, relaxationfactor) + + if debug: print 'smoothing callback returned new actor' + + if debug: print 'removing old actor' + targetVTKWindow.ren.RemoveActor(targetVTKWindow.curr_actor) + if debug: print 'adding new actor' + targetVTKWindow.ren.AddActor(newActor) + + if debug: print 'added new actor - changing curr actor pointer' + targetVTKWindow.curr_actor = newActor #lets keep track + + if debug: print 'trying to update ' + # run mesh update to reset the color map (smoothing "messes" this up) + updateMeshData(theMeshInstance, [], mainWindowUI, the_TCPserver) + if debug: print 'update completed' + + #return success + return 0 + + +def updateMeshData(theMeshInstance, commandArgs, mainWindowUI, the_TCPserver): + + # here the base mesh is already loaded and we are simply updating with the + # current View settings in from the vista session WITH THE COLOR VALUES FROM + # VISTA - i.e. do not go through a lookuptable + + #lets try to get the apt window + try: + targetVTKWindow = mainWindowUI.vtkInstances[mainWindowUI.vtkDict[theMeshInstance]] + except: + print ('No mesh instance with id:%s currently available - may need a re-synch' %theMeshInstance) + #return error + return 1 + + + # lets show the correct tab + mainWindowUI.tabWidget.setCurrentIndex(int(mainWindowUI.vtkDict[theMeshInstance])) #zero index + #mainWindowUI.tabWidget.repaint() + mainWindowUI.tabWidget.update() + + + #lets get the original data + the_polyData = targetVTKWindow.curr_polydata + the_mapper = targetVTKWindow.curr_mapper + + #first get all the data we are expecting from the server + ## NB this assumes that the order of sending by the server is + # 1) r_vector - red component + # 2) g_vector - blue component + # 3) b_vector - green component + # 4) a_vector - aplha component + + + if debug: + print('received request for UPDATE DIRECT mesh with Args:') + print(commandArgs) + + if len(commandArgs) != 0 : #new data has come from MATLAB so recompute + + # load the surfaces data + rVecArgs = commandArgs[0].strip().split(',') + r_vec = unpackData(rVecArgs[1], int(rVecArgs[2]), the_TCPserver) + r_vec = array(r_vec,'uint8') + if debug: print(r_vec) + + gVecArgs = commandArgs[1].strip().split(',') + g_vec = unpackData(gVecArgs[1], int(gVecArgs[2]), the_TCPserver) + g_vec = array(g_vec,'uint8') + + bVecArgs = commandArgs[2].strip().split(',') + b_vec = unpackData(bVecArgs[1], int(bVecArgs[2]), the_TCPserver) + b_vec = array(b_vec,'uint8') + + aVecArgs = commandArgs[3].strip().split(',') + a_vec = unpackData(aVecArgs[1], int(aVecArgs[2]), the_TCPserver) + a_vec = array(a_vec,'uint8') + + if debug: + print(len(r_vec)) + print(len(g_vec)) + print(len(b_vec)) + print(len(a_vec)) + + #combine into numpy array + colorDat = squeeze(array(squeeze([r_vec,g_vec,b_vec,a_vec]),'B',order='F').transpose()) + + # convert this to a VTK unsigned char array + vtkColorArray = numpy_support.numpy_to_vtk(colorDat,0) + + # keep a "deep" copy - this is to workaround some artifacts generated + # by vtk algorithms (e.g. smoothing) that also smooth the color data + # on the surface and then automatically update the inherited color map + # - we allow vtk to do this but then overwrite the recomptued color + # map AFTER the algorithms have run + + deepCopyScalars = vtk.vtkUnsignedCharArray() + deepCopyScalars.DeepCopy(vtkColorArray) + targetVTKWindow.curr_scalars = deepCopyScalars + + #TODO - this may have impact on later processing - investigate + + else: + # no new data from MATLAB, probably just an internal re-draw call + # after something like smoothing - just grab the current deep + # copy of the required scalars + vtkColorArray = targetVTKWindow.curr_scalars + + + # OK - we have the data - let's update the mesh + newActor = VTK_updateMesh(targetVTKWindow, vtkColorArray, mainWindowUI) + + targetVTKWindow.ren.AddActor(newActor) + targetVTKWindow.ren.RemoveActor(targetVTKWindow.curr_actor) + + targetVTKWindow.curr_actor = newActor #lets keep track + targetVTKWindow.ren.Render() + targetVTKWindow.Render() + print('success with direct mesh update routine') + + #return success + return 0 + + + +## -------------------------------------------------------------------------------- +# test example animation +def rotateMeshAnimation(currVTKInstance, commandArgs, mainWindowUI, the_TCPserver): + + #rotation args + rotations = commandArgs[0].strip().split(',') + rotations = unpackData(rotations[1], int(rotations[2]), the_TCPserver) + + if debug: print(rotations) + + targetVTKWindow = mainWindowUI.vtkInstances[int(currVTKInstance)] #NB zero indexing + + camera = targetVTKWindow.ren.GetActiveCamera() + if debug: print(camera) + + for i in range(len(rotations)): + camera.Azimuth(rotations[i]) + #targetVTKWindow.ren.Render() + targetVTKWindow.iren.Render() + + time.sleep(0.02) + + the_TCPserver.socket.write(str('send useful message back here TODO')) + +## -------------------------------------------------------------------------------- + + + + diff --git a/mrMeshPy/mp_setupVTKWindow.py b/mrMeshPy/mp_setupVTKWindow.py new file mode 100644 index 00000000..7cb4070e --- /dev/null +++ b/mrMeshPy/mp_setupVTKWindow.py @@ -0,0 +1,110 @@ +#!/usr/bin/python + +''' +This module sets up a VTK window instance for each new +mesh instance sent from mrVista + +Andre' Gouws 2017 +''' + + +from PyQt5 import QtCore, QtGui, QtNetwork, QtWidgets +from QVTKRenderWindowInteractor import QVTKRenderWindowInteractor + +import vtk + +debug = False + + +# add a new VTK window as an extra tab in the Main window's tab widget +def mrMeshVTKWindow(parentUI, theMeshInstance, data): + _translate = QtCore.QCoreApplication.translate + + #create a new vtkWidget, appending to the list of exisitng widgets + if debug: print(parentUI) + newVTKWidgetInstance = QVTKRenderWindowInteractor(parentUI.centralwidget) + parentUI.vtkInstances.append(newVTKWidgetInstance) + + currMeshCount = (len(parentUI.vtkInstances)) + + # we have to make a new tab and layout for it to go to. + parentUI.newTab = QtWidgets.QWidget() + sizePolicy = QtWidgets.QSizePolicy(QtWidgets.QSizePolicy.MinimumExpanding, QtWidgets.QSizePolicy.MinimumExpanding) + sizePolicy.setHorizontalStretch(0) + sizePolicy.setVerticalStretch(0) + sizePolicy.setHeightForWidth(parentUI.newTab.sizePolicy().hasHeightForWidth()) + parentUI.newTab.setSizePolicy(sizePolicy) + parentUI.newTab.setObjectName("tab%s" %currMeshCount) + parentUI.gridLayoutTabNew = QtWidgets.QGridLayout(parentUI.newTab) + parentUI.gridLayoutTabNew.setContentsMargins(0, 0, 0, 0) + parentUI.gridLayoutTabNew.setObjectName("gridLayoutTab%s" %currMeshCount) + parentUI.gridLayoutVTKWinNew = QtWidgets.QGridLayout() + parentUI.gridLayoutVTKWinNew.setObjectName("gridLayoutVTKWin%s" %currMeshCount) + parentUI.gridLayoutTabNew.addLayout(parentUI.gridLayoutVTKWinNew, 0, 0, 1, 1) + + parentUI.gridLayoutVTKWinNew.addWidget(newVTKWidgetInstance, 0, 0, 1, 1) + parentUI.tabWidget.addTab(parentUI.newTab, "Mesh-%s" %theMeshInstance) + + # have to set the resize policy now + sizePolicy = QtWidgets.QSizePolicy(QtWidgets.QSizePolicy.MinimumExpanding, QtWidgets.QSizePolicy.MinimumExpanding) + parentUI.vtkInstances[-1].setSizePolicy(sizePolicy) + parentUI.vtkInstances[-1].ren = vtk.vtkRenderer() + parentUI.vtkInstances[-1].GetRenderWindow().AddRenderer(parentUI.vtkInstances[-1].ren) + parentUI.vtkInstances[-1].iren = parentUI.vtkInstances[-1].GetRenderWindow().GetInteractor() + parentUI.gridLayoutVTKWinNew.addWidget(parentUI.vtkInstances[-1], 0, 0, 1, 1) + + parentUI.tabWidget.setCurrentIndex(len(parentUI.vtkInstances)-1) #zero index + parentUI.tabWidget.update() + + style = vtk.vtkInteractorStyleTrackballCamera() + parentUI.vtkInstances[-1].SetInteractorStyle(style) + + # flip the camera + parentUI.vtkInstances[-1].ren.GetActiveCamera().SetViewUp(0,-1,0) + + if data == 'debug': + loadTestVTKWindow(parentUI) + else: + initializeEmptyVTKWindowInstance(parentUI) + print('Finished setting up new VTK window instance ...') + + + +# just set up an empty window +def initializeEmptyVTKWindowInstance(parentUI): + parentUI.vtkInstances[-1].ren.ResetCamera() + parentUI.vtkInstances[-1].iren.Start() + + +#debug testing content for VTK window instance +def loadTestVTKWindow(parentUI): + + cube = vtk.vtkCubeSource() + cube.SetXLength(200) + cube.SetYLength(200) + cube.SetZLength(200) + cube.Update() + cm = vtk.vtkPolyDataMapper() + cm.SetInputConnection(cube.GetOutputPort()) + ca = vtk.vtkActor() + ca.SetMapper(cm) + parentUI.vtkInstances[-1].ren.AddActor(ca) + + axesActor = vtk.vtkAnnotatedCubeActor(); + axesActor.SetXPlusFaceText('R') + axesActor.SetXMinusFaceText('L') + axesActor.SetYMinusFaceText('H') + axesActor.SetYPlusFaceText('F') + axesActor.SetZMinusFaceText('P') + axesActor.SetZPlusFaceText('A') + axesActor.GetTextEdgesProperty().SetColor(1,1,0) + axesActor.GetTextEdgesProperty().SetLineWidth(2) + axesActor.GetCubeProperty().SetColor(0,0,1) + axes = vtk.vtkOrientationMarkerWidget() + axes.SetOrientationMarker(axesActor) + axes.SetInteractor(parentUI.vtkInstances[-1].iren) + + parentUI.vtkInstances[-1].ren.ResetCamera() + parentUI.vtkInstances[-1].iren.Start() + + diff --git a/mrMeshPy/mp_unpackIncomingData.py b/mrMeshPy/mp_unpackIncomingData.py new file mode 100644 index 00000000..709082b0 --- /dev/null +++ b/mrMeshPy/mp_unpackIncomingData.py @@ -0,0 +1,96 @@ +#!/usr/bin/python + +''' +This module interprets the commands sent from mrVista output to +mrMeshPy (python) through the mrMeshPyServer .. + +Matlab sends a string in one transaction giving a command to the +visualisation module. This command either performs an explicit +function in the viewer (e.g. rotate the camera 90 degrees) or the +commnd describes the configuration/content of a large data chunk to +will be sent in the subsequent transaction so that we know how to +unpack the data, and how to process it (e.g. 70,000 floating point +numbers which are scalar values to show as an amplitude map. + +Each command string is interpreted by the pm_commandInterpret module. + +N.B. - currently command strings have a maximum length of 1024 bytes. + +Commands are specifically ordered, semi-colon seperated strings which are +unpacked to describe what the user is trying to do / send from matlab. +Commands have a MINIMUM LENGTH of 6 arguments and have the following +structure and item order (zero-indexed) + +0 - "cmd" -- always this, identifies it as a cmd :) +1 - "drawOnly" or "newData" - allows an immediate action or tells + the program the content of the next TCP transaction +2 - "None" if no data is incoming, "i" integer, "d" double", "f" float etc +3 - integer -- the number of samples we are sending +4 - commandName - should match a command in mp_Commands file +5 - theMeshInstance - integer pointing to the the mesh window that we + want to operate on +6 onwards - commandArgs - a list of comma-separated pairs of arguments + to characterise the processing of the incoming data + blob or apply some settings to the viewport - + CAN BE EMPTY but must be set to [] + + +Andre' Gouws 2017 +''' + +import struct + +#local modules +from mp_Commands import * + + +debug = False + + + +def unpackData(incomingDataType, incomingNumberOfSamples, the_TCPserver): + + try: + # lets determine the size and nature of the incoming data + + # .. and how many bytes will that be? + incomingBytes = struct.calcsize(incomingDataType)*incomingNumberOfSamples + + # report + print('Expecting %i bytes from %i samples of %s dataType ... waiting ...' %(incomingBytes, incomingNumberOfSamples, incomingDataType)) + + ##### Read number of data bytes ## TODO revise the weird behaviour here + ####the_TCPserver.socket.waitForReadyRead() + + # report started + print("socket open") + + # a placeholder that we append the incoming string to + dataReceived = '' + + # lets keep track of how many bytes we have actually received + bytesReceived = 0 + + # .. start reading ## TODO - 10ms timeout here OK? + while bytesReceived < incomingBytes: + the_TCPserver.socket.waitForReadyRead(10) + if debug: print('bytesReceived', bytesReceived) + tmp = the_TCPserver.socket.read(incomingBytes-bytesReceived) + bytesReceived += len(tmp) + dataReceived += tmp + + print('read completed ... ') + + # now unpack the data to + dataObject = struct.unpack(incomingNumberOfSamples*incomingDataType, dataReceived) + print('unpacked %i data points.. ' %incomingNumberOfSamples) + print('done loading incoming data!') + + if debug: print(dataObject[-10]) + + return dataObject + + except Exception as e: + print("Trouble unpacking data") + + diff --git a/mrMeshPy/mrMeshPy.py b/mrMeshPy/mrMeshPy.py new file mode 100644 index 00000000..a082b313 --- /dev/null +++ b/mrMeshPy/mrMeshPy.py @@ -0,0 +1,94 @@ +#!/usr/bin/python + +''' +The main mrMeshPy qapplication window built with Qt5 + +Andre' Gouws 2017 +''' + + +import sys +import vtk + +# some local modules +from mrMeshPyQtTCPServer import mrMeshPyQtTCPServer +from mp_MainWindow import Ui_MrMeshMainWindow +from mp_MenuCallbacks import Ui_setupMenuBarCallbacks +from mp_setupVTKWindow import mrMeshVTKWindow + +# for the gui +from PyQt5 import QtCore, QtGui, QtNetwork, QtWidgets + + +#for debugging/testing +debug = False +test = False + + +class mrMeshPyMainWindow(QtWidgets.QMainWindow): + + def __init__(self, parent = None): + + # set up the initial user interface window on Qt + QtWidgets.QMainWindow.__init__(self, parent) + self.ui = Ui_MrMeshMainWindow() + self.ui.setupUi(self) + + + # set up a tcp server instance (Qt TCP server) + self.server = mrMeshPyQtTCPServer(self) + self.ui.TCP_server = self.server ## this puts the service instance into the scope of the user interface (ui) + + + # TODO - hard code TCP port 9999 for now #next available may be better + if not self.server.listen(QtNetwork.QHostAddress('127.0.0.1'), 9999): + print('Error starting TCP instance on requested port - quitting.') + return + print('Running TCP instance on port %d' % self.server.serverPort()) + + + # set up the menu bar in the UI window with its callbacks + Ui_setupMenuBarCallbacks(self.ui) + + + # lets keep track of how many vtk instances are open in a list that is in the scope of the ui + self.ui.vtkInstances = [] #empty for now + # and we are going to use a dictionary to keep track of unique instances + # the unique string generated by matlab will point at an index to the vtkInstances list #TODO is this sensible? + self.ui.vtkDict = {} + + + self.ui.statusbar.showMessage(' ... Waiting for matlab to send a mesh ...') + + + + # ---Ready to go unless testing/debugging --------------------------------- + + # some debug options + if debug: + # --------------- + # add a vtk window! + mrMeshVTKWindow(self.ui, 'debug') + self.ui.statusbar.showMessage(' ... in debug mode; loaded test VTK window ...') + + def testMenuPrint(self): + print('test menu print message') + + + # some other test options + if test: + from testRender import mp_launchVTKWindow + mp_launchVTKWindow(self.ui) + + # ------------------------------------------------------------------------- + + + + +if __name__ == '__main__': + + app = QtWidgets.QApplication(sys.argv) + window = mrMeshPyMainWindow() + window.show() + sys.exit(app.exec_()) + diff --git a/mrMeshPy/mrMeshPyQtTCPServer.py b/mrMeshPy/mrMeshPyQtTCPServer.py new file mode 100644 index 00000000..09839df1 --- /dev/null +++ b/mrMeshPy/mrMeshPyQtTCPServer.py @@ -0,0 +1,123 @@ + #!/usr/bin/python + +''' +A TCP server in QT to handle mrVista throughput to mrMeshPy + +Mark Hymers & Andre' Gouws 2017 +''' + + +from PyQt5 import QtCore, QtGui, QtNetwork, QtWidgets + +#local modules +from mp_Commands import run_mp_command + +debug = True + +class mrMeshPyQtTCPServer(QtNetwork.QTcpServer): + def __init__(self, parent): + super(mrMeshPyQtTCPServer, self).__init__(parent) + + # start up the server checking for incoming data + self.newConnection.connect(self.startReceiveCommand) + + # place the parent object (mrMeshPyMainWindow) in the scope of the server + self.mainWindow = parent + + # set up a counter for testing later #TODO remove? + self.counter = 0 + + + + def startReceiveCommand(self): + + self.mainWindow.ui.lockedForProcessing = 1 #TODO remove? + + self.socket = self.nextPendingConnection() + print("Server waiting for incoming data ...") + # Wait until we've read the command + # TODO Add timeout? + self.socket.waitForReadyRead(1000) + + # Parse command #TODO we assume a command will never be longer than 1024 bytes + command = self.socket.read(1024).strip() + + # commands have multiple description strings that are ;separeated - unpack here + incomingData = command.strip().split(';') + + # for testing + if debug: + print('incomingData') + print(incomingData) + + # lets start unpacking the incoming data packet + # we expect anything coming in to first have a command (cmd) + # which is then followed by some specifically sized data blobs + + if incomingData[0] != 'cmd': + #something has gone wrong - stop here + print('error! - command string received but not recognised:') + + else: + # we got a command + print('Full incoming command is:') + print(incomingData) + + # fields 2,3,and 4 (index 1,2,3) in the command string are placeholders -- ignore for now #TODO + + # the 5th field (zero index = 4) is the actual, sensibly-named command string e.g. paint_my_car + commandName = incomingData[4] + + # field 6 (index 5) -- the incoming command should always reference a particular mesh instance + # Major change here - no longer an index - will be a unique string identifier related to clock + theMeshInstance = incomingData[5] + + + # field 7 onwards (may be multiple) are extra arguments specific to a particluar command + + # now check if there are any extra arguments (e.g. smoothing. threshold etc) and add to a list + commandArgs = [] #placeholder + for i in range(6,len(incomingData)): + commandArgs.append(incomingData[i]) + + # report the incoming command + print('.. attempting to run command: %s\n with args: %s \n on mesh window %s ..' %(commandName, str(commandArgs), str(theMeshInstance))) + + #and try to run it + run_mp_command(commandName, commandArgs, theMeshInstance, self.mainWindow.ui, self) + +# try: +# run_mp_command(commandName, commandArgs, theMeshInstance, self.mainWindow.ui, self) +# +# except Exception as e: +# print("Bad command") +# self.socket.write(str('cmd error')) + + + + ## TODO: generate some feedback output + +#----------------- + + +class Dialog(QtWidgets.QDialog): + def __init__(self, parent=None): + super(Dialog, self).__init__(parent) + + self.server = mrMeshPyQtTCPServer() + if not self.server.listen(QtNetwork.QHostAddress("127.0.0.1"), 9995): + # if not self.server.listen(): #next available ## TODO may be better + print("Error") + return + + print("Running on %d" % self.server.serverPort()) + + + +if __name__ == '__main__': + import sys + app = QtWidgets.QApplication(sys.argv) + dialog = Dialog() + dialog.show() + sys.exit(dialog.exec_()) +