Python myVTKPythonLibrary.createFloatArray函数代码示例

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 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 computeHelixTransverseSheetAngles(
        farray_eRR,
        farray_eCC,
        farray_eLL,
        farray_eF,
        farray_eS,
        farray_eN,
        use_new_definition=False,
        verbose=0):

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

    n_tuples = farray_eRR.GetNumberOfTuples()
    assert (farray_eCC.GetNumberOfTuples() == n_tuples)
    assert (farray_eLL.GetNumberOfTuples() == n_tuples)
    assert (farray_eF.GetNumberOfTuples() == n_tuples)
    assert (farray_eS.GetNumberOfTuples() == n_tuples)
    assert (farray_eN.GetNumberOfTuples() == n_tuples)

    farray_angle_helix = myVTK.createFloatArray("angle_helix", 1, n_tuples)
    farray_angle_trans = myVTK.createFloatArray("angle_trans", 1, n_tuples)
    farray_angle_sheet = myVTK.createFloatArray("angle_sheet", 1, n_tuples)

    eRR = numpy.empty(3)
    eCC = numpy.empty(3)
    eLL = numpy.empty(3)
    eF = numpy.empty(3)
    for k_tuple in xrange(n_tuples):
        farray_eRR.GetTuple(k_tuple, eRR)
        farray_eCC.GetTuple(k_tuple, eCC)
        farray_eLL.GetTuple(k_tuple, eLL)

        farray_eF.GetTuple(k_tuple, eF)
        eF -= numpy.dot(eF, eRR) * eRR
        eF /= numpy.linalg.norm(eF)
        angle_helix = math.copysign(1., numpy.dot(eF, eCC)) * math.asin(min(1., max(-1., numpy.dot(eF, eLL)))) * (180./math.pi)
        farray_angle_helix.SetTuple1(k_tuple, angle_helix)

        eF  = numpy.array(farray_eF.GetTuple(k_tuple))
        eF -= numpy.dot(eF, eLL) * eLL
        eF /= numpy.linalg.norm(eF)
        angle_trans = math.copysign(-1., numpy.dot(eF, eCC)) * math.asin(min(1., max(-1., numpy.dot(eF, eRR)))) * (180./math.pi)
        farray_angle_trans.SetTuple1(k_tuple, angle_trans)

        #if (use_new_definition):
            #assert 0, "TODO"
        #else:
            #assert 0, "TODO"

    return (farray_angle_helix,
            farray_angle_trans,
            farray_angle_sheet)

def createArray(
        array_name,
        n_components=1,
        n_tuples=0,
        array_type="float",
        init_to_zero=0,
        verbose=1):

    assert (type(array_type) in [type, str]), "array_type must be a type or a str. Aborting."

    if (type(array_type) is type):
        assert (array_type in [float, int]), "if a type, array_type must be equal to float or int. Aborting."

        if (array_type == float):
            return myVTK.createFloatArray(
                       name=array_name,
                       n_components=n_components,
                       n_tuples=n_tuples,
                       init_to_zero=init_to_zero,
                       verbose=verbose-1)
        elif (array_type == int):
            return myVTK.createIntArray(
                       name=array_name,
                       n_components=n_components,
                       n_tuples=n_tuples,
                       init_to_zero=init_to_zero,
                       verbose=verbose-1)
    elif (type(array_type) is str):
        assert (array_type in ["double", "float", "int", "short"]), "if a str, array_type must be equal to double, float, int or short. Aborting."

        if (array_type == "float") or (array_type == "double"):
            return myVTK.createFloatArray(
                       name=array_name,
                       n_components=n_components,
                       n_tuples=n_tuples,
                       init_to_zero=init_to_zero,
                       verbose=verbose-1)
        elif (array_type == "int"):
            return myVTK.createIntArray(
                       name=array_name,
                       n_components=n_components,
                       n_tuples=n_tuples,
                       init_to_zero=init_to_zero,
                       verbose=verbose-1)
        elif (array_type == "short"):
            return myVTK.createShortArray(
                       name=array_name,
                       n_components=n_components,
                       n_tuples=n_tuples,
                       init_to_zero=init_to_zero,
                       verbose=verbose-1)

