/* R5: fixed-address micromirror modulation -> projection lens -> folding mirror.
   Source pixels provide exposure-averaged ON duty. No inverse-CDF target aiming.
   Positions and normals are shared with the visual renderer. Geometry is cached
   until lens/aperture/tilt changes; each decoded frame is forward-splatted. */
#include <stdint.h>
#ifdef NATIVE
#include <math.h>
#define EX
#else
#define EX __attribute__((visibility("default")))
extern double sin(double),cos(double);
#endif
#define MAXN (720*720)
#define MAXR (1080*1080)
#define SAMPLES 5
static int fw,fh,rw,rh,N,color=0,isolate=-1,built=0,live=0,rebuilds=0;
static double ow,oh,ww,wh,focus=0,aperture=.50,tilt=0,mass=0,rms=0,spot=.18;
static unsigned char input[MAXN*4],output[MAXR*4],valid[MAXN*SAMPLES];
static float centers[MAXN*3],normals[MAXN*3],sources[MAXN*3],collimators[MAXN*3],lenses[MAXN*3],folds[MAXN*3],hits[MAXN*2],duty[MAXN];
static float mapx[MAXN*SAMPLES],mapy[MAXN*SAMPLES];
static float energy[3][MAXR],temp[MAXR];
static const double S=2.2,F0=2.0,YF=2.75,SZ=-21.45,FC=.66,DC=3.0;
static const double sn=.20791169081775931,cs=.9781476007338057;
static const double si=.4067366430758002,ci=.9135454576426009;
static const double slopes[5][2]={{0,0},{.024,0},{-.024,0},{0,.024},{0,-.024}};
static const double sub[5][2]={{0,0},{.27,.21},{-.23,-.18},{.21,-.26},{-.25,.24}};
static double sq(double a){return a*a;}
static double root(double x){return __builtin_sqrt(x);}
static double cl(double x,double a,double b){return x<a?a:x>b?b:x;}
static void unit(double *x,double *y,double *z){double l=root(sq(*x)+sq(*y)+sq(*z));*x/=l;*y/=l;*z/=l;}
static void reflect(double ix,double iy,double iz,double nx,double ny,double nz,double *rx,double *ry,double *rz){double d=ix*nx+iy*ny+iz*nz;*rx=ix-2*d*nx;*ry=iy-2*d*ny;*rz=iz-2*d*nz;}
/* Independent surfaces determine the landing point; image coordinates are only
   used to place the physical object pixel, never as a destination to draw. */
