ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/cv/bsdfinterp.c
(Generate patch)

Comparing ray/src/cv/bsdfinterp.c (file contents):
Revision 2.2 by greg, Sat Oct 20 07:02:00 2012 UTC vs.
Revision 2.3 by greg, Tue Oct 23 05:10:42 2012 UTC

# Line 13 | Line 13 | static const char RCSid[] = "$Id$";
13   #include <string.h>
14   #include <math.h>
15   #include "bsdfrep.h"
16                                /* migration edges drawn in raster fashion */
17 MIGRATION               *mig_grid[GRIDRES][GRIDRES];
16  
19 #ifdef DEBUG
20 #include "random.h"
21 #include "bmpfile.h"
22 /* Hash pointer to byte value (must return 0 for NULL) */
23 static int
24 byte_hash(const void *p)
25 {
26        size_t  h = (size_t)p;
27        h ^= (size_t)p >> 8;
28        h ^= (size_t)p >> 16;
29        h ^= (size_t)p >> 24;
30        return(h & 0xff);
31 }
32 /* Write out BMP image showing edges */
33 static void
34 write_edge_image(const char *fname)
35 {
36        BMPHeader       *hdr = BMPmappedHeader(GRIDRES, GRIDRES, 0, 256);
37        BMPWriter       *wtr;
38        int             i, j;
39
40        fprintf(stderr, "Writing incident mesh drawing to '%s'\n", fname);
41        hdr->compr = BI_RLE8;
42        for (i = 256; --i; ) {                  /* assign random color map */
43                hdr->palette[i].r = random() & 0xff;
44                hdr->palette[i].g = random() & 0xff;
45                hdr->palette[i].b = random() & 0xff;
46                                                /* reject dark colors */
47                i += (hdr->palette[i].r + hdr->palette[i].g +
48                                                hdr->palette[i].b < 128);
49        }
50        hdr->palette[0].r = hdr->palette[0].g = hdr->palette[0].b = 0;
51                                                /* open output */
52        wtr = BMPopenOutputFile(fname, hdr);
53        if (wtr == NULL) {
54                free(hdr);
55                return;
56        }
57        for (i = 0; i < GRIDRES; i++) {         /* write scanlines */
58                for (j = 0; j < GRIDRES; j++)
59                        wtr->scanline[j] = byte_hash(mig_grid[i][j]);
60                if (BMPwriteScanline(wtr) != BIR_OK)
61                        break;
62        }
63        BMPcloseOutput(wtr);                    /* close & clean up */
64 }
65 #endif
66
67 /* Draw edge list into mig_grid array */
68 void
69 draw_edges(void)
70 {
71        int             nnull = 0, ntot = 0;
72        MIGRATION       *ej;
73        int             p0[2], p1[2];
74
75        memset(mig_grid, 0, sizeof(mig_grid));
76        for (ej = mig_list; ej != NULL; ej = ej->next) {
77                ++ntot;
78                pos_from_vec(p0, ej->rbfv[0]->invec);
79                pos_from_vec(p1, ej->rbfv[1]->invec);
80                if ((p0[0] == p1[0]) & (p0[1] == p1[1])) {
81                        ++nnull;
82                        mig_grid[p0[0]][p0[1]] = ej;
83                        continue;
84                }
85                if (abs(p1[0]-p0[0]) > abs(p1[1]-p0[1])) {
86                        const int       xstep = 2*(p1[0] > p0[0]) - 1;
87                        const double    ystep = (double)((p1[1]-p0[1])*xstep) /
88                                                        (double)(p1[0]-p0[0]);
89                        int             x;
90                        double          y;
91                        for (x = p0[0], y = p0[1]+.5; x != p1[0];
92                                                x += xstep, y += ystep)
93                                mig_grid[x][(int)y] = ej;
94                        mig_grid[x][(int)y] = ej;
95                } else {
96                        const int       ystep = 2*(p1[1] > p0[1]) - 1;
97                        const double    xstep = (double)((p1[0]-p0[0])*ystep) /
98                                                        (double)(p1[1]-p0[1]);
99                        int             y;
100                        double          x;
101                        for (y = p0[1], x = p0[0]+.5; y != p1[1];
102                                                y += ystep, x += xstep)
103                                mig_grid[(int)x][y] = ej;
104                        mig_grid[(int)x][y] = ej;
105                }
106        }
107        if (nnull)
108                fprintf(stderr, "Warning: %d of %d edges are null\n",
109                                nnull, ntot);
110 #ifdef DEBUG
111        write_edge_image("bsdf_edges.bmp");
112 #endif
113 }
114
115 /* Identify enclosing triangle for this position (flood fill raster check) */
116 static int
117 identify_tri(MIGRATION *miga[3], unsigned char vmap[GRIDRES][(GRIDRES+7)/8],
118                        int px, int py)
119 {
120        const int       btest = 1<<(py&07);
121
122        if (vmap[px][py>>3] & btest)            /* already visited here? */
123                return(1);
124                                                /* else mark it */
125        vmap[px][py>>3] |= btest;
126
127        if (mig_grid[px][py] != NULL) {         /* are we on an edge? */
128                int     i;
129                for (i = 0; i < 3; i++) {
130                        if (miga[i] == mig_grid[px][py])
131                                return(1);
132                        if (miga[i] != NULL)
133                                continue;
134                        miga[i] = mig_grid[px][py];
135                        return(1);
136                }
137                return(0);                      /* outside triangle! */
138        }
139                                                /* check neighbors (flood) */
140        if (px > 0 && !identify_tri(miga, vmap, px-1, py))
141                return(0);
142        if (px < GRIDRES-1 && !identify_tri(miga, vmap, px+1, py))
143                return(0);
144        if (py > 0 && !identify_tri(miga, vmap, px, py-1))
145                return(0);
146        if (py < GRIDRES-1 && !identify_tri(miga, vmap, px, py+1))
147                return(0);
148        return(1);                              /* this neighborhood done */
149 }
150
17   /* Insert vertex in ordered list */
18   static void
19   insert_vert(RBFNODE **vlist, RBFNODE *v)
# Line 219 | Line 85 | order_triangle(MIGRATION *miga[3])
85          return(1);
86   }
87  
88 + /* Determine if we are close enough to the given edge */
89 + static int
90 + on_edge(const MIGRATION *ej, const FVECT ivec)
91 + {
92 +        double  cos_a = DOT(ej->rbfv[0]->invec, ivec);
93 +        double  cos_b = DOT(ej->rbfv[1]->invec, ivec);
94 +        double  cos_c = DOT(ej->rbfv[0]->invec, ej->rbfv[1]->invec);
95 +        double  cos_aplusb = cos_a*cos_b -
96 +                                sqrt((1.-cos_a*cos_a)*(1.-cos_b*cos_b));
97 +
98 +        return(cos_aplusb - cos_c < .01);
99 + }
100 +
101 + /* Determine if we are inside the given triangle */
102 + static int
103 + in_tri(const RBFNODE *v1, const RBFNODE *v2, const RBFNODE *v3, const FVECT p)
104 + {
105 +        FVECT   vc;
106 +        int     sgn1, sgn2, sgn3;
107 +                                        /* signed volume test */
108 +        VCROSS(vc, v1->invec, v2->invec);
109 +        sgn1 = (DOT(p, vc) > 0);
110 +        VCROSS(vc, v2->invec, v3->invec);
111 +        sgn2 = (DOT(p, vc) > 0);
112 +        if (sgn1 != sgn2)
113 +                return(0);
114 +        VCROSS(vc, v3->invec, v1->invec);
115 +        sgn3 = (DOT(p, vc) > 0);
116 +        return(sgn2 == sgn3);
117 + }
118 +
119 + /* Compute intersection with the given position over remaining mesh */
120 + static int
121 + in_mesh(MIGRATION *miga[3], unsigned char *emap, int nedges,
122 +                        const FVECT ivec, MIGRATION *mig)
123 + {
124 +        MIGRATION       *ej1, *ej2;
125 +        RBFNODE         *tv;
126 +        int             ejndx;
127 +                                                /* check visitation record */
128 +        if (mig->rbfv[0]->ord > mig->rbfv[1]->ord)
129 +                ejndx = mig->rbfv[1]->ord + (nedges-1)*mig->rbfv[0]->ord;
130 +        else
131 +                ejndx = mig->rbfv[0]->ord + (nedges-1)*mig->rbfv[1]->ord;
132 +        if (emap[ejndx>>3] & 1<<(ejndx&07))     /* tested already? */
133 +                return(0);
134 +        emap[ejndx>>3] |= 1<<(ejndx&07);        /* else mark & test it */
135 +        if (on_edge(mig, ivec)) {
136 +                miga[0] = mig;                  /* close enough to edge */
137 +                return(1);
138 +        }
139 +                                                /* do triangles either side */
140 +        for (ej1 = mig->rbfv[0]->ejl; ej1 != NULL;
141 +                                ej1 = nextedge(mig->rbfv[0],ej1)) {
142 +                if (ej1 == mig)
143 +                        continue;
144 +                tv = opp_rbf(mig->rbfv[0],ej1);
145 +                for (ej2 = tv->ejl; ej2 != NULL; ej2 = nextedge(tv,ej2))
146 +                        if (opp_rbf(tv,ej2) == mig->rbfv[1]) {
147 +                                if (in_mesh(miga, emap, nedges, ivec, ej1))
148 +                                        return(1);
149 +                                if (in_mesh(miga, emap, nedges, ivec, ej2))
150 +                                        return(1);
151 +                                if (in_tri(mig->rbfv[0], mig->rbfv[1],
152 +                                                tv, ivec)) {
153 +                                        miga[0] = mig;
154 +                                        miga[1] = ej1;
155 +                                        miga[2] = ej2;
156 +                                        return(1);
157 +                                }
158 +                        }
159 +        }
160 +        return(0);
161 + }
162 +
163   /* Find edge(s) for interpolating the given vector, applying symmetry */
164   int
165   get_interp(MIGRATION *miga[3], FVECT invec)
# Line 242 | Line 183 | get_interp(MIGRATION *miga[3], FVECT invec)
183                  return(-1);                     /* outside range! */
184          }
185          {                                       /* else use triangle mesh */
186 <                const int       sym = use_symmetry(invec);
187 <                unsigned char   floodmap[GRIDRES][(GRIDRES+7)/8];
188 <                int             pstart[2];
189 <                RBFNODE         *vother;
190 <                MIGRATION       *ej;
191 <                int             i;
192 <
193 <                pos_from_vec(pstart, invec);
194 <                memset(floodmap, 0, sizeof(floodmap));
195 <                                                /* call flooding function */
196 <                if (!identify_tri(miga, floodmap, pstart[0], pstart[1]))
197 <                        return(-1);             /* outside mesh */
198 <                if ((miga[0] == NULL) | (miga[2] == NULL))
199 <                        return(-1);             /* should never happen */
200 <                if (miga[1] == NULL)
201 <                        return(sym);            /* on edge */
202 <                                                /* verify triangle */
203 <                if (!order_triangle(miga)) {
186 >                int             sym = use_symmetry(invec);
187 >                int             nedges = 0;
188 >                MIGRATION       *mep;
189 >                unsigned char   *emap;
190 >                                                /* clear visitation map */
191 >                for (mep = mig_list; mep != NULL; mep = mep->next)
192 >                        ++nedges;
193 >                emap = (unsigned char *)calloc((nedges*(nedges-1) + 7)>>3, 1);
194 >                if (emap == NULL) {
195 >                        fprintf(stderr, "%s: Out of memory in get_interp()\n",
196 >                                        progname);
197 >                        exit(1);
198 >                }
199 >                                                /* identify intersection  */
200 >                if (!in_mesh(miga, emap, nedges, invec, mig_list))
201 >                        sym = -1;               /* outside mesh */
202 >                else if (miga[1] != NULL &&
203 >                                (miga[2] == NULL || !order_triangle(miga))) {
204   #ifdef DEBUG
205                          fputs("Munged triangle in get_interp()\n", stderr);
206   #endif
207 <                        vother = NULL;          /* find triangle from edge */
267 <                        for (i = 3; i--; ) {
268 <                            RBFNODE     *tpair[2];
269 <                            if (get_triangles(tpair, miga[i]) &&
270 <                                        (vother = tpair[ is_rev_tri(
271 <                                                        miga[i]->rbfv[0]->invec,
272 <                                                        miga[i]->rbfv[1]->invec,
273 <                                                        invec) ]) != NULL)
274 <                                        break;
275 <                        }
276 <                        if (vother == NULL) {   /* couldn't find 3rd vertex */
277 < #ifdef DEBUG
278 <                                fputs("No triangle in get_interp()\n", stderr);
279 < #endif
280 <                                return(-1);
281 <                        }
282 <                                                /* reassign other two edges */
283 <                        for (ej = vother->ejl; ej != NULL;
284 <                                                ej = nextedge(vother,ej)) {
285 <                                RBFNODE *vorig = opp_rbf(vother,ej);
286 <                                if (vorig == miga[i]->rbfv[0])
287 <                                        miga[(i+1)%3] = ej;
288 <                                else if (vorig == miga[i]->rbfv[1])
289 <                                        miga[(i+2)%3] = ej;
290 <                        }
291 <                        if (!order_triangle(miga)) {
292 < #ifdef DEBUG
293 <                                fputs("Bad triangle in get_interp()\n", stderr);
294 < #endif
295 <                                return(-1);
296 <                        }
207 >                        sym = -1;
208                  }
209 +                free(emap);
210                  return(sym);                    /* return in standard order */
211          }
212   }

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines