//  view.c
//
//  These routines transform three-dimensional objects into two-dimensional screen images.
//  The view may be either a plane projection or a perspective from a vantage point.
//
//  by John D. de Boer

#include "common.h"
#include "cstring.h"
#include "cyber.h"
#include "disk.h"
#include "floating.h"
#include "frame.h"
#include "graph.h"
#include "vector.h"
#include "view.h"


//  Allocation

  projection *newprojection(real X, real Y, real Z, real Mag, rect R) {
projection *P; //initzoom();
P=pointer(sizeof(projection)); if (!P) return NULL;
setprojangles(P,0,0); P->persp=0; P->terr=0; P->str=0; P->stretch=1.0;
setvector(&(P->look),X,Y,Z); P->mag=Mag; P->maglim=0;
P->rads= 40.0 / fiftyseven; P->fov=anglestring(P->rads,0);
P->port=R; resizeproj(P,R);
P->zout=NULL; P->zmin=NULL; P->shlf=NULL; P->shup=NULL;
P->shbk=NULL; P->shfw=NULL; P->shdn=NULL; P->shrt=NULL;
return P; }

  void oldprojection(projection *P) {
if (!P) return; if (P->fov) old(P->fov); old(P); }


//  Setting projection properties

//  Note:  The line of nodes between world coordinates and screen coordinates is a scanline.  (Or at least it was
//  on a CRT.)  With the change from vertical screen coordinate increasing downward with QuickDraw to it increasing
//  upward with Quartz 2D, the screen coordinate system is now left handed.  The "real-world" coordinates are
//  x=away, y=leftward, z=up when facing zero azimuth and zero elevation.  The screen coordinates are now x=right,
//  y=up, z=in.  That has led to the requirements for many procedures in this file to contain a statement which
//  inverts the vertical screen coordinate.  In the long run it might make more sense to rotate so that it is
//  the depth coordinate that has to be inverted at the end.

  static void setrotation(projection *P) {
//matrix *R; real SA,CA,SE,CE;
if (!P) return;
eulermatrix(P->azim + pibytwo,pibytwo - P->elev,pi,&(P->rot));
return;
/*
R=&(P->rot);
SA=sin(P->azim); CA=cos(P->azim); SE=sin(P->elev); CE=cos(P->elev);
R->xx=      SA; R->xy=     -CA; R->xz= 0.0;  //  projection to horizontal screen coordinate
R->yx= CA * SE; R->yy= SA * SE; R->yz= -CE;  //  projection to vertical screen coordinate
R->zx= CA * CE; R->zy= SA * CE; R->zz=  SE;  //  projection to depth coordinate
//  These used to be coded directly for reason of speed.  I can probably ditch this now.
*/
}

  void resizeproj(projection *P, rect R) {
if (!P) return; P->port=R; P->width=width(R);
P->org=middle(R);
P->pix= P->width / P->rads; }

  void projmode(projection *P, bool Persp) {
if (!P) return; P->persp=Persp; }

  void setprojangles(projection *P, ControlRef Azim, ControlRef Elev) {
int AZ,EL; if (!P) return;
if (Azim) { AZ=GetControlValue(Azim); if (!P->terr) AZ= 360 - AZ;
  P->azim= AZ / fiftyseven; }
if (Elev) { EL=GetControlValue(Elev); if (P->terr) EL-=90; else EL= 90 - EL;
  P->elev= EL / fiftyseven; }
setrotation(P); }

  void setscrollangles(projection *P, ControlRef Azim, ControlRef Elev) {
real Az,El; if (!P) return;
Az= P->terr ? P->azim : (twopi - P->azim);
El= P->terr ? (pibytwo + P->elev) : (pibytwo - P->elev);
if (Azim) SetControlValue(Azim,mod(Az * fiftyseven,360));
if (Elev) SetControlValue(Elev,between(El * fiftyseven,0,180));
setrotation(P); }

  void setazimelev(projection *P, real Azim, real Elev) {
P->elev=fbetween(Elev,-pibytwo,pibytwo); P->azim=modtwopisym(Azim); setrotation(P); }

  void setviewangles(projection *P, real Theta, real Phi) {
setazimelev(P,Phi,pibytwo - Theta); }

  void setdirection(projection *P, const vector *To) {
vector V; real Theta; if (!P || !To) return;
vectordifference(To,&(P->look),&V);
cartesiantopolar(&V,NULL,&Theta,&(P->azim)); P->elev= pibytwo - Theta;
setrotation(P); }

  void setlook(projection *P, const vector *V) {
if (!P || !V) return; P->look=*V; }

  void setlookpoint(projection *P, real X, real Y, real Z) {
if (!P) return; setvector(&(P->look),X,Y,Z); }

  void shiftlook(projection *P, const vector *V) {
vector W; if (!P || !V) return; W=*V;
if (P->str) W.z/=P->stretch; vectorincr(&(P->look),&W); }

  void shiftlookpoint(projection *P, real DX, real DY, real DZ) {
vector V; if (!P) return; setvector(&V,DX,DY,DZ); shiftlook(P,&V); }

  void shiftproj(projection *P, float DH, float DV, float DD) {
vector Del; if (!P) return; setvector(&Del,DH,-DV,DD);
scalevector(&Del,1.0 / P->mag); unrotatevector(&(P->rot),&Del); shiftlook(P,&Del); }

  void verticalstretch(projection *P, real Factor) {
if (!P) return; P->str=1; P->stretch=Factor; }

  static void checkmaglimits(projection *P) {
if (!P || !P->maglim) return; P->mag=fbetween(P->mag,P->mag1,P->mag2); }

  void setmaglimits(projection *P, real Lower, real Upper) {
if (!P) return; P->mag1=Lower; P->mag2=Upper; P->maglim=1; checkmaglimits(P); }

  void setmagnification(projection *P, real Mag) {
if (!P) return; P->mag=Mag; checkmaglimits(P); }

  void zoommagnification(projection *P, real Zoom) {
if (!P) return; if (!P->persp) { P->mag*=Zoom; checkmaglimits(P); return; }
P->pix*=Zoom; P->rads= P->width / P->pix;
if (P->rads>twopi) { P->rads=twopi; P->pix= P->width / twopi; }
P->wide= P->rads > 0.4; old(P->fov); P->fov=anglestring(P->rads,0); }


