-- File: SVMappingsImpl.mesa
-- Last edited by Bier on December 19, 1982 9:52 pm
-- Author: Eric Bier on October 19, 1982 1:09 pm
-- Contents: A set of bidirectional mappings between 3d surfaces and a 2d image plane

DIRECTORY
 Graphics,
 Matrix3d,
 RealFns,
 Rope,
 SV2d,
 SVAngle,
 SVArtwork,
 SVDraw,
 SVMappings,
 SVPolygon2d,
 SVVector2d,
 SVVector3d;

SVMappingsImpl: PROGRAM
IMPORTS RealFns, SVAngle, SVDraw, SVPolygon2d
EXPORTS SVMappings =
BEGIN

Point2d: TYPE = SV2d.Point2d;
Point3d: TYPE = Matrix3d.Point3d;
Polygon: TYPE = SVPolygon2d.Polygon;
Vector: TYPE = SVVector3d.Vector;
Vector2d: TYPE = SVVector2d.Vector2d;

Artwork: TYPE = REF ArtworkObj;
ArtworkObj: TYPE = SVArtwork.ArtworkObj;

PDisk: TYPE = REF PDiskObj;
PDiskObj: TYPE = SVMappings.PDiskObj;

Tube: TYPE = REF TubeObj;
TubeObj: TYPE = SVMappings.TubeObj;

Box: TYPE = REF BoxObj;
BoxObj: TYPE = SVMappings.BoxObj;

Disk: TYPE = REF DiskObj;
DiskObj: TYPE = SVMappings.DiskObj;

Pi: REAL = 3.14159;
PiUnder180: REAL = 180.0/Pi;
PiOver180: REAL = Pi/180.0;

CartesianToPolar: PUBLIC PROC [x, y: REAL] RETURNS [r, theta: REAL] = {
 r ← RealFns.SqRt[x*x + y*y];
 theta ← SVAngle.ArcTan[y,x];
 };

PolarToCartesian: PUBLIC PROC [r, theta: REAL] RETURNS [x, y: REAL] = {
 x ← r*RealFns.CosDeg[theta];
 y ← r*RealFns.SinDeg[theta];
 };

CreatePDisk: PUBLIC PROC [r: REAL, origin: Point2d] RETURNS [pdisk: PDisk] = {
 pdisk ← NEW[PDiskObj ← [r, origin]];
 };

DrawPDisk: PUBLIC PROC [dc: Graphics.Context, pdisk: PDisk] = {
 SVDraw.Circle[dc, pdisk.origin[1], pdisk.origin[2], pdisk.r];
 };

CreateTube: PUBLIC PROC [R: REAL, H: REAL] RETURNS [tube: Tube] = {
 tubePoly: Polygon ← SVPolygon2d.CreatePoly[4]; -- a poly of [theta, h] pairs.
 tubePoly ← SVPolygon2d.AddPolyPoint[tubePoly, [-180, H/2]];
 tubePoly ← SVPolygon2d.AddPolyPoint[tubePoly, [180, H/2]];
 tubePoly ← SVPolygon2d.AddPolyPoint[tubePoly, [180, -H/2]];
 tubePoly ← SVPolygon2d.AddPolyPoint[tubePoly, [-180, -H/2]];
 tube ← NEW[TubeObj ← [R, H, tubePoly]];
 };

