Merge pull request #93 from elias-soltani/path1/sidederivatives

Path1/sidederivatives

src/scaffoldmaker/meshtypes/meshtype_1d_path1.py

@@ -11,6 +11,9 @@
 from opencmiss.zinc.node import Node
 from scaffoldmaker.meshtypes.scaffold_base import Scaffold_base
 from scaffoldmaker.utils.interpolation import DerivativeScalingMode, smoothCubicHermiteDerivativesLine
+from scaffoldmaker.utils import vector
+from opencmiss.zinc.result import RESULT_OK
+
 
 class MeshType_1d_path1(Scaffold_base):
     '''
@@ -24,6 +27,8 @@ def getName():
     def getDefaultOptions(parameterSetName='Default'):
         return {
             'Coordinate dimensions' : 3,
+            'D2 derivatives': False,
+            'D3 derivatives': False,
             'Length' : 1.0,
             'Number of elements' : 1
         }
@@ -32,19 +37,24 @@ def getDefaultOptions(parameterSetName='Default'):
     def getOrderedOptionNames():
         return [
             'Coordinate dimensions',
+            'D2 derivatives',
+            'D3 derivatives',
             'Length',
             'Number of elements'
         ]
 
     @staticmethod
     def checkOptions(options):
+        dependentChanges = False
         if (options['Coordinate dimensions'] < 1) :
             options['Coordinate dimensions'] = 1
         elif (options['Coordinate dimensions'] > 3) :
             options['Coordinate dimensions'] = 3
         if (options['Number of elements'] < 1) :
             options['Number of elements'] = 1
 
+        return dependentChanges
+
     @classmethod
     def generateBaseMesh(cls, region, options):
         """
@@ -53,6 +63,8 @@ def generateBaseMesh(cls, region, options):
         :return: [] empty list of AnnotationGroup
         """
         coordinateDimensions = options['Coordinate dimensions']
+        d2derivatives = options['D2 derivatives']
+        d3derivatives = options['D3 derivatives']
         length = options['Length']
         elementsCount = options['Number of elements']
 
@@ -69,22 +81,34 @@ def generateBaseMesh(cls, region, options):
         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)
-        nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D2_DS1DS2, 1)
+        if d2derivatives:
+            nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS2, 1)
+            nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D2_DS1DS2, 1)
+        if d3derivatives:
+            nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D_DS3, 1)
+            nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_D2_DS1DS3, 1)
 
         nodeIdentifier = 1
         x = [ 0.0, 0.0, 0.0 ]
         dx_ds1 = [ length/elementsCount, 0.0, 0.0 ]
-        dx_ds2 = [ 0.0, 1.0, 0.0 ]
-        d2x_ds1ds2 = [ 0.0, 0.0, 0.0 ]
+        if d2derivatives:
+            dx_ds2 = [ 0.0, 1.0, 0.0 ]
+            d2x_ds1ds2 = [ 0.0, 0.0, 0.0 ]
+        if d3derivatives:
+            dx_ds3 = [ 0.0, 0.0, 1.0 ]
+            d2x_ds1ds3 = [ 0.0, 0.0, 0.0 ]
         for n in range(elementsCount + 1):
             x[0] = length*n/elementsCount
             node = nodes.createNode(nodeIdentifier, nodetemplate)
             cache.setNode(node)
             coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, x)
             coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, dx_ds1)
-            coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, dx_ds2)
-            coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS2, 1, d2x_ds1ds2)
+            if d2derivatives:
+                coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, dx_ds2)
+                coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS2, 1, d2x_ds1ds2)
+            if d3derivatives:
+                coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS3, 1, dx_ds3)
+                coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS3, 1, d2x_ds1ds3)
             nodeIdentifier = nodeIdentifier + 1
 
         #################
@@ -108,26 +132,58 @@ def generateBaseMesh(cls, region, options):
 
     @classmethod
     def smoothPath(cls, region, options, mode : DerivativeScalingMode):
-        x, d1 = extractPathParametersFromRegion(region)[0:2]
+        x, d1 = extractPathParametersFromRegion(region, [Node.VALUE_LABEL_VALUE, Node.VALUE_LABEL_D_DS1])
         d1 = smoothCubicHermiteDerivativesLine(x, d1, magnitudeScalingMode=mode)
         setPathParameters(region, [ Node.VALUE_LABEL_VALUE, Node.VALUE_LABEL_D_DS1 ], [ x, d1 ])
 
+    @classmethod
+    def makeD2Normal(cls, region, options):
+        if not options['D2 derivatives']:
+            return
+        d1, d2 = extractPathParametersFromRegion(region, [Node.VALUE_LABEL_D_DS1, Node.VALUE_LABEL_D_DS2])
+        for c in range(len(d1)):
+            td2 = vector.vectorRejection(d2[c], d1[c])
+            d2[c] = vector.setMagnitude(td2, vector.magnitude(d2[c]))
+        setPathParameters(region, [Node.VALUE_LABEL_D_DS2], [d2])
+
+    @classmethod
+    def makeD3Normal(cls, region, options):
+        if not options['D3 derivatives']:
+            return
+        if options['D2 derivatives']:
+            d1, d2, d3 = extractPathParametersFromRegion(region, [Node.VALUE_LABEL_D_DS1, Node.VALUE_LABEL_D_DS2,
+                                                                  Node.VALUE_LABEL_D_DS3])
+            for c in range(len(d1)):
+                d3[c] = vector.setMagnitude(vector.crossproduct3(d1[c], d2[c]), vector.magnitude(d3[c]))
+            setPathParameters(region, [Node.VALUE_LABEL_D_DS3], [d3])
+        else:
+            d1, d3 = extractPathParametersFromRegion(region, [Node.VALUE_LABEL_D_DS1, Node.VALUE_LABEL_D_DS3])
+            for c in range(len(d1)):
+                td3 = vector.vectorRejection(d3[c], d1[c])
+                d3[c] = vector.setMagnitude(td3, vector.magnitude(d3[c]))
+            setPathParameters(region, [Node.VALUE_LABEL_D_DS3], [d3])
+
+
     @classmethod
     def getInteractiveFunctions(cls):
         """
         Supply client with functions for smoothing path parameters.
         """
         return [
-            ("Smooth d1 arithmetic", lambda region, options: cls.smoothPath(region, options, DerivativeScalingMode.ARITHMETIC_MEAN)),
-            ("Smooth d1 harmonic", lambda region, options: cls.smoothPath(region, options, DerivativeScalingMode.HARMONIC_MEAN)) ]
+            ("Smooth D1 arithmetic", lambda region, options: cls.smoothPath(region, options, DerivativeScalingMode.ARITHMETIC_MEAN)),
+            ("Smooth D1 harmonic", lambda region, options: cls.smoothPath(region, options, DerivativeScalingMode.HARMONIC_MEAN)),
+            ("Make D2 normal", lambda region, options: cls.makeD2Normal(region, options)),
+            ("Make D3 normal", lambda region, options: cls.makeD3Normal(region, options))
+        ]
 
 
-def extractPathParametersFromRegion(region, groupName=None):
+def extractPathParametersFromRegion(region, valueLabels, groupName=None):
     '''
     Returns parameters of all nodes in region in identifier order.
     Assumes nodes in region have field coordinates (1 to 3 components).
     Currently limited to nodes with exactly value, d_ds1, d_ds2, d2_ds12,
     same as path 1 scaffold.
