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

Implement phase #68

Merged
merged 7 commits into from
Feb 21, 2020
Merged

Implement phase #68

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
1bda1fb
Remove redundant parameters and arguments
mlin865 Jan 2, 2020
9f61c69
Add colonsection and allow phase selection to start along elements
mlin865 Jan 8, 2020
72aef74
Allow variation in start phase
mlin865 Jan 15, 2020
45bc1ef
Merge remote-tracking branch 'abi/master' into phase
mlin865 Jan 15, 2020
9518818
Add assertion to ensure start phase is zero for varying radius and tc…
mlin865 Jan 19, 2020
5bdb1b0
Implement phase for colon segment with varying radius
mlin865 Feb 13, 2020
de3c728
Merge remote-tracking branch 'abi/master' into phase
mlin865 Feb 13, 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
97 changes: 54 additions & 43 deletions src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,6 @@
from scaffoldmaker.utils.meshrefinement import MeshRefinement
from scaffoldmaker.utils import interpolation as interp
from scaffoldmaker.utils import tubemesh
from scaffoldmaker.utils import vector
from scaffoldmaker.utils.zinc_utils import exnodeStringFromNodeValues
from opencmiss.zinc.node import Node

Expand Down Expand Up @@ -138,7 +137,20 @@ class MeshType_3d_colon1(Scaffold_base):
[ [ 32.0, -49.3, 24.4 ], [ -51.1, -50.1, -1.2 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ],
[ [ -51.1, -45.3, 15.7 ], [ -38.0, 52.2, -16.2 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ],
[ [ -122.9, 124.2, -36.0 ], [ -21.3, 64.5, -18.3 ], [ 0.0, 0.0, 5.0 ], [ 0.0, 0.0, 0.5 ] ] ] )
} )
} ),
'Pig 2': ScaffoldPackage(MeshType_1d_path1, {
'scaffoldSettings': {
'Coordinate dimensions': 3,
'Length': 90.0,
'Number of elements': 3
},
'meshEdits': exnodeStringFromNodeValues(
[Node.VALUE_LABEL_VALUE, Node.VALUE_LABEL_D_DS1, Node.VALUE_LABEL_D_DS2, Node.VALUE_LABEL_D2_DS1DS2], [
[ [ 0.0, 0.0, 0.0 ], [ 30.0, 0.0, 0.0 ], [ 0.0, 1.0, 0.0 ], [ 0.0, 0.0, 0.0 ] ],
[ [ 30.0, 0.0, 0.0 ], [ 30.0, 0.0, 0.0 ], [ 0.0, 1.0, 0.0 ], [ 0.0, 0.0, 0.0 ] ],
[ [ 60.0, 0.0, 0.0 ], [ 30.0, 0.0, 0.0 ], [ 0.0, 1.0, 0.0 ], [ 0.0, 0.0, 0.0 ] ],
[ [ 90.0, 0.0, 0.0 ], [ 30.0, 0.0, 0.0 ], [ 0.0, 1.0, 0.0 ], [ 0.0, 0.0, 0.0 ] ] ] )
} ),
}

@staticmethod
Expand All @@ -153,7 +165,8 @@ def getParameterSetNames():
'Human 2',
'Mouse 1',
'Mouse 2',
'Pig 1']
'Pig 1',
'Pig 2']

@classmethod
def getDefaultOptions(cls, parameterSetName='Default'):
Expand All @@ -163,8 +176,10 @@ def getDefaultOptions(cls, parameterSetName='Default'):
centralPathOption = cls.centralPathDefaultScaffoldPackages['Mouse 1']
elif 'Mouse 2' in parameterSetName:
centralPathOption = cls.centralPathDefaultScaffoldPackages['Mouse 2']
elif 'Pig' in parameterSetName:
elif 'Pig 1' in parameterSetName:
centralPathOption = cls.centralPathDefaultScaffoldPackages['Pig 1']
elif 'Pig 2' in parameterSetName:
centralPathOption = cls.centralPathDefaultScaffoldPackages['Pig 2']
else:
centralPathOption = cls.centralPathDefaultScaffoldPackages['Human 1']
if 'Mouse' in parameterSetName:
Expand All @@ -177,6 +192,7 @@ def getDefaultOptions(cls, parameterSetName='Default'):
'Central path' : copy.deepcopy(centralPathOption),
'Segment profile' : segmentProfileOption,
'Number of segments': 30,
'Start phase': 0.0,
'Proximal length': 420.0,
'Transverse length': 460.0,
'Distal length': 620.0,
Expand Down Expand Up @@ -210,7 +226,7 @@ def getDefaultOptions(cls, parameterSetName='Default'):
options['Transverse-distal tenia coli width'] = 1.0
options['Distal inner radius'] = 0.7
options['Distal tenia coli width'] = 1.0
elif 'Pig' in parameterSetName:
elif 'Pig 1' in parameterSetName:
options['Number of segments'] = 120
options['Proximal length'] = 2610.0
options['Transverse length'] = 200.0
Expand All @@ -223,6 +239,19 @@ def getDefaultOptions(cls, parameterSetName='Default'):
options['Transverse-distal tenia coli width'] = 5.0
options['Distal inner radius'] = 8.0
options['Distal tenia coli width'] = 5.0
elif 'Pig 2' in parameterSetName:
options['Number of segments'] = 3
options['Proximal length'] = 30.0
options['Transverse length'] = 30.0
options['Distal length'] = 30.0
options['Proximal inner radius'] = 16.0
options['Proximal tenia coli width'] = 5.0
options['Proximal-transverse inner radius'] = 16.0
options['Proximal-transverse tenia coli width'] = 5.0
options['Transverse-distal inner radius'] = 16.0
options['Transverse-distal tenia coli width'] = 5.0
options['Distal inner radius'] = 16.0
options['Distal tenia coli width'] = 5.0
return options

