<html><head><meta name="color-scheme" content="light dark"></head><body><pre style="word-wrap: break-word; white-space: pre-wrap;">
#Import necessary libraries
import salome
import MEDCoupling
from MEDLoader import MEDLoader
# from MEDLoader import MEDLoader
import MEDCoupling as mc
from MEDCouplingRemapper import *
import SimpleITK as sitk
import HOMARD
import os
import unittest
import SMESH, SALOMEDS
from salome.smesh import smeshBuilder
import numpy as np

homard = salome.lcc.FindOrLoadComponent('FactoryServer', 'HOMARD')
homard.SetCurrentStudy(salome.myStudy)


def CreateImageMesh(imagePath):

    # Import the labeled image
    reader2 = sitk.ImageFileReader()
    reader2.SetFileName(imagePath)
    img_label = reader2.Execute()

    img_origin = np.array(img_label.GetOrigin())
    img_size = img_label.GetSize()
    img_spacing = np.array(img_label.GetSpacing())

    # Can only handle one image orientation (RAI code) right now (LPI - Left Posterior Inferior)
    if img_label.GetDirection()[0] != -1.0 or img_label.GetDirection()[4] != -1.0 or img_label.GetDirection()[8] != 1.0:
        print img_label.GetDirection()
        print "+++++++++++++++ Image is not in LPI orientation, please re-orient and try again. +++++++++++++++"
        return

    img_l = sitk.GetArrayFromImage(img_label)

    def build_vector(dim):
        arr = mc.DataArrayDouble(img_size[dim] + 1)
        arr.iota(0)
        if dim == 0:
            arr *= img_spacing[dim]
            arr += -img_origin[dim] - 0.5 * img_spacing[dim]
            # arr *= -1
        elif dim == 1:
            arr *= img_spacing[dim]
            arr += -img_origin[dim] - 0.5 * img_spacing[dim]
            # arr *= -1
        elif dim == 2:
            arr *= img_spacing[dim]
            arr += img_origin[dim] - 0.5 * img_spacing[dim]
        return arr

    # # Create the source mesh matching the image dimensions
    x_dim = build_vector(0)
    y_dim = build_vector(1)
    z_dim = build_vector(2)

    srcMesh = mc.MEDCouplingCMesh()
    srcMesh.setCoords(x_dim, y_dim, z_dim)
    srcMesh = srcMesh.buildUnstructured()
    srcMesh.setName("Mesh_1")

    return srcMesh, img_l, img_size

