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

Network tube #240

Merged
merged 42 commits into from
Nov 23, 2023
Merged
Changes from all commits
Commits
Show all changes
42 commits
Select commit Hold shift + click to select a range
d274d50
Add tube bifurcation scaffold and utilities
rchristie Aug 6, 2023
5561a12
Merge remote-tracking branch 'abi/main' into network_tube
rchristie Aug 20, 2023
46f4865
Get intersection point between 2 surfaces
rchristie Aug 23, 2023
ba65c96
Merge branch 'main' into network_tube
rchristie Aug 27, 2023
035d1ea
Add TrackSurface method to get intersection curve
rchristie Aug 28, 2023
761a431
Improve TrackSurface get intersection curve
rchristie Aug 30, 2023
5c543e0
Compute curve nearest and intersection points
rchristie Sep 4, 2023
318025a
Find nearest/intersection of curve and TrackSurface
rchristie Sep 4, 2023
82c987d
Sample curves between start, end locations
rchristie Sep 5, 2023
b0ac138
Limit tube coordinates to start, end surface intersections
rchristie Sep 5, 2023
44318ec
Resample path tube coordinates between surfaces
rchristie Sep 8, 2023
f907c0d
Fix cross derivatives in tube network
rchristie Sep 8, 2023
6cca7cd
Use target element aspect ratio
rchristie Sep 8, 2023
ee1a47b
Improve increments on track surfaces
rchristie Sep 13, 2023
d4bef3d
Fix surface boundary intersections
rchristie Oct 10, 2023
3653c89
Fix getting intersection curve between aligned surfaces
rchristie Oct 10, 2023
10ad980
Fix surface intersection curve boundaries
rchristie Oct 12, 2023
bb75fa1
Improve tube intersection tests
rchristie Oct 12, 2023
ec10965
Output surfaces and curves to groups
rchristie Oct 12, 2023
f344568
Improve surface intersections
rchristie Oct 16, 2023
8f6b220
Pass all surface intersection tests
rchristie Oct 16, 2023
dbfc3df
Handle no progress on surface intersection curves
rchristie Oct 17, 2023
845e409
Get trim surfaces working
rchristie Oct 17, 2023
db0523e
Create bifurcation nodes and elements
rchristie Oct 26, 2023
2fcbc25
Improve tube bifurcations
rchristie Nov 5, 2023
76790fb
Find surface intersection point in more cases
rchristie Nov 5, 2023
e82b013
Improve curve intersection at boundaries
rchristie Nov 6, 2023
0f488a6
Minor intersection curve boundary improvements
rchristie Nov 7, 2023
0ea1eae
Fix converging bifurcation
rchristie Nov 7, 2023
9e3b91c
Fix more bifurcation cases
rchristie Nov 8, 2023
2e58917
Add sphere cube network layout
rchristie Nov 9, 2023
0956685
Transfer annotation groups from layout to tube network
rchristie Nov 10, 2023
f4fa495
Add inner coordinates to network layout
rchristie Nov 14, 2023
590d118
Add 3D tube network
rchristie Nov 20, 2023
5872124
Merge remote-tracking branch 'abi/main' into network_tube
rchristie Nov 20, 2023
57e2162
Increment version number
rchristie Nov 20, 2023
bfb8149
Use common nodes on bifurcation-bifurcation connections
rchristie Nov 21, 2023
8721d5b
Fix assign inner coordinates for selected elements
rchristie Nov 21, 2023
13b8285
Fix tube-tube segment connections
rchristie Nov 21, 2023
5e2ba2e
Merge scale/offset into assign coordinates
rchristie Nov 22, 2023
b2d35d6
Review fixes, pep8 fixes
rchristie Nov 22, 2023
e4638f4
Fix uterus calls to bifurcation function
rchristie Nov 23, 2023
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
2 changes: 1 addition & 1 deletion setup.py
Original file line number Diff line number Diff line change
@@ -35,7 +35,7 @@ def readfile(filename, split=False):

