Merge remote-tracking branch 'abi/main' into DRG

src/scaffoldmaker/meshtypes/meshtype_3d_esophagus1.py

@@ -36,11 +36,11 @@ def getDefaultNetworkLayoutScaffoldPackage(cls, parameterSetName):
                 ['coordinates'],
                 [Node.VALUE_LABEL_VALUE, Node.VALUE_LABEL_D_DS1, Node.VALUE_LABEL_D_DS2, Node.VALUE_LABEL_D2_DS1DS2,
                  Node.VALUE_LABEL_D_DS3, Node.VALUE_LABEL_D2_DS1DS3], [[
-                    (1, [[0.394, -100.872, 1402.818], [-0.035, 12.367, -48.020], [8.730, -0.526, -0.142], [0.613, -0.153, -0.037], [-0.272, -4.224, -1.088], [-0.169, -1.491, -0.564]]),
-                    (2, [[0.520, -86.043, 1340.066], [0.501, 16.682, -77.602], [9.142, -0.799, -0.113], [0.212, -0.392, 0.096], [-0.465, -5.159, -1.112], [-0.215, -0.377, 0.515]]),
-                    (3, [[1.368, -67.733, 1247.932], [0.235, -3.685, -89.672], [9.061, -1.366, 0.080], [-0.833, -0.231, 0.187], [-0.714, -4.722, 0.192], [-0.167, 0.445, 1.659]]),
-                    (4, [[0.361, -91.057, 1165.531], [-2.499, -24.560, -49.102], [7.540, -1.290, 0.261], [-0.809, 1.514, 2.095], [-0.806, -4.269, 2.176], [0.001, 0.896, 0.910]]),
-                    (5, [[11.750, -111.874, 1127.887], [7.636, -5.715, -7.930], [5.678, 1.265, 4.556], [-8.397, 13.092, 24.878], [-0.708, -3.530, 1.862], [-0.807, -7.995, 7.596]])
+                    (1, [[0.394,-100.872,1402.818], [-0.035,12.367,-48.020], [8.730,-0.526,-0.142], [0.645,-0.149,-0.084], [-0.272,-4.224,-1.088], [-0.172,-1.936,0.474]]),
+                    (2, [[0.520,-86.043,1340.066], [0.501,16.682,-77.602], [9.142,-0.799,-0.113], [0.176,-0.353,-0.027], [-0.465,-5.159,-1.112], [-0.217,0.101,-0.821]]),
+                    (3, [[1.368,-67.733,1247.932], [0.235,-3.685,-89.672], [9.061,-1.366,0.080], [-0.816,-0.266,0.913], [-0.714,-4.722,0.192], [-0.191,0.388,4.146]]),
+                    (4, [[0.361,-91.057,1165.531], [-2.499,-24.560,-49.102], [7.540,-1.290,0.261], [-1.321,2.063,3.911], [-0.806,-4.269,2.176], [0.019,0.649,0.615]]),
+                    (5, [[8.225,-110.042,1130.494], [17.674,-15.065,-22.956], [6.092,0.940,4.074], [1.229,3.363,5.229], [-0.669,-3.559,1.821], [0.256,-3.730,2.087]])
                 ]]),
 
             'userAnnotationGroups': [

src/scaffoldmaker/meshtypes/meshtype_3d_ostium2.py

@@ -1005,7 +1005,7 @@ def generateOstiumMesh(region, options, trackSurface, centralPath, startNodeIden
             forceMidLinearXi3=not useCubicHermiteThroughVesselWall,
             maxStartThickness=vesselWallThickness, maxEndThickness=vesselWallThickness,
             elementsCountRadial=elementsCountAlong, meshGroups=rowMeshGroups, wallAnnotationGroups=wallAnnotationGroups,
-            coordinates=coordinates)
+            coordinates=coordinates, rescaleStartDerivatives=True)
 
         # To visualise central path
         # nodeIdentifierLine = nodeIdentifier