+    :param valueLabels: List of parameters required as list of node value labels. e.g. [Node.VALUE_LABEL_VALUE, Node.VALUE_LABEL_D_DS1].
     :param groupName: Optional name of Zinc group to get parameters from.
     :return: cx, cd1, cd2, cd12 (all padded with zeroes to 3 components)
     '''
@@ -136,10 +192,9 @@ def extractPathParametersFromRegion(region, groupName=None):
     componentsCount = coordinates.getNumberOfComponents()
     assert componentsCount in [ 1, 2, 3 ], 'extractPathParametersFromRegion.  Invalid coordinates number of components'
     cache = fieldmodule.createFieldcache()
-    cx = []
-    cd1 = []
-    cd2 = []
-    cd12 = []
+
+    valueLabelsCount = len(valueLabels)
+    returnValues = [[] for i in range(valueLabelsCount)]
     nodes = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES)
     if groupName:
         group = fieldmodule.findFieldByName(groupName).castGroup()
@@ -152,21 +207,14 @@ def extractPathParametersFromRegion(region, groupName=None):
     node = nodeIter.next()
     while node.isValid():
         cache.setNode(node)
-        result, x  = coordinates.getNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, componentsCount)
-        result, d1 = coordinates.getNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS1, 1, componentsCount)
-        result, d2 = coordinates.getNodeParameters(cache, -1, Node.VALUE_LABEL_D_DS2, 1, componentsCount)
-        result, d12 = coordinates.getNodeParameters(cache, -1, Node.VALUE_LABEL_D2_DS1DS2, 1, componentsCount)
-        for c in range(componentsCount, 3):
-            x.append(0.0)
-            d1.append(0.0)
-            d2.append(0.0)
-            d12.append(0.0)
-        cx.append(x)
-        cd1.append(d1)
-        cd2.append(d2)
-        cd12.append(d12)
+        for i in range(valueLabelsCount):
+            result, values = coordinates.getNodeParameters(cache, -1, valueLabels[i], 1, componentsCount)
+            for c in range(componentsCount, 3):
+                values.append(0.0)
+            returnValues[i].append(values)
         node = nodeIter.next()
-    return cx, cd1, cd2, cd12
+
+    return returnValues
 
 
 def setPathParameters(region, nodeValueLabels, nodeValues):

src/scaffoldmaker/meshtypes/meshtype_3d_cecum1.py

@@ -41,6 +41,7 @@ class MeshType_3d_cecum1(Scaffold_base):
         'Pig 1': ScaffoldPackage(MeshType_1d_path1, {
             'scaffoldSettings': {
                 'Coordinate dimensions': 3,
+                'D2 derivatives': True,
                 'Length': 120.0,
                 'Number of elements': 3
             },
@@ -307,7 +308,9 @@ def generateBaseMesh(cls, region, options):
         # Central path
         tmpRegion = region.createRegion()
         centralPath.generate(tmpRegion)
-        cx, cd1, cd2, cd12 = extractPathParametersFromRegion(tmpRegion)
+        cx, cd1, cd2, cd12 = extractPathParametersFromRegion(tmpRegion,
+                                                             [Node.VALUE_LABEL_VALUE, Node.VALUE_LABEL_D_DS1,
+                                                              Node.VALUE_LABEL_D_DS2, Node.VALUE_LABEL_D2_DS1DS2])
         # for i in range(len(cx)):
         #     print(i, '[', cx[i], ',', cd1[i], ',', cd2[i], ',', cd12[i], '],')
         del tmpRegion

src/scaffoldmaker/meshtypes/meshtype_3d_colon1.py

@@ -30,6 +30,7 @@ class MeshType_3d_colon1(Scaffold_base):
         'Human 1' : ScaffoldPackage(MeshType_1d_path1, {
             'scaffoldSettings' : {
                 'Coordinate dimensions' : 3,
+                'D2 derivatives': True,
                 'Length' : 1.0,
                 'Number of elements' : 8
                 },
@@ -48,6 +49,7 @@ class MeshType_3d_colon1(Scaffold_base):
         'Human 2' : ScaffoldPackage(MeshType_1d_path1, {
             'scaffoldSettings' : {
                 'Coordinate dimensions' : 3,
+                'D2 derivatives': True,
                 'Length' : 1.0,
                 'Number of elements' : 8
                 },
@@ -66,6 +68,7 @@ class MeshType_3d_colon1(Scaffold_base):
         'Mouse 1' : ScaffoldPackage(MeshType_1d_path1, {
             'scaffoldSettings' : {
                 'Coordinate dimensions' : 3,
+                'D2 derivatives': True,
                 'Length' : 1.0,
                 'Number of elements' : 7
                 },
@@ -83,6 +86,7 @@ class MeshType_3d_colon1(Scaffold_base):
         'Mouse 2' : ScaffoldPackage(MeshType_1d_path1, {
             'scaffoldSettings' : {
                 'Coordinate dimensions' : 3,
+                'D2 derivatives': True,
                 'Length' : 1.0,
                 'Number of elements' : 4
                 },
@@ -97,6 +101,7 @@ class MeshType_3d_colon1(Scaffold_base):
         'Pig 1' : ScaffoldPackage(MeshType_1d_path1, {
             'scaffoldSettings' : {
                 'Coordinate dimensions' : 3,
+                'D2 derivatives': True,
                 'Length' : 1.0,
                 'Number of elements' : 36
                 },
@@ -143,6 +148,7 @@ class MeshType_3d_colon1(Scaffold_base):
         'Pig 2': ScaffoldPackage(MeshType_1d_path1, {
             'scaffoldSettings': {
                 'Coordinate dimensions': 3,
+                'D2 derivatives': True,
                 'Length': 90.0,
                 'Number of elements': 3
             },
@@ -391,7 +397,9 @@ def generateBaseMesh(cls, region, options):
         # Central path
         tmpRegion = region.createRegion()
         centralPath.generate(tmpRegion)
-        cx, cd1, cd2, cd12 = extractPathParametersFromRegion(tmpRegion)
+        cx, cd1, cd2, cd12 = extractPathParametersFromRegion(tmpRegion,
+                                                             [Node.VALUE_LABEL_VALUE, Node.VALUE_LABEL_D_DS1,
+                                                              Node.VALUE_LABEL_D_DS2, Node.VALUE_LABEL_D2_DS1DS2])
         # for i in range(len(cx)):
         #     print(i, '[', cx[i], ',', cd1[i], ',', cd2[i], ',', cd12[i], '],')
         del tmpRegion

src/scaffoldmaker/meshtypes/meshtype_3d_smallintestine1.py

@@ -29,6 +29,7 @@ class MeshType_3d_smallintestine1(Scaffold_base):
     centralPathDefaultScaffoldPackages = {
         'Mouse 1' : ScaffoldPackage(MeshType_1d_path1, {
             'scaffoldSettings' : {
+                'D2 derivatives': True,
                 'Coordinate dimensions' : 3,
                 'Length' : 1.0,
                 'Number of elements' : 45
@@ -231,7 +232,9 @@ def generateBaseMesh(cls, region, options):
         # Central path
         tmpRegion = region.createRegion()
         centralPath.generate(tmpRegion)
-        cx, cd1, cd2, cd12 = extractPathParametersFromRegion(tmpRegion)
+        cx, cd1, cd2, cd12 = extractPathParametersFromRegion(tmpRegion,
+                                                             [Node.VALUE_LABEL_VALUE, Node.VALUE_LABEL_D_DS1,
+                                                              Node.VALUE_LABEL_D_DS2, Node.VALUE_LABEL_D2_DS1DS2])
         # for i in range(len(cx)):
         #     print(i, '[', cx[i], ',', cd1[i], ',', cd2[i],',', cd12[i], '],')
         del tmpRegion

src/scaffoldmaker/utils/vector.py

@@ -42,3 +42,49 @@ def addVectors(v1,v2,s1=1.0,s2=1.0):
     :return: Vector s1*v1+s2*v2
     '''
     return [(s1 * v1[c] + s2 * v2[c]) for c in range(len(v1))]
