diff --git a/src/scaffoldmaker/annotation/body_terms.py b/src/scaffoldmaker/annotation/body_terms.py index 80de066b..03c32809 100644 --- a/src/scaffoldmaker/annotation/body_terms.py +++ b/src/scaffoldmaker/annotation/body_terms.py @@ -10,8 +10,18 @@ ("abdominopelvic cavity", "UBERON:0035819"), ("upper limb", "UBERON:0001460"), ("left upper limb", "UBERON:8300002", "FMA:7186"), + ("left shoulder", ""), + ("left brachium", ""), + ("left elbow", ""), + ("left antebrachium", ""), + ("left hand", ""), ("left upper limb skin epidermis outer surface", "ILX:0796504"), ("right upper limb", "UBERON:8300001", "FMA:7185"), + ("right shoulder", ""), + ("right brachium", ""), + ("right elbow", ""), + ("right antebrachium", ""), + ("right hand", ""), ("right upper limb skin epidermis outer surface", "ILX:0796503"), ("body", "UBERON:0000468", "ILX:0101370"), ("core", ""), @@ -21,11 +31,18 @@ ("head core", ""), ("diaphragm", "UBERON:0001103", "ILX:0103194"), ("hand", "ILX:0104885", "FMA:9712"), + ("hip", ""), ("left", ""), ("lower limb", "UBERON:0000978"), ("left lower limb", "UBERON:8300004", "FMA:24981"), + ("left upper leg", ""), + ("left lower leg", ""), + ("left foot", ""), ("left lower limb skin epidermis outer surface", "ILX:0796506"), - ("right lower limb ", "UBERON:8300003", "FMA:24980"), + ("right lower limb", "UBERON:8300003", "FMA:24980"), + ("right upper leg", ""), + ("right lower leg", ""), + ("right foot", ""), ("right lower limb skin epidermis outer surface", "ILX:0796505"), ("foot", "ILX:0745450", "FMA:9664"), ("neck", "UBERON:0000974", "ILX:0733967"), @@ -38,7 +55,34 @@ ("thoracic cavity", "UBERON:0002224"), ("thoracic cavity boundary surface", "ILX:0796508"), ("thorax", "ILX:0742178"), - ("ventral", "") + ("ventral", ""), + # kinematic tree markers + ('pelvis', ""), + ('femur_r', ""), + ('tibia_r', ""), + ('talus_r', ""), + ('calcn_r', ""), + ('toes_r', ""), + ('femur_l', ""), + ('tibia_l', ""), + ('talus_l', ""), + ('calcn_l', ""), + ('toes_l', ""), + ('lumbar_body', ""), + ('thorax_top', ""), + ('head_marker', ""), + ('scapula_r', ""), + ('humerus_r', ""), + ('ulna_r', ""), + ('radius_r', ""), + ('hand_r', ""), + ('scapula_l', ""), + ('humerus_l', ""), + ('ulna_l', ""), + ('radius_l', ""), + ('hand_l', ""), + + ] def get_body_term(name : str): @@ -51,3 +95,12 @@ def get_body_term(name : str): if name in term: return ( term[0], term[1] ) raise NameError("Body annotation term '" + name + "' not found.") + +def marker_name_in_terms(name: str): + """ + Check if term exists in approved marker terms + """ + for term in body_terms: + if name in term: + return True + return False diff --git a/src/scaffoldmaker/meshtypes/meshtype_3d_wholebody2.py b/src/scaffoldmaker/meshtypes/meshtype_3d_wholebody2.py index 374d7d9e..a8e55578 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_3d_wholebody2.py +++ b/src/scaffoldmaker/meshtypes/meshtype_3d_wholebody2.py @@ -1,13 +1,18 @@ """ Generates a 3D body coordinates using tube network mesh. """ -from cmlibs.maths.vectorops import add, cross, mult, set_magnitude, sub -from cmlibs.utils.zinc.field import Field, find_or_create_field_coordinates +from cmlibs.maths.vectorops import add, cross, mult, set_magnitude, sub, magnitude, axis_angle_to_rotation_matrix, matrix_vector_mult, matrix_mult, dot, angle +from cmlibs.utils.zinc.general import ChangeManager +from cmlibs.utils.zinc.field import Field, find_or_create_field_coordinates, find_or_create_field_group, find_or_create_field_stored_string, find_or_create_field_finite_element +from cmlibs.utils.zinc.finiteelement import get_maximum_node_identifier +from cmlibs.utils.zinc.region import copy_fitting_data from cmlibs.zinc.element import Element from cmlibs.zinc.node import Node +from scaffoldfitter.fitter import Fitter as GeometryFitter +from scaffoldfitter.fitterstepalign import FitterStepAlign from scaffoldmaker.annotation.annotationgroup import ( - AnnotationGroup, findOrCreateAnnotationGroupForTerm, getAnnotationGroupForTerm) -from scaffoldmaker.annotation.body_terms import get_body_term + AnnotationGroup, findOrCreateAnnotationGroupForTerm, getAnnotationGroupForTerm, evaluateAnnotationMarkerNearestMeshLocation) +from scaffoldmaker.annotation.body_terms import get_body_term, marker_name_in_terms from scaffoldmaker.meshtypes.meshtype_1d_network_layout1 import MeshType_1d_network_layout1 from scaffoldmaker.meshtypes.scaffold_base import Scaffold_base from scaffoldmaker.scaffoldpackage import ScaffoldPackage @@ -16,8 +21,12 @@ sampleCubicHermiteCurvesSmooth, smoothCubicHermiteDerivativesLine) from scaffoldmaker.utils.networkmesh import NetworkMesh from scaffoldmaker.utils.tubenetworkmesh import BodyTubeNetworkMeshBuilder, TubeNetworkMeshGenerateData +from scaffoldmaker.utils.human_network_layout import constructNetworkLayoutStructure, humanElementCounts +from scaffoldmaker.utils.zinc_utils import generate_datapoints, find_or_create_field_zero_fibres +# from scaffoldmaker.utils.human_network_layout import import math - +import logging +logger = logging.getLogger(__name__) class MeshType_1d_human_body_network_layout1(MeshType_1d_network_layout1): """ @@ -36,47 +45,58 @@ def getParameterSetNames(cls): def getDefaultOptions(cls, parameterSetName="Default"): options = {} options["Base parameter set"] = parameterSetName - options["Structure"] = ( - "1-2-3-4," - "4-5-6.1," - "6.2-14-15-16-17-18-19,19-20," - "6.3-21-22-23-24-25-26,26-27," - "6.1-7-8-9," - "9-10-11-12-13.1," - "13.2-28-29-30-31-32,32-33-34," - "13.3-35-36-37-38-39,39-40-41") + options["Structure"] = constructNetworkLayoutStructure(humanElementCounts) options["Define inner coordinates"] = True options["Head depth"] = 2.0 options["Head length"] = 2.2 options["Head width"] = 2.0 options["Neck length"] = 1.3 options["Shoulder drop"] = 1.0 - options["Shoulder width"] = 4.5 - options["Arm lateral angle degrees"] = 10.0 + options["Shoulder width"] = 5.0 + options["Left shoulder flexion degrees"] = 0.0 + options["Right shoulder flexion degrees"] = 0.0 + options["Left shoulder abduction degrees"] = 10.0 + options["Right shoulder abduction degrees"] = 10.0 + options["Left elbow flexion degrees"] = 0.0 + options["Right elbow flexion degrees"] = 0.0 + options["Left arm rotation degrees"] = 0.0 + options["Right arm rotation degrees"] = 0.0 options["Arm length"] = 7.5 options["Arm top diameter"] = 1.0 options["Arm twist angle degrees"] = 0.0 options["Wrist thickness"] = 0.5 options["Wrist width"] = 0.7 + options["Left wrist flexion degrees"] = 0.0 + options["Right wrist flexion degrees"] = 0.0 + options["Left wrist deviation degrees"] = 0.0 + options["Right wrist deviation degrees"] = 0.0 options["Hand length"] = 2.0 options["Hand thickness"] = 0.2 options["Hand width"] = 1.0 - options["Thorax length"] = 2.5 + options["Thorax length"] = 2.5 options["Abdomen length"] = 3.0 options["Torso depth"] = 2.5 options["Torso width"] = 3.2 - options["Pelvis drop"] = 1.5 + options["Pelvis drop"] = 0.8 options["Pelvis width"] = 2.0 - options["Leg lateral angle degrees"] = 10.0 - options["Leg length"] = 10.0 + options["Left leg abduction degrees"] = 10.0 + options["Right leg abduction degrees"] = 10.0 + options["Left hip flexion degrees"] = 0.0 + options["Right hip flexion degrees"] = 0.0 + options["Leg length"] = 11.0 options["Leg top diameter"] = 2.0 options["Leg bottom diameter"] = 0.7 + options["Left knee flexion degrees"] = 0.0 + options["Right knee flexion degrees"] = 0.0 + options["Left ankle flexion degrees"] = 90.0 + options["Right ankle flexion degrees"] = 90.0 options["Foot height"] = 1.25 - options["Foot length"] = 2.5 + options["Foot length"] = 2.0 options["Foot thickness"] = 0.3 options["Foot width"] = 1.0 options["Inner proportion default"] = 0.7 options["Inner proportion head"] = 0.35 + options['scale'] = [1,1,1] return options @classmethod @@ -88,7 +108,18 @@ def getOrderedOptionNames(cls): "Neck length", "Shoulder drop", "Shoulder width", - "Arm lateral angle degrees", + "Left shoulder abduction degrees", + "Right shoulder abduction degrees", + "Left shoulder flexion degrees", + "Right shoulder flexion degrees", + "Left arm rotation degrees", + "Right arm rotation degrees", + "Left elbow flexion degrees", + "Right elbow flexion degrees", + "Left wrist flexion degrees", + "Right wrist flexion degrees", + "Left wrist deviation degrees", + "Right wrist deviation degrees", "Arm length", "Arm top diameter", "Arm twist angle degrees", @@ -103,14 +134,21 @@ def getOrderedOptionNames(cls): "Torso width", "Pelvis drop", "Pelvis width", - "Leg lateral angle degrees", + "Left leg abduction degrees", + "Right leg abduction degrees", + "Left hip flexion degrees", + "Right hip flexion degrees", "Leg length", "Leg top diameter", "Leg bottom diameter", + "Left knee flexion degrees", + "Right knee flexion degrees", "Foot height", "Foot length", "Foot thickness", "Foot width", + "Left ankle flexion degrees", + "Right ankle flexion degrees", "Inner proportion default", "Inner proportion head" ] @@ -157,9 +195,25 @@ def checkOptions(cls, options): elif options[key] > 0.9: options[key] = 0.9 for key, angleRange in { - "Arm lateral angle degrees": (-60.0, 200.0), + "Left shoulder abduction degrees": (-100.0, 180.0), + "Right shoulder abduction degrees": (-100.0, 180.0), + "Left shoulder flexion degrees": (-60.0, 200.0), + "Right shoulder flexion degrees": (-60.0, 200.0), + "Left elbow flexion degrees": (0.0, 120.0), + "Right elbow flexion degrees": (0.0, 120.0), + "Left arm rotation degrees": (-90.0, 90.0), + "Right arm rotation degrees": (-90.0, 90.0), + # "Left wrist flexion degrees": (-30.0, 30.0), + # "Right wrist flexion degrees": (-30.0, 30.0), + "Left hip flexion degrees": (0.0, 150.0), + "Right hip flexion degrees": (0.0, 150.0), + "Left knee flexion degrees": (0.0, 140.0), + "Right knee flexion degrees": (0.0, 140.0), + "Left ankle flexion degrees": (60.0, 140.0), + "Right ankle flexion degrees": (60.0, 140.0), "Arm twist angle degrees": (-90.0, 90.0), - "Leg lateral angle degrees": (-20.0, 60.0) + "Left leg abduction degrees": (-20.0, 60.0), + "Right leg abduction degrees": (-20.0, 60.0) }.items(): if options[key] < angleRange[0]: options[key] = angleRange[0] @@ -183,7 +237,18 @@ def generateBaseMesh(cls, region, options): neckLength = options["Neck length"] shoulderDrop = options["Shoulder drop"] halfShoulderWidth = 0.5 * options["Shoulder width"] - armAngleRadians = math.radians(options["Arm lateral angle degrees"]) + shoulderLeftFlexionRadians = math.radians(options["Left shoulder flexion degrees"]) + shoulderRightFlexionRadians = math.radians(options["Right shoulder flexion degrees"]) + armLeftAngleRadians = math.radians(options["Left shoulder abduction degrees"]) + armRightAngleRadians = math.radians(options["Right shoulder abduction degrees"]) + elbowLeftFlexionRadians = math.radians(options["Left elbow flexion degrees"]) + elbowRightFlexionRadians = math.radians(options["Right elbow flexion degrees"]) + armLeftRotationRadians = math.radians(options["Left arm rotation degrees"]) + armRightRotationRadians = math.radians(options["Right arm rotation degrees"]) + wristLeftFlexionRadians = math.radians(options["Left wrist flexion degrees"]) + wristRightFlexionRadians = math.radians(options["Right wrist flexion degrees"]) + wristLeftAbductionRadians = math.radians(options["Left wrist deviation degrees"]) + wristRightAbductionRadians = math.radians(options["Right wrist deviation degrees"]) armLength = options["Arm length"] armTopRadius = 0.5 * options["Arm top diameter"] armTwistAngleRadians = math.radians(options["Arm twist angle degrees"]) @@ -198,53 +263,80 @@ def generateBaseMesh(cls, region, options): halfTorsoWidth = 0.5 * options["Torso width"] pelvisDrop = options["Pelvis drop"] halfPelvisWidth = 0.5 * options["Pelvis width"] - legAngleRadians = math.radians(options["Leg lateral angle degrees"]) + leftLegAbductionRadians = math.radians(options["Left leg abduction degrees"]) + rightLegAbductionRadians = math.radians(options["Right leg abduction degrees"]) + hipLeftFlexionRadians = math.radians(options["Left hip flexion degrees"]) + hipRightFlexionRadians = math.radians(options["Right hip flexion degrees"]) legLength = options["Leg length"] legTopRadius = 0.5 * options["Leg top diameter"] legBottomRadius = 0.5 * options["Leg bottom diameter"] + kneeLeftFlexionRadians = math.radians(options["Left knee flexion degrees"]) + kneeRightFlexionRadians = math.radians(options["Right knee flexion degrees"]) + ankleLeftFlexionRadians = math.radians( options["Left ankle flexion degrees"]) + ankleRightFlexionRadians = math.radians(options["Right ankle flexion degrees"]) footHeight = options["Foot height"] footLength = options["Foot length"] halfFootThickness = 0.5 * options["Foot thickness"] halfFootWidth = 0.5 * options["Foot width"] innerProportionDefault = options["Inner proportion default"] innerProportionHead = options["Inner proportion head"] - + # Store coordinates for kinematic tree markers + options['Kinematic tree'] = {} networkMesh = NetworkMesh(structure) networkMesh.create1DLayoutMesh(region) - fieldmodule = region.getFieldmodule() mesh = fieldmodule.findMeshByDimension(1) - # set up element annotations bodyGroup = AnnotationGroup(region, get_body_term("body")) headGroup = AnnotationGroup(region, get_body_term("head")) neckGroup = AnnotationGroup(region, get_body_term("neck")) armGroup = AnnotationGroup(region, get_body_term("upper limb")) - armToHandGroup = AnnotationGroup(region, ("arm to hand", "")) + # armToHandGroup = AnnotationGroup(region, ("arm to hand", "")) leftArmGroup = AnnotationGroup(region, get_body_term("left upper limb")) + leftShoulderGroup = AnnotationGroup(region, get_body_term("left shoulder")) + leftBrachiumGroup = AnnotationGroup(region, get_body_term("left brachium")) + leftAntebrachiumGroup = AnnotationGroup(region, get_body_term("left antebrachium")) + # leftElbowGroup = AnnotationGroup(region, get_body_term("left elbow")) + leftHandGroup = AnnotationGroup(region, get_body_term("left hand")) rightArmGroup = AnnotationGroup(region, get_body_term("right upper limb")) + rightShoulderGroup = AnnotationGroup(region, get_body_term("right shoulder")) + rightBrachiumGroup = AnnotationGroup(region, get_body_term("right brachium")) + rightAntebrachiumGroup = AnnotationGroup(region, get_body_term("right antebrachium")) + # rightElbowGroup = AnnotationGroup(region, get_body_term("right elbow")) + rightHandGroup = AnnotationGroup(region, get_body_term("right hand")) handGroup = AnnotationGroup(region, get_body_term("hand")) thoraxGroup = AnnotationGroup(region, get_body_term("thorax")) abdomenGroup = AnnotationGroup(region, get_body_term("abdomen")) + hipGroup = AnnotationGroup(region, get_body_term("hip")) legGroup = AnnotationGroup(region, get_body_term("lower limb")) legToFootGroup = AnnotationGroup(region, ("leg to foot", "")) leftLegGroup = AnnotationGroup(region, get_body_term("left lower limb")) - rightLegGroup = AnnotationGroup(region, get_body_term("right lower limb ")) + leftUpperLegGroup = AnnotationGroup(region, get_body_term("left upper leg")) + leftLowerLegGroup = AnnotationGroup(region, get_body_term("left lower leg")) + leftFootGroup = AnnotationGroup(region, get_body_term("left foot")) + rightLegGroup = AnnotationGroup(region, get_body_term("right lower limb")) + rightUpperLegGroup = AnnotationGroup(region, get_body_term("right upper leg")) + rightLowerLegGroup = AnnotationGroup(region, get_body_term("right lower leg")) + rightFootGroup = AnnotationGroup(region, get_body_term("right foot")) footGroup = AnnotationGroup(region, get_body_term("foot")) annotationGroups = [bodyGroup, headGroup, neckGroup, - armGroup, armToHandGroup, leftArmGroup, rightArmGroup, handGroup, - thoraxGroup, abdomenGroup, - legGroup, legToFootGroup, leftLegGroup, rightLegGroup, footGroup] + thoraxGroup, abdomenGroup, hipGroup, + leftShoulderGroup, leftBrachiumGroup, leftAntebrachiumGroup, leftHandGroup, + rightShoulderGroup, rightBrachiumGroup, rightAntebrachiumGroup, rightHandGroup, + leftLegGroup, leftUpperLegGroup, leftLowerLegGroup, leftFootGroup, + rightLegGroup, rightUpperLegGroup, rightLowerLegGroup, rightFootGroup, + armGroup, leftArmGroup, rightArmGroup, handGroup, + legGroup, footGroup] bodyMeshGroup = bodyGroup.getMeshGroup(mesh) elementIdentifier = 1 - headElementsCount = 3 + headElementsCount = humanElementCounts['headElementsCount'] meshGroups = [bodyMeshGroup, headGroup.getMeshGroup(mesh)] for e in range(headElementsCount): element = mesh.findElementByIdentifier(elementIdentifier) for meshGroup in meshGroups: meshGroup.addElement(element) elementIdentifier += 1 - neckElementsCount = 2 + neckElementsCount = humanElementCounts['neckElementsCount'] meshGroups = [bodyMeshGroup, neckGroup.getMeshGroup(mesh)] for e in range(neckElementsCount): element = mesh.findElementByIdentifier(elementIdentifier) @@ -253,53 +345,113 @@ def generateBaseMesh(cls, region, options): elementIdentifier += 1 left = 0 right = 1 - armToHandElementsCount = 6 - handElementsCount = 1 + shoulderElementsCount = humanElementCounts['shoulderElementsCount'] + brachiumElementsCount = humanElementCounts['brachiumElementsCount'] + antebrachiumElementsCount = humanElementCounts['antebrachiumElementsCount'] + handElementsCount = humanElementCounts['handElementsCount'] + armToHandElementsCount = shoulderElementsCount + brachiumElementsCount + antebrachiumElementsCount armMeshGroup = armGroup.getMeshGroup(mesh) - armToHandMeshGroup = armToHandGroup.getMeshGroup(mesh) + # armToHandMeshGroup = armToHandGroup.getMeshGroup(mesh) handMeshGroup = handGroup.getMeshGroup(mesh) for side in (left, right): sideArmGroup = leftArmGroup if (side == left) else rightArmGroup - meshGroups = [bodyMeshGroup, armMeshGroup, armToHandMeshGroup, sideArmGroup.getMeshGroup(mesh)] - for e in range(armToHandElementsCount): + sideShoulderGroup = leftShoulderGroup if (side == left) else rightShoulderGroup + sideBrachiumGroup = leftBrachiumGroup if (side == left) else rightBrachiumGroup + sideAntebrachiumGroup = leftAntebrachiumGroup if (side == left) else rightAntebrachiumGroup + # sideElbowGroup = leftElbowGroup if (side == left) else rightElbowGroup + sideHandGroup = leftHandGroup if (side == left) else rightHandGroup + # Setup shoulder elements + meshGroups = [bodyMeshGroup, + armMeshGroup, + sideArmGroup.getMeshGroup(mesh), sideShoulderGroup.getMeshGroup(mesh)] + for e in range(shoulderElementsCount): element = mesh.findElementByIdentifier(elementIdentifier) for meshGroup in meshGroups: meshGroup.addElement(element) elementIdentifier += 1 - meshGroups = [bodyMeshGroup, armMeshGroup, handMeshGroup, sideArmGroup.getMeshGroup(mesh)] + # Setup brachium elements + meshGroups = [bodyMeshGroup, + armMeshGroup, + sideArmGroup.getMeshGroup(mesh), sideBrachiumGroup.getMeshGroup(mesh)] + for e in range(brachiumElementsCount): + element = mesh.findElementByIdentifier(elementIdentifier) + for meshGroup in meshGroups: + meshGroup.addElement(element) + elementIdentifier += 1 + # Setup antebrachium elements + meshGroups = [bodyMeshGroup, + armMeshGroup, + sideArmGroup.getMeshGroup(mesh), sideAntebrachiumGroup.getMeshGroup(mesh)] + for e in range(antebrachiumElementsCount): + element = mesh.findElementByIdentifier(elementIdentifier) + for meshGroup in meshGroups: + meshGroup.addElement(element) + elementIdentifier += 1 + # Setup hand elements + meshGroups = [bodyMeshGroup, + armMeshGroup, + handMeshGroup, sideHandGroup.getMeshGroup(mesh)] for e in range(handElementsCount): element = mesh.findElementByIdentifier(elementIdentifier) for meshGroup in meshGroups: meshGroup.addElement(element) elementIdentifier += 1 - thoraxElementsCount = 3 - abdomenElementsCount = 4 + thoraxElementsCount = humanElementCounts['thoraxElementsCount'] + abdomenElementsCount = humanElementCounts['abdomenElementsCount'] + # Setup thorax elements meshGroups = [bodyMeshGroup, thoraxGroup.getMeshGroup(mesh)] for e in range(thoraxElementsCount): element = mesh.findElementByIdentifier(elementIdentifier) for meshGroup in meshGroups: meshGroup.addElement(element) elementIdentifier += 1 + # Setup abdomen elements meshGroups = [bodyMeshGroup, abdomenGroup.getMeshGroup(mesh)] for e in range(abdomenElementsCount): element = mesh.findElementByIdentifier(elementIdentifier) for meshGroup in meshGroups: meshGroup.addElement(element) elementIdentifier += 1 - legToFootElementsCount = 5 - footElementsCount = 2 + hipElementsCount = humanElementCounts['hipElementsCount'] + upperLegElementsCount = humanElementCounts['upperLegElementsCount'] + lowerLegElementsCount = humanElementCounts['lowerLegElementsCount'] + footElementsCount = humanElementCounts['footElementsCount'] + legToFootElementsCount = hipElementsCount + upperLegElementsCount + lowerLegElementsCount legMeshGroup = legGroup.getMeshGroup(mesh) + hipMeshGroup = hipGroup.getMeshGroup(mesh) legToFootMeshGroup = legToFootGroup.getMeshGroup(mesh) footMeshGroup = footGroup.getMeshGroup(mesh) for side in (left, right): sideLegGroup = leftLegGroup if (side == left) else rightLegGroup - meshGroups = [bodyMeshGroup, legMeshGroup, legToFootMeshGroup, sideLegGroup.getMeshGroup(mesh)] - for e in range(legToFootElementsCount): + sideUpperLegGroup = leftUpperLegGroup if (side == left) else rightUpperLegGroup + sideLowerLegGroup = leftLowerLegGroup if (side == left) else rightLowerLegGroup + sideFootGroup = leftFootGroup if (side == left) else rightFootGroup + # Hip + meshGroups = [bodyMeshGroup, legMeshGroup, hipMeshGroup, + sideLegGroup.getMeshGroup(mesh), sideUpperLegGroup.getMeshGroup(mesh)] + for e in range(hipElementsCount): + element = mesh.findElementByIdentifier(elementIdentifier) + for meshGroup in meshGroups: + meshGroup.addElement(element) + elementIdentifier += 1 + # Upper leg + meshGroups = [bodyMeshGroup, legMeshGroup, + sideLegGroup.getMeshGroup(mesh), sideUpperLegGroup.getMeshGroup(mesh)] + for e in range(upperLegElementsCount): + element = mesh.findElementByIdentifier(elementIdentifier) + for meshGroup in meshGroups: + meshGroup.addElement(element) + elementIdentifier += 1 + # Lower leg + meshGroups = [bodyMeshGroup, legMeshGroup, + sideLegGroup.getMeshGroup(mesh), sideLowerLegGroup.getMeshGroup(mesh)] + for e in range(lowerLegElementsCount): element = mesh.findElementByIdentifier(elementIdentifier) for meshGroup in meshGroups: meshGroup.addElement(element) elementIdentifier += 1 - meshGroups = [bodyMeshGroup, legMeshGroup, footMeshGroup, sideLegGroup.getMeshGroup(mesh)] + # Foot + meshGroups = [bodyMeshGroup, legMeshGroup, footMeshGroup, sideFootGroup.getMeshGroup(mesh)] for e in range(footElementsCount): element = mesh.findElementByIdentifier(elementIdentifier) for meshGroup in meshGroups: @@ -343,10 +495,12 @@ def generateBaseMesh(cls, region, options): setNodeFieldParameters(innerCoordinates, fieldcache, x, d1, id2, id3) nodeIdentifier += 1 armJunctionNodeIdentifier = nodeIdentifier + thoraxScale = thoraxLength / thoraxElementsCount thoraxStartX = headLength + neckLength sx = [thoraxStartX, 0.0, 0.0] + options['Kinematic tree']['thorax_top'] = sx for i in range(thoraxElementsCount): node = nodes.findNodeByIdentifier(nodeIdentifier) fieldcache.setNode(node) @@ -370,7 +524,7 @@ def generateBaseMesh(cls, region, options): setNodeFieldParameters(coordinates, fieldcache, x, d1, d2, d3, d12) setNodeFieldParameters(innerCoordinates, fieldcache, x, d1, id2, id3, id12) nodeIdentifier += 1 - + abdomenScale = abdomenLength / abdomenElementsCount d2 = [0.0, halfTorsoWidth, 0.0] d3 = [0.0, 0.0, halfTorsoDepth] @@ -387,37 +541,44 @@ def generateBaseMesh(cls, region, options): nodeIdentifier += 1 legJunctionNodeIdentifier = nodeIdentifier - 1 px = [abdomenStartX + abdomenLength, 0.0, 0.0] - + # options['Kinematic tree']['lumbar_body'] = px # arms - # rotate shoulder with arm, pivoting about shoulder drop below arm junction on network - # this has the realistic effect of shoulders becoming narrower with higher angles - # initial shoulder rotation with arm is negligible, hence: - shoulderRotationFactor = 1.0 - math.cos(0.5 * armAngleRadians) - # assume shoulder drop is half shrug distance to get limiting shoulder angle for 180 degree arm rotation - shoulderLimitAngleRadians = math.asin(1.5 * shoulderDrop / halfShoulderWidth) - shoulderAngleRadians = shoulderRotationFactor * shoulderLimitAngleRadians - armStartX = thoraxStartX + shoulderDrop - halfShoulderWidth * math.sin(shoulderAngleRadians) - nonHandArmLength = armLength - handLength - armScale = nonHandArmLength / (armToHandElementsCount - 2) # 2 == shoulder elements count - d12_mag = (halfWristThickness - armTopRadius) / (armToHandElementsCount - 2) - d13_mag = (halfWristWidth - armTopRadius) / (armToHandElementsCount - 2) for side in (left, right): - armAngle = armAngleRadians if (side == left) else -armAngleRadians + # Shoulder rotation + # rotate shoulder with arm, pivoting about shoulder drop below arm junction on network + # this has the realistic effect of shoulders becoming narrower with higher angles + # initial shoulder rotation with arm is negligible, hence: + armAbductionRadians = armLeftAngleRadians if (side == left) else armRightAngleRadians + shoulderRotationFactor = 1.0 - math.cos(0.5 * armAbductionRadians) + # assume shoulder drop is half shrug distance to get limiting shoulder angle for 180 degree arm rotation + shoulderLimitAngleRadians = math.asin(1.5 * shoulderDrop / halfShoulderWidth) + shoulderAngleRadians = shoulderRotationFactor * shoulderLimitAngleRadians + nonHandArmLength = armLength - handLength + # shoulderElementsCount = 2 + armScale = nonHandArmLength / (armToHandElementsCount - shoulderElementsCount) + d12_mag = (halfWristThickness - armTopRadius) / (armToHandElementsCount - shoulderElementsCount) + d13_mag = (halfWristWidth - armTopRadius) / (armToHandElementsCount - shoulderElementsCount) + armAngle = armAbductionRadians if (side == left) else -armAbductionRadians cosArmAngle = math.cos(armAngle) sinArmAngle = math.sin(armAngle) + armStartX = thoraxStartX + shoulderDrop - halfShoulderWidth * math.sin(shoulderAngleRadians) armStartY = (halfShoulderWidth if (side == left) else -halfShoulderWidth) * math.cos(shoulderAngleRadians) - x = [armStartX, armStartY, 0.0] + armStart = [armStartX, armStartY, 0.0] + x = armStart armDirn = [cosArmAngle, sinArmAngle, 0.0] armSide = [-sinArmAngle, cosArmAngle, 0.0] armFront = cross(armDirn, armSide) d1 = mult(armDirn, armScale) - # set leg versions 2 (left) and 3 (right) on leg junction node, and intermediate shoulder node + # set arm versions 2 (left) and 3 (right) on arm junction node, and intermediate shoulder node sd1 = interpolateLagrangeHermiteDerivative(sx, x, d1, 0.0) nx, nd1 = sampleCubicHermiteCurvesSmooth([sx, x], [sd1, d1], 2, derivativeMagnitudeEnd=armScale)[0:2] arcLengths = [getCubicHermiteArcLength(nx[i], nd1[i], nx[i + 1], nd1[i + 1]) for i in range(2)] sd2_list = [] sd3_list = [] sNodeIdentifiers = [] + side_label = 'l' if (side == left) else 'r' + options['Kinematic tree']['humerus_' + side_label] = nx[1] + # Upper shoulder nodes for i in range(2): sNodeIdentifiers.append(nodeIdentifier if (i > 0) else armJunctionNodeIdentifier) node = nodes.findNodeByIdentifier(sNodeIdentifiers[-1]) @@ -439,8 +600,70 @@ def generateBaseMesh(cls, region, options): nodeIdentifier += 1 setNodeFieldVersionDerivatives(coordinates, fieldcache, version, sd1, sd2, sd3) setNodeFieldVersionDerivatives(innerCoordinates, fieldcache, version, sd1, sid2, sid3) - sd2_list.append([-armTopRadius * sinArmAngle, armTopRadius * cosArmAngle, 0.0]) - sd3_list.append([0.0, 0.0, armTopRadius]) + # Arm twist + elementTwistAngle = ((armTwistAngleRadians if (side == left) else -armTwistAngleRadians) / + (armToHandElementsCount - 3)) + # Shoulder rotation + #Reset flexion angle w.r.t. current node + shoulderFlexionRadians = shoulderLeftFlexionRadians if (side == left) else shoulderRightFlexionRadians + ventralFlexion = True + if shoulderFlexionRadians < 0: + ventralFlexion = False + shoulderFlexionRadians = -shoulderFlexionRadians + #Reset abdution angle w.r.t. current node + shoulderAbductionRadians = (1 if (side == left) else -1)*angle([1,0,0], sd1) - armAngle + upwardAbduction = True + if shoulderAbductionRadians < 0: + upwardAbduction = False + shoulderAbductionRadians = -shoulderAbductionRadians + # Calculate magntiudes for joint node + i = 0 + xi = i / (armToHandElementsCount - 2) + halfWidth = xi * halfWristWidth + (1.0 - xi) * armTopRadius + halfThickness = xi * halfWristThickness + (1.0 - xi) * armTopRadius + upperShoulderDir = [nx[1], sd1, sd2, sd3, nd1[1], 0, 0] + armRotationRadians = armLeftRotationRadians if (side == left) else armRightRotationRadians + x1 = [0,1,0] + for i in range(1, 4): + upperShoulderDir[i] = set_magnitude(upperShoulderDir[i], 1) #normalize frame + shoulderDir, armDir = getJointRotationFrames( \ + shoulderAbductionRadians, shoulderFlexionRadians, upperShoulderDir, upwardAbduction, ventralFlexion) + shoulderDir[0] = armStart + shoulderDir[4:] = [armScale, halfThickness, halfWidth] + rotationCoeff = 0.2 + x, shoulderd2_mag, shoulderd3_mag, shoulderd12_mag, shoulderd13_mag = getJointRotationPosition( + shoulderAbductionRadians, shoulderFlexionRadians, upperShoulderDir, \ + shoulderDir, armDir, rotationCoeff, upwardAbduction, ventralFlexion + ) + armRotationRadians = armLeftRotationRadians if (side == left) else armRightRotationRadians + jointAbductionMatrix = axis_angle_to_rotation_matrix(mult(shoulderDir[1], 1), armRotationRadians) + for i in range(2, 4): + shoulderDir[i] = matrix_vector_mult(jointAbductionMatrix, shoulderDir[i]) + armDir[i] = matrix_vector_mult(jointAbductionMatrix, armDir[i]) + shoulderDirn, shoulderSide, shoulderFront = shoulderDir[1:4] + d1 = mult(shoulderDirn, armScale) + d2 = mult(shoulderSide, shoulderd2_mag) + d3 = mult(shoulderFront, shoulderd3_mag) + d12 = add( + mult(shoulderSide, d12_mag), + mult(shoulderDirn, shoulderd12_mag) + ) + d13 = add( + mult(shoulderFront, d13_mag), + mult(shoulderDirn, shoulderd13_mag) + ) + id2 = mult(d2, innerProportionDefault) + id3 = mult(d3, innerProportionDefault) + id12 = mult(d12, innerProportionDefault) + id13 = mult(d13, innerProportionDefault) + node = nodes.findNodeByIdentifier(nodeIdentifier) + fieldcache.setNode(node) + setNodeFieldParameters(coordinates, fieldcache, x, d1, d2, d3, d12, d13) + setNodeFieldParameters(innerCoordinates, fieldcache, x, d1, id2, id3, id12, id13) + nodeIdentifier += 1 + sd2_list.append(d2) + sd3_list.append(d3) + # Adjusting d12 and 13 for the upper Shoulder node for i in range(2): node = nodes.findNodeByIdentifier(sNodeIdentifiers[i]) fieldcache.setNode(node) @@ -453,20 +676,25 @@ def generateBaseMesh(cls, region, options): sid13 = mult(sd13, innerProportionDefault) innerCoordinates.setNodeParameters(fieldcache, -1, Node.VALUE_LABEL_D2_DS1DS2, version, sid12) innerCoordinates.setNodeParameters(fieldcache, -1, Node.VALUE_LABEL_D2_DS1DS3, version, sid13) - # main part of arm to wrist - elementTwistAngle = ((armTwistAngleRadians if (side == left) else -armTwistAngleRadians) / - (armToHandElementsCount - 3)) - for i in range(armToHandElementsCount - 1): + options['Kinematic tree']['ulna_' + side_label] = x + # Initial position for arm node + armDirn, armSide, armFront = armDir[1:4] + d1 = set_magnitude(armDirn, armScale) + armStart = add(x, d1) + armDir[0] = armStart + armStart = getDistalNodePosition(shoulderAbductionRadians, shoulderFlexionRadians, + upperShoulderDir, shoulderDir, armDir, 0.2) + # Setting brachium coordinates + j = 0 + for i in range(1, brachiumElementsCount): xi = i / (armToHandElementsCount - 2) node = nodes.findNodeByIdentifier(nodeIdentifier) fieldcache.setNode(node) - x = [armStartX + d1[0] * i, armStartY + d1[1] * i, d1[2] * i] + x = add(armStart, mult(d1, j)) halfThickness = xi * halfWristThickness + (1.0 - xi) * armTopRadius halfWidth = xi * halfWristWidth + (1.0 - xi) * armTopRadius if i == 0: twistAngle = 0.0 - elif i == (armToHandElementsCount - 2): - twistAngle = armTwistAngleRadians if (side == left) else -armTwistAngleRadians else: twistAngle = -0.5 * elementTwistAngle + elementTwistAngle * i if twistAngle == 0.0: @@ -483,7 +711,104 @@ def generateBaseMesh(cls, region, options): mult(armSide, halfWidth * sinTwistAngle)) d12 = set_magnitude(d2, d12_mag) d13 = set_magnitude(d3, d13_mag) - if i < (armToHandElementsCount - 2): + id2 = mult(d2, innerProportionDefault) + id3 = mult(d3, innerProportionDefault) + id12 = mult(d12, innerProportionDefault) + id13 = mult(d13, innerProportionDefault) + setNodeFieldParameters(coordinates, fieldcache, x, d1, d2, d3, d12, d13) + setNodeFieldParameters(innerCoordinates, fieldcache, x, d1, id2, id3, id12, id13) + nodeIdentifier += 1 + j += 1 + # Elbow + i += 1 + twistAngle = -0.5 * elementTwistAngle + elementTwistAngle * (i) + if twistAngle == 0.0: + d2 = armSide + d3 = armFront + else: + cosTwistAngle = math.cos(twistAngle) + sinTwistAngle = math.sin(twistAngle) + d2 = sub(mult(armSide, cosTwistAngle), + mult(armFront, sinTwistAngle)) + d3 = add(mult(armFront, cosTwistAngle), + mult(armSide, sinTwistAngle)) + armSide = d2 + armFront = d3 + # Updating frame of reference wrt flexion angle (using d2 as rotation axis) + elbowFlexionRadians = elbowLeftFlexionRadians if (side == left) else elbowRightFlexionRadians + elbowAbductionRadians = 0 + ventralFlexion = True + if shoulderFlexionRadians < 0: + ventralFlexion = False + shoulderFlexionRadians = -shoulderFlexionRadians + armDir[0] = x + xi = i / (armToHandElementsCount - 2) + halfThickness = xi * halfWristThickness + (1.0 - xi) * armTopRadius + halfWidth = xi * halfWristWidth + (1.0 - xi) * armTopRadius + rotationCoeff = 0.25 + elbowDir, antebrachiumDir = getJointRotationFrames( \ + 0, elbowFlexionRadians, armDir, upwardAbduction, ventralFlexion) + elbowDir[0] = add(x, d1) + elbowDir[4:] = [armScale, halfThickness, halfWidth] + x, elbowd2_mag, elbowd3_mag, elbowd12_mag, elbowd13_mag = getJointRotationPosition( + 0, elbowFlexionRadians, armDir, elbowDir, antebrachiumDir, rotationCoeff, upwardAbduction, ventralFlexion) + elbowDirn, elbowSide, elbowFront = elbowDir[1:4] + d1 = mult(elbowDirn, armScale) + d2 = mult(elbowSide, elbowd2_mag) + d3 = mult(elbowFront, elbowd3_mag) + d12 = add( + mult(elbowSide, d12_mag), + mult(elbowDirn, elbowd12_mag) + ) + d13 = add( + mult(elbowFront, d13_mag), + mult(elbowDirn, elbowd13_mag) + ) + id2 = mult(d2, innerProportionDefault) + id3 = mult(d3, innerProportionDefault) + id12 = mult(d12, innerProportionDefault) + id13 = mult(d13, innerProportionDefault) + node = nodes.findNodeByIdentifier(nodeIdentifier) + fieldcache.setNode(node) + setNodeFieldParameters(coordinates, fieldcache, x, d1, d2, d3, d12, d13) + setNodeFieldParameters(innerCoordinates, fieldcache, x, d1, id2, id3, id12, id13) + nodeIdentifier += 1 + options['Kinematic tree']['ulna_' + side_label] = x + # Antebrachium nodes starts after the elbow node + # antebrachiumStart = jointPositions[-1] + antebrachiumDirn, antebrachiumSide, antebrachiumFront = antebrachiumDir[1:4] + d1 = set_magnitude(antebrachiumDirn, armScale) + antebrachiumStart = add(x, d1) + antebrachiumDir[0] = antebrachiumStart + antebrachiumStart = getDistalNodePosition(elbowAbductionRadians, elbowFlexionRadians, + armDir, elbowDir, antebrachiumDir, rotationCoeff) + j=0 + for i in range(brachiumElementsCount + 1, armToHandElementsCount - 2): + xi = (i) / (armToHandElementsCount - 2) + node = nodes.findNodeByIdentifier(nodeIdentifier) + fieldcache.setNode(node) + x = add(antebrachiumStart, mult(d1, j)) + halfThickness = xi * halfWristThickness + (1.0 - xi) * armTopRadius + halfWidth = xi * halfWristWidth + (1.0 - xi) * armTopRadius + if i == 0: + twistAngle = armTwistAngleRadians if (side == left) else -armTwistAngleRadians + else: + twistAngle = -0.5 * elementTwistAngle + elementTwistAngle * (i) + if twistAngle == 0.0: + d2 = mult(antebrachiumSide, halfThickness) + d3 = mult(antebrachiumFront, halfWidth) + d12 = mult(antebrachiumSide, d12_mag) + d13 = mult(antebrachiumFront, d13_mag) + else: + cosTwistAngle = math.cos(twistAngle) + sinTwistAngle = math.sin(twistAngle) + d2 = sub(mult(antebrachiumSide, halfThickness * cosTwistAngle), + mult(antebrachiumFront, halfThickness * sinTwistAngle)) + d3 = add(mult(antebrachiumFront, halfWidth * cosTwistAngle), + mult(antebrachiumSide, halfWidth * sinTwistAngle)) + d12 = set_magnitude(d2, d12_mag) + d13 = set_magnitude(d3, d13_mag) + if i < (antebrachiumElementsCount - 1): d12 = add(d12, set_magnitude(d3, -halfThickness * elementTwistAngle)) d13 = add(d13, set_magnitude(d2, halfWidth * elementTwistAngle)) id2 = mult(d2, innerProportionDefault) @@ -493,40 +818,104 @@ def generateBaseMesh(cls, region, options): setNodeFieldParameters(coordinates, fieldcache, x, d1, d2, d3, d12, d13) setNodeFieldParameters(innerCoordinates, fieldcache, x, d1, id2, id3, id12, id13) nodeIdentifier += 1 - # hand + j += 1 + # Wrist flexion + i += 1 + twistAngle = -0.5 * elementTwistAngle + elementTwistAngle * (i) + if twistAngle == 0.0: + d2 = antebrachiumSide + d3 = antebrachiumFront + else: + cosTwistAngle = math.cos(twistAngle) + sinTwistAngle = math.sin(twistAngle) + d2 = sub(mult(antebrachiumSide, cosTwistAngle), + mult(antebrachiumFront, sinTwistAngle)) + d3 = add(mult(antebrachiumFront, cosTwistAngle), + mult(antebrachiumSide, sinTwistAngle)) + antebrachiumSide = d2 + antebrachiumFront = d3 + wristFlexionRadians = wristLeftFlexionRadians if (side == left) else wristRightFlexionRadians + wristAbductionRadians = wristLeftAbductionRadians if (side == left) else -wristRightAbductionRadians + ventralFlexion = True + if wristFlexionRadians < 0: + ventralFlexion = False + wristFlexionRadians = -wristFlexionRadians + upwardAbduction = True + if wristAbductionRadians < 0: + upwardAbduction = False + wristAbductionRadians = -wristAbductionRadians + antebrachiumDir[0] = x + xi = i / (armToHandElementsCount - 2) + halfThickness = halfWristThickness + halfWidth = halfWristWidth + rotationCoeff = 0.3 + wristDir, handDir = getJointRotationFrames( + wristAbductionRadians, wristFlexionRadians, antebrachiumDir, upwardAbduction, ventralFlexion + ) + wristDir[0] = add(x, d1) + wristDir[4:] = [armScale, halfThickness, halfWidth] + x, wristd2_mag, wristd3_mag, wristd12_mag, wristd13_mag = getJointRotationPosition( + wristAbductionRadians, wristFlexionRadians, antebrachiumDir, wristDir, handDir, rotationCoeff, upwardAbduction, ventralFlexion + ) + wristDirn, wristSide, wristFront = wristDir[1:4] + handDirn, handSide, handFront = handDir[1:4] + d1 = mult(handDirn, armScale) + d2 = mult(wristSide, wristd2_mag) + d3 = mult(wristFront, wristd3_mag) + d12 = add( + mult(wristSide, d12_mag), + mult(handDirn, wristd12_mag) + ) + d13 = add( + mult(wristFront, d13_mag), + mult(handDirn, wristd13_mag) + ) + id2 = mult(d2, innerProportionDefault) + id3 = mult(d3, innerProportionDefault) + id12 = mult(d12, innerProportionDefault) + id13 = mult(d13, innerProportionDefault) + node = nodes.findNodeByIdentifier(nodeIdentifier) + fieldcache.setNode(node) + setNodeFieldParameters(coordinates, fieldcache, x, d1, d2, d3, d12, d13) + setNodeFieldParameters(innerCoordinates, fieldcache, x, d1, id2, id3, id12, id13) + nodeIdentifier += 1 + options['Kinematic tree']['hand_' + side_label] = x + + d1 = set_magnitude(handDirn, armScale) + handStart = add(x, mult(handDirn, handLength)) + handDir[0] = handStart + # Hand assert handElementsCount == 1 node = nodes.findNodeByIdentifier(nodeIdentifier) fieldcache.setNode(node) - hx = [armStartX + armLength * cosArmAngle, armStartY + armLength * sinArmAngle, 0.0] - hd1 = computeCubicHermiteEndDerivative(x, d1, hx, d1) + hd1 = computeCubicHermiteEndDerivative(x, d1, handStart, d1) twistAngle = armTwistAngleRadians if (side == left) else -armTwistAngleRadians - if twistAngle == 0.0: - hd2 = set_magnitude(d2, halfHandThickness) - hd3 = [0.0, 0.0, halfHandWidth] + if twistAngle >= 0.0: + hd2 = set_magnitude(handSide, halfHandThickness) + hd3 = set_magnitude(handFront, halfHandWidth) else: cosTwistAngle = math.cos(twistAngle) sinTwistAngle = math.sin(twistAngle) - hd2 = sub(mult(armSide, halfHandThickness * cosTwistAngle), - mult(armFront, halfHandThickness * sinTwistAngle)) - hd3 = add(mult(armFront, halfHandWidth * cosTwistAngle), - mult(armSide, halfHandWidth * sinTwistAngle)) + hd2 = sub(mult(handSide, halfHandThickness * cosTwistAngle), + mult(handFront, halfHandThickness * sinTwistAngle)) + hd3 = add(mult(handFront, halfHandWidth * cosTwistAngle), + mult(handSide, halfHandWidth * sinTwistAngle)) hid2 = mult(hd2, innerProportionDefault) hid3 = mult(hd3, innerProportionDefault) - setNodeFieldParameters(coordinates, fieldcache, hx, hd1, hd2, hd3) - setNodeFieldParameters(innerCoordinates, fieldcache, hx, hd1, hid2, hid3) + setNodeFieldParameters(coordinates, fieldcache, handStart, hd1, hd2, hd3) + setNodeFieldParameters(innerCoordinates, fieldcache, handStart, hd1, hid2, hid3) nodeIdentifier += 1 - # legs legStartX = abdomenStartX + abdomenLength + pelvisDrop nonFootLegLength = legLength - footHeight - legScale = nonFootLegLength / (legToFootElementsCount - 1) - d12_mag = (legBottomRadius - legTopRadius) / (armToHandElementsCount - 2) - d13_mag = (legBottomRadius - legTopRadius) / (armToHandElementsCount - 2) - + legScale = nonFootLegLength / (legToFootElementsCount - 1) + d12_mag = (legBottomRadius - legTopRadius) / (legToFootElementsCount) + d13_mag = (legBottomRadius - legTopRadius) / (legToFootElementsCount) pd3 = [0.0, 0.0, 0.5 * legTopRadius + 0.5 * halfTorsoDepth] pid3 = mult(pd3, innerProportionDefault) for side in (left, right): - legAngle = legAngleRadians if (side == left) else -legAngleRadians + side_label = 'l' if (side == left) else 'r' + legAngle = leftLegAbductionRadians if (side == left) else -rightLegAbductionRadians cosLegAngle = math.cos(legAngle) sinLegAngle = math.sin(legAngle) legStartY = halfPelvisWidth if (side == left) else -halfPelvisWidth @@ -552,47 +941,223 @@ def generateBaseMesh(cls, region, options): id12 = mult(d12, innerProportionDefault) d13 = [0.0, 0.0, d13_mag] id13 = mult(d13, innerProportionDefault) - # main part of leg to ankle - for i in range(legToFootElementsCount): + # options['Kinematic tree']['femur_' + side_label] = x + # Upper leg + for i in range(hipElementsCount-1): xi = i / legToFootElementsCount node = nodes.findNodeByIdentifier(nodeIdentifier) fieldcache.setNode(node) - x = [legStartX + d1[0] * i, legStartY + d1[1] * i, d1[2] * i] + x = add(legStart, mult(d1, i)) radius = xi * legBottomRadius + (1.0 - xi) * legTopRadius d2 = mult(legSide, radius) - d3 = [0.0, 0.0, radius] + d3 = mult(legFront, radius) + d13 = mult(legFront, d13_mag) + d12 = mult(legSide, d12_mag) + id2 = mult(d2, innerProportionDefault) + id3 = mult(d3, innerProportionDefault) + id13 = mult(d13, innerProportionDefault) + id12 = mult(d12, innerProportionDefault) + setNodeFieldParameters(coordinates, fieldcache, x, d1, d2, d3, d12, d13) + setNodeFieldParameters(innerCoordinates, fieldcache, x, d1, id2, id3, id12, id13) + nodeIdentifier += 1 + # Frontal hip flexion + # Updating frame of reference wrt flexion angle (using d2 as rotation axis) + hipFlexionRadians = hipLeftFlexionRadians if (side == left) else hipRightFlexionRadians + ventralFlexion = True + if hipFlexionRadians < 0: + ventralFlexion = False + hipFlexionRadians = -hipFlexionRadians + i += 1 + xi = i / legToFootElementsCount + radius = xi * legBottomRadius + (1.0 - xi) * legTopRadius + legDir = [x, legDirn, legSide, legFront, legScale, 0, 0] + rotationCoeff = 0.18 + hipDir, upperLegDir = getJointRotationFrames( \ + 0, hipFlexionRadians, legDir, True, ventralFlexion) + hipDir[0] = add(x, d1) + hipDir[4:] = [legScale, radius, radius] + x, hipd2_mag, hipd3_mag, hipd12_mag, hipd13_mag = getJointRotationPosition( + 0, hipFlexionRadians, legDir, hipDir, upperLegDir, rotationCoeff + ) + hipDirn, hipSide, hipFront = hipDir[1:4] + d1 = mult(hipDirn, legScale) + d2 = mult(hipSide, hipd2_mag) + d3 = mult(hipFront, 0.9*hipd3_mag) + d12 = mult(hipSide, d12_mag) + d13 = add( + mult(hipFront, d13_mag), + mult(hipDirn, 0.7*hipd13_mag) + ) + id2 = mult(d2, innerProportionDefault) + id3 = mult(d3, innerProportionDefault) + id12 = mult(d12, innerProportionDefault) + id13 = mult(d13, innerProportionDefault) + node = nodes.findNodeByIdentifier(nodeIdentifier) + fieldcache.setNode(node) + setNodeFieldParameters(coordinates, fieldcache, x, d1, d2, d3, d12, d13) + setNodeFieldParameters(innerCoordinates, fieldcache, x, d1, id2, id3, id12, id13) + nodeIdentifier += 1 + options['Kinematic tree']['femur_' + side_label] = x + upperLegDirn, upperLegSide, upperLegFront = upperLegDir[1:4] + d1 = set_magnitude(upperLegDirn, legScale) + upperLegStart = add(x, d1) + upperLegDir[0] = upperLegStart + upperLegStart = getDistalNodePosition(0, hipFlexionRadians, + legDir, hipDir, upperLegDir, rotationCoeff) + # Rest of upper leg + j = 0 + for i in range(hipElementsCount, hipElementsCount+upperLegElementsCount-1): + xi = i / legToFootElementsCount + node = nodes.findNodeByIdentifier(nodeIdentifier) + fieldcache.setNode(node) + x = add(upperLegStart, mult(d1, j)) + radius = xi * legBottomRadius + (1.0 - xi) * legTopRadius + d2 = mult(upperLegSide, radius) + d3 = mult(upperLegFront, radius) + d12 = mult(upperLegSide, d12_mag) + d13 = mult(upperLegFront, d13_mag) + id2 = mult(d2, innerProportionDefault) + id3 = mult(d3, innerProportionDefault) + id12 = set_magnitude(d12, d12_mag) + id13 = set_magnitude(d13, d13_mag) + setNodeFieldParameters(coordinates, fieldcache, x, d1, d2, d3, d12, d13) + setNodeFieldParameters(innerCoordinates, fieldcache, x, d1, id2, id3, id12, id13) + nodeIdentifier += 1 + j += 1 + # knee + kneeFlexionRadians = kneeLeftFlexionRadians if (side == left) else kneeRightFlexionRadians + ventralFlexion = False + if shoulderFlexionRadians < 0: + ventralFlexion = True + kneeFlexionRadians = -kneeFlexionRadians + # # Set coordiantes for joint node + i += 1 + xi = i / legToFootElementsCount + radius = xi * legBottomRadius + (1.0 - xi) * legTopRadius + rotationCoeff = 0.25 + kneeDir, lowerLegDir = getJointRotationFrames( + 0, kneeFlexionRadians, upperLegDir, True, ventralFlexion + ) + kneeDir[0] = add(x, d1) + kneeDir[4:] = [legScale, radius, radius] + x, kneed2_mag, kneed3_mag, kneed12_mag, kneed13_mag = getJointRotationPosition( + 0, kneeFlexionRadians, upperLegDir, kneeDir, lowerLegDir, rotationCoeff, True, ventralFlexion) + kneeDirn, kneeSide, kneeFront = kneeDir[1:4] + d1 = mult(kneeDirn, legScale) + d2 = mult(kneeSide, kneed2_mag) + d3 = mult(kneeFront, kneed3_mag) + d12 = mult(kneeSide, d12_mag) + d13 = add( + set_magnitude(d3, d13_mag), + set_magnitude(d1, kneed13_mag) + ) + id2 = mult(d2, innerProportionDefault) + id3 = mult(d3, innerProportionDefault) + id12 = mult(d12, innerProportionDefault) + id13 = mult(d13, innerProportionDefault) + node = nodes.findNodeByIdentifier(nodeIdentifier) + fieldcache.setNode(node) + setNodeFieldParameters(coordinates, fieldcache, x, d1, d2, d3, d12, d13) + setNodeFieldParameters(innerCoordinates, fieldcache, x, d1, id2, id3, id12, id13) + nodeIdentifier += 1 + options['Kinematic tree']['tibia_' + side_label] = x + # Lower leg + lowerLegDirn, lowerLegSide, lowerLegFront = lowerLegDir[1:4] + d1 = set_magnitude(lowerLegDirn, legScale) + lowerLegStart = add(x, d1) + lowerLegDir[0] = lowerLegStart + lowerLegStart = getDistalNodePosition(0, kneeFlexionRadians, + upperLegDir, kneeDir, lowerLegDir, rotationCoeff) + j = 0 + for i in range(hipElementsCount+upperLegElementsCount, legToFootElementsCount-1): + xi = i / legToFootElementsCount + node = nodes.findNodeByIdentifier(nodeIdentifier) + fieldcache.setNode(node) + x = add(lowerLegStart, mult(d1, j)) + radius = xi * legBottomRadius + (1.0 - xi) * legTopRadius + d2 = set_magnitude(lowerLegSide, radius) + d3 = set_magnitude(lowerLegFront, radius) + d12 = set_magnitude(d2, d12_mag) + d13 = set_magnitude(d3, d13_mag) id2 = mult(d2, innerProportionDefault) id3 = mult(d3, innerProportionDefault) + id12 = set_magnitude(d12, d12_mag) + id13 = set_magnitude(d13, d13_mag) setNodeFieldParameters(coordinates, fieldcache, x, d1, d2, d3, d12, d13) setNodeFieldParameters(innerCoordinates, fieldcache, x, d1, id2, id3, id12, id13) nodeIdentifier += 1 + j+=1 # foot - fx = [x, - add(add(legStart, mult(legDirn, legLength - 1.5 * halfFootThickness)), - [0.0, 0.0, legBottomRadius]), - add(add(legStart, mult(legDirn, legLength - halfFootThickness)), - [0.0, 0.0, footLength - legBottomRadius])] - fd1 = smoothCubicHermiteDerivativesLine( - fx, [d1, [0.0, 0.0, 0.5 * footLength], [0.0, 0.0, 0.5 * footLength]], - fixAllDirections=True, fixStartDerivative=True) - fd2 = [d2, mult(legSide, halfFootWidth), mult(legSide, halfFootWidth)] - fd3 = [d3, - set_magnitude(sub(legFront, legDirn), - math.sqrt(2.0 * halfFootThickness * halfFootThickness) + legBottomRadius), - set_magnitude(cross(fd1[2], fd2[2]), halfFootThickness)] - fd12 = sub(fd2[2], fd2[1]) - fd13 = sub(fd3[2], fd3[1]) - fid12 = mult(fd12, innerProportionDefault) - fid13 = mult(fd13, innerProportionDefault) - for i in range(1, 3): + ankleFlexionRadians = ankleLeftFlexionRadians if (side == left) else ankleRightFlexionRadians + ventralFlexion = True + if ankleFlexionRadians < 0: + ventralFlexion = False + ankleFlexionRadians = -ankleFlexionRadians + i += 1 + radius = math.sqrt(2)*legBottomRadius + rotationCoeff = 0.4 + ankleDir, footDir = getJointRotationFrames( + 0, ankleFlexionRadians, lowerLegDir, True, ventralFlexion + ) + ankleDir[0] = add(x, d1) + ankleDir[4:] = [legScale, radius, radius] + [x, ankled2_mag, ankled3_mag, ankled12_mag, ankled13_mag] = getJointRotationPosition( + 0, ankleFlexionRadians, lowerLegDir, ankleDir, footDir, rotationCoeff, True, ventralFlexion + ) + ankleDirn, ankleSide, ankleFront = ankleDir[1:4] + footDirn, footSide, footFront = footDir[1:4] + d1 = mult(footDirn, legScale) + d2 = set_magnitude(ankleSide, ankled2_mag) + d3 = set_magnitude(ankleFront, ankled3_mag) + d12 = set_magnitude(ankleSide, d12_mag) + d13 = add( + set_magnitude(ankleFront, d13_mag), + set_magnitude(footDirn, ankled13_mag) + ) + id2 = mult(d2, innerProportionDefault) + id3 = mult(d3, innerProportionDefault) + id12 = mult(d12, innerProportionDefault) + id13 = mult(d13, innerProportionDefault) + node = nodes.findNodeByIdentifier(nodeIdentifier) + fieldcache.setNode(node) + setNodeFieldParameters(coordinates, fieldcache, x, d1, d2, d3, d12, d13) + setNodeFieldParameters(innerCoordinates, fieldcache, x, d1, id2, id3, id12, id13) + nodeIdentifier += 1 + # Foot end nodes + d1 = mult(footDirn, footLength) + footStart = add(x, d1) + footDir[0] = footStart + j = 0 + for i in range(footElementsCount): node = nodes.findNodeByIdentifier(nodeIdentifier) fieldcache.setNode(node) - setNodeFieldParameters(coordinates, fieldcache, fx[i], fd1[i], fd2[i], fd3[i], fd12, fd13) - fid2 = mult(fd2[i], innerProportionDefault) - fid3 = mult(fd3[i], innerProportionDefault) - setNodeFieldParameters(innerCoordinates, fieldcache, fx[i], fd1[i], fid2, fid3, fid12, fid13) + x = add(footStart, mult(d1, j)) + d2 = set_magnitude(footSide, halfFootWidth) + d3 = set_magnitude(footFront, halfFootThickness) + d12 = set_magnitude(d2, d12_mag) + d13 = sub(d3, set_magnitude(ankleFront, ankled3_mag)) + id2 = mult(d2, innerProportionDefault) + id3 = mult(d3, innerProportionDefault) + id12 = set_magnitude(d12, d12_mag) + id13 = set_magnitude(d13, d13_mag) + setNodeFieldParameters(coordinates, fieldcache, x, d1, d2, d3, d12, d13) + setNodeFieldParameters(innerCoordinates, fieldcache, x, d1, id2, id3, id12, id13) nodeIdentifier += 1 + j+=1 + options['Kinematic tree']['toes_' + side_label] = x + # Kinematic tree markers + nodes = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) + node_identifier = max(1, get_maximum_node_identifier(nodes) + 1) + coordinates = find_or_create_field_coordinates(fieldmodule) + stickman_markers = options['Kinematic tree'] + for marker_name, marker_position in stickman_markers.items(): + marker_group = findOrCreateAnnotationGroupForTerm( + annotationGroups, region, (marker_name, ""), isMarker=True + ) + marker_group.createMarkerNode( + node_identifier, coordinates, marker_position + ) return annotationGroups, networkMesh @classmethod @@ -605,8 +1170,63 @@ def getInteractiveFunctions(cls): if interactiveFunction[0] == "Edit structure...": interactiveFunctions.remove(interactiveFunction) break + interactiveFunctions = interactiveFunctions + [\ + ("Align to body markers", + {"Load transformation parameters into Settings": True}, + lambda region, options, networkMesh, functionOptions, editGroupName: + cls.alignNetworkLayoutToMarkers(region, options, networkMesh, functionOptions, editGroupName) + ), + ] return interactiveFunctions + @classmethod + def alignNetworkLayoutToMarkers(cls, region, options, networkMesh, functionOptions, editGroupName): + + # Load the data into the current region + data_region = region.getParent().findChildByName('data') + if not data_region.isValid(): + logger.warning('Missing input data') + return None + data_fieldmodule = data_region.getFieldmodule() + with ChangeManager(data_fieldmodule): + copy_fitting_data(region, data_region) + del data_region + # Setup fitting routine + fieldmodule = region.getFieldmodule() + fitter = GeometryFitter(region=region) + fitter.getInitialFitterStepConfig() + # Manually run the _loadModel routines + fitter._discoverModelCoordinatesField() + fitter._discoverModelFitGroup() + # + zero_fibres = find_or_create_field_zero_fibres(fieldmodule) + fitter.setFibreField(zero_fibres) + del zero_fibres + # + fitter._discoverFlattenGroup() + fitter.defineCommonMeshFields() + # Manually run the _loadData routines + fitter._discoverDataCoordinatesField() + fitter._discoverMarkerGroup() + # Stuff + fitter.defineCommonMeshFields() + fitter.defineDataProjectionFields() + fitter.initializeFit() + # Call Align step + fit1 = FitterStepAlign() + fitter.addFitterStep(fit1) + fit1.setAlignMarkers(True) + fit1._doAutoAlign() + # Pass the graphical transformation settings into options + options['rotation'] = fit1._rotation + options['scale'] = fit1._scale + options['translation'] = fit1._translation + del fit1 + return True, False + + + + class MeshType_3d_wholebody2(Scaffold_base): """ @@ -636,10 +1256,14 @@ def getDefaultOptions(cls, parameterSetName="Default"): options["Number of elements along neck"] = 1 options["Number of elements along thorax"] = 2 options["Number of elements along abdomen"] = 2 - options["Number of elements along arm to hand"] = 5 + options["Number of elements along shoulder"] = 2 + options["Number of elements along brachium"] = 3 + options["Number of elements along antebrachium"] = 2 options["Number of elements along hand"] = 1 - options["Number of elements along leg to foot"] = 4 - options["Number of elements along foot"] = 2 + options["Number of elements along hip"] = 2 + options["Number of elements along upper leg"] = 3 + options["Number of elements along lower leg"] = 3 + options["Number of elements along foot"] = 1 options["Number of elements around head"] = 12 options["Number of elements around torso"] = 12 options["Number of elements around arm"] = 8 @@ -654,9 +1278,12 @@ def getDefaultOptions(cls, parameterSetName="Default"): options["Number of elements along neck"] = 2 options["Number of elements along thorax"] = 3 options["Number of elements along abdomen"] = 3 - options["Number of elements along arm to hand"] = 6 + options["Number of elements along shoulder"] = 2 + options["Number of elements along brachium"] = 3 + options["Number of elements along antebrachium"] = 3 options["Number of elements along hand"] = 1 - options["Number of elements along leg to foot"] = 6 + options["Number of elements along upper leg"] = 2 + options["Number of elements along lower leg"] = 2 options["Number of elements along foot"] = 2 options["Number of elements around head"] = 16 options["Number of elements around torso"] = 16 @@ -666,9 +1293,12 @@ def getDefaultOptions(cls, parameterSetName="Default"): options["Number of elements along neck"] = 2 options["Number of elements along thorax"] = 4 options["Number of elements along abdomen"] = 4 - options["Number of elements along arm to hand"] = 8 + options["Number of elements along shoulder"] = 2 + options["Number of elements along brachium"] = 3 + options["Number of elements along antebrachium"] = 4 options["Number of elements along hand"] = 2 - options["Number of elements along leg to foot"] = 8 + options["Number of elements along upper leg"] = 3 + options["Number of elements along lower leg"] = 2 options["Number of elements along foot"] = 3 options["Number of elements around head"] = 20 options["Number of elements around torso"] = 20 @@ -687,9 +1317,13 @@ def getOrderedOptionNames(cls): "Number of elements along neck", "Number of elements along thorax", "Number of elements along abdomen", - "Number of elements along arm to hand", + "Number of elements along shoulder", + "Number of elements along brachium", + "Number of elements along antebrachium", "Number of elements along hand", - "Number of elements along leg to foot", + "Number of elements along hip", + "Number of elements along upper leg", + "Number of elements along lower leg", "Number of elements along foot", "Number of elements around head", "Number of elements around torso", @@ -735,9 +1369,11 @@ def checkOptions(cls, options): "Number of elements along neck", "Number of elements along thorax", "Number of elements along abdomen", - "Number of elements along arm to hand", + "Number of elements along brachium", + "Number of elements along antebrachium", "Number of elements along hand", - "Number of elements along leg to foot", + "Number of elements along upper leg", + "Number of elements along lower leg", "Number of elements along foot" ]: if options[key] < 1: @@ -790,9 +1426,13 @@ def generateBaseMesh(cls, region, options): elementsCountAlongNeck = options["Number of elements along neck"] elementsCountAlongThorax = options["Number of elements along thorax"] elementsCountAlongAbdomen = options["Number of elements along abdomen"] - elementsCountAlongArmToHand = options["Number of elements along arm to hand"] + elementsCountAlongShoulder = options["Number of elements along shoulder"] + elementsCountAlongBrachium = options["Number of elements along brachium"] + elementsCountAlongAntebrachium = options["Number of elements along antebrachium"] elementsCountAlongHand = options["Number of elements along hand"] - elementsCountAlongLegToFoot = options["Number of elements along leg to foot"] + elementsCountAlongHip = options["Number of elements along hip"] + elementsCountAlongUpperLeg = options["Number of elements along upper leg"] + elementsCountAlongLowerLeg = options["Number of elements along lower leg"] elementsCountAlongFoot = options["Number of elements along foot"] elementsCountAroundHead = options["Number of elements around head"] elementsCountAroundTorso = options["Number of elements around torso"] @@ -828,14 +1468,26 @@ def generateBaseMesh(cls, region, options): alongCount = elementsCountAlongAbdomen aroundCount = elementsCountAroundTorso coreBoundaryScalingMode = 2 - elif "arm to hand" in name: - alongCount = elementsCountAlongArmToHand + elif "shoulder" in name: + alongCount = elementsCountAlongShoulder + aroundCount = elementsCountAroundArm + elif " brachium" in name: + alongCount = elementsCountAlongBrachium + aroundCount = elementsCountAroundArm + elif " antebrachium" in name: + alongCount = elementsCountAlongAntebrachium aroundCount = elementsCountAroundArm elif "hand" in name: alongCount = elementsCountAlongHand aroundCount = elementsCountAroundArm - elif "leg to foot" in name: - alongCount = elementsCountAlongLegToFoot + elif "hip" in name: + alongCount = elementsCountAlongHip + aroundCount = elementsCountAroundLeg + elif "upper leg" in name: + alongCount = elementsCountAlongUpperLeg + aroundCount = elementsCountAroundLeg + elif "lower leg" in name: + alongCount = elementsCountAlongLowerLeg aroundCount = elementsCountAroundLeg elif "foot" in name: alongCount = elementsCountAlongFoot @@ -868,9 +1520,8 @@ def generateBaseMesh(cls, region, options): isShowTrimSurfaces=options["Show trim surfaces"]) tubeNetworkMeshBuilder.generateMesh(generateData) annotationGroups = generateData.getAnnotationGroups() - + fieldmodule = region.getFieldmodule() if isCore: - fieldmodule = region.getFieldmodule() mesh = fieldmodule.findMeshByDimension(meshDimension) thoraxGroup = getAnnotationGroupForTerm(annotationGroups, get_body_term("thorax")) abdomenGroup = getAnnotationGroupForTerm(annotationGroups, get_body_term("abdomen")) @@ -886,6 +1537,18 @@ def generateBaseMesh(cls, region, options): is_abdominal_cavity = fieldmodule.createFieldAnd(abdomenGroup.getGroup(), coreGroup.getGroup()) abdominalCavityGroup.getMeshGroup(mesh).addElementsConditional(is_abdominal_cavity) + # Kinematic tree markers + nodes = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) + node_identifier = max(1, get_maximum_node_identifier(nodes) + 1) + coordinates = find_or_create_field_coordinates(fieldmodule) + stickman_markers = networkLayout._scaffoldSettings['Kinematic tree'] + for marker_name, marker_position in stickman_markers.items(): + marker_group = findOrCreateAnnotationGroupForTerm( + annotationGroups, region, (marker_name, ""), isMarker=True + ) + marker_group.createMarkerNode( + node_identifier, coordinates, marker_position + ) return annotationGroups, None @classmethod @@ -893,6 +1556,7 @@ def defineFaceAnnotations(cls, region, options, annotationGroups): """ Add face annotation groups from the highest dimension mesh. Must have defined faces and added subelements for highest dimension groups. + :param region: Zinc region containing model. :param options: Dict containing options. See getDefaultOptions(). :param annotationGroups: List of annotation groups for top-level elements. @@ -928,7 +1592,7 @@ def defineFaceAnnotations(cls, region, options, annotationGroups): annotationGroups, region, get_body_term("left lower limb skin epidermis outer surface")) leftLegSkinGroup.getMeshGroup(mesh2d).addElementsConditional( fieldmodule.createFieldAnd(leftLegGroup.getGroup(), is_exterior)) - rightLegGroup = getAnnotationGroupForTerm(annotationGroups, get_body_term("right lower limb ")) + rightLegGroup = getAnnotationGroupForTerm(annotationGroups, get_body_term("right lower limb")) rightLegSkinGroup = findOrCreateAnnotationGroupForTerm( annotationGroups, region, get_body_term("right lower limb skin epidermis outer surface")) rightLegSkinGroup.getMeshGroup(mesh2d).addElementsConditional( @@ -981,6 +1645,7 @@ def defineFaceAnnotations(cls, region, options, annotationGroups): def setNodeFieldParameters(field, fieldcache, x, d1, d2, d3, d12=None, d13=None): """ Assign node field parameters x, d1, d2, d3 of field. + :param field: Field parameters to assign. :param fieldcache: Fieldcache with node set. :param x: Parameters to set for Node.VALUE_LABEL_VALUE. @@ -1004,6 +1669,7 @@ def setNodeFieldParameters(field, fieldcache, x, d1, d2, d3, d12=None, d13=None) def setNodeFieldVersionDerivatives(field, fieldcache, version, d1, d2, d3, d12=None, d13=None): """ Assign node field parameters d1, d2, d3 of field. + :param field: Field to assign parameters of. :param fieldcache: Fieldcache with node set. :param version: Version of d1, d2, d3 >= 1. @@ -1021,3 +1687,103 @@ def setNodeFieldVersionDerivatives(field, fieldcache, version, d1, d2, d3, d12=N field.setNodeParameters(fieldcache, -1, Node.VALUE_LABEL_D2_DS1DS2, version, d12) if d13: field.setNodeParameters(fieldcache, -1, Node.VALUE_LABEL_D2_DS1DS3, version, d13) + + +def getJointRotationFrames(jointAbductionRadians, jointFlexionRadians, proximalDir, upwardAbduction=True, ventralFlexion=True): + proximalDirn, proximalSide, proximalFront = proximalDir[1:4] + ventral = 1 if upwardAbduction else -1 + jointAbductionMatrix = axis_angle_to_rotation_matrix(mult(proximalFront, -1*ventral), jointAbductionRadians) + jointHalfAbductionMatrix = axis_angle_to_rotation_matrix(mult(proximalFront, -1*ventral), jointAbductionRadians/2) + ventral = 1 if ventralFlexion else -1 + jointFlexionMatrix = axis_angle_to_rotation_matrix(mult(proximalSide, -1*ventral), jointFlexionRadians) + jointHalfFlexionMatrix = axis_angle_to_rotation_matrix(mult(proximalSide, -1*ventral), jointFlexionRadians/2) + jointRotationMatrix = matrix_mult(jointFlexionMatrix,jointAbductionMatrix) + jointHalfRotationMatrix = matrix_mult(jointHalfFlexionMatrix, jointHalfAbductionMatrix) + # Joint directions (frame is rotated by half the abduction angle) + jointDirn = matrix_vector_mult(jointHalfRotationMatrix, proximalDirn) + jointSide = matrix_vector_mult(jointHalfRotationMatrix, proximalSide) + jointFront = matrix_vector_mult(jointHalfRotationMatrix, proximalFront) + # Distal directions + distalDirn = matrix_vector_mult(jointRotationMatrix, proximalDirn) + distalSide = matrix_vector_mult(jointRotationMatrix, proximalSide) + distalFront = matrix_vector_mult(jointRotationMatrix, proximalFront) + jointDir = [0, jointDirn, jointSide, jointFront, 0, 0, 0] + distalDir = [0, distalDirn, distalSide, distalFront, 0, 0, 0] + return jointDir, distalDir + + +def getJointRotationPosition(jointAbductionRadians, jointFlexionRadians, \ + proximalDir, jointDir, distalDir, rotationCoeff, upwardAbduction=True, ventralFlexion=True): + ventral = 1 if ventralFlexion else -1 + upward = 1 if upwardAbduction else -1 + proximalNodePosition, proximalDirn, proximalSide, proximalFront = proximalDir[0:4] + jointNodePosition, jointDirn, jointSide, jointFront = jointDir[0:4] + distalNodePosition, distalDirn, distalSide, distalFront = distalDir[0:4] + # Flexion node adjustment + jointAngleRadians = math.pi - jointFlexionRadians + frontScale = jointDir[6] + rotationRotFactor = 1*math.sin(jointAngleRadians) + jointRotFactor = 1/math.sin(jointAngleRadians/2) + d13RotFactor = math.sqrt(2)*math.tan(jointFlexionRadians/2) + rotDisplacementFactor = rotationCoeff*frontScale*rotationRotFactor + jointNodePosition = adjustJointNodePosition(proximalNodePosition, jointNodePosition, \ + jointFront, distalDirn, rotDisplacementFactor, ventralFlexion) + d3_mag = frontScale*(jointRotFactor-(rotationCoeff*rotationRotFactor)) + d13_mag = -1*ventral*frontScale*d13RotFactor + # Abduction node adjustment + jointAngleRadians = math.pi - jointAbductionRadians + sideScale = jointDir[5] + jointSide = mult(jointSide, -1) + rotationRotFactor = 1*math.sin(jointAngleRadians) + jointRotFactor = 1/math.sin(jointAngleRadians/2) + d12RotFactor = math.sqrt(2)*math.tan(jointAbductionRadians/2) + rotDisplacementFactor = rotationCoeff*sideScale*rotationRotFactor + jointNodePosition = adjustJointNodePosition(proximalNodePosition, jointNodePosition, \ + jointSide, distalDirn, rotDisplacementFactor, upwardAbduction) + d2_mag = sideScale*(jointRotFactor-(rotationCoeff*rotationRotFactor)) + d12_mag = upward*sideScale*d12RotFactor + return [jointNodePosition, d2_mag, d3_mag, d12_mag, d13_mag] + +def getDistalNodePosition(jointAbductionRadians, jointFlexionRadians, proximalDir, jointDir, distalDir, rotationCoeff): + proximalNodePosition, proximalDirn, proximalSide, proximalFront = proximalDir[0:4] + jointNodePosition, jointDirn, jointSide, jointFront = jointDir[0:4] + distalNodePosition, distalDirn, distalSide, distalFront = distalDir[0:4] + # Flexion + jointFrontScale = jointDir[6] + jointAngleRadians = math.pi - jointFlexionRadians + flexionRotFactor = 1*math.sin(jointAngleRadians) + rotDisplacementFactor = rotationCoeff*jointFrontScale*flexionRotFactor + distalNodePosition = adjustJointNodePosition(jointNodePosition, distalNodePosition, \ + jointDirn, proximalDirn, rotDisplacementFactor) + # Abduction + jointSideScale = jointDir[5] + jointAngleRadians = math.pi - jointAbductionRadians + abductionRotFactor = 1*math.sin(jointAngleRadians) + rotDisplacementFactor = rotationCoeff*jointSideScale*abductionRotFactor + distalNodePosition = adjustJointNodePosition(jointNodePosition, distalNodePosition, \ + jointDirn, proximalDirn, rotDisplacementFactor) + return distalNodePosition + +def adjustJointNodePosition(proximalNodePosition, nodePosition, d1, d2, dispFactor, ventral = True): + """ + Displaces a node towards the center of the tube network, + while preserving the distance between the current node and the previous (proximal) node + + :param proximalNodePosition: Position of the node previous to nodePosition + :param nodePosition: Original position of the node to be rotated. + :param lenScale: Distance between the joint node and the proximal/distal node. + :param d1: Typically jointFront or or jointSide. + :param d2: Typically distalDirn or jointDirn. + :param dispFactor: Measure of displacement of the node away from the corner. + """ + initialDir = sub(nodePosition, proximalNodePosition) + lenScale = magnitude(initialDir) + ventral = 1 if (ventral) else -1 + jointAdjustDir = add( + set_magnitude(d1, ventral*dispFactor), + set_magnitude(d2, dispFactor), + ) + jointAdjustDir = add(initialDir, jointAdjustDir) + jointAdjustDir = set_magnitude(jointAdjustDir, lenScale) + adjustedNodePosition = add(proximalNodePosition, jointAdjustDir) + return adjustedNodePosition diff --git a/src/scaffoldmaker/scaffolds.py b/src/scaffoldmaker/scaffolds.py index f6bd8a5b..c5f157ca 100644 --- a/src/scaffoldmaker/scaffolds.py +++ b/src/scaffoldmaker/scaffolds.py @@ -68,7 +68,6 @@ class Scaffolds(object): - _allScaffoldTypes = [ MeshType_1d_bifurcationtree1, MeshType_1d_uterus_network_layout1, diff --git a/src/scaffoldmaker/utils/human_network_layout.py b/src/scaffoldmaker/utils/human_network_layout.py new file mode 100644 index 00000000..596f150b --- /dev/null +++ b/src/scaffoldmaker/utils/human_network_layout.py @@ -0,0 +1,147 @@ +#%% +# Number of elements per segment. Used to calculate the number of nodes per segment. +humanElementCounts = { + 'headElementsCount': 3, + 'neckElementsCount': 2, + 'shoulderElementsCount': 2, + 'brachiumElementsCount': 3, + 'antebrachiumElementsCount': 3, + 'handElementsCount': 1, + 'thoraxElementsCount': 3, + 'abdomenElementsCount': 4, + 'hipElementsCount': 2, + 'upperLegElementsCount': 4, + 'lowerLegElementsCount': 3, + 'footElementsCount': 1 +} + + +def createSegment(nodeCount:int, networkLayout:str, nodeIdentifier:int, endSegment=False, version=0): + """ + Construct a segment of the human network node + + :param nodeCount: Number of nodes to add. + :param networkLayout: String containing the current network layout. + :param nodeIdentifier: Integer denoting the current node. + :param endSegment: If true, adds a comma at the end of the segment. + :param version: If > 0, adds version number on the last node of the segment. + :return networklayout: String containing the updated layout. + :return nodeIdentifier: The updated nodeIdentifier after adding the segment. + """ + for i in range(nodeCount): + networkLayout = networkLayout + str(nodeIdentifier) + if i == nodeCount - 1: + if endSegment: + if version == 0: + segmentConnector = ',' + else: + segmentConnector = '.' + str(version) + ',' + networkLayout = networkLayout + segmentConnector + else: + networkLayout = networkLayout + '-' + nodeIdentifier += 1 + else: + networkLayout = networkLayout + '-' + nodeIdentifier += 1 + return networkLayout, nodeIdentifier + +def constructNetworkLayoutStructure(humanElementCounts:dict): + """ + Construct the network layout of the human wholebody scaffold. + The network layout consists of the following segments: + head, neck, thorax, abdomen, right/left arm, right/left leg. + Arms are subdivided into brachium, antebrachium and hand. + Legs are subdivided into upper leg, lower leg and foot. + + :param humanElementCounts: Dictionary containing the number of elements + corresponding to each segment. + :return humanNetworkLayout: String containing the network layout + """ + nodeIdentifier = 1 + # Head + headNetworkLayout = str(nodeIdentifier) + '-' + nodeIdentifier += 1 + headNetworkLayout, nodeIdentifier = createSegment( + humanElementCounts['headElementsCount'], + headNetworkLayout, nodeIdentifier, endSegment=True) + # Neck + necknNetworkLayout = str(nodeIdentifier) + '-' + nodeIdentifier += 1 + necknNetworkLayout, nodeIdentifier = createSegment( + humanElementCounts['neckElementsCount'], + necknNetworkLayout, nodeIdentifier, endSegment=True, version=1) + neckJointNode = nodeIdentifier + # Thorax + thoraxNetworkLayout = str(nodeIdentifier) + '.1-' + nodeIdentifier += 1 + thoraxNetworkLayout, nodeIdentifier = createSegment( + humanElementCounts['thoraxElementsCount'], + thoraxNetworkLayout, nodeIdentifier, endSegment=True) + # Abdomen + abdomenNetworkLayout = str(nodeIdentifier) + '-' + nodeIdentifier += 1 + abdomenNetworkLayout, nodeIdentifier = createSegment( + humanElementCounts['abdomenElementsCount'], + abdomenNetworkLayout, nodeIdentifier, endSegment=True, version=1) + pelvisNodeJoint = nodeIdentifier + # Arms + arms = [] + for i in range(2): + version = 2 if (i == 0) else 3 #Left is 2, right is 3 + armNetworkLayout = str(neckJointNode) + '.' + str(version) + '-' + nodeIdentifier += 1 + # Shoulder + armNetworkLayout, nodeIdentifier = createSegment( + humanElementCounts['shoulderElementsCount'], + armNetworkLayout, nodeIdentifier, endSegment=True) + # Brachium + armNetworkLayout = armNetworkLayout + str(nodeIdentifier) + '-' + nodeIdentifier += 1 + armNetworkLayout, nodeIdentifier = createSegment( + humanElementCounts['brachiumElementsCount'], + armNetworkLayout, nodeIdentifier, endSegment=True) + # Antebrachium + armNetworkLayout = armNetworkLayout + str(nodeIdentifier) + '-' + nodeIdentifier += 1 + armNetworkLayout, nodeIdentifier = createSegment( + humanElementCounts['antebrachiumElementsCount'], + armNetworkLayout, nodeIdentifier, endSegment=True) + # Hand + armNetworkLayout = armNetworkLayout + str(nodeIdentifier) + '-' + nodeIdentifier += 1 + armNetworkLayout, nodeIdentifier = createSegment( + humanElementCounts['handElementsCount'], + armNetworkLayout, nodeIdentifier, endSegment=True) + arms.append(armNetworkLayout) + #Legs + legs = [] + for i in range(2): + version = 2 if (i == 0) else 3 #Left is 2, right is 3 + legNetworkLayout = str(pelvisNodeJoint) + '.' + str(version) + '-' + nodeIdentifier += 1 + # Hip + legNetworkLayout, nodeIdentifier = createSegment( + humanElementCounts['hipElementsCount'], + legNetworkLayout, nodeIdentifier, endSegment=True) + # Upper leg + legNetworkLayout = legNetworkLayout + str(nodeIdentifier) + '-' + nodeIdentifier += 1 + legNetworkLayout, nodeIdentifier = createSegment( + humanElementCounts['upperLegElementsCount'], + legNetworkLayout, nodeIdentifier, endSegment=True) + # Lower leg + legNetworkLayout = legNetworkLayout + str(nodeIdentifier) + '-' + nodeIdentifier += 1 + legNetworkLayout, nodeIdentifier = createSegment( + humanElementCounts['lowerLegElementsCount'], + legNetworkLayout, nodeIdentifier, endSegment=True) + # Foot + legNetworkLayout = legNetworkLayout + str(nodeIdentifier) + '-' + nodeIdentifier += 1 + legNetworkLayout, nodeIdentifier = createSegment( + humanElementCounts['footElementsCount'], + legNetworkLayout, nodeIdentifier, endSegment=True) + legs.append(legNetworkLayout) + humanNetworkLayout = headNetworkLayout + necknNetworkLayout + arms[0] + arms[1] + thoraxNetworkLayout + abdomenNetworkLayout + legs[0] + legs[1] + humanNetworkLayout = humanNetworkLayout[:-1] #Remove a comma at the end + return humanNetworkLayout \ No newline at end of file diff --git a/src/scaffoldmaker/utils/networkmesh.py b/src/scaffoldmaker/utils/networkmesh.py index 039c7476..efc9b546 100644 --- a/src/scaffoldmaker/utils/networkmesh.py +++ b/src/scaffoldmaker/utils/networkmesh.py @@ -168,7 +168,7 @@ def isPatch(self): """ return self._isPatch - def split(self, splitNetworkNode): + def split(self, splitNetworkNode, isPatch): """ Split segment to finish at splitNetworkNode, returning remainder as a new NetworkSegment. :param splitNetworkNode: NetworkNode to split segment at. @@ -176,7 +176,7 @@ def split(self, splitNetworkNode): """ index = self._networkNodes.index(splitNetworkNode, 1, -1) # throws exception if not an interior node splitNetworkNode.setInteriorSegment(None) - nextSegment = NetworkSegment(self._networkNodes[index:], self._nodeVersions[index:]) + nextSegment = NetworkSegment(self._networkNodes[index:], self._nodeVersions[index:], isPatch) self._networkNodes = self._networkNodes[:index + 1] self._nodeVersions = self._nodeVersions[:index + 1] return nextSegment @@ -250,7 +250,7 @@ def build(self, structureString): if networkNode: interiorSegment = networkNode.getInteriorSegment() if interiorSegment: - nextSegment = interiorSegment.split(networkNode) + nextSegment = interiorSegment.split(networkNode, isPatch) index = self._networkSegments.index(interiorSegment) + 1 self._networkSegments.insert(index, nextSegment) else: diff --git a/tests/test_wholebody2.py b/tests/test_wholebody2.py index d97e33ac..45017ceb 100644 --- a/tests/test_wholebody2.py +++ b/tests/test_wholebody2.py @@ -34,15 +34,18 @@ def test_wholebody2_core(self): parameterSetNames = scaffold.getParameterSetNames() self.assertEqual(parameterSetNames, ["Default", "Human 1 Coarse", "Human 1 Medium", "Human 1 Fine"]) options = scaffold.getDefaultOptions("Human 1 Coarse") - self.assertEqual(19, len(options)) + self.assertEqual(23, len(options)) self.assertEqual(2, options["Number of elements along head"]) self.assertEqual(1, options["Number of elements along neck"]) self.assertEqual(2, options["Number of elements along thorax"]) self.assertEqual(2, options["Number of elements along abdomen"]) - self.assertEqual(5, options["Number of elements along arm to hand"]) + self.assertEqual(3, options["Number of elements along brachium"]) + self.assertEqual(2, options["Number of elements along antebrachium"]) self.assertEqual(1, options["Number of elements along hand"]) - self.assertEqual(4, options["Number of elements along leg to foot"]) - self.assertEqual(2, options["Number of elements along foot"]) + self.assertEqual(2, options["Number of elements along hip"]) + self.assertEqual(3, options["Number of elements along upper leg"]) + self.assertEqual(3, options["Number of elements along lower leg"]) + self.assertEqual(1, options["Number of elements along foot"]) self.assertEqual(12, options["Number of elements around head"]) self.assertEqual(12, options["Number of elements around torso"]) self.assertEqual(8, options["Number of elements around arm"]) @@ -57,18 +60,18 @@ def test_wholebody2_core(self): region = context.getDefaultRegion() self.assertTrue(region.isValid()) annotationGroups = scaffold.generateMesh(region, options)[0] - self.assertEqual(32, len(annotationGroups)) + self.assertEqual(58, len(annotationGroups)) fieldmodule = region.getFieldmodule() self.assertEqual(RESULT_OK, fieldmodule.defineAllFaces()) mesh3d = fieldmodule.findMeshByDimension(3) - self.assertEqual(704, mesh3d.getSize()) + self.assertEqual(904, mesh3d.getSize()) mesh2d = fieldmodule.findMeshByDimension(2) - self.assertEqual(2306, mesh2d.getSize()) + self.assertEqual(2946, mesh2d.getSize()) mesh1d = fieldmodule.findMeshByDimension(1) - self.assertEqual(2533, mesh1d.getSize()) + self.assertEqual(3223, mesh1d.getSize()) nodes = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) - self.assertEqual(932, nodes.getSize()) + self.assertEqual(1195, nodes.getSize()) datapoints = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_DATAPOINTS) self.assertEqual(0, datapoints.getSize()) @@ -77,8 +80,8 @@ def test_wholebody2_core(self): self.assertTrue(coordinates.isValid()) minimums, maximums = evaluateFieldNodesetRange(coordinates, nodes) tol = 1.0E-4 - assertAlmostEqualList(self, minimums, [0.0, -3.650833433150559, -1.25], tol) - assertAlmostEqualList(self, maximums, [20.48332368641937, 3.650833433150559, 2.15], tol) + assertAlmostEqualList(self, minimums, [0.0, -3.936660189011623, -1.25], tol) + assertAlmostEqualList(self, maximums, [19.725896216328984, 3.936660189011623, 2.354767205067949], tol) with ChangeManager(fieldmodule): one = fieldmodule.createFieldConstant(1.0) @@ -96,17 +99,17 @@ def test_wholebody2_core(self): result, surfaceArea = surfaceAreaField.evaluateReal(fieldcache, 1) self.assertEqual(result, RESULT_OK) - self.assertAlmostEqual(volume, 98.41184838749453, delta=tol) - self.assertAlmostEqual(surfaceArea, 228.9017722286146, delta=tol) + self.assertAlmostEqual(volume, 93.46891042325106, delta=tol) + self.assertAlmostEqual(surfaceArea, 224.0370925379395, delta=tol) # check some annotation groups: expectedSizes3d = { - "abdominal cavity": (40, 10.074522362520469), - "core": (428, 45.535080392793354), - "head": (64, 6.909618374858558), - "thoracic cavity": (40, 6.974341918899202), - "shell": (276, 52.878054197629744) + "abdominal cavity": (40, 9.250133376720713), + "core": (548, 42.72181578233019), + "head": (64, 6.909618374858556), + "thoracic cavity": (40, 6.858627377211818), + "shell": (356, 50.74698672407125) } for name in expectedSizes3d: term = get_body_term(name) @@ -122,14 +125,14 @@ def test_wholebody2_core(self): self.assertAlmostEqual(volume, expectedSizes3d[name][1], delta=tol) expectedSizes2d = { - "abdominal cavity boundary surface": (64, 27.428203203724674), + "abdominal cavity boundary surface": (64, 25.710221456528195), "diaphragm": (20, 3.0778936559347208), - "left upper limb skin epidermis outer surface": (68, 22.433540462588258), - "left lower limb skin epidermis outer surface": (68, 55.24949927903832), - "right upper limb skin epidermis outer surface": (68, 22.433540462588258), - "right lower limb skin epidermis outer surface": (68, 55.24949927903832), - "skin epidermis outer surface": (388, 228.9017722286146), - "thoracic cavity boundary surface": (64, 20.2627556651649) + "left upper limb skin epidermis outer surface": (56, 18.880989688740083), + "left lower limb skin epidermis outer surface": (64, 48.27786366921722), + "right upper limb skin epidermis outer surface": (56, 18.880989688740083), + "right lower limb skin epidermis outer surface": (64, 48.27786366921722), + "skin epidermis outer surface": (468, 224.0370925379395), + "thoracic cavity boundary surface": (64, 19.93080798484531) } for name in expectedSizes2d: term = get_body_term(name) @@ -145,7 +148,7 @@ def test_wholebody2_core(self): self.assertAlmostEqual(surfaceArea, expectedSizes2d[name][1], delta=tol) expectedSizes1d = { - "spinal cord": (8, 10.85369421002332) + "spinal cord": (8, 10.154987441354445) } for name in expectedSizes1d: term = get_body_term(name) @@ -172,18 +175,18 @@ def test_wholebody2_tube(self): region = context.getDefaultRegion() self.assertTrue(region.isValid()) annotationGroups = scaffold.generateMesh(region, options)[0] - self.assertEqual(24, len(annotationGroups)) + self.assertEqual(50, len(annotationGroups)) fieldmodule = region.getFieldmodule() self.assertEqual(RESULT_OK, fieldmodule.defineAllFaces()) mesh3d = fieldmodule.findMeshByDimension(3) - self.assertEqual(276, mesh3d.getSize()) + self.assertEqual(356, mesh3d.getSize()) mesh2d = fieldmodule.findMeshByDimension(2) - self.assertEqual(1126, mesh2d.getSize()) + self.assertEqual(1446, mesh2d.getSize()) mesh1d = fieldmodule.findMeshByDimension(1) - self.assertEqual(1443, mesh1d.getSize()) + self.assertEqual(1843, mesh1d.getSize()) nodes = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) - self.assertEqual(590, nodes.getSize()) + self.assertEqual(763, nodes.getSize()) datapoints = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_DATAPOINTS) self.assertEqual(0, datapoints.getSize()) @@ -192,8 +195,8 @@ def test_wholebody2_tube(self): self.assertTrue(coordinates.isValid()) minimums, maximums = evaluateFieldNodesetRange(coordinates, nodes) tol = 1.0E-4 - assertAlmostEqualList(self, minimums, [0.0, -3.650833433150559, -1.25], tol) - assertAlmostEqualList(self, maximums, [20.48332368641937, 3.650833433150559, 2.15], tol) + assertAlmostEqualList(self, minimums, [0.0, -3.936660189011623, -1.25], tol) + assertAlmostEqualList(self, maximums, [19.725896216328984, 3.936660189011623, 2.354767205067949], tol) with ChangeManager(fieldmodule): one = fieldmodule.createFieldConstant(1.0) @@ -220,13 +223,13 @@ def test_wholebody2_tube(self): result, innerSurfaceArea = innerSurfaceAreaField.evaluateReal(fieldcache, 1) self.assertEqual(result, RESULT_OK) - self.assertAlmostEqual(volume, 52.87781598884186, delta=tol) - self.assertAlmostEqual(outerSurfaceArea, 224.9456647093771, delta=tol) - self.assertAlmostEqual(innerSurfaceArea, 155.4114878443358, delta=tol) + self.assertAlmostEqual(volume, 50.747766334772045, delta=tol) + self.assertAlmostEqual(outerSurfaceArea, 220.08098501870208, delta=tol) + self.assertAlmostEqual(innerSurfaceArea, 151.26319776665184, delta=tol) # check some annotationGroups: expectedSizes2d = { - "skin epidermis outer surface": (320, 228.11749988635236) + "skin epidermis outer surface": (400, 223.25282019567732) } for name in expectedSizes2d: term = get_body_term(name)