def filterUGridIntoPData(
        ugrid,
        only_trianlges=False,
        verbose=0):

    myVTK.myPrint(verbose, "*** filterUGridIntoPData ***")

    filter_geometry = vtk.vtkGeometryFilter()
    if (vtk.vtkVersion.GetVTKMajorVersion() >= 6):
        filter_geometry.SetInputData(ugrid)
    else:
        filter_geometry.SetInput(ugrid)
    filter_geometry.Update()
    pdata = filter_geometry.GetOutput()

    if (only_trianlges):
        filter_triangle = vtk.vtkTriangleFilter()
        if (vtk.vtkVersion.GetVTKMajorVersion() >= 6):
            filter_triangle.SetInputData(pdata)
        else:
            filter_triangle.SetInput(pdata)
        filter_triangle.Update()
        pdata = filter_triangle.GetOutput()

    return pdata

def computeABPointsFromBoundsAndCenter(
        mesh,
        AB=[0,0,1],
        verbose=0):

    myVTK.myPrint(verbose, "*** computeABPointsFromBoundsAndCenter ***")

    C = numpy.array(mesh.GetCenter())
    #print "C ="+str(C)

    bounds = mesh.GetBounds()
    diag = numpy.array([bounds[1]-bounds[0], bounds[3]-bounds[2], bounds[5]-bounds[4]])
    AB = numpy.array(AB)
    AB = abs(numpy.dot(diag, AB)) * AB
    #print "bounds ="+str(bounds)
    #print "diag ="+str(diag)
    #print "AB ="+str(AB)

    point_A = C - AB/2
    point_B = C + AB/2
    #print "point_A ="+str(point_A)
    #print "point_B ="+str(point_B)

    points_AB = vtk.vtkPoints()
    points_AB.InsertNextPoint(point_A)
    points_AB.InsertNextPoint(point_B)
    #print points_AB

    return points_AB

def clipSurfacesForFullLVMesh(endo, epi, verbose=1):

    myVTK.myPrint(verbose, "*** clipSurfacesForFullLVMesh ***")

    endo_implicit_distance = vtk.vtkImplicitPolyDataDistance()
    endo_implicit_distance.SetInput(endo)

    epi_implicit_distance = vtk.vtkImplicitPolyDataDistance()
    epi_implicit_distance.SetInput(epi)

    epi_clip = vtk.vtkClipPolyData()
    epi_clip.SetInputData(epi)
    epi_clip.SetClipFunction(endo_implicit_distance)
    epi_clip.GenerateClippedOutputOn()
    epi_clip.Update()
    clipped_epi = epi_clip.GetOutput(0)
    clipped_valve = epi_clip.GetOutput(1)

    endo_clip = vtk.vtkClipPolyData()
    endo_clip.SetInputData(endo)
    endo_clip.SetClipFunction(epi_implicit_distance)
    endo_clip.InsideOutOn()
    endo_clip.Update()
    clipped_endo = endo_clip.GetOutput(0)

    return (clipped_endo, clipped_epi, clipped_valve)

def computeRegionsForBiV(points, pdata_endLV, pdata_endRV, pdata_epi, verbose=0):

    myVTK.myPrint(verbose, "*** computeRegionsForBiV ***")

    myVTK.myPrint(verbose, "Initializing cell locators...")

    (cell_locator_endLV, closest_point_endLV, generic_cell, cellId_endLV, subId, dist_endLV) = myVTK.getCellLocator(
        mesh=pdata_endLV, verbose=verbose - 1
    )
    (cell_locator_endRV, closest_point_endRV, generic_cell, cellId_endRV, subId, dist_endRV) = myVTK.getCellLocator(
        mesh=pdata_endRV, verbose=verbose - 1
    )
    (cell_locator_epi, closest_point_epi, generic_cell, cellId_epi, subId, dist_epi) = myVTK.getCellLocator(
        mesh=pdata_epi, verbose=verbose - 1
    )

    n_points = points.GetNumberOfPoints()

    iarray_region = myVTK.createIntArray("region_id", 1, n_points)

    for k_point in range(n_points):
        point = numpy.array(points.GetPoint(k_point))
        cell_locator_endLV.FindClosestPoint(point, closest_point_endLV, generic_cell, cellId_endLV, subId, dist_endLV)
        cell_locator_endRV.FindClosestPoint(point, closest_point_endRV, generic_cell, cellId_endRV, subId, dist_endRV)
        cell_locator_epi.FindClosestPoint(point, closest_point_epi, generic_cell, cellId_epi, subId, dist_epi)

        if dist_endRV == max(dist_endLV, dist_endRV, dist_epi):
            iarray_region.SetTuple(k_point, [0])
        elif dist_epi == max(dist_endLV, dist_endRV, dist_epi):
            iarray_region.SetTuple(k_point, [1])
        elif dist_endLV == max(dist_endLV, dist_endRV, dist_epi):
            iarray_region.SetTuple(k_point, [2])

    return iarray_region

