How to inspect and reverse face normals?

[gymclean]

I am trying to inspect the normals of some faces with pythonocc-core (Version 7.8.1) using

from OCC.Core.BRep import BRep_Tool
from OCC.Core.GeomLProp import GeomLProp_SLProps

def get_normal(face, u, v, order=1, resolution=1e-6):
    brep_surf = BRep_Tool.Surface(face)
    props = GeomLProp_SLProps(brep_surf, u, v, order, resolution)
    return props.Normal()

and to reverse the normal of a face am trying face.Reverse(), but I get the same result for the normal before and after face.Reverse() so I expect I am not properly understanding how to inspect or manipulate face normals. My main question is "What is the correct way to inspect and reverse face normals?"

As a test case I am using a prism-shaped body constructed in FreeCAD 1.0.0 and exported to STEP using the AP214 International Standard scheme.

On a side note, I find that the face normals obtained in the FreeCAD console by selecting faces and calling face.normalAt() are as I expected -- that is, facing outward -- but in the STEP some of the normals point inward (per my method of inspection, at least) which is what led me to the question of how to reverse normals in the first place. Should I expect the normals in STEP to match those reported from within FreeCAD? Is there a way to make it so?

A script detailing my efforts follows. The output I get from the script (which contains the STEP text) is

$ python test_normals.py 
The normals I expect to find:
   (0,0,-1) for the square base
   (+sqrt(3)/2, 0, 0.5) = (+0.866, 0, 0.5) for the square face in the +x half-space
   (-sqrt(3)/2, 0, 0.5) = (-0.866, 0, 0.5) for the square face in the -x half-space
   (0, +1, 0) for the triangular face located at y = +5 mm.
   (0, -1, 0) for the triangular face located at y = -5 mm.
Face normals using BRep_Tool.Surface(face) and GeomLProp_SLProps:
  (0.000, 0.000, 1.000) [@ uv=(0.000, 0.000)]
  (-0.866, -0.000, -0.500) [@ uv=(0.000, 0.000)]
  (0.866, -0.000, -0.500) [@ uv=(0.000, 0.000)]
  (0.000, 1.000, -0.000) [@ uv=(0.000, 0.000)]
  (0.000, 1.000, -0.000) [@ uv=(0.000, 0.000)]

Face normals using BRep_Adaptor.Surface(face) and GeomLProp_SLProps:
  (0.000, 0.000, 1.000) [@ uv=(0.000, 0.000)]
  (-0.866, -0.000, -0.500) [@ uv=(0.000, 0.000)]
  (0.866, -0.000, -0.500) [@ uv=(0.000, 0.000)]
  (0.000, 1.000, -0.000) [@ uv=(0.000, 0.000)]
  (0.000, 1.000, -0.000) [@ uv=(0.000, 0.000)]

Checking normal of base before and after reversing the face.
  (0.000, 1.000, -0.000) [@ uv=(0.000, 0.000)] (before)
  (0.000, 1.000, -0.000) [@ uv=(0.000, 0.000)] (after)

FreeCAD construction

Test script:

"""Inspect face normals of a single prism-shaped body created in FreeCAD.

Compare face normals as reported by FreeCAD and pythonocc.

Try reversing the normal of a face.
"""

import os.path

from OCC.Core.BRep import BRep_Tool
from OCC.Core.BRepAdaptor import BRepAdaptor_Surface
from OCC.Core.BRepTools import breptools
from OCC.Core.GeomAbs import GeomAbs_Plane
from OCC.Core.GeomLProp import GeomLProp_SLProps
from OCC.Core.STEPControl import STEPControl_Reader
from OCC.Core.TopAbs import TopAbs_FACE
from OCC.Core.TopExp import TopExp_Explorer

