import vtk
import os, wx
from numpy import *
from struct import unpack, calcsize
from dv3dObjectsToRender import *
from dv3dPropertyFrames import *
from dv3dFreeSurferLoad import *
from common_functions import align_all_planes


def Generate_surface_from_loaded_volume(threshold, data_vol, frame_parent, data_max, data_source_path, reduction_factor, my_colour=(0.5,0.5,0.5), my_transp = 1, my_tree_group = None, overlay_or_not = 0, connectivity=0):
    
    if overlay_or_not == 0:
        a,b,c = extract_rough_surface(float(threshold), data_vol, data_max, my_colour, my_transp, reduction_factor, connectivity)
        AddVTKObjectWithAttributes(a,\
                    frame_parent,\
                    'vtkPolyDataSurfaceFromLoadedVolume',\
                    'surface at %.2f' %float(threshold),\
                    None,\
                    None,\
                    0,\
                    1,\
                    1,\
                    'structural',\
                    None,\
                    data_source_path,\
                    1,\
                    'SurfaceProperties',\
                    None)
    else:
        a,b,c = extract_rough_surface_for_overlays(float(threshold), data_vol, data_max, my_colour, my_transp, reduction_factor, connectivity)
        AddVTKObjectWithAttributes(a,\
                    frame_parent,\
                    'vtkPolyDataSurfaceFromLoadedVolume',\
                    'group surface at %.2f' %float(threshold),\
                    None,\
                    None,\
                    0,\
                    1,\
                    1,\
                    'structural',\
                    None,\
                    data_source_path,\
                    1,\
                    'SurfaceProperties',\
                    None)        
    
    

    ###TODO URGENT FIX
    ###frame_parent.ortho_window.ortho_ren1.AddActor(frame_parent.ListOfObjects[-1])
    ###frame_parent.ortho_window.ortho_ren2.AddActor(frame_parent.ListOfObjects[-1])
    ###frame_parent.ortho_window.ortho_ren3.AddActor(frame_parent.ListOfObjects[-1])
    frame_parent.ren.AddActor(frame_parent.ListOfObjects[-1])
    
    new_surface = frame_parent.tree.AppendItem(my_tree_group, frame_parent.ListOfObjects[-1].my_label)
    #add the newly created treeitem as an atribute of the plane_group
    frame_parent.ListOfObjects[-1].my_treeitem = new_surface
    frame_parent.tree.Refresh()
    
    # .. create a properties window for the plane set ...
    # we pass the parent frame and the newly created object as args
    frame = MyPropertyFrame(frame_parent, frame_parent.ListOfObjects[-1], 0)
    frame.Show(True)
    
    # ... and apply it to the relevant list object
    frame_parent.ListOfObjects[-1].my_property_frame_instance = frame
    
    align_all_planes(frame_parent)
    
    frame_parent.ListOfObjects[-1].GetMapper().ScalarVisibilityOff()
    frame_parent.ListOfObjects[-1].GetMapper().Update()
    frame_parent.ListOfObjects[-1].GetProperty().SetColor(my_colour)
    frame_parent.ListOfObjects[-1].GetProperty().SetOpacity(my_transp)
    
    frame_parent.ListOfObjects[-1].Modified()
    

