Skip to content

Commit

Permalink
Merge pull request #267 from rchristie/body-shape
Browse files Browse the repository at this point in the history 
Improve body shape for fitting
  • Loading branch information
@rchristie
rchristie authored Oct 29, 2024
2 parents 9ad9121 + 57fe52f commit 2665486
Show file tree
Hide file tree
Showing 2 changed files with 87 additions and 55 deletions.
96 changes: 64 additions & 32 deletions src/scaffoldmaker/meshtypes/meshtype_3d_wholebody2.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -50,15 +50,16 @@ def getDefaultOptions(cls, parameterSetName="Default"):
options["Head length"] = 2.2
options["Head width"] = 2.0
options["Neck length"] = 1.3
options["Shoulder drop"] = 0.7
options["Shoulder drop"] = 1.0
options["Shoulder width"] = 4.5
options["Arm lateral angle degrees"] = 10.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["Hand length"] = 1.5
options["Hand thickness"] = 0.3
options["Hand length"] = 2.0
options["Hand thickness"] = 0.2
options["Hand width"] = 1.0
options["Thorax length"] = 2.5
options["Abdomen length"] = 3.0
Expand All @@ -70,7 +71,7 @@ def getDefaultOptions(cls, parameterSetName="Default"):
options["Leg length"] = 10.0
options["Leg top diameter"] = 2.0
options["Leg bottom diameter"] = 0.7
options["Foot height"] = 1.0
options["Foot height"] = 1.25
options["Foot length"] = 2.5
options["Foot thickness"] = 0.3
options["Foot width"] = 1.0
Expand All @@ -90,6 +91,7 @@ def getOrderedOptionNames(cls):
"Arm lateral angle degrees",
"Arm length",
"Arm top diameter",
"Arm twist angle degrees",
"Wrist thickness",
"Wrist width",
"Hand length",
Expand Down Expand Up @@ -154,20 +156,15 @@ def checkOptions(cls, options):
options[key] = 0.1
elif options[key] > 0.9:
options[key] = 0.9
for key in [
"Arm lateral angle degrees"
]:
if options[key] < -60.0:
options[key] = -60.0
elif options[key] > 200.0:
options[key] = 200.0
for key in [
"Leg lateral angle degrees"
]:
if options[key] < -20.0:
options[key] = -20.0
elif options[key] > 60.0:
options[key] = 60.0
for key, angleRange in {
"Arm lateral angle degrees": (-60.0, 200.0),
"Arm twist angle degrees": (-90.0, 90.0),
"Leg lateral angle degrees": (-20.0, 60.0)
}.items():
if options[key] < angleRange[0]:
options[key] = angleRange[0]
elif options[key] > angleRange[1]:
options[key] = angleRange[1]
return dependentChanges

