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.3 by greg, Tue Oct 23 05:10:42 2012 UTC vs.
Revision 2.9 by greg, Fri Dec 14 22:41:44 2012 UTC

# Line 51 | Line 51 | order_triangle(MIGRATION *miga[3])
51                  insert_vert(vert, miga[i]->rbfv[1]);
52          }
53                                                  /* should be just 3 vertices */
54 <        if ((vert[3] == NULL) | (vert[4] != NULL))
54 >        if ((vert[2] == NULL) | (vert[3] != NULL))
55                  return(0);
56                                                  /* identify edge 0 */
57          for (i = 3; i--; )
# Line 85 | Line 85 | order_triangle(MIGRATION *miga[3])
85          return(1);
86   }
87  
88 < /* Determine if we are close enough to the given edge */
88 > /* Determine if we are close enough to an 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));
92 >        double  cos_a, cos_b, cos_c, cos_aplusb;
93 >                                        /* use triangle inequality */
94 >        cos_a = DOT(ej->rbfv[0]->invec, ivec);
95 >        if (cos_a <= 0)
96 >                return(0);
97  
98 <        return(cos_aplusb - cos_c < .01);
98 >        cos_b = DOT(ej->rbfv[1]->invec, ivec);
99 >        if (cos_b <= 0)
100 >                return(0);
101 >
102 >        cos_aplusb = cos_a*cos_b - sqrt((1.-cos_a*cos_a)*(1.-cos_b*cos_b));
103 >        if (cos_aplusb <= 0)
104 >                return(0);
105 >
106 >        cos_c = DOT(ej->rbfv[0]->invec, ej->rbfv[1]->invec);
107 >
108 >        return(cos_c - cos_aplusb < .001);
109   }
110  
111   /* Determine if we are inside the given triangle */
# Line 116 | Line 126 | in_tri(const RBFNODE *v1, const RBFNODE *v2, const RBF
126          return(sgn2 == sgn3);
127   }
128  
129 < /* Compute intersection with the given position over remaining mesh */
129 > /* Test (and set) bitmap for edge */
130   static int
131 < in_mesh(MIGRATION *miga[3], unsigned char *emap, int nedges,
122 <                        const FVECT ivec, MIGRATION *mig)
131 > check_edge(unsigned char *emap, int nedges, const MIGRATION *mig, int mark)
132   {
133 <        MIGRATION       *ej1, *ej2;
134 <        RBFNODE         *tv;
126 <        int             ejndx;
127 <                                                /* check visitation record */
133 >        int     ejndx, bit2check;
134 >
135          if (mig->rbfv[0]->ord > mig->rbfv[1]->ord)
136                  ejndx = mig->rbfv[1]->ord + (nedges-1)*mig->rbfv[0]->ord;
137          else
138                  ejndx = mig->rbfv[0]->ord + (nedges-1)*mig->rbfv[1]->ord;
139 <        if (emap[ejndx>>3] & 1<<(ejndx&07))     /* tested already? */
139 >
140 >        bit2check = 1<<(ejndx&07);
141 >
142 >        if (emap[ejndx>>3] & bit2check)
143                  return(0);
144 <        emap[ejndx>>3] |= 1<<(ejndx&07);        /* else mark & test it */
144 >        if (mark)
145 >                emap[ejndx>>3] |= bit2check;
146 >        return(1);
147 > }
148 >
149 > /* Compute intersection with the given position over remaining mesh */
150 > static int
151 > in_mesh(MIGRATION *miga[3], unsigned char *emap, int nedges,
152 >                        const FVECT ivec, MIGRATION *mig)
153 > {
154 >        RBFNODE         *tv[2];
155 >        MIGRATION       *sej[2], *dej[2];
156 >        int             i;
157 >                                                /* check visitation record */
158 >        if (!check_edge(emap, nedges, mig, 1))
159 >                return(0);
160          if (on_edge(mig, ivec)) {
161                  miga[0] = mig;                  /* close enough to edge */
162                  return(1);
163          }
164 <                                                /* do triangles either side */
165 <        for (ej1 = mig->rbfv[0]->ejl; ej1 != NULL;
166 <                                ej1 = nextedge(mig->rbfv[0],ej1)) {
167 <                if (ej1 == mig)
164 >        if (!get_triangles(tv, mig))            /* do triangles either side? */
165 >                return(0);
166 >        for (i = 2; i--; ) {                    /* identify edges to check */
167 >                MIGRATION       *ej;
168 >                sej[i] = dej[i] = NULL;
169 >                if (tv[i] == NULL)
170                          continue;
171 <                tv = opp_rbf(mig->rbfv[0],ej1);
172 <                for (ej2 = tv->ejl; ej2 != NULL; ej2 = nextedge(tv,ej2))
173 <                        if (opp_rbf(tv,ej2) == mig->rbfv[1]) {
174 <                                if (in_mesh(miga, emap, nedges, ivec, ej1))
175 <                                        return(1);
176 <                                if (in_mesh(miga, emap, nedges, ivec, ej2))
177 <                                        return(1);
178 <                                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 <                                }
171 >                for (ej = tv[i]->ejl; ej != NULL; ej = nextedge(tv[i],ej)) {
172 >                        RBFNODE *rbfop = opp_rbf(tv[i],ej);
173 >                        if (rbfop == mig->rbfv[0]) {
174 >                                if (check_edge(emap, nedges, ej, 0))
175 >                                        sej[i] = ej;
176 >                        } else if (rbfop == mig->rbfv[1]) {
177 >                                if (check_edge(emap, nedges, ej, 0))
178 >                                        dej[i] = ej;
179                          }
180 +                }
181          }
182 <        return(0);
182 >        for (i = 2; i--; ) {                    /* check triangles just once */
183 >                if (sej[i] != NULL && in_mesh(miga, emap, nedges, ivec, sej[i]))
184 >                        return(1);
185 >                if (dej[i] != NULL && in_mesh(miga, emap, nedges, ivec, dej[i]))
186 >                        return(1);
187 >                if ((sej[i] == NULL) | (dej[i] == NULL))
188 >                        continue;
189 >                if (in_tri(mig->rbfv[0], mig->rbfv[1], tv[i], ivec)) {
190 >                        miga[0] = mig;
191 >                        miga[1] = sej[i];
192 >                        miga[2] = dej[i];
193 >                        return(1);
194 >                }
195 >        }
196 >        return(0);                              /* not near this edge */
197   }
198  
199   /* Find edge(s) for interpolating the given vector, applying symmetry */
# Line 175 | Line 211 | get_interp(MIGRATION *miga[3], FVECT invec)
211                                                          input_orient*invec[2]) {
212                                  for (miga[0] = rbf->ejl; miga[0] != NULL;
213                                                  miga[0] = nextedge(rbf,miga[0]))
214 <                                        if (opp_rbf(rbf,miga[0]) == rbf->next)
214 >                                        if (opp_rbf(rbf,miga[0]) == rbf->next) {
215 >                                                double  nf = 1.-rbf->invec[2]*rbf->invec[2];
216 >                                                if (nf > FTINY) {
217 >                                                        nf = sqrt((1.-invec[2]*invec[2])/nf);
218 >                                                        invec[0] = nf*rbf->invec[0];
219 >                                                        invec[1] = nf*rbf->invec[1];
220 >                                                }
221                                                  return(0);
222 +                                        }
223                                  break;
224                          }
225                  }
# Line 197 | Line 240 | get_interp(MIGRATION *miga[3], FVECT invec)
240                          exit(1);
241                  }
242                                                  /* identify intersection  */
243 <                if (!in_mesh(miga, emap, nedges, invec, mig_list))
243 >                if (!in_mesh(miga, emap, nedges, invec, mig_list)) {
244 > #ifdef DEBUG
245 >                        fprintf(stderr,
246 >                        "Incident angle (%.1f,%.1f) deg. outside mesh\n",
247 >                                        get_theta180(invec), get_phi360(invec));
248 > #endif
249                          sym = -1;               /* outside mesh */
250 <                else if (miga[1] != NULL &&
250 >                } else if (miga[1] != NULL &&
251                                  (miga[2] == NULL || !order_triangle(miga))) {
252   #ifdef DEBUG
253                          fputs("Munged triangle in get_interp()\n", stderr);
# Line 220 | Line 268 | e_advect_rbf(const MIGRATION *mig, const FVECT invec)
268          double          t, full_dist;
269                                                  /* get relative position */
270          t = acos(DOT(invec, mig->rbfv[0]->invec));
271 <        if (t < M_PI/GRIDRES) {                 /* near first DSF */
271 >        if (t < M_PI/grid_res) {                /* near first DSF */
272                  n = sizeof(RBFNODE) + sizeof(RBFVAL)*(mig->rbfv[0]->nrbf-1);
273                  rbf = (RBFNODE *)malloc(n);
274                  if (rbf == NULL)
275                          goto memerr;
276                  memcpy(rbf, mig->rbfv[0], n);   /* just duplicate */
277 +                rbf->next = NULL; rbf->ejl = NULL;
278                  return(rbf);
279          }
280          full_dist = acos(DOT(mig->rbfv[0]->invec, mig->rbfv[1]->invec));
281 <        if (t > full_dist-M_PI/GRIDRES) {       /* near second DSF */
281 >        if (t > full_dist-M_PI/grid_res) {      /* near second DSF */
282                  n = sizeof(RBFNODE) + sizeof(RBFVAL)*(mig->rbfv[1]->nrbf-1);
283                  rbf = (RBFNODE *)malloc(n);
284                  if (rbf == NULL)
285                          goto memerr;
286                  memcpy(rbf, mig->rbfv[1], n);   /* just duplicate */
287 +                rbf->next = NULL; rbf->ejl = NULL;
288                  return(rbf);
289          }
290          t /= full_dist;
# Line 305 | Line 355 | advect_rbf(const FVECT invec)
355                  return(NULL);
356          if (miga[1] == NULL) {                  /* advect along edge? */
357                  rbf = e_advect_rbf(miga[0], sivec);
358 <                rev_rbf_symmetry(rbf, sym);
358 >                if (single_plane_incident)
359 >                        rotate_rbf(rbf, invec);
360 >                else
361 >                        rev_rbf_symmetry(rbf, sym);
362                  return(rbf);
363          }
364   #ifdef DEBUG
# Line 332 | Line 385 | advect_rbf(const FVECT invec)
385              for (j = 0; j < mtx_ncols(miga[0]); j++)
386                  for (k = (mtx_coef(miga[0],i,j) > FTINY) *
387                                          mtx_ncols(miga[2]); k--; )
388 <                        n += (mtx_coef(miga[2],i,k) > FTINY &&
388 >                        n += (mtx_coef(miga[2],i,k) > FTINY ||
389                                  mtx_coef(miga[1],j,k) > FTINY);
390   #ifdef DEBUG
391          fprintf(stderr, "Input RBFs have %d, %d, %d nodes -> output has %d\n",
# Line 373 | Line 426 | advect_rbf(const FVECT invec)
426                      float               mc = mtx_coef(miga[2],i,k);
427                      const RBFVAL        *rbf2k;
428                      double              rad2k;
376 –                    FVECT               vout;
429                      int                 pos[2];
430 <                    if ((mb <= FTINY) | (mc <= FTINY))
430 >                    if ((mb <= FTINY) & (mc <= FTINY))
431                          continue;
432                      rbf2k = &miga[2]->rbfv[1]->rbfa[k];
433                      rbf->rbfa[n].peak = w0i * ma * (mb*mbfact + mc*mcfact);
434                      rad2k = R2ANG(rbf2k->crad);
435                      rbf->rbfa[n].crad = ANG2R(sqrt(srad2 + t*rad2k*rad2k));
436                      ovec_from_pos(v2, rbf2k->gx, rbf2k->gy);
437 <                    geodesic(vout, v1, v2, t, GEOD_REL);
438 <                    pos_from_vec(pos, vout);
437 >                    geodesic(v2, v1, v2, t, GEOD_REL);
438 >                    pos_from_vec(pos, v2);
439                      rbf->rbfa[n].gx = pos[0];
440                      rbf->rbfa[n].gy = pos[1];
441                      ++n;

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines