brought back x-mirror editing, though it's currently buggy. also made tesselation...
[blender.git] / source / blender / editors / mesh / editbmesh_bvh.c
1  /* $Id:
2  *
3  * ***** BEGIN GPL LICENSE BLOCK *****
4  *
5  * This program is free software; you can redistribute it and/or
6  * modify it under the terms of the GNU General Public License
7  * as published by the Free Software Foundation; either version 2
8  * of the License, or (at your option) any later version.
9  *
10  * This program is distributed in the hope that it will be useful,
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13  * GNU General Public License for more details.
14  *
15  * You should have received a copy of the GNU General Public License
16  * along with this program; if not, write to the Free Software Foundation,
17  * Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
18  *
19  * The Original Code is Copyright (C) 2004 by Blender Foundation.
20  * All rights reserved.
21  *
22  * The Original Code is: all of this file.
23  *
24  * Contributor(s): Joseph Eagar
25  *
26  * ***** END GPL LICENSE BLOCK *****
27  */
28 #include <stdlib.h>
29 #include <stdarg.h>
30 #include <string.h>
31 #include <math.h>
32 #include <float.h>
33
34 #include "MEM_guardedalloc.h"
35 #include "PIL_time.h"
36
37 #include "BLO_sys_types.h" // for intptr_t support
38
39 #include "DNA_mesh_types.h"
40 #include "DNA_material_types.h"
41 #include "DNA_meshdata_types.h"
42 #include "DNA_modifier_types.h"
43 #include "DNA_object_types.h"
44 #include "DNA_scene_types.h"
45 #include "DNA_screen_types.h"
46 #include "DNA_view3d_types.h"
47 #include "DNA_key_types.h"
48 #include "DNA_windowmanager_types.h"
49
50 #include "RNA_types.h"
51 #include "RNA_define.h"
52 #include "RNA_access.h"
53
54 #include "BLI_blenlib.h"
55 #include "BLI_math.h"
56 #include "BLI_editVert.h"
57 #include "BLI_rand.h"
58 #include "BLI_ghash.h"
59 #include "BLI_linklist.h"
60 #include "BLI_heap.h"
61 #include "BLI_array.h"
62 #include "BLI_kdopbvh.h"
63
64 #include "BKE_context.h"
65 #include "BKE_customdata.h"
66 #include "BKE_depsgraph.h"
67 #include "BKE_global.h"
68 #include "BKE_library.h"
69 #include "BKE_mesh.h"
70 #include "BKE_object.h"
71 #include "BKE_utildefines.h"
72 #include "BKE_bmesh.h"
73 #include "BKE_report.h"
74 #include "BKE_tessmesh.h"
75
76 #include "BIF_gl.h"
77 #include "BIF_glutil.h"
78
79 #include "WM_api.h"
80 #include "WM_types.h"
81
82 #include "ED_mesh.h"
83 #include "ED_view3d.h"
84 #include "ED_util.h"
85 #include "ED_screen.h"
86 #include "ED_transform.h"
87
88 #include "UI_interface.h"
89
90 #include "mesh_intern.h"
91 #include "bmesh.h"
92
93 #define IN_EDITMESHBVH
94 #include "editbmesh_bvh.h"
95
96 typedef struct BMBVHTree {
97         BMEditMesh *em;
98         BMesh *bm;
99         BVHTree *tree;
100         float epsilon;
101         float maxdist; //for nearest point search
102
103         /*stuff for topological vert search*/
104         BMVert *v, *curv;
105         GHash *gh;
106         float curw, curd;
107         float co[3];
108         int curtag;
109 } BMBVHTree;
110
111 BMBVHTree *BMBVH_NewBVH(BMEditMesh *em)
112 {
113         BMBVHTree *tree = MEM_callocN(sizeof(*tree), "BMBVHTree");
114         float cos[3][3];
115         int i;
116
117         BMEdit_RecalcTesselation(em);
118
119         tree->em = em;
120         tree->bm = em->bm;
121         tree->epsilon = FLT_EPSILON*2.0f;
122
123         tree->tree = BLI_bvhtree_new(em->tottri, tree->epsilon, 8, 8);
124
125         for (i=0; i<em->tottri; i++) {
126                 VECCOPY(cos[0], em->looptris[i][0]->v->co);
127                 VECCOPY(cos[1], em->looptris[i][1]->v->co);
128                 VECCOPY(cos[2], em->looptris[i][2]->v->co);
129
130                 BLI_bvhtree_insert(tree->tree, i, (float*)cos, 3);
131         }
132         
133         BLI_bvhtree_balance(tree->tree);
134         
135         return tree;
136 }
137
138 void BMBVH_FreeBVH(BMBVHTree *tree)
139 {
140         BLI_bvhtree_free(tree->tree);
141         MEM_freeN(tree);
142 }
143
144 /*taken from bvhutils.c*/
145 static float ray_tri_intersection(const BVHTreeRay *ray, const float m_dist, float *v0, 
146                                   float *v1, float *v2, float *uv, float e)
147 {
148         float dist;
149 #if 0
150         float vv1[3], vv2[3], vv3[3], cent[3];
151
152         /*expand triangle by an epsilon.  this is probably a really stupid
153           way of doing it, but I'm too tired to do better work.*/
154         VECCOPY(vv1, v0);
155         VECCOPY(vv2, v1);
156         VECCOPY(vv3, v2);
157
158         add_v3_v3v3(cent, vv1, vv2);
159         add_v3_v3v3(cent, cent, vv3);
160         mul_v3_fl(cent, 1.0f/3.0f);
161
162         sub_v3_v3v3(vv1, vv1, cent);
163         sub_v3_v3v3(vv2, vv2, cent);
164         sub_v3_v3v3(vv3, vv3, cent);
165
166         mul_v3_fl(vv1, 1.0f + e);
167         mul_v3_fl(vv2, 1.0f + e);
168         mul_v3_fl(vv3, 1.0f + e);
169
170         add_v3_v3v3(vv1, vv1, cent);
171         add_v3_v3v3(vv2, vv2, cent);
172         add_v3_v3v3(vv3, vv3, cent);
173
174         if(isect_ray_tri_v3((float*)ray->origin, (float*)ray->direction, vv1, vv2, vv3, &dist, uv))
175                 return dist;
176 #else
177         if(isect_ray_tri_v3((float*)ray->origin, (float*)ray->direction, v0, v1, v2, &dist, uv))
178                 return dist;
179 #endif
180
181         return FLT_MAX;
182 }
183
184 static void raycallback(void *userdata, int index, const BVHTreeRay *ray, BVHTreeRayHit *hit)
185 {
186         BMBVHTree *tree = userdata;
187         BMLoop **ls = tree->em->looptris[index];
188         float dist, uv[2], co1[3], co2[3], co3[3];
189
190         dist = ray_tri_intersection(ray, hit->dist, ls[0]->v->co, ls[1]->v->co,
191                                     ls[2]->v->co, uv, tree->epsilon);
192         if (dist < hit->dist) {
193                 hit->dist = dist;
194                 hit->index = index;
195                 
196                 VECCOPY(hit->no, ls[0]->v->no);
197
198                 VECCOPY(co1, ls[0]->v->co);
199                 VECCOPY(co2, ls[1]->v->co);
200                 VECCOPY(co3, ls[2]->v->co);
201
202                 mul_v3_fl(co1, uv[0]);
203                 mul_v3_fl(co2, uv[1]);
204                 mul_v3_fl(co3, 1.0f-uv[0]-uv[1]);
205
206                 add_v3_v3v3(hit->co, co1, co2);
207                 add_v3_v3v3(hit->co, hit->co, co3);
208         }
209 }
210
211 BMFace *BMBVH_RayCast(BMBVHTree *tree, float *co, float *dir, float *hitout)
212 {
213         BVHTreeRayHit hit;
214
215         hit.dist = FLT_MAX;
216         hit.index = -1;
217
218         BLI_bvhtree_ray_cast(tree->tree, co, dir, FLT_MAX, &hit, raycallback, tree);
219         if (hit.dist != FLT_MAX && hit.index != -1) {
220                 if (hitout) {
221                         VECCOPY(hitout, hit.co);
222                 }
223
224                 return tree->em->looptris[hit.index][0]->f;
225         }
226
227         return NULL;
228 }
229
230
231 static void vertsearchcallback(void *userdata, int index, const float *co, BVHTreeNearest *hit)
232 {
233         BMBVHTree *tree = userdata;
234         BMLoop **ls = tree->em->looptris[index];
235         float dist, maxdist, v[3];
236         int i;
237
238         maxdist = tree->maxdist;
239
240         for (i=0; i<3; i++) {
241                 sub_v3_v3v3(v, hit->co, ls[i]->v->co);
242
243                 dist = sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]);
244                 if (dist < hit->dist && dist < maxdist) {
245                         VECCOPY(hit->co, ls[i]->v->co);
246                         VECCOPY(hit->no, ls[i]->v->no);
247                         hit->dist = dist;
248                 }
249         }
250 }
251
252 BMVert *BMBVH_FindClosestVert(BMBVHTree *tree, float *co, float maxdist)
253 {
254         BVHTreeNearest hit;
255
256         VECCOPY(hit.co, co);
257         hit.dist = maxdist*5;
258         hit.index = -1;
259
260         tree->maxdist = maxdist;
261
262         BLI_bvhtree_find_nearest(tree->tree, co, &hit, vertsearchcallback, tree);
263         if (hit.dist != FLT_MAX && hit.index != -1) {
264                 BMLoop **ls = tree->em->looptris[hit.index];
265                 float dist, curdist = tree->maxdist, v[3];
266                 int cur=0, i;
267
268                 maxdist = tree->maxdist;
269
270                 for (i=0; i<3; i++) {
271                         sub_v3_v3v3(v, hit.co, ls[i]->v->co);
272
273                         dist = sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]);
274                         if (dist < curdist) {
275                                 cur = i;
276                                 curdist = dist;
277                         }
278                 }
279
280                 return ls[i]->v;
281         }
282
283         return NULL;
284 }
285
286 typedef struct walklist {
287         BMVert *v;
288         int valence;
289         int depth;
290         float w, r;
291         int totwalked;
292
293         /*state data*/
294         BMVert *lastv;
295         BMLoop *curl, *firstl;
296         BMEdge *cure;
297 } walklist;
298
299
300 static short winding(float *v1, float *v2, float *v3)
301 /* is v3 to the right of v1-v2 ? With exception: v3==v1 || v3==v2 */
302 {
303         double inp;
304
305         //inp= (v2[cox]-v1[cox])*(v1[coy]-v3[coy]) +(v1[coy]-v2[coy])*(v1[cox]-v3[cox]);
306         inp= (v2[0]-v1[0])*(v1[1]-v3[1]) +(v1[1]-v2[1])*(v1[0]-v3[0]);
307
308         if(inp<0.0) return 0;
309         else if(inp==0) {
310                 if(v1[0]==v3[0] && v1[1]==v3[1]) return 0;
311                 if(v2[0]==v3[0] && v2[1]==v3[1]) return 0;
312         }
313         return 1;
314 }
315
316 static float topo_compare(BMesh *bm, BMVert *v1, BMVert *v2, int tag)
317 {
318         BMIter iter1, iter2;
319         BMEdge *e1, *e2, *cure1 = NULL, *cure2 = NULL;
320         BMLoop *l1, *l2;
321         BMVert *lastv1, *lastv2;
322         GHash *gh;
323         walklist *stack1=NULL, *stack2=NULL;
324         BLI_array_declare(stack1);
325         BLI_array_declare(stack2);
326         float vec1[3], vec2[3], minangle=FLT_MAX, w;
327         int lvl=1;
328         static int maxlevel = 3;
329
330         /*ok.  see how similar v is to v2, based on topological similaritys in the local
331           topological neighborhood*/
332
333         /*step 1: find two edges, one that contains v and one that contains v2, with the
334           smallest angle between the two edges*/
335
336         BM_ITER(e1, &iter1, bm, BM_EDGES_OF_VERT, v1) {
337                 BM_ITER(e2, &iter2, bm, BM_EDGES_OF_VERT, v2) {
338                         float angle;
339                         
340                         if (e1->v1 == e2->v1 || e1->v2 == e2->v2 || e1->v1 == e2->v2 || e1->v2 == e2->v1)
341                                 continue;
342
343                         sub_v3_v3v3(vec1, BM_OtherEdgeVert(e1, v1)->co, v1->co);
344                         sub_v3_v3v3(vec2, BM_OtherEdgeVert(e2, v2)->co, v2->co);
345
346                         angle = fabs(angle_v3v3(vec1, vec2));
347
348                         if (angle < minangle) {
349                                 minangle = angle;
350                                 cure1 = e1;
351                                 cure2 = e2;
352                         }
353                 }
354         }
355
356         if (!cure1 || !cure1->loop || !cure2->loop) {
357                 /*just return 1.0 in this case*/
358                 return 1.0f;
359         }
360
361         /*assumtions
362         
363           we assume a 2-manifold mesh here.  if at any time this isn't the case,
364           e.g. a hole or an edge with more then 2 faces around it, we um ignore
365           that edge I guess, and try to make the algorithm go around as necassary.*/
366
367         l1 = cure1->loop;
368         l2 = cure2->loop;
369
370         lastv1 = l1->v == v1 ? ((BMLoop*)l1->head.next)->v : ((BMLoop*)l1->head.prev)->v;
371         lastv2 = l2->v == v2 ? ((BMLoop*)l2->head.next)->v : ((BMLoop*)l2->head.prev)->v;
372
373         /*we can only provide meaningful comparisons if v1 and v2 have the same valence*/
374         if (BM_Vert_EdgeCount(v1) != BM_Vert_EdgeCount(v2))
375                 return 1.0f; /*full mismatch*/
376
377         gh = BLI_ghash_new(BLI_ghashutil_ptrhash, BLI_ghashutil_ptrcmp);
378
379 #define SPUSH(s, d, vt, lv, e)\
380         if (BLI_array_count(s) <= lvl) BLI_array_growone(s);\
381         memset((s+lvl), 0, sizeof(*s));\
382         s[lvl].depth = d;\
383         s[lvl].v = vt;\
384         s[lvl].cure = e;\
385         s[lvl].lastv = lv;\
386         s[lvl].valence = BM_Vert_EdgeCount(vt);\
387
388         lvl = 0;
389
390         SPUSH(stack1, 0, v1, lastv1, cure1);
391         SPUSH(stack2, 0, v2, lastv2, cure2);
392
393         BLI_srand( BLI_rand() ); /* random seed */
394
395         lvl = 1;
396         while (lvl) {
397                 int term = 0;
398                 walklist *s1 = stack1 + lvl - 1, *s2 = stack2 + lvl - 1;
399
400                 /*pop from the stack*/
401                 lvl--;
402
403                 if (s1->curl && s1->curl->e == s1->cure)
404                         term = 1;
405                 if (s2->curl && s2->curl->e == s2->cure)
406                         term = 1;
407
408                 /*find next case to do*/
409                 if (!s1->curl)
410                         s1->curl = s1->cure->loop;
411                 if (!s2->curl) {
412                         float no1[3], no2[3], angle;
413                         int wind1, wind2;
414                         
415                         s2->curl = s2->cure->loop;
416
417                         /*find which of two possible faces to use*/
418                         l1 = BM_OtherFaceLoop(s1->curl->e, s1->curl->f, s1->lastv);
419                         l2 = BM_OtherFaceLoop(s2->curl->e, s2->curl->f, s2->lastv);
420
421                         if (l1->v == s2->lastv) {
422                                 l1 = (BMLoop*) l1->head.next;
423                                 if (l1->v == s2->v)
424                                         l1 = (BMLoop*) l1->head.prev->prev;
425                         } else if (l1->v == s2->v) {
426                                 l1 = (BMLoop*) l1->head.next;
427                                 if (l1->v == s2->lastv)
428                                         l1 = (BMLoop*) l1->head.prev->prev;
429                         }
430
431                         if (l2->v == s2->lastv) {
432                                 l2 = (BMLoop*) l2->head.next;
433                                 if (l2->v == s2->v)
434                                         l2 = (BMLoop*) l2->head.prev->prev;
435                         } else if (l2->v == s2->v) {
436                                 l2 = (BMLoop*) l2->head.next;
437                                 if (l2->v == s2->lastv)
438                                         l2 = (BMLoop*) l2->head.prev->prev;
439                         }
440
441                         wind1 = winding(s1->v->co, s1->lastv->co, l1->v->co);
442
443                         wind2 = winding(s2->v->co, s2->lastv->co, l2->v->co);
444                         
445                         /*if angle between the two adjacent faces is greater then 90 degrees,
446                           we need to flip wind2*/
447                         l1 = l2;
448                         l2 = s2->curl->radial.next->data;
449                         l2 = BM_OtherFaceLoop(l2->e, l2->f, s2->lastv);
450                         
451                         if (l2->v == s2->lastv) {
452                                 l2 = (BMLoop*) l2->head.next;
453                                 if (l2->v == s2->v)
454                                         l2 = (BMLoop*) l2->head.prev->prev;
455                         } else if (l2->v == s2->v) {
456                                 l2 = (BMLoop*) l2->head.next;
457                                 if (l2->v == s2->lastv)
458                                         l2 = (BMLoop*) l2->head.prev->prev;
459                         }
460
461                         normal_tri_v3(no1, s2->v->co, s2->lastv->co, l1->v->co);
462                         normal_tri_v3(no2, s2->v->co, s2->lastv->co, l2->v->co);
463                         
464                         /*enforce identical winding as no1*/
465                         mul_v3_fl(no2, -1.0);
466
467                         angle = angle_v3v3(no1, no2);
468                         if (angle > M_PI/2 - FLT_EPSILON*2)
469                                 wind2 = !wind2;
470
471                         if (wind1 == wind2)
472                                 s2->curl = s2->curl->radial.next->data;
473                 }
474
475                 /*handle termination cases of having already looped through all child
476                   nodes, or the valence mismatching between v1 and v2, or we hit max
477                   recursion depth*/
478                 term |= s1->valence != s2->valence || lvl+1 > maxlevel;
479                 term |= s1->curl->radial.next->data == (BMLoop*)l1;
480                 term |= s2->curl->radial.next->data == (BMLoop*)l2;
481
482                 if (!term) {
483                         lastv1 = s1->v;
484                         lastv2 = s2->v;
485                         v1 = BM_OtherEdgeVert(s1->curl->e, lastv1);
486                         v2 = BM_OtherEdgeVert(s2->curl->e, lastv2);
487                         
488                         e1 = s1->curl->e;
489                         e2 = s2->curl->e;
490
491                         if (!BLI_ghash_haskey(gh, v1) && !BLI_ghash_haskey(gh, v2)) {
492                                 /*repush the current stack item*/
493                                 lvl++;
494                                 
495                                 //if (maxlevel % 2 == 0) {
496                                         BLI_ghash_insert(gh, v1, NULL);
497                                         BLI_ghash_insert(gh, v2, NULL);
498                                 //}
499
500                                 /*now push the child node*/
501                                 SPUSH(stack1, lvl, v1, lastv1, e1);
502                                 SPUSH(stack2, lvl, v2, lastv2, e2);
503
504                                 lvl++;
505
506                                 s1 = stack1 + lvl - 2;
507                                 s2 = stack2 + lvl - 2;
508                         }
509
510                         s1->curl = s1->curl->v == s1->v ? (BMLoop*) s1->curl->head.prev : (BMLoop*) s1->curl->head.next;
511                         s2->curl = s2->curl->v == s2->v ? (BMLoop*) s2->curl->head.prev : (BMLoop*) s2->curl->head.next;
512                 
513                         s1->curl = (BMLoop*) s1->curl->radial.next->data;
514                         s2->curl = (BMLoop*) s2->curl->radial.next->data;
515                 }
516
517 #define WADD(stack, s)\
518                 if (lvl) {/*silly attempt to make this non-commutative: randomize\
519                               how much this particular weight adds to the total*/\
520                         stack[lvl-1].r += r;\
521                         s->w *= r;\
522                         stack[lvl-1].totwalked++;\
523                         stack[lvl-1].w += s->w;\
524                 }
525
526                 /*if no next case to do, update parent weight*/
527                 if (term) {
528                         float r = 0.8f + BLI_frand()*0.2f - FLT_EPSILON;
529
530                         if (s1->totwalked) {
531                                 s1->w /= s1->r;
532                         } else
533                                 s1->w = s1->valence == s2->valence ? 1.0f : 0.0f;
534
535                         WADD(stack1, s1);
536
537                         if (s2->totwalked) {
538                                 s2->w /= s2->r;
539                         } else
540                                 s2->w = s1->valence == s2->valence ? 1.0f : 0.0f;
541                         
542                         WADD(stack2, s2);
543
544                         /*apply additional penalty to weight mismatch*/
545                         if (s2->w != s1->w)
546                                 s2->w *= 0.8f;
547                 }
548         }
549
550         w = (stack1[0].w + stack2[0].w)*0.5f;
551
552         BLI_array_free(stack1);
553         BLI_array_free(stack2);
554
555         BLI_ghash_free(gh, NULL, NULL);
556
557         return 1.0f - w;
558 }
559
560 static void vertsearchcallback_topo(void *userdata, int index, const float *co, BVHTreeNearest *hit)
561 {
562         BMBVHTree *tree = userdata;
563         BMLoop **ls = tree->em->looptris[index];
564         int i;
565         float dist, maxdist, vec[3], w;
566
567         maxdist = tree->maxdist;
568
569         for (i=0; i<3; i++) {
570                 float dis;
571
572                 if (BLI_ghash_haskey(tree->gh, ls[i]->v))
573                         continue;
574
575                 sub_v3_v3v3(vec, tree->co, ls[i]->v->co);
576                 dis = dot_v3v3(vec, vec);
577
578                 w = topo_compare(tree->em->bm, tree->v, ls[i]->v, tree->curtag++);
579                 
580                 if (w < tree->curw-FLT_EPSILON*4) {
581                         tree->curw = w;
582                         tree->curv = ls[i]->v;
583                    
584                         sub_v3_v3v3(vec, tree->co, ls[i]->v->co);
585                         tree->curd = dot_v3v3(vec, vec);
586
587                         /*we deliberately check for equality using (smallest possible float)*4 
588                           comparison factor, to always prefer distance in cases of verts really
589                           close to each other*/
590                 } else if (fabs(tree->curw - w) < FLT_EPSILON*4) {
591                         /*if w is equal to hitex->curw, sort by distance*/
592                         sub_v3_v3v3(vec, tree->co, ls[i]->v->co);
593                         dis = dot_v3v3(vec, vec);
594
595                         if (dis < tree->curd) {
596                                 tree->curd = dis;
597                                 tree->curv = ls[i]->v;
598                         }
599                 }
600
601                 BLI_ghash_insert(tree->gh, ls[i]->v, NULL);
602         }
603 }
604
605 BMVert *BMBVH_FindClosestVertTopo(BMBVHTree *tree, float *co, float maxdist, BMVert *sourcev)
606 {
607         BVHTreeNearest hit;
608         BMVert *v;
609         BMIter iter;
610
611         memset(&hit, 0, sizeof(hit));
612
613         VECCOPY(hit.co, co);
614         VECCOPY(tree->co, co);
615         hit.index = -1;
616         hit.dist = maxdist;
617
618         tree->curw = FLT_MAX;
619         tree->curd = FLT_MAX;
620         tree->curv = NULL;
621         tree->curtag = 1;
622
623         tree->gh = BLI_ghash_new(BLI_ghashutil_ptrhash, BLI_ghashutil_ptrcmp);
624
625         tree->maxdist = maxdist;
626         tree->v = sourcev;
627
628         BLI_bvhtree_find_nearest(tree->tree, co, &hit, vertsearchcallback_topo, tree);
629         
630         BLI_ghash_free(tree->gh, NULL, NULL);
631         tree->gh = NULL;
632
633         return tree->curv;
634 }
635
636
637 #if 0 //BMESH_TODO: not implemented yet
638 int BMBVH_VertVisible(BMBVHTree *tree, BMEdge *e, RegionView3D *r3d)
639 {
640
641 }
642 #endif
643
644 static BMFace *edge_ray_cast(BMBVHTree *tree, float *co, float *dir, float *hitout, BMEdge *e)
645 {
646         BMFace *f = BMBVH_RayCast(tree, co, dir, hitout);
647         
648         if (f && BM_Edge_In_Face(f, e))
649                 return NULL;
650
651         return f;
652 }
653
654 int BMBVH_EdgeVisible(BMBVHTree *tree, BMEdge *e, RegionView3D *r3d, Object *obedit)
655 {
656         BMFace *f;
657         float co1[3], co2[3], co3[3], dir1[4], dir2[4], dir3[4];
658         float origin[3], invmat[4][4];
659         float epsilon = 0.01f; 
660         
661         VECCOPY(origin, r3d->viewinv[3]);
662         invert_m4_m4(invmat, obedit->obmat);
663         mul_m4_v3(invmat, origin);
664
665         VECCOPY(co1, e->v1->co);
666         add_v3_v3v3(co2, e->v1->co, e->v2->co);
667         mul_v3_fl(co2, 0.5f);
668         VECCOPY(co3, e->v2->co);
669         
670         /*ok, idea is to generate rays going from the camera origin to the 
671           three points on the edge (v1, mid, v2)*/
672         sub_v3_v3v3(dir1, origin, co1);
673         sub_v3_v3v3(dir2, origin, co2);
674         sub_v3_v3v3(dir3, origin, co3);
675         
676         normalize_v3(dir1);
677         normalize_v3(dir2);
678         normalize_v3(dir3);
679
680         mul_v3_fl(dir1, epsilon);
681         mul_v3_fl(dir2, epsilon);
682         mul_v3_fl(dir3, epsilon);
683         
684         /*offset coordinates slightly along view vectors, to avoid
685           hitting the faces that own the edge.*/
686         add_v3_v3v3(co1, co1, dir1);
687         add_v3_v3v3(co2, co2, dir2);
688         add_v3_v3v3(co3, co3, dir3);
689
690         normalize_v3(dir1);
691         normalize_v3(dir2);
692         normalize_v3(dir3);
693
694         /*do three samplings: left, middle, right*/
695         f = edge_ray_cast(tree, co1, dir1, NULL, e);
696         if (f && !edge_ray_cast(tree, co2, dir2, NULL, e))
697                 return 1;
698         else if (f && !edge_ray_cast(tree, co3, dir3, NULL, e))
699                 return 1;
700         else if (!f)
701                 return 1;
702
703         return 0;
704 }