CreateBox: PUBLIC PROC [x, y, z: REAL] RETURNS [box: Box] = {
 box ← NEW[BoxObj];
 box.x ← x;
 box.halfX ← x/2;
 box.y ← y;
 box.halfY ← y/2;
 box.threeHalvesY ← 3*box.halfY;
 box.z ← z;
 box.halfZ ← z/2.0;
FOR i: NAT IN [1..6] DO
  box.faces[i] ← SVPolygon2d.CreatePoly[4]; -- a poly of [x,y] pairs.
ENDLOOP;
-- The eight vertex points are named mnemonically. luf = left, up, front, rdb = right, down, back, etc. Initialize these now.

 box.faces[1] ← SVPolygon2d.AddPolyPoint[box.faces[1], [-(box.x+box.z), -box.halfY]];
 box.faces[1] ← SVPolygon2d.AddPolyPoint[box.faces[1], [-(box.x+box.z), box.halfY]];
 box.faces[1] ← SVPolygon2d.AddPolyPoint[box.faces[1], [-box.x, box.halfY]];
 box.faces[1] ← SVPolygon2d.AddPolyPoint[box.faces[1], [-box.x, -box.halfY]];

 box.faces[2] ← SVPolygon2d.AddPolyPoint[box.faces[2], [-box.x, -box.halfY]];
 box.faces[2] ← SVPolygon2d.AddPolyPoint[box.faces[2], [-box.x, box.halfY]];
 box.faces[2] ← SVPolygon2d.AddPolyPoint[box.faces[2], [0, box.halfY]];
 box.faces[2] ← SVPolygon2d.AddPolyPoint[box.faces[2], [0, -box.halfY]];

 box.faces[3] ← SVPolygon2d.AddPolyPoint[box.faces[3], [0, -box.halfY]];
 box.faces[3] ← SVPolygon2d.AddPolyPoint[box.faces[3], [0, box.halfY]];
 box.faces[3] ← SVPolygon2d.AddPolyPoint[box.faces[3], [box.z, box.halfY]];
 box.faces[3] ← SVPolygon2d.AddPolyPoint[box.faces[3], [box.z, -box.halfY]];

 box.faces[4] ← SVPolygon2d.AddPolyPoint[box.faces[4], [box.z, -box.halfY]];
 box.faces[4] ← SVPolygon2d.AddPolyPoint[box.faces[4], [box.z, box.halfY]];
 box.faces[4] ← SVPolygon2d.AddPolyPoint[box.faces[4], [(box.x+box.z), box.halfY]];
 box.faces[4] ← SVPolygon2d.AddPolyPoint[box.faces[4], [(box.x+box.z), -box.halfY]];

 box.faces[5] ← SVPolygon2d.AddPolyPoint[box.faces[5], [-box.x, box.halfY]];
 box.faces[5] ← SVPolygon2d.AddPolyPoint[box.faces[5], [-box.x, box.halfY + box.z]];
 box.faces[5] ← SVPolygon2d.AddPolyPoint[box.faces[5], [0, box.halfY + box.z]];
 box.faces[5] ← SVPolygon2d.AddPolyPoint[box.faces[5], [0, box.halfY]];

 box.faces[6] ← SVPolygon2d.AddPolyPoint[box.faces[6], [-box.x, -(box.halfY + box.z)]];
 box.faces[6] ← SVPolygon2d.AddPolyPoint[box.faces[6], [-box.x, -box.halfY]];
 box.faces[6] ← SVPolygon2d.AddPolyPoint[box.faces[6], [0, -box.halfY]];
 box.faces[6] ← SVPolygon2d.AddPolyPoint[box.faces[6], [0, -(box.halfY + box.z)]];

 };

Point3dToTubePoint: PUBLIC PROC [point3d: Point3d] RETURNS [tubePoint: Point2d] = {
-- given a point in space, we project it onto a vertical tube by dropping a perpendicular from the point to the y-axis. h = y and theta = ArcTan[-z/x] measuring angles in the xz plane as counter-clockwise (as seen from the positive y axis) from the positive x axis.
 tubePoint[2] ← point3d[2];
 tubePoint[1] ← SVAngle.ArcTan[-point3d[3], point3d[1]];
 };

TubePointToPoint3d: PUBLIC PROC [tubePoint: Point2d, assembly: REF ANY] RETURNS [point3d: Point3d] = {};
-- Given a point on a tube, we find a point on an object by casting a ray directly at the y-axis and finding out where it hits the object.


ImageToTube: PUBLIC PROC [imagePoint: Point2d, data: REF ANY] RETURNS [surfacePoint: Point2d, onSurf: BOOL] = {
-- imagePoint[2] just maps to h
-- imagePoint[1] maps to r*theta where r is the radius of the cylinder. ie theta ← u/r;
-- theta is surfacePoint[1]
 tube: Tube ← NARROW[data];
 surfacePoint[2] ← imagePoint[2];
 surfacePoint[1] ← (imagePoint[1]/tube.R)*PiUnder180;
 onSurf ← TRUE;
 };