def readAbaqusStressFromDAT(
        data_filename,
        verbose=1):

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

    s_array = myVTK.createFloatArray("", 6)

    data_file = open(data_filename, 'r')
    context = ""
    k_cell = 0
    for line in data_file:
        if (context == "reading stresses"):
            #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."
                s_list = [float(splitted_line[k]) for k in xrange(3,9)]
                s_array.InsertNextTuple(s_list)
                k_cell += 1

        if (line == "    ELEMENT  PT FOOT-       S11         S22         S33         S12         S13         S23     \n"):
            context = "reading stresses"

    data_file.close()

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

    return s_array

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 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 computeHelixAngles(
        farray_eRR,
        farray_eCC,
        farray_eLL,
        farray_eF,
        verbose=0):

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

    n_tuples = farray_eRR.GetNumberOfTuples()
    assert (farray_eCC.GetNumberOfTuples() == n_tuples)
    assert (farray_eLL.GetNumberOfTuples() == n_tuples)
    assert (farray_eF.GetNumberOfTuples() == n_tuples)

    farray_angle_helix = myVTK.createFloatArray("angle_helix", 1, n_tuples)

    eRR = numpy.empty(3)
    eCC = numpy.empty(3)
    eLL = numpy.empty(3)
    eF = numpy.empty(3)
    for k_tuple in xrange(n_tuples):
        farray_eRR.GetTuple(k_tuple, eRR)
        farray_eCC.GetTuple(k_tuple, eCC)
        farray_eLL.GetTuple(k_tuple, eLL)
        farray_eF.GetTuple(k_tuple, eF)
        eF -= numpy.dot(eF, eRR) * eRR
        eF /= numpy.linalg.norm(eF)
        helix_angle = math.copysign(1., numpy.dot(eF, eCC)) * math.asin(min(1., max(-1., numpy.dot(eF, eLL)))) * (180./math.pi)
        farray_angle_helix.SetTuple1(k_tuple, helix_angle)

    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 computeSystolicStrainsFromEndDiastolicAndEndSystolicStates(
        farray_F_dia,
        farray_F_sys,
        verbose=0):

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

    n_tuples = farray_F_dia.GetNumberOfTuples()
    assert (farray_F_sys.GetNumberOfTuples() == n_tuples)

    farray_E_dia = myVTK.createFloatArray('E_dia', 6, n_tuples)
    farray_E_sys = myVTK.createFloatArray('E_sys', 6, n_tuples)
    farray_F_num = myVTK.createFloatArray('F_num', 6, n_tuples)
    farray_E_num = myVTK.createFloatArray('E_num', 6, n_tuples)

    I = numpy.eye(3)
    E_vec = numpy.empty(6)
    for k_tuple in xrange(n_tuples):
        F_dia = numpy.reshape(farray_F_dia.GetTuple(k_tuple), (3,3), order='C')
        F_sys = numpy.reshape(farray_F_sys.GetTuple(k_tuple), (3,3), order='C')
        #print 'F_dia =', F_dia
        #print 'F_sys =', F_sys

        C = numpy.dot(numpy.transpose(F_dia), F_dia)
        E = (C - I)/2
        mat_sym33_to_vec_col6(E, E_vec)
        farray_E_dia.SetTuple(k_tuple, E_vec)

        C = numpy.dot(numpy.transpose(F_sys), F_sys)
        E = (C - I)/2
        mat_sym33_to_vec_col6(E, E_vec)
        farray_E_sys.SetTuple(k_tuple, E_vec)

        F = numpy.dot(F_sys, numpy.linalg.inv(F_dia))
        farray_F_num.SetTuple(k_tuple, numpy.reshape(F, 9, order='C'))
        #print 'F =', F

        C = numpy.dot(numpy.transpose(F), F)
        E = (C - I)/2
        mat_sym33_to_vec_col6(E, E_vec)
        farray_E_num.SetTuple(k_tuple, E_vec)

    return (farray_E_dia,
            farray_E_sys,
            farray_F_num,
            farray_E_num)

