Finding surface normal for 3D region at a specific point

Question

I would like to find the surface normal for a point on a 3D filled shape in Mathematica.

I know how to calculate the normal of a parametric surface using the cross product but this method will not work for a shape like Cone[] or Ball[].

1. Is there some sort of RegionNormal option? There is an option to find VertexNormals here, but this is something to with shading and seems unhelpful.
2. Is there a method I can use to convert the region into a parametric expression and use the normal cross product method?

The plan is to take an arbitrary line and find the angle of intersection between the line and the surface of the shape.

Answer 1

(* put inequality into u ≤ 0 form, return just u *)

standardize[a_ <= b_] := a - b;
standardize[a_ >= b_] := b - a;
regnormal[reg_, {x_, y_, z_}] := Module[{impl},
   impl = LogicalExpand@ Simplify[RegionMember[reg, {x, y, z}], {x, y, z} ∈ Reals];
   If[Head@impl === Or,
    impl = List @@ impl,
    impl = List@impl];
   impl = Replace[impl, {Verbatim[And][a___] :> {a}, e_ :> {e}}, 1];
   Piecewise[
    Flatten[
     Function[{component},
       Table[{
         D[standardize[component[[i]]], {{x, y, z}}],
         Simplify[
          (And @@ Drop[component, {i}] /. {LessEqual -> Less, GreaterEqual -> Greater}) &&
           (component[[i]] /. {LessEqual -> Equal, GreaterEqual -> Equal}),
          TransformationFunctions -> {Automatic, 
            Reduce[#, {}, Reals] &}]
         }, {i, Length@component}]
       ] /@ impl,
     1],
    Indeterminate]
   ];

Examples:

regnormal[Cone[{{0, 0, 0}, {1, 1, 1}}, 1/2], {x, y, z}]

regnormal[Ball[{1, 2, 3}, 4], {x, y, z}]

regnormal[RegionUnion[Ball[], Cone[{{0, 0, 0}, {1, 1, 1}}, 1/2]], {x, y, z}]

regnormal[Cylinder[{{1, 1, 1}, {2, 3, 1}}], {x, y, z}]

It assumes that the RegionMember expression can be computed (which is not always the case) and that it will be a union (via Or) of intersections (via And). It also assumes that the RegionMember expression includes the boundary. Thus, it is probably not very robust, but it handles the OP's examples.

Also, if this is used in numerical applications, which seems to be the case for the OP, one should worry about the exact conditions in the Piecewise expressions returned. It's unlikely the numerical calculations will be accurate enough to satisfy Equal. So either change the conditions or possible change Internal`$EqualTolerance:

Block[{Internal`$EqualTolerance = Log10[2.^28]}, (* ~single-precision FP equality *)
 <evaluate regnormal[...] expression>
 ]

Answer 2

Wether this is enough to warrant an "answer" is debatable and it relies on MichaelE2's work, but I felt it was helpful to share.

Using Michael E2's solution, we can plot and clearly see the normals for 3D shapes.

numberofpoints = 50; 
pts = RandomPoint[RegionBoundary[shape], numberofpoints];
normals = regnormal[shape, {x, y, z}]
surfaces = Length[normals[[1]]];
magnitude = 0.1;

pointonsurface = ConstantArray[0, surfaces];

lines[{a_, b_}] :=  {magnitude*normals[[1, a, 1]] + b, 
    b} /. {x -> b[[1]], y ->  b[[2]], z -> b[[3]]};

For[i = 1, i <= surfaces, i++, 

 pointonsurface[[i]] = 
  Table[ {If[
     Evaluate[
       normals[[1, i, 2]] /. {x -> pts[[j, 1]], y ->  pts[[j, 2]], 
         z -> pts[[j, 3]]}] == True, i], pts[[j]]}, {j, 1, 
    numberofpoints, 1}];

 pointonsurface[[i]]  = 
  DeleteCases[pointonsurface[[i]], {a_, b_} /; a == Null];
 pointonsurface[[i]] = lines /@ pointonsurface[[i]];

 ]

normallines = Flatten[pointonsurface, 1];

Graphics3D[{shape, {Red, Point[pts]}, Line[normallines] }, 
 Boxed -> False, Axes -> True, PlotRange -> All]

So, for Cone[]

And for Cuboid[]

And for Sphere[]

Answer 3

If you know the equation of your surface, you can do this. For instance a circle at the origin of radius 5.

f[x_,y_,z_]=x^2+y^2+z^2-5^2;

grad[x_,y_,z_]=Grad[f[x,y,z],{x,y,z}]

The unit normal vector on that surface is

normal[x_,y_,z_]=Simplify[grad[x,y,z]/Sqrt[grad[x,y,z].grad[x,y,z]]]

Simple case

normal[0,0,5]
(*{0,0,1}*)

normal[2,3,2 Sqrt[3]]
(*{2/5,3/5,(2 Sqrt[3])/5}*)

The hard part is getting the surface into the equation.

Answer 4

I just came across the UnitNormal[] function that might be applicable: https://resources.wolframcloud.com/FunctionRepository/resources/UnitNormal

Also, SurfaceData[surface,"NormalVector"] might work: https://reference.wolfram.com/language/ref/SurfaceData.html