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.

Sign up for GitHub

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

Annotation #92

Merged
merged 13 commits into from
Sep 24, 2020
Merged

Annotation #92

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
13 commits
Select commit Hold shift + click to select a range
fede36b
Adding markers and parameter set for species-specific datasets (mouse)
s-fong Sep 10, 2020
d68d9fc
Merge remote-tracking branch 'main/master' into annotation
s-fong Sep 10, 2020
07cbe01
With face groups for fitting
s-fong Sep 14, 2020
796f692
Tidy up
s-fong Sep 17, 2020
d00c971
Merge remote-tracking branch 'main/master' into annotation
s-fong Sep 17, 2020
42ee4ae
Update stellate_terms
s-fong Sep 17, 2020
ee71a0d
Add soma markers
s-fong Sep 21, 2020
545d30a
With marker positions as proportion of arms, and calculated outside e…
s-fong Sep 22, 2020
d5d709b
Remove Mouse 2
s-fong Sep 23, 2020
6ef7604
Style fixes.
s-fong Sep 23, 2020
632c038
Specify 'stellate' or 'cervicothoracic ganglion' depending on choice …
s-fong Sep 23, 2020
8df9408
Rename wholeScaffold to stellate. Improve efficiency.
s-fong Sep 24, 2020
ee7db4f
Remove unecessary lines.
s-fong Sep 24, 2020
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
19 changes: 19 additions & 0 deletions src/scaffoldmaker/annotation/stellate_terms.py
View file
Open in desktop
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
Open in desktop
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)