#246 breaks `face_normals(coordinates(mesh), faces(mesh))`

[KeitaNakamura]

Hi, thank you for providing a great package.
After updating to v0.5.6 (specifically due to #246), the behavior of face_normals(coordinates(mesh), faces(mesh)) changed, as this no longer calls orthogonal_vector for Triangle. I’m wondering if this might be unintended behavior?

[ffreyer]

The Triangle method is only a fast path. Not hitting that should not result in different normals.

#245 updated orthogonal_vector for larger faces, which could cause different result but those should be a better approximation of the face normal than before.

Do you have an example for when it changes?

[KeitaNakamura]

Thank you for your reply. I didn't check it in detail, but for example, this STL file, which is used on Wikipedia, generates some normal vectors with NaN values:

julia> using FileIO, MeshIO, LinearAlgebra

julia> mesh = FileIO.load("Stanford Bunny -- Digitized! - 88208/files/Bunny.stl")
Mesh{3, Float32, GeometryBasics.NgonFace{3, GeometryBasics.OffsetInteger{-1, UInt32}}}
    faces: 112402
    vertex position: 337206
    vertex normal: 337206

julia> normals = values(face_normals(coordinates(mesh), faces(mesh)));

julia> face_indices = findall(n->any(!isfinite,n), normals)
37-element Vector{Int64}:
   3563
   5514
   8240
   9483
   9738
  10255
  12977
  14917
  25745
  26774
  32270
  40324
  40572
  41620
  42117
  44864
  45078
  45158
  45222
  45229
  46717
  46796
  46948
  50143
  52202
  63950
  65056
  65910
  76008
  81433
  85390
  92728
  92745
  95490
  95642
  99348
 108959

julia> normals[face_indices[1]]
3-element Vec{3, Float32} with indices SOneTo(3):
 NaN
 NaN
 NaN

julia> tri = mesh[face_indices[1]]
Triangle(Float32[-10.180012, 9.928894, 63.42825], Float32[-10.180012, 9.92804, 63.424137], Float32[-10.181978, 9.92804, 63.42825])

julia> normalize(cross(tri[2]-tri[1], tri[3]-tri[1]))
3-element Point{3, Float32} with indices SOneTo(3):
 -0.39148504
  0.9009399
 -0.18720828

julia> normals_from_tri = map(tri -> GeometryBasics.orthogonal_vector(Vec3f, tri), mesh); # check for all triangles

julia> findall(n->any(!isfinite,n), normals_from_tri)
Int64[]

[ffreyer]

Hmm, I'm not sure if this is a lucky case for the fast path or the new summation being unstable. Mathematically it's the same as far as I can tell. I guess effectively it's swapping the order of a-b and cross() which could be a precision problem?

In your example the direction vectors contain 0s which happen to 0 out the subtraction in the cross product. That's why I'm wondering if it's just a lucky case.

t = Triangle(
    Point3f(-10.180012, 9.928894, 63.42825), 
    Point3f(-10.180012, 9.92804, 63.424137), 
    Point3f(-10.181978, 9.92804, 63.42825)
)
a, b, c = t
julia> a-b
3-element Point{3, Float32} with indices SOneTo(3):
 0.0
 0.0008544922
 0.0041122437

julia> c-b
3-element Point{3, Float32} with indices SOneTo(3):
 -0.0019664764
  0.0
  0.0041122437