src/scaffoldmaker/meshtypes/meshtype_3d_stomach1.py

@@ -123,21 +123,21 @@ def getDefaultNetworkLayoutScaffoldPackage(cls, parameterSetName):
             'meshEdits': exnode_string_from_nodeset_field_parameters(
                 ['coordinates'],
                 [Node.VALUE_LABEL_VALUE, Node.VALUE_LABEL_D_DS1, Node.VALUE_LABEL_D_DS2, Node.VALUE_LABEL_D2_DS1DS2, Node.VALUE_LABEL_D_DS3, Node.VALUE_LABEL_D2_DS1DS3], [[
-                (1, [[11.750,-111.874,1127.887], [7.636,-5.715,-7.930], [5.678,1.265,4.556], [-8.397,13.092,24.878], [-0.708,-3.530,1.862], [-0.807,-7.995,7.596]]),
-                (2, [[23.789,-117.922,1120.040], [26.354,-6.724,-6.404], [4.223,6.205,10.864], [10.037,1.800,8.968], [-1.192,-11.215,6.869], [-2.926,-13.889,10.204]]),
-                (3, [[63.704,-120.094,1123.374], [[0.500,-9.138,-13.405],[50.106,1.267,11.056]], [[37.742,-3.477,3.778],[-2.509,7.605,10.499]], [[3.190,-0.590,-0.290],[3.190,-0.590,-0.290]], [[-5.452,-34.121,23.056],[-1.379,-10.790,7.486]], [[-0.100,-1.220,0.550],[-0.100,-1.220,0.550]]]),
-                (4, [[61.247,-99.931,1152.681], [0.346,-2.728,-3.873], [11.320,-0.365,1.269], [15.701,-3.029,2.832], [-0.653,-5.931,4.119], [-5.756,-22.456,14.946]]),
-                (5, [[61.743,-103.510,1147.760], [0.413,-3.592,-5.311], [24.159,-2.387,3.493], [10.039,-1.491,0.928], [-2.982,-15.339,10.142], [-2.384,-12.784,8.034]]),
-                (6, [[62.381,-107.527,1141.785], [0.249,-4.737,-7.073], [30.559,-2.973,3.067], [5.720,-0.448,-0.071], [-3.839,-23.420,15.550], [-0.206,-5.098,2.599]]),
-                (7, [[62.800,-113.150,1133.665], [0.116,-6.546,-9.651], [35.408,-3.630,2.888], [4.541,-0.488,0.325], [-4.677,-29.659,20.061], [-0.383,-3.012,1.219]]),
-                (8, [[64.339,-131.197,1107.233], [0.086,-11.682,-16.915], [39.201,-3.705,2.758], [0.490,-2.470,-4.030], [-5.108,-35.712,24.638], [0.260,0.760,0.250]]),
-                (9, [[62.912,-143.954,1088.811], [-5.216,-12.408,-17.967], [34.623,-8.161,-4.415], [-4.950,-4.940,-8.800], [-4.917,-34.532,25.275], [0.770,0.610,1.050]]),
-                (10, [[53.361,-155.397,1072.006], [-15.833,-9.775,-15.486], [25.117,-13.916,-16.896], [-14.190,-3.410,-7.580], [-2.355,-30.712,21.794], [2.030,4.350,-1.240]]),
-                (11, [[32.110,-162.230,1059.680], [-22.173,-3.681,-8.142], [10.067,-16.126,-20.126], [-15.400,0.330,-0.620], [-2.559,-23.629,17.653], [1.200,8.620,-5.080]]),
-                (12, [[10.560,-162.970,1055.650], [-20.956,2.164,-0.908], [-0.696,-13.976,-17.241], [-8.750,2.830,4.760], [-2.502,-18.048,14.732], [-1.000,7.490,-3.260]]),
-                (13, [[-8.740,-158.280,1057.630], [-17.269,8.027,2.588], [-6.253,-10.164,-10.200], [-2.630,4.430,6.290], [-3.442,-11.913,13.981], [-1.380,7.080,-2.660]]),
-                (14, [[-23.260,-147.620,1060.640], [-11.787,11.369,1.720], [-5.807,-5.321,-4.621], [-0.700,3.020,1.710], [-3.207,-4.764,9.515], [-2.080,5.830,-2.460]]),
-                (15, [[-32.281,-136.261,1061.249], [-6.125,11.114,-0.493], [-5.735,-3.374,-4.823], [-2.070,-0.220,-4.520], [-4.630,-2.238,7.071], [-4.400,2.440,0.570]]) 
+                (1, [[8.225,-110.042,1130.494], [17.674,-15.065,-22.956], [6.092,0.940,4.074], [1.229,3.363,5.229], [-0.669,-3.559,1.821], [0.256,-3.730,2.087]]),
+                (2, [[27.268,-119.049,1120.363], [29.345,-6.635,-5.589], [3.536,6.355,11.021], [-4.563,2.025,2.177], [-1.229,-11.215,6.861], [1.619,-4.227,2.277]]),
+                (3, [[63.704,-120.094,1123.374], [[0.842,-9.011,-13.227],[41.742,4.358,11.135]], [[37.708,-3.017,4.455],[-3.481,6.664,10.443]], [[1.955,-0.619,-0.927],[1.246,-0.837,-1.425]], [[-5.447,-34.193,22.949],[-0.678,-11.218,6.933]], [[-0.269,-2.561,2.982],[-1.558,4.051,-2.388]]]),
+                (4, [[61.247,-99.931,1152.681], [0.427,-3.349,-4.397], [11.314,-0.313,1.337], [16.186,-3.105,3.004], [-0.672,-5.774,4.333], [-3.234,-9.056,7.808]]),
+                (5, [[61.743,-103.510,1147.760], [0.565,-3.807,-5.444], [24.118,-2.068,3.951], [9.485,-0.404,2.240], [-2.996,-15.211,10.328], [-1.475,-9.388,4.780]]),
+                (6, [[62.381,-107.527,1141.785], [0.559,-4.813,-7.053], [30.505,-2.333,4.010], [5.777,-1.924,-2.194], [-3.842,-23.362,15.637], [-0.957,-6.170,6.414]]),
+                (7, [[62.800,-113.150,1133.665], [0.642,-6.292,-9.204], [35.350,-2.675,4.295], [3.358,1.888,3.501], [-4.676,-29.711,19.984], [-0.685,-6.454,2.135]]),
+                (8, [[64.339,-131.197,1107.233], [-0.327,-11.961,-17.331], [39.197,-4.146,2.121], [-1.413,0.507,0.138], [-5.109,-35.662,24.710], [0.316,-0.518,0.671]]),
+                (9, [[62.912,-143.954,1088.811], [-5.607,-12.396,-18.044], [34.482,-8.485,-4.886], [-6.065,-6.479,-11.671], [-4.921,-34.543,25.260], [1.721,3.665,0.276]]),
+                (10, [[53.361,-155.397,1072.006], [-15.827,-9.773,-15.483], [25.119,-13.915,-16.894], [-19.643,-7.564,-13.701], [-2.356,-30.712,21.794], [-0.170,4.599,-5.192]]),
+                (11, [[32.110,-162.230,1059.680], [-22.173,-3.681,-8.142], [10.068,-16.126,-20.126], [-12.488,0.413,0.532], [-2.559,-23.629,17.653], [1.057,6.496,-3.159]]),
+                (12, [[10.560,-162.970,1055.650], [-20.956,2.164,-0.908], [-0.697,-13.976,-17.241], [-8.999,3.378,5.491], [-2.502,-18.048,14.732], [1.647,5.538,-1.785]]),
+                (13, [[-8.740,-158.280,1057.630], [-17.269,8.027,2.588], [-6.253,-10.164,-10.200], [-3.575,5.519,6.887], [-3.442,-11.913,13.981], [-3.280,8.069,-2.019]]),
+                (14, [[-23.260,-147.620,1060.640], [-11.776,11.361,1.718], [-5.807,-5.321,-4.621], [1.528,2.286,2.595], [-3.207,-4.763,9.515], [-0.837,5.004,-3.547]]),
+                (15, [[-32.281,-136.261,1061.249], [-8.825,12.725,0.137], [-5.485,-3.752,-4.834], [-0.599,1.952,0.267], [-4.191,-2.983,7.070], [0.644,-1.559,-1.898]])
                 ]]),
 
             'userAnnotationGroups': [
@@ -496,7 +496,7 @@ def getDefaultNetworkLayoutScaffoldPackage(cls, parameterSetName):
                 }]
         })
 
