Implement ellipse shape

[adam-urbanczyk]

cf. GC_MakeArcOfEllipse

Related to #260

[bernhard-42]

I also missed ellipses in one project and I am happy to give it a try.
Some first test are promising, however, I have an issue that the ellipse center is slightly translated if major_radius != minor-radius (which actually is the point about ellipses ...) I could further look into it if you are interested

[bernhard-42]

Found a possible reason for the shifted center for ellipses with major_radius > minor_radius:

from OCC.Core.gp import gp_Ax1, gp_Ax2, gp_Elips, gp_Pnt, gp_Dir
from OCC.Core.GC import GC_MakeArcOfCircle, GC_MakeArcOfEllipse
from OCC.Core.BRepBuilderAPI import BRepBuilderAPI_MakeEdge
from OCC.Core.GProp import GProp_GProps
from OCC.Core.BRepGProp import brepgprop_LinearProperties, brepgprop_SurfaceProperties

pnt = gp_Pnt(0.0, 0.0, 0.0)
dir = gp_Dir(0.0, 0.0, 1.0)

ellipse_gp = gp_Elips(gp_Ax2(pnt, dir), 4000, 40)
ellipse = BRepBuilderAPI_MakeEdge(ellipse_gp).Edge()

l_gp = ellipse_gp.Location()
print("ellipse_gp: %9.6f %9.6f %9.6f" % (l_gp.X(), l_gp.Y(), l_gp.Z()))

Properties = GProp_GProps()
brepgprop_SurfaceProperties(ellipse, Properties)
com = Properties.CentreOfMass()
print("ellipse (surface): %9.6f %9.6f %9.6f" % (com.X(), com.Y(), com.Z()))

Properties = GProp_GProps()
brepgprop_LinearProperties(ellipse, Properties)
com = Properties.CentreOfMass()
print("ellipse (linear): %9.6f %9.6f %9.6f" % (com.X(), com.Y(), com.Z()))

from cadquery.occ_impl.geom import BoundBox
bb = BoundBox._fromTopoDS(ellipse, tol=1e-8)
print("BoundingBox: %9.6f %9.6f %9.6f" % bb.center.toTuple())

results in

ellipse_gp:  0.000000  0.000000  0.000000
ellipse (surface):  0.000000  0.000000  0.000000
ellipse (linear): -81.295066  0.000000  0.000000
BoundingBox:  0.000001 -0.000009  0.000000

That means, the "linear center of mass" and the "surface center of mass" are not identical even though the edge is symmetrical along x-axis and y-axis.
The bounding box center is correct given a tolerance.

Shape.centerOfMass uses shape_properties_LUT[obj.wrapped.ShapeType()] which for ta.TopAbs_Edge returns brepgprop_LinearProperties and as such the 'shifted' center appears.

If I create a circle with ellipse_gp = gp_Elips(gp_Ax2(pnt, dir), 40, 40) then the result is:

ellipse_gp:  0.000000  0.000000  0.000000
ellipse (surface):  0.000000  0.000000  0.000000
ellipse (linear):  0.000000  0.000000  0.000000
BoundingBox:  0.000000  0.000000  0.000000

That's why it worked for circles.

I tested with OCC 7.4, same behavior.

Could it be that for closed edges one should use brepgprop_SurfacePropertiesto calculate the center of mass?

[bernhard-42]

Opened an issue in pythonocc tpaviot/pythonocc-core#747

[bernhard-42]

Opened an issue at OCCT mantis (https://tracker.dev.opencascade.org/view.php?id=31290)

Answer:

This seems to be a bug in calculation of linear properties of an ellipse.

Workaround is to create a face of that ellipse and to compute surface properties of that face.

The computed center of such face is (0, -3.91083e-15, 0).

@adam-urbanczyk @jmwright What shall we do?

1. Should we change Shape.centerOfMass to create a face first, when a closed edge is detected?

    @staticmethod
    def centerOfMass(obj):
        """
        Calculates the 'mass' of an object.
        """
        # start workaround
        # TODO remove as soon as center of mass for ellipses is fixed in OCC
        if obj.wrapped.ShapeType() == ta.TopAbs_WIRE and obj.wrapped.Closed():
            return Face.makeFromWires(obj).Center()
        # end workaround
  
        Properties = GProp_GProps()
        calc_function = shape_properties_LUT[obj.wrapped.ShapeType()]

        if calc_function:
            calc_function(obj.wrapped, Properties)
            return Vector(Properties.CentreOfMass())
        else:
            raise NotImplementedError

2. Edge.Center -> create wire and use Shape.Center

    def Center(self):

        # Workaround
        # TODO remove as soon as center of mass for ellipses is fixed in OCC
        # Properties = GProp_GProps()
        # brepgprop_LinearProperties(self.wrapped,
        #                            Properties)

        # return Vector(Properties.CentreOfMass())

        w = Wire.assembleEdges([self])
        return w.Center()

3. What about Shape._center_of_mass? It is defined but never used and it only uses brepgprop_VolumeProperties without ensuring that shape is actually a solid.

Or would you like to postpone ellipses until the bug is fixed in OCCT?

[bernhard-42]

Interesting findings:

1. For an ellipse approximation

def ellipse(maj_r, min_r, alpha):
    return (maj_r * math.cos(alpha), min_r * math.sin(alpha))

alphas = np.linspace(0, 2*math.pi, 3601)
points = [ellipse(400, 10, alpha) for alpha in alphas]

obj = cq.Workplane("XY").polyline(points).close()

the center of mass is correct. So it is not a matter of the shape dimensions or the shape form, but of the ellipse itself.

2. Using gp_Elips leads to extrem slow performance:

obj = cq.Workplane("XY").ellipse(400, 10)
show(obj)

takes 24sec in jupyter-cadquery.

[adam-urbanczyk]

@bernhard-42 I would not change center of mass for now. Note that center is mostly calculated on faces and not on edges.

[bernhard-42]

Good for me.
Given that the test whether the center is correct for an ellipse (similar to _make_circle) fails, are you interested in a PR with ellipse or shall I drop it?

[adam-urbanczyk]

PR is always welcome and the testing strategy can be adjusted.

[bernhard-42]

As an extension (doesn't work with the above PR) I could also support rotation_angle being an array. It needs to be same length as number of points on the stack. We then could do:

from math import sin, cos, pi
import numpy as np

r, r1, r2 = 80, 30, 10
angles = np.linspace(0, 360, 12, endpoint=False).tolist()
points = [(r * cos(angle), r * sin(angle)) for angle in np.deg2rad(angles)]
`
(cq.Workplane("XY")
     .circle(r+r1+10).extrude(5)
     .pushPoints(points)
     .ellipse(r1, r2, angle1=45, angle2=-45, closed=True, rotation_angle=angles)
     .cutThruAll()
)

This would involve adding a keyword argument argsArray to Workplane.each and Workplane.eachpoint. It can easily be done keeping the current functionality, but not sure, whether you want to change these functions

[adam-urbanczyk]

Thanks @bernhard-42 , there is already #118. I'd prefer to go for a generic solution and handle the rotations at the pushPoints level.

[adam-urbanczyk]

Solved by #265