static int trace(int i,int k,double *hx,double *hy,int store){
 int col=i%fw,row=i/fw,j=i*3;
 double du=sub[k][0]*ow/fw,dv=sub[k][1]*oh/fh;
 double px=(.5-(col+.5)/fw)*ow+du*cs,py=du*sn,pz=(.5-(row+.5)/fh)*oh+dv;
 double nx=-sin(.20943951023931953+tilt),ny=cos(.20943951023931953+tilt),nz=0;
 if(store){centers[j]=px;centers[j+1]=py;centers[j+2]=pz;normals[j]=nx;normals[j+1]=ny;normals[j+2]=nz;
   px=centers[j];py=centers[j+1];pz=centers[j+2];nx=normals[j];ny=normals[j+1];nz=normals[j+2];unit(&nx,&ny,&nz);}
 double ox=slopes[k][0],oy=1,oz=slopes[k][1];unit(&ox,&oy,&oz);
 double ix,iy,iz;reflect(ox,oy,oz,-sn,cs,0,&ix,&iy,&iz);
 /* Collimator normal Ic=(si,-ci,0). Every sub-source gives a near-parallel
    illumination plane wave, not a ray from a miniature video screen. */
 double dot=ix*si-iy*ci,ell=(px*si-py*ci+DC)/dot;
 double qx=px-ix*ell,qy=py-iy*ell,qz=pz-iz*ell;
 double cx=-si*DC,cy=ci*DC;
 double tx=ix/dot-si,ty=iy/dot+ci,tz=iz/dot;
 double sx=cx-si*FC-FC*tx,sy=cy+ci*FC-FC*ty,sz=-FC*tz;
 if(store){sources[j]=sx;sources[j+1]=sy;sources[j+2]=sz;collimators[j]=qx;collimators[j+1]=qy;collimators[j+2]=qz;}
 double rx,ry,rz;reflect(ix,iy,iz,nx,ny,nz,&rx,&ry,&rz);
 if(ry<=.001)return 0;
 double d=(S-py)/ry,lx=px+rx*d,lz=pz+rz*d;
 if(store){lenses[j]=lx;lenses[j+1]=S;lenses[j+2]=lz;}
 if(sq(lx)+sq(lz)>sq(aperture))return 0;
 /* Paraxial thin lens: transverse slope -= height/f. */
 double bx=rx/ry-lx/(F0+focus),bz=rz/ry-lz/(F0+focus);
 /* Exact intersection with the 45 degree folding plane y+z=YF. */
 double fd=(YF-S-lz)/(1+bz),fx=lx+bx*fd,fy=S+fd,fz=lz+bz*fd;
 if(store){folds[j]=fx;folds[j+1]=fy;folds[j+2]=fz;}
 if(fd<=0||sq(fx)+sq(fy-YF)+sq(fz)>.50*.50)return 0;
 /* Reflection (bx,1,bz) -> (bx,-bz,-1), intersect screen z=SZ. */
 double sd=fz-SZ;*hx=fx+bx*sd;*hy=fy-bz*sd;
 if(store){hits[i*2]=*hx;hits[i*2+1]=*hy;}
 return sd>0;
}
static void rebuild(void){
 rebuilds++;
 double e=0;live=0;
 for(int i=0;i<N;i++)for(int k=0;k<SAMPLES;k++){
  double x=0,y=0;int ok=trace(i,k,&x,&y,k==0);int a=i*SAMPLES+k;valid[a]=ok;
  mapx[a]=(x/ww+.5)*rw-.5;mapy[a]=(.5-(y-YF)/wh)*rh-.5;
  if(k==0&&ok){double ex=((i%fw+.5)/fw-.5)*ww,ey=YF+(.5-(i/fw+.5)/fh)*wh;e+=sq(x-ex)+sq(y-ey);live++;}
 }
 rms=root(e/(live?live:1));built=1;
}
EX int configure(int w,int h,int outw,int outh,double objectWidth){
 if(w<2||h<2||w>720||h>720||outw<2||outh<2||outw>1080||outh>1080)return -1;
 fw=w;fh=h;rw=outw;rh=outh;N=w*h;ow=objectWidth;oh=ow*h/w;ww=ow*10;wh=oh*10;
 focus=0;aperture=.50;tilt=0;isolate=-1;spot=.18;built=0;rebuild();return 0;
}
EX void set_optics(double f,double a,double t){f=cl(f,-.35,.35);a=cl(a,.2,1.4);t=cl(t,-.12,.12);if(f!=focus||a!=aperture||t!=tilt){focus=f;aperture=a;tilt=t;built=0;rebuild();}}
EX void set_isolate(int i){isolate=(i>=0&&i<N)?i:-1;}
EX void set_spot(double s){spot=cl(s,0,.28);}
EX void set_frame(int followColor){
 double sum=0;int maxColor=0;
 for(int i=0;i<N;i++){int p=i*4,r=input[p],g=input[p+1],b=input[p+2];double l=(.2126*r+.7152*g+.0722*b)/255.;duty[i]=l;sum+=l;
   int d=r>g?r-g:g-r,ee=b>g?b-g:g-b;if(d>maxColor)maxColor=d;if(ee>maxColor)maxColor=ee;}
 mass=sum/N;color=followColor&&maxColor>5;
}
static void blur(float *a){
 if(spot<=0)return;double mid=1-2*spot;
 for(int y=0;y<rh;y++)for(int x=0;x<rw;x++){int p=y*rw+x;temp[p]=a[p]*mid+spot*(a[y*rw+(x?x-1:x)]+a[y*rw+(x+1<rw?x+1:x)]);}
 for(int y=0;y<rh;y++)for(int x=0;x<rw;x++){int p=y*rw+x;a[p]=temp[p]*mid+spot*(temp[(y?y-1:y)*rw+x]+temp[(y+1<rh?y+1:y)*rw+x]);}
}
EX void solve(double power){
 if(!built)rebuild();int chans=color?3:1;int size=rw*rh;
 for(int c=0;c<chans;c++)for(int p=0;p<size;p++)energy[c][p]=0;
 double scale=(double)size/N/SAMPLES*cl(power,0,1);
 for(int i=0;i<N;i++){
   if((isolate>=0&&i!=isolate)||duty[i]<=0)continue;
   double rgb[3];if(color){rgb[0]=input[i*4]*scale;rgb[1]=input[i*4+1]*scale;rgb[2]=input[i*4+2]*scale;}else rgb[0]=duty[i]*255*scale;
   for(int k=0;k<SAMPLES;k++){int a=i*SAMPLES+k;if(!valid[a])continue;double x=mapx[a],y=mapy[a];int x0=(int)__builtin_floor(x),y0=(int)__builtin_floor(y);double ax=x-x0,ay=y-y0;
     /* Low-resolution detector bins integrate each pixel footprint. No
        bilinear tent leaks light into a neighbouring OFF address. */
     if(rw==fw && rh==fh){int X=(int)__builtin_floor(x+.5),Y=(int)__builtin_floor(y+.5);if(X>=0&&X<rw&&Y>=0&&Y<rh)for(int c=0;c<chans;c++)energy[c][Y*rw+X]+=rgb[c];continue;}
     for(int yy=0;yy<2;yy++){int row=y0+yy;if(row<0||row>=rh)continue;double wy=yy?ay:1-ay;for(int xx=0;xx<2;xx++){int col=x0+xx;if(col<0||col>=rw)continue;double w=wy*(xx?ax:1-ax);int p=row*rw+col;for(int c=0;c<chans;c++)energy[c][p]+=w*rgb[c];}}
   }
 }
 for(int c=0;c<chans;c++)blur(energy[c]);
 for(int p=0;p<size;p++){output[p*4]=(unsigned char)cl(energy[0][p]+.5,0,255);output[p*4+1]=color?(unsigned char)cl(energy[1][p]+.5,0,255):output[p*4];output[p*4+2]=color?(unsigned char)cl(energy[2][p]+.5,0,255):output[p*4];output[p*4+3]=255;}
}
EX double get_rms(void){return rms;}EX double get_mass(void){return mass;}EX int get_color(void){return color;}EX int get_live(void){return live;}
#define PTR(name) EX uintptr_t name##_ptr(void){return (uintptr_t)name;}
PTR(input) PTR(output) PTR(centers) PTR(normals) PTR(sources) PTR(collimators) PTR(lenses) PTR(folds) PTR(hits) PTR(duty) PTR(mapx) PTR(mapy) PTR(valid)