//  Projection query

  bool projatnadir(const projection *P) {
if (!P) return 0; return P->elev + pibytwo < 0.01; }


//  Saving projections to files

  ulong lengthofprojection(const projection *P) {
if (!P) return 0; return sizeof(projection) + lengthofstring(P->fov); }

  void writeprojection(const projection *P) {
if (!P) return; writebytes(P,sizeof(projection)); writestring(P->fov); }

  void readprojection(projection *P) {
if (!P) return; readbytes(P,sizeof(projection)); P->fov=readstring(); }


//  Rendering

  static long distfromrect(point P, rect R) {
long DH,DV;
DH=gtr(P.x - right(R),left(R) - P.x); if (DH<0) DH=0;
DV=gtr(P.y - top(R),bottom(R) - P.y); if (DV<0) DV=0;
return gtr(DH,DV); }

//  return values:
//    0: failure
//    1: success
//    2: okay, but too close to edge of 360¡ fish-eye view

  static int projvect(const projection *P, bool Dir, const vector *V, vector *R, real *Mag) {
real Rho2,Rho,Q; if (Dir) *R=*V; else vectordifference(V,&(P->look),R);
if (P->str) R->z*=P->stretch; rotatevector(&(P->rot),R);
if (!P->persp) { *Mag=P->mag; return 1; }
if (!P->wide) { if (R->z<=0.0) return 0; *Mag= P->pix / R->z; return 1; }
Rho2= (R->x * R->x) + (R->y * R->y); Rho=sqrt(Rho2); Q=atan2(Rho,R->z);
if (Q<0.01) *Mag= P->pix / R->z; else *Mag= P->pix * Q / Rho;
R->z=sqrt(Rho2 + (R->z * R->z)); if (Dir && (Q > pi * 0.9)) return 2; return 1; }

  static bool projpoint(const projection *P, bool Dir, const vector *V, image *D) {
vector R; real Mag; int RES; if (!P || !V || !D) return 0;
if (!(RES=projvect(P,Dir,V,&R,&Mag))) return 0;
R.y=-R.y;  //  Quartz
D->pt.x= P->org.x + (R.x * Mag); if (fabs(D->pt.x)>30000.0) return 0;  //  required?
D->pt.y= P->org.y + (R.y * Mag); if (fabs(D->pt.y)>30000.0) return 0;
D->mag=Mag; D->depth=R.z; D->miss=distfromrect(D->pt,P->port); return RES==1; }

  bool projectpoint(const projection *P, const vector *V, image *D) {
D->good=projpoint(P,0,V,D); return D->good; }

  bool projectdirection(const projection *P, const vector *V, image *D) {
if (!P->persp) { D->good=0; return 0; }
D->good=projpoint(P,1,V,D); return D->good; }

  static bool reverseplaneproj(const projection *P, real H, real V, vector *Ve) {
Ve->x= (H - P->org.x) / (real)P->mag; Ve->y= (V - P->org.y) / (real)P->mag; Ve->z=0.0;
Ve->y=-Ve->y;
unrotatevector(&(P->rot),Ve); vectorincr(Ve,&(P->look)); return 1; }

  static bool reverseperspective(const projection *P, real Horz, real Vert, vector *U) {
real X,Y,R,Beta; X= (Horz - P->org.x) / P->pix; Y= (Vert - P->org.y) / P->pix;
Y=-Y;
R=sqrt((X * X) + (Y * Y)); if (R>=pi) return 0;
Beta=atan2(Y,X); polartocartesian(1.0,R,Beta,U);
unrotatevector(&(P->rot),U); normalise(U); return 1; }

  bool reverseprojection(const projection *P, point Q, vector *V) {
if (!P || !V) return 0;
return P->persp ? reverseperspective(P,Q.x,Q.y,V) : reverseplaneproj(P,Q.x,Q.y,V); }

  bool reversevector(const projection *P, vector *Q, vector *V) {
if (!P || !Q || !V) return 0;
return P->persp ? reverseperspective(P,Q->x,Q->y,V) : reverseplaneproj(P,Q->x,Q->y,V); }


