Merge pull request #263 from rchristie/sample-tube

Improve tube mesh core sampling

src/scaffoldmaker/meshtypes/meshtype_1d_network_layout1.py

@@ -186,8 +186,8 @@ def _defineInnerCoordinates(cls, region, coordinates, options, networkMesh):
             "To field": {"coordinates": False, "inner coordinates": True},
             "From field": {"coordinates": True, "inner coordinates": False},
             "Mode": {"Scale": True, "Offset": False},
-            "D2 value": 0.5,
-            "D3 value": 0.5}
+            "D2 value": 0.8,
+            "D3 value": 0.8}
         cls.assignCoordinates(region, options, networkMesh, functionOptions, editGroupName=None)
 
     @classmethod

src/scaffoldmaker/meshtypes/meshtype_2d_tubenetwork1.py

@@ -24,8 +24,8 @@ def getParameterSetNames(cls):
     def getDefaultOptions(cls, parameterSetName="Default"):
         options = {
             "Network layout": ScaffoldPackage(MeshType_1d_network_layout1, defaultParameterSetName=parameterSetName),
-            "Elements count around": 8,
-            "Annotation elements counts around": [0],
+            "Number of elements around": 8,
+            "Annotation numbers of elements around": [0],
             "Target element density along longest segment": 4.0,
             "Show trim surfaces": False
         }
@@ -35,8 +35,8 @@ def getDefaultOptions(cls, parameterSetName="Default"):
     def getOrderedOptionNames():
         return [
             "Network layout",
-            "Elements count around",
-            "Annotation elements counts around",
+            "Number of elements around",
+            "Annotation numbers of elements around",
             "Target element density along longest segment",
             "Show trim surfaces"
         ]
@@ -74,12 +74,12 @@ def getOptionScaffoldPackage(cls, optionName, scaffoldType, parameterSetName=Non
     def checkOptions(cls, options):
         if not options["Network layout"].getScaffoldType() in cls.getOptionValidScaffoldTypes("Network layout"):
             options["Network layout"] = cls.getOptionScaffoldPackage("Network layout", MeshType_1d_network_layout1)
-        elementsCountAround = options["Elements count around"]
-        if options["Elements count around"] < 4:
-            options["Elements count around"] = 4
-        annotationElementsCountsAround = options["Annotation elements counts around"]
+        dependentChanges = False
+        if options["Number of elements around"] < 4:
+            options["Number of elements around"] = 4
+        annotationElementsCountsAround = options["Annotation numbers of elements around"]
         if len(annotationElementsCountsAround) == 0:
-            options["Annotation elements count around"] = [0]
+            options["Annotation numbers of elements around"] = [0]
         else:
             for i in range(len(annotationElementsCountsAround)):
                 if annotationElementsCountsAround[i] <= 0:
@@ -88,7 +88,6 @@ def checkOptions(cls, options):
                     annotationElementsCountsAround[i] = 4
         if options["Target element density along longest segment"] < 1.0:
             options["Target element density along longest segment"] = 1.0
-        dependentChanges = False
         return dependentChanges
 
     @classmethod
@@ -107,10 +106,10 @@ def generateBaseMesh(cls, region, options):
         tubeNetworkMeshBuilder = TubeNetworkMeshBuilder(
             networkMesh,
             targetElementDensityAlongLongestSegment=options["Target element density along longest segment"],
-            defaultElementsCountAround=options["Elements count around"],
+            defaultElementsCountAround=options["Number of elements around"],
             elementsCountThroughWall=1,
             layoutAnnotationGroups=networkLayout.getAnnotationGroups(),
-            annotationElementsCountsAround=options["Annotation elements counts around"])
+            annotationElementsCountsAround=options["Annotation numbers of elements around"])
         tubeNetworkMeshBuilder.build()
         generateData = TubeNetworkMeshGenerateData(
             region, 2,

src/scaffoldmaker/meshtypes/meshtype_3d_wholebody2.py

@@ -198,11 +198,11 @@ def getDefaultNetworkLayoutScaffoldPackage(cls, parameterSetName):
 
 class MeshType_3d_wholebody2(Scaffold_base):
     """
-    Generates a 3-D hermite bifurcating tube network, with linear basis through wall.
+    Generates a 3-D hermite bifurcating tube network with core representing the human body.
     """
 
-    @staticmethod
-    def getName():
+    @classmethod
+    def getName(cls):
         return "3D Whole Body 2"
 
     @classmethod
@@ -225,38 +225,29 @@ def getDefaultOptions(cls, parameterSetName="Default"):
             "Number of elements around torso": 12,
             "Number of elements around arm": 8,
             "Number of elements around leg": 8,
-            "Number of elements through wall": 1,
+            "Number of elements through shell": 1,
             "Target element density along longest segment": 5.0,
             "Show trim surfaces": False,
             "Use Core": True,
-            "Number of elements across head core major": 6,
-            "Number of elements across torso core major": 6,
-            "Number of elements across arm core major": 4,
-            "Number of elements across leg core major": 4,
+            "Number of elements across core box minor": 2,
             "Number of elements across core transition": 1
         }
 
         if "Human 2 Medium" in parameterSetName:
             options["Number of elements around head"] = 16
             options["Number of elements around torso"] = 16
