Skip to content

Commit

Permalink
Merge pull request #92 from s-fong/annotation
Browse files Browse the repository at this point in the history 
Annotation
  • Loading branch information
@rchristie
rchristie authored Sep 24, 2020
2 parents 059ad2f + ee7db4f commit c0ecb19
Show file tree
Hide file tree
Showing 2 changed files with 198 additions and 36 deletions.
19 changes: 19 additions & 0 deletions src/scaffoldmaker/annotation/stellate_terms.py
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,19 @@
"""
Common resource for stellate annotation terms.
"""

# convention: preferred name, preferred id, followed by any other ids and alternative names
stellate_terms = [
( "cervicothoracic ganglion", "UBERON:2441", "FMA:6469", "ILX:733799")
]

def get_stellate_term(name : str):
"""
Find term by matching name to any identifier held for a term.
Raise exception if name not found.
:return ( preferred name, preferred id )
"""
for term in stellate_terms:
if name in term:
return ( term[0], term[1] )
raise NameError("Stellate annotation term '" + name + "' not found.")
215 changes: 179 additions & 36 deletions src/scaffoldmaker/meshtypes/meshtype_3d_stellate1.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,10 +4,12 @@

from __future__ import division
import math
from opencmiss.utils.zinc.field import findOrCreateFieldCoordinates
from opencmiss.utils.zinc.field import findOrCreateFieldCoordinates, findOrCreateFieldGroup, findOrCreateFieldNodeGroup, findOrCreateFieldStoredMeshLocation, findOrCreateFieldStoredString
from opencmiss.zinc.element import Element
from opencmiss.zinc.field import Field
from opencmiss.zinc.node import Node
from scaffoldmaker.annotation.annotationgroup import AnnotationGroup, findOrCreateAnnotationGroupForTerm, getAnnotationGroupForTerm
from scaffoldmaker.annotation.stellate_terms import get_stellate_term
from scaffoldmaker.meshtypes.scaffold_base import Scaffold_base
from scaffoldmaker.utils.eftfactory_bicubichermitelinear import eftfactory_bicubichermitelinear
from scaffoldmaker.utils.meshrefinement import MeshRefinement
Expand All @@ -25,22 +27,39 @@ def getName():
return '3D Stellate 1'

@staticmethod
def getDefaultOptions(parameterSetName='Default'):
return {
'Numbers of elements along arms' : [4,2,2],
'Element length along arm' : 1.0,
'Element width across arm' : 0.5,
'Element thickness' : 0.5,
'Refine' : False,
'Refine number of elements along arm' : 1,
'Refine number of elements across arm' : 1,
'Refine number of elements through thickness' : 1
}
def getParameterSetNames():
return [
'Default',
'Mouse cervicothoracic ganglion 1']

@classmethod
def getDefaultOptions(cls, parameterSetName='Default'):
options = {}
options['Base parameter set'] = parameterSetName

isMouse = 'Mouse' in parameterSetName

if isMouse:
options['Numbers of elements along arms'] = [4,2,2]
else:
options['Numbers of elements along arms'] = [4,2,2]
options['Element width central'] = 0.8
options['Element length along arm'] = 0.8
options['Element width across arm'] = 0.3
options['Element thickness'] = 0.1
options['Refine'] = False
options['Refine number of elements along arm'] = 1
options['Refine number of elements across arm'] = 1
options['Refine number of elements through thickness'] = 1
return options


@staticmethod
def getOrderedOptionNames():

return [
'Numbers of elements along arms',
'Element width central',
'Element length along arm',
'Element width across arm',
'Element thickness',
Expand All @@ -50,8 +69,8 @@ def getOrderedOptionNames():
'Refine number of elements through thickness'
]

