home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1994 January / usenetsourcesnewsgroupsinfomagicjanuary1994.iso / sources / misc / volume21 / rayshade / part16 < prev    next >
Encoding:
Text File  |  1991-07-21  |  37.7 KB  |  1,500 lines
  1. Newsgroups: comp.sources.misc
  2. From: Rayshade Construction Co. <rayshade@weedeater.math.YALE.EDU>
  3. Subject:  v21i019:  rayshade - A raytracing package for UNIX, Part16/19
  4. Message-ID: <1991Jul21.223347.3819@sparky.IMD.Sterling.COM>
  5. X-Md4-Signature: efda1f5851c6f429b69764464f3e1055
  6. Date: Sun, 21 Jul 1991 22:33:47 GMT
  7. Approved: kent@sparky.imd.sterling.com
  8.  
  9. Submitted-by: Rayshade Construction Co. <rayshade@weedeater.math.YALE.EDU>
  10. Posting-number: Volume 21, Issue 19
  11. Archive-name: rayshade/part16
  12. Environment: UNIX, !16BIT
  13.  
  14. #! /bin/sh
  15. # This is a shell archive.  Remove anything before this line, then unpack
  16. # it by saving it into a file and typing "sh file".  To overwrite existing
  17. # files, type "sh file -c".  You can also feed this as standard input via
  18. # unshar, or by typing "sh <file", e.g..  If this archive is complete, you
  19. # will see the following message at the end:
  20. #        "End of archive 16 (of 19)."
  21. # Contents:  libray/libobj/hf.c raypaint/render.c
  22. # Wrapped by kolb@woody on Wed Jul 17 17:56:55 1991
  23. PATH=/bin:/usr/bin:/usr/ucb ; export PATH
  24. if test -f 'libray/libobj/hf.c' -a "${1}" != "-c" ; then 
  25.   echo shar: Will not clobber existing file \"'libray/libobj/hf.c'\"
  26. else
  27. echo shar: Extracting \"'libray/libobj/hf.c'\" \(17333 characters\)
  28. sed "s/^X//" >'libray/libobj/hf.c' <<'END_OF_FILE'
  29. X/*
  30. X * hf.c
  31. X *
  32. X * Copyright (C) 1989, 1991, Craig E. Kolb
  33. X * All rights reserved.
  34. X *
  35. X * This software may be freely copied, modified, and redistributed
  36. X * provided that this copyright notice is preserved on all copies.
  37. X *
  38. X * You may not distribute this software, in whole or in part, as part of
  39. X * any commercial product without the express consent of the authors.
  40. X *
  41. X * There is no warranty or other guarantee of fitness of this software
  42. X * for any purpose.  It is provided solely "as is".
  43. X *
  44. X * $Id: hf.c,v 4.0 91/07/17 14:38:15 kolb Exp Locker: kolb $
  45. X *
  46. X * $Log:    hf.c,v $
  47. X * Revision 4.0  91/07/17  14:38:15  kolb
  48. X * Initial version.
  49. X * 
  50. X */
  51. X#include "geom.h"
  52. X#include "hf.h"
  53. X
  54. Xstatic Methods *iHfMethods = NULL;
  55. Xstatic char hfName[] = "heighfield";
  56. X
  57. Xstatic void integrate_grid(), QueueTri();
  58. Xstatic int DDA2D(), CheckCell();
  59. Xstatic Float intHftri();
  60. Xstatic float minalt(), maxalt();
  61. X
  62. Xtypedef struct {
  63. X    int stepX, stepY;
  64. X    Float tDX, tDY;
  65. X    float minz, maxz;
  66. X    int outX, outY;
  67. X    Vector cp, pDX, pDY;
  68. X} Trav2D;
  69. X
  70. XhfTri *CreateHfTriangle(), *GetQueuedTri();
  71. X
  72. Xunsigned long HFTests, HFHits;
  73. X
  74. XHf *
  75. XHfCreate(filename)
  76. Xchar *filename;
  77. X{
  78. X    Hf *hf;
  79. X    FILE *fp;
  80. X    float val, *maxptr, *minptr;
  81. X    int i, j;
  82. X
  83. X    fp = fopen(filename, "r");
  84. X    if (fp == (FILE *)NULL) {
  85. X        RLerror(RL_ABORT, "Cannot open heightfield file \"%s\".",
  86. X                filename);
  87. X        return (Hf *)NULL;
  88. X    }
  89. X
  90. X    hf = (Hf *)Malloc(sizeof(Hf));
  91. X    /*
  92. X     * Make the following an option someday.
  93. X     */
  94. X    hf->BestSize = BESTSIZE;
  95. X    /*
  96. X     * Store the inverse for faster computation.
  97. X     */
  98. X    hf->iBestSize = 1. / (float)hf->BestSize;
  99. X    /*
  100. X     * Get HF size.
  101. X     */
  102. X    if (fread((char *)&hf->size, sizeof(int), 1, fp) == 0) {
  103. X        RLerror(RL_ABORT, "Cannot read height field size.");
  104. X        return (Hf *)NULL;
  105. X    }
  106. X        
  107. X    hf->data = (float **)share_malloc(hf->size * sizeof(float *));
  108. X    for (i = 0; i < hf->size; i++) {
  109. X        hf->data[i] = (float *)share_malloc(hf->size * sizeof(float));
  110. X        /*
  111. X         * Read in row of HF data.
  112. X         */
  113. X        if (fread((char *)hf->data[i],sizeof(float),hf->size,fp)
  114. X            != hf->size) {
  115. X            RLerror(RL_ABORT, "Not enough heightfield data.");
  116. X            return (Hf *)NULL;
  117. X        }
  118. X        for (j = 0; j < hf->size; j++) {
  119. X            val = hf->data[i][j];
  120. X            if (val <= HF_UNSET) {
  121. X                hf->data[i][j] = HF_UNSET;
  122. X                /*
  123. X                 * Don't include the point in min/max
  124. X                 * calculations.
  125. X                 */
  126. X                continue;
  127. X            }
  128. X            if (val > hf->maxz)
  129. X                hf->maxz = val;
  130. X            if (val < hf->minz)
  131. X                hf->minz = val;
  132. X        }
  133. X    }
  134. X    (void)fclose(fp);
  135. X    /*
  136. X     * Allocate levels of grid.  hf->levels = log base BestSize of hf->size
  137. X     */
  138. X    for (i = hf->size, hf->levels = 0; i > hf->BestSize; i /= hf->BestSize,
  139. X                hf->levels++)
  140. X            ;
  141. X    hf->levels++;
  142. X    hf->qsize = CACHESIZE;
  143. X    hf->q = (hfTri **)Calloc((unsigned)hf->qsize, sizeof(hfTri *));
  144. X    hf->qtail = 0;
  145. X
  146. X    hf->lsize = (int *)share_malloc(hf->levels * sizeof(int));
  147. X    hf->spacing = (float *)share_malloc(hf->levels * sizeof(float));
  148. X    hf->boundsmax = (float ***)share_malloc(hf->levels * sizeof(float **));
  149. X    hf->boundsmin = (float ***)share_malloc(hf->levels * sizeof(float **));
  150. X
  151. X    hf->spacing[0] = hf->size -1;
  152. X    hf->lsize[0] = (int)hf->spacing[0];
  153. X    hf->boundsmax[0] = (float **)share_malloc(hf->lsize[0]*sizeof(float *));
  154. X    hf->boundsmin[0] = (float **)share_malloc(hf->lsize[0]*sizeof(float *));
  155. X    /*
  156. X     * Compute initial bounding boxes
  157. X     */
  158. X    for (i = 0; i < hf->lsize[0]; i++) {
  159. X        hf->boundsmax[0][i]=(float *)share_malloc(hf->lsize[0]*sizeof(float));
  160. X        hf->boundsmin[0][i]=(float *)share_malloc(hf->lsize[0]*sizeof(float));
  161. X        maxptr = hf->boundsmax[0][i];
  162. X        minptr = hf->boundsmin[0][i];
  163. X        for (j = 0; j < hf->lsize[0]; j++) {
  164. X            *maxptr++ = maxalt(i, j, hf->data) + EPSILON;
  165. X            *minptr++ = minalt(i, j, hf->data) - EPSILON;
  166. X        }
  167. X    }
  168. X
  169. X    for (i = 1; i < hf->levels; i++) {
  170. X        hf->spacing[i] = hf->spacing[i-1] * hf->iBestSize;
  171. X        hf->lsize[i] = (int)hf->spacing[i];
  172. X        if ((Float)hf->lsize[i] != hf->spacing[i])
  173. X            hf->lsize[i]++;
  174. X        hf->boundsmax[i]=(float **)share_malloc(hf->lsize[i]*sizeof(float *));
  175. X        hf->boundsmin[i]=(float **)share_malloc(hf->lsize[i]*sizeof(float *));
  176. X        for (j = 0; j < hf->lsize[i]; j++) {
  177. X            hf->boundsmax[i][j] = (float *)share_malloc(hf->lsize[i] *
  178. X                            sizeof(float));
  179. X            hf->boundsmin[i][j] = (float *)share_malloc(hf->lsize[i] *
  180. X                            sizeof(float));
  181. X        }
  182. X        integrate_grid(hf, i);
  183. X    }
  184. X
  185. X    hf->boundbox[LOW][X] = hf->boundbox[LOW][Y] = 0;
  186. X    hf->boundbox[HIGH][X] = hf->boundbox[HIGH][Y] = 1;
  187. X    hf->boundbox[LOW][Z] = hf->minz;
  188. X    hf->boundbox[HIGH][Z] = hf->maxz;
  189. X
  190. X    return hf;
  191. X}
  192. X
  193. XMethods *
  194. XHfMethods()
  195. X{
  196. X    if (iHfMethods == (Methods *)NULL) {
  197. X        iHfMethods = MethodsCreate();
  198. X        iHfMethods->create = (GeomCreateFunc *)HfCreate;
  199. X        iHfMethods->methods = HfMethods;
  200. X        iHfMethods->name = HfName;
  201. X        iHfMethods->intersect = HfIntersect;
  202. X        iHfMethods->normal = HfNormal;
  203. X        iHfMethods->uv = HfUV;
  204. X        iHfMethods->bounds = HfBounds;
  205. X        iHfMethods->stats = HfStats;
  206. X        iHfMethods->checkbounds = TRUE;
  207. X        iHfMethods->closed = FALSE;
  208. X    }
  209. X    return iHfMethods;
  210. X}
  211. X
  212. X/*
  213. X * Intersect ray with height field.
  214. X */
  215. Xint
  216. XHfIntersect(hf, ray, mindist, maxdist)
  217. XHf *hf;
  218. XRay *ray;
  219. XFloat mindist, *maxdist;
  220. X{
  221. X    Vector hitpos;
  222. X    Float offset;
  223. X    Trav2D trav;
  224. X
  225. X    HFTests++;
  226. X
  227. X    /*
  228. X     * Find where we hit the hf cube.
  229. X     */
  230. X    VecAddScaled(ray->pos, mindist, ray->dir, &hitpos);
  231. X    if (OutOfBounds(&hitpos, hf->boundbox)) {
  232. X        offset = *maxdist;
  233. X        if (!BoundsIntersect(ray, hf->boundbox, mindist, &offset))
  234. X            return FALSE;
  235. X        hitpos.x = ray->pos.x + ray->dir.x * offset;
  236. X        hitpos.y = ray->pos.y + ray->dir.y * offset;
  237. X        hitpos.z = ray->pos.z + ray->dir.z * offset;
  238. X    } else
  239. X        hitpos = ray->pos;
  240. X    /*
  241. X     * Find out in which cell "hitpoint" is.
  242. X     */
  243. X    if (equal(hitpos.x, 1.))
  244. X        hitpos.x -= EPSILON;
  245. X    if (equal(hitpos.y, 1.))
  246. X        hitpos.y -= EPSILON;
  247. X
  248. X    if (ray->dir.x < 0.) {
  249. X        trav.stepX = trav.outX = -1;
  250. X        trav.tDX = -1. / (ray->dir.x * hf->spacing[hf->levels -1]);
  251. X    } else if (ray->dir.x > 0.) {
  252. X        trav.stepX = 1;
  253. X        trav.outX = hf->lsize[hf->levels -1];
  254. X        /*
  255. X         * (1./size) / ray
  256. X         */
  257. X        trav.tDX = 1. / (ray->dir.x * hf->spacing[hf->levels -1]);
  258. X    }
  259. X
  260. X    if (ray->dir.y < 0.) {
  261. X        trav.stepY = trav.outY = -1;
  262. X        trav.tDY = -1. / (ray->dir.y * hf->spacing[hf->levels -1]);
  263. X    } else if (ray->dir.y > 0.) {
  264. X        trav.stepY = 1;
  265. X        trav.outY = hf->lsize[hf->levels -1];
  266. X        trav.tDY = 1. / (ray->dir.y * hf->spacing[hf->levels -1]);
  267. X    }
  268. X
  269. X    trav.pDX.x = ray->dir.x * trav.tDX;
  270. X    trav.pDX.y = ray->dir.y * trav.tDX;
  271. X    trav.pDX.z = ray->dir.z * trav.tDX;
  272. X    trav.pDY.x = ray->dir.x * trav.tDY;
  273. X    trav.pDY.y = ray->dir.y * trav.tDY;
  274. X    trav.pDY.z = ray->dir.z * trav.tDY;
  275. X
  276. X    trav.cp = hitpos;
  277. X    trav.minz = hf->minz;
  278. X    trav.maxz = hf->maxz;
  279. X    if (DDA2D(hf, &ray->pos, &ray->dir, hf->levels -1, &trav, maxdist)) {
  280. X        HFHits++;
  281. X        return TRUE;
  282. X    }
  283. X    return FALSE;
  284. X}
  285. X
  286. X/*
  287. X * Traverse the grid using a modified DDA algorithm.  If the extent of
  288. X * the ray over a cell intersects the bounding volume defined by the
  289. X * four corners of the cell, either recurse or perform ray/surface
  290. X * intersection test.
  291. X */
  292. Xstatic int
  293. XDDA2D(hf, pos, ray, level, trav, maxdist)
  294. XHf *hf;
  295. XVector *pos, *ray;
  296. Xint level;
  297. XTrav2D *trav;
  298. XFloat *maxdist;
  299. X{
  300. X    int x, y, size, posZ;
  301. X    float **boundsmin, **boundsmax, spacing;
  302. X    Float tX, tY;
  303. X    Trav2D newtrav;
  304. X    Vector nxp, nyp;
  305. X
  306. X    size = hf->lsize[level];
  307. X    spacing = hf->spacing[level];
  308. X
  309. X    posZ = (ray->z > 0.);
  310. X
  311. X    x = trav->cp.x * hf->spacing[level];
  312. X    if (x == size)
  313. X        x--;
  314. X    y = trav->cp.y * hf->spacing[level];
  315. X    if (y == size)
  316. X        y--;
  317. X    boundsmax = hf->boundsmax[level];
  318. X    boundsmin = hf->boundsmin[level];
  319. X
  320. X    if (trav->outX > size) trav->outX = size;
  321. X    if (trav->outY > size) trav->outY = size;
  322. X    if (trav->outX < 0) trav->outX = -1;
  323. X    if (trav->outY < 0) trav->outY = -1;
  324. X
  325. X    if (ray->x < 0.) {
  326. X        tX = (x /spacing - trav->cp.x) / ray->x;
  327. X    } else if (ray->x > 0.)
  328. X        tX = ((x+1)/spacing - trav->cp.x) / ray->x;
  329. X    else
  330. X        tX = FAR_AWAY;
  331. X    if (ray->y < 0.) {
  332. X        tY = (y /spacing - trav->cp.y) / ray->y;
  333. X    } else if (ray->y > 0.)
  334. X        tY = ((y+1)/spacing - trav->cp.y) / ray->y;
  335. X    else
  336. X        tY = FAR_AWAY;
  337. X
  338. X    nxp.x = trav->cp.x + tX * ray->x;
  339. X    nxp.y = trav->cp.y + tX * ray->y;
  340. X    nxp.z = trav->cp.z + tX * ray->z;
  341. X
  342. X    nyp.x = trav->cp.x + tY * ray->x;
  343. X    nyp.y = trav->cp.y + tY * ray->y;
  344. X    nyp.z = trav->cp.z + tY * ray->z;
  345. X
  346. X    do {
  347. X        if (tX < tY) {
  348. X            if ((posZ && trav->cp.z <= boundsmax[y][x] &&
  349. X                 nxp.z >= boundsmin[y][x]) ||
  350. X                (!posZ && trav->cp.z >= boundsmin[y][x] &&
  351. X                 nxp.z <= boundsmax[y][x])) {
  352. X                if (level) {
  353. X                    /*
  354. X                     * Recurse -- compute constants
  355. X                     * needed for next level.
  356. X                     * Nicely enough, this just
  357. X                     * involves a few multiplications.
  358. X                     */
  359. X                    newtrav = *trav;
  360. X                    newtrav.tDX *= hf->iBestSize;
  361. X                    newtrav.tDY *= hf->iBestSize;
  362. X                    newtrav.maxz = boundsmax[y][x];
  363. X                    newtrav.minz = boundsmin[y][x];
  364. X                    if (ray->x < 0.)
  365. X                        newtrav.outX=hf->BestSize*x-1;
  366. X                    else
  367. X                        newtrav.outX=hf->BestSize*(x+1);
  368. X                    if (ray->y < 0.)
  369. X                        newtrav.outY=hf->BestSize*y-1;
  370. X                    else
  371. X                        newtrav.outY=hf->BestSize*(y+1);
  372. X                    newtrav.pDX.x *= hf->iBestSize;
  373. X                    newtrav.pDX.y *= hf->iBestSize;
  374. X                    newtrav.pDX.z *= hf->iBestSize;
  375. X                    newtrav.pDY.x *= hf->iBestSize;
  376. X                    newtrav.pDY.y *= hf->iBestSize;
  377. X                    newtrav.pDY.z *= hf->iBestSize;
  378. X                    if (DDA2D(hf,pos,ray,level-1,
  379. X                        &newtrav, maxdist))
  380. X                        return TRUE;
  381. X                } else if (CheckCell(x,y,hf,ray,pos,maxdist))
  382. X                    return TRUE;
  383. X            }
  384. X            x += trav->stepX;        /* Move in X */
  385. X            if (*maxdist < tX || x == trav->outX)
  386. X                /* If outside, quit */
  387. X                return FALSE;
  388. X            tX += trav->tDX;    /* Update position on ray */
  389. X            trav->cp = nxp;        /* cur pos gets next pos */
  390. X            nxp.x += trav->pDX.x;    /* Compute next pos */
  391. X            nxp.y += trav->pDX.y;
  392. X            nxp.z += trav->pDX.z;
  393. X        } else {
  394. X            if ((posZ && trav->cp.z <= boundsmax[y][x] &&
  395. X                 nyp.z >= boundsmin[y][x]) ||
  396. X                (!posZ && trav->cp.z >= boundsmin[y][x] &&
  397. X                 nyp.z <= boundsmax[y][x])) {
  398. X                if (level) {
  399. X                    /* Recurse */
  400. X                    newtrav = *trav;
  401. X                    newtrav.tDX *= hf->iBestSize;
  402. X                    newtrav.tDY *= hf->iBestSize;
  403. X                    newtrav.maxz = boundsmax[y][x];
  404. X                    newtrav.minz = boundsmin[y][x];
  405. X                    if (ray->x < 0.)
  406. X                        newtrav.outX=hf->BestSize*x-1;
  407. X                    else
  408. X                        newtrav.outX=hf->BestSize*(x+1);
  409. X                    if (ray->y < 0.)
  410. X                        newtrav.outY=hf->BestSize*y-1;
  411. X                    else
  412. X                        newtrav.outY=hf->BestSize*(y+1);
  413. X                    newtrav.pDX.x *= hf->iBestSize;
  414. X                    newtrav.pDX.y *= hf->iBestSize;
  415. X                    newtrav.pDX.z *= hf->iBestSize;
  416. X                    newtrav.pDY.x *= hf->iBestSize;
  417. X                    newtrav.pDY.y *= hf->iBestSize;
  418. X                    newtrav.pDY.z *= hf->iBestSize;
  419. X                    if (DDA2D(hf,pos,ray,level-1,
  420. X                            &newtrav, maxdist))
  421. X                        return TRUE;
  422. X                } else if (CheckCell(x,y,hf,ray,pos,maxdist))
  423. X                    return TRUE;
  424. X            }
  425. X            y += trav->stepY;
  426. X            if (*maxdist < tY || y == trav->outY)
  427. X                return FALSE;
  428. X            tY += trav->tDY;
  429. X            trav->cp = nyp;
  430. X            nyp.x += trav->pDY.x;
  431. X            nyp.y += trav->pDY.y;
  432. X            nyp.z += trav->pDY.z;
  433. X        }
  434. X    } while ((trav->cp.x <= 1. && trav->cp.y <= 1.) &&
  435. X         ((posZ && trav->cp.z <= trav->maxz) ||
  436. X         (!posZ && trav->cp.z >= trav->minz)));
  437. X
  438. X        /*
  439. X         * while ((we're inside the horizontal bounding box)
  440. X         *        (usually caught by outX & outY, but
  441. X         *         it's possible to go "too far" due to
  442. X         *         the fact that our levels of grids do
  443. X         *         not "nest" exactly if gridsize%BestSize != 0)
  444. X         *      and
  445. X         *      ((if ray->z is positive and we haven't gone through
  446. X         *       the upper bounding plane) or
  447. X         *      (if ray->z is negative and we haven't gone through
  448. X         *       the lower bounding plane)));
  449. X         */
  450. X
  451. X    return FALSE;
  452. X}
  453. X
  454. X/*
  455. X * Check for ray/cell intersection
  456. X */
  457. Xstatic int
  458. XCheckCell(x, y, hf, ray, pos, maxdist)
  459. Xint x, y;
  460. XHf *hf;
  461. XVector *ray, *pos;
  462. XFloat *maxdist;
  463. X{
  464. X    hfTri *tri1, *tri2;
  465. X    Float d1, d2;
  466. X
  467. X    d1 = d2 = FAR_AWAY;
  468. X
  469. X    if (tri1 = CreateHfTriangle(hf, x, y, x+1, y, x, y+1, TRI1))
  470. X        d1 = intHftri(ray, pos, tri1);
  471. X    if (tri2 = CreateHfTriangle(hf, x+1, y, x+1, y+1, x, y+1, TRI2))
  472. X        d2 = intHftri(ray, pos, tri2);
  473. X
  474. X    if (d1 == FAR_AWAY && d2 == FAR_AWAY)
  475. X        return FALSE;
  476. X
  477. X    if (d1 < d2) {
  478. X        if (d1 < *maxdist) {
  479. X            hf->hittri = *tri1;
  480. X            *maxdist = d1;
  481. X            return TRUE;
  482. X        }
  483. X        return FALSE;
  484. X    }
  485. X
  486. X    if (d2 < *maxdist) {
  487. X        hf->hittri = *tri2;
  488. X        *maxdist = d2;
  489. X        return TRUE;
  490. X    }
  491. X    return FALSE;
  492. X}
  493. X
  494. Xstatic hfTri *
  495. XCreateHfTriangle(hf, x1, y1, x2, y2, x3, y3, which)
  496. XHf *hf;
  497. Xint x1, y1, x2, y2, x3, y3, which;
  498. X{
  499. X    hfTri *tri;
  500. X    Float xid, yid;
  501. X    Vector tmp1, tmp2;
  502. X
  503. X    /*
  504. X     * Don't use triangles with "unset" vertices.
  505. X     */
  506. X    if (hf->data[y1][x1] == HF_UNSET ||
  507. X        hf->data[y2][x2] == HF_UNSET ||
  508. X        hf->data[y3][x3] == HF_UNSET)
  509. X        return (hfTri *)0;
  510. X
  511. X    xid = (Float)x1 / (Float)(hf->size -1);
  512. X    yid = (Float)y1 / (Float)(hf->size -1);
  513. X
  514. X    if ((tri = GetQueuedTri(hf, xid, yid, which)) != (hfTri *)0)
  515. X        return tri;
  516. X
  517. X    tri = (hfTri *)Malloc(sizeof(hfTri));
  518. X
  519. X    tri->type = which;
  520. X    tri->v1.x = xid;
  521. X    tri->v1.y = yid;
  522. X    tri->v1.z = hf->data[y1][x1];
  523. X    tri->v2.x = (Float)x2 / (Float)(hf->size-1);
  524. X    tri->v2.y = (Float)y2 / (Float)(hf->size-1);
  525. X    tri->v2.z = hf->data[y2][x2];
  526. X    tri->v3.x = (Float)x3 / (Float)(hf->size-1);
  527. X    tri->v3.y = (Float)y3 / (Float)(hf->size-1);
  528. X    tri->v3.z = hf->data[y3][x3];
  529. X
  530. X    tmp1.x = tri->v2.x - tri->v1.x;
  531. X    tmp1.y = tri->v2.y - tri->v1.y;
  532. X    tmp1.z = tri->v2.z - tri->v1.z;
  533. X    tmp2.x = tri->v3.x - tri->v1.x;
  534. X    tmp2.y = tri->v3.y - tri->v1.y;
  535. X    tmp2.z = tri->v3.z - tri->v1.z;
  536. X
  537. X    (void)VecNormCross(&tmp1, &tmp2, &tri->norm);
  538. X
  539. X    tri->d = -dotp(&tri->v1, &tri->norm);
  540. X
  541. X    QueueTri(hf, tri);
  542. X    return tri;
  543. X}
  544. X
  545. X/*
  546. X * Intersect ray with right isoscoles triangle, the hypotenuse of which
  547. X * has slope of -1.
  548. X */
  549. Xstatic Float
  550. XintHftri(ray, pos, tri)
  551. XhfTri *tri;
  552. XVector *pos, *ray;
  553. X{
  554. X    Float u, v, dist, xpos, ypos;
  555. X
  556. X    u = dotp(&tri->norm, pos) + tri->d;
  557. X    v = dotp(&tri->norm, ray);
  558. X
  559. X    if ((u <= 0. || v > -EPSILON) && (u >= 0. && v < EPSILON))
  560. X        return FAR_AWAY;
  561. X
  562. X    dist = -u / v;
  563. X
  564. X    if (dist < EPSILON)
  565. X        return FAR_AWAY;
  566. X
  567. X    xpos = pos->x + dist*ray->x;
  568. X    ypos = pos->y + dist*ray->y;
  569. X
  570. X    if (tri->type == TRI1 && xpos >= tri->v1.x && ypos >= tri->v1.y &&
  571. X        xpos + ypos <= tri->v2.x + tri->v2.y)
  572. X        return dist;
  573. X    if (tri->type == TRI2 && xpos <= tri->v2.x && ypos <= tri->v2.y &&
  574. X        xpos + ypos >= tri->v1.x + tri->v1.y)
  575. X        return dist;
  576. X    return FAR_AWAY;
  577. X}
  578. X
  579. X/*
  580. X * Compute normal to height field.
  581. X */
  582. X/*ARGSUSED*/
  583. Xint
  584. XHfNormal(hf, pos, nrm, gnrm)
  585. XHf *hf;
  586. XVector *pos, *nrm, *gnrm;
  587. X{
  588. X    *gnrm = *nrm = hf->hittri.norm;
  589. X    return FALSE;
  590. X}
  591. X
  592. X/*ARGSUSED*/
  593. Xvoid
  594. XHfUV(hf, pos, norm, uv, dpdu, dpdv)
  595. XHf *hf;
  596. XVector *pos, *norm, *dpdu, *dpdv;
  597. XVec2d *uv;
  598. X{
  599. X    uv->u = pos->x;
  600. X    uv->v = pos->y;
  601. X    if (dpdu) {
  602. X        dpdu->x = 1.;
  603. X        dpdu->y = dpdv->z = 0.;
  604. X        dpdv->x = dpdv->z = 0.;
  605. X        dpdv->y = 1.;
  606. X    }
  607. X}
  608. X
  609. X/*
  610. X * Compute heightfield bounding box.
  611. X */
  612. Xvoid
  613. XHfBounds(hf, bounds)
  614. XHf *hf;
  615. XFloat bounds[2][3];
  616. X{
  617. X    /*
  618. X     * By default, height fields are centered at (0.5, 0.5, 0.)
  619. X     */
  620. X    bounds[LOW][X] = bounds[LOW][Y] = 0;
  621. X    bounds[HIGH][X] = bounds[HIGH][Y] = 1;
  622. X    bounds[LOW][Z] = hf->minz;
  623. X    bounds[HIGH][Z] = hf->maxz;
  624. X}
  625. X
  626. X/*
  627. X * Build min/max altitude value arrays for the given grid level.
  628. X */
  629. Xstatic void
  630. Xintegrate_grid(hf, level)
  631. XHf *hf;
  632. Xint level;
  633. X{
  634. X    int i, j, k, l, ii, ji;
  635. X    float max_alt, min_alt;
  636. X    float **maxinto, **mininto, **frommax, **frommin, *minptr, *maxptr;
  637. X    int insize, fromsize, fact;
  638. X
  639. X    maxinto = hf->boundsmax[level];
  640. X    mininto = hf->boundsmin[level];
  641. X    insize = hf->lsize[level];
  642. X    frommax = hf->boundsmax[level-1];
  643. X    frommin = hf->boundsmin[level-1];
  644. X    fact = hf->BestSize;
  645. X    fromsize = hf->lsize[level-1];
  646. X
  647. X    ii = 0;
  648. X
  649. X    for (i = 0; i < insize; i++) {
  650. X        ji = 0;
  651. X        for (j = 0; j < insize; j++) {
  652. X            max_alt = HF_UNSET;
  653. X            min_alt = -HF_UNSET;
  654. X            for (k = 0; k <= fact; k++) {
  655. X                if (ii+k >= fromsize)
  656. X                    continue;
  657. X                maxptr = &frommax[ii+k][ji];
  658. X                minptr = &frommin[ii+k][ji];
  659. X                for (l = 0; l <= fact; l++,maxptr++,minptr++) {
  660. X                    if (ji+l >= fromsize)
  661. X                        continue;
  662. X                    if (*maxptr > max_alt)
  663. X                        max_alt = *maxptr;
  664. X                    if (*minptr < min_alt)
  665. X                        min_alt = *minptr;
  666. X                }
  667. X            }
  668. X            maxinto[i][j] = max_alt + EPSILON;
  669. X            mininto[i][j] = min_alt - EPSILON;
  670. X            ji += fact;
  671. X        }
  672. X        ii += fact;
  673. X    }
  674. X}
  675. X
  676. X/*
  677. X * Place the given triangle in the triangle cache.
  678. X */
  679. Xstatic void
  680. XQueueTri(hf, tri)
  681. XHf *hf;
  682. XhfTri *tri;
  683. X{
  684. X    if (hf->q[hf->qtail])        /* Free old triangle data */
  685. X        free((voidstar)hf->q[hf->qtail]);
  686. X    hf->q[hf->qtail] = tri;        /* Put on tail */
  687. X    hf->qtail = (hf->qtail + 1) % hf->qsize;    /* Increment tail */
  688. X}
  689. X
  690. X/*
  691. X * Search list of cached trianges to see if this triangle has been
  692. X * cached.  If so, return a pointer to it.  If not, return null pointer.
  693. X */
  694. Xstatic hfTri *
  695. XGetQueuedTri(hf, x, y, flag)
  696. XHf *hf;
  697. XFloat x, y;
  698. Xint flag;
  699. X{
  700. X    register int i;
  701. X    register hfTri **tmp;
  702. X
  703. X    for (i = 0, tmp = hf->q; i < hf->qsize; i++, tmp++) {
  704. X        if (*tmp && (*tmp)->v1.x == x && (*tmp)->v1.y == y &&
  705. X            (*tmp)->type == flag)
  706. X            return *tmp;    /* vertices & flag match, return it */
  707. X    }
  708. X
  709. X    return (hfTri *)0;
  710. X}
  711. X
  712. X/*
  713. X * Return maximum height of cell indexed by y,x.  This could be done
  714. X * as a macro, but many C compliers will choke on it.
  715. X */
  716. Xstatic float
  717. Xminalt(y,x,data)
  718. Xint x, y;
  719. Xfloat **data;
  720. X{
  721. X    float  min_alt;
  722. X
  723. X    min_alt = min(data[y][x], data[y+1][x]);
  724. X    min_alt = min(min_alt, data[y][x+1]);
  725. X    min_alt = min(min_alt, data[y+1][x+1]);
  726. X    return min_alt;
  727. X}
  728. X
  729. X/*
  730. X * Return maximum cell height, as above.
  731. X */
  732. Xstatic float
  733. Xmaxalt(y,x,data)
  734. Xint x, y;
  735. Xfloat **data;
  736. X{
  737. X    float  max_alt;
  738. X
  739. X    max_alt = max(data[y][x], data[y+1][x]);
  740. X    max_alt = max(max_alt, data[y][x+1]);
  741. X    max_alt = max(max_alt, data[y+1][x+1]);
  742. X    return max_alt;
  743. X}
  744. X
  745. Xchar *
  746. XHfName()
  747. X{
  748. X    return hfName;
  749. X}
  750. X
  751. Xvoid
  752. XHfStats(tests, hits)
  753. Xunsigned long *tests, *hits;
  754. X{
  755. X    *tests = HFTests;
  756. X    *hits = HFHits;
  757. X}
  758. X
  759. Xvoid
  760. XHfMethodRegister(meth)
  761. XUserMethodType meth;
  762. X{
  763. X    if (iHfMethods)
  764. X        iHfMethods->user = meth;
  765. X}
  766. END_OF_FILE
  767. if test 17333 -ne `wc -c <'libray/libobj/hf.c'`; then
  768.     echo shar: \"'libray/libobj/hf.c'\" unpacked with wrong size!
  769. fi
  770. # end of 'libray/libobj/hf.c'
  771. fi
  772. if test -f 'raypaint/render.c' -a "${1}" != "-c" ; then 
  773.   echo shar: Will not clobber existing file \"'raypaint/render.c'\"
  774. else
  775. echo shar: Extracting \"'raypaint/render.c'\" \(17440 characters\)
  776. sed "s/^X//" >'raypaint/render.c' <<'END_OF_FILE'
  777. X/*
  778. X * render.c
  779. X *
  780. X * Copyright (C) 1989, 1991 Craig E. Kolb, Rod G. Bogart
  781. X *
  782. X * This software may be freely copied, modified, and redistributed
  783. X * provided that this copyright notice is preserved on all copies.
  784. X *
  785. X * You may not distribute this software, in whole or in part, as part of
  786. X * any commercial product without the express consent of the authors.
  787. X * 
  788. X * There is no warranty or other guarantee of fitness of this software
  789. X * for any purpose.  It is provided solely "as is".
  790. X *
  791. X * $Id: render.c,v 4.0 91/07/17 17:37:09 kolb Exp Locker: kolb $
  792. X *
  793. X * $Log:    render.c,v $
  794. X * Revision 4.0  91/07/17  17:37:09  kolb
  795. X * Initial version.
  796. X * 
  797. X */
  798. X
  799. X#include "rayshade.h"
  800. X#include "libcommon/sampling.h"
  801. X#include "libsurf/atmosphere.h"
  802. X#include "viewing.h"
  803. X#include "options.h"
  804. X#include "stats.h"
  805. X#include "picture.h"
  806. X
  807. X/*
  808. X * "Dither matrices" used to encode the 'number' of a ray that samples a
  809. X * particular portion of a pixel.  Hand-coding is ugly, but...
  810. X */
  811. Xstatic int        *SampleNumbers;
  812. Xstatic int OneSample[1] =     {0};
  813. Xstatic int TwoSamples[4] =    {0, 2,
  814. X                 3, 1};
  815. Xstatic int ThreeSamples[9] =    {0, 2, 7,
  816. X                 6, 5, 1,
  817. X                 3, 8, 4};
  818. Xstatic int FourSamples[16] =    { 0,  8,  2, 10,
  819. X                 12,  4, 14,  6,
  820. X                  3, 11,  1,  9,
  821. X                 15,  7, 13,  5};
  822. Xstatic int FiveSamples[25] =    { 0,  8, 23, 17,  2,
  823. X                 19, 12,  4, 20, 15,
  824. X                  3, 21, 16,  9,  6,
  825. X                 14, 10, 24,  1, 13,
  826. X                 22,  7, 18, 11,  5};
  827. Xstatic int SixSamples[36] =    { 6, 32,  3, 34, 35,  1,
  828. X                  7, 11, 27, 28,  8, 30,
  829. X                 24, 14, 16, 15, 23, 19,
  830. X                 13, 20, 22, 21, 17, 18,
  831. X                 25, 29, 10,  9, 26, 12,
  832. X                 36,  5, 33,  4,  2, 31};
  833. Xstatic int SevenSamples[49] =    {22, 47, 16, 41, 10, 35,  4,
  834. X                  5, 23, 48, 17, 42, 11, 29,
  835. X                 30,  6, 24, 49, 18, 36, 12,
  836. X                 13, 31,  7, 25, 43, 19, 37,
  837. X                 38, 14, 32,  1, 26, 44, 20,
  838. X                 21, 39,  8, 33,  2, 27, 45,
  839. X                 46, 15, 40,  9, 34,  3, 28};
  840. Xstatic int EightSamples[64] =    { 8, 58, 59,  5,  4, 62, 63,  1,
  841. X                 49, 15, 14, 52, 53, 11, 10, 56,
  842. X                 41, 23, 22, 44, 45, 19, 18, 48,
  843. X                 32, 34, 35, 29, 28, 38, 39, 25,
  844. X                 40, 26, 27, 37, 36, 30, 31, 33,
  845. X                 17, 47, 46, 20, 21, 43, 42, 24,
  846. X                  9, 55, 54, 12, 13, 51, 50, 16,
  847. X                 64,  2,  3, 61, 60,  6,  7, 57};
  848. X#define RFAC    0.299
  849. X#define GFAC    0.587
  850. X#define BFAC    0.114
  851. X
  852. X#define NOT_CLOSED 0
  853. X#define CLOSED_PARTIAL   1
  854. X#define CLOSED_SUPER 2
  855. X/*
  856. X * If a region has area < MINAREA, it is considered to be "closed",
  857. X * (a permanent leaf).  Regions that meet this criterion
  858. X * are drawn pixel-by-pixel rather.
  859. X */
  860. X#define MINAREA        9
  861. X
  862. X#define SQ_AREA(s)    (((s)->xsize+1)*((s)->ysize+1))
  863. X
  864. X#define PRIORITY(s)    ((s)->var * SQ_AREA(s))
  865. X
  866. X#define INTENSITY(p)    ((RFAC*(p)[0] + GFAC*(p)[1] + BFAC*(p)[2])/255.)
  867. X
  868. X#define OVERLAPS_RECT(s) (!Rectmode || \
  869. X                ((s)->xpos <= Rectx1 && \
  870. X                 (s)->ypos <= Recty1 && \
  871. X                 (s)->xpos+(s)->xsize >= Rectx0 && \
  872. X                 (s)->ypos+(s)->ysize >= Recty0))
  873. X
  874. Xtypedef unsigned char RGB[3];
  875. X
  876. Xstatic RGB **Image;
  877. Xstatic char **SampleMap;
  878. X
  879. X/*
  880. X * Sample square
  881. X */
  882. Xtypedef struct SSquare {
  883. X    short xpos, ypos, xsize, ysize;
  884. X    short depth;
  885. X    short leaf, closed;
  886. X    float mean, var, prio;
  887. X    struct SSquare *child[4], *parent;
  888. X} SSquare;
  889. X
  890. XSSquare *SSquares, *SSquareCreate(), *SSquareInstall(), *SSquareSelect(),
  891. X    *SSquareFetchAtMouse();
  892. X
  893. XFloat SSquareComputeLeafVar();
  894. X
  895. XRay    TopRay;                /* Top-level ray. */
  896. Xint    Rectmode = FALSE,
  897. X    Rectx0, Recty0, Rectx1, Recty1;
  898. Xint    SuperSampleMode = 0;
  899. X#define SSCLOSED (SuperSampleMode + 1)
  900. X
  901. XRender(argc, argv)
  902. Xint argc;
  903. Xchar **argv;
  904. X{
  905. X    /*
  906. X     * Do an adaptive trace, displaying results in a
  907. X     * window as we go.
  908. X     */
  909. X    SSquare *cursq;
  910. X    Pixel *pixelbuf;
  911. X    int y, x;
  912. X
  913. X    /*
  914. X     * The top-level ray TopRay always has as its origin the
  915. X     * eye position and as its medium NULL, indicating that it
  916. X     * is passing through a medium with index of refraction
  917. X     * equal to DefIndex.
  918. X     */
  919. X    TopRay.pos = Camera.pos;
  920. X    TopRay.media = (Medium *)0;
  921. X    TopRay.depth = 0;
  922. X    /*
  923. X     * Doesn't handle motion blur yet.
  924. X     */
  925. X    TopRay.time = Options.framestart;
  926. X
  927. X    GraphicsInit(Screen.xsize, Screen.ysize, "rayview");
  928. X    /*
  929. X     * Allocate array of samples.
  930. X     */
  931. X    Image=(RGB **)Malloc(Screen.ysize*sizeof(RGB *));
  932. X    SampleMap = (char **)Malloc(Screen.ysize * sizeof(char *));
  933. X    for (y = 0; y < Screen.ysize; y++) {
  934. X        Image[y]=(RGB *)Calloc(Screen.xsize, sizeof(RGB));
  935. X        SampleMap[y] = (char *)Calloc(Screen.xsize, sizeof(char));
  936. X    }
  937. X    switch (Sampling.sidesamples) {
  938. X        case 1:
  939. X            SampleNumbers = OneSample;
  940. X            break;
  941. X        case 2:
  942. X            SampleNumbers = TwoSamples;
  943. X            break;
  944. X        case 3:
  945. X            SampleNumbers = ThreeSamples;
  946. X            break;
  947. X        case 4:
  948. X            SampleNumbers = FourSamples;
  949. X            break;
  950. X        case 5:
  951. X            SampleNumbers = FiveSamples;
  952. X            break;
  953. X        case 6:
  954. X            SampleNumbers = SixSamples;
  955. X            break;
  956. X        case 7:
  957. X            SampleNumbers = SevenSamples;
  958. X            break;
  959. X        case 8:
  960. X            SampleNumbers = EightSamples;
  961. X            break;
  962. X        default:
  963. X            RLerror(RL_PANIC,
  964. X                "Sorry, %d rays/pixel not supported.\n",
  965. X                    Sampling.totsamples);
  966. X    }
  967. X    /*
  968. X     * Take initial four samples
  969. X     */
  970. X    SSquareSample(0, 0, FALSE);
  971. X    SSquareSample(Screen.xsize -1, 0, FALSE);
  972. X    SSquareSample(Screen.xsize -1, Screen.ysize -1, FALSE);
  973. X    SSquareSample(0, Screen.ysize -1, FALSE);
  974. X    SSquares = SSquareInstall(0, 0, Screen.xsize -1, Screen.ysize -1,
  975. X                  0, (SSquare *) NULL);
  976. X
  977. X    for (y = 1; y <= 5; y++) {
  978. X        /*
  979. X         * Create and draw every region at depth y.
  980. X         */
  981. X        SSquareDivideToDepth(SSquares, y);
  982. X    }
  983. X        
  984. X
  985. X    while ((cursq = SSquareSelect(SSquares)) != (SSquare *)NULL) {
  986. X        SSquareDivide(cursq);
  987. X        if (GraphicsRedraw())
  988. X            SSquareDraw(SSquares);
  989. X        if (GraphicsMiddleMouseEvent()) 
  990. X            SSetRectMode();
  991. X    }
  992. X
  993. X    /*
  994. X     * Finished the image; write to image file.
  995. X     */
  996. X    pixelbuf = (Pixel *)Malloc(Screen.xsize * sizeof(Pixel));
  997. X    PictureStart(argv);
  998. X    for (y = 0; y < Screen.ysize; y++) {
  999. X        /*
  1000. X         * This is really disgusting.
  1001. X         */
  1002. X        for (x = 0; x < Screen.xsize; x++) {
  1003. X            pixelbuf[x].r = Image[y][x][0] / 255.;
  1004. X            pixelbuf[x].g = Image[y][x][1] / 255.;
  1005. X            pixelbuf[x].b = Image[y][x][2] / 255.;
  1006. X            pixelbuf[x].alpha = SampleMap[y][x];
  1007. X        }
  1008. X        PictureWriteLine(pixelbuf);
  1009. X    }
  1010. X    PictureEnd();
  1011. X    free((voidstar)pixelbuf);
  1012. X}
  1013. X
  1014. XFloat
  1015. XSampleTime(sampnum)
  1016. Xint sampnum;
  1017. X{
  1018. X    Float window, jitter = 0.0, res;
  1019. X
  1020. X    if (Options.shutterspeed <= 0.)
  1021. X        return Options.framestart;
  1022. X    if (Options.jitter)
  1023. X        jitter = nrand();
  1024. X    window = Options.shutterspeed / Sampling.totsamples;
  1025. X    res = Options.framestart + window * (sampnum + jitter);
  1026. X    TimeSet(res);
  1027. X    return res;
  1028. X}
  1029. X
  1030. XSSetRectMode()
  1031. X{
  1032. X    int x1,y1,x2,y2;
  1033. X
  1034. X    if (Rectmode) {
  1035. X        Rectmode = FALSE;
  1036. X        RecomputePriority(SSquares);
  1037. X    }
  1038. X    GraphicsGetMousePos(&x1, &y1);
  1039. X    while (GraphicsMiddleMouseEvent())
  1040. X        ;
  1041. X    GraphicsGetMousePos(&x2, &y2);
  1042. X    if (x1 < x2) {
  1043. X        Rectx0 = (x1 < 0) ? 0 : x1;
  1044. X        Rectx1 = (x2 >= Screen.xsize) ? Screen.xsize - 1 : x2;
  1045. X    } else {
  1046. X        Rectx0 = (x2 < 0) ? 0 : x2;
  1047. X        Rectx1 = (x1 >= Screen.xsize) ? Screen.xsize - 1 : x1;
  1048. X    } if (y1 < y2) {
  1049. X        Recty0 = (y1 < 0) ? 0 : y1;
  1050. X        Recty1 = (y2 >= Screen.ysize) ? Screen.ysize - 1 : y2;
  1051. X    } else {
  1052. X        Recty0 = (y2 < 0) ? 0 : y2;
  1053. X        Recty1 = (y1 >= Screen.ysize) ? Screen.ysize - 1 : y1;
  1054. X    }
  1055. X    Rectmode = TRUE;
  1056. X    /* setup current rect */
  1057. X    (void)RecomputePriority(SSquares);
  1058. X}
  1059. X
  1060. XRecomputePriority(sq)
  1061. XSSquare *sq;
  1062. X{
  1063. X    Float maxp;
  1064. X
  1065. X    if (!OVERLAPS_RECT(sq)) {
  1066. X        sq->closed = SSCLOSED;
  1067. X        return FALSE;
  1068. X    }
  1069. X
  1070. X    if (sq->leaf) {
  1071. X        if (SQ_AREA(sq) >= MINAREA)
  1072. X            sq->closed = NOT_CLOSED;
  1073. X        return TRUE;
  1074. X    }
  1075. X    maxp = 0.;
  1076. X    if (RecomputePriority(sq->child[0]))
  1077. X        maxp = max(maxp, sq->child[0]->prio);
  1078. X    if (RecomputePriority(sq->child[1]))
  1079. X        maxp = max(maxp, sq->child[1]->prio);
  1080. X    if (RecomputePriority(sq->child[2]))
  1081. X        maxp = max(maxp, sq->child[2]->prio);
  1082. X    if (RecomputePriority(sq->child[3]))
  1083. X        maxp = max(maxp, sq->child[3]->prio);
  1084. X    sq->prio = maxp;
  1085. X#if 0
  1086. X    if ((sq->child[0]->closed == CLOSED_SUPER) &&
  1087. X        (sq->child[1]->closed == CLOSED_SUPER) &&
  1088. X        (sq->child[2]->closed == CLOSED_SUPER) &&
  1089. X        (sq->child[3]->closed == CLOSED_SUPER))
  1090. X        sq->closed = CLOSED_SUPER;
  1091. X    else if (sq->child[0]->closed && sq->child[1]->closed &&
  1092. X         sq->child[2]->closed && sq->child[3]->closed)
  1093. X        sq->closed = CLOSED_PARTIAL;
  1094. X    else
  1095. X        sq->closed = NOT_CLOSED;
  1096. X#else
  1097. X    if ((sq->child[0]->closed >= SSCLOSED) &&
  1098. X        (sq->child[1]->closed >= SSCLOSED) &&
  1099. X        (sq->child[2]->closed >= SSCLOSED) &&
  1100. X        (sq->child[3]->closed >= SSCLOSED))
  1101. X        sq->closed = SSCLOSED;
  1102. X    else
  1103. X        sq->closed = NOT_CLOSED;
  1104. X#endif
  1105. X    return TRUE;
  1106. X}
  1107. X
  1108. XSSquareSample(x, y, supersample)
  1109. Xint x, y, supersample;
  1110. X{
  1111. X    Float upos, vpos, u, v;
  1112. X    int xx, yy, xp, sampnum, usamp, vsamp;
  1113. X    Pixel ctmp;
  1114. X    Pixel p;
  1115. X    extern unsigned char correct();
  1116. X
  1117. X    if (SampleMap[y][x] >= 128 + supersample)
  1118. X        /* already a sample there */
  1119. X        return;
  1120. X    SampleMap[y][x] = 128 + supersample;
  1121. X    if (supersample) {
  1122. X        p.r = p.g = p.b = p.alpha = 0;
  1123. X        sampnum = 0;
  1124. X        xp = x + Screen.minx;
  1125. X        vpos = Screen.miny + y - 0.5*Sampling.filterwidth;
  1126. X        for (yy = 0; yy < Sampling.sidesamples; yy++,
  1127. X             vpos += Sampling.filterdelta) {
  1128. X            upos = xp - 0.5*Sampling.filterwidth;
  1129. X            for (xx = 0; xx < Sampling.sidesamples; xx++,
  1130. X                 upos += Sampling.filterdelta) {
  1131. X                if (Options.jitter) {
  1132. X                    u = upos + nrand()*Sampling.filterdelta;
  1133. X                    v = vpos + nrand()*Sampling.filterdelta;
  1134. X                } else {
  1135. X                    u = upos;
  1136. X                    v = vpos;
  1137. X                }
  1138. X                TopRay.time = SampleTime(SampleNumbers[sampnum]);
  1139. X                SampleScreen(u, v, &TopRay, &ctmp,
  1140. X                         SampleNumbers[sampnum]);
  1141. X                p.r += ctmp.r*Sampling.filter[xx][yy];
  1142. X                p.g += ctmp.g*Sampling.filter[xx][yy];
  1143. X                p.b += ctmp.b*Sampling.filter[xx][yy];
  1144. X                if (++sampnum == Sampling.totsamples)
  1145. X                    sampnum = 0;
  1146. X            }
  1147. X        }
  1148. X    }
  1149. X    else {
  1150. X        sampnum = nrand() * Sampling.totsamples;
  1151. X        usamp = sampnum % Sampling.sidesamples;
  1152. X        vsamp = sampnum / Sampling.sidesamples;
  1153. X
  1154. X        vpos = Screen.miny + y - 0.5*Sampling.filterwidth
  1155. X            + vsamp * Sampling.filterdelta;
  1156. X        upos = x + Screen.minx - 0.5*Sampling.filterwidth +
  1157. X                usamp*Sampling.filterdelta;
  1158. X        if (Options.jitter) {
  1159. X            vpos += nrand()*Sampling.filterdelta;
  1160. X            upos += nrand()*Sampling.filterdelta;
  1161. X        }
  1162. X        TopRay.time = SampleTime(SampleNumbers[sampnum]);
  1163. X        SampleScreen(upos, vpos, &TopRay, &p, SampleNumbers[sampnum]);
  1164. X    }
  1165. X    Image[y][x][0] = CORRECT(p.r);
  1166. X    Image[y][x][1] = CORRECT(p.g);
  1167. X    Image[y][x][2] = CORRECT(p.b);
  1168. X}
  1169. X
  1170. XSSquare *
  1171. XSSquareCreate(xp, yp, xs, ys, d, parent)
  1172. Xint xp, yp, xs, ys, d;
  1173. XSSquare *parent;
  1174. X{
  1175. X    SSquare *new;
  1176. X    Float i1, i2, i3, i4;
  1177. X
  1178. X    new = (SSquare *)Calloc(1, sizeof(SSquare));
  1179. X    new->xpos = xp; new->ypos = yp;
  1180. X    new->xsize = xs; new->ysize = ys;
  1181. X    new->depth = d;
  1182. X    new->parent = parent;
  1183. X    i1 = INTENSITY(Image[new->ypos][new->xpos]);
  1184. X    i2 = INTENSITY(Image[new->ypos+new->ysize][new->xpos]);
  1185. X    i3 = INTENSITY(Image[new->ypos+new->ysize][new->xpos+new->xsize]);
  1186. X    i4 = INTENSITY(Image[new->ypos][new->xpos+new->xsize]);
  1187. X    new->mean = 0.25 * (i1+i2+i3+i4);
  1188. X    if (SQ_AREA(new) < MINAREA) {
  1189. X        new->prio = 0;
  1190. X        new->closed = SSCLOSED;
  1191. X    } else {
  1192. X        new->var = SSquareComputeLeafVar(new, i1, i2, i3, i4);
  1193. X        new->prio = PRIORITY(new);
  1194. X        new->closed = NOT_CLOSED;
  1195. X    }
  1196. X    new->leaf = TRUE;
  1197. X    return new;
  1198. X}
  1199. X
  1200. XFloat
  1201. XSSquareComputeLeafVar(sq, i1, i2, i3, i4)
  1202. XSSquare *sq;
  1203. XFloat i1, i2, i3, i4;
  1204. X{
  1205. X    Float res, diff;
  1206. X
  1207. X    diff = i1 - sq->mean;
  1208. X    res = diff*diff;
  1209. X    diff = i2 - sq->mean;
  1210. X    res += diff*diff;
  1211. X    diff = i3 - sq->mean;
  1212. X    res += diff*diff;
  1213. X    diff = i4 - sq->mean;
  1214. X    return res + diff*diff;
  1215. X}
  1216. X
  1217. XSSquareDivideToDepth(sq, d)
  1218. XSSquare *sq;
  1219. Xint d;
  1220. X{
  1221. X    if (sq->depth == d)
  1222. X        return;
  1223. X    if (sq->leaf)
  1224. X        SSquareDivide(sq);
  1225. X    SSquareDivideToDepth(sq->child[0], d);
  1226. X    SSquareDivideToDepth(sq->child[1], d);
  1227. X    SSquareDivideToDepth(sq->child[2], d);
  1228. X    SSquareDivideToDepth(sq->child[3], d);
  1229. X}
  1230. X
  1231. XSSquareDivide(sq)
  1232. XSSquare *sq;
  1233. X{
  1234. X    int lowx, lowy, midx, midy, hix, hiy;
  1235. X    int newxsize, newysize, ndepth, supersample;
  1236. X    /*
  1237. X     * Divide the square into fourths by tracing 12
  1238. X     * new samples if necessary.
  1239. X     */
  1240. X    newxsize = sq->xsize / 2;
  1241. X    newysize = sq->ysize / 2;
  1242. X    lowx = sq->xpos; lowy = sq->ypos;
  1243. X    midx = sq->xpos + newxsize;
  1244. X    midy = sq->ypos + newysize;
  1245. X    hix  = sq->xpos + sq->xsize;
  1246. X    hiy  = sq->ypos + sq->ysize;
  1247. X    ndepth = sq->depth + 1;
  1248. X    /* create new samples */
  1249. X    supersample = FALSE;
  1250. X    SSquareSample(midx, lowy, supersample);
  1251. X    SSquareSample(lowx, midy, supersample);
  1252. X    SSquareSample(midx, midy, supersample);
  1253. X    SSquareSample(hix,  midy, supersample);
  1254. X    SSquareSample(midx, hiy, supersample);
  1255. X#ifdef SHARED_EDGES
  1256. X    /* create and draw new squares */
  1257. X    sq->child[0] = SSquareInstall(lowx,lowy,newxsize,newysize,ndepth,sq);
  1258. X    sq->child[1] = SSquareInstall(midx, lowy, sq->xsize - newxsize,
  1259. X            newysize, ndepth,sq);
  1260. X    sq->child[2] = SSquareInstall(lowx, midy, newxsize,
  1261. X            sq->ysize - newysize, ndepth,sq);
  1262. X    sq->child[3] = SSquareInstall(midx, midy, sq->xsize - newxsize,
  1263. X             sq->ysize - newysize, ndepth,sq);
  1264. X#else
  1265. X    /*
  1266. X     *  draw additional samples in order to subdivide such that
  1267. X     * edges of regions do not overlap
  1268. X     */
  1269. X    SSquareSample(midx +1, lowy, supersample);
  1270. X    SSquareSample(midx +1, midy, supersample);
  1271. X    SSquareSample(lowx, midy +1, supersample);
  1272. X    SSquareSample(midx, midy +1, supersample);
  1273. X    SSquareSample(midx +1, midy +1, supersample);
  1274. X    SSquareSample(hix, midy +1, supersample);
  1275. X    SSquareSample(midx +1, hiy, supersample);
  1276. X
  1277. X    /* create and draw new squares */
  1278. X    sq->child[0] = SSquareInstall(lowx,lowy,
  1279. X                newxsize,newysize,ndepth,sq);
  1280. X    sq->child[1] = SSquareInstall(midx+1, lowy, sq->xsize - newxsize -1,
  1281. X            newysize, ndepth,sq);
  1282. X    sq->child[2] = SSquareInstall(lowx, midy+1, newxsize,
  1283. X            sq->ysize - newysize-1, ndepth,sq);
  1284. X    sq->child[3] = SSquareInstall(midx+1, midy+1, sq->xsize - newxsize-1,
  1285. X             sq->ysize - newysize-1, ndepth,sq);
  1286. X#endif
  1287. X    sq->leaf = FALSE;
  1288. X    /*
  1289. X     * Recompute parent's mean and variance.
  1290. X     */
  1291. X    if (OVERLAPS_RECT(sq))
  1292. X        SSquareRecomputeStats(sq);
  1293. X}
  1294. X
  1295. XSSquareRecomputeStats(sq)
  1296. XSSquare *sq;
  1297. X{
  1298. X    Float maxp;
  1299. X    int in[4], closeflag;
  1300. X
  1301. X    in[0] = OVERLAPS_RECT(sq->child[0]);
  1302. X    in[1] = OVERLAPS_RECT(sq->child[1]);
  1303. X    in[2] = OVERLAPS_RECT(sq->child[2]);
  1304. X    in[3] = OVERLAPS_RECT(sq->child[3]);
  1305. X
  1306. X    if ((in[0] && (sq->child[0]->closed < SSCLOSED)) ||
  1307. X        (in[1] && (sq->child[1]->closed < SSCLOSED)) ||
  1308. X        (in[2] && (sq->child[2]->closed < SSCLOSED)) ||
  1309. X        (in[3] && (sq->child[3]->closed < SSCLOSED))) {
  1310. X        maxp = 0.;
  1311. X        if (in[0])
  1312. X            maxp = max(maxp, sq->child[0]->prio);
  1313. X        if (in[1])
  1314. X            maxp = max(maxp, sq->child[1]->prio);
  1315. X        if (in[2])
  1316. X            maxp = max(maxp, sq->child[2]->prio);
  1317. X        if (in[3])
  1318. X            maxp = max(maxp, sq->child[3]->prio);
  1319. X        sq->closed = NOT_CLOSED;
  1320. X        sq->prio = maxp;
  1321. X    } else if ((sq->child[0]->closed == CLOSED_SUPER) &&
  1322. X           (sq->child[1]->closed == CLOSED_SUPER) &&
  1323. X           (sq->child[2]->closed == CLOSED_SUPER) &&
  1324. X           (sq->child[3]->closed == CLOSED_SUPER)) {
  1325. X        sq->closed = CLOSED_SUPER;
  1326. X        sq->prio = 0;
  1327. X#if 0
  1328. X    } else if ((!in[0] || sq->child[0]->closed >= SSCLOSED) &&
  1329. X           (!in[1] || sq->child[1]->closed >= SSCLOSED) &&
  1330. X           (!in[2] || sq->child[2]->closed >= SSCLOSED) &&
  1331. X           (!in[3] || sq->child[3]->closed >= SSCLOSED)) {
  1332. X        sq->closed = SSCLOSED;
  1333. X        sq->prio = 0;
  1334. X#endif
  1335. X    } else {
  1336. X        sq->closed = SSCLOSED;
  1337. X        sq->prio = 0;
  1338. X    }
  1339. X    if (sq->parent)
  1340. X        SSquareRecomputeStats(sq->parent);
  1341. X}
  1342. X
  1343. XSSquare *
  1344. XSSquareInstall(xp, yp, xs, ys, d, parent)
  1345. Xint xp, yp, xs, ys, d;
  1346. XSSquare *parent;
  1347. X{
  1348. X    SSquare *new;
  1349. X
  1350. X    new = SSquareCreate(xp, yp, xs, ys, d, parent);
  1351. X    SSquareDraw(new);
  1352. X    return new;
  1353. X}
  1354. X
  1355. XSSquare *
  1356. XSSquareSelect(list)
  1357. XSSquare *list;
  1358. X{
  1359. X    int i;
  1360. X    SSquare *res, *which;
  1361. X
  1362. X    /*
  1363. X     * If mousebutton is pressed,
  1364. X     * find the square in which the mouse is located and
  1365. X     * return that if not a leaf (pixel-sized square).
  1366. X     */
  1367. X    if (GraphicsLeftMouseEvent()) {
  1368. X        SuperSampleMode = 0;
  1369. X        if ((res = SSquareFetchAtMouse(list)) != (SSquare *)NULL)
  1370. X            return res;
  1371. X    }
  1372. X    else if (GraphicsRightMouseEvent()) {
  1373. X        SuperSampleMode = 1;
  1374. X        if ((res = SSquareFetchAtMouse(list)) != (SSquare *)NULL) {
  1375. X            return res;
  1376. X        }
  1377. X    }
  1378. X    if (list->closed >= SSCLOSED) {
  1379. X        if (Rectmode) {
  1380. X            Rectmode = FALSE;
  1381. X            RecomputePriority(SSquares);
  1382. X            return SSquareSelect(list);
  1383. X        }
  1384. X        return (SSquare *)NULL;
  1385. X    }
  1386. X    /*
  1387. X     * Otherwise, find the square with the greatest
  1388. X     * 'priority'.
  1389. X     */
  1390. X    res = list;
  1391. X    while (res && !res->leaf) {
  1392. X        which = (SSquare *)NULL;
  1393. X        for (i = 0; i < 4; i++) {
  1394. X            if ((res->child[i]->closed < SSCLOSED) &&
  1395. X                OVERLAPS_RECT(res->child[i])) {
  1396. X                which = res->child[i];
  1397. X                break;
  1398. X            }
  1399. X        }
  1400. X        while (++i < 4) {
  1401. X            if ((res->child[i]->closed < SSCLOSED) &&
  1402. X                which->prio < res->child[i]->prio &&
  1403. X                OVERLAPS_RECT(res->child[i]))
  1404. X                which = res->child[i];
  1405. X        }
  1406. X        res = which;
  1407. X    }
  1408. X    return res;
  1409. X}
  1410. X
  1411. XSSquare *
  1412. XSSquareFetchAtMouse(list)
  1413. XSSquare *list;
  1414. X{
  1415. X    SSquare *res;
  1416. X    int x, y;
  1417. X
  1418. X    /*
  1419. X     * Get mouse position.
  1420. X     */
  1421. X    GraphicsGetMousePos(&x, &y);
  1422. X    res = list;
  1423. X    while (!res->leaf && (res->closed < SSCLOSED)) { 
  1424. X        /*
  1425. X         * Find in which child the mouse is located.
  1426. X         */
  1427. X        if (x < res->child[1]->xpos) {
  1428. X            if (y < res->child[2]->ypos)
  1429. X                res = res->child[0];
  1430. X            else
  1431. X                res = res->child[2];
  1432. X        } else if (y < res->child[3]->ypos)
  1433. X            res = res->child[1];
  1434. X        else
  1435. X            res = res->child[3];
  1436. X    }
  1437. X    if (res->closed >= SSCLOSED)
  1438. X        return (SSquare *)NULL;
  1439. X    return res;
  1440. X}
  1441. X
  1442. XSSquareDraw(sq)
  1443. XSSquare *sq;
  1444. X{
  1445. X    if (SQ_AREA(sq) >= MINAREA)
  1446. X        GraphicsDrawRectangle(sq->xpos, sq->ypos, sq->xsize, sq->ysize,
  1447. X            Image[sq->ypos][sq->xpos],
  1448. X            Image[sq->ypos][sq->xpos+sq->xsize],
  1449. X            Image[sq->ypos+sq->ysize][sq->xpos+sq->xsize],
  1450. X            Image[sq->ypos+sq->ysize][sq->xpos]);
  1451. X    else
  1452. X        DrawPixels(sq->xpos, sq->ypos, sq->xsize, sq->ysize);
  1453. X    if (!sq->leaf) {
  1454. X        SSquareDraw(sq->child[0]);
  1455. X        SSquareDraw(sq->child[1]);
  1456. X        SSquareDraw(sq->child[2]);
  1457. X        SSquareDraw(sq->child[3]);
  1458. X    }
  1459. X}
  1460. X
  1461. XDrawPixels(xp, yp, xs, ys)
  1462. Xint xp, yp, xs, ys;
  1463. X{
  1464. X    int x, y, xi, yi;
  1465. X
  1466. X    yi = yp;
  1467. X    for (y = 0; y <= ys; y++, yi++) {
  1468. X        xi = xp;
  1469. X        for (x = 0; x <= xs; x++, xi++) {
  1470. X            SSquareSample(xi, yi, SuperSampleMode);
  1471. X            GraphicsDrawPixel(xi, yi, Image[yi][xi]);
  1472. X        }
  1473. X    }
  1474. X}
  1475. END_OF_FILE
  1476. if test 17440 -ne `wc -c <'raypaint/render.c'`; then
  1477.     echo shar: \"'raypaint/render.c'\" unpacked with wrong size!
  1478. fi
  1479. # end of 'raypaint/render.c'
  1480. fi
  1481. echo shar: End of archive 16 \(of 19\).
  1482. cp /dev/null ark16isdone
  1483. MISSING=""
  1484. for I in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 ; do
  1485.     if test ! -f ark${I}isdone ; then
  1486.     MISSING="${MISSING} ${I}"
  1487.     fi
  1488. done
  1489. if test "${MISSING}" = "" ; then
  1490.     echo You have unpacked all 19 archives.
  1491.     rm -f ark[1-9]isdone ark[1-9][0-9]isdone
  1492. else
  1493.     echo You still need to unpack the following archives:
  1494.     echo "        " ${MISSING}
  1495. fi
  1496. ##  End of shell archive.
  1497. exit 0
  1498.  
  1499. exit 0 # Just in case...
  1500.