@staticmethod
Expand All @@ -231,6 +260,7 @@ def getOrderedOptionNames():
'Central path',
'Segment profile',
'Number of segments',
'Start phase',
'Proximal length',
'Transverse length',
'Distal length',
Expand Down Expand Up @@ -321,6 +351,7 @@ def generateBaseMesh(region, options):
centralPath = options['Central path']
segmentProfile = options['Segment profile']
segmentCount = options['Number of segments']
startPhase = options['Start phase'] % 360.0
proximalLength = options['Proximal length']
transverseLength = options['Transverse length']
distalLength = options['Distal length']
Expand All @@ -332,7 +363,6 @@ def generateBaseMesh(region, options):
transverseDistalTCWidth = options['Transverse-distal tenia coli width']
distalInnerRadius = options['Distal inner radius']
distalTCWidth = options['Distal tenia coli width']
segmentScaffoldType = segmentProfile.getScaffoldType()
segmentSettings = segmentProfile.getScaffoldSettings()

elementsCountAroundTC = segmentSettings['Number of elements around tenia coli']
Expand All @@ -347,10 +377,6 @@ def generateBaseMesh(region, options):

elementsCountAlongSegment = segmentSettings['Number of elements along segment']
elementsCountThroughWall = segmentSettings['Number of elements through wall']
startRadius = segmentSettings['Start inner radius']
startRadiusDerivative = segmentSettings['Start inner radius derivative']
endRadius = segmentSettings['End inner radius']
endRadiusDerivative = segmentSettings['End inner radius derivative']
wallThickness = segmentSettings['Wall thickness']
useCrossDerivatives = segmentSettings['Use cross derivatives']
useCubicHermiteThroughWall = not(segmentSettings['Use linear through wall'])
Expand All @@ -364,10 +390,7 @@ def generateBaseMesh(region, options):
centralPath.generate(tmpRegion)
cx, cd1, cd2, cd12 = extractPathParametersFromRegion(tmpRegion)
# for i in range(len(cx)):
# print('cx = ', i+1, cx[i])
# print('cd1 = ', i+1, cd1[i])
# print('cd2 = ', i+1, cd2[i])
# print('cd12 = ', i+1, cd12[i])
# print(i, '[', cx[i], ',', cd1[i], ',', cd2[i], ',', cd12[i], '],')
del tmpRegion

# find arclength of colon
Expand All @@ -388,11 +411,16 @@ def generateBaseMesh(region, options):

# Generate variation of radius & tc width along length
lengthList = [0.0, proximalLength, proximalLength + transverseLength, length]
innerRadiusList = [proximalInnerRadius, proximalTransverseInnerRadius, transverseDistalInnerRadius, distalInnerRadius]
innerRadiusSegmentList, dInnerRadiusSegmentList = interp.sampleParameterAlongLine(lengthList, innerRadiusList, segmentCount)
innerRadiusList = [proximalInnerRadius, proximalTransverseInnerRadius,
transverseDistalInnerRadius, distalInnerRadius]
innerRadiusAlongElementList, dInnerRadiusAlongElementList = interp.sampleParameterAlongLine(lengthList,
innerRadiusList,
elementsCountAlong)

tcWidthList = [proximalTCWidth, proximalTransverseTCWidth, transverseDistalTCWidth, distalTCWidth]
tcWidthSegmentList, dTCWidthSegmentList = interp.sampleParameterAlongLine(lengthList, tcWidthList, segmentCount)
tcWidthAlongElementList, dTCWidthAlongElementList = interp.sampleParameterAlongLine(lengthList,
tcWidthList,
elementsCountAlong)