def computeFractionalAnisotropy(
        farray_e1,
        farray_e2,
        farray_e3,
        verbose=1):

    myVTK.myPrint(verbose, "*** computeFractionalAnisotropy ***")

    n_tuples = farray_e1.GetNumberOfTuples()

    farray_FA   = myVTK.createFloatArray("FA"   , 1, n_tuples)
    farray_FA12 = myVTK.createFloatArray("FA_12", 1, n_tuples)
    farray_FA23 = myVTK.createFloatArray("FA_23", 1, n_tuples)

    for k_tuple in xrange(n_tuples):
        e1 = farray_e1.GetTuple1(k_tuple)
        e2 = farray_e2.GetTuple1(k_tuple)
        e3 = farray_e3.GetTuple1(k_tuple)
        FA   = ((e1-e2)**2+(e1-e3)**2+(e2-e3)**2)**(0.5) / (2*(e1**2+e2**2+e3**2))**(0.5)
        FA12 = ((e1-e2)**2)**(0.5) / (e1**2+e2**2)**(0.5)
        FA23 = ((e2-e3)**2)**(0.5) / (e2**2+e3**2)**(0.5)

        farray_FA.SetTuple1(k_tuple, FA)
        farray_FA12.SetTuple1(k_tuple, FA12)
        farray_FA23.SetTuple1(k_tuple, FA23)

    return (farray_FA,
            farray_FA12,
            farray_FA23)

def thresholdUGrid(
        ugrid_mesh,
        field_support,
        field_name,
        threshold_value,
        threshold_by_upper_or_lower,
        verbose=0):

    myVTK.myPrint(verbose, "*** thresholdUGrid ***")

    threshold = vtk.vtkThreshold()
    if (vtk.vtkVersion.GetVTKMajorVersion() >= 6):
        threshold.SetInputData(ugrid_mesh)
    else:
        threshold.SetInput(ugrid_mesh)
    if (field_support == "points"):
        association = vtk.vtkDataObject.FIELD_ASSOCIATION_POINTS
    elif (field_support == "cells"):
        association = vtk.vtkDataObject.FIELD_ASSOCIATION_CELLS
    threshold.SetInputArrayToProcess(0, 0, 0, association, field_name)
    if (threshold_by_upper_or_lower == "upper"):
        threshold.ThresholdByUpper(threshold_value)
    elif (threshold_by_upper_or_lower == "lower"):
        threshold.ThresholdByLower(threshold_value)
    threshold.Update()
    ugrid_thresholded_mesh = threshold.GetOutput()
    return ugrid_thresholded_mesh

def readPData(
        filename,
        verbose=1):

    myVTK.myPrint(verbose, "*** readPData: " + filename + " ***")

    if ('vtk' in filename):
        pdata_reader = vtk.vtkPolyDataReader()
    elif ('vtp' in filename):
        pdata_reader = vtk.vtkXMLPolyDataReader()
    else:
        assert 0, "File must be .vtk or .vtp. Aborting."

    assert (os.path.isfile(filename)), "Wrong filename. Aborting."

    pdata_reader.SetFileName(filename)
    pdata_reader.Update()
    pdata = pdata_reader.GetOutput()

    if (verbose):
        print "n_points = " + str(pdata.GetNumberOfPoints())
        print "n_verts = " + str(pdata.GetNumberOfVerts())
        print "n_lines = " + str(pdata.GetNumberOfLines())
        print "n_polys = " + str(pdata.GetNumberOfPolys())
        print "n_strips = " + str(pdata.GetNumberOfStrips())

    return pdata

