# Python script transalted from mrReadMrM.m in mrVista
# Andre Gouws 2009

from struct import unpack, calcsize
from numpy import *
import vtk

def LoadMrMFile(frame_parent, filename=None):

    #f = open('./left_full_bob.MrM','r')
    f = open('./sm_co20_test_right.MrM','r')
    #f = open('/groups/Projects/P1026/data/Anatomy/CM/Left/3DMeshes/CM_Left.MrM','r')
    
    
    signature = f.read(11) #filetype header TODO check if coorect type and exceptions
    print signature
    
    (flags, nStrips, nTriangles) = unpack('<lll', f.read(calcsize('<lll')))
    print flags, nStrips, nTriangles
    
    bounds = array(unpack('<6f', f.read(calcsize('<6f')))).reshape(3,2).transpose()
    print bounds
    
    (strip_off, strip_vertex_off, tri_vertex_off) = unpack('<lll', f.read(calcsize('<lll')))
    print strip_off, strip_vertex_off, tri_vertex_off
    
    #stripList = transpose(fread(MrMfile, [2,nStrips], 'long')); #orignal matlab
    nToGet = 2*nStrips
    stripList = array(unpack('<%dl' %nToGet, f.read(calcsize('<%dl' %nToGet))))
    stripList = stripList.reshape((nStrips,2))
    stripList = stripList.transpose()
    #stripList = ((array(unpack('<%dl' %nToGet, f.read(calcsize('<%dl' %nToGet))))).reshape((2,nStrips))).transpose()
    
    
    #nStrippedVerts = stripList(nStrips,1)+stripList(nStrips,2); #orignal matlab
    nStrippedVerts = stripList[0,nStrips-1]+stripList[1,nStrips-1]
    
    triangleOffset = nStrippedVerts + 1
    
    
    #sVertices = transpose(fread(MrMfile, [7,nStrippedVerts+3*nTriangles], 'float')); #orignal matlab
    nVertexValuesToGet = 7*(nStrippedVerts+3*nTriangles)
    sVertices = array(unpack('<%df' %nVertexValuesToGet, f.read(calcsize('<%df' %nVertexValuesToGet))))
    sVertices = sVertices.reshape(nStrippedVerts+3*nTriangles, 7)
    
    vertices = sVertices[:,0:3]
    normals = sVertices[:,3:6]
    
    
    if signature=='mr3DMesh_v2':
           
        GrayFile = f.read(60) #fgets(MrMfile, 60);
        GrayToAdd = unpack('<i' , f.read(calcsize('<i'))) #fread(MrMfile,1,'int');
        EncodeCurvature = unpack('<i' , f.read(calcsize('<i'))) #fread(MrMfile,1,'int');
        voxSize = unpack('<3f' , f.read(calcsize('<3f'))) #fread(MrMfile, [1,3], 'float');
        IsoSurfacer =  unpack('<i' , f.read(calcsize('<i'))) #fread(MrMfile,1,'int');
        
        # Overlays
        #
        OverlayROIs = unpack('<i' , f.read(calcsize('<i'))) #fread(MrMfile,1,'int');
        ClipToROIs = unpack('<i' , f.read(calcsize('<i'))) #fread(MrMfile,1,'int');
        OverlayFunctionals = unpack('<i' , f.read(calcsize('<i'))) #fread(MrMfile,1,'int');
        OverlayCuts = unpack('<i' , f.read(calcsize('<i'))) #fread(MrMfile,1,'int');
        ColorModDepth = unpack('<f' , f.read(calcsize('<f'))) #fread(MrMfile,1,'float');
        
        # Smoothing 1
        #
        smth1_Convergence = unpack('<f' , f.read(calcsize('<f'))) #fread(MrMfile,1,'float');
        smth1_NumberOfIterations = unpack('<i' , f.read(calcsize('<i'))) #fread(MrMfile,1,'int');
        smth1_RelaxationFactor = unpack('<f' , f.read(calcsize('<f'))) #fread(MrMfile,1,'float');
        smth1_FeatureAngle = unpack('<f' , f.read(calcsize('<f'))) #fread(MrMfile,1,'float');
        smth1_EdgeAngle = unpack('<f' , f.read(calcsize('<f'))) #fread(MrMfile,1,'float');
        smth1_FeatureEdgeSmoothing = unpack('<i' , f.read(calcsize('<i'))) #fread(MrMfile,1,'int');
        smth1_BoundarySmoothing = unpack('<i' , f.read(calcsize('<i'))) #fread(MrMfile,1,'int');
        smth1_Skip = unpack('<i' , f.read(calcsize('<i'))) #fread(MrMfile,1,'int');
        
        # Decimation
        #
        dec1_InitialFeatureAngle = unpack('<f' , f.read(calcsize('<f'))) #fread(MrMfile,1,'float');
        dec1_FeatureAngleIncrement = unpack('<f' , f.read(calcsize('<f'))) #fread(MrMfile,1,'float');
        dec1_MaximumFeatureAngle = unpack('<f' , f.read(calcsize('<f'))) #fread(MrMfile,1,'float');
        dec1_PreserveEdges = unpack('<i' , f.read(calcsize('<i'))) #fread(MrMfile,1,'int');
        dec1_BoundaryVertexDeletion = unpack('<i' , f.read(calcsize('<i'))) #fread(MrMfile,1,'int');
        dec1_InitialError = unpack('<f' , f.read(calcsize('<f'))) #fread(MrMfile,1,'float');
        dec1_ErrorIncrement = unpack('<f' , f.read(calcsize('<f'))) #fread(MrMfile,1,'float');
        dec1_MaximumError = unpack('<f' , f.read(calcsize('<f'))) #fread(MrMfile,1,'float');
        dec1_MaximumIterations = unpack('<i' , f.read(calcsize('<i'))) #fread(MrMfile,1,'int');
        dec1_MaximumSubIterations = unpack('<i' , f.read(calcsize('<i'))) #fread(MrMfile,1,'int');
        dec1_AspectRatio = unpack('<f' , f.read(calcsize('<f'))) #fread(MrMfile,1,'float');
        dec1_Degree = unpack('<i' , f.read(calcsize('<i'))) #fread(MrMfile,1,'int');
        dec1_PreserveTopology = unpack('<i' , f.read(calcsize('<i'))) #fread(MrMfile,1,'int');
        dec1_TgtNPolys = unpack('<i' , f.read(calcsize('<i'))) #fread(MrMfile,1,'int');
        dec1_Skip = unpack('<i' , f.read(calcsize('<i'))) #fread(MrMfile,1,'int');
        
        # Smoothing 2
        #
        smth2_Convergence = unpack('<f' , f.read(calcsize('<f'))) #fread(MrMfile,1,'float');
        smth2_NumberOfIterations = unpack('<i' , f.read(calcsize('<i'))) #fread(MrMfile,1,'int');
        smth2_RelaxationFactor = unpack('<f' , f.read(calcsize('<f'))) #fread(MrMfile,1,'float');
        smth2_FeatureAngle = unpack('<f' , f.read(calcsize('<f'))) #fread(MrMfile,1,'float');
        smth2_EdgeAngle = unpack('<f' , f.read(calcsize('<f'))) #fread(MrMfile,1,'float');
        smth2_FeatureEdgeSmoothing = unpack('<i' , f.read(calcsize('<i'))) #fread(MrMfile,1,'int');
        smth2_BoundarySmoothing = unpack('<i' , f.read(calcsize('<i'))) #fread(MrMfile,1,'int');
        smth2_Skip = unpack('<i' , f.read(calcsize('<i'))) #fread(MrMfile,1,'int');
        
        # Output
        #
        #garbage = fgets(MrMfile, 1);
        OutputFile = f.read(60)#fgets(MrMfile, 60);
        SaveOutput = unpack('<i' , f.read(calcsize('<i'))) #fread(MrMfile,1,'int');
        ViewOutput = unpack('<i' , f.read(calcsize('<i'))) #fread(MrMfile,1,'int');
        
        # Scale the coordinates into real physical units.
        #mesh.vertices = mrScaleCoords(mesh.vertices, mesh.parameters.structural.voxSize); #original matlab
        vertices[:,0] = vertices[:,0]/voxSize[0]
        vertices[:,1] = vertices[:,1]/voxSize[1]
        vertices[:,2] = vertices[:,2]/voxSize[2]
        #mesh.bounds = mrScaleCoords(mesh.bounds, mesh.parameters.structural.voxSize); #original matlab
        bounds[:,0] = bounds[:,0]/voxSize[0]
        bounds[:,1] = bounds[:,1]/voxSize[1]
        bounds[:,2] = bounds[:,2]/voxSize[2]    
    else:
        print '*** This mesh was made with mrGray 3.0. Use mrGray 4.1 or above for mrFlatMesh ****'
    #original matlab comment: Added 073001 ARW - now returns the original voxels size as saved out by mrGray.
    #voxelSize=voxSize
    
    
    
    f.close()
    
    #? after closing the file we re-open it and extract some more values? - the RGBA
    
    #f = open('./left_full_bob.MrM','r')
    f = open('./sm_co20_test_right.MrM','r')
    garbage = f.read(11)
    nToGet = 2*(nStrips+6)
    garbage = array(unpack('<%dl' %nToGet, f.read(calcsize('<%dl' %nToGet))))
    nToGet = 28*(nStrippedVerts+3*nTriangles)
    rgba = (array(unpack('<%dB' %nToGet, f.read(calcsize('<%dB' %nToGet))))).reshape(nStrippedVerts+3*nTriangles, 28)[:,24:]
    
    #return(vertices, rgba)
    #w = (vertices[triangleOffset-1:triangleOffset-1+nTriangles*3,:]).round() #-1 for zero indexing
    
    w = vertices
    
    print w[0:7,:]
        
    points = vtk.vtkPoints()
    for i in range(w.shape[0]):
        points.InsertPoint(i,w[i,0],w[i,1],w[i,2]) #ar - YNIC_FS
    
    triangles = vtk.vtkCellArray()
    
    
    initial_scalars = vtk.vtkFloatArray()
    
    #f = open('/Applications/MATLAB71/my_gray.out','r')
    #w = f.readlines()
    #my_gray = empty(len(w))
    #for i in range(len(w)):
    #    my_gray[i] = float(w[i].strip())
    #
    for i in range(w.shape[0]):
        initial_scalars.InsertNextValue(rgba[i,0]/255.0)

    #for i in range(len(w)):
    #    initial_scalars.InsertNextValue(my_gray[i]/255.0)
    #
    #for i in range(len(w)):
    #    if my_gray[i] < 128.0:
    #        initial_scalars.InsertNextValue(0.25)
        #else:
        #    initial_scalars.InsertNextValue(0.75)
    
    for i in range(nTriangles):
        triangles.InsertNextCell(3)
        triangles.InsertCellPoint(triangleOffset+3*i-4)
        triangles.InsertCellPoint(triangleOffset+3*i-3)
        triangles.InsertCellPoint(triangleOffset+3*i-2)
        if i == 1:
            print triangleOffset+3*i-4,triangleOffset+3*i-3,triangleOffset+3*i-2
        
        
    strips = vtk.vtkCellArray()
    for j in range(nStrips):
        curr_strip_num_points = stripList[1,j]
        start_offset = stripList[0,j]
        strips.InsertNextCell(curr_strip_num_points)
        for i in range(curr_strip_num_points):
            strips.InsertCellPoint(start_offset+i)


    # Assemble as PolyData
    polyData = vtk.vtkPolyData()
    polyData.SetPoints(points)
    polyData.SetPolys(triangles)
    polyData.SetStrips(strips)
    
    #return polyData
    
    tf = vtk.vtkTriangleFilter()
    tf.SetInput(polyData)
    tf.Update()

    #sm = vtk.vtkSmoothPolyDataFilter()
    #sm.SetNumberOfIterations(50)
    #sm.SetInput(tf.GetOutput())
    #sm.Update()

    pd_ScalarValues = polyData.GetPointData()
    pd_ScalarValues.SetScalars(initial_scalars)

    pdn = vtk.vtkPolyDataNormals()
    pdn.SetInput(tf.GetOutput())
    pdn.FlipNormalsOn()
    pdn.Update()

    t = vtk.vtkPolyDataMapper()
    t.SetInput(pdn.GetOutput())
    t.SetScalarRange(0.0,1.0)
    t.ScalarVisibilityOn()

    y = vtk.vtkLookupTable()
    y.SetNumberOfTableValues(256)
    y.SetHueRange(1.0,0.0)
    y.SetValueRange(0.25,0.75)
    y.SetSaturationRange(1.0,1.0)
    y.SetAlphaRange(1.0,1.0)
    y.SetTableRange(0.0,1.0)
    y.Build()

    #new_LUT.Build()

    t.SetLookupTable(y)
    t.SetColorModeToMapScalars()

    
    y = vtk.vtkActor()
    y.SetMapper(t)
    frame_parent.ren.AddActor(y)
    frame_parent.widget.Render()
     
    #recon = vtk.vtkSurfaceReconstructionFilter()
    #recon.SetInput(polyData)
    #
    #cf = vtk.vtkContourFilter()
    #cf.SetInput(recon.GetOutput())
    #cf.SetValue(0, 0.0)
    #
    #rev = vtk.vtkReverseSense()
    #rev.SetInput(cf.GetOutput())
    #rev.ReverseCellsOn()
    #rev.ReverseNormalsOn()
    #
    #t2 = vtk.vtkPolyDataMapper()
    #t2.SetInput(rev.GetOutput())
    #
    #y2 = vtk.vtkActor()
    #y2 .SetMapper(t2)
    #frame_parent.ren.AddActor(y2)
    #frame_parent.widget.Render()
    
    
    
    