-def getDefaultOstiumSettings():
+def getDefaultOstiumSettings(parameterSetName):
     """
     Generate list of default options for ostium.
     """
@@ -516,6 +516,8 @@ def getDefaultOstiumSettings():
                 'Refine number of elements around': 4,
                 'Refine number of elements along': 4,
                 'Refine number of elements through wall': 1}
+    if 'Human 2' in parameterSetName:
+        options['Number of elements along'] = 2
 
     return options
 
@@ -627,6 +629,7 @@ def getDefaultOptions(cls, parameterSetName='Default'):
             options['Longitudinal muscle layer relative thickness'] = 0.25
             options['Limiting ridge'] = False
 
+        options['Base parameter set'] = parameterSetName
         return options
 
     @staticmethod
@@ -1301,6 +1304,7 @@ def createStomachMesh3d(region, fm, coordinates, stomachTermsAlong, allAnnotatio
     :return allAnnotationGroups, nextNodeIdentifier, nextElementIdentifier, elementsAlongSections, nodeIdxDistal,
     xDistal, d1Distal, d2Distal, d3Distal, arclengthDuodenumCP, xPrev, d2Prev
     """
+    parameterSetName = options['Base parameter set']
     elementsCountAroundEso = options['Number of elements around esophagus']
     elementsCountAroundDuod = options['Number of elements around duodenum']
     elementsCountAlong = options['Number of elements along']
@@ -1312,7 +1316,7 @@ def createStomachMesh3d(region, fm, coordinates, stomachTermsAlong, allAnnotatio
     useCrossDerivatives = False
     useCubicHermiteThroughWall = not (options['Use linear through wall'])
 
-    ostiumOptions = getDefaultOstiumSettings()
+    ostiumOptions = getDefaultOstiumSettings(parameterSetName)
     GEJSettings = updateOstiumOptions(options, ostiumOptions)
     elementsAlongEsophagus = GEJSettings['Number of elements along']
     elementsThroughEsophagusWall = GEJSettings['Number of elements through wall']
@@ -1322,6 +1326,7 @@ def createStomachMesh3d(region, fm, coordinates, stomachTermsAlong, allAnnotatio
     GEJSettings['Use linear through ostium wall'] = options['Use linear through wall']
     GEJSettings['Use linear through vessel wall'] = options['Use linear through wall']
 
+    stomachStartElementIdentifier = elementIdentifier
     elementsCountAcrossCardia = 1
     cardiaDiameterFactor = 1.4  # scale to ostium diameter
     sf = (cardiaDiameterFactor - 1) * ostiumRadius
@@ -2429,7 +2434,7 @@ def createStomachMesh3d(region, fm, coordinates, stomachTermsAlong, allAnnotatio
     # Create elements
     fundusMucosaElementIdentifiers = []
     elementIdxMat = []
-    n = 0
+    n = stomachStartElementIdentifier - 1
     for n2 in range(elementsAlongEsophagus):
         elementIdxThroughWall = []
         for n3 in range(elementsThroughEsophagusWall):

src/scaffoldmaker/utils/zinc_utils.py

@@ -13,7 +13,6 @@
 from cmlibs.zinc.result import RESULT_OK
 from scaffoldmaker.utils import interpolation as interp
 
-
 def interpolateNodesCubicHermite(cache, coordinates, xi, normal_scale,
         node1, derivative1, scale1, cross_derivative1, cross_scale1,
         node2, derivative2, scale2, cross_derivative2, cross_scale2):
@@ -755,3 +754,36 @@ def group_add_connected_elements(group: FieldGroup, other_mesh_group: MeshGroup)
                     del elementid_to_nodeids[elementid]
                     break
         group.setSubelementHandlingMode(old_subelement_mode)
+
+def blendDerivativeBetweenDifferentSizeElements(coordinates, cache, nodes, prevNodeID, transitionNodeID, nextNodeID):
+    """
+    Blends derivatives for nodes between two elements of different lengths. Calculates an end derivative using
+    upstream coordinates and derivative, and a start derivative using downstream coordinates and derivative. The
+    derivative at the junction is then rescaled with the harmonic mean of start and end derivative magnitudes.
+    :param coordinates: coordinate field
+    :param cache: field cache
+    :param nodes: nodeset
+    :param prevNodeID: node identifier of upstream node
+    :param transitionNodeID: node identifier of node between the two elements of different lengths
+    :param nextNodeID: node identifier of downstream node
+    """
+    node = nodes.findNodeByIdentifier(prevNodeID)
+    cache.setNode(node)
+    xPrev = coordinates.getNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, 3)[1]
+    d2Prev = coordinates.getNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, 3)[1]
+
+    node = nodes.findNodeByIdentifier(nextNodeID)
+    cache.setNode(node)
+    xNext = coordinates.getNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, 3)[1]
+    d2Next = coordinates.getNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, 3)[1]
+
+    node = nodes.findNodeByIdentifier(transitionNodeID)
+    cache.setNode(node)
+    xTransition = coordinates.getNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, 3)[1]
+    d2Transition = coordinates.getNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, 3)[1]
+
+    d2_in = interp.computeCubicHermiteEndDerivative(xPrev, d2Prev, xTransition, d2Transition)
+    d2_out = interp.computeCubicHermiteStartDerivative(xTransition, d2Transition, xNext, d2Next)
+    d2MeanMag = 2.0 / ((1.0 / magnitude(d2_in)) + (1.0 / magnitude(d2_out)))
+    d2Mean = set_magnitude(d2Transition, d2MeanMag)
+    coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, d2Mean)

tests/test_cecum.py

@@ -51,21 +51,21 @@ def test_cecum1(self):
         fieldmodule = region.getFieldmodule()
         self.assertEqual(RESULT_OK, fieldmodule.defineAllFaces())
         mesh3d = fieldmodule.findMeshByDimension(3)
-        self.assertEqual(460, mesh3d.getSize())
+        self.assertEqual(436, mesh3d.getSize())
         mesh2d = fieldmodule.findMeshByDimension(2)
-        self.assertEqual(1738, mesh2d.getSize())
+        self.assertEqual(1642, mesh2d.getSize())
         mesh1d = fieldmodule.findMeshByDimension(1)
-        self.assertEqual(2102, mesh1d.getSize())
+        self.assertEqual(1982, mesh1d.getSize())
         nodes = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES)
-        self.assertEqual(824, nodes.getSize())
+        self.assertEqual(776, nodes.getSize())
         datapoints = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_DATAPOINTS)
         self.assertEqual(0, datapoints.getSize())
 
         coordinates = fieldmodule.findFieldByName("coordinates").castFiniteElement()
         self.assertTrue(coordinates.isValid())
         minimums, maximums = evaluateFieldNodesetRange(coordinates, nodes)
-        assertAlmostEqualList(self, minimums, [-112.4222871639696, -146.3433620526202, 852.5876977230726], 1.0E-6)
-        assertAlmostEqualList(self, maximums, [-54.14347619948218, -77.56, 899.9973429272325], 1.0E-6)
+        assertAlmostEqualList(self, minimums, [-112.41968159644387, -146.34422797225153, 852.6082677676069], 1.0E-6)
+        assertAlmostEqualList(self, maximums, [-54.22370705290374, -77.56, 899.9973429272325], 1.0E-6)
 
         with ChangeManager(fieldmodule):
             one = fieldmodule.createFieldConstant(1.0)
@@ -77,15 +77,15 @@ def test_cecum1(self):
         fieldcache = fieldmodule.createFieldcache()
         result, surfaceArea = surfaceAreaField.evaluateReal(fieldcache, 1)
         self.assertEqual(result, RESULT_OK)
-        self.assertAlmostEqual(surfaceArea, 8554.513081780715, delta=1.0E-6)
+        self.assertAlmostEqual(surfaceArea, 8546.983090282285, delta=1.0E-6)
         result, volume = volumeField.evaluateReal(fieldcache, 1)
         self.assertEqual(result, RESULT_OK)
-        self.assertAlmostEqual(volume, 13809.9960233912, delta=1.0E-6)
+        self.assertAlmostEqual(volume, 13790.25181377472, delta=1.0E-6)
 
         # check some annotationGroups:
         expectedSizes3d = {
-            "caecum": 460,
-            "ileum": 24,
+            "caecum": 436,
+            "ileum": 12,
             "ileocecal junction": 12
         }
 

tests/test_esophagus.py

@@ -58,14 +58,14 @@ def test_esophagus1(self):
         coordinates = fieldmodule.findFieldByName("coordinates").castFiniteElement()
         self.assertTrue(coordinates.isValid())
         minimums, maximums = evaluateFieldNodesetRange(coordinates, nodes)
-        assertAlmostEqualList(self, minimums, [-8.645606404075409, -115.4469482741834, 1123.3296869517246], 1.0E-6)
-        assertAlmostEqualList(self, maximums, [17.428066853931824, -63.05159903323634, 1403.9055789269714], 1.0E-6)
+        assertAlmostEqualList(self, minimums, [-8.651939179594637, -113.60141744677131, 1126.4189460716025], 1.0E-6)
+        assertAlmostEqualList(self, maximums, [14.317229342198544, -62.96131235320371, 1403.9055789269714], 1.0E-6)
 
         flatCoordinates = fieldmodule.findFieldByName("flat coordinates").castFiniteElement()
         self.assertTrue(flatCoordinates.isValid())
         minimums, maximums = evaluateFieldNodesetRange(flatCoordinates, nodes)
-        assertAlmostEqualList(self, minimums, [-2.023676644962851, 0.0, 0.0], 1.0E-6)
-        assertAlmostEqualList(self, maximums, [38.56398979676183, 289.9673186862639, 1.2], 1.0E-6)
+        assertAlmostEqualList(self, minimums, [-2.0052709915460802, 0.0, 0.0], 1.0E-6)
+        assertAlmostEqualList(self, maximums, [38.55018555669926, 286.2207860825292, 1.2], 1.0E-6)
 
         esophagusCoordinates = fieldmodule.findFieldByName("esophagus coordinates").castFiniteElement()
         minimums, maximums = evaluateFieldNodesetRange(esophagusCoordinates, nodes)
@@ -82,10 +82,10 @@ def test_esophagus1(self):
         fieldcache = fieldmodule.createFieldcache()
         result, surfaceArea = surfaceAreaField.evaluateReal(fieldcache, 1)
         self.assertEqual(result, RESULT_OK)
-        self.assertAlmostEqual(surfaceArea, 12544.147466204333, delta=1.0E-6)
+        self.assertAlmostEqual(surfaceArea, 12404.496450194649, delta=1.0E-6)
         result, volume = volumeField.evaluateReal(fieldcache, 1)
         self.assertEqual(result, RESULT_OK)
-        self.assertAlmostEqual(volume, 13787.433711373638, delta=1.0E-6)
+        self.assertAlmostEqual(volume, 13637.443048857438, delta=1.0E-6)
 
         # check some annotationGroups:
         expectedSizes3d = {