def generateWarpedMesh(
        mesh_folder,
        mesh_basename,
        images,
        structure,
        deformation,
        evolution,
        verbose=0):

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

    mesh = myVTK.readUGrid(
        filename=mesh_folder+"/"+mesh_basename+".vtk",
        verbose=verbose-1)
    n_points = mesh.GetNumberOfPoints()
    n_cells = mesh.GetNumberOfCells()

    if os.path.exists(mesh_folder+"/"+mesh_basename+"-WithLocalBasis.vtk"):
        ref_mesh = myVTK.readUGrid(
            filename=mesh_folder+"/"+mesh_basename+"-WithLocalBasis.vtk",
            verbose=verbose-1)
    else:
        ref_mesh = None

    farray_disp = myVTK.createFloatArray(
        name="displacement",
        n_components=3,
        n_tuples=n_points,
        verbose=verbose-1)
    mesh.GetPointData().AddArray(farray_disp)

    mapping = Mapping(images, structure, deformation, evolution)

    X = numpy.empty(3)
    x = numpy.empty(3)
    U = numpy.empty(3)
    if ("zfill" not in images.keys()):
        images["zfill"] = len(str(images["n_frames"]))
    for k_frame in xrange(images["n_frames"]):
        t = images["T"]*float(k_frame)/(images["n_frames"]-1) if (images["n_frames"]>1) else 0.
        mapping.init_t(t)

        for k_point in xrange(n_points):
            mesh.GetPoint(k_point, X)
            mapping.x(X, x)
            U = x - X
            farray_disp.SetTuple(k_point, U)

        myVTK.computeStrainsFromDisplacements(
            mesh=mesh,
            disp_array_name="displacement",
            ref_mesh=ref_mesh,
            verbose=verbose-1)

        myVTK.writeUGrid(
            ugrid=mesh,
            filename=mesh_folder+"/"+mesh_basename+"_"+str(k_frame).zfill(images["zfill"])+".vtk",
            verbose=verbose-1)

def computeCartesianCoordinates(
        points,
        verbose=1):

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

    n_points = points.GetNumberOfPoints()

    [xmin, xmax, ymin, ymax, zmin, zmax] = points.GetBounds()
    dx = xmax-xmin
    dy = ymax-ymin
    dz = zmax-zmin
    if (verbose >= 2): print "xmin = "+str(xmin)
    if (verbose >= 2): print "xmax = "+str(xmax)
    if (verbose >= 2): print "dx = "+str(dx)
    if (verbose >= 2): print "ymin = "+str(ymin)
    if (verbose >= 2): print "ymax = "+str(ymax)
    if (verbose >= 2): print "dy = "+str(dy)
    if (verbose >= 2): print "zmin = "+str(zmin)
    if (verbose >= 2): print "zmax = "+str(zmax)
    if (verbose >= 2): print "dz = "+str(dz)

    farray_xx = myVTK.createFloatArray("xx", 1, n_points)
    farray_yy = myVTK.createFloatArray("yy", 1, n_points)
    farray_zz = myVTK.createFloatArray("zz", 1, n_points)

    point = numpy.empty(3)
    for k_point in xrange(n_points):
        if (verbose >= 2): print "k_point = "+str(k_point)

        points.GetPoint(k_point, point)
        if (verbose >= 2): print "point = "+str(point)

        xx = (point[0] - xmin) / dx
        yy = (point[1] - ymin) / dy
        zz = (point[2] - zmin) / dz

        farray_xx.SetTuple1(k_point, xx)
        farray_yy.SetTuple1(k_point, yy)
        farray_zz.SetTuple1(k_point, zz)

    return (farray_xx,
            farray_yy,
            farray_zz)

def computeCartesianCoordinates(
        points,
        verbose=1):

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

    n_points = points.GetNumberOfPoints()

    [xmin, xmax, ymin, ymax, zmin, zmax] = points.GetBounds()
    dx = xmax-xmin
    dy = ymax-ymin
    dz = zmax-zmin
    if (verbose >= 2): print "xmin = " + str(xmin)
    if (verbose >= 2): print "xmax = " + str(xmax)
    if (verbose >= 2): print "dx = " + str(dx)
    if (verbose >= 2): print "ymin = " + str(ymin)
    if (verbose >= 2): print "ymax = " + str(ymax)
    if (verbose >= 2): print "dy = " + str(dy)
    if (verbose >= 2): print "zmin = " + str(zmin)
    if (verbose >= 2): print "zmax = " + str(zmax)
    if (verbose >= 2): print "dz = " + str(dz)

    farray_norm_x = myVTK.createFloatArray("norm_x", 1, n_points)
    farray_norm_y = myVTK.createFloatArray("norm_y", 1, n_points)
    farray_norm_z = myVTK.createFloatArray("norm_z", 1, n_points)

    for k_point in xrange(n_points):
        if (verbose >= 2): print "k_point = " + str(k_point)

        point = points.GetPoint(k_point)
        if (verbose >= 2): print "point = " + str(point)

        norm_x = (point[0] - xmin) / dx
        norm_y = (point[1] - ymin) / dy
        norm_z = (point[2] - zmin) / dz

        farray_norm_x.InsertTuple(k_point, [norm_x])
        farray_norm_y.InsertTuple(k_point, [norm_y])
        farray_norm_z.InsertTuple(k_point, [norm_z])

    return (farray_norm_x,
            farray_norm_y,
            farray_norm_z)