xExtrude = []
d1Extrude = []
Expand All @@ -404,48 +432,31 @@ def generateBaseMesh(region, options):
region, elementsCountAroundTC, elementsCountAroundHaustrum, elementsCountAlongSegment,
tcCount, segmentLengthEndDerivativeFactor, segmentLengthMidDerivativeFactor,
segmentLength, wallThickness, cornerInnerRadiusFactor, haustrumInnerRadiusFactor,
innerRadiusSegmentList, dInnerRadiusSegmentList, tcWidthSegmentList, dTCWidthSegmentList)
innerRadiusAlongElementList, dInnerRadiusAlongElementList, tcWidthAlongElementList,
startPhase)

for nSegment in range(segmentCount):
# Create inner points
xInner, d1Inner, d2Inner, transitElementList, segmentAxis, annotationGroups, annotationArray = \
colonSegmentTubeMeshInnerPoints.getColonSegmentTubeMeshInnerPoints(nSegment)
xInner, d1Inner, d2Inner, transitElementList, segmentAxis, annotationGroups, annotationArray, \
faceMidPointsZ = colonSegmentTubeMeshInnerPoints.getColonSegmentTubeMeshInnerPoints(nSegment)

# Warp segment points
xWarpedList, d1WarpedList, d2WarpedList, d3WarpedUnitList = tubemesh.warpSegmentPoints(
xInner, d1Inner, d2Inner, segmentAxis, segmentLength, sx, sd1, sd2,
elementsCountAround, elementsCountAlongSegment, nSegment)
elementsCountAround, elementsCountAlongSegment, nSegment, faceMidPointsZ)

# Store points along length
xExtrude = xExtrude + (xWarpedList if nSegment == 0 else xWarpedList[elementsCountAround:])
d1Extrude = d1Extrude + (d1WarpedList if nSegment == 0 else d1WarpedList[elementsCountAround:])

# Smooth d2 for nodes between segments and recalculate d3
if nSegment == 0:
d2Extrude = d2Extrude + (d2WarpedList[:-elementsCountAround])
d3UnitExtrude = d3UnitExtrude + (d3WarpedUnitList[:-elementsCountAround])
else:
xSecondFace = xWarpedList[elementsCountAround:elementsCountAround*2]
d1SecondFace = d1WarpedList[elementsCountAround:elementsCountAround*2]
d2SecondFace = d2WarpedList[elementsCountAround:elementsCountAround*2]
for n1 in range(elementsCountAround):
nx = [xLastTwoFaces[n1], xLastTwoFaces[n1 + elementsCountAround], xSecondFace[n1]]
nd2 = [d2LastTwoFaces[n1], d2LastTwoFaces[n1 + elementsCountAround], d2SecondFace[n1]]
d2 = interp.smoothCubicHermiteDerivativesLine(nx, nd2, fixStartDerivative = True, fixEndDerivative = True)[1]
d2Extrude.append(d2)
d3Unit = vector.normalise(vector.crossproduct3(vector.normalise(d1LastTwoFaces[n1 + elementsCountAround]), vector.normalise(d2)))
d3UnitExtrude.append(d3Unit)
d2Extrude = d2Extrude + (d2WarpedList[elementsCountAround:-elementsCountAround] if nSegment < segmentCount - 1 else d2WarpedList[elementsCountAround:])
d3UnitExtrude = d3UnitExtrude + (d3WarpedUnitList[elementsCountAround:-elementsCountAround] if nSegment < segmentCount - 1 else d3WarpedUnitList[elementsCountAround:])
xLastTwoFaces = xWarpedList[-elementsCountAround*2:]
d1LastTwoFaces = d1WarpedList[-elementsCountAround*2:]
d2LastTwoFaces = d2WarpedList[-elementsCountAround*2:]
d2Extrude = d2Extrude + (d2WarpedList if nSegment == 0 else d2WarpedList[elementsCountAround:])
d3UnitExtrude = d3UnitExtrude + (
d3WarpedUnitList if nSegment == 0 else d3WarpedUnitList[elementsCountAround:])

contractedWallThicknessList = colonSegmentTubeMeshInnerPoints.getContractedWallThicknessList()

# Create coordinates and derivatives
xList, d1List, d2List, d3List, curvatureList = tubemesh.getCoordinatesFromInner(xExtrude, d1Extrude,
d2Extrude, d3UnitExtrude, sx, contractedWallThicknessList,
d2Extrude, d3UnitExtrude, contractedWallThicknessList,
elementsCountAround, elementsCountAlong, elementsCountThroughWall, transitElementList)

relaxedLengthList, xiList = colonSegmentTubeMeshInnerPoints.getRelaxedLengthAndXiList()
Expand Down
Loading