diff --git a/src/scaffoldmaker/annotation/bone_terms.py b/src/scaffoldmaker/annotation/bone_terms.py new file mode 100644 index 00000000..97620011 --- /dev/null +++ b/src/scaffoldmaker/annotation/bone_terms.py @@ -0,0 +1,33 @@ +""" +Common resource for testing annotation terms. +""" + +# convention: preferred name, preferred id, followed by any other ids and alternative names +bone_terms = [ + #General bones + ("Cancellous bone", "ILX:0491986"), + ("Cortical bone", "ILX:0745808"), + # Geometric markers for bones + ("Olecranon", "UBERON:0006810","ILX:0735843"), + ("Coronoid process of ulna", "UBERON:0010994","ILX:0735423"), + ("Greater tubercle of radius","ILX:0745221"), + ("Medial epicondyle of radius","ILX:0748959"), + ("Lateral epicondyle of radius", "ILX:0748301"), + ("Intercondylar eminence of tibia", "ILX:0743525"), + ("Inferior articular surface of tibia", "ILX:0748099"), + ("Medial surface of tibia", "ILX:0747893"), + ("Lateral surface of tibia", "ILX:0742857"), + +] + + +def get_bone_term(name : str): + """ + Find term by matching name to any identifier held for a term. + Raise exception if name not found. + :return ( preferred name, preferred id ) + """ + for term in bone_terms: + if name in term: + return ( term[0], term[1] ) + raise NameError("Bone annotation term '" + name + "' not found.") diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_bone1.py b/src/scaffoldmaker/meshtypes/meshtype_3d_bone1.py index fd93ffcd..ec26ca8f 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_bone1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_bone1.py @@ -7,7 +7,8 @@ import copy -from cmlibs.utils.zinc.field import findOrCreateFieldCoordinates +from cmlibs.utils.zinc.field import findOrCreateFieldCoordinates, findOrCreateFieldGroup, findOrCreateFieldStoredString, \ + findOrCreateFieldStoredMeshLocation from cmlibs.utils.zinc.finiteelement import getMaximumNodeIdentifier, getMaximumElementIdentifier from cmlibs.zinc.field import Field from cmlibs.zinc.node import Node @@ -22,7 +23,7 @@ from scaffoldmaker.utils.zinc_utils import exnode_string_from_nodeset_field_parameters from scaffoldmaker.utils.derivativemoothing import DerivativeSmoothing from scaffoldmaker.annotation.annotationgroup import AnnotationGroup, findOrCreateAnnotationGroupForTerm - +from scaffoldmaker.annotation.bone_terms import get_bone_term class MeshType_3d_bone1 (Scaffold_base): """ @@ -30,7 +31,7 @@ class MeshType_3d_bone1 (Scaffold_base): with variable numbers of elements in major, minor, shell and axial directions. """ centralPathDefaultScaffoldPackages = { - 'Cylinder 1': ScaffoldPackage(MeshType_1d_path1, { + 'Bone 1': ScaffoldPackage(MeshType_1d_path1, { 'scaffoldSettings': { 'Coordinate dimensions': 3, 'D2 derivatives': True, @@ -46,21 +47,99 @@ class MeshType_3d_bone1 (Scaffold_base): (3, [[0.0, 0.0, 2.0], [0.0, 0.0, 1.0], [1.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 0.0]]), (4, [[0.0, 0.0, 3.0], [0.0, 0.0, 1.0], [1.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 0.0]]) ]) - }) + }), + 'Ulna 1': ScaffoldPackage(MeshType_1d_path1, { + 'scaffoldSettings': { + 'Coordinate dimensions': 3, + 'D2 derivatives': True, + 'D3 derivatives': True, + 'Length': 3.0, + 'Number of elements': 3 + }, + 'meshEdits': exnode_string_from_nodeset_field_parameters( + [Node.VALUE_LABEL_VALUE, Node.VALUE_LABEL_D_DS1, Node.VALUE_LABEL_D_DS2, Node.VALUE_LABEL_D2_DS1DS2, + Node.VALUE_LABEL_D_DS3, Node.VALUE_LABEL_D2_DS1DS3], [ + (1, [[ 5.356588e-01, 1.344625e+00,-1.057989e+00], [-2.179484e-01,-6.019707e-01, 8.118885e-01], [ 1.000000e+00, 0.000000e+00, 0.000000e+00], [ 0.000000e+00, 0.000000e+00, 0.000000e+00], [ 0.000000e+00, 1.000000e+00, 0.000000e+00], [ 0.000000e+00, 0.000000e+00, 0.000000e+00]]), + (2, [[ 5.342651e-01, 2.042545e-01, 7.154451e-01], [ 0.000000e+00, 0.000000e+00, 1.000000e+00], [ 1.000000e+00, 0.000000e+00, 0.000000e+00], [ 0.000000e+00, 0.000000e+00, 0.000000e+00], [ 0.000000e+00, 1.000000e+00, 0.000000e+00], [ 0.000000e+00, 0.000000e+00, 0.000000e+00]]), + (3, [[ 4.256270e-01,-7.595472e-01, 2.423991e+00], [ 0.000000e+00, 0.000000e+00, 1.000000e+00], [ 1.000000e+00, 0.000000e+00, 0.000000e+00], [ 0.000000e+00, 0.000000e+00, 0.000000e+00], [ 0.000000e+00, 1.000000e+00, 0.000000e+00], [ 0.000000e+00, 0.000000e+00, 0.000000e+00]]), + (4, [[ 8.636938e-01,-2.271568e+00, 3.615462e+00], [ 3.065797e-01,-1.318752e+00, 8.066436e-01], [ 1.000000e+00, 0.000000e+00, 0.000000e+00], [ 0.000000e+00, 0.000000e+00, 0.000000e+00], [ 0.000000e+00, 1.000000e+00, 0.000000e+00], [ 0.000000e+00, 0.000000e+00, 0.000000e+00]]) + ]) + }), + 'Radius 1': ScaffoldPackage(MeshType_1d_path1, { + 'scaffoldSettings': { + 'Coordinate dimensions': 3, + 'D2 derivatives': True, + 'D3 derivatives': True, + 'Length': 3.0, + 'Number of elements': 3 + }, + 'meshEdits': exnode_string_from_nodeset_field_parameters( + [Node.VALUE_LABEL_VALUE, Node.VALUE_LABEL_D_DS1, Node.VALUE_LABEL_D_DS2, Node.VALUE_LABEL_D2_DS1DS2, + Node.VALUE_LABEL_D_DS3, Node.VALUE_LABEL_D2_DS1DS3], [ + (1, [[-1.046980e+00,-6.643628e-01,-1.123687e+00], [ 0.000000e+00, 0.000000e+00, 1.000000e+00], [ 2.606645e-01, 3.753671e-02, 1.551701e-02], [ 0.000000e+00, 0.000000e+00, 0.000000e+00], [-6.214456e-02, 2.659024e-01, 0.000000e+00], [ 0.000000e+00, 0.000000e+00, 0.000000e+00]]), + (2, [[-9.877670e-01,-7.028367e-01,-4.943481e-01], [ 8.919182e-02,-5.293884e-02, 7.619579e-01], [ 2.444035e-01, 9.618276e-03, 1.657146e-02], [ 0.000000e+00, 0.000000e+00, 0.000000e+00], [-4.259481e-02, 2.119403e-01, 3.984784e-03], [ 0.000000e+00, 0.000000e+00, 0.000000e+00]]), + (3, [[-9.341713e-01,-6.580989e-01, 1.537421e-01], [ 1.589471e-01,-5.608755e-02, 8.047271e-01], [ 2.196727e-01, 2.198117e-02, 0.000000e+00], [ 0.000000e+00, 0.000000e+00, 0.000000e+00], [-2.156681e-02, 2.155317e-01, 0.000000e+00], [ 0.000000e+00, 0.000000e+00, 0.000000e+00]]), + (4, [[-8.515940e-01,-6.576398e-01, 6.944961e-01], [ 6.039198e-03, 5.546106e-02, 2.692810e-01], [ 2.174544e-01, 6.164221e-02, 0.000000e+00], [ 0.000000e+00, 0.000000e+00, 0.000000e+00], [-6.445803e-02, 2.273877e-01, 0.000000e+00], [ 0.000000e+00, 0.000000e+00, 0.000000e+00]]) + ]) + }), + 'Tibia 1': ScaffoldPackage(MeshType_1d_path1, { + 'scaffoldSettings': { + 'Coordinate dimensions': 3, + 'D2 derivatives': True, + 'D3 derivatives': True, + 'Length': 3.0, + 'Number of elements': 3 + }, + 'meshEdits': exnode_string_from_nodeset_field_parameters( + [Node.VALUE_LABEL_VALUE, Node.VALUE_LABEL_D_DS1, Node.VALUE_LABEL_D_DS2, Node.VALUE_LABEL_D2_DS1DS2, + Node.VALUE_LABEL_D_DS3, Node.VALUE_LABEL_D2_DS1DS3], [ + (1, [[ 5.478874e-02,-2.031468e-02,-9.043539e+00], [ 1.707559e-01,-5.112576e-02, 5.609759e+00], [ 7.943801e-01,-7.964818e-02,-7.322865e-03], [ 0.000000e+00, 0.000000e+00, 0.000000e+00], [ 9.087504e-02, 1.136004e+00, 7.587069e-03], [ 0.000000e+00, 0.000000e+00, 0.000000e+00]]), + (2, [[ 6.058305e-02,-9.069170e-02,-3.895807e+00], [-1.109146e-01, 5.163967e-02, 5.480916e+00], [ 5.653736e-01,-1.230512e-01, 1.260054e-02], [ 0.000000e+00, 0.000000e+00, 0.000000e+00], [ 1.389491e-01, 6.380904e-01,-3.200067e-03], [ 0.000000e+00, 0.000000e+00, 0.000000e+00]]), + (3, [[-1.877598e-01, 1.107110e-02, 1.514034e+00], [-4.509590e-02, 6.698397e-02, 5.530504e+00], [ 7.153731e-01,-1.353970e-01, 7.473066e-03], [ 0.000000e+00, 0.000000e+00, 0.000000e+00], [ 1.388781e-01, 7.333377e-01,-7.749571e-03], [ 0.000000e+00, 0.000000e+00, 0.000000e+00]]), + (4, [[-2.093406e-02, 4.163051e-02, 7.155221e+00], [ 3.783990e-01,-5.859759e-03, 5.746579e+00], [ 1.346095e+00,-1.457858e-01,-8.878590e-02], [ 0.000000e+00, 0.000000e+00, 0.000000e+00], [ 2.337453e-01, 2.166286e+00,-1.318264e-02], [ 0.000000e+00, 0.000000e+00, 0.000000e+00]]) + ]) + }), + } @staticmethod def getName(): return '3D Bone 1' + @staticmethod + def getParameterSetNames(): + return [ + 'Default', + 'Bone 1', + 'Ulna 1', + 'Radius 1', + 'Tibia 1', + ] @classmethod def getDefaultOptions(cls, parameterSetName='Default'): - centralPathOption = cls.centralPathDefaultScaffoldPackages['Cylinder 1'] + + + if parameterSetName == 'Default': + parameterSetName = 'Bone 1' + + if 'Bone 1' in parameterSetName: + centralPathOption = cls.centralPathDefaultScaffoldPackages['Bone 1'] + elif 'Ulna 1' in parameterSetName: + centralPathOption = cls.centralPathDefaultScaffoldPackages['Ulna 1'] + elif 'Radius 1' in parameterSetName: + centralPathOption = cls.centralPathDefaultScaffoldPackages['Radius 1'] + elif 'Tibia 1' in parameterSetName: + centralPathOption = cls.centralPathDefaultScaffoldPackages['Tibia 1'] + + + + + options = { 'Central path': copy.deepcopy(centralPathOption), - 'Number of elements across major': 6, - 'Number of elements across minor': 6, - 'Number of elements across shell': 1, + 'Number of elements across major': 8, + 'Number of elements across minor': 8, + 'Number of elements across shell': 2, 'Number of elements across transition': 1, 'Number of elements along': 4, 'Shell element thickness proportion': 1.0, @@ -74,6 +153,8 @@ def getDefaultOptions(cls, parameterSetName='Default'): 'Lower scale': 1.0, 'Lower scale_Z': 1.0 } + + options['Base parameter set'] = parameterSetName return options @staticmethod @@ -86,6 +167,7 @@ def getOrderedOptionNames(): 'Number of elements across transition', 'Number of elements along', 'Shell element thickness proportion', + # 'Lower half', 'Refine', 'Refine number of elements across major', 'Refine number of elements along', @@ -132,28 +214,39 @@ def checkOptions(cls, options): options['Central path'] = cls.getOptionScaffoldPackage('Central path', MeshType_1d_path1) dependentChanges = False - if options['Number of elements across major'] < 4: - options['Number of elements across major'] = 4 + if options['Number of elements across major'] < 6: + options['Number of elements across major'] = 6 if options['Number of elements across major'] % 2: options['Number of elements across major'] += 1 - if options['Number of elements across minor'] < 4: - options['Number of elements across minor'] = 4 - if options['Number of elements across minor'] % 2: - options['Number of elements across minor'] += 1 + #if options['Number of elements across minor'] < 6: + # options['Number of elements across minor'] = 6 + #if options['Number of elements across minor'] % 2: + # options['Number of elements across minor'] += 1 + options['Number of elements across minor'] = options['Number of elements across major' ] + #if options['Number of elements across minor'] != options['Number of elements across major']: + # options['number of elements across minor'] = options['number of elements across major'] if options['Number of elements along'] < 1: options['Number of elements along'] = 1 if options['Number of elements across transition'] < 1: options['Number of elements across transition'] = 1 + Rcrit = min(options['Number of elements across major']-4, options['Number of elements across minor']-4)//2 if options['Number of elements across shell'] + options['Number of elements across transition'] - 1 > Rcrit: dependentChanges = True options['Number of elements across shell'] = Rcrit options['Number of elements across transition'] = 1 + if options['Number of elements across major'] == 8: + options['Number of elements across shell'] = 2 if options['Shell element thickness proportion'] < 0.15: options['Shell element thickness proportion'] = 1.0 + parameterSetName = options['Base parameter set'] + isRadius = 'Radius 1' in parameterSetName + isUlna = 'Ulna 1' in parameterSetName + isTibia = 'Tibia 1' in parameterSetName + return dependentChanges @staticmethod @@ -164,7 +257,10 @@ def generateBaseMesh(region, options): :param options: Dict containing options. See getDefaultOptions(). :return: None """ - + parameterSetName = options['Base parameter set'] + isRadius = 'Radius 1' in parameterSetName + isUlna = 'Ulna 1' in parameterSetName + isTibia = 'Tibia 1' in parameterSetName centralPath = options['Central path'] full = not options['Lower half'] elementsCountAcrossMajor = options['Number of elements across major'] @@ -195,6 +291,7 @@ def generateBaseMesh(region, options): nodes = fm.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) mesh = fm.findMeshByDimension(3) + cache = fm.createFieldcache() sphere_shape = SphereShape.SPHERE_SHAPE_FULL sphere_base = cylinder1._ellipses[0] @@ -218,6 +315,7 @@ def generateBaseMesh(region, options): # get hemisphere nodes from both cylinder end and top of the sphere and mix them hemisphere._boxDerivatives = sphere1._shield3D._boxDerivatives + hemisphere._boxMapping = sphere1._shield3D._boxMapping hemisphere._box_deriv_mapping = sphere1._shield3D._box_deriv_mapping hemisphere._element_needs_scale_factor = sphere1._shield3D._element_needs_scale_factor @@ -253,7 +351,6 @@ def generateBaseMesh(region, options): if cylinder1._shield.nodeId[0][n2c][n1c]: nodesIdCylinderProximalEnd.append(cylinder1._shield.nodeId[0][n2c][n1c]) - # ****************************************************************************************************************************** # generate hemisphere extra nodes. rangeOfRequiredElements = [[0, elementsCountAcross[0]], [0, elementsCountAcross[1]], [0, int(0.5 * elementsCountAcross[2]) - 1]] @@ -278,6 +375,7 @@ def generateBaseMesh(region, options): hemisphere.nodeId[elementsCountAcross[2] // 2][n2][n1] = nodesIdCylinderProximalEnd[count] count += 1 + # generate hemisphere elements. rangeOfRequiredElements = [[0, elementsCountAcross[0]], [0, elementsCountAcross[1]], [0, int(0.5 * elementsCountAcross[2])]] @@ -378,7 +476,155 @@ def generateBaseMesh(region, options): rangeOfRequiredElements) - annotationGroup = [] + + + + cancellousGroup = AnnotationGroup(region, get_bone_term("Cancellous bone")) + + cancellousMeshGroup = cancellousGroup.getMeshGroup(mesh) + + + corticalGroup = AnnotationGroup(region, get_bone_term("Cortical bone")) + + corticalMeshGroup = corticalGroup.getMeshGroup(mesh) + annotationGroup = [cancellousGroup,corticalGroup] + + + array2 = [1,2,7,14,15,22,27,28] + array1 = [1,2,5,10,11,16,19,20] + + More_elements = array1 if elementsCountAcrossMajor == 6 else array2 if elementsCountAcrossMajor == 8 else None + addition = elementsCountAcrossMajor-2 if elementsCountAcrossMajor == 6 else elementsCountAcrossMajor if elementsCountAcrossMajor == 8 else None + + + + number_of_elements_per_layer = max(More_elements) #(elementsCountAcrossMajor-1)*4 + + outer_element_trunk = [x+number_of_elements_per_layer*i for i in range(0,elementsCountAlong) for x in More_elements] + print(outer_element_trunk) + max_value = max(outer_element_trunk) + dome1=[max_value + i for i in range(1,5)] + print(dome1) + dome1_final = dome1 + [x+max(dome1)+addition for x in More_elements] + print(dome1_final) + max_value = max(dome1_final) + + dome2 = [x + max_value for x in More_elements] + max_value = max(dome2) + print(dome2) + dome2_final = dome2 + [max_value + i+addition for i in range(1,5)] + print(dome2_final) + final_outer = outer_element_trunk+dome1_final+dome2_final + + + final_inner = [] + + for number in range (1,max(final_outer)-1): + if number not in final_outer: + final_inner.append(number) + + + + for key in final_inner: + element = mesh.findElementByIdentifier(int(key)) + cancellousMeshGroup.addElement(element) + + for key in final_outer: + element = mesh.findElementByIdentifier(int(key)) + corticalMeshGroup.addElement(element) + + top_elem = max(final_outer) + bottom_elem = number_of_elements_per_layer*elementsCountAlong+1 + lateral_elem = int(number_of_elements_per_layer*(elementsCountAlong/2-1)+More_elements[3]) + medial_elem = int(number_of_elements_per_layer*(elementsCountAlong/2-1)+More_elements[4]) + + + # markers with element number and xi position + allMarkersRadius = {"Greater tuberosity of radius": {"elementID": top_elem, "xi": [1.0, 0.0, 0.5]}, + "Lateral epicondyle of radius": {"elementID": bottom_elem, "xi": [1.0, 0.0, 1.0]}, + "Medial epicondyle of radius": {"elementID": medial_elem, "xi": [0.0, 0.0, 0.0]}, + "Epiphysial plate of radius": {"elementID": lateral_elem, "xi": [0.0, 0.0, 0.0]}, + } + + allMarkersUlna = {"Greater tuberosity of radius": {"elementID": top_elem, "xi": [1.0, 0.0, 0.5]}, + "Lateral epicondyle of radius": {"elementID": bottom_elem, "xi": [1.0, 0.0, 1.0]}, + } + + allMarkersTibia = {"Intercondylar eminence of tibia": {"elementID": top_elem, "xi": [0.0,1.0,0.0]}, + "Inferior articular surface of tibia": {"elementID": bottom_elem, "xi": [1.0,0.0,1.0]}, + "Medial surface of tibia": {"elementID": medial_elem, "xi": [1.0,1.0,1.0]}, + "Lateral surface of tibia": {"elementID": lateral_elem, "xi":[1.0,1.0,1.0]}, + + } + + markerGroup = findOrCreateFieldGroup(fm, "marker") + markerName = findOrCreateFieldStoredString(fm, name="marker_name") + markerLocation = findOrCreateFieldStoredMeshLocation(fm, mesh, name="marker_location") + + markerPoints = markerGroup.getOrCreateNodesetGroup(nodes) + markerTemplateInternal = nodes.createNodetemplate() + markerTemplateInternal.defineField(markerName) + markerTemplateInternal.defineField(markerLocation) + nodes = fm.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) + nodeIdentifier = max(1, getMaximumNodeIdentifier(nodes) + 1)+10 + + #annotation fiducial point + if isRadius: + for key in allMarkersRadius: + xi = allMarkersRadius[key]["xi"] + addMarker = {"name": key, "xi": allMarkersRadius[key]["xi"]} + markerPoint = markerPoints.createNode(nodeIdentifier, markerTemplateInternal) + nodeIdentifier += 1 + cache.setNode(markerPoint) + markerName.assignString(cache, addMarker["name"]) + elementID = allMarkersRadius[key]["elementID"] + element = mesh.findElementByIdentifier(elementID) + markerLocation.assignMeshLocation(cache, element, addMarker["xi"]) + + elif isUlna: + for key in allMarkersUlna: + xi = allMarkersUlna[key]["xi"] + addMarker = {"name": key, "xi": allMarkersUlna[key]["xi"]} + markerPoint = markerPoints.createNode(nodeIdentifier, markerTemplateInternal) + nodeIdentifier += 1 + cache.setNode(markerPoint) + markerName.assignString(cache, addMarker["name"]) + elementID = allMarkersRadius[key]["elementID"] + element = mesh.findElementByIdentifier(elementID) + markerLocation.assignMeshLocation(cache, element, addMarker["xi"]) + + elif isTibia: + for key in allMarkersTibia: + xi = allMarkersTibia[key]["xi"] + addMarker = {"name": key, "xi": allMarkersTibia[key]["xi"]} + markerPoint = markerPoints.createNode(nodeIdentifier, markerTemplateInternal) + nodeIdentifier += 1 + cache.setNode(markerPoint) + markerName.assignString(cache, addMarker["name"]) + elementID = allMarkersTibia[key]["elementID"] + element = mesh.findElementByIdentifier(elementID) + markerLocation.assignMeshLocation(cache, element, addMarker["xi"]) + + +# # AnnotationGroup.createMarkerNode(startNodeIdentifier=1, materialCoordinatesField: FieldFiniteElement = None, materialCoordinates = None, element = None, xi = [0.0, 0.0, 0.0]) +# +# +# # markerTermNameBoneCoordinatesMap = { +# # 'tibial tuberosity': [-5.076472492200136e+01, -4.592226612078402e+01, -1.261953033704384e+03], +# # } +# # annotationGroups = [] +# # nodes = fm.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) +# # nodeIdentifier = max(1, getMaximumNodeIdentifier(nodes) + 1) +# # print(nodeIdentifier) +# # for termName, boneCoordinatesValues in markerTermNameBoneCoordinatesMap.items(): +# # annotationGroup = findOrCreateAnnotationGroupForTerm(annotationGroups, region, ('tibia point', 'FM:63')) +# # annotationGroup.createMarkerNode(nodeIdentifier, coordinates, boneCoordinatesValues) +# # nodeIdentifier += 1 + + # smoothing = DerivativeSmoothing(region, coordinates) + # smoothing.smooth(True) + + #annotationGroup = [] return annotationGroup, None @classmethod