From 83a8d4dbdd9742beda544de6452db112d6fdf81e Mon Sep 17 00:00:00 2001 From: Richard Christie Date: Tue, 21 Jul 2020 16:10:46 +1200 Subject: [PATCH 1/2] Add start of tube bifurcation tree scaffold --- .../meshtypes/meshtype_1d_bifurcationtree1.py | 10 +- .../meshtype_2d_tubebifurcationtree1.py | 149 ++++++++++++++++++ src/scaffoldmaker/scaffolds.py | 2 + src/scaffoldmaker/utils/bifurcation.py | 43 ++++- 4 files changed, 201 insertions(+), 3 deletions(-) create mode 100644 src/scaffoldmaker/meshtypes/meshtype_2d_tubebifurcationtree1.py diff --git a/src/scaffoldmaker/meshtypes/meshtype_1d_bifurcationtree1.py b/src/scaffoldmaker/meshtypes/meshtype_1d_bifurcationtree1.py index 58e6c911..d502dd7c 100644 --- a/src/scaffoldmaker/meshtypes/meshtype_1d_bifurcationtree1.py +++ b/src/scaffoldmaker/meshtypes/meshtype_1d_bifurcationtree1.py @@ -134,6 +134,14 @@ def addChild(self, childTreeNode, d1=None, r=None): self._r.append(r) self._children.append(childTreeNode) + def getChild(self, childIndex): + ''' + :param childIndex: Index starting at 0 into child tree node list. + Return x1, d1, r1, x2, d2, r2 + ''' + assert 0 <= childIndex < len(self._children) + return self._children[childIndex] + def getChildCurve(self, childIndex): ''' :param childIndex: Index starting at 0 into child tree node list. @@ -223,7 +231,7 @@ def generateZincModel(self, region, nextNodeIdentifier=1, nextElementIdentifier= self._radius = findOrCreateFieldFiniteElement(self._fieldmodule, "radius", components_count=1, managed=True) self._fieldcache = self._fieldmodule.createFieldcache() parentNode = None - self._generateZincModelTree(self._rootNode, parentNode, nextNodeIdentifier, nextElementIdentifier) + nextNodeIdentifier, nextElementIdentifier = self._generateZincModelTree(self._rootNode, parentNode, nextNodeIdentifier, nextElementIdentifier) return nextNodeIdentifier, nextElementIdentifier def _getZincNodetemplate(self, d1VersionsCount, rVersionsCount): diff --git a/src/scaffoldmaker/meshtypes/meshtype_2d_tubebifurcationtree1.py b/src/scaffoldmaker/meshtypes/meshtype_2d_tubebifurcationtree1.py new file mode 100644 index 00000000..7692b131 --- /dev/null +++ b/src/scaffoldmaker/meshtypes/meshtype_2d_tubebifurcationtree1.py @@ -0,0 +1,149 @@ +""" +Generates a 2-D tube bifurcation mesh. +""" + +from __future__ import division +import math +from opencmiss.utils.maths.vectorops import cross, mult, normalize, sub +from opencmiss.utils.zinc.field import findOrCreateFieldCoordinates +from opencmiss.zinc.field import Field +from opencmiss.zinc.node import Node +from scaffoldmaker.meshtypes.meshtype_1d_bifurcationtree1 import MeshType_1d_bifurcationtree1 +from scaffoldmaker.meshtypes.scaffold_base import Scaffold_base +from scaffoldmaker.scaffoldpackage import ScaffoldPackage +from scaffoldmaker.utils.bifurcation import get_curve_circle_points, \ + make_tube_bifurcation_points, make_tube_bifurcation_elements_2d +from scaffoldmaker.utils.geometry import createCirclePoints +from scaffoldmaker.utils.interpolation import getCubicHermiteArcLength + + +class MeshType_2d_tubebifurcationtree1(Scaffold_base): + ''' + Generates a 2-D tube bifurcation tree mesh. + ''' + + @staticmethod + def getName(): + return '2D Tube Bifurcation Tree 1' + + @staticmethod + def getParameterSetNames(): + return ['Default'] + + @staticmethod + def getDefaultOptions(parameterSetName='Default'): + options = {} + options['Bifurcation tree'] = ScaffoldPackage(MeshType_1d_bifurcationtree1) + options['Number of elements around root'] = 6 + options['Maximum element length'] = 0.5 + return options + + @staticmethod + def getOrderedOptionNames(): + return [ + 'Bifurcation tree', + 'Number of elements around root', + 'Maximum element length'] + + @classmethod + def getOptionValidScaffoldTypes(cls, optionName): + if optionName == 'Bifurcation tree': + return [ MeshType_1d_bifurcationtree1 ] + return [] + + @classmethod + def getOptionScaffoldTypeParameterSetNames(cls, optionName, scaffoldType): + assert scaffoldType in cls.getOptionValidScaffoldTypes(optionName), cls.__name__ + '.getOptionScaffoldTypeParameterSetNames. ' + \ + 'Invalid option \'' + optionName + '\' scaffold type ' + scaffoldType.getName() + return scaffoldType.getParameterSetNames() + + @classmethod + def getOptionScaffoldPackage(cls, optionName, scaffoldType, parameterSetName=None): + ''' + :param parameterSetName: Name of valid parameter set for option Scaffold, or None for default. + :return: ScaffoldPackage. + ''' + if parameterSetName: + assert parameterSetName in cls.getOptionScaffoldTypeParameterSetNames(optionName, scaffoldType), \ + 'Invalid parameter set ' + str(parameterSetName) + ' for scaffold ' + str(scaffoldType.getName()) + ' in option ' + str(optionName) + ' of scaffold ' + cls.getName() + if optionName == 'Bifurcation tree': + if not parameterSetName: + parameterSetName = MeshType_1d_bifurcationtree1.getParameterSetNames()[0] + return ScaffoldPackage(MeshType_1d_bifurcationtree1, defaultParameterSetName=parameterSetName) + assert False, cls.__name__ + '.getOptionScaffoldPackage: Option ' + optionName + ' is not a scaffold' + + @staticmethod + def checkOptions(options): + ''' + :return: True if dependent options changed, otherwise False. This + happens where two or more options must change together to be valid. + ''' + dependentChanges = False + for key in [ + 'Number of elements around root']: + if options[key] < 4: + options[key] = 4 + if options['Maximum element length'] <= 0.0: + options['Maximum element length'] = 1.0 + return dependentChanges + + @classmethod + def generateBaseMesh(cls, region, options): + """ + Generate the base bicubic Hermite mesh. + :param region: Zinc region to define model in. Must be empty. + :param options: Dict containing options. See getDefaultOptions(). + :return: list of AnnotationGroup + """ + bifurcationTreeScaffold = options['Bifurcation tree'] + maxElementLength = options['Maximum element length'] + elementsCountAroundRoot = options['Number of elements around root'] + useCrossDerivatives = False + + fm = region.getFieldmodule() + coordinates = findOrCreateFieldCoordinates(fm) + cache = fm.createFieldcache() + + ############## + # Create nodes + ############## + + bifurcationTree = bifurcationTreeScaffold.getScaffoldType().generateBifurcationTree(bifurcationTreeScaffold.getScaffoldSettings()) + nodeIdentifier = bifurcationTree.generateZincModel(region)[0] + + nodes = fm.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) + nodetemplate = nodes.createNodetemplate() + nodetemplate.defineField(coordinates) + nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_VALUE, 1) + nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS1, 1) + nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS2, 1) + + rootNode = bifurcationTree.getRootNode() + child1 = rootNode.getChild(0) + + x1, xd1, r1, x2, xd2, r2 = child1.getChildCurve(0) + rd1 = rd2 = (r2 - r1) + + curveLength = getCubicHermiteArcLength(x1, xd1, x2, xd2) + elementsCount = max(2, math.ceil(curveLength/maxElementLength)) + elementLength = curveLength/elementsCount + side = [ 1.0, 0.0, 0.0 ] + for i in range(elementsCount + 1): + xi = i/elementsCount + x, d1, d2 = get_curve_circle_points(x1, xd1, x2, xd2, r1, rd1, r2, rd2, xi, elementLength, side, elementsCountAroundRoot) + for n in range(elementsCountAroundRoot): + #print('node',nodeIdentifier,i,n,x[n],d1[n],d2[n]) + node = nodes.createNode(nodeIdentifier, nodetemplate) + cache.setNode(node) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, x [n]) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, d1[n]) + coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, d2[n]) + nodeIdentifier = nodeIdentifier + 1 + + ################# + # Create elements + ################# + + elementIdentifier = 1 + + return [] diff --git a/src/scaffoldmaker/scaffolds.py b/src/scaffoldmaker/scaffolds.py index bcc0443d..ef1b3c0c 100644 --- a/src/scaffoldmaker/scaffolds.py +++ b/src/scaffoldmaker/scaffolds.py @@ -10,6 +10,7 @@ from scaffoldmaker.meshtypes.meshtype_2d_sphere1 import MeshType_2d_sphere1 from scaffoldmaker.meshtypes.meshtype_2d_tube1 import MeshType_2d_tube1 from scaffoldmaker.meshtypes.meshtype_2d_tubebifurcation1 import MeshType_2d_tubebifurcation1 +from scaffoldmaker.meshtypes.meshtype_2d_tubebifurcationtree1 import MeshType_2d_tubebifurcationtree1 from scaffoldmaker.meshtypes.meshtype_3d_bladder1 import MeshType_3d_bladder1 from scaffoldmaker.meshtypes.meshtype_3d_box1 import MeshType_3d_box1 from scaffoldmaker.meshtypes.meshtype_3d_boxhole1 import MeshType_3d_boxhole1 @@ -48,6 +49,7 @@ def __init__(self): MeshType_2d_sphere1, MeshType_2d_tube1, MeshType_2d_tubebifurcation1, + MeshType_2d_tubebifurcationtree1, MeshType_3d_bladder1, MeshType_3d_box1, MeshType_3d_boxhole1, diff --git a/src/scaffoldmaker/utils/bifurcation.py b/src/scaffoldmaker/utils/bifurcation.py index 581af52e..a1199f3f 100644 --- a/src/scaffoldmaker/utils/bifurcation.py +++ b/src/scaffoldmaker/utils/bifurcation.py @@ -3,15 +3,54 @@ """ from __future__ import division -from opencmiss.utils.maths.vectorops import add, cross, mult, normalize, sub +from opencmiss.utils.maths.vectorops import add, cross, dot, magnitude, mult, normalize, sub from opencmiss.zinc.element import Element, Elementbasis from opencmiss.zinc.field import Field from opencmiss.zinc.node import Node from scaffoldmaker.utils.eft_utils import remapEftNodeValueLabel, scaleEftNodeValueLabels, setEftScaleFactorIds +from scaffoldmaker.utils.geometry import createCirclePoints from scaffoldmaker.utils.interpolation import DerivativeScalingMode, interpolateCubicHermite, interpolateCubicHermiteDerivative, \ - smoothCubicHermiteDerivativesLine, interpolateLagrangeHermiteDerivative + interpolateCubicHermiteSecondDerivative, smoothCubicHermiteDerivativesLine, interpolateLagrangeHermiteDerivative +def get_curve_circle_points(x1, xd1, x2, xd2, r1, rd1, r2, rd2, xi, dmag, side, elementsCountAround): + ''' + :param dmag: Magnitude of derivative on curve. + :param side: Vector in side direction of first node around. + Need not be unit or exactly normal to curve at xi. + :return: x[], d1[] around, d2[] along + ''' + cx = interpolateCubicHermite(x1, xd1, x2, xd2, xi) + cxd = interpolateCubicHermiteDerivative(x1, xd1, x2, xd2, xi) + mag_cxd = magnitude(cxd) + cxd2 = interpolateCubicHermiteSecondDerivative(x1, xd1, x2, xd2, xi) + mag_cxd2 = magnitude(cxd2) + r = interpolateCubicHermite([ r1 ], [ rd1 ], [ r2 ], [ rd2 ], xi)[0] + rd = interpolateCubicHermiteDerivative([ r1 ], [ rd1 ], [ r2 ], [ rd2 ], xi)[0] + axis1 = normalize(cxd) + axis3 = normalize(cross(axis1, side)) + axis2 = cross(axis3, axis1) + x, d1 = createCirclePoints(cx, mult(axis2, r), mult(axis3, r), elementsCountAround) + curvatureVector = mult(cxd2, 1.0/(mag_cxd*mag_cxd)) + d2 = [] + radialGrowth = rd/(mag_cxd*r) + for e in range(elementsCountAround): + radialVector = sub(x[e], cx) + dmagFinal = dmag*(1.0 - dot(radialVector, curvatureVector)) + # add curvature and radius change components: + d2.append(add(mult(cxd, dmagFinal/mag_cxd), mult(radialVector, dmagFinal*radialGrowth))) + return x, d1, d2 + + +def track_curve_side_axis(x1, d1, x2, d2, sideStart, xiStart, xiEnd): + ''' + Get side vector normal to curve at xiEnd for smoothest transition from + sideStart at xiStart. + :param xi, d1, x2, d2: + :param sideStart: Unit vector normal to curve + ''' + pass + def get_bifurcation_triple_point(p1x, p1d, p2x, p2d, p3x, p3d): ''' Get coordinates and derivatives of triple point between p1, p2 and p3 with derivatives. From 4d2222fc2be76759057f441cb5c8ff820740f7ab Mon Sep 17 00:00:00 2001 From: Richard Christie Date: Fri, 12 Feb 2021 10:29:17 +1300 Subject: [PATCH 2/2] Disable bifurcation tree 2D scaffold --- src/scaffoldmaker/scaffolds.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/scaffoldmaker/scaffolds.py b/src/scaffoldmaker/scaffolds.py index 24513ba7..0993d621 100644 --- a/src/scaffoldmaker/scaffolds.py +++ b/src/scaffoldmaker/scaffolds.py @@ -55,7 +55,7 @@ def __init__(self): MeshType_2d_sphere1, MeshType_2d_tube1, MeshType_2d_tubebifurcation1, - MeshType_2d_tubebifurcationtree1, + #MeshType_2d_tubebifurcationtree1, MeshType_3d_bladder1, MeshType_3d_bladderurethra1, MeshType_3d_box1,