def computeSystolicStrains(
        farray_F_dia,
        farray_F_sys,
        verbose=1):

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

    n_tuples = farray_F_dia.GetNumberOfTuples()

    farray_E_dia = myVTK.createFloatArray('E_dia', 6, n_tuples)
    farray_E_sys = myVTK.createFloatArray('E_sys', 6, n_tuples)
    farray_F_num = myVTK.createFloatArray('F_num', 6, n_tuples)
    farray_E_num = myVTK.createFloatArray('E_num', 6, n_tuples)

    for k_tuple in xrange(n_tuples):
        F_dia = numpy.reshape(farray_F_dia.GetTuple(k_tuple), (3,3), order='C')
        F_sys = numpy.reshape(farray_F_sys.GetTuple(k_tuple), (3,3), order='C')
        #print 'F_dia =', F_dia
        #print 'F_sys =', F_sys

        C = numpy.dot(numpy.transpose(F_dia), F_dia)
        E = (C - numpy.eye(3))/2
        farray_E_dia.InsertTuple(k_tuple, mat_sym_to_vec_col(E))

        C = numpy.dot(numpy.transpose(F_sys), F_sys)
        E = (C - numpy.eye(3))/2
        farray_E_sys.InsertTuple(k_tuple, mat_sym_to_vec_col(E))

        F = numpy.dot(F_sys, numpy.linalg.inv(F_dia))
        farray_F_num.InsertTuple(k_tuple, numpy.reshape(F, 9, order='C'))
        #print 'F =', F

        C = numpy.dot(numpy.transpose(F), F)
        E = (C - numpy.eye(3))/2
        farray_E_num.InsertTuple(k_tuple, mat_sym_to_vec_col(E))

    return (farray_E_dia,
            farray_E_sys,
            farray_F_num,
            farray_E_num)

def rotateSymmetricMatrix(
        old_array,
        old_array_storage="vec",
        in_vecs=None,
        out_vecs=None,
        verbose=1):

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

    n_tuples = old_array.GetNumberOfTuples()
    new_array = myVTK.createFloatArray("", 6, n_tuples)

    for k_tuple in xrange(n_tuples):
        old_matrix = old_array.GetTuple(k_tuple)
        if (old_array_storage == "vec"):
            old_matrix = vec_col_to_mat_sym(old_matrix)
        elif (old_array_storage == "Cmat"):
            old_matrix = numpy.reshape(old_matrix, (3,3), order='C')
        elif (old_array_storage == "Fmat"):
            old_matrix = numpy.reshape(old_matrix, (3,3), order='F')

        if (in_vecs == None):
            in_R = numpy.eye(3)
        else:
            in_R = numpy.transpose(numpy.array([in_vecs[0].GetTuple(k_tuple),
                                                in_vecs[1].GetTuple(k_tuple),
                                                in_vecs[2].GetTuple(k_tuple)]))

        if (out_vecs == None):
            out_R = numpy.eye(3)
        else:
            out_R = numpy.transpose(numpy.array([out_vecs[0].GetTuple(k_tuple),
                                                 out_vecs[1].GetTuple(k_tuple),
                                                 out_vecs[2].GetTuple(k_tuple)]))

        R = numpy.dot(numpy.transpose(in_R), out_R)

        new_matrix = numpy.dot(numpy.dot(numpy.transpose(R), old_matrix), R)

        if (old_array_storage == "vec"):
            new_matrix = mat_sym_to_vec_col(new_matrix)
        elif (old_array_storage == "Cmat"):
            new_matrix = numpy.reshape(new_matrix, 9, order='C')
        elif (old_array_storage == "Fmat"):
            new_matrix = numpy.reshape(new_matrix, 9, order='F')

        new_array.InsertTuple(k_tuple, new_matrix)

    return new_array