setup(
name="scaffoldmaker",
version="0.11.0",
version="0.12.0",
description="Python client for generating anatomical scaffolds using Zinc",
long_description="\n".join(readme) + source_license,
long_description_content_type="text/x-rst",
316 changes: 298 additions & 18 deletions src/scaffoldmaker/meshtypes/meshtype_1d_network_layout1.py

Large diffs are not rendered by default.

10 changes: 6 additions & 4 deletions src/scaffoldmaker/meshtypes/meshtype_1d_path1.py
Original file line number Diff line number Diff line change
@@ -131,7 +131,7 @@ def generateBaseMesh(cls, region, options):
return [], None

@classmethod
def makeSideDerivativesNormal(cls, region, options, functionOptions, editGroupName):
def makeSideDerivativesNormal(cls, region, options, constructionObject, functionOptions, editGroupName):
fieldmodule = region.getFieldmodule()
coordinates = fieldmodule.findFieldByName("coordinates").castFiniteElement()
nodeset = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES)
@@ -143,7 +143,7 @@ def makeSideDerivativesNormal(cls, region, options, functionOptions, editGroupNa
return False, True # settings not changed, nodes changed

@classmethod
def smoothSideCrossDerivatives(cls, region, options, functionOptions, editGroupName):
def smoothSideCrossDerivatives(cls, region, options, constructionObject, functionOptions, editGroupName):
smoothD12 = options['D2 derivatives'] and functionOptions['Smooth D12']
smoothD13 = options['D3 derivatives'] and functionOptions['Smooth D13']
if not (smoothD12 or smoothD13):
@@ -182,9 +182,11 @@ def getInteractiveFunctions(cls):
("Make side derivatives normal...",
{ 'Make D2 normal': True,
'Make D3 normal': True },
lambda region, options, functionOptions, editGroupName: cls.makeSideDerivativesNormal(region, options, functionOptions, editGroupName)),
lambda region, options, constructionObject, functionOptions, editGroupName:
cls.makeSideDerivativesNormal(region, options, constructionObject, functionOptions, editGroupName)),
("Smooth side cross derivatives...",
{ 'Smooth D12' : True,
'Smooth D13' : True },
lambda region, options, functionOptions, editGroupName: cls.smoothSideCrossDerivatives(region, options, functionOptions, editGroupName))
lambda region, options, constructionObject, functionOptions, editGroupName:
cls.smoothSideCrossDerivatives(region, options, constructionObject, functionOptions, editGroupName))
]
4 changes: 2 additions & 2 deletions src/scaffoldmaker/meshtypes/meshtype_2d_tubebifurcation1.py
Original file line number Diff line number Diff line change
@@ -84,7 +84,7 @@ def checkOptions(options):
c2Count = options['Number of elements around child 2']
if (paCount + c1Count + c2Count) % 2:
c2Count += 1
pac1Count, pac2Count, c1c2Count = get_tube_bifurcation_connection_elements_counts(paCount, c1Count, c2Count)
pac1Count, c1c2Count, pac2Count = get_tube_bifurcation_connection_elements_counts([paCount, c1Count, c2Count])
if pac1Count < 2:
c2Count -= 2*(2 - pac1Count)
elif pac2Count < 2:
@@ -126,7 +126,7 @@ def generateBaseMesh(cls, region, options):
c2Centre = [ child2Length*math.sin(child2AngleRadians), 0.0, child2Length*math.cos(child2AngleRadians) ]
c12 = sub(c1Centre, c2Centre)

pac1Count, pac2Count, c1c2Count = get_tube_bifurcation_connection_elements_counts(paCount, c1Count, c2Count)
pac1Count, c1c2Count, pac2Count = get_tube_bifurcation_connection_elements_counts([paCount, c1Count, c2Count])

# parent ring
paAxis3 = [ 0.0, 0.0, parentLength ]
122 changes: 122 additions & 0 deletions src/scaffoldmaker/meshtypes/meshtype_2d_tubenetwork1.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,122 @@
"""
Generates a 2-D Hermite bifurcating tube network.
"""

import copy

from cmlibs.utils.zinc.field import find_or_create_field_coordinates
from cmlibs.utils.zinc.finiteelement import get_maximum_element_identifier, get_maximum_node_identifier
from cmlibs.zinc.field import Field
from scaffoldmaker.meshtypes.meshtype_1d_network_layout1 import MeshType_1d_network_layout1
from scaffoldmaker.meshtypes.scaffold_base import Scaffold_base
from scaffoldmaker.scaffoldpackage import ScaffoldPackage
from scaffoldmaker.utils.bifurcation import generateTubeBifurcationTree


class MeshType_2d_tubenetwork1(Scaffold_base):
"""
Generates a 2-D hermite bifurcating tube network.
"""

@staticmethod
def getName():
return "2D Tube Network 1"

@classmethod
def getParameterSetNames(cls):
return MeshType_1d_network_layout1.getParameterSetNames()

@classmethod
def getDefaultOptions(cls, parameterSetName="Default"):
options = {
"Network layout": ScaffoldPackage(MeshType_1d_network_layout1, defaultParameterSetName=parameterSetName),
"Elements count around": 8,
"Target element aspect ratio": 2.0,
"Serendipity": True
}
return options

@staticmethod
def getOrderedOptionNames():
return [
"Network layout",
"Elements count around",
"Target element aspect ratio",
"Serendipity"
]