TubeToImage: PUBLIC PROC [surfacePoint: Point2d, data: REF ANY] RETURNS [imagePoint: Point2d, onSurf: BOOL] = {
-- u ← (tube.R*theta)*PiOver180;
-- v ← h;
 tube: Tube ← NARROW[data];
 imagePoint[1] ← (tube.R*surfacePoint[1])*PiOver180;
 imagePoint[2] ← surfacePoint[2];
 onSurf ← TRUE;
 };

ImageToTopDisk: PUBLIC PROC [u,v: REAL, disk: Disk] RETURNS [r, theta: REAL] = {
-- Basically this is just a conversion to cartesian coordinates and a translation
 [r, theta] ← CartesianToPolar[u,v-disk.R-disk.halfH];
 };

TopDiskToImage: PUBLIC PROC [r, theta: REAL, disk: Disk] RETURNS [u,v: REAL] = {
 [u,v] ← PolarToCartesian[r, theta];
 v ← v + disk.R + disk.halfH;
 };

ImageDiskToTopDiskBoundary: PUBLIC PROC [pdisk: PDisk, disk: Disk] RETURNS [R, halfH: REAL] = {
R ← disk.R;
 halfH ← disk.halfH;
 };

TopDiskBoundaryToImageDisk: PUBLIC PROC [R, halfH: REAL, disk: Disk] RETURNS [pdisk: PDisk] = {
 pdisk ← NEW[PDiskObj ← [R, [0, halfH + R]]];
 };


ImageToBottomDisk: PUBLIC PROC [u,v: REAL, disk: Disk] RETURNS [r, theta: REAL] = {
-- Basically this is just a conversion to cartesian coordinates
 [r, theta] ← CartesianToPolar[u,v+disk.R+disk.halfH];
 };

BottomDiskToImage: PUBLIC PROC [r, theta: REAL, disk: Disk] RETURNS [u,v: REAL] = {
 [u,v] ← PolarToCartesian[r, theta];
 v ← v - disk.R - disk.halfH;
 };

ImageDiskToBottomDiskBoundary: PUBLIC PROC [pdisk: PDisk, disk: Disk] RETURNS [R, halfH: REAL] = {
R ← disk.R;
 halfH ← disk.halfH;
 };

BottomDiskBoundaryToImageDisk: PUBLIC PROC [R, halfH: REAL, disk: Disk] RETURNS [pdisk: PDisk] = {
 pdisk ← NEW[PDiskObj ← [R, [0, -halfH -R]]];
 };
Point3dToDiskPoint: PUBLIC PROC [point3d: Point3d] RETURNS [r, theta: REAL] = {
-- given a point in space, we project it onto a horizontal disk by dropping a perpendicular to the disk. This point of intersection differs from point3d in y value only. Hence, this mapping is just a conversion from Cartesian to polar coordinates.
  [r, theta] ← CartesianToPolar[point3d[1], point3d[2]];
  };

DiskPointToPoint3d: PUBLIC PROC [r, theta: REAL, assembly: REF ANY] RETURNS [point3d: Point3d] = {};


-- Box Mapping
-- Imagine an unfolded box as follows:
-- There are four rectangles numbered left to right as 1..4.
-- Two more rectangles numbered 5 and 6 are above and below rectangle 2 respectively.
-- Rectangle 2 maps to the front (pos z) side of the box. Square 6 maps to the bottom (neg y) side. All other rectangles map accordingly.
-- An [x,y] pair maps to a surface point as follows:
-- If a box has height sy, width sx, and depth sz then the correspondences are:
-- IF -sy/2 <= y <= sy/2 AND
-- -(sz+sx) <= x < -sx THEN Rectangle 1.
-- -sx <= x <= 0 THEN Rectangle 2.
-- 0 < x <= sz THEN Rectangle 3.
-- sz < x <= sz + sx THEN Rectangle 4.
-- IF -3sy/2 <= y < -sy/2 AND -sx <= x <= 0 THEN Rectangle 6.
-- IF sy/2 < y <= 3sy/2 AND -sx <= x <= 0 THEN Rectangle 5.
-- Otherwise, out of box.