def Load_FS_File(self, the_parent_frame):
    my_filetypes = '*' #the restiction for the file choice dialog
    surface_data_to_load = ChooseFile(self, my_filetypes)    
    
    print surface_data_to_load
    
    if surface_data_to_load[-4:] == '.asc':
            vertices, triangles = read_fs_ascii(surface_data_to_load)
    else:
            vertices, triangles = read_fs_surf(surface_data_to_load, 4)
            
    polyData = get_surf_polyData(vertices, triangles)
    


    #######TODOTODOTODOTODO



    deci = vtk.vtkDecimatePro() 
    deci.SetInput(polyData)
    deci.PreserveTopologyOn()
    deci.SetTargetReduction(0.5) #target decimation factor -> 0.9 = 90% reduction in number of vertices
    deci.Update()

    smooth = vtk.vtkSmoothPolyDataFilter()
    smooth.SetNumberOfIterations(10)
    smooth.SetInput(deci.GetOutput())
    
    #surfNormals = vtk.vtkPolyDataNormals() #required for correct display 
    #surfNormals.SetInput(smooth.GetOutput())
    #surfNormals.SetFeatureAngle(60.0)

    #conn = vtk.vtkPolyDataConnectivityFilter() #actually connect the dots now     
    #conn.SetInput(surfNormals.GetOutput())
    ##conn.AddSpecifiedRegion(0)
    #conn.SetExtractionModeToLargestRegion() #show only the largest connected surface
    #conn.Update()
    #
    # since the process above generates surfaces with (typically) 600k vertices,
    # this option allows us to downsample the surface structure. This is useful
    # for speeding up rendering (we spend a bitmore time here ~10 seconds)
    # TODO - option to display either the decimated / full surface

    
    #######TODOTODOTODOTODOTODOTODO
    
    
    
    origMapper = vtk.vtkPolyDataMapper()
    origMapper.ImmediateModeRenderingOn()
    #origMapper.SetInput(polyData)
    origMapper.SetInput(deci.GetOutput())
    origMapper.Update()
    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)

    #the_parent_frame.ren.AddActor(origActor)
    
    AddVTKObjectWithAttributes(origActor,\
                            the_parent_frame,\
                            'vtkPolyDataSurfaceFromOFF',\
                            'surface',\
                            None,\
                            None,\
                            0,\
                            1,\
                            1,\
                            'structural',\
                            None,\
                            surface_data_to_load,\
                            1,\
                            'SurfaceProperties',\
                            None)
    
    ###TODO URGENT FIX
    ###the_parent_frame.ortho_window.ortho_ren1.AddActor(the_parent_frame.ListOfObjects[-1])
    ###the_parent_frame.ortho_window.ortho_ren2.AddActor(the_parent_frame.ListOfObjects[-1])
    ###the_parent_frame.ortho_window.ortho_ren3.AddActor(the_parent_frame.ListOfObjects[-1])
    the_parent_frame.ren.AddActor(the_parent_frame.ListOfObjects[-1])
    
        
    new_surface = the_parent_frame.tree.AppendItem(the_parent_frame.tree_structural, the_parent_frame.ListOfObjects[-1].my_label)
    #add the newly created treeitem as an atribute of the plane_group
    the_parent_frame.ListOfObjects[-1].my_treeitem = new_surface
    the_parent_frame.tree.Refresh()
    
    # .. create a properties window for the plane set ...
    # we pass the parent frame and the newly created object as args
    frame = MyPropertyFrame(the_parent_frame, the_parent_frame.ListOfObjects[-1], 0)
    frame.Show(True)
    
    # ... and apply it to the relevant list object
    the_parent_frame.ListOfObjects[-1].my_property_frame_instance = frame
    
    align_all_planes(the_parent_frame)
    

    
    