+
+
+def scalarProjection(v1, v2):
+    """
+    :return: Scalar projection of v1 onto v2.
+    """
+    return dotproduct(v1, normalise(v2))
+
+
+def vectorProjection(v1, v2):
+    """
+    Calculate vector projection of v1 on v2
+    :return: A projection vector.
+    """
+    s1 = scalarProjection(v1, v2)
+    return scalarProduct(s1, normalise(v2))
+
+
+def vectorRejection(v1, v2):
+    """
+    Calculate vector rejection of v1 on v2
+    :return: A rejection vector.
+    """
+    v1p = vectorProjection(v1, v2)
+    return addVectors(v1, v1p, 1.0, -1.0)
+
+
+def scalarProduct(s, v):
+    """
+    Calculate s * v
+    :param s: Scalar.
+    :param v: Vector.
+    :return:
+    """
+    return [s * v[c] for c in range(len(v))]
+
+
+def parallelVectors(v1, v2):
+    """
+    :return: True if the vectors are parallel.
+    """
+    assert (len(v2) == len(v1)), 'Vectors lengths are not the same.'
+    TOL = 1.0e-6/2.0
+    if magnitude(crossproduct3(v1, v2)) < TOL * (magnitude(v1)+magnitude(v2)):
+        return True
+    return False