IndicesOfMax: PRIVATE PROC [v: Vector] RETURNS [indices: ARRAY[1..3] OF NAT, negArray: ARRAY[1..3] OF BOOL] = {
-- bubble sort
 temp: REAL;
 tempIndex: NAT;
IF v[1] < 0 THEN {negArray[1] ← TRUE; v[1] ← -v[1]}
  ELSE negArray[1] ← FALSE;
IF v[2] < 0 THEN {negArray[2] ← TRUE; v[2] ← -v[2]}
  ELSE negArray[2] ← FALSE;
IF v[3] < 0 THEN {negArray[3] ← TRUE; v[3] ← -v[3]}
  ELSE negArray[3] ← FALSE;
 indices ← [1,2,3];
IF v[1] < v[2] THEN
  {temp ← v[1]; v[1] ← v[2]; v[2] ← temp;
   tempIndex ← indices[1]; indices[1] ← indices[2]; indices[2] ← tempIndex};
IF v[2] < v[3] THEN {temp ← v[2]; v[2] ← v[3]; v[3] ← temp;
   tempIndex ← indices[2]; indices[2] ← indices[3]; indices[3] ← tempIndex};
IF v[1] < v[2] THEN {temp ← v[1]; v[1] ← v[2]; v[2] ← temp;
   tempIndex ← indices[1]; indices[1] ← indices[2]; indices[2] ← tempIndex};
 };

BoxOrthogonalOInverse: PUBLIC PROC [point3d: Point3d, normal: Vector, box: Box] RETURNS [boxPoint: Point2d, hit: BOOL] = {
-- given a point in space on some surface, we project it onto an aligned box by finding the aligned direction which the normal most closely follows and dropping a surface normal from the surface of the box which this normal points to, to the point. We then decode this point in three space to an unfolded box point (See the comment in SVMappings on boxes).
 indices: ARRAY[1..3] OF NAT;
 negArray: ARRAY[1..3] OF BOOL;
 [indices, negArray] ← IndicesOfMax[normal];
-- indices is a permutation. 3,1,2 for instance means that normal is most closely aligned with the z axis, then the x axis, then the yaxis. negArray[indices[1]] tells whether it is most closely aligned with the positive or negative part of the axis. Together they define the 6 sides of a box as required.
-- in each case, make sure that the calculated box point is within the expected rectangle. This will not be the case in general if the simple surface box is smaller than the object being painted.
SELECT indices[1] FROM
 1 => IF negArray[1] THEN {-- rectangle 1 of box = left side
   -- y and z are of interest
   IF point3d[2] < -box.halfY OR point3d[2] > box.halfY OR point3d[3] < -box.halfZ OR point3d[3] > box.halfZ THEN RETURN [[0,0], FALSE];
   boxPoint[2] ← point3d[2];
   boxPoint[1] ← -(box.x+box.z) + point3d[3] + box.halfZ;
   }
   ELSE {-- rectangle 3 of box = right side
   -- y and z are of interest
   IF point3d[2] < -box.halfY OR point3d[2] > box.halfY OR point3d[3] < -box.halfZ OR point3d[3] > box.halfZ THEN RETURN [[0,0], FALSE];
   boxPoint[2] ← point3d[2];
   boxPoint[1] ← -point3d[3] + box.halfZ;
   };
 2 => IF negArray[2] THEN {-- rectangle 6 of box = down side
   -- x and z are of interest
   IF point3d[1] < -box.halfX OR point3d[1] > box.halfX OR point3d[3] < -box.halfZ OR point3d[3] > box.halfZ THEN RETURN [[0,0], FALSE];
   boxPoint[1] ← point3d[1] + -box.halfX;
   boxPoint[2] ← point3d[3] -box.halfY -box.halfZ;
   }
   ELSE {-- rectangle 5 of box = up side
   IF point3d[1] < -box.halfX OR point3d[1] > box.halfX OR point3d[3] < -box.halfZ OR point3d[3] > box.halfZ THEN RETURN [[0,0], FALSE];
   boxPoint[1] ← point3d[1] + -box.halfX;
   boxPoint[2] ← -point3d[3] + box.halfY + box.halfZ;
   };
 3 => IF negArray[3] THEN {-- rectangle 4 of box = back
   -- x and y are of interest
   IF point3d[1] < -box.halfX OR point3d[1] > box.halfX OR point3d[2] < -box.halfY OR point3d[2] > box.halfY THEN RETURN [[0,0], FALSE];
   boxPoint[2] ← point3d[2];
   boxPoint[1] ← box.z - point3d[1] + box.halfX;
   }
   ELSE {-- rectangle 2 of box = front
   IF point3d[1] < -box.halfX OR point3d[1] > box.halfX OR point3d[2] < -box.halfY OR point3d[2] > box.halfY THEN RETURN [[0,0], FALSE];
   boxPoint[2] ← point3d[2];
   boxPoint[1] ← point3d[1] - box.halfX;
   };
ENDCASE => ERROR;
 hit ← TRUE;
 };