def MapImageToOriginalMesh(case_dir, imagePath, trgVolPath, MedName):

    srcMesh, img_l, img_size = CreateImageMesh(imagePath)

    # Import the target tetrahedral mesh
    trgMesh = MEDLoader.ReadUMeshFromFile(trgVolPath, 0)

    # Interpolate cells to cells
    remap = MEDCouplingRemapper()
    remap.prepare(srcMesh,trgMesh,"P0P0")

    # Source field construction
    srcField = mc.MEDCouplingFieldDouble(mc.ON_CELLS, mc.ONE_TIME)
    srcField.setMesh(srcMesh)

    array = mc.DataArrayDouble()
    array.alloc(srcField.getMesh().getNumberOfCells(),4)  # Implicitely fieldOnNodes will be a 4 component field.
    array.setInfoOnComponents(['Skin', 'Fat', 'Muscle', 'Bone'])

    for i in range(0,srcField.getMesh().getNumberOfCells(), 1):
        val = int(img_l[i//(img_size[0]*img_size[1])][i%(img_size[0]*img_size[1])//(img_size[0])][i%img_size[0]])
        if val == 4: #Skin
            array[i] = [1,0,0,0]
        elif val == 3: #Fat
            array[i] = [0,1,0,0]
        elif val == 2: #Muscle
            array[i] = [0,0,1,0]
        elif val == 1: #Bone
            array[i] = [0,0,0,1]
        else:
            array[i] = [0,0,0,0]

    srcField.setArray(array)

    # Transfer field
    # srcField.setNature(mc.Integral) # Extensive remapping - Integral (summation) of field
    srcField.setNature(mc.ConservativeVolumic) # Intensive remapping - Average of field
    trgFieldCV = remap.transferField(srcField, 0)

    srcField.setName("FIELD_src")
    trgFieldCV.setName("FIELD_trg")

    MEDLoader.WriteUMesh(os.path.join(case_dir, 'image_1.med'), srcMesh, True)
    MEDLoader.WriteFieldUsingAlreadyWrittenMesh(os.path.join(case_dir, 'image_1.med'), srcField)

    MEDLoader.WriteFieldUsingAlreadyWrittenMesh(os.path.join(case_dir, MedName), trgFieldCV)

    return trgFieldCV


def AdaptMeshFromImage(case_dir, MedName, option, i):

    srcMesh = MEDLoader.ReadUMeshFromFile(os.path.join(case_dir, 'image_1.med'), 0)
    srcField = MEDLoader.ReadFieldCell(os.path.join(case_dir, 'image_1.med'), "Mesh_1", 0, 'FIELD_src', -1, -1)


    Cases = homard.GetAllCasesName()
    if len(Cases) == 0:
        # Creation of the case
        Case_1 = homard.CreateCase(os.path.split(case_dir)[1], "Flesh", os.path.join(case_dir, MedName))
        Case_1.SetDirName(case_dir)
        Case_1.SetConfType(0)
        Case_1.SetExtType(0)
    else:
        Case_1 = homard.GetIteration("Iter_"+str(i-1))

    if option == 'FieldRefineNeighbor':
        # Creation of the hypothesis Hypo_1
        Hypo_1 = homard.CreateHypothesis("Hypo_"+str(i+1))
        Hypo_1.SetField("FIELD_trg")
        Hypo_1.SetUseField(1)  # 0: value by element (default), 1: jump between an element and neighbors
        Hypo_1.SetUseComp(0)  # 0: L2 norm (default), 1: infinite norm, 2: relative value (only 1 component)
        # Hypo_1.AddComp("Skin")
        Hypo_1.AddComp("Fat")
        Hypo_1.AddComp("Muscle")
        Hypo_1.SetRefinThr(1, 0.1)
        Hypo_1.SetUnRefThr(1, 0.001)

    elif option == 'FieldRefineValue':
        # Source field construction
        trgField = mc.MEDCouplingFieldDouble(mc.ON_CELLS, mc.ONE_TIME)
        trgField.setMesh(MEDLoader.ReadUMeshFromFile(os.path.join(case_dir, MedName[0:-5]+str(i)+'.med'), 0))

        array = mc.DataArrayDouble()
        array.alloc(trgField.getMesh().getNumberOfCells(), 1)  # Implicitely fieldOnNodes will be a 4 component field.
        array.setInfoOnComponents(['Max_comp'])

        Field_final = MEDLoader.ReadFieldCell(os.path.join(case_dir, MedName[0:-5]+str(i)+'.med'), "Flesh", 0, 'FIELD_trg', -1, -1).getArray()

        for ii in range(0, trgField.getMesh().getNumberOfCells(), 1):
            array[ii] = 1 - np.max(Field_final.getTuple(ii)[1:]/np.sum(Field_final.getTuple(ii)[1:]))

        trgField.setArray(array)
        trgField.setName("FIELD_trg_max")
        print MedName
        MEDLoader.WriteFieldUsingAlreadyWrittenMesh(os.path.join(case_dir, MedName[0:-5]+str(i)+'.med'), trgField)

        # Refine mesh using the value of the element
        Hypo_1 = homard.CreateHypothesis("Hypo_"+str(i+1))
        Hypo_1.SetField("FIELD_trg_max")
        Hypo_1.SetUseField(0)  # 0: value by element (default), 1: jump between an element and neighbors
        Hypo_1.SetUseComp(2)   # 0: L2 norm (default), 1: infinite norm, 2: relative value (only 1 component)
        Hypo_1.AddComp("Max_comp")
        Hypo_1.SetRefinThr(1, 0.3)
        Hypo_1.SetUnRefThr(1, 0.2)

    elif option == 'UniformRefine':
        # Refine Mesh uniformly
        Hypo_1 = homard.CreateHypothesis("Hypo_"+str(i+1))
        Hypo_1.SetUnifRefinUnRef(1)

    elif option == 'UnrefineNeighbor':
        # Creation of the hypothesis Hypo_1
        Hypo_1 = homard.CreateHypothesis("Hypo_"+str(i+1))
        Hypo_1.SetField("FIELD_trg")
        Hypo_1.SetUseField(1)  # 0: value by element (default), 1: jump between an element and neighbors
        Hypo_1.SetUseComp(0)  # 0: L2 norm (default), 1: infinite norm, 2: relative value (only 1 component)
        # Hypo_1.AddComp("Skin")
        Hypo_1.AddComp("Fat")
        Hypo_1.AddComp("Muscle")
        Hypo_1.SetUnRefThr(1, 0.1)

    else:
        "Not a recognized option: Choose 'FieldRefineNeighbor', 'FieldRefineValue', or 'UniformRefine'"


    Hypo_1.SetDiamMin(0.5)

    Iter_1 = Case_1.NextIteration("Iter_"+str(i))
    Iter_1.AssociateHypo("Hypo_"+str(i+1))
    Iter_1.SetMeshName('Flesh')
    Iter_1.SetMeshFile(os.path.join(case_dir, MedName[0:-5]+str(i+1)+'.med'))
    Iter_1.SetFieldFile(os.path.join(case_dir, MedName[0:-5]+str(i)+'.med'))
    codret = Iter_1.Compute(1, 1)

    # Map field on final iteration
    # Import the target tetrahedral mesh
    trgMesh = MEDLoader.ReadUMeshFromFile(os.path.join(case_dir, MedName[0:-5] + str(i+1) + '.med'), 0)

    # Interpolate cells to cells
    remap = MEDCouplingRemapper()
    remap.prepare(srcMesh, trgMesh, "P0P0")

    # Transfer field
    # srcField.setNature(mc.Integral) # Extensive remapping - Integral (summation) of field
    srcField.setNature(mc.ConservativeVolumic)  # Intensive remapping - Average of field
    trgFieldCV = remap.transferField(srcField, 0)

    trgFieldCV.setName("FIELD_trg")

    MEDLoader.WriteFieldUsingAlreadyWrittenMesh(os.path.join(case_dir, MedName[0:-5] + str(i+1) + '.med'), trgFieldCV)


def create_groups_from_field(case_dir, model_dir, MedName, i):

    smesh = smeshBuilder.New(salome.myStudy)

    ([Mesh_final], status) = smesh.CreateMeshesFromMED(os.path.join(case_dir, MedName[0:-5] + str(i) + '.med'))
    mesh = MEDLoader.ReadUMeshFromFile(os.path.join(case_dir, MedName[0:-5] + str(i) + '.med'), 0)
    Field_final = MEDLoader.ReadFieldCell(os.path.join(case_dir, MedName[0:-5] + str(i) + '.med'), "Flesh", 0, 'FIELD_trg', -1, -1)
    Field_final.setMesh(mesh)

    for group in Mesh_final.GetGroups():
        if 'Fat' in group.GetName() or 'Muscle' in group.GetName() or 'Bone' in group.GetName():
            Mesh_final.RemoveGroup(group)

    elem_bkg = []
    elem_skin = []
    elem_fat = []
    elem_muscle = []
    elem_bone = []

    nFaces = Mesh_final.NbFaces()

    for n, f in enumerate(Field_final.getArray()):
        if f[0] == 0 and f[1] == 0 and f[2] == 0 and f[3] == 0:
            elem_bkg.append(n + 1 + nFaces)
        else:
            f1 = np.array(f)
            tissue_idx = np.argmax(f1[1:])
            # if tissue_idx == 0:
            #     elem_skin.append(n + 1 + nFaces)
            if tissue_idx == 0:
                elem_fat.append(n + 1 + nFaces)
            elif tissue_idx == 1:
                elem_muscle.append(n + 1 + nFaces)
            elif tissue_idx == 2:
                elem_bone.append(n + 1 + nFaces)

    # S = Mesh_final.MakeGroupByIds("Skin", SMESH.VOLUME, elem_skin)
    F = Mesh_final.MakeGroupByIds("Fat", SMESH.VOLUME, elem_fat)
    M = Mesh_final.MakeGroupByIds("Muscle", SMESH.VOLUME, elem_muscle)
    B = Mesh_final.MakeGroupByIds("Bone", SMESH.VOLUME, elem_bone)
    # Mesh_final.MakeGroupByIds("Background", SMESH.VOLUME, elem_bkg)

    Mesh_final.ExportMED(os.path.join(case_dir, MedName[0:-5] + str(i) + '.med'))
    MEDLoader.WriteFieldUsingAlreadyWrittenMesh(os.path.join(case_dir, MedName[0:-5] + str(i) + '.med'), Field_final)
    try:
        Max_field = MEDLoader.ReadFieldCell(os.path.join(case_dir, MedName[0:-5] + str(i) + '.med'), "Flesh", 0, 'FIELD_trg_max', -1, -1)
        Max_field.setMesh(mesh)
        MEDLoader.WriteFieldUsingAlreadyWrittenMesh(os.path.join(case_dir, MedName[0:-5] + str(i) + '.med'), Max_field)
    except:
        pass
</pre></body></html>