import vtk
from vtk.util.colors import tomato, banana
from struct import unpack, calcsize #from surface
from numpy import *
import scipy as Sci
import scipy.linalg
import sys
import os
from subprocess import Popen, PIPE


#set globals 
x_const = 0
y_const = 0
z_const = 0


def get_path(path):
    return os.path.join(os.path.split(__file__)[0], path)

# --------- File handling ------------

def my_file_dialog(processing_type):
    p = Popen('pythonw %s %s' % (get_path('ynicDV3D_wxdialogs.py'), processing_type), stdout=PIPE, shell=True)
    p.wait() 
    (so, se) = p.communicate()
    return so.strip()

### IMPORT FIBER DATA #################################################################
#____________________________________________________________________________________

plotData = []

def import_fiber_data(fiber_data_to_load):
    global x_const, y_const, z_constz, plotData
    #get the offsets of the original volume
    dti_ref_filename = my_file_dialog(7)
    if dti_ref_filename == '':
        print 'Load cancelled' 
        return
    
    mat_filename = my_file_dialog(8)
    if mat_filename == '':
        print 'Load cancelled' 
        return
        
    p = os.popen('avwhd %s' %dti_ref_filename)
    d = p.readlines()
    p.close()
    
    orig_offsetList=[]
    
    for c in d:
        s = c.split()
        try:
            if s[0][0:8] == 'qto_xyz:':
                orig_offsetList.append(float(s[4]))
        except:
            pass
    
    #print orig_offsetList
    
    coreg_corections = []
    
    f = open(mat_filename)
    d = f.readlines()
    f.close()
    
    for c in d:
        s = c.split()
        try:
            coreg_corections.append(float(s[3]))
        except:
            pass
            
    print coreg_corections
    
    manual_correction=[0.0,0.0,0.0]
    
    #co-registration is never perfect .. and for some reason the human brain seems to be a
    #slightly more advanced processing toolbox than the average pc / algorithm .. as such 
    #you may feel you can improve on the gross co-registration by tweaking the placement of
    #fiber bundles with these parameters
    manual_correction[0] = float(my_file_dialog(16))
    if manual_correction[0] == '':
        print 'Load cancelled' 
        return
    manual_correction[1] = float(my_file_dialog(15))
    if manual_correction[1] == '':
        print 'Load cancelled' 
        return
    manual_correction[2] = float(my_file_dialog(14))
    if manual_correction[1] == '':
        print 'Load cancelled' 
        return
        
    x_const = coreg_corections[0] + orig_offsetList[0] + manual_correction[0]
    y_const = coreg_corections[1] - orig_offsetList[1] + manual_correction[1]
    z_const = coreg_corections[2] - orig_offsetList[2] + manual_correction[2]
    
    lineData = []
    plotData = []
    f = open(fiber_data_to_load)
    s = f.readlines()
    f.close()
    for x in s:
        lineData.append([float(y) for y in x.split()])

    for t in range(len(lineData)):
        currentLineData = []
        points_per_line = int(lineData[t][0])
        for m in range(points_per_line):
            currentLineData.append([lineData[t][m+1],lineData[t][m+points_per_line+1],lineData[t][m+points_per_line*2+1]])
        
        plotData.append(currentLineData)

    inputPoints = vtk.vtkPoints()
    
    # Set up some empty arrays to store our calculated / extracted values in   - One spline for each direction.
    aSplineX = vtk.vtkCardinalSpline()
    aSplineY = vtk.vtkCardinalSpline()
    aSplineZ = vtk.vtkCardinalSpline()
    
    # List of points for each line to be drawn through
    glyphList = []
    
    # List of line profiles to be drawn through the points
    profileList = []
             
               
    # --------- START THE PROCESSING ------------
    #----------Add the fiber data ---------------------------
    
    fibres = []
    
    for j in range(len(plotData)):
        
        # Total number of points.
        numberOfInputPoints = int(len(plotData[j]))
        
        #print numberOfInputPoints
        
        for i in range(numberOfInputPoints):
        
            # Extract values for each point
            x = plotData[j][i][0]+x_const   #add half FOV of DTI set to .mat transform of t1Flair to t1
            y = plotData[j][i][1]+y_const
            z = plotData[j][i][2]+z_const
            
            # Add the corresponding coordinates to the splines.
            # aSplineX will interpolate the x values of the points
            # aSplineY will interpolate the y values of the points
            # aSplineZ will interpolate the z values of the points
            aSplineX.AddPoint(i, x)
            aSplineY.AddPoint(i, y)
            aSplineZ.AddPoint(i, z)
            inputPoints.InsertPoint(i, x, y, z)
        
        # The following section will create glyphs for the pivot points
        # in order to make the effect of the spline more clear.
        
        # Create a polydata to be glyphed.
        inputData = vtk.vtkPolyData()
        inputData.SetPoints(inputPoints)
        
        # Use sphere as glyph source.
        balls = vtk.vtkSphereSource()
        balls.SetRadius(10.01)
        balls.SetPhiResolution(10)
        balls.SetThetaResolution(10)
        
        glyphPoints = vtk.vtkGlyph3D()
        glyphPoints.SetInput(inputData)
        glyphPoints.SetSource(balls.GetOutput())
        
        glyphMapper = vtk.vtkPolyDataMapper()
        glyphMapper.SetInput(glyphPoints.GetOutput())
        
        glyphList.append(vtk.vtkActor())
        glyphList[j].SetMapper(glyphMapper)
        glyphList[j].GetProperty().SetDiffuseColor(tomato)
        glyphList[j].GetProperty().SetSpecular(.3)
        glyphList[j].GetProperty().SetSpecularPower(30)
        
        # Generate the polyline for the spline.
        points = vtk.vtkPoints()
        profileData = vtk.vtkPolyData()
        
        # Number of points on the spline
        numberOfOutputPoints = 25
        
        # Interpolate x, y and z by using the three spline filters and
        # create new points
        for i in range(0, numberOfOutputPoints):
            t = (numberOfInputPoints-1.0)/(numberOfOutputPoints-1.0)*i
            points.InsertPoint(i, aSplineX.Evaluate(t), aSplineY.Evaluate(t),
                               aSplineZ.Evaluate(t))
        
        # Create the polyline.
        lines = vtk.vtkCellArray()
        lines.InsertNextCell(numberOfOutputPoints)
        for i in range(0, numberOfOutputPoints):
            lines.InsertCellPoint(i)
         
        profileData.SetPoints(points)
        profileData.SetLines(lines)
        
        # Add thickness to the resulting line.
        profileTubes = vtk.vtkTubeFilter()
        profileTubes.SetNumberOfSides(1)
        profileTubes.SetInput(profileData)
        profileTubes.SetRadius(0.1)
        #profileTubes = vtk.vtkSplineFilter()
        #profileTubes.SetInput(profileData)
        
        profileMapper = vtk.vtkPolyDataMapper()
        profileMapper.SetInput(profileTubes.GetOutput())
        
    #    glyphList.append(vtk.vtkActor())
        profileList.append(vtk.vtkActor())
        profileList[j].SetMapper(profileMapper)
        profileList[j].GetProperty().SetDiffuseColor(banana)
        profileList[j].GetProperty().SetSpecular(.3)
        profileList[j].GetProperty().SetSpecularPower(30)
        profileList[j].VisibilityOff()
        profileList[j].SetPickable(0)
                        
    return (profileList, lineData, plotData, x_const, y_const, z_const)