def test(stepfile):
    """Check normals in STEP file

    Args:
        stepfile: a filename or STEP content as a string
    """

    print("""The normals I expect to find:
   (0,0,-1) for the square base
   (+sqrt(3)/2, 0, 0.5) = (+0.866, 0, 0.5) for the square face in the +x half-space
   (-sqrt(3)/2, 0, 0.5) = (-0.866, 0, 0.5) for the square face in the -x half-space
   (0, +1, 0) for the triangular face located at y = +5 mm.
   (0, -1, 0) for the triangular face located at y = -5 mm.""")

    # Read the shape from STEP.
    step_reader = STEPControl_Reader()
    if os.path.isfile(stepfile):
        step_reader.ReadFile(stepfile)
    else:
        step_reader.ReadStream("STEP file", stepfile)
    step_reader.TransferRoot()
    shape = step_reader.Shape()

    # Get the faces. There should be five.
    faces = []
    explorer = TopExp_Explorer(shape, TopAbs_FACE)
    while explorer.More():
        faces.append(explorer.Current())
        explorer.Next()

    # Get the normal direction for each face using OCCT.

    def get_normal(face, u, v, order=1, resolution=1e-6):
        brep_surf = BRep_Tool.Surface(face)
        props = GeomLProp_SLProps(brep_surf, u, v, order, resolution)
        return props.Normal()

    def get_uv_midpoint(face):
        umin, umax, vmin, vmax = breptools.UVBounds(face)
        u = (umin + umax) / 2
        v = (vmin + vmax) / 2
        return u, v

    #get_eval_point = get_uv_midpoint
    get_eval_point = lambda face: (0.0, 0.0)

    print(f"Face normals using BRep_Tool.Surface(face) and GeomLProp_SLProps:")
    for face in faces:
        # Get uv coordinates at which to evaluate the normal.
        u, v = get_eval_point(face)

        # Verify that the coordinates are within the uv bounds of the face.
        umin, umax, vmin, vmax = breptools.UVBounds(face)
        assert (umin <= u <= umax)
        assert (vmin <= v <= vmax)

        # Evaluate and report the face normal.
        uvstr = f"({u:0.3f}, {v:0.3f})"
        normal = get_normal(face, u, v)
        nx, ny, nz = normal.X(), normal.Y(), normal.Z()
        print(f"  ({nx:0.3f}, {ny:0.3f}, {nz:0.3f}) [@ uv={uvstr}]")
    print()

    print(f"Face normals using BRep_Adaptor.Surface(face) and GeomLProp_SLProps:")
    for face in faces:
        surf = BRepAdaptor_Surface(face, True)
        surf_type = surf.GetType()
        if surf_type == GeomAbs_Plane:
            gp_pln = surf.Plane()
            location = gp_pln.Location()
            normal = gp_pln.Axis().Direction() 
            nx, ny, nz = normal.X(), normal.Y(), normal.Z()
            print(f"  ({nx:0.3f}, {ny:0.3f}, {nz:0.3f}) [@ uv={uvstr}]")
        else:
            print("Skipping non-planar face.")
    print()
            
    # Try flipping the normal of the first face (the base of the prism.)
    print("Checking normal of base before and after reversing the face.")
    base = faces[0]
    n_before = get_normal(face, 0, 0)
    nx, ny, nz = n_before.X(), n_before.Y(), n_before.Z()
    print(f"  ({nx:0.3f}, {ny:0.3f}, {nz:0.3f}) [@ uv={uvstr}] (before)")
    base.Reverse()
    n_after = get_normal(face, 0, 0)
    nx, ny, nz = n_after.X(), n_after.Y(), n_after.Z()
    print(f"  ({nx:0.3f}, {ny:0.3f}, {nz:0.3f}) [@ uv={uvstr}] (after)")

# Below is the content of a STEP file containing a prism. It was constructed
# by drawing a sketch of an equilateral triangle with sides of length 10 mm in
# the X-Z plane with the base horizontal and the apex on the z axis at z=10 mm
# and by padding the triangle symmetrically about the X-Z plane by 5 mm in
# each direction (+y, -y). The base of the prism is a square whose sides
# have length 10 mm, and whose center lies at x=0, y=0, and
#   z = [1 - sqrt(3)/2] * (10 mm) = 1.340 mm.
# If the face normals all point outward they should have values
#   (0,0,-1) for the square base
#   (+sqrt(3)/2, 0, 0.5) = (+0.866, 0, 0.5) for the square face in the +x half-space
#   (-sqrt(3)/2, 0, 0.5) = (-0.866, 0, 0.5) for the square face in the -x half-space
#   (0, +1, 0) for the triangular face located at y = +5 mm.
#   (0, -1, 0) for the triangular face located at y = -5 mm.
# The prism body was exported to STEP format from FreeCAD 1.0.0 using the
# AP214 International Standard with "Export single object placement" enabled.