@classmethod
def getOptionValidScaffoldTypes(cls, optionName):
if optionName == "Network layout":
return [MeshType_1d_network_layout1]
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() # use the defaults from the network layout scaffold

@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 == "Network layout":
if not parameterSetName:
parameterSetName = "Default"
return ScaffoldPackage(scaffoldType, defaultParameterSetName=parameterSetName)
assert False, cls.__name__ + ".getOptionScaffoldPackage: Option " + optionName + " is not a scaffold"

@classmethod
def checkOptions(cls, options):
if not options["Network layout"].getScaffoldType() in cls.getOptionValidScaffoldTypes("Network layout"):
options["Network layout"] = cls.getOptionScaffoldPackage("Network layout")
elementsCountAround = options["Elements count around"]
options["Elements count around"] = max(4, elementsCountAround + (elementsCountAround % 2))
if options["Target element aspect ratio"] < 0.01:
options["Target element aspect ratio"] = 0.01
dependentChanges = False
return dependentChanges

@staticmethod
def generateBaseMesh(region, options):
"""
Generate the base hermite-bilinear mesh. See also generateMesh().
:param region: Zinc region to define model in. Must be empty.
:param options: Dict containing options. See getDefaultOptions().
:return: list of AnnotationGroup, None
"""
elementsCountAround = options["Elements count around"]
networkLayout = options["Network layout"]
targetElementAspectRatio = options["Target element aspect ratio"]
serendipity = options["Serendipity"]

layoutRegion = region.createRegion()
networkLayout.generate(layoutRegion) # ask scaffold to generate to get user-edited parameters
layoutAnnotationGroups = networkLayout.getAnnotationGroups()

fieldmodule = region.getFieldmodule()
nodes = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES)
coordinates = find_or_create_field_coordinates(fieldmodule)
nodeIdentifier = max(get_maximum_node_identifier(nodes), 0) + 1
mesh = fieldmodule.findMeshByDimension(2)
elementIdentifier = max(get_maximum_element_identifier(mesh), 0) + 1

networkMesh = networkLayout.getConstructionObject()

try:
nodeIdentifier, elementIdentifier, annotationGroups = generateTubeBifurcationTree(
networkMesh, region, coordinates, nodeIdentifier, elementIdentifier,
elementsCountAround, targetElementAspectRatio, 1, layoutAnnotationGroups, serendipity=serendipity)
except Exception as e:
print(e, "\nException occurred while generating tube network: Please edit network layout")
return [], None

return annotationGroups, None
8 changes: 4 additions & 4 deletions src/scaffoldmaker/meshtypes/meshtype_3d_bladderurethra1.py
Original file line number Diff line number Diff line change
@@ -1480,12 +1480,12 @@ def generateUreterInlets(region, nodes, mesh, ureterDefaultOptions, elementsCoun
startProportions2.append(trackSurfaceUreter2.getProportion(o2_Positions[n]))

endProportions1 = []
elementsAroundTrackSurface1 = trackSurfaceUreter1.elementsCount1
elementsAlongTrackSurface1 = trackSurfaceUreter1.elementsCount2
elementsAroundTrackSurface1 = trackSurfaceUreter1._elementsCount1
elementsAlongTrackSurface1 = trackSurfaceUreter1._elementsCount2

endProportions2 = []
elementsAroundTrackSurface2 = trackSurfaceUreter2.elementsCount1
elementsAlongTrackSurface2 = trackSurfaceUreter2.elementsCount2
elementsAroundTrackSurface2 = trackSurfaceUreter2._elementsCount1
elementsAlongTrackSurface2 = trackSurfaceUreter2._elementsCount2

firstIdxAround1 = ureterElementPositionAround + (0 if ureter1Position.xi1 > 0.5 else -1)
firstIdxAlong = ureterElementPositionDown - (0 if ureter1Position.xi2 > 0.5 else 1)
4 changes: 2 additions & 2 deletions src/scaffoldmaker/meshtypes/meshtype_3d_cecum1.py
Original file line number Diff line number Diff line change
@@ -956,8 +956,8 @@ def getElementIdxOfOstiumBoundary(centrePosition, trackSurfaceOstium, ostiumDiam
:return: element indices on the left, right, bottom and top boundaries around tracksurface.
"""

elementsAroundTrackSurface = trackSurfaceOstium.elementsCount1
elementsAlongTrackSurface = trackSurfaceOstium.elementsCount2
elementsAroundTrackSurface = trackSurfaceOstium._elementsCount1
elementsAlongTrackSurface = trackSurfaceOstium._elementsCount2
ei1 = centrePosition.e1
ei2 = centrePosition.e2
xi1 = centrePosition.xi1
Loading