[OpenSCAD] Reimplement surface() to match OpenSCAD

The original implementation of the surface() function used a simple
B-spline representation for the surface, which generated degenerate
surface with several of OpenSCAD's demo input files. This commit
modifies the algorithm to generate a discrete surface identical to the
one generated within OpenSCAD itself. It also adds several units tests
to identify future regressions.

Note that PNG input is not yet supported for the surface() function.

src/Mod/OpenSCAD/CMakeLists.txt

@@ -43,6 +43,9 @@ SET(OpenSCADTestsFiles_SRCS
     OpenSCADTest/data/CSG.csg
     OpenSCADTest/data/Cube.stl
     OpenSCADTest/data/Square.dxf
+    OpenSCADTest/data/Surface.dat
+    OpenSCADTest/data/Surface2.dat
+    OpenSCADTest/data/Surface.png
 )
 
 SET(OpenSCADTests_ALL

src/Mod/OpenSCAD/OpenSCADFeatures.py

@@ -502,33 +502,139 @@ def execute(self,fp):
             raise ValueError
 
 def makeSurfaceVolume(filename):
-    import FreeCAD,Part
-    f1=open(filename)
-    coords=[]
-    miny=1
-    for line in f1.readlines():
-        sline=line.strip()
-        if sline and not sline.startswith('#'):
-            ycoord=len(coords)
-            lcoords=[]
-            for xcoord, num in enumerate(sline.split()):
-                fnum=float(num)
-                lcoords.append(FreeCAD.Vector(float(xcoord),float(ycoord),fnum))
-                miny=min(fnum,miny)
-            coords.append(lcoords)
-    s=Part.BSplineSurface()
-    s.interpolate(coords)
-    plane=Part.makePlane(len(coords[0])-1,len(coords)-1,FreeCAD.Vector(0,0,miny-1))
-    l1=Part.makeLine(plane.Vertexes[0].Point,s.value(0,0))
-    l2=Part.makeLine(plane.Vertexes[1].Point,s.value(1,0))
-    l3=Part.makeLine(plane.Vertexes[2].Point,s.value(0,1))
-    l4=Part.makeLine(plane.Vertexes[3].Point,s.value(1,1))
-    f0=plane.Faces[0]
-    f0.reverse()
-    f1=Part.Face(Part.Wire([plane.Edges[0],l1.Edges[0],s.vIso(0).toShape(),l2.Edges[0]]))
-    f2=Part.Face(Part.Wire([plane.Edges[1],l3.Edges[0],s.uIso(0).toShape(),l1.Edges[0]]))
-    f3=Part.Face(Part.Wire([plane.Edges[2],l4.Edges[0],s.vIso(1).toShape(),l3.Edges[0]]))
-    f4=Part.Face(Part.Wire([plane.Edges[3],l2.Edges[0],s.uIso(1).toShape(),l4.Edges[0]]))
-    f5=s.toShape().Faces[0]
-    solid=Part.Solid(Part.Shell([f0,f1,f2,f3,f4,f5]))
-    return solid,(len(coords[0])-1)/2.0,(len(coords)-1)/2.0
+    import FreeCAD,Part,sys
+    with open(filename) as f1:
+        coords = []
+        min_z = sys.float_info.max
+        for line in f1.readlines():
+            sline=line.strip()
+            if sline and not sline.startswith('#'):
+                ycoord=len(coords)
+                lcoords=[]
+                for xcoord, num in enumerate(sline.split()):
+                    fnum=float(num)
+                    lcoords.append(FreeCAD.Vector(float(xcoord),float(ycoord),fnum))
+                    min_z = min(fnum,min_z)
+                coords.append(lcoords)
+
+        num_rows = len(coords)
+        num_cols = len(coords[0])
+
+        # OpenSCAD does not spline this surface, so neither do we: just create a bunch of faces, 
+        # using four triangles per quadrilateral
+        faces = []
+        for row in range(num_rows-1):
+            for col in range(num_cols-1):
+                a = coords[row+0][col+0]
+                b = coords[row+0][col+1]
+                c = coords[row+1][col+1]
+                d = coords[row+1][col+0]
+                centroid = 0.25 * (a + b + c + d)
+                ab = Part.makeLine(a,b)
+                bc = Part.makeLine(b,c)
+                cd = Part.makeLine(c,d)
+                da = Part.makeLine(d,a)
+
+                diag_a = Part.makeLine(a, centroid)
+                diag_b = Part.makeLine(b, centroid)
+                diag_c = Part.makeLine(c, centroid)
+                diag_d = Part.makeLine(d, centroid)
+
+                wire1 = Part.Wire([ab,diag_a,diag_b])
+                wire2 = Part.Wire([bc,diag_b,diag_c])
+                wire3 = Part.Wire([cd,diag_c,diag_d])
+                wire4 = Part.Wire([da,diag_d,diag_a])
+
+                try:
+                    face = Part.Face(wire1)
+                    faces.append(face)
+                    face = Part.Face(wire2)
+                    faces.append(face)
+                    face = Part.Face(wire3)
+                    faces.append(face)
+                    face = Part.Face(wire4)
+                    faces.append(face)
+                except Exception:
+                    print ("Failed to create the face from {},{},{},{}".format(coords[row+0][col+0],\
+                        coords[row+0][col+1],coords[row+1][col+1],coords[row+1][col+0]))
+
+        last_row = num_rows-1
+        last_col = num_cols-1
+
+        # Create the face to close off the y-min border: OpenSCAD places the lower surface of the shell
+        # at 1 unit below the lowest coordinate in the surface
+        lines = []
+        corner1 = FreeCAD.Vector(coords[0][0].x, coords[0][0].y, min_z-1)
+        lines.append (Part.makeLine(corner1,coords[0][0]))
+        for col in range(num_cols-1):
+            a = coords[0][col]
+            b = coords[0][col+1]
+            lines.append (Part.makeLine(a, b))
+        corner2 = FreeCAD.Vector(coords[0][last_col].x, coords[0][last_col].y, min_z-1)
+        lines.append (Part.makeLine(corner2,coords[0][last_col]))
+        lines.append (Part.makeLine(corner1,corner2))
+        wire = Part.Wire(lines)
+        face = Part.Face(wire)
+        faces.append(face)
+
+        # Create the face to close off the y-max border
+        lines = []
+        corner1 = FreeCAD.Vector(coords[last_row][0].x, coords[last_row][0].y, min_z-1)
+        lines.append (Part.makeLine(corner1,coords[last_row][0]))
+        for col in range(num_cols-1):
+            a = coords[last_row][col]
+            b = coords[last_row][col+1]
+            lines.append (Part.makeLine(a, b))
+        corner2 = FreeCAD.Vector(coords[last_row][last_col].x, coords[last_row][last_col].y, min_z-1)
+        lines.append (Part.makeLine(corner2,coords[last_row][last_col]))
+        lines.append (Part.makeLine(corner1,corner2))
+        wire = Part.Wire(lines)
+        face = Part.Face(wire)
+        faces.append(face)
+
+         # Create the face to close off the x-min border
+        lines = []
+        corner1 = FreeCAD.Vector(coords[0][0].x, coords[0][0].y, min_z-1)
+        lines.append (Part.makeLine(corner1,coords[0][0]))
+        for row in range(num_rows-1):
+            a = coords[row][0]
+            b = coords[row+1][0]
+            lines.append (Part.makeLine(a, b))
+        corner2 = FreeCAD.Vector(coords[last_row][0].x, coords[last_row][0].y, min_z-1)
+        lines.append (Part.makeLine(corner2,coords[last_row][0]))
+        lines.append (Part.makeLine(corner1,corner2))
+        wire = Part.Wire(lines)
+        face = Part.Face(wire)
+        faces.append(face)
+
+         # Create the face to close off the x-max border
+        lines = []
+        corner1 = FreeCAD.Vector(coords[0][last_col].x, coords[0][last_col].y, min_z-1)
+        lines.append (Part.makeLine(corner1,coords[0][last_col]))
+        for row in range(num_rows-1):
+            a = coords[row][last_col]
+            b = coords[row+1][last_col]
+            lines.append (Part.makeLine(a, b))
+        corner2 = FreeCAD.Vector(coords[last_row][last_col].x, coords[last_row][last_col].y, min_z-1)
+        lines.append (Part.makeLine(corner2,coords[last_row][last_col]))
+        lines.append (Part.makeLine(corner1,corner2))
+        wire = Part.Wire(lines)
+        face = Part.Face(wire)
+        faces.append(face)
+
+        # Create a bottom surface to close off the shell
+        a = FreeCAD.Vector(coords[0][0].x, coords[0][0].y, min_z-1)
+        b = FreeCAD.Vector(coords[0][last_col].x, coords[0][last_col].y, min_z-1)
+        c = FreeCAD.Vector(coords[last_row][last_col].x, coords[last_row][last_col].y, min_z-1)
+        d = FreeCAD.Vector(coords[last_row][0].x, coords[last_row][0].y, min_z-1)
+        ab = Part.makeLine(a,b)
+        bc = Part.makeLine(b,c)
+        cd = Part.makeLine(c,d)
+        da = Part.makeLine(d,a)
+        wire = Part.Wire([ab,bc,cd,da])
+        face = Part.Face(wire)
+        faces.append(face)
+
+        s = Part.Shell(faces)
+        solid = Part.Solid(s)
+        return solid,last_col,last_row

src/Mod/OpenSCAD/OpenSCADTest/app/test_importCSG.py

@@ -313,7 +313,38 @@ def test_import_resize(self):
         FreeCAD.closeDocument(doc.Name)
 
     def test_import_surface(self):
-        pass
+        testfile = join(self.test_dir, "Surface.dat").replace('\\','/')
+        doc = self.utility_create_scad(f"surface(file = \"{testfile}\", center = true, convexity = 5);", "surface_simple_dat")
+        object = doc.ActiveObject
+        self.assertTrue (object is not None)
+        self.assertAlmostEqual (object.Shape.Volume, 275.000000, 6)
+        self.assertAlmostEqual (object.Shape.BoundBox.XMin, -4.5, 6)
+        self.assertAlmostEqual (object.Shape.BoundBox.XMax, 4.5, 6)
+        self.assertAlmostEqual (object.Shape.BoundBox.YMin, -4.5, 6)
+        self.assertAlmostEqual (object.Shape.BoundBox.YMax, 4.5, 6)
+        FreeCAD.closeDocument(doc.Name)
+
+        testfile = join(self.test_dir, "Surface.dat").replace('\\','/')
+        doc = self.utility_create_scad(f"surface(file = \"{testfile}\", convexity = 5);", "surface_uncentered_dat")
+        object = doc.ActiveObject
+        self.assertTrue (object is not None)
+        self.assertAlmostEqual (object.Shape.Volume, 275.000000, 6)
+        self.assertAlmostEqual (object.Shape.BoundBox.XMin, 0, 6)
+        self.assertAlmostEqual (object.Shape.BoundBox.XMax, 9, 6)
+        self.assertAlmostEqual (object.Shape.BoundBox.YMin, 0, 6)
+        self.assertAlmostEqual (object.Shape.BoundBox.YMax, 9, 6)
+        FreeCAD.closeDocument(doc.Name)
+
+        testfile = join(self.test_dir, "Surface2.dat").replace('\\','/')
+        doc = self.utility_create_scad(f"surface(file = \"{testfile}\", center = true, convexity = 5);", "surface_rectangular_dat")
+        object = doc.ActiveObject
+        self.assertTrue (object is not None)
+        self.assertAlmostEqual (object.Shape.Volume, 24.5500000, 6)
+        self.assertAlmostEqual (object.Shape.BoundBox.XMin, -2, 6)
+        self.assertAlmostEqual (object.Shape.BoundBox.XMax, 2, 6)
+        self.assertAlmostEqual (object.Shape.BoundBox.YMin, -1.5, 6)
+        self.assertAlmostEqual (object.Shape.BoundBox.YMax, 1.5, 6)
+        FreeCAD.closeDocument(doc.Name)
 
     def test_import_projection(self):
         pass

src/Mod/OpenSCAD/OpenSCADTest/data/Surface.dat

@@ -0,0 +1,10 @@
+10 9 8 7 6 5 5 5 5 5 
+9 8 7 6 6 4 3 2 1 0 
+8 7 6 6 4 3 2 1 0 0
+7 6 6 4 3 2 1 0 0 0
+6 6 4 3 2 1 1 0 0 0
+6 6 3 2 1 1 1 0 0 0
+6 6 2 1 1 1 1 0 0 0
+6 6 1 0 0 0 0 0 0 0
+3 1 0 0 0 0 0 0 0 0
+3 0 0 0 0 0 0 0 0 0

src/Mod/OpenSCAD/OpenSCADTest/data/Surface2.dat

@@ -0,0 +1,5 @@
+# Example comment
+1.0 1.1 1.3 1.7 2.5
+1.1 1.3 1.7 2.5 2.7
+1.3 1.7 2.5 2.7 2.9
+1.7 2.5 2.7 2.9 3.0