def addSyntheticHelixAngles(
        ugrid,
        helix_angle_end,
        helix_angle_epi,
        type_of_support="cell",
        verbose=1):

    myVTK.myPrint(verbose, "*** addSyntheticHelixAngles ***")

    if (type_of_support == "cell"):
        ugrid_data = ugrid.GetCellData()
    elif (type_of_support == "point"):
        ugrid_data = ugrid.GetPointData()

    farray_rr = ugrid_data.GetArray("rr")
    farray_angle_helix = ugrid_data.GetArray("angle_helix")

    farray_angle_helix = computeSyntheticHelixAngles(
        farray_rr=farray_rr,
        helix_angle_end=helix_angle_end,
        helix_angle_epi=helix_angle_epi,
        farray_angle_helix=farray_angle_helix,
        verbose=verbose-1)

    ugrid_data.AddArray(farray_angle_helix)

    return farray_angle_helix

def computeSyntheticHelixAngles(
        farray_rr,
        helix_angle_end,
        helix_angle_epi,
        farray_angle_helix=None,
        verbose=1):

    myVTK.myPrint(verbose, "*** computeSyntheticHelixAngles ***")

    n_cells = farray_rr.GetNumberOfTuples()

    if (farray_angle_helix is None):
        farray_angle_helix = myVTK.createFloatArray(
            name="angle_helix",
            n_components=1,
            n_tuples=n_cells)

    for k_cell in xrange(n_cells):
        rr = farray_rr.GetTuple1(k_cell)

        helix_angle_in_degrees = (1.-rr) * helix_angle_end \
                               +     rr  * helix_angle_epi
        farray_angle_helix.SetTuple1(
            k_cell,
            helix_angle_in_degrees)

    return farray_angle_helix

def computeSectorsForLV(
        farray_rr,
        farray_cc,
        farray_ll,
        n_r=1,
        n_c=1,
        n_l=1,
        iarray_part_id=None,
        verbose=0):

    myVTK.myPrint(verbose, "*** computeSectorsForLV ***")

    n_cells = farray_rr.GetNumberOfTuples()

    iarray_sector = myVTK.createIntArray("sector_id", 1, n_cells)

    for k_cell in range(n_cells):
        if (iarray_part_id is not None) and (int(iarray_part_id.GetTuple1(k_cell)) > 0):
            sector_id = -1

        else:
            rr = farray_rr.GetTuple1(k_cell)
            cc = farray_cc.GetTuple1(k_cell)
            ll = farray_ll.GetTuple1(k_cell)

            k_r = int(rr*n_r/1.000001)
            k_c = int(cc*n_c/1.000001)
            k_l = int((1.-ll)*n_l/1.000001)

            sector_id = k_l * n_c * n_r + k_c * n_r + k_r

        iarray_sector.SetTuple1(k_cell, sector_id)

    return iarray_sector

def computeMeanStddevAngles(
        angles,
        angles_in_degrees=True,
        angles_in_pm_pi=True,
        verbose=1):

    myVTK.myPrint(verbose, "*** computeMeanStddevAngles ***")

    if (angles_in_degrees):
        if (angles_in_pm_pi):
            mean = math.atan2(numpy.mean([numpy.sin(2*numpy.array(angles)*numpy.pi/180)]),
                              numpy.mean([numpy.cos(2*numpy.array(angles)*numpy.pi/180)]))*180/math.pi/2
        else:
            mean = math.atan2(numpy.mean([numpy.sin(numpy.array(angles)*numpy.pi/180)]),
                              numpy.mean([numpy.cos(numpy.array(angles)*numpy.pi/180)]))*180/math.pi

        stddev = numpy.sqrt(numpy.mean(((((numpy.array(angles)-mean)+90)%180)-90)**2))
    else:
        if (angles_in_pm_pi):
            mean = math.atan2(numpy.mean([numpy.sin(2*numpy.array(angles))]),
                              numpy.mean([numpy.cos(2*numpy.array(angles))]))/2
        else:
            mean = math.atan2(numpy.mean([numpy.sin(numpy.array(angles))]),
                              numpy.mean([numpy.cos(numpy.array(angles))]))

        stddev = numpy.sqrt(numpy.mean(((((numpy.array(angles)-mean)+math.pi/2)%math.pi)-math.pi/2)**2))

    return (mean, stddev)

def computeImageGradient(
        image=None,
        image_filename=None,
        image_dimensionality=None,
        verbose=0):

    myVTK.myPrint(verbose, "*** computeImageGradient ***")

    image = myVTK.initImage(image, image_filename, verbose-1)

    if (image_dimensionality is None):
        image_dimensionality = myVTK.computeImageDimensionality(
            image=image,
            verbose=verbose-1)

    image_gradient = vtk.vtkImageGradient()
    if (vtk.vtkVersion.GetVTKMajorVersion() >= 6):
        image_gradient.SetInputData(image)
    else:
        image_gradient.SetInput(image)
    image_gradient.SetDimensionality(image_dimensionality)
    image_gradient.Update()
    image_w_grad = image_gradient.GetOutput()

    return image_w_grad