BoxPointToUnfoldedBoxPoint: PROC [boxPoint: Point3d, rect: NAT, pos: BOOL, box: Box] RETURNS [unfoldedBoxPoint: Point2d, hit: BOOL] = {
SELECT rect FROM
 1 => IF NOT pos THEN {-- rectangle 1 of box = left side
   -- y and z are of interest
   IF boxPoint[2] < -box.halfY OR boxPoint[2] > box.halfY OR boxPoint[3] < -box.halfZ OR boxPoint[3] > box.halfZ THEN RETURN [[0,0], FALSE];
   unfoldedBoxPoint[2] ← boxPoint[2];
   unfoldedBoxPoint[1] ← -(box.x+box.z) + boxPoint[3] + box.halfZ;
   }
   ELSE {-- rectangle 3 of box = right side
   -- y and z are of interest
   IF boxPoint[2] < -box.halfY OR boxPoint[2] > box.halfY OR boxPoint[3] < -box.halfZ OR boxPoint[3] > box.halfZ THEN RETURN [[0,0], FALSE];
   unfoldedBoxPoint[2] ← boxPoint[2];
   unfoldedBoxPoint[1] ← -boxPoint[3] + box.halfZ;
   };
 2 => IF NOT pos THEN {-- rectangle 6 of box = down side
   -- x and z are of interest
   IF boxPoint[1] < -box.halfX OR boxPoint[1] > box.halfX OR boxPoint[3] < -box.halfZ OR boxPoint[3] > box.halfZ THEN RETURN [[0,0], FALSE];
   unfoldedBoxPoint[1] ← boxPoint[1] + -box.halfX;
   unfoldedBoxPoint[2] ← boxPoint[3] -box.halfY -box.halfZ;
   }
   ELSE {-- rectangle 5 of box = up side
   IF boxPoint[1] < -box.halfX OR boxPoint[1] > box.halfX OR boxPoint[3] < -box.halfZ OR boxPoint[3] > box.halfZ THEN RETURN [[0,0], FALSE];
   unfoldedBoxPoint[1] ← boxPoint[1] + -box.halfX;
   unfoldedBoxPoint[2] ← -boxPoint[3] + box.halfY + box.halfZ;
   };
 3 => IF NOT pos THEN {-- rectangle 4 of box = back
   -- x and y are of interest
   IF boxPoint[1] < -box.halfX OR boxPoint[1] > box.halfX OR boxPoint[2] < -box.halfY OR boxPoint[2] > box.halfY THEN RETURN [[0,0], FALSE];
   unfoldedBoxPoint[2] ← boxPoint[2];
   unfoldedBoxPoint[1] ← box.z - boxPoint[1] + box.halfX;
   }
   ELSE {-- rectangle 2 of box = front
   IF boxPoint[1] < -box.halfX OR boxPoint[1] > box.halfX OR boxPoint[2] < -box.halfY OR boxPoint[2] > box.halfY THEN RETURN [[0,0], FALSE];
   unfoldedBoxPoint[2] ← boxPoint[2];
   unfoldedBoxPoint[1] ← boxPoint[1] - box.halfX;
   };
ENDCASE => ERROR;
 hit ← TRUE;
 };


