--- ray/src/common/image.c 2008/12/12 22:05:38 2.35 +++ ray/src/common/image.c 2022/01/13 00:26:09 2.54 @@ -1,5 +1,5 @@ #ifndef lint -static const char RCSid[] = "$Id: image.c,v 2.35 2008/12/12 22:05:38 greg Exp $"; +static const char RCSid[] = "$Id: image.c,v 2.54 2022/01/13 00:26:09 greg Exp $"; #endif /* * image.c - routines for image generation. @@ -15,11 +15,6 @@ static const char RCSid[] = "$Id: image.c,v 2.35 2008/ #include "paths.h" #include "view.h" - -#define FEQ(x,y) (fabs((x)-(y)) <= FTINY) -#define VEQ(v,w) (FEQ((v)[0],(w)[0]) && FEQ((v)[1],(w)[1]) \ - && FEQ((v)[2],(w)[2])) - VIEW stdview = STDVIEW; /* default view parameters */ static gethfunc gethview; @@ -33,7 +28,8 @@ VIEW *v static char ill_horiz[] = "illegal horizontal view size"; static char ill_vert[] = "illegal vertical view size"; - if (v->vaft < -FTINY || (v->vaft > FTINY && v->vaft <= v->vfore)) + if ((v->vfore < -FTINY) | (v->vaft < -FTINY) || + (v->vaft > FTINY) & (v->vaft <= v->vfore)) return("illegal fore/aft clipping plane"); if (v->vdist <= FTINY) @@ -67,8 +63,8 @@ VIEW *v return(ill_horiz); if (v->vert >= 180.0-FTINY) return(ill_vert); - v->hn2 = 2.0 * tan(v->horiz*(PI/180.0/2.0)); - v->vn2 = 2.0 * tan(v->vert*(PI/180.0/2.0)); + v->hn2 = 2.0 * tan(v->horiz*(PI/360.)); + v->vn2 = 2.0 * tan(v->vert*(PI/360.)); break; case VT_CYL: /* cylindrical panorama */ if (v->horiz > 360.0+FTINY) @@ -76,7 +72,7 @@ VIEW *v if (v->vert >= 180.0-FTINY) return(ill_vert); v->hn2 = v->horiz * (PI/180.0); - v->vn2 = 2.0 * tan(v->vert*(PI/180.0/2.0)); + v->vn2 = 2.0 * tan(v->vert*(PI/360.)); break; case VT_ANG: /* angular fisheye */ if (v->horiz > 360.0+FTINY) @@ -91,18 +87,18 @@ VIEW *v return(ill_horiz); if (v->vert > 180.0+FTINY) return(ill_vert); - v->hn2 = 2.0 * sin(v->horiz*(PI/180.0/2.0)); - v->vn2 = 2.0 * sin(v->vert*(PI/180.0/2.0)); + v->hn2 = 2.0 * sin(v->horiz*(PI/360.)); + v->vn2 = 2.0 * sin(v->vert*(PI/360.)); break; case VT_PLS: /* planispheric fisheye */ if (v->horiz >= 360.0-FTINY) return(ill_horiz); if (v->vert >= 360.0-FTINY) return(ill_vert); - v->hn2 = 2.*sin(v->horiz*(PI/180.0/2.0)) / - (1.0 + cos(v->horiz*(PI/180.0/2.0))); - v->vn2 = 2.*sin(v->vert*(PI/180.0/2.0)) / - (1.0 + cos(v->vert*(PI/180.0/2.0))); + v->hn2 = 2.*sin(v->horiz*(PI/360.)) / + (1.0 + cos(v->horiz*(PI/360.))); + v->vn2 = 2.*sin(v->vert*(PI/360.)) / + (1.0 + cos(v->vert*(PI/360.))); break; default: return("unknown view type"); @@ -124,6 +120,79 @@ VIEW *v } +char * +cropview( /* crop a view to the indicated bounds */ +VIEW *v, +double x0, +double y0, +double x1, +double y1 +) +{ + static char ill_hemi[] = "illegal crop for hemispherical view"; + double d; + /* order crop extrema */ + if (x0 > x1) { d=x0; x0=x1; x1=d; } + if (y0 > y1) { d=y0; y0=y1; y1=d; } + + if ((x1-x0 <= FTINY) | (y1-y0 <= FTINY)) + return("zero crop area"); + + d = x1 - x0; /* adjust horizontal size? */ + if (!FABSEQ(d, 1.)) + switch (v->type) { + case VT_PER: + v->horiz = 360./PI*atan( d*tan(PI/360.*v->horiz) ); + break; + case VT_PAR: + case VT_ANG: + case VT_CYL: + v->horiz *= d; + break; + case VT_HEM: + d *= sin(PI/360.*v->horiz); + if (d > 1.) + return(ill_hemi); + v->horiz = 360./PI*asin( d ); + break; + case VT_PLS: + d *= sin(PI/360.*v->horiz) / + (1. + cos(PI/360.*v->horiz)); + v->horiz = 360./PI*acos( (1. - d*d) / (1. + d*d) ); + break; + } + + d = y1 - y0; /* adjust vertical size? */ + if (!FABSEQ(d, 1.)) + switch (v->type) { + case VT_PER: + case VT_CYL: + v->vert = 360./PI*atan( d*tan(PI/360.*v->vert) ); + break; + case VT_PAR: + case VT_ANG: + v->vert *= d; + break; + case VT_HEM: + d *= sin(PI/360.*v->vert); + if (d > 1.) + return(ill_hemi); + v->vert = 360./PI*asin( d ); + break; + case VT_PLS: + d *= sin(PI/360.*v->vert) / + (1. + cos(PI/360.*v->vert)); + v->vert = 360./PI*acos( (1. - d*d) / (1. + d*d) ); + break; + } + /* adjust offsets */ + v->hoff = ((x0 + x1)*.5 - .5 + v->hoff) / (x1 - x0); + v->voff = ((y0 + y1)*.5 - .5 + v->voff) / (y1 - y0); + + return(setview(v)); /* final error checks & set-up */ +} + + void normaspect( /* fix pixel aspect or resolution */ double va, /* view aspect ratio */ @@ -169,9 +238,7 @@ double y direc[0] = v->vdir[0] + x*v->hvec[0] + y*v->vvec[0]; direc[1] = v->vdir[1] + x*v->hvec[1] + y*v->vvec[1]; direc[2] = v->vdir[2] + x*v->hvec[2] + y*v->vvec[2]; - orig[0] = v->vp[0] + v->vfore*direc[0]; - orig[1] = v->vp[1] + v->vfore*direc[1]; - orig[2] = v->vp[2] + v->vfore*direc[2]; + VSUM(orig, v->vp, direc, v->vfore); d = normalize(direc); return(v->vaft > FTINY ? (v->vaft - v->vfore)*d : 0.0); case VT_HEM: /* hemispherical fisheye */ @@ -182,9 +249,7 @@ double y direc[0] = z*v->vdir[0] + x*v->hvec[0] + y*v->vvec[0]; direc[1] = z*v->vdir[1] + x*v->hvec[1] + y*v->vvec[1]; direc[2] = z*v->vdir[2] + x*v->hvec[2] + y*v->vvec[2]; - orig[0] = v->vp[0] + v->vfore*direc[0]; - orig[1] = v->vp[1] + v->vfore*direc[1]; - orig[2] = v->vp[2] + v->vfore*direc[2]; + VSUM(orig, v->vp, direc, v->vfore); return(v->vaft > FTINY ? v->vaft - v->vfore : 0.0); case VT_CYL: /* cylindrical panorama */ d = x * v->horiz * (PI/180.0); @@ -193,9 +258,7 @@ double y direc[0] = z*v->vdir[0] + x*v->hvec[0] + y*v->vvec[0]; direc[1] = z*v->vdir[1] + x*v->hvec[1] + y*v->vvec[1]; direc[2] = z*v->vdir[2] + x*v->hvec[2] + y*v->vvec[2]; - orig[0] = v->vp[0] + v->vfore*direc[0]; - orig[1] = v->vp[1] + v->vfore*direc[1]; - orig[2] = v->vp[2] + v->vfore*direc[2]; + VSUM(orig, v->vp, direc, v->vfore); d = normalize(direc); return(v->vaft > FTINY ? (v->vaft - v->vfore)*d : 0.0); case VT_ANG: /* angular fisheye */ @@ -212,37 +275,33 @@ double y direc[0] = z*v->vdir[0] + x*v->hvec[0] + y*v->vvec[0]; direc[1] = z*v->vdir[1] + x*v->hvec[1] + y*v->vvec[1]; direc[2] = z*v->vdir[2] + x*v->hvec[2] + y*v->vvec[2]; - orig[0] = v->vp[0] + v->vfore*direc[0]; - orig[1] = v->vp[1] + v->vfore*direc[1]; - orig[2] = v->vp[2] + v->vfore*direc[2]; + VSUM(orig, v->vp, direc, v->vfore); return(v->vaft > FTINY ? v->vaft - v->vfore : 0.0); case VT_PLS: /* planispheric fisheye */ x *= sqrt(v->hn2); y *= sqrt(v->vn2); d = x*x + y*y; z = (1. - d)/(1. + d); - d = d <= FTINY*FTINY ? PI : sqrt((1.0 - z*z)/d); - x *= d; - y *= d; + x *= (1. + z); + y *= (1. + z); direc[0] = z*v->vdir[0] + x*v->hvec[0] + y*v->vvec[0]; direc[1] = z*v->vdir[1] + x*v->hvec[1] + y*v->vvec[1]; direc[2] = z*v->vdir[2] + x*v->hvec[2] + y*v->vvec[2]; - orig[0] = v->vp[0] + v->vfore*direc[0]; - orig[1] = v->vp[1] + v->vfore*direc[1]; - orig[2] = v->vp[2] + v->vfore*direc[2]; + VSUM(orig, v->vp, direc, v->vfore); return(v->vaft > FTINY ? v->vaft - v->vfore : 0.0); } return(-1.0); } -void +int viewloc( /* find image location for point */ FVECT ip, VIEW *v, FVECT p -) +) /* Use VL_* flags to interpret return value */ { + int rflags = VL_GOOD; double d, d2; FVECT disp; @@ -254,17 +313,19 @@ FVECT p break; case VT_PER: /* perspective view */ d = DOT(disp,v->vdir); + if ((v->vaft > FTINY) & (d >= v->vaft)) + rflags |= VL_BEYOND; ip[2] = VLEN(disp); - if (d < 0.0) { /* fold pyramid */ + if (d < -FTINY) { /* fold pyramid */ ip[2] = -ip[2]; d = -d; - } - if (d > FTINY) { - d = 1.0/d; - disp[0] *= d; - disp[1] *= d; - disp[2] *= d; - } + } else if (d <= FTINY) + return(VL_BAD); /* at infinite edge */ + d = 1.0/d; + disp[0] *= d; + disp[1] *= d; + disp[2] *= d; + if (ip[2] < 0.0) d = -d; ip[2] *= (1.0 - v->vfore*d); break; case VT_HEM: /* hemispherical fisheye */ @@ -279,42 +340,57 @@ FVECT p d = DOT(disp,v->hvec); d2 = DOT(disp,v->vdir); ip[0] = 180.0/PI * atan2(d,d2) / v->horiz + 0.5 - v->hoff; - d = 1.0/sqrt(d*d + d2*d2); + d2 = d*d + d2*d2; + if (d2 <= FTINY*FTINY) + return(VL_BAD); /* at pole */ + if ((v->vaft > FTINY) & (d2 >= v->vaft*v->vaft)) + rflags |= VL_BEYOND; + d = 1.0/sqrt(d2); ip[1] = DOT(disp,v->vvec)*d/v->vn2 + 0.5 - v->voff; ip[2] = VLEN(disp); ip[2] *= (1.0 - v->vfore*d); - return; + goto gotall; case VT_ANG: /* angular fisheye */ ip[0] = 0.5 - v->hoff; ip[1] = 0.5 - v->voff; ip[2] = normalize(disp) - v->vfore; d = DOT(disp,v->vdir); if (d >= 1.0-FTINY) - return; + goto gotall; if (d <= -(1.0-FTINY)) { ip[0] += 180.0/v->horiz; - return; + goto gotall; } d = (180.0/PI)*acos(d) / sqrt(1.0 - d*d); ip[0] += DOT(disp,v->hvec)*d/v->horiz; ip[1] += DOT(disp,v->vvec)*d/v->vert; - return; + goto gotall; case VT_PLS: /* planispheric fisheye */ ip[0] = 0.5 - v->hoff; ip[1] = 0.5 - v->voff; ip[2] = normalize(disp) - v->vfore; d = DOT(disp,v->vdir); if (d >= 1.0-FTINY) - return; + goto gotall; if (d <= -(1.0-FTINY)) - return; /* really an error */ - d = sqrt(1.0 - d*d) / (1.0 + d); - ip[0] += DOT(disp,v->hvec)*d/sqrt(v->hn2); - ip[1] += DOT(disp,v->vvec)*d/sqrt(v->vn2); - return; + return(VL_BAD); + ip[0] += DOT(disp,v->hvec)/((1. + d)*sqrt(v->hn2)); + ip[1] += DOT(disp,v->vvec)/((1. + d)*sqrt(v->vn2)); + goto gotall; + default: + return(VL_BAD); } ip[0] = DOT(disp,v->hvec)/v->hn2 + 0.5 - v->hoff; ip[1] = DOT(disp,v->vvec)/v->vn2 + 0.5 - v->voff; +gotall: /* add appropriate return flags */ + if (ip[2] <= 0.0) + rflags |= VL_BEHIND; + else if ((v->type != VT_PER) & (v->type != VT_CYL)) + rflags |= VL_BEYOND*((v->vaft > FTINY) & + (ip[2] >= v->vaft - v->vfore)); + rflags |= VL_OUTSIDE*((0.0 >= ip[0]) | (ip[0] >= 1.0) | + (0.0 >= ip[1]) | (ip[1] >= 1.0)); + return(rflags); } @@ -326,7 +402,7 @@ int px, int py ) { - register int x, y; + int x, y; if (rp->rt & YMAJOR) { x = px; @@ -352,10 +428,11 @@ double lx, double ly ) { - register int x, y; + int x, y; - x = lx * rp->xr; - y = ly * rp->yr; + x = (int)(lx*rp->xr + .5 - (lx < 0.0)); + y = (int)(ly*rp->yr + .5 - (ly < 0.0)); + if (rp->rt & XDECR) x = rp->xr-1 - x; if (rp->rt & YDECR) @@ -377,14 +454,14 @@ int ac, char *av[] ) { -#define check(c,l) if ((av[0][c]&&av[0][c]!=' ') || \ +#define check(c,l) if ((av[0][c]&&!isspace(av[0][c])) || \ badarg(ac-1,av+1,l)) return(-1) if (ac <= 0 || av[0][0] != '-' || av[0][1] != 'v') return(-1); switch (av[0][2]) { case 't': /* type */ - if (!av[0][3] || av[0][3]==' ') + if (!av[0][3] || isspace(av[0][3])) return(-1); check(4,""); v->type = av[0][3]; @@ -498,51 +575,51 @@ VIEW *vp ) { static char vwstr[128]; - register char *cp = vwstr; + char *cp = vwstr; *cp = '\0'; if (vp->type != stdview.type) { sprintf(cp, " -vt%c", vp->type); cp += strlen(cp); } - if (!VEQ(vp->vp,stdview.vp)) { + if (!VABSEQ(vp->vp,stdview.vp)) { sprintf(cp, " -vp %.6g %.6g %.6g", vp->vp[0], vp->vp[1], vp->vp[2]); cp += strlen(cp); } - if (!FEQ(vp->vdist,stdview.vdist) || !VEQ(vp->vdir,stdview.vdir)) { + if (!FABSEQ(vp->vdist,stdview.vdist) || !VABSEQ(vp->vdir,stdview.vdir)) { sprintf(cp, " -vd %.6g %.6g %.6g", vp->vdir[0]*vp->vdist, vp->vdir[1]*vp->vdist, vp->vdir[2]*vp->vdist); cp += strlen(cp); } - if (!VEQ(vp->vup,stdview.vup)) { + if (!VABSEQ(vp->vup,stdview.vup)) { sprintf(cp, " -vu %.6g %.6g %.6g", vp->vup[0], vp->vup[1], vp->vup[2]); cp += strlen(cp); } - if (!FEQ(vp->horiz,stdview.horiz)) { + if (!FABSEQ(vp->horiz,stdview.horiz)) { sprintf(cp, " -vh %.6g", vp->horiz); cp += strlen(cp); } - if (!FEQ(vp->vert,stdview.vert)) { + if (!FABSEQ(vp->vert,stdview.vert)) { sprintf(cp, " -vv %.6g", vp->vert); cp += strlen(cp); } - if (!FEQ(vp->vfore,stdview.vfore)) { + if (!FABSEQ(vp->vfore,stdview.vfore)) { sprintf(cp, " -vo %.6g", vp->vfore); cp += strlen(cp); } - if (!FEQ(vp->vaft,stdview.vaft)) { + if (!FABSEQ(vp->vaft,stdview.vaft)) { sprintf(cp, " -va %.6g", vp->vaft); cp += strlen(cp); } - if (!FEQ(vp->hoff,stdview.hoff)) { + if (!FABSEQ(vp->hoff,stdview.hoff)) { sprintf(cp, " -vs %.6g", vp->hoff); cp += strlen(cp); } - if (!FEQ(vp->voff,stdview.voff)) { + if (!FABSEQ(vp->voff,stdview.voff)) { sprintf(cp, " -vl %.6g", vp->voff); cp += strlen(cp); } @@ -557,8 +634,8 @@ char *s { static char *altname[]={NULL,VIEWSTR,"rpict","rview","rvu","rpiece","pinterp",NULL}; extern char *progname; - register char *cp; - register char **an; + char *cp; + char **an; /* add program name to list */ if (altname[0] == NULL) { for (cp = progname; *cp; cp++) @@ -569,7 +646,7 @@ char *s } /* skip leading path */ cp = s; - while (*cp && *cp != ' ') + while (*cp && !isspace(*cp)) cp++; while (cp > s && !ISDIRSEP(cp[-1])) cp--; @@ -621,7 +698,8 @@ RESOLU *rp if (rp != NULL && !fgetsresolu(rp, fp)) mvs.ok = 0; - fclose(fp); + if (fp != stdin) + fclose(fp); return(mvs.ok); }