def mulArrays(
    array1,
    array2,
    array3=None,
    verbose=0):

    myVTK.myPrint(verbose, "*** mulArrays ***")

    n_components = array1.GetNumberOfComponents()
    assert (array2.GetNumberOfComponents() == n_components)

    n_tuples = array1.GetNumberOfTuples()
    assert (array2.GetNumberOfTuples() == n_tuples)

    array_type = type(array1.GetTuple(0)[0])
    assert (array_type in [int, float])
    assert (type(array2.GetTuple(0)[0]) is array_type)

    if (array3 is None):
        array3 = myVTK.createArray(
            name="",
            n_components=n_components,
            n_tuples=n_tuples,
            array_type=array_type)
    else:
        assert (array3.GetNumberOfComponents() == n_components)
        assert (array3.GetNumberOfTuples() == n_tuples)
        assert (type(array3.GetTuple(0)[0]) is array_type)

    for k_tuple in xrange(n_tuples):
        array3.SetTuple(
            k_tuple,
            numpy.array(array1.GetTuple(k_tuple)) * numpy.array(array2.GetTuple(k_tuple)))

    return array3

def readAbaqusDeformationGradientsFromDAT(
        data_filename,
        verbose=0):

    myVTK.myPrint(verbose, "*** readAbaqusDeformationGradientsFromDAT: "+data_filename+" ***")

    farray_F = myVTK.createFloatArray("F", 9)

    data_file = open(data_filename, 'r')
    context = ""
    k_cell = 0
    for line in data_file:
        if (context == "reading deformation gradients"):
            #print line
            if ("MAXIMUM" in line):
                context = ""
                continue
            if ("OR" in line):
                splitted_line = line.split()
                assert (int(splitted_line[0]) == k_cell+1), "Wrong element number. Aborting."
                F_list = [float(splitted_line[ 3]), float(splitted_line[ 6]), float(splitted_line[7]),
                          float(splitted_line[ 9]), float(splitted_line[ 4]), float(splitted_line[8]),
                          float(splitted_line[10]), float(splitted_line[11]), float(splitted_line[5])]
                farray_F.InsertNextTuple(F_list)
                k_cell += 1

        if (line == "    ELEMENT  PT FOOT-       DG11        DG22        DG33        DG12        DG13        DG23        DG21        DG31        DG32    \n"):
            context = "reading deformation gradients"

    data_file.close()

    myVTK.myPrint(verbose-1, "n_tuples = "+str(farray_F.GetNumberOfTuples()))

    return farray_F

def addCartesianCoordinates(
        ugrid,
        verbose=1):

    myVTK.myPrint(verbose, "*** addCartesianCoordinates ***")

    points = ugrid.GetPoints()
    (farray_xx,
     farray_yy,
     farray_zz) = computeCartesianCoordinates(
        points=points,
        verbose=verbose-1)

    ugrid.GetPointData().AddArray(farray_xx)
    ugrid.GetPointData().AddArray(farray_yy)
    ugrid.GetPointData().AddArray(farray_zz)

    cell_centers = myVTK.getCellCenters(
        mesh=ugrid,
        verbose=verbose-1)
    (farray_xx,
     farray_yy,
     farray_zz) = computeCartesianCoordinates(
        points=cell_centers,
        verbose=verbose-1)

    ugrid.GetCellData().AddArray(farray_xx)
    ugrid.GetCellData().AddArray(farray_yy)
    ugrid.GetCellData().AddArray(farray_zz)

def writeFiberOrientationFileForAbaqus(
        mesh,
        filename,
        eF_field_name="eF",
        eS_field_name="eS",
        sep=", ",
        verbose=0):

    myVTK.myPrint(verbose, "*** writeFiberOrientationFileForAbaqus ***")

    orientation_file = open(filename, "w")
    orientation_file.write(", 1., 0., 0., 0., 1., 0."+"\n")

    n_cells = mesh.GetNumberOfCells()

    eF_array = mesh.GetCellData().GetArray(eF_field_name)
    eS_array = mesh.GetCellData().GetArray(eS_field_name)
    eF = numpy.empty(3)
    eS = numpy.empty(3)

    for k_cell in xrange(n_cells):
        eF_array.GetTuple(k_cell, eF)
        eS_array.GetTuple(k_cell, eS)

        line = str(k_cell+1)
        for k in xrange(3): line += sep + str(eF[k])
        for k in xrange(3): line += sep + str(eS[k])
        line += "\n"
        orientation_file.write(line)

    orientation_file.close()