def Load_OFF_File(self, the_parent_frame, my_file = None):
##def Load_OFF_File(the_parent_frame):
        
    if my_file is None:
        my_filetypes = '*' #the restiction for the file choice dialog
        surface_data_to_load = ChooseFile(self, my_filetypes)
    else:
        surface_data_to_load = my_file
    
    #surface_data_to_load = './DV3D_examples/cortex.off'   
    
    f=open(str(surface_data_to_load),'rb')
    d=f.readline() #ignore first line
       
    dim = unpack('>lll', f.read(calcsize('lll'))) #get dims
    print dim
    #p = fromfile(f, float32, dim[0]*3) 
    #vertices = reshape(p, (- 1,3))
    
    
    #TODO unpack('>%dfff' %dim[0]) etc to reove loops and speed up load
    vertices = []
    for i in range(dim[0]):
        vert = unpack('>fff', f.read(calcsize('fff')))
        vertices.append(vert)    
    
    triangles = []
    for i in range(dim[1]):
        tri = unpack('>lllll', f.read(calcsize('lllll')))
        triangles.append(tri)    
    
    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]),float(vertices[i][1]),float(vertices[i][2]))
    pts.Reset()

    # loop here for multiple parts
    
    for j in range(nTriangles):
        pts.Reset()
        
        for i in range(1,len(triangles[j])-1):
        
            pts.InsertNextId(triangles[j][i])
        
        newPolys.InsertNextCell(pts)
        
    output = vtk.vtkPolyData()
    output.SetPolys(newPolys)
    output.SetPoints(newPts)

    origMapper = vtk.vtkPolyDataMapper()
    origMapper.ImmediateModeRenderingOn()
    origMapper.SetInput(output)
    origMapper.Update()
    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)

    #the_parent_frame.ren.AddActor(origActor)
    
    AddVTKObjectWithAttributes(origActor,\
                            the_parent_frame,\
                            'vtkPolyDataSurfaceFromOFF',\
                            'surface',\
                            None,\
                            None,\
                            0,\
                            1,\
                            1,\
                            'structural',\
                            None,\
                            surface_data_to_load,\
                            1,\
                            'SurfaceProperties',\
                            None)
    
    ###TODO URGENT FIX
    ###the_parent_frame.ortho_window.ortho_ren1.AddActor(the_parent_frame.ListOfObjects[-1])
    ###the_parent_frame.ortho_window.ortho_ren2.AddActor(the_parent_frame.ListOfObjects[-1])
    ###the_parent_frame.ortho_window.ortho_ren3.AddActor(the_parent_frame.ListOfObjects[-1])
    the_parent_frame.ren.AddActor(the_parent_frame.ListOfObjects[-1])
    
        
    new_surface = the_parent_frame.tree.AppendItem(the_parent_frame.tree_structural, the_parent_frame.ListOfObjects[-1].my_label)
    #add the newly created treeitem as an atribute of the plane_group
    the_parent_frame.ListOfObjects[-1].my_treeitem = new_surface
    the_parent_frame.tree.Refresh()
    
    # .. create a properties window for the plane set ...
    # we pass the parent frame and the newly created object as args
    frame = MyPropertyFrame(the_parent_frame, the_parent_frame.ListOfObjects[-1], 0)
    frame.Show(True)
    
    # ... and apply it to the relevant list object
    the_parent_frame.ListOfObjects[-1].my_property_frame_instance = frame
    the_parent_frame.status_message = 'OFF Surface load completed.'
    
    align_all_planes(the_parent_frame)
    

def ChooseFile(self, my_filetypes):
    dlg = wx.FileDialog(self, "Choose file:", '', "", my_filetypes, wx.OPEN)
    if dlg.ShowModal() == wx.ID_OK:
        path = dlg.GetPath()
        return path


#------------------------------------------------------------------------------
#generate a rough surface fropm scalar data in the underlying
# volume dataset - can be used to extract a rough cortex / 3d-stat-blob
def extract_rough_surface(value, vol, data_max, my_colour, my_transp, reduction_factor, connectivity):
    
    #print reduction_factor
    brain_volume = vol.GetOutput() #get the underlying volume
    brain_volume.Update() #actually load the data

    shrink = vtk.vtkImageShrink3D()
    shrink.SetShrinkFactors(reduction_factor, reduction_factor, reduction_factor)
    shrink.SetInput(brain_volume)
    shrink.AveragingOn()
      
    isoSurfaceExtractor = vtk.vtkContourFilter() #used to generate a suface by connecting scalars of equal value
    isoSurfaceExtractor.UseScalarTreeOn()
    isoSurfaceExtractor.SetInput(shrink.GetOutput())
    isoSurfaceExtractor.SetValue(0, value) #sets the 'range' of the dots to connect
    isoSurfaceExtractor.ComputeGradientsOn()
    isoSurfaceExtractor.UseScalarTreeOn()
    
    smooth = vtk.vtkSmoothPolyDataFilter()
    smooth.SetNumberOfIterations(500)
    smooth.SetInput(isoSurfaceExtractor.GetOutput())
    
    surfNormals = vtk.vtkPolyDataNormals() #required for correct display 
    surfNormals.SetInput(smooth.GetOutput())
    surfNormals.SetFeatureAngle(60.0)

    conn = vtk.vtkPolyDataConnectivityFilter() #actually connect the dots now     
    conn.SetInput(surfNormals.GetOutput())
    #conn.AddSpecifiedRegion(0)
    if connectivity == 1:
        conn.SetExtractionModeToLargestRegion() #show only the largest connected surface
    else:
        conn.SetExtractionModeToAllRegions()
    conn.Update()
    
    # since the process above generates surfaces with (typically) 600k vertices,
    # this option allows us to downsample the surface structure. This is useful
    # for speeding up rendering (we spend a bitmore time here ~10 seconds)
    # TODO - option to display either the decimated / full surface
    deci = vtk.vtkDecimatePro() 
    deci.SetInput(conn.GetOutput())
    deci.PreserveTopologyOn()
    deci.SetTargetReduction(0.99) #target decimation factor -> 0.9 = 90% reduction in number of vertices
    deci.Update()

    #now pass our surface to the tool that draws(or 'maps') the data to the
    #screen
    surfMapper = vtk.vtkPolyDataMapper()
    surfMapper.ImmediateModeRenderingOn()
    surfMapper.SetInput(deci.GetOutput())
    surfMapper.ScalarVisibilityOff() 
    surfMapper.Update()
    surf = vtk.vtkActor()
    surf.GetProperty().SetLineWidth(5.005) 
    surf.SetMapper(surfMapper)
    
    surf.GetProperty().SetColor(my_colour)
    surf.GetProperty().SetOpacity(my_transp)

    return surf,deci,conn