VectorHitsSide: PROC [v: Vector2d, width, height: REAL] RETURNS [vertOrHoriz: NAT, positive: BOOL] = {
-- given a two dimensional vector v whose tail is on the origin and a rectangle of dimensions width by height centered on the origin and aligned with the axes, which side of the rectangle does v point to?
 yOverX: REAL;
IF v[1] = 0 THEN {vertOrHoriz ← 2; positive ← v[2] > 0; RETURN};
 yOverX ← v[2]/v[1];
IF v[1] < 0 THEN {
  IF v[2] > 0 THEN {
   IF -height/width <= yOverX AND yOverX <= 0 THEN RETURN[1, FALSE]
   ELSE RETURN[2, TRUE];
   }
  ELSE {
   IF 0 <= yOverX AND yOverX <= height/width THEN RETURN[1, FALSE]
   ELSE RETURN[2, FALSE];
   }
  }
ELSE {
  IF v[2] > 0 THEN {
   IF 0 <= yOverX AND yOverX <= height/width THEN RETURN[1, TRUE]
   ELSE RETURN[2, TRUE];
   }
  ELSE {
   IF -height/width <= yOverX AND yOverX <= 0 THEN RETURN[1, TRUE]
   ELSE RETURN[2, FALSE];
   }
  };
 }; -- end of VectorHitsSide

RayHitsRect: PROC [v: Vector, rect: NAT, positive: BOOL, box: Box] RETURNS [point3d: Point3d] = {
-- rect will be in [1..3]. Along with positive we describe a face of the box. For rect = 1, postive = FALSE for example, we describe the negative x face (face 1 in SVMappings). Its plane equation is x = -box.x. v, with its tail at the origin will forms a ray R(t) = t*v. R hits the plane when t*v[1] = -box.x. t = -box.x/v[1]. This is the point [-box.x, v[2]*(-box.x)/v[1], v[3]*(-box.x)/v[1]]. Similarly for other planes.
 k, t: REAL;
SELECT rect FROM
 1 => {k ← box.halfX; IF NOT positive THEN k ← -k;
  IF v[1] = 0 THEN ERROR; -- this situation should be filtered out by VectorHitsSide.
  t ← k/v[1];
  point3d ← [k, v[2]*t, v[3]*t];
  };
 2 => {k ← box.halfY; IF NOT positive THEN k ← -k;
  IF v[2] = 0 THEN ERROR; -- this situation should be filtered out by VectorHitsSide.
  t ← k/v[2];
  point3d ← [v[1]*t, k, v[3]*t];
  };
 3 => {k ← box.halfZ; IF NOT positive THEN k ← -k;
  IF v[3] = 0 THEN ERROR; -- this situation should be filtered out by VectorHitsSide.
  t ← k/v[3];
  point3d ← [v[1]*t, v[2]*t, k];
  };
ENDCASE;
};
  
BoxRadialOInverse: PUBLIC PROC [point3d: Point3d, normal: Vector, box: Box] RETURNS [boxPoint: Point2d, hit: BOOL] = {
-- a mapping from point3d to a box point along a ray passing from the origin of the box through point3d. We must figure out which face of the box is hit and then perform a ray-plane intersection.
-- Project the vector from origin to point3d onto xy and zy planes. Find out which rectangle it is nearest in both cases and infer which face it points at.
 xyVert, zyVert, zxVert: NAT;
 xyPos, zyPos, zxPos: BOOL;
 boxPoint3d: Point3d;
 xScore, yScore, zScore: NAT;
 xScore ← yScore ← zScore ← 0;
 [xyVert, xyPos] ← VectorHitsSide[[point3d[1], point3d[2]], box.x, box.z];
IF xyVert = 1 THEN xScore ← xScore + 1 ELSE yScore ← yScore + 1;
 [zyVert, zyPos] ← VectorHitsSide[[point3d[3], point3d[2]], box.z, box.y];
IF zyVert = 1 THEN zScore ← zScore + 1 ELSE yScore ← yScore + 1;
 [zxVert, zxPos] ← VectorHitsSide[[point3d[3], point3d[1]], box.z, box.x];
IF zxVert = 1 THEN zScore ← zScore + 1 ELSE xScore ← xScore + 1;
-- who wins 2 out of three?
SELECT TRUE FROM
 xScore = 2 => {boxPoint3d ← RayHitsRect[point3d, 1, xyPos, box];
     [boxPoint, hit] ← BoxPointToUnfoldedBoxPoint[boxPoint3d, 1, xyPos, box];
     };
 yScore = 2 => {boxPoint3d ← RayHitsRect[point3d, 2, xyPos, box];
      [boxPoint, hit] ← BoxPointToUnfoldedBoxPoint[boxPoint3d, 2, xyPos, box];
     };
 zScore = 2 => {boxPoint3d ← RayHitsRect[point3d, 3, zxPos, box];
     [boxPoint, hit] ← BoxPointToUnfoldedBoxPoint[boxPoint3d, 3, zxPos, box];
     };
ENDCASE => ERROR;
 };
    

