How to project a 2D region from one plane to another?

Question

I have one of a general set of 2D regions (Disk, Rectangle, Polygon, etc.), as well as 3D origin points and normal vectors defining two reference planes/coordinate systems. I want to find a series of 2D transformations that I can apply using TransformedRegion[] to project a region defined in one plane to the other.

I was able to fairly easily project "down" from the first plane to the second, where the projected region is always smaller than the original. This is akin to shining a flashlight through the object in the first plane when the beam is normal to the second plane. However, what I need to do is find the projection as if the beam is normal to the first plane, where the "shadow" of the region on the second plane can become arbitrarily large. The function below works for basic axial transformations and rotations of each plane, but the result distorts and rotates substantially when the two plane normal vectors are skewed from each other.

FindProjectedRegion[region_, origin1_, normal1_, origin2_, normal2_] :=
  Module[{rot, trans, normal2ref, origin2ref, rmat, tmat},
  (*Get transformations to go from global coordinates to local \
coordinates for plane 1*)
  rot = RotationTransform[{normal1, {0, 0, 1}}];
  trans = TranslationTransform[-origin1];
  (*Get the normal vector and origin of plane 2 in the coordinates of \
plane 1*)
  normal2ref = rot[normal2];
  origin2ref = rot[trans[origin2]];
  (*Construct 2D transformations to project the region.*)
  rmat = TransformationFunction[
    Inverse[Drop[
      RotationTransform[{{0, 0, 1}, normal2ref}][[1]], {3}, {3}]]];
  tmat = TransformationFunction[
    Drop[TranslationTransform[-origin2ref][[1]], {3}, {3}]];
  TransformedRegion[TransformedRegion[region, tmat], rmat]]

The following Manipulate shows what I want this to do, giving control over the parameters for the x,y,z coordinates of the origin of each plane as well as the two spherical angles for their normal vectors. The rectangle region is projected correctly from the red plane 1 to the blue plane 2, except for when the phi angle of planes 1 and 2 are different.

LocalToGlobal[obj_, n_, p_] := 
  GeometricTransformation[
   GeometricTransformation[obj, RotationTransform[{{0, 0, 1}, n}]], 
   TranslationTransform[p]];
Manipulate[Block[{n1, n2, p1, p2, xl, reg},
  n1 = {Cos[ph1] Sin[th1], Sin[ph1] Sin[th1], Cos[th1]};
  n2 = {Cos[ph2] Sin[th2], Sin[ph2] Sin[th2], Cos[th2]};
  p1 = {px1, py1, pz1}; p2 = {px2, py2, pz2};
  reg = FindProjectedRegion[Rectangle[{-0.5, -1}, {0.5, 1}], p1, n1, 
    p2, n2];
  Graphics3D[
   (*Projection of the region as an extended cube*)
   {{Purple, Opacity[0.1], 
     LocalToGlobal[Cuboid[{-0.5, -1, -20}, {0.5, 1, 20}], n1, p1],
     (*Planes 1 and 2*)
     Red, 
     LocalToGlobal[Cuboid[{-20, -20, -0.001}, {20, 20, 0.001}], n1, 
      p1],
     Blue, 
     LocalToGlobal[Cuboid[{-20, -20, -0.001}, {20, 20, 0.001}], n2, 
      p2]},
    (* Plane plot windowed by the projected region *)
    {Opacity[1], 
     LocalToGlobal[
      Plot3D[0, {x, y} \[Element] reg, 
        PlotRange -> {{-1, 1}, {-1, 1}, {-3, 3}}, ImageSize -> 250, 
        Mesh -> False][[1]], n2, p2]}},
   ViewPoint -> {-2, 0.5, 2},
   ViewVertical -> {0, 1, 0},
   Boxed -> False,
   AxesLabel -> Table[Style[i, Bold, Black, 13], {i, {"x", "y", "z"}}],
    Axes -> True,
   PlotRange -> {{-3, 3}, {-3, 3}, {-3, 3}}, ImageSize -> 300]
  ], {{th1, 0.0}, 0, Pi}, {{ph1, Pi/2}, 0, 2*Pi}, {{th2, Pi/4}, 0, 
  Pi}, {{ph2, Pi/2}, 0, 2*Pi}, {{px1, 0}, -1, 1}, {{py1, 0}, -1, 
  1}, {{pz1, 0}, -1, 1}, {{px2, 0}, -1, 1}, {{py2, 0}, -1, 
  1}, {{pz2, 0}, -1, 1}, ControlPlacement -> Left]

Answer 1

I assume you already figured how to go from 2D in a plane to 3D and back.

There remains the problem of projection of a 3D point: pos along some unit direction: dir onto a plane through the origin specified by a unit normal: normal.

First, dir and norm must not be perpendicular. That is dir.normal != 0

Then, the 3D projection: proj of pos onto the plane can be written as p plus some vector in the direction of dir:

proj= pos + lam dir

the parameter lam is obtained from proj.normal==0;

(pos + lam dir).normal==0
lam== - pos.normal /dir.normal
proj== pos- pos.normal /dir.normal dir

We may create a function to do this that also takes care of normalization:

projection[pos_, dir0_, normal0_] := 
  Module[{dir = Normalize[dir0], normal = Normalize[normal0]},
   If[dir . normal == 0, Print["dir and normal are perpendicular"]; 
    Return[]];
   pos - pos . normal/dir . normal dir
   ];

As a test we project a point {0,0,1} onto the plane through the origin with normal {-1,0,1}.:

projection[{0, 0, 1}, {1, 0, 0}, {-1, 0, 1}]
(*{1, 0, 1}*)