#------------------------------------------------------------------------------
#generate a rough surface fropm scalar data in the underlying
# volume dataset - can be used to extract a rough cortex / 3d-stat-blob
def extract_rough_surface_for_overlays(value, vol, data_max, my_colour, my_transp, reduction_factor, connectivity):
    
    brain_volume = vol.GetOutput() #get the underlying volume
    brain_volume.Update() #actually load the data

    isoSurfaceExtractor = vtk.vtkContourFilter() #used to generate a suface by connecting scalars of equal value
    isoSurfaceExtractor.SetInput(brain_volume)
    isoSurfaceExtractor.SetValue(0, value) #sets the 'range' of the dots to connect
    isoSurfaceExtractor.ComputeGradientsOn()
    isoSurfaceExtractor.UseScalarTreeOn()
    
    surfNormals = vtk.vtkPolyDataNormals() #required for correct display 
    surfNormals.SetInput(isoSurfaceExtractor.GetOutput())
    surfNormals.SetFeatureAngle(60.0)

    conn = vtk.vtkPolyDataConnectivityFilter() #actually connect the dots now     
    conn.SetInput(surfNormals.GetOutput())
    #conn.AddSpecifiedRegion(0)
    if connectivity == 1:
        conn.SetExtractionModeToLargestRegion() #show only the largest connected surface
    else:
        conn.SetExtractionModeToAllRegions()
        conn.Update()
    
    # since the process above generates surfaces with (typically) 600k vertices,
    # this option allows us to downsample the surface structure. This is useful
    # for speeding up rendering (we spend a bitmore time here ~10 seconds)
    # TODO - option to display either the decimated / full surface
    deci = vtk.vtkDecimatePro() 
    deci.SetInput(conn.GetOutput())
    deci.PreserveTopologyOff()
    deci.SetTargetReduction(0.001) #target decimation factor -> 0.9 = 90% reduction in number of vertices
    deci.Update()

    #now pass our surface to the tool that draws(or 'maps') the data to the
    #screen
    surfMapper = vtk.vtkPolyDataMapper()
    surfMapper.ImmediateModeRenderingOn()
    surfMapper.SetInput(deci.GetOutput())
    surfMapper.ScalarVisibilityOff() 
    surfMapper.Update()
    
    #propertyBack = vtk.vtkProperty()
    #propertyBack.SetColor(my_colour)
    #propertyBack.SetOpacity(my_transp)
    
    surf = vtk.vtkActor()
    surf.GetProperty().SetLineWidth(0.005) 
    surf.SetMapper(surfMapper)
    surf.GetProperty().SetColor(my_colour)
    surf.GetProperty().SetOpacity(my_transp)

    #surf.SetBackfaceProperty(propertyBack)
    surfMapper.ScalarVisibilityOff() 
    surfMapper.Update()

    return surf,deci,conn