### IMPORT SURFACE DATA #################################################################
#____________________________________________________________________________________

def import_cortex_data(surface_data_to_load):
    global vertices, triangles

    f=open(surface_data_to_load,'rb')
    d=f.readline()                                                                                                                                    
    dim = unpack('>lll', f.read(calcsize('lll')))
    p = fromfile(f, float32, dim[0]*3) 
    vertices = reshape(p, (- 1,3))
        
    t = fromfile(f, int, dim[1]*5)
    tmp = reshape(t, (-1,5))
    triangles = tmp[:,0:4].copy()
    f.close()
    nVertices = dim[0]
    nTriangles = dim[1]
    
    pts = vtk.vtkIdList()
    newPolys = vtk.vtkCellArray()
    newPts = vtk.vtkPoints()

    # generate parts
    for i in range(nVertices):
        #newPts.InsertPoint(i, float(vertices[i][0])+191,float(vertices[i][1])-15,float(vertices[i][2])-85) #.. hack for piers
        newPts.InsertPoint(i, float(vertices[i][0]),float(vertices[i][1]),float(vertices[i][2])) #.. hack for piers

    pts.Reset()

    # loop here for multiple parts
    
    for j in range(nTriangles):
        pts.Reset()
        
        for i in range(1,len(triangles[j])):
        
            pts.InsertNextId(triangles[j][i])
        
        newPolys.InsertNextCell(pts)
        
    output = vtk.vtkPolyData()
    output.SetPolys(newPolys)
    output.SetPoints(newPts)
    
    origMapper = vtk.vtkPolyDataMapper()
    origMapper.SetInput(output)
    origActor = vtk.vtkLODActor()
    origActor.SetMapper(origMapper)
    origActor.VisibilityOn()
    
    
    # --- Add the Surface Actor with a set opacity and Colour
    origActor.GetProperty().SetOpacity(1.0)
    origActor.GetProperty().SetColor(0.8, 0.0, 0.0)
        
    return(origActor,output)

#----------END Add a cortical surface ---------------------------
######## END CORTEX -- #############################################################