ImageToBox: PUBLIC PROC [imagePoint: Point2d, data: REF ANY] RETURNS [surfacePoint: Point2d, onSurf: BOOL] = {
-- data will be a Box. It is a ref any so we can use this mapping in ImagePolyToSurfacePoly below
-- surfacePoint is a [x, y] pair (See comment on Box Mapping above)
-- since I use unfolded coordinates to refer to points on the surface of a box, this mapping is the identity. However, we can check for illegal values.
 box: Box ← NARROW[data];
IF -box.halfY <= imagePoint[2] AND imagePoint[2] <= box.halfY
  AND -(box.x+box.z) <= imagePoint[1] AND imagePoint[1] <= (box.x+box.z)
THEN onSurf ← TRUE
ELSE IF -(box.halfY+box.z) <= imagePoint[2] AND imagePoint[2] <= -box.halfY
   THEN onSurf ← TRUE
   ELSE IF box.halfY <= imagePoint[2] AND imagePoint[2] <= box.halfY+box.z
      THEN onSurf ← TRUE
      ELSE onSurf ← FALSE;
 surfacePoint ← imagePoint;
 }; -- end of ImageToBox

BoxToImage: PUBLIC PROC [surfacePoint: Point2d, data: REF ANY] RETURNS [imagePoint: Point2d, onSurf: BOOL] = {
-- data will be a Box. It is a ref any so we can use this mapping in SurfacePolyToImagePoly below
-- surfacePoint is a [x, y] pair (See comment on Box Mapping above)
-- since I use unfolded coordinates to refer to points on the surface of a box, this mapping is the identity. However, we can check for illegal values.
 box: Box ← NARROW[data];
IF -box.halfY <= surfacePoint[2] AND surfacePoint[2] <= box.halfY
  AND -(box.x+box.z) <= surfacePoint[1] AND surfacePoint[1] <= (box.x+box.z)
THEN onSurf ← TRUE
ELSE IF -(box.halfY+box.z) <= surfacePoint[2] AND surfacePoint[2] <= -box.halfY
   THEN onSurf ← TRUE
   ELSE IF box.halfY <= surfacePoint[2] AND surfacePoint[2] <= box.halfY + box.z
      THEN onSurf ← TRUE
      ELSE onSurf ← FALSE;
 imagePoint ← surfacePoint;
 }; -- end of BoxToImage

ImagePolyToSurfacePoly: PUBLIC PROC [imagePoly: Polygon, map: PROC [imagePoint: Point2d, data: REF ANY] RETURNS [surfacePoint: Point2d, onSurf: BOOL], data: REF ANY] RETURNS [surfacePoly: Polygon] = {
 surfacePoint: Point2d;
 onSurf: BOOL;
-- data will probably be a REF <Surface> for some type of 3d surface.
 surfacePoly ← SVPolygon2d.CreatePoly[imagePoly.len];
FOR i: NAT IN [0..imagePoly.len) DO
  [surfacePoint, onSurf] ← map[imagePoly[i], data];
  surfacePoly ← SVPolygon2d.AddPolyPoint[surfacePoly, surfacePoint];
ENDLOOP;
 };

SurfacePolyToImagePoly: PUBLIC PROC [surfacePoly: Polygon, map: PROC [surfacePoint: Point2d, data: REF ANY] RETURNS [imagePoint: Point2d, onSurf: BOOL], data: REF ANY] RETURNS [imagePoly: Polygon] = {
 imagePoint: Point2d;
 onSurf: BOOL;
-- data will probably be a REF <Surface> for some type of 3d surface
 imagePoly ← SVPolygon2d.CreatePoly[surfacePoly.len];
FOR i: NAT IN [0..surfacePoly.len) DO
  [imagePoint, onSurf] ← map[surfacePoly[i], data];
  imagePoly ← SVPolygon2d.AddPolyPoint[imagePoly, imagePoint];
ENDLOOP;
 };

END.