@classmethod
Expand All @@ -189,6 +186,7 @@ def generateBaseMesh(cls, region, options):
armAngleRadians = math.radians(options["Arm lateral angle degrees"])
armLength = options["Arm length"]
armTopRadius = 0.5 * options["Arm top diameter"]
armTwistAngleRadians = math.radians(options["Arm twist angle degrees"])
halfWristThickness = 0.5 * options["Wrist thickness"]
halfWristWidth = 0.5 * options["Wrist width"]
handLength = options["Hand length"]
Expand Down Expand Up @@ -396,22 +394,23 @@ def generateBaseMesh(cls, region, options):
# 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(2.0 * shoulderDrop / halfShoulderWidth)
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)
hd3 = [0.0, 0.0, halfHandWidth]
hid3 = mult(hd3, innerProportionDefault)
for side in (left, right):
armAngle = armAngleRadians if (side == left) else -armAngleRadians
cosArmAngle = math.cos(armAngle)
sinArmAngle = math.sin(armAngle)
armStartY = (halfShoulderWidth if (side == left) else -halfShoulderWidth) * math.cos(shoulderAngleRadians)
x = [armStartX, armStartY, 0.0]
d1 = [armScale * cosArmAngle, armScale * sinArmAngle, 0.0]
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
sd1 = interpolateLagrangeHermiteDerivative(sx, x, d1, 0.0)
nx, nd1 = sampleCubicHermiteCurvesSmooth([sx, x], [sd1, d1], 2, derivativeMagnitudeEnd=armScale)[0:2]
Expand Down Expand Up @@ -454,22 +453,43 @@ 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)
d12 = [-d12_mag * sinArmAngle, d12_mag * cosArmAngle, 0.0]
id12 = mult(d12, innerProportionDefault)
d13 = [0.0, 0.0, d13_mag]
id13 = mult(d13, innerProportionDefault)
# main part of arm to wrist
elementTwistAngle = ((armTwistAngleRadians if (side == left) else -armTwistAngleRadians) /
(armToHandElementsCount - 3))
for i in range(armToHandElementsCount - 1):
xi = i / (armToHandElementsCount - 2)
node = nodes.findNodeByIdentifier(nodeIdentifier)
fieldcache.setNode(node)
x = [armStartX + d1[0] * i, armStartY + d1[1] * i, d1[2] * i]
halfThickness = xi * halfWristThickness + (1.0 - xi) * armTopRadius
halfWidth = xi * halfWristWidth + (1.0 - xi) * armTopRadius
d2 = [-halfThickness * sinArmAngle, halfThickness * cosArmAngle, 0.0]
d3 = [0.0, 0.0, halfWidth]
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:
d2 = mult(armSide, halfThickness)
d3 = mult(armFront, halfWidth)
d12 = mult(armSide, d12_mag)
d13 = mult(armFront, d13_mag)
else:
cosTwistAngle = math.cos(twistAngle)
sinTwistAngle = math.sin(twistAngle)
d2 = sub(mult(armSide, halfThickness * cosTwistAngle),
mult(armFront, halfThickness * sinTwistAngle))
d3 = add(mult(armFront, halfWidth * cosTwistAngle),
mult(armSide, halfWidth * sinTwistAngle))
d12 = set_magnitude(d2, d12_mag)
d13 = set_magnitude(d3, d13_mag)
if i < (armToHandElementsCount - 2):
d12 = add(d12, set_magnitude(d3, -halfThickness * elementTwistAngle))
d13 = add(d13, set_magnitude(d2, halfWidth * elementTwistAngle))
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
Expand All @@ -479,8 +499,19 @@ def generateBaseMesh(cls, region, options):
fieldcache.setNode(node)
hx = [armStartX + armLength * cosArmAngle, armStartY + armLength * sinArmAngle, 0.0]
hd1 = computeCubicHermiteEndDerivative(x, d1, hx, d1)
hd2 = set_magnitude(d2, halfHandThickness)
twistAngle = armTwistAngleRadians if (side == left) else -armTwistAngleRadians
if twistAngle == 0.0:
hd2 = set_magnitude(d2, halfHandThickness)
hd3 = [0.0, 0.0, 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))
hid2 = mult(hd2, innerProportionDefault)
hid3 = mult(hd3, innerProportionDefault)
setNodeFieldParameters(coordinates, fieldcache, hx, hd1, hd2, hd3)
setNodeFieldParameters(innerCoordinates, fieldcache, hx, hd1, hid2, hid3)
nodeIdentifier += 1
Expand Down Expand Up @@ -537,16 +568,17 @@ def generateBaseMesh(cls, region, options):
nodeIdentifier += 1
# foot
fx = [x,
add(add(legStart, mult(legDirn, legLength - halfFootThickness)),
[0.0, 0.0, 0.25 * footLength + legBottomRadius]),
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), legBottomRadius),
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])
Expand Down
46 changes: 23 additions & 23 deletions tests/test_wholebody2.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -69,8 +69,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.7000751482231564, -1.25], tol)
assertAlmostEqualList(self, maximums, [20.437483381451223, 3.7000751482231564, 2.15], tol)
assertAlmostEqualList(self, minimums, [0.0, -3.650833433150559, -1.25], tol)
assertAlmostEqualList(self, maximums, [20.48332368641937, 3.650833433150559, 2.15], tol)

with ChangeManager(fieldmodule):
one = fieldmodule.createFieldConstant(1.0)
Expand All @@ -88,17 +88,17 @@ def test_wholebody2_core(self):
result, surfaceArea = surfaceAreaField.evaluateReal(fieldcache, 1)
self.assertEqual(result, RESULT_OK)