-            options["Number of elements around arm"] = 12
+            options["Number of elements around arm"] = 8
             options["Number of elements around leg"] = 12
-
-            options["Number of elements across head core major"] = 8
-            options["Number of elements across torso core major"] = 8
-            options["Number of elements across arm core major"] = 6
-            options["Number of elements across leg core major"] = 6
-
+            options["Target element density along longest segment"] = 8.0
+            options["Number of elements across core box minor"] = 2
         elif "Human 2 Fine" in parameterSetName:
             options["Number of elements around head"] = 24
             options["Number of elements around torso"] = 24
-            options["Number of elements around arm"] = 20
-            options["Number of elements around leg"] = 20
-
-            options["Number of elements across head core major"] = 10
-            options["Number of elements across torso core major"] = 10
-            options["Number of elements across arm core major"] = 8
-            options["Number of elements across leg core major"] = 8
+            options["Number of elements around arm"] = 12
+            options["Number of elements around leg"] = 16
+            options["Number of elements through shell"] = 1
+            options["Target element density along longest segment"] = 10.0
+            options["Number of elements across core box minor"] = 4
 
         return options
 
@@ -268,14 +259,12 @@ def getOrderedOptionNames(cls):
             "Number of elements around torso",
             "Number of elements around arm",
             "Number of elements around leg",
-            "Number of elements through wall",
+            "Number of elements through shell",
             "Target element density along longest segment",
             "Show trim surfaces",
             "Use Core",
-            "Number of elements across head core major",
-            "Number of elements across torso core major",
-            "Number of elements across arm core major",
-            "Number of elements across leg core major"]
+            "Number of elements across core box minor",
+            "Number of elements across core transition"]
         return optionNames
 
     @classmethod