STEPFILE = """ISO-10303-21;
HEADER;
FILE_DESCRIPTION(('FreeCAD Model'),'2;1');
FILE_NAME('Open CASCADE Shape Model','2025-04-21T08:35:16',(''),(''),
  'Open CASCADE STEP processor 7.7','FreeCAD','Unknown');
FILE_SCHEMA(('AUTOMOTIVE_DESIGN { 1 0 10303 214 1 1 1 1 }'));
ENDSEC;
DATA;
#1 = APPLICATION_PROTOCOL_DEFINITION('international standard',
  'automotive_design',2000,#2);
#2 = APPLICATION_CONTEXT(
  'core data for automotive mechanical design processes');
#3 = SHAPE_DEFINITION_REPRESENTATION(#4,#10);
#4 = PRODUCT_DEFINITION_SHAPE('','',#5);
#5 = PRODUCT_DEFINITION('design','',#6,#9);
#6 = PRODUCT_DEFINITION_FORMATION('','',#7);
#7 = PRODUCT('prism','prism','',(#8));
#8 = PRODUCT_CONTEXT('',#2,'mechanical');
#9 = PRODUCT_DEFINITION_CONTEXT('part definition',#2,'design');
#10 = ADVANCED_BREP_SHAPE_REPRESENTATION('',(#11,#15),#132);
#11 = AXIS2_PLACEMENT_3D('',#12,#13,#14);
#12 = CARTESIAN_POINT('',(0.,0.,0.));
#13 = DIRECTION('',(0.,0.,1.));
#14 = DIRECTION('',(1.,0.,-0.));
#15 = MANIFOLD_SOLID_BREP('',#16);
#16 = CLOSED_SHELL('',(#17,#57,#88,#110,#121));
#17 = ADVANCED_FACE('',(#18),#52,.F.);
#18 = FACE_BOUND('',#19,.F.);
#19 = EDGE_LOOP('',(#20,#30,#38,#46));
#20 = ORIENTED_EDGE('',*,*,#21,.T.);
#21 = EDGE_CURVE('',#22,#24,#26,.T.);
#22 = VERTEX_POINT('',#23);
#23 = CARTESIAN_POINT('',(-5.,5.,1.339745962156));
#24 = VERTEX_POINT('',#25);
#25 = CARTESIAN_POINT('',(-5.,-5.,1.339745962156));
#26 = LINE('',#27,#28);
#27 = CARTESIAN_POINT('',(-5.,5.,1.339745962156));
#28 = VECTOR('',#29,1.);
#29 = DIRECTION('',(0.,-1.,2.22044604925E-16));
#30 = ORIENTED_EDGE('',*,*,#31,.T.);
#31 = EDGE_CURVE('',#24,#32,#34,.T.);
#32 = VERTEX_POINT('',#33);
#33 = CARTESIAN_POINT('',(5.,-5.,1.339745962156));
#34 = LINE('',#35,#36);
#35 = CARTESIAN_POINT('',(-5.,-5.,1.339745962156));
#36 = VECTOR('',#37,1.);
#37 = DIRECTION('',(1.,0.,0.));
#38 = ORIENTED_EDGE('',*,*,#39,.F.);
#39 = EDGE_CURVE('',#40,#32,#42,.T.);
#40 = VERTEX_POINT('',#41);
#41 = CARTESIAN_POINT('',(5.,5.,1.339745962156));
#42 = LINE('',#43,#44);
#43 = CARTESIAN_POINT('',(5.,5.,1.339745962156));
#44 = VECTOR('',#45,1.);
#45 = DIRECTION('',(0.,-1.,2.22044604925E-16));
#46 = ORIENTED_EDGE('',*,*,#47,.F.);
#47 = EDGE_CURVE('',#22,#40,#48,.T.);
#48 = LINE('',#49,#50);
#49 = CARTESIAN_POINT('',(-5.,5.,1.339745962156));
#50 = VECTOR('',#51,1.);
#51 = DIRECTION('',(1.,0.,0.));
#52 = PLANE('',#53);
#53 = AXIS2_PLACEMENT_3D('',#54,#55,#56);
#54 = CARTESIAN_POINT('',(-5.,5.,1.339745962156));
#55 = DIRECTION('',(0.,2.22044604925E-16,1.));
#56 = DIRECTION('',(1.,0.,0.));
#57 = ADVANCED_FACE('',(#58),#83,.F.);
#58 = FACE_BOUND('',#59,.F.);
#59 = EDGE_LOOP('',(#60,#61,#69,#77));
#60 = ORIENTED_EDGE('',*,*,#39,.T.);
#61 = ORIENTED_EDGE('',*,*,#62,.T.);
#62 = EDGE_CURVE('',#32,#63,#65,.T.);
#63 = VERTEX_POINT('',#64);
#64 = CARTESIAN_POINT('',(0.,-5.,10.));
#65 = LINE('',#66,#67);
#66 = CARTESIAN_POINT('',(5.,-5.,1.339745962156));
#67 = VECTOR('',#68,1.);
#68 = DIRECTION('',(-0.5,1.922962686384E-16,0.866025403784));
#69 = ORIENTED_EDGE('',*,*,#70,.F.);
#70 = EDGE_CURVE('',#71,#63,#73,.T.);
#71 = VERTEX_POINT('',#72);
#72 = CARTESIAN_POINT('',(0.,5.,10.));
#73 = LINE('',#74,#75);
#74 = CARTESIAN_POINT('',(0.,5.,10.));
#75 = VECTOR('',#76,1.);
#76 = DIRECTION('',(0.,-1.,2.22044604925E-16));
#77 = ORIENTED_EDGE('',*,*,#78,.F.);
#78 = EDGE_CURVE('',#40,#71,#79,.T.);
#79 = LINE('',#80,#81);
#80 = CARTESIAN_POINT('',(5.,5.,1.339745962156));
#81 = VECTOR('',#82,1.);
#82 = DIRECTION('',(-0.5,1.922962686384E-16,0.866025403784));
#83 = PLANE('',#84);
#84 = AXIS2_PLACEMENT_3D('',#85,#86,#87);
#85 = CARTESIAN_POINT('',(5.,5.,1.339745962156));
#86 = DIRECTION('',(-0.866025403784,-1.110223024625E-16,-0.5));
#87 = DIRECTION('',(-0.5,1.922962686384E-16,0.866025403784));
#88 = ADVANCED_FACE('',(#89),#105,.F.);
#89 = FACE_BOUND('',#90,.F.);
#90 = EDGE_LOOP('',(#91,#92,#98,#99));
#91 = ORIENTED_EDGE('',*,*,#70,.T.);
#92 = ORIENTED_EDGE('',*,*,#93,.T.);
#93 = EDGE_CURVE('',#63,#24,#94,.T.);
#94 = LINE('',#95,#96);
#95 = CARTESIAN_POINT('',(0.,-5.,10.));
#96 = VECTOR('',#97,1.);
#97 = DIRECTION('',(-0.5,-1.922962686384E-16,-0.866025403784));
#98 = ORIENTED_EDGE('',*,*,#21,.F.);
#99 = ORIENTED_EDGE('',*,*,#100,.F.);
#100 = EDGE_CURVE('',#71,#22,#101,.T.);
#101 = LINE('',#102,#103);
#102 = CARTESIAN_POINT('',(0.,5.,10.));
#103 = VECTOR('',#104,1.);
#104 = DIRECTION('',(-0.5,-1.922962686384E-16,-0.866025403784));
#105 = PLANE('',#106);
#106 = AXIS2_PLACEMENT_3D('',#107,#108,#109);
#107 = CARTESIAN_POINT('',(0.,5.,10.));
#108 = DIRECTION('',(0.866025403784,-1.110223024625E-16,-0.5));
#109 = DIRECTION('',(-0.5,-1.922962686384E-16,-0.866025403784));
#110 = ADVANCED_FACE('',(#111),#116,.T.);
#111 = FACE_BOUND('',#112,.F.);
#112 = EDGE_LOOP('',(#113,#114,#115));
#113 = ORIENTED_EDGE('',*,*,#47,.T.);
#114 = ORIENTED_EDGE('',*,*,#78,.T.);
#115 = ORIENTED_EDGE('',*,*,#100,.T.);
#116 = PLANE('',#117);
#117 = AXIS2_PLACEMENT_3D('',#118,#119,#120);
#118 = CARTESIAN_POINT('',(7.105427357601E-16,5.,4.226497308104));
#119 = DIRECTION('',(0.,1.,-4.930380657631E-32));
#120 = DIRECTION('',(0.,4.930380657631E-32,1.));
#121 = ADVANCED_FACE('',(#122),#127,.F.);
#122 = FACE_BOUND('',#123,.T.);
#123 = EDGE_LOOP('',(#124,#125,#126));
#124 = ORIENTED_EDGE('',*,*,#31,.T.);
#125 = ORIENTED_EDGE('',*,*,#62,.T.);
#126 = ORIENTED_EDGE('',*,*,#93,.T.);
#127 = PLANE('',#128);
#128 = AXIS2_PLACEMENT_3D('',#129,#130,#131);
#129 = CARTESIAN_POINT('',(7.105427357601E-16,-5.,4.226497308104));
#130 = DIRECTION('',(0.,1.,-4.930380657631E-32));
#131 = DIRECTION('',(0.,4.930380657631E-32,1.));
#132 = ( GEOMETRIC_REPRESENTATION_CONTEXT(3) 
GLOBAL_UNCERTAINTY_ASSIGNED_CONTEXT((#136)) GLOBAL_UNIT_ASSIGNED_CONTEXT
((#133,#134,#135)) REPRESENTATION_CONTEXT('Context #1',
  '3D Context with UNIT and UNCERTAINTY') );
#133 = ( LENGTH_UNIT() NAMED_UNIT(*) SI_UNIT(.MILLI.,.METRE.) );
#134 = ( NAMED_UNIT(*) PLANE_ANGLE_UNIT() SI_UNIT($,.RADIAN.) );
#135 = ( NAMED_UNIT(*) SI_UNIT($,.STERADIAN.) SOLID_ANGLE_UNIT() );
#136 = UNCERTAINTY_MEASURE_WITH_UNIT(LENGTH_MEASURE(1.E-07),#133,
  'distance_accuracy_value','confusion accuracy');
#137 = PRODUCT_RELATED_PRODUCT_CATEGORY('part',$,(#7));
#138 = MECHANICAL_DESIGN_GEOMETRIC_PRESENTATION_REPRESENTATION('',(#139)
  ,#132);
#139 = STYLED_ITEM('color',(#140),#15);
#140 = PRESENTATION_STYLE_ASSIGNMENT((#141,#147));
#141 = SURFACE_STYLE_USAGE(.BOTH.,#142);
#142 = SURFACE_SIDE_STYLE('',(#143));
#143 = SURFACE_STYLE_FILL_AREA(#144);
#144 = FILL_AREA_STYLE('',(#145));
#145 = FILL_AREA_STYLE_COLOUR('',#146);
#146 = COLOUR_RGB('',0.800000010877,0.800000010877,0.90196104065);
#147 = CURVE_STYLE('',#148,POSITIVE_LENGTH_MEASURE(0.1),#149);
#148 = DRAUGHTING_PRE_DEFINED_CURVE_FONT('continuous');
#149 = COLOUR_RGB('',9.803921802644E-02,9.803921802644E-02,
  9.803921802644E-02);
ENDSEC;
END-ISO-10303-21;
"""

def write_step_file(path):
    with open(path, 'w') as fh:
        fh.write(STEPFILE)

#fn = "prism.step"
#write_step_file(fn)
#test(fn)

test(STEPFILE)

[gymclean]

It seems that if I reverse the normal of each face (or not) according to whether its face.Orientation() is reversed, I get the results I was after.

from OCC.Core.BRep import BRep_Tool
from OCC.Core.GeomLProp import GeomLProp_SLProps
from OCC.Core.TopAbs import TopAbs_REVERSED

def get_normal(face, u, v, order=1, resolution=1e-6):
    brep_surf = BRep_Tool.Surface(face)
    props = GeomLProp_SLProps(brep_surf, u, v, order, resolution)
    if face.Orientation() == TopAbs_REVERSED: 
        return props.Normal().Reversed()
    else:
        return props.Normal()