self.assertAlmostEqual(volume, 98.99087587225652, delta=tol)
self.assertAlmostEqual(surfaceArea, 229.8973868830034, delta=tol)
self.assertAlmostEqual(volume, 98.41184838749453, delta=tol)
self.assertAlmostEqual(surfaceArea, 228.9017722286146, delta=tol)

# check some annotation groups:

expectedSizes3d = {
"abdominal cavity": (40, 10.081327011840981),
"core": (428, 45.78886468970665),
"abdominal cavity": (40, 10.074522362520469),
"core": (428, 45.535080392793354),
"head": (64, 6.909618374858558),
"thoracic cavity": (40, 7.135491643165788),
"shell": (276, 53.20466827197639)
"thoracic cavity": (40, 6.974341918899202),
"shell": (276, 52.878054197629744)
}
for name in expectedSizes3d:
term = get_body_term(name)
Expand All @@ -114,14 +114,14 @@ def test_wholebody2_core(self):
self.assertAlmostEqual(volume, expectedSizes3d[name][1], delta=tol)

expectedSizes2d = {
"abdominal cavity boundary": (64, 27.46017763836879),
"abdominal cavity boundary": (64, 27.428203203724674),
"diaphragm": (20, 3.0778936559347208),
"left arm skin epidermis": (68, 22.627795339108015),
"left leg skin epidermis": (68, 55.21582811667045),
"right arm skin epidermis": (68, 22.62779533911023),
"right leg skin epidermis": (68, 55.21582811667045),
"skin epidermis": (388, 229.8973868830034),
"thoracic cavity boundary": (64, 20.606925296069125)
"left arm skin epidermis": (68, 22.433540462588258),
"left leg skin epidermis": (68, 55.24949927903832),
"right arm skin epidermis": (68, 22.433540462588258),
"right leg skin epidermis": (68, 55.24949927903832),
"skin epidermis": (388, 228.9017722286146),
"thoracic cavity boundary": (64, 20.2627556651649)
}
for name in expectedSizes2d:
term = get_body_term(name)
Expand All @@ -137,7 +137,7 @@ def test_wholebody2_core(self):
self.assertAlmostEqual(surfaceArea, expectedSizes2d[name][1], delta=tol)

expectedSizes1d = {
"spinal cord": (8, 10.856804626156244)
"spinal cord": (8, 10.85369421002332)
}
for name in expectedSizes1d:
term = get_body_term(name)
Expand Down Expand Up @@ -183,9 +183,9 @@ def test_wholebody2_tube(self):
coordinates = fieldmodule.findFieldByName("coordinates").castFiniteElement()
self.assertTrue(coordinates.isValid())
minimums, maximums = evaluateFieldNodesetRange(coordinates, nodes)
tol = 1.0E-6
assertAlmostEqualList(self, minimums, [0.0, -3.7000751482231564, -1.25], tol)
assertAlmostEqualList(self, maximums, [20.437483381451223, 3.7000751482231564, 2.15], tol)
tol = 1.0E-4
assertAlmostEqualList(self, minimums, [0.0, -3.650833433150559, -1.25], tol)
assertAlmostEqualList(self, maximums, [20.48332368641937, 3.650833433150559, 2.15], tol)

with ChangeManager(fieldmodule):
one = fieldmodule.createFieldConstant(1.0)
Expand All @@ -212,13 +212,13 @@ def test_wholebody2_tube(self):
result, innerSurfaceArea = innerSurfaceAreaField.evaluateReal(fieldcache, 1)
self.assertEqual(result, RESULT_OK)

self.assertAlmostEqual(volume, 53.20377108353156, delta=tol)
self.assertAlmostEqual(outerSurfaceArea, 225.79241553960935, delta=tol)
self.assertAlmostEqual(innerSurfaceArea, 155.88335089354402, delta=tol)
self.assertAlmostEqual(volume, 52.87781598884186, delta=tol)
self.assertAlmostEqual(outerSurfaceArea, 224.9456647093771, delta=tol)
self.assertAlmostEqual(innerSurfaceArea, 155.4114878443358, delta=tol)

# check some annotationGroups:
expectedSizes2d = {
"skin epidermis": (320, 229.04017126690428)
"skin epidermis": (320, 228.11749988635236)
}
for name in expectedSizes2d:
term = get_body_term(name)
Expand Down

0 comments on commit 2665486

Please sign in to comment.