//  Rotation-solving utility used by solveperspdrag() and solvespheredrag()

  static void setrot(matrix *R, real Azim, real Elev) {
Elev=fbetween(Elev,-pibytwo,pibytwo); Azim=modtwopisym(Azim);
eulermatrix(Azim + pibytwo,pibytwo - Elev,pi,R); }

  static void solverotation(vector W, vector U, bool Globe, real *Az, real *El) {
real A,B,Ang1,Ang2,UZ1,UZ2,Rho2,Rho,Elev,Azim,SaveUZ; matrix Rot; if (!Az || !El) return;
//  W is the point on the sky or globe that is being dragged
//  U is the rotated (screen-sense) point to which the user has dragged it
//  W and U are both normalised, normally
if (Globe) {
  Rho2= (U.x * U.x) + (U.y * U.y);
  if (Rho2<1.0) U.z=-sqrt(1.0 - Rho2); else { U.z=-0.0; Rho=sqrt(Rho2); U.x/=Rho; U.y/=Rho; } }
SaveUZ=U.z;  //  used for discriminating Azim options later
//  Problem:  find angles that make a rotation matrix that maps W into U
//  U.x= (     SA * W.x) + (    -CA * W.y);                #1
//  U.y= (CA * SE * W.x) + (SA * SE * W.y) + (-CE * W.z);  #2
//  U.z= (CA * CE * W.x) + (SA * CE * W.y) + ( SE * W.z);  #3
//  #1 says that |U.x| must be at most sqrt((W.x * W.x) + (W.y * W.y))
Rho=sqrt((W.x * W.x) + (W.y * W.y)); U.x=fbetween(U.x,-Rho,Rho);
//  (CE * #2) - (SE * #3) gives:
//  (CE * U.y) - (SE * U.z) = -W.z;  #4
//  There will be two solutions in general.
//  I want the one that gives the smallest Elev angle (equator closer to centre of view).
Rho=sqrt((U.y * U.y) + (U.z * U.z)); A= -W.z / Rho;
if (A>=1.0) A=0.0; else if (A<=-1.0) A=pi; else A=acos(A);
B=atan2(U.z,U.y); Ang1= A - B; Ang2=-(A + B);
Elev= (Ang1 * Ang1 < Ang2 * Ang2) ? Ang1 : Ang2;
//  Now, using #1,
Rho=sqrt((W.x * W.x) + (W.y * W.y));
A= -U.x / Rho; A=acos(fbetween(A,-1.0,1.0));
B=atan2(W.x,W.y); Ang1= A - B; Ang2=-(A + B);
//  The solution I want has Azim angle that makes U.z the most like the argument.
//  Is that cheating?  There may be a more appropriate condition to check.
setrot(&Rot,Ang1,Elev); U=W; rotatevector(&Rot,&U); UZ1=U.z;
setrot(&Rot,Ang2,Elev); U=W; rotatevector(&Rot,&U); UZ2=U.z;
Azim= (fabs(SaveUZ - UZ1) < fabs(SaveUZ - UZ2)) ? Ang1 : Ang2;
//  In that case of the globe, it is the one that gives U.z the closest to the user.
if (Globe) Azim= (UZ1 < UZ2) ? Ang1 : Ang2;
*Az=Azim; *El=Elev; }