@@ -311,8 +300,10 @@ def getOptionScaffoldPackage(cls, optionName, scaffoldType, parameterSetName=Non
 
     @classmethod
     def checkOptions(cls, options):
+        dependentChanges = False
         if not options["Network layout"].getScaffoldType() in cls.getOptionValidScaffoldTypes("Network layout"):
             options["Network layout"] = cls.getOptionScaffoldPackage('Network layout', MeshType_1d_network_layout1)
+        minElementsCountAround = None
         for key in [
             "Number of elements around head",
             "Number of elements around torso",
@@ -321,23 +312,32 @@ def checkOptions(cls, options):
         ]:
             if options[key] < 8:
                 options[key] = 8
+            elif options[key] % 4:
+                options[key] += 4 - (options[key] % 4)
+            if (minElementsCountAround is None) or (options[key] < minElementsCountAround):
+                minElementsCountAround = options[key]
 
-        for key in [
-            "Number of elements across head core major",
-            "Number of elements across torso core major",
-            "Number of elements across arm core major",
-            "Number of elements across leg core major"
-        ]:
-            if options[key] < 4:
-                options[key] = 4
-
-        if options["Number of elements through wall"] < 0:
-            options["Number of elements through wall"] = 1
+        if options["Number of elements through shell"] < 0:
+            options["Number of elements through shell"] = 1
 
         if options["Target element density along longest segment"] < 1.0:
             options["Target element density along longest segment"] = 1.0
 
-        dependentChanges = False
+        if options["Number of elements across core transition"] < 1:
+            options["Number of elements across core transition"] = 1
+
+        maxElementsCountCoreBoxMinor = minElementsCountAround // 2 - 2
+        for key in [
+            "Number of elements across core box minor"
+        ]:
+            if options[key] < 2:
+                options[key] = 2
+            elif options[key] > maxElementsCountCoreBoxMinor:
+                options[key] = maxElementsCountCoreBoxMinor
+                dependentChanges = True
+            elif options[key] % 2:
+                options[key] += options[key] % 2
+
         return dependentChanges
 
     @classmethod
@@ -348,39 +348,42 @@ def generateBaseMesh(cls, region, options):
         :param options: Dict containing options. See getDefaultOptions().
         :return: list of AnnotationGroup, None
         """
-        parameterSetName = options['Base parameter set']
-        isHuman = parameterSetName in ["Default", "Human 1", "Human 2 Coarse", "Human 2 Medium", "Human 2 Fine"]
+        # parameterSetName = options['Base parameter set']
+        # isHuman = parameterSetName in ["Default", "Human 1", "Human 2 Coarse", "Human 2 Medium", "Human 2 Fine"]
 
         layoutRegion = region.createRegion()
         networkLayout = options["Network layout"]
         networkLayout.generate(layoutRegion)  # ask scaffold to generate to get user-edited parameters
         layoutAnnotationGroups = networkLayout.getAnnotationGroups()
         networkMesh = networkLayout.getConstructionObject()
 
-        annotationElementsCountsAround = []
-        annotationElementsCountsAcross = []
+        coreBoxMinorCount = options["Number of elements across core box minor"]
+        coreTransitionCount = options['Number of elements across core transition']
+        annotationAroundCounts = []
+        # implementation currently uses major count including transition
+        annotationCoreMajorCounts = []
         for layoutAnnotationGroup in layoutAnnotationGroups:
-            elementsCountAround = 0
-            elementsCountAcrossMajor = 0
+            aroundCount = 0
+            coreMajorCount = 0
             name = layoutAnnotationGroup.getName()
             if name in ["head", "torso", "arm", "leg"]:
-                elementsCountAround = options["Number of elements around " + name]
-                elementsCountAcrossMajor = options["Number of elements across " + name + " core major"]
-            annotationElementsCountsAround.append(elementsCountAround)
-            annotationElementsCountsAcross.append(elementsCountAcrossMajor)
+                aroundCount = options["Number of elements around " + name]
+                coreMajorCount = aroundCount // 2 - coreBoxMinorCount + 2 * coreTransitionCount
+            annotationAroundCounts.append(aroundCount)
+            annotationCoreMajorCounts.append(coreMajorCount)
 
         isCore = options["Use Core"]
 
         tubeNetworkMeshBuilder = TubeNetworkMeshBuilder(
             networkMesh,
             targetElementDensityAlongLongestSegment=options["Target element density along longest segment"],
             defaultElementsCountAround=options["Number of elements around head"],
-            elementsCountThroughWall=options["Number of elements through wall"],
+            elementsCountThroughWall=options["Number of elements through shell"],
             layoutAnnotationGroups=layoutAnnotationGroups,
-            annotationElementsCountsAround=annotationElementsCountsAround,
-            defaultElementsCountAcrossMajor=options['Number of elements across head core major'],
-            elementsCountTransition=options['Number of elements across core transition'],
-            annotationElementsCountsAcrossMajor=annotationElementsCountsAcross,
+            annotationElementsCountsAround=annotationAroundCounts,
+            defaultElementsCountAcrossMajor=annotationCoreMajorCounts[-1],
+            elementsCountTransition=coreTransitionCount,
+            annotationElementsCountsAcrossMajor=annotationCoreMajorCounts,
             isCore=isCore)
 
         tubeNetworkMeshBuilder.build()

src/scaffoldmaker/utils/interpolation.py

@@ -138,6 +138,44 @@ def computeCubicHermiteDerivativeScaling(v1, d1, v2, d2):
     print('computeCubicHermiteDerivativeScaling:  Max iters reached:', iters, ' mag', mag, 'arc', arcLength)
     return scaling
 
+def computeCubicHermiteStartDerivative(v1, d1_in, v2, d2):
+    """
+    Compute scaled d1 which makes their sum of d1 and d2 magnitudes twice the arc length.
+    :param d1_in: Original start derivative.
+    :return: Scaled d1
+    """
+    d2_mag = magnitude(d2)
+    d1 = set_magnitude(d1_in, 0.5 * d2_mag)
+    for iters in range(100):
+        arcLength = getCubicHermiteArcLength(v1, d1, v2, d2)
+        d1_mag = 2.0 * arcLength - d2_mag
+        d1 = set_magnitude(d1_in, d1_mag)
+        if math.fabs(2.0 * arcLength - d1_mag - d2_mag) < (1.0E-6 * arcLength):
+            break
+    else:
+        print('computeCubicHermiteStartDerivative:  Max iters reached:', iters, ' mag d1', d1_mag, 'mag d2', d2_mag,
+              'arc length', arcLength)
+    return d1
+
+def computeCubicHermiteEndDerivative(v1, d1, v2, d2_in):
+    """
+    Compute scaled d2 which makes their sum of d1 and d2 magnitudes twice the arc length.
+    :param d2_in: Original end derivative.
+    :return: Scaled d2
+    """
+    d1_mag = magnitude(d1)
+    d2 = set_magnitude(d2_in, 0.5 * d1_mag)
+    for iters in range(100):
+        arcLength = getCubicHermiteArcLength(v1, d1, v2, d2)
+        d2_mag = 2.0 * arcLength - d1_mag
+        d2 = set_magnitude(d2_in, d2_mag)
+        if math.fabs(2.0 * arcLength - d1_mag - d2_mag) < (1.0E-6 * arcLength):
+            break
+    else:
+        print('computeCubicHermiteEndDerivative:  Max iters reached:', iters, ' mag d1', d1_mag, 'mag d2', d2_mag,
+              'arc length', arcLength)
+    return d2
+
 def getCubicHermiteArcLength(v1, d1, v2, d2):
     """
     Note this is approximate.
@@ -725,11 +763,17 @@ def smoothCubicHermiteDerivativesLine(nx, nd1,
     componentRange = range(componentsCount)
     if elementsCount == 1:
         # special cases for one element
+        if fixStartDerivative and fixEndDerivative:
+            return nd1
         if not (fixStartDerivative or fixEndDerivative or fixStartDirection or fixEndDirection or fixAllDirections):
             # straight line
             delta = [ (nx[1][c] - nx[0][c]) for c in componentRange ]
             return [ delta, copy.deepcopy(delta) ]
-        if fixAllDirections or (fixStartDirection and fixEndDirection):
+        if fixAllDirections or ((fixStartDirection or fixStartDerivative) and (fixEndDirection or fixEndDerivative)):
+            if fixStartDerivative:
+                return [nd1[0], computeCubicHermiteEndDerivative(nx[0], nd1[0], nx[1], nd1[1])]
+            if fixEndDerivative:
+                return [computeCubicHermiteStartDerivative(nx[0], nd1[0], nx[1], nd1[1]), nd1[1]]
             # fixed directions, equal magnitude
             arcLength = computeCubicHermiteArcLength(nx[0], nd1[0], nx[1], nd1[1], rescaleDerivatives=True)
             return [ set_magnitude(nd1[0], arcLength), set_magnitude(nd1[1], arcLength) ]