EX int get_rebuilds(void){return rebuilds;}
/* Presentation-only continuous ray raster. Optical rays have already been
   computed above; this does NOT create or modify the receiver buffer. */
#ifndef NATIVE
extern float expf(float);
#endif
static float beam_segments[20000*5],beam_acc[MAXR];
static unsigned char beam_pixels[MAXR*4];
EX uintptr_t beam_segments_ptr(void){return (uintptr_t)beam_segments;}
EX uintptr_t beam_pixels_ptr(void){return (uintptr_t)beam_pixels;}
EX void draw_beams(int count,int w,int h,double exposure){
 if(count<0||count>20000||w<2||h<2||w*h>MAXR)return;
 for(int i=0;i<w*h;i++)beam_acc[i]=0;
 for(int k=0;k<count;k++){float *s=beam_segments+k*5,dx=s[2]-s[0],dy=s[3]-s[1];float ax=dx<0?-dx:dx,ay=dy<0?-dy:dy;int n=(int)__builtin_ceil(ax>ay?ax:ay);if(n<1)n=1;if(n>10000)continue;
  for(int j=0;j<=n;j++){float x=s[0]+dx*j/n,y=s[1]+dy*j/n;int ix=(int)__builtin_floor(x),iy=(int)__builtin_floor(y);float a=x-ix,b=y-iy;
   for(int yy=0;yy<2;yy++)for(int xx=0;xx<2;xx++){int X=ix+xx,Y=iy+yy;if(X>=0&&X<w&&Y>=0&&Y<h)beam_acc[Y*w+X]+=s[4]*(xx?a:1-a)*(yy?b:1-b);}
  }
 }
 for(int i=0;i<w*h;i++){unsigned char g=(unsigned char)(255*(1-expf(-beam_acc[i]*exposure))+.5f);beam_pixels[i*4]=g;beam_pixels[i*4+1]=g;beam_pixels[i*4+2]=g;beam_pixels[i*4+3]=255;}
}