//  Dragging a perspective view

  void solveperspdrag(projection *Proj, const vector *Where, point P) {
vector U,W; real Beta,Rho,Elev,Azim; if (!Proj || !Where) return;
W=*Where; U.x= (P.x - Proj->org.x) / Proj->pix; U.y= (P.y - Proj->org.y) / Proj->pix;
U.y=-U.y;  //  Quartz
Rho= sqrt((U.x * U.x) + (U.y * U.y));
if (Rho>pi) { U.x*= pi / Rho; U.y*= pi / Rho; Rho=pi; }
Beta=atan2(U.y,U.x); polartocartesian(1.0,Rho,Beta,&U);
solverotation(W,U,0,&Azim,&Elev); setazimelev(Proj,Azim,Elev); }


//  These procedures are for projecting onto an opaque unit spherical surface centred
//  in the view.

  bool projontosphere(const projection *Proj, const vector *V, image *D) {
if (!Proj || !V || !D) return 0;
projectpoint(Proj,V,D); if (D->depth>0.0) D->good=0; return D->good; }

  bool reverseprojontosphere(const projection *Proj, point P, vector *U) {
real Rho2; if (!Proj || !U) return 0;
U->x= (P.x - Proj->org.x) / Proj->mag; U->y= (P.y - Proj->org.y) / Proj->mag;
U->y=-U->y;  //  Quartz
Rho2= (U->x * U->x) + (U->y * U->y); if (Rho2>=1.0) return 0;
U->z=-sqrt(1.0 - Rho2); /*negative depth choice*/ if (isnan(U->z)) return 0;
unrotatevector(&(Proj->rot),U); normalise(U); return 1; }

  void solvespheredrag(projection *Proj, const vector *Where, point P) {
vector U,W; real A,B,Ang1,Ang2,UZ1,UZ2,Rho,Elev,Azim,SaveUZ; if (!Proj || !Where) return;
W=*Where; U.x= (P.x - Proj->org.x) / Proj->mag; U.y= (P.y - Proj->org.y) / Proj->mag;
U.y=-U.y;  //  Quartz
Rho=sqrt((U.x * U.x) + (U.y * U.y));
if (Rho>1.0) { U.x/=Rho; U.y/=Rho; }

solverotation(W,U,1,&Azim,&Elev);
setazimelev(Proj,Azim,Elev);
return;

//  problem:  find angles that make a rotation matrix that maps W into U
//  U.x= (     SA * W.x) + (    -CA * W.y);                #1
//  U.y= (CA * SE * W.x) + (SA * SE * W.y) + (-CE * W.z);  #2
//  U.z= (CA * CE * W.x) + (SA * CE * W.y) + ( SE * W.z);  #3
//  #1 says that |U.x| must be at most sqrt((W.x * W.x) + (W.y * W.y))
Rho=sqrt((W.x * W.x) + (W.y * W.y)); U.x=fbetween(U.x,-Rho,Rho);
U.z= 1.0 - (U.x * U.x) - (U.y * U.y); if (U.z<=0.0) U.z=-0.0; else U.z=-sqrt(U.z);
SaveUZ=U.z;  //  used for discriminating Azim options later
//  (CE * #2) - (SE * #3) gives:
//  (CE * U.y) - (SE * U.z) = -W.z;  #4
//  There will be two solutions in general.
//  I want the one that gives the smallest Elev angle (equator closer to centre of view).
Rho=sqrt((U.y * U.y) + (U.z * U.z)); A= -W.z / Rho;
if (A>=1.0) A=0.0; else if (A<=-1.0) A=pi; else A=acos(A);
B=atan2(U.z,U.y); Ang1= A - B; Ang2=-(A + B);
Elev= (Ang1 * Ang1 < Ang2 * Ang2) ? Ang1 : Ang2;
//  Now, using #1,
Rho=sqrt((W.x * W.x) + (W.y * W.y));
A= -U.x / Rho; A=acos(fbetween(A,-1.0,1.0));
B=atan2(W.x,W.y); Ang1= A - B; Ang2=-(A + B);
//  The solution I want has Azim angle that makes U.z the most negative (?)
setazimelev(Proj,Ang1,Elev); U=W; rotatevector(&(Proj->rot),&U); UZ1=U.z;
setazimelev(Proj,Ang2,Elev); U=W; rotatevector(&(Proj->rot),&U); UZ2=U.z;
Azim= (UZ1 < UZ2) ? Ang1 : Ang2;
Azim= (fabs(SaveUZ - UZ1) < fabs(SaveUZ - UZ2)) ? Ang1 : Ang2; //test
setazimelev(Proj,Azim,Elev); }