tests/test_colon.py

@@ -26,6 +26,7 @@ def test_colon1(self):
             'Test line': ScaffoldPackage(MeshType_1d_path1, {
                 'scaffoldSettings': {
                     'Coordinate dimensions': 3,
+                    'D2 derivatives':True,
                     'Length': 1.0,
                     'Number of elements': 1
                 },
@@ -83,7 +84,7 @@ def test_colon1(self):
 
         tmpRegion = region.createRegion()
         centralPath.generate(tmpRegion)
-        cx = extractPathParametersFromRegion(tmpRegion)[0]
+        cx = extractPathParametersFromRegion(tmpRegion, [Node.VALUE_LABEL_VALUE])[0]
         self.assertEqual(2, len(cx))
         assertAlmostEqualList(self, cx[0], [ 163.7, -25.2, 12.2 ], 1.0E-6)
         assertAlmostEqualList(self, cx[1], [ 117.2, 32.8, -2.6 ], 1.0E-6)

tests/test_smallintestine.py

@@ -24,6 +24,7 @@ def test_smallintestine1(self):
         centralPathDefaultScaffoldPackages = {
             'Test line': ScaffoldPackage(MeshType_1d_path1, {
                 'scaffoldSettings': {
+                    'D2 derivatives': True,
                     'Coordinate dimensions': 3,
                     'Length': 1.0,
                     'Number of elements': 3
@@ -62,7 +63,7 @@ def test_smallintestine1(self):
 
         tmpRegion = region.createRegion()
         centralPath.generate(tmpRegion)
-        cx = extractPathParametersFromRegion(tmpRegion)[0]
+        cx = extractPathParametersFromRegion(tmpRegion, [Node.VALUE_LABEL_VALUE])[0]
         self.assertEqual(4, len(cx))
         assertAlmostEqualList(self, cx[0], [  -2.3, 18.5, -4.4 ], 1.0E-6)
         assertAlmostEqualList(self, cx[2], [ -18.3, 12.6, -1.5 ], 1.0E-6)