@staticmethod
def checkOptions(options):
@classmethod
def checkOptions(cls, options):
for key in [
'Refine number of elements along arm',
'Refine number of elements across arm',
Expand All @@ -60,6 +79,7 @@ def checkOptions(options):
options[key] = 1

for key in [
'Element width central',
'Element length along arm',
'Element width across arm',
'Element thickness']:
Expand All @@ -85,7 +105,11 @@ def generateBaseMesh(cls, region, options):
:param options: Dict containing options. See getDefaultOptions().
:return: None
"""
isDefault = 'Default' in options['Base parameter set']
isMouse = 'Mouse' in options['Base parameter set']

armCount = 3
elementLengthCentral = options['Element width central']
elementLengths = [options['Element length along arm'],
options['Element width across arm'],
options['Element thickness']]
Expand All @@ -95,9 +119,9 @@ def generateBaseMesh(cls, region, options):
useCrossDerivatives = False

fm = region.getFieldmodule()
nodes = fm.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES)
coordinates = findOrCreateFieldCoordinates(fm)

nodes = fm.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES)
nodetemplate = nodes.createNodetemplate()
nodetemplate.defineField(coordinates)
nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_VALUE, 1)
Expand All @@ -108,27 +132,86 @@ def generateBaseMesh(cls, region, options):
mesh = fm.findMeshByDimension(3)
cache = fm.createFieldcache()

markerGroup = findOrCreateFieldGroup(fm, "marker")
markerName = findOrCreateFieldStoredString(fm, name="marker_name")
markerLocation = findOrCreateFieldStoredMeshLocation(fm, mesh, name="marker_location")

markerPoints = findOrCreateFieldNodeGroup(markerGroup, nodes).getNodesetGroup()
markerTemplateInternal = nodes.createNodetemplate()
markerTemplateInternal.defineField(markerName)
markerTemplateInternal.defineField(markerLocation)

# markers with element number and xi position
allMarkers = {}
if isMouse:
xProportion = {}
xProportion['ICN'] = 0.9
xProportion['VA'] = 0.9
xProportion['DA'] = 0.9
xProportion['C8'] = 0.9
xProportion['T1'] = 0.25
xProportion['T2'] = 0.5
xProportion['T3'] = 0.75
xProportion['TST'] = 1
armNumber = {}
armNumber['ICN'] = 2
armNumber['VA'] = 2
armNumber['DA'] = 3
armNumber['C8'] = 3
armNumber['T1'] = 1
armNumber['T2'] = 1
armNumber['T3'] = 1
armNumber['TST'] = 1
nerveAbbrev = list(xProportion.keys())
elementIndex = {}
xi1 = {}
for nerve in nerveAbbrev:
elementIndex[nerve] = int(xProportion[nerve] * elementsCount1[armNumber[nerve]-1])
xi1[nerve] = 1 if xProportion[nerve] == 1 else xProportion[nerve] * elementsCount1[armNumber[nerve]-1] - elementIndex[nerve]
elementIndex[nerve] += 1 if xProportion[nerve] < 1 else 0
j = 10

allMarkers = { "Inferior cardiac nerve" : {"elementID": elementIndex['ICN']+2*elementsCount1[0], "xi": [xi1['ICN'], 0.0, 0.5]},
"Ventral ansa subclavia" : {"elementID": elementIndex['VA']+2*elementsCount1[0]+elementsCount1[1], "xi": [xi1['VA'], 1.0, 0.5]},
"Dorsal ansa subclavia" : {"elementID": elementIndex['DA']+2*(elementsCount1[0]+elementsCount1[1]), "xi": [xi1['DA'], 0.0, 0.5]},
"Cervical spinal nerve 8" : {"elementID": elementIndex['C8']+2*(elementsCount1[0]+elementsCount1[1])+elementsCount1[2], "xi": [xi1['C8'], 1.0, 0.5]},
"Thoracic spinal nerve 1" : {"elementID": elementIndex['T1'], "xi": [xi1['T1'], 0.0, 0.5]},
"Thoracic spinal nerve 2" : {"elementID": elementIndex['T2'], "xi": [xi1['T2'], 0.0, 0.5]},
"Thoracic spinal nerve 3" : {"elementID": elementIndex['T3'], "xi": [xi1['T3'], 0.0, 0.5]},
"Thoracic sympathetic nerve trunk" : {"elementID": elementIndex['TST'], "xi": [xi1['TST'], 1.0, 0.5]},
}

# arm group annotations for user
armTerms, _ = getAutomaticArmFaceTerms(armCount)
armGroups = [AnnotationGroup(region, armTerm) for armTerm in armTerms]
stellateTerm = get_stellate_term("cervicothoracic ganglion") if isMouse else ("stellate", None)
stellateGroup = AnnotationGroup(region, stellateTerm)
annotationGroups = [stellateGroup] + armGroups

armMeshGroups = [a.getMeshGroup(mesh) for a in armGroups]
stellateMeshGroup = stellateGroup.getMeshGroup(mesh)

# Create nodes
numNodesPerArm = [0]
halfArmArcAngleRadians = math.pi / armCount
dipMultiplier = 1 #elementLengths[1] * 1.2 * 1.5

dipMultiplier = 1
nodeIdentifier = 1
minArmAngle = 2 * math.pi / armCount
halfArmArcAngleRadians = minArmAngle/2
xx = []
xds1 = []
xds2 = []
x_in_nodes = []
for na in range(armCount):
elementsCount_i = [elementsCount1[na], elementsCount2, elementsCount3]
x, ds1, ds2, nWheelEdge = createArm(halfArmArcAngleRadians, elementLengths, minArmAngle * na, minArmAngle, elementsCount_i, dipMultiplier, armCount, na)
x, ds1, ds2, nWheelEdge = createArm(halfArmArcAngleRadians, elementLengths, elementLengthCentral, elementsCount_i, dipMultiplier, armCount, na)
for ix in range(len(x)):
if na == 0 or ix not in nWheelEdge:
node = nodes.createNode(nodeIdentifier, nodetemplate)
cache.setNode(node)
coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, x[ix])
coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, ds1[ix])
coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, ds2[ix])

nodeIdentifier += 1
x_in_nodes.append(x[ix])
numNodesPerArm.append(len(x))
Expand Down Expand Up @@ -173,13 +256,12 @@ def generateBaseMesh(cls, region, options):
nwPrev[1], no2 + em - 1 + nplUq,
nCentre[1], bni + (4 * em) - 2]
else:
# nplPrev = int(numNodesPerArm[na]/2) - elementsCount1[na-1] - 2
nplPrev = int(numNodesPerArm[na]/2) - 2
no2 = elementsCount1[na]-1
no3 = int(numNodesPerArm[na+1]/2) - 3
nwPrev = [cumNumNodesPerArm[na-1] + 2*(elementsCount1[na-1]),
cumNumNodesPerArm[na-1] + 2*(elementsCount1[na-1]) + nplPrev]
start = cumNumNodesPerArm[na] - 3 # -4 + 1
start = cumNumNodesPerArm[na] - 3
nodeIdentifiers = [nwPrev[0], start,
nCentre[0], start + no2,
nwPrev[1], start + no3,
Expand Down Expand Up @@ -346,7 +428,6 @@ def generateBaseMesh(cls, region, options):
# rounded ends of arms. Collapse xi2 at xi1 = 1
eft1 = bicubichermitelinear.createEftNoCrossDerivatives()
remapEftNodeValueLabel(eft1, [2, 4, 6, 8], Node.VALUE_LABEL_D_DS2, [])
# remapEftNodeValueLabel(eft1, [2, 5], Node.VALUE_LABEL_D_DS2, [])
if e2 == 0:
remapEftNodeValueLabel(eft1, [2, 6], Node.VALUE_LABEL_D_DS1,
[(Node.VALUE_LABEL_D_DS2, [])])
Expand All @@ -367,8 +448,29 @@ def generateBaseMesh(cls, region, options):
element = mesh.createElement(elementIdentifier, elementtemplate1)
result = element.setNodesByIdentifier(eft1, nodeIdentifiers)
result3 = element.setScaleFactors(eft1, scalefactors) if scalefactors else None

# add to meshGroup
stellateMeshGroup.addElement(element)
armMeshGroups[na].addElement(element)

elementIdentifier += 1
return []

# annotation fiducial points
if isMouse:
for key in allMarkers:

xi = allMarkers[key]["xi"]
addMarker = {"name": key, "xi": allMarkers[key]["xi"]}

markerPoint = markerPoints.createNode(nodeIdentifier, markerTemplateInternal)
nodeIdentifier += 1
cache.setNode(markerPoint)
markerName.assignString(cache, addMarker["name"])
elementID = allMarkers[key]["elementID"]
element = mesh.findElementByIdentifier(elementID)
markerLocation.assignMeshLocation(cache, element, addMarker["xi"])

return annotationGroups

@classmethod
def refineMesh(cls, meshrefinement, options):
Expand All @@ -384,16 +486,47 @@ def refineMesh(cls, meshrefinement, options):
meshrefinement.refineAllElementsCubeStandard3d(refineElementsCount1, refineElementsCount2, refineElementsCount3)


def createArm(halfArmArcAngleRadians, elementLengths, armAngle, armAngleConst, elementsCount, dipMultiplier, armCount, armIndex):
@classmethod
def defineFaceAnnotations(cls, region, options, annotationGroups):
"""
Add point annotation groups from the 1D mesh.
:param region: Zinc region containing model.
:param options: Dict containing options. See getDefaultOptions().
:param annotationGroups: List of annotation groups for top-level elements.
New point annotation groups are appended to this list.
"""
# create groups
fm = region.getFieldmodule()
armCount = len(options['Numbers of elements along arms'])
mesh2d = fm.findMeshByDimension(2)
is_exterior = fm.createFieldIsExterior()
is_exterior_face_xi2_0 = fm.createFieldAnd(is_exterior, fm.createFieldIsOnFace(Element.FACE_TYPE_XI2_0))
is_exterior_face_xi2_1 = fm.createFieldAnd(is_exterior, fm.createFieldIsOnFace(Element.FACE_TYPE_XI2_1))
armTerms, faceTerms = getAutomaticArmFaceTerms(armCount)

armGroups = [getAnnotationGroupForTerm(annotationGroups, armTerm) for armTerm in armTerms]
isArm =[armGroup.getFieldElementGroup(mesh2d) for armGroup in armGroups]
for arm in range(armCount):
is_face = fm.createFieldOr(
fm.createFieldAnd(isArm[arm - 1], is_exterior_face_xi2_1),
fm.createFieldAnd(isArm[arm], is_exterior_face_xi2_0))
faceGroup = findOrCreateAnnotationGroupForTerm(annotationGroups, region, faceTerms[arm - 1])
faceGroup.getMeshGroup(mesh2d).addElementsConditional(is_face)

def getAutomaticArmFaceTerms(armCount):
armTerms = [("stellate arm %d"%(i), None) for i in range(1,armCount+1)]
faceTerms = [("stellate face %d-%d" % (i, (i % armCount) + 1), None) for i in range(1, armCount + 1)]
return armTerms, faceTerms

def createArm(halfArmArcAngleRadians, elementLengths, elementLengthCentral, elementsCount, dipMultiplier, armCount, armIndex):
"""
Create single arm unit.
Base length of element is 1.
Direction: anticlockwise.
Minimum arm length is 2 elements.
:param halfArmArcAngleRadians: angle arising from base node (rad)
:param elementLengths: [Element length along arm, half Element width across arm, Element thickness]
:param armAngle: angle from x axis of the arm about origin (rad)
:param armAngleConst: minimum angle from x axis of the first arm about origin (rad)
:param elementLengthCentral: Radial element length about central node.
:param elementsCount: list of numbers of elements along arm length, across width and through thickness directions
:param dipMultiplier: factor that wheel nodes protrude by, relative to unit length
:param armCount: number of arms in body
Expand All @@ -408,11 +541,14 @@ def createArm(halfArmArcAngleRadians, elementLengths, armAngle, armAngleConst, e
elementsCount1, elementsCount2, elementsCount3 = elementsCount
[elementLength, elementWidth, elementHeight] = elementLengths

shorterArmEnd = True
armAngle = 2*halfArmArcAngleRadians*armIndex
armAngleConst = 2*halfArmArcAngleRadians

xnodes_ds1 = []
xnodes_ds2 = []
dx_ds1 = [elementLength, 0.0, 0.0]
dx_ds2 = [0.0, elementWidth, 0.0]
wheelDvMult = [0.5, 0.75]

nodes_per_layer = (elementsCount1 + 1) * (elementsCount2 + 1) - 2 # discount 2 armEnd corners
x = []
Expand All @@ -427,14 +563,14 @@ def createArm(halfArmArcAngleRadians, elementLengths, armAngle, armAngleConst, e
else:
rmVertexNodes = nCentre + [0, nodes_per_layer,
2* (elementsCount1+1) - 1, 2* (elementsCount1+1) - 1 + nodes_per_layer ]
dcent = [elementLength * math.cos(armAngleConst / 2), elementLength * math.sin(armAngleConst / 2), 0.0]
dvertex = [elementLength * dipMultiplier * math.cos(halfArmArcAngleRadians),
elementLength * dipMultiplier * math.sin(halfArmArcAngleRadians)]
dcent = [elementLengthCentral * math.cos(armAngleConst / 2), elementLengthCentral * math.sin(armAngleConst / 2), 0.0]
dvertex = [elementLengthCentral * dipMultiplier * math.cos(halfArmArcAngleRadians),
elementLengthCentral * dipMultiplier * math.sin(halfArmArcAngleRadians)]

dipLength = 0.5*(elementLength + elementWidth)*dipMultiplier
dipLength = 0.5*(elementLengthCentral + elementWidth)*dipMultiplier
dvertex = [dipLength * math.cos(halfArmArcAngleRadians),
dipLength * math.sin(halfArmArcAngleRadians)]
dipMag = 2*dipLength - elementLength
dipMag = 2*dipLength - elementLengthCentral

nid = 0
for e3 in range(elementsCount3 + 1):
Expand All @@ -452,8 +588,11 @@ def createArm(halfArmArcAngleRadians, elementLengths, armAngle, armAngleConst, e
x1 = dvertex[0]
x2 = dvertex[1]
else:
x1 = (elementLength*e1)
x2 = ((e2 - 1) * elementWidth )
if e1 == elementsCount1 and shorterArmEnd:
x1 = 0.5*(elementLength+elementWidth) + elementLength*(e1-1)
else:
x1 = elementLength*e1
x2 = (e2 - 1) * elementWidth
x.append([x1, x2, x3])

# DERIVATIVES
Expand All @@ -466,8 +605,11 @@ def createArm(halfArmArcAngleRadians, elementLengths, armAngle, armAngleConst, e
ds2 = [dcent[0], -dcent[1], dcent[2]] if (e2 == 0) else dcent
ds2 = setMagnitude(ds2, dipMag)
ds1 = rotateAboutZAxis(ds2, -math.pi / 2)
elif e1 == elementsCount1 and e2 == elementsCount2-1:
elif e1 == elementsCount1 and e2 == elementsCount2-1: # armEnd
ds1 = [elementWidth,0,0]
ds2 = [0, 1 * (elementLength + elementWidth), 0]
if shorterArmEnd:
ds2 = [0, 0.5 * (elementLength + elementWidth), 0]
xnodes_ds1.append(ds1)
xnodes_ds2.append(ds2)
nid += 1
Expand Down Expand Up @@ -497,4 +639,5 @@ def createArm(halfArmArcAngleRadians, elementLengths, armAngle, armAngleConst, e
xnodes_ds1[j] = rotateAboutZAxis(xnodes_ds1[j], armAngle)
xnodes_ds2[j] = rotateAboutZAxis(xnodes_ds2[j], armAngle)

return (x, xnodes_ds1, xnodes_ds2, rmVertexNodes)
return (x, xnodes_ds1, xnodes_ds2, rmVertexNodes)

0 comments on commit c0ecb19

Please sign in to comment.