Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Improve body shape for fitting #267

Merged
merged 2 commits into from
Oct 29, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
96 changes: 64 additions & 32 deletions src/scaffoldmaker/meshtypes/meshtype_3d_wholebody2.py
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
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