//

  bool clicktovector(const mouseclick *M, vector *V) {
window *W; projection *P; if (!M) return 0;
W=M->wind; if (!W) return 0; P=W->proj; if (!P) return 0;
if (!projatnadir(P) || !pointinrect(M->point,W->page)) { beep(); return 0; }
return reverseprojection(P,M->point,V); }

  bool imagerotation(const projection *P, const image *D, matrix *R) {
real X,Y,Q,N; if (!P || !D || !R) return 0;
if (!P->persp) { *R=P->rot; return 1; }
X= (D->pt.x - P->org.x) / P->pix; Y= (D->pt.y - P->org.y) / P->pix;
Q=sqrt((X * X) + (Y * Y)); N=atan2(X,Y);
eulermatrix(-N,-Q,N,R); matrixproduct(R,&(P->rot),R); return 1; }

  void vectmove(window *W, const projection *P, const vector *V) {
image D; if (!W || !P || !V) return; projectpoint(P,V,&D); if (D.good) move(W,D.pt.x,D.pt.y); }

  void vectmoveto(window *W, const projection *P, const vector *V) {
image D; if (!W || !P || !V) return; projectpoint(P,V,&D); if (D.good) moveto(W,D.pt.x,D.pt.y); }

  void vectline(window *W, const projection *P, const vector *V) {
image D; if (!W || !P || !V) return; projectpoint(P,V,&D); if (D.good) line(W,D.pt.x,D.pt.y); }

  void vectlineto(window *W, const projection *P, const vector *V) {
image D; if (!W || !P || !V) return; projectpoint(P,V,&D); if (D.good) lineto(W,D.pt.x,D.pt.y); }

  void projmove(window *W, const projection *P, real X, real Y, real Z) {
vector V; V.x=X; V.y=Y; V.z=Z; vectmove(W,P,&V); }

  void projmoveto(window *W, const projection *P, real X, real Y, real Z) {
vector V; V.x=X; V.y=Y; V.z=Z; vectmoveto(W,P,&V); }

  void projline(window *W, const projection *P, real X, real Y, real Z) {
vector V; V.x=X; V.y=Y; V.z=Z; vectline(W,P,&V); }

  void projlineto(window *W, const projection *P, real X, real Y, real Z) {
vector V; V.x=X; V.y=Y; V.z=Z; vectlineto(W,P,&V); }

  void movetolineto(window *W, const image *P, bool *PenDown) {
if (!P || !PenDown) return; if (!P->good) { *PenDown=0; return; }
if (*PenDown) lineto(W,P->pt.x,P->pt.y); else moveto(W,P->pt.x,P->pt.y); *PenDown=1; }

/*  void addpentorect(CGRect *R, image *D, float Width) {
CGRect Pen; if (!R || !D) return;
setpenrect(&Pen,D->pt,Width); *R=CGRectUnion(*R,Pen); }*/

  void imagepoly(window *W, const image *Im, int N) {
int i; if (!W || !Im) return;
movept(W,Im[0].pt); for (i=1;i<N;i++) linept(W,Im[i].pt); linept(W,Im[0].pt); }