def clipSurfacesForCutLVMesh(
        endo,
        epi,
        height,
        verbose=1):

    myVTK.myPrint(verbose, "*** clipSurfacesForCutLVMesh ***")

    plane = vtk.vtkPlane()
    plane.SetNormal(0,0,-1)
    plane.SetOrigin(0,0,height)

    clip = vtk.vtkClipPolyData()
    clip.SetClipFunction(plane)
    clip.SetInputData(endo)
    clip.Update()
    clipped_endo = clip.GetOutput(0)

    clip = vtk.vtkClipPolyData()
    clip.SetClipFunction(plane)
    clip.SetInputData(epi)
    clip.Update()
    clipped_epi = clip.GetOutput(0)

    return (clipped_endo,
            clipped_epi)

def readDynaDeformationGradients(
        mesh,
        hystory_files_basename,
        array_name,
        verbose=1):

    myVTK.myPrint(verbose, "*** readDynaDeformationGradients ***")

    n_cells = mesh.GetNumberOfCells()

    history_files_names = [hystory_files_basename + '.history#' + str(num) for num in xrange(11,20)]

    F_list = [[0. for k_component in xrange(9)] for k_cell in xrange(n_cells)]

    for k_component in xrange(9):
        history_file = open(history_files_names[k_component], 'r')
        for line in history_file:
            if line.startswith('*') or line.startswith('$'): continue
            line = line.split()
            F_list[int(line[0])-1][k_component] = float(line[1])
        history_file.close()

    F_array = myVTK.createFloatArray(array_name, 9, n_cells)

    for k_cell in xrange(n_cells):
        F_array.InsertTuple(k_cell, F_list[k_cell])

    myVTK.myPrint(verbose, "n_tuples = " + str(F_array.GetNumberOfTuples()))

    mesh.GetCellData().AddArray(F_array)

def readAbaqusFibersFromINP(
        filename,
        verbose=1):

    myVTK.myPrint(verbose, "*** readAbaqusFibersFromINP: " + filename + " ***")

    eF_array = myVTK.createFloatArray('eF', 3)
    eS_array = myVTK.createFloatArray('eS', 3)
    eN_array = myVTK.createFloatArray('eN', 3)

    file = open(filename, 'r')
    file.readline()

    for line in file:
        line = line.split(', ')
        #print line

        eF = [float(item) for item in line[1:4]]
        eS = [float(item) for item in line[4:7]]
        eN = numpy.cross(eF,eS)
        #print "eF =", eF
        #print "eS =", eS
        #print "eN =", eN

        eF_array.InsertNextTuple(eF)
        eS_array.InsertNextTuple(eS)
        eN_array.InsertNextTuple(eN)

    file.close()

    myVTK.myPrint(verbose, "n_tuples = " + str(eF_array.GetNumberOfTuples()))

    return (eF_array,
            eS_array,
            eN_array)

def addSyntheticHelixAngles2(
        ugrid,
        angles_end,
        angles_epi,
        type_of_support="cell",
        sigma=0,
        verbose=0):

    myVTK.myPrint(verbose, "*** addSyntheticHelixAngles2 ***")

    if (type_of_support == "cell"):
        ugrid_data = ugrid.GetCellData()
    elif (type_of_support == "point"):
        ugrid_data = ugrid.GetPointData()

    farray_rr = ugrid_data.GetArray("rr")
    farray_cc = ugrid_data.GetArray("cc")
    farray_ll = ugrid_data.GetArray("ll")

    farray_angle_helix = ugrid_data.GetArray("angle_helix")

    farray_angle_helix = computeSyntheticHelixAngles2(
        farray_rr=farray_rr,
        farray_cc=farray_cc,
        farray_ll=farray_ll,
        angles_end=angles_end,
        angles_epi=angles_epi,
        sigma=sigma,
        farray_angle_helix=farray_angle_helix,
        verbose=verbose-1)

    ugrid_data.AddArray(farray_angle_helix)

    return farray_angle_helix