Blender  V2.59
bvhutils.c
Go to the documentation of this file.
00001 /*
00002  *
00003  * $Id: bvhutils.c 37787 2011-06-24 05:34:03Z campbellbarton $
00004  *
00005  * ***** BEGIN GPL LICENSE BLOCK *****
00006  *
00007  * This program is free software; you can redistribute it and/or
00008  * modify it under the terms of the GNU General Public License
00009  * as published by the Free Software Foundation; either version 2
00010  * of the License, or (at your option) any later version.
00011  *
00012  * This program is distributed in the hope that it will be useful,
00013  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00014  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00015  * GNU General Public License for more details.
00016  *
00017  * You should have received a copy of the GNU General Public License
00018  * along with this program; if not, write to the Free Software Foundation,
00019  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
00020  *
00021  * The Original Code is Copyright (C) Blender Foundation.
00022  * All rights reserved.
00023  *
00024  * The Original Code is: all of this file.
00025  *
00026  * Contributor(s): Andr Pinto.
00027  *
00028  * ***** END GPL LICENSE BLOCK *****
00029  */
00030 
00035 #include <stdio.h>
00036 #include <string.h>
00037 #include <math.h>
00038 #include <assert.h>
00039 
00040 #include "DNA_meshdata_types.h"
00041 
00042 #include "BLI_editVert.h"
00043 #include "BLI_utildefines.h"
00044 
00045 #include "BKE_DerivedMesh.h"
00046 
00047 
00048 #include "BLI_math.h"
00049 #include "MEM_guardedalloc.h"
00050 
00051 /* Math stuff for ray casting on mesh faces and for nearest surface */
00052 
00053 static float ray_tri_intersection(const BVHTreeRay *ray, const float UNUSED(m_dist), const float *v0, const float *v1, const float *v2)
00054 {
00055         float dist;
00056 
00057         if(isect_ray_tri_v3((float*)ray->origin, (float*)ray->direction, (float*)v0, (float*)v1, (float*)v2, &dist, NULL))
00058                 return dist;
00059 
00060         return FLT_MAX;
00061 }
00062 
00063 static float sphereray_tri_intersection(const BVHTreeRay *ray, float radius, const float m_dist, const float *v0, const float *v1, const float *v2)
00064 {
00065         
00066         float idist;
00067         float p1[3];
00068         float plane_normal[3], hit_point[3];
00069 
00070         normal_tri_v3( plane_normal,(float*)v0, (float*)v1, (float*)v2);
00071 
00072         VECADDFAC( p1, ray->origin, ray->direction, m_dist);
00073         if(isect_sweeping_sphere_tri_v3((float*)ray->origin, p1, radius, (float*)v0, (float*)v1, (float*)v2, &idist, hit_point))
00074         {
00075                 return idist * m_dist;
00076         }
00077 
00078         return FLT_MAX;
00079 }
00080 
00081 
00082 /*
00083  * Function adapted from David Eberly's distance tools (LGPL)
00084  * http://www.geometrictools.com/LibFoundation/Distance/Distance.html
00085  */
00086 static float nearest_point_in_tri_surface(const float *v0,const float *v1,const float *v2,const float *p, int *v, int *e, float *nearest )
00087 {
00088         float diff[3];
00089         float e0[3];
00090         float e1[3];
00091         float A00;
00092         float A01;
00093         float A11;
00094         float B0;
00095         float B1;
00096         float C;
00097         float Det;
00098         float S;
00099         float T;
00100         float sqrDist;
00101         int lv = -1, le = -1;
00102         
00103         VECSUB(diff, v0, p);
00104         VECSUB(e0, v1, v0);
00105         VECSUB(e1, v2, v0);
00106         
00107         A00 = INPR ( e0, e0 );
00108         A01 = INPR( e0, e1 );
00109         A11 = INPR ( e1, e1 );
00110         B0 = INPR( diff, e0 );
00111         B1 = INPR( diff, e1 );
00112         C = INPR( diff, diff );
00113         Det = fabs( A00 * A11 - A01 * A01 );
00114         S = A01 * B1 - A11 * B0;
00115         T = A01 * B0 - A00 * B1;
00116 
00117         if ( S + T <= Det )
00118         {
00119                 if ( S < 0.0f )
00120                 {
00121                         if ( T < 0.0f )  // Region 4
00122                         {
00123                                 if ( B0 < 0.0f )
00124                                 {
00125                                         T = 0.0f;
00126                                         if ( -B0 >= A00 )
00127                                         {
00128                                                 S = (float)1.0;
00129                                                 sqrDist = A00 + 2.0f * B0 + C;
00130                                                 lv = 1;
00131                                         }
00132                                         else
00133                                         {
00134                                                 if(fabsf(A00) > FLT_EPSILON)
00135                                                         S = -B0/A00;
00136                                                 else
00137                                                         S = 0.0f;
00138                                                 sqrDist = B0 * S + C;
00139                                                 le = 0;
00140                                         }
00141                                 }
00142                                 else
00143                                 {
00144                                         S = 0.0f;
00145                                         if ( B1 >= 0.0f )
00146                                         {
00147                                                 T = 0.0f;
00148                                                 sqrDist = C;
00149                                                 lv = 0;
00150                                         }
00151                                         else if ( -B1 >= A11 )
00152                                         {
00153                                                 T = 1.0f;
00154                                                 sqrDist = A11 + 2.0f * B1 + C;
00155                                                 lv = 2;
00156                                         }
00157                                         else
00158                                         {
00159                                                 if(fabsf(A11) > FLT_EPSILON)
00160                                                         T = -B1 / A11;
00161                                                 else
00162                                                         T = 0.0f;
00163                                                 sqrDist = B1 * T + C;
00164                                                 le = 1;
00165                                         }
00166                                 }
00167                         }
00168                         else  // Region 3
00169                         {
00170                                 S = 0.0f;
00171                                 if ( B1 >= 0.0f )
00172                                 {
00173                                         T = 0.0f;
00174                                         sqrDist = C;
00175                                         lv = 0;
00176                                 }
00177                                 else if ( -B1 >= A11 )
00178                                 {
00179                                         T = 1.0f;
00180                                         sqrDist = A11 + 2.0f * B1 + C;
00181                                         lv = 2;
00182                                 }
00183                                 else
00184                                 {
00185                                         if(fabsf(A11) > FLT_EPSILON)
00186                                                 T = -B1 / A11;
00187                                         else
00188                                                 T = 0.0;
00189                                         sqrDist = B1 * T + C;
00190                                         le = 1;
00191                                 }
00192                         }
00193                 }
00194                 else if ( T < 0.0f )  // Region 5
00195                 {
00196                         T = 0.0f;
00197                         if ( B0 >= 0.0f )
00198                         {
00199                                 S = 0.0f;
00200                                 sqrDist = C;
00201                                 lv = 0;
00202                         }
00203                         else if ( -B0 >= A00 )
00204                         {
00205                                 S = 1.0f;
00206                                 sqrDist = A00 + 2.0f * B0 + C;
00207                                 lv = 1;
00208                         }
00209                         else
00210                         {
00211                                 if(fabsf(A00) > FLT_EPSILON)
00212                                         S = -B0 / A00;
00213                                 else
00214                                         S = 0.0f;
00215                                 sqrDist = B0 * S + C;
00216                                 le = 0;
00217                         }
00218                 }
00219                 else  // Region 0
00220                 {
00221                         // Minimum at interior lv
00222                         float invDet;
00223                         if(fabsf(Det) > FLT_EPSILON)
00224                                 invDet = 1.0f / Det;
00225                         else
00226                                 invDet = 0.0f;
00227                         S *= invDet;
00228                         T *= invDet;
00229                         sqrDist = S * ( A00 * S + A01 * T + 2.0f * B0) +
00230                                         T * ( A01 * S + A11 * T + 2.0f * B1 ) + C;
00231                 }
00232         }
00233         else
00234         {
00235                 float tmp0, tmp1, numer, denom;
00236 
00237                 if ( S < 0.0f )  // Region 2
00238                 {
00239                         tmp0 = A01 + B0;
00240                         tmp1 = A11 + B1;
00241                         if ( tmp1 > tmp0 )
00242                         {
00243                                 numer = tmp1 - tmp0;
00244                                 denom = A00 - 2.0f * A01 + A11;
00245                                 if ( numer >= denom )
00246                                 {
00247                                         S = 1.0f;
00248                                         T = 0.0f;
00249                                         sqrDist = A00 + 2.0f * B0 + C;
00250                                         lv = 1;
00251                                 }
00252                                 else
00253                                 {
00254                                         if(fabsf(denom) > FLT_EPSILON)
00255                                                 S = numer / denom;
00256                                         else
00257                                                 S = 0.0f;
00258                                         T = 1.0f - S;
00259                                         sqrDist = S * ( A00 * S + A01 * T + 2.0f * B0 ) +
00260                                                         T * ( A01 * S + A11 * T + 2.0f * B1 ) + C;
00261                                         le = 2;
00262                                 }
00263                         }
00264                         else
00265                         {
00266                                 S = 0.0f;
00267                                 if ( tmp1 <= 0.0f )
00268                                 {
00269                                         T = 1.0f;
00270                                         sqrDist = A11 + 2.0f * B1 + C;
00271                                         lv = 2;
00272                                 }
00273                                 else if ( B1 >= 0.0f )
00274                                 {
00275                                         T = 0.0f;
00276                                         sqrDist = C;
00277                                         lv = 0;
00278                                 }
00279                                 else
00280                                 {
00281                                         if(fabsf(A11) > FLT_EPSILON)
00282                                                 T = -B1 / A11;
00283                                         else
00284                                                 T = 0.0f;
00285                                         sqrDist = B1 * T + C;
00286                                         le = 1;
00287                                 }
00288                         }
00289                 }
00290                 else if ( T < 0.0f )  // Region 6
00291                 {
00292                         tmp0 = A01 + B1;
00293                         tmp1 = A00 + B0;
00294                         if ( tmp1 > tmp0 )
00295                         {
00296                                 numer = tmp1 - tmp0;
00297                                 denom = A00 - 2.0f * A01 + A11;
00298                                 if ( numer >= denom )
00299                                 {
00300                                         T = 1.0f;
00301                                         S = 0.0f;
00302                                         sqrDist = A11 + 2.0f * B1 + C;
00303                                         lv = 2;
00304                                 }
00305                                 else
00306                                 {
00307                                         if(fabsf(denom) > FLT_EPSILON)
00308                                                 T = numer / denom;
00309                                         else
00310                                                 T = 0.0f;
00311                                         S = 1.0f - T;
00312                                         sqrDist = S * ( A00 * S + A01 * T + 2.0f * B0 ) +
00313                                                         T * ( A01 * S + A11 * T + 2.0f * B1 ) + C;
00314                                         le = 2;
00315                                 }
00316                         }
00317                         else
00318                         {
00319                                 T = 0.0f;
00320                                 if ( tmp1 <= 0.0f )
00321                                 {
00322                                         S = 1.0f;
00323                                         sqrDist = A00 + 2.0f * B0 + C;
00324                                         lv = 1;
00325                                 }
00326                                 else if ( B0 >= 0.0f )
00327                                 {
00328                                         S = 0.0f;
00329                                         sqrDist = C;
00330                                         lv = 0;
00331                                 }
00332                                 else
00333                                 {
00334                                         if(fabsf(A00) > FLT_EPSILON)
00335                                                 S = -B0 / A00;
00336                                         else
00337                                                 S = 0.0f;
00338                                         sqrDist = B0 * S + C;
00339                                         le = 0;
00340                                 }
00341                         }
00342                 }
00343                 else  // Region 1
00344                 {
00345                         numer = A11 + B1 - A01 - B0;
00346                         if ( numer <= 0.0f )
00347                         {
00348                                 S = 0.0f;
00349                                 T = 1.0f;
00350                                 sqrDist = A11 + 2.0f * B1 + C;
00351                                 lv = 2;
00352                         }
00353                         else
00354                         {
00355                                 denom = A00 - 2.0f * A01 + A11;
00356                                 if ( numer >= denom )
00357                                 {
00358                                         S = 1.0f;
00359                                         T = 0.0f;
00360                                         sqrDist = A00 + 2.0f * B0 + C;
00361                                         lv = 1;
00362                                 }
00363                                 else
00364                                 {
00365                                         if(fabsf(denom) > FLT_EPSILON)
00366                                                 S = numer / denom;
00367                                         else
00368                                                 S = 0.0f;
00369                                         T = 1.0f - S;
00370                                         sqrDist = S * ( A00 * S + A01 * T + 2.0f * B0 ) +
00371                                                         T * ( A01 * S + A11 * T + 2.0f * B1 ) + C;
00372                                         le = 2;
00373                                 }
00374                         }
00375                 }
00376         }
00377 
00378         // Account for numerical round-off error
00379         if ( sqrDist < FLT_EPSILON )
00380                 sqrDist = 0.0f;
00381         
00382         {
00383                 float w[3], x[3], y[3], z[3];
00384                 VECCOPY(w, v0);
00385                 VECCOPY(x, e0);
00386                 mul_v3_fl(x, S);
00387                 VECCOPY(y, e1);
00388                 mul_v3_fl(y, T);
00389                 VECADD(z, w, x);
00390                 VECADD(z, z, y);
00391                 //VECSUB(d, p, z);
00392                 VECCOPY(nearest, z);
00393                 // d = p - ( v0 + S * e0 + T * e1 );
00394         }
00395         *v = lv;
00396         *e = le;
00397 
00398         return sqrDist;
00399 }
00400 
00401 
00402 /*
00403  * BVH from meshs callbacks
00404  */
00405 
00406 // Callback to bvh tree nearest point. The tree must bust have been built using bvhtree_from_mesh_faces.
00407 // userdata must be a BVHMeshCallbackUserdata built from the same mesh as the tree.
00408 static void mesh_faces_nearest_point(void *userdata, int index, const float *co, BVHTreeNearest *nearest)
00409 {
00410         const BVHTreeFromMesh *data = (BVHTreeFromMesh*) userdata;
00411         MVert *vert     = data->vert;
00412         MFace *face = data->face + index;
00413 
00414         float *t0, *t1, *t2, *t3;
00415         t0 = vert[ face->v1 ].co;
00416         t1 = vert[ face->v2 ].co;
00417         t2 = vert[ face->v3 ].co;
00418         t3 = face->v4 ? vert[ face->v4].co : NULL;
00419 
00420         
00421         do
00422         {       
00423                 float nearest_tmp[3], dist;
00424                 int vertex, edge;
00425                 
00426                 dist = nearest_point_in_tri_surface(t0, t1, t2, co, &vertex, &edge, nearest_tmp);
00427                 if(dist < nearest->dist)
00428                 {
00429                         nearest->index = index;
00430                         nearest->dist = dist;
00431                         VECCOPY(nearest->co, nearest_tmp);
00432                         normal_tri_v3( nearest->no,t0, t1, t2);
00433                 }
00434 
00435                 t1 = t2;
00436                 t2 = t3;
00437                 t3 = NULL;
00438 
00439         } while(t2);
00440 }
00441 
00442 // Callback to bvh tree raycast. The tree must bust have been built using bvhtree_from_mesh_faces.
00443 // userdata must be a BVHMeshCallbackUserdata built from the same mesh as the tree.
00444 static void mesh_faces_spherecast(void *userdata, int index, const BVHTreeRay *ray, BVHTreeRayHit *hit)
00445 {
00446         const BVHTreeFromMesh *data = (BVHTreeFromMesh*) userdata;
00447         MVert *vert     = data->vert;
00448         MFace *face = data->face + index;
00449 
00450         float *t0, *t1, *t2, *t3;
00451         t0 = vert[ face->v1 ].co;
00452         t1 = vert[ face->v2 ].co;
00453         t2 = vert[ face->v3 ].co;
00454         t3 = face->v4 ? vert[ face->v4].co : NULL;
00455 
00456         
00457         do
00458         {       
00459                 float dist;
00460                 if(data->sphere_radius == 0.0f)
00461                         dist = ray_tri_intersection(ray, hit->dist, t0, t1, t2);
00462                 else
00463                         dist = sphereray_tri_intersection(ray, data->sphere_radius, hit->dist, t0, t1, t2);
00464 
00465                 if(dist >= 0 && dist < hit->dist)
00466                 {
00467                         hit->index = index;
00468                         hit->dist = dist;
00469                         VECADDFAC(hit->co, ray->origin, ray->direction, dist);
00470 
00471                         normal_tri_v3( hit->no,t0, t1, t2);
00472                 }
00473 
00474                 t1 = t2;
00475                 t2 = t3;
00476                 t3 = NULL;
00477 
00478         } while(t2);
00479 }
00480 
00481 // Callback to bvh tree nearest point. The tree must bust have been built using bvhtree_from_mesh_edges.
00482 // userdata must be a BVHMeshCallbackUserdata built from the same mesh as the tree.
00483 static void mesh_edges_nearest_point(void *userdata, int index, const float *co, BVHTreeNearest *nearest)
00484 {
00485         const BVHTreeFromMesh *data = (BVHTreeFromMesh*) userdata;
00486         MVert *vert     = data->vert;
00487         MEdge *edge = data->edge + index;
00488         float nearest_tmp[3], dist;
00489 
00490         float *t0, *t1;
00491         t0 = vert[ edge->v1 ].co;
00492         t1 = vert[ edge->v2 ].co;
00493         
00494         // NOTE: casts to "float*" here are due to co being "const float*"
00495         closest_to_line_segment_v3(nearest_tmp, (float*)co, t0, t1);
00496         dist = len_squared_v3v3(nearest_tmp, (float*)co);
00497         
00498         if(dist < nearest->dist)
00499         {
00500                 nearest->index = index;
00501                 nearest->dist = dist;
00502                 VECCOPY(nearest->co, nearest_tmp);
00503                 sub_v3_v3v3(nearest->no, t0, t1);
00504                 normalize_v3(nearest->no);
00505         }
00506 }
00507 
00508 /*
00509  * BVH builders
00510  */
00511 // Builds a bvh tree.. where nodes are the vertexs of the given mesh
00512 BVHTree* bvhtree_from_mesh_verts(BVHTreeFromMesh *data, DerivedMesh *mesh, float epsilon, int tree_type, int axis)
00513 {
00514         BVHTree *tree = bvhcache_find(&mesh->bvhCache, BVHTREE_FROM_VERTICES);
00515 
00516         //Not in cache
00517         if(tree == NULL)
00518         {
00519                 int i;
00520                 int numVerts= mesh->getNumVerts(mesh);
00521                 MVert *vert     = mesh->getVertDataArray(mesh, CD_MVERT);
00522 
00523                 if(vert != NULL)
00524                 {
00525                         tree = BLI_bvhtree_new(numVerts, epsilon, tree_type, axis);
00526 
00527                         if(tree != NULL)
00528                         {
00529                                 for(i = 0; i < numVerts; i++)
00530                                         BLI_bvhtree_insert(tree, i, vert[i].co, 1);
00531 
00532                                 BLI_bvhtree_balance(tree);
00533 
00534                                 //Save on cache for later use
00535 //                              printf("BVHTree built and saved on cache\n");
00536                                 bvhcache_insert(&mesh->bvhCache, tree, BVHTREE_FROM_VERTICES);
00537                         }
00538                 }
00539         }
00540         else
00541         {
00542 //              printf("BVHTree is already build, using cached tree\n");
00543         }
00544 
00545 
00546         //Setup BVHTreeFromMesh
00547         memset(data, 0, sizeof(*data));
00548         data->tree = tree;
00549 
00550         if(data->tree)
00551         {
00552                 data->cached = TRUE;
00553 
00554                 //a NULL nearest callback works fine
00555                 //remeber the min distance to point is the same as the min distance to BV of point
00556                 data->nearest_callback = NULL;
00557                 data->raycast_callback = NULL;
00558 
00559                 data->mesh = mesh;
00560                 data->vert = mesh->getVertDataArray(mesh, CD_MVERT);
00561                 data->face = mesh->getFaceDataArray(mesh, CD_MFACE);
00562 
00563                 data->sphere_radius = epsilon;
00564         }
00565 
00566         return data->tree;
00567 }
00568 
00569 // Builds a bvh tree.. where nodes are the faces of the given mesh.
00570 BVHTree* bvhtree_from_mesh_faces(BVHTreeFromMesh *data, DerivedMesh *mesh, float epsilon, int tree_type, int axis)
00571 {
00572         BVHTree *tree = bvhcache_find(&mesh->bvhCache, BVHTREE_FROM_FACES);
00573 
00574         //Not in cache
00575         if(tree == NULL)
00576         {
00577                 int i;
00578                 int numFaces= mesh->getNumFaces(mesh);
00579                 MVert *vert     = mesh->getVertDataArray(mesh, CD_MVERT);
00580                 MFace *face = mesh->getFaceDataArray(mesh, CD_MFACE);
00581 
00582                 if(vert != NULL && face != NULL)
00583                 {
00584                         /* Create a bvh-tree of the given target */
00585                         tree = BLI_bvhtree_new(numFaces, epsilon, tree_type, axis);
00586                         if(tree != NULL)
00587                         {
00588                                 /* XXX, for snap only, em & dm are assumed to be aligned, since dm is the em's cage */
00589                                 EditMesh *em= data->em_evil;
00590                                 if(em) {
00591                                         EditFace *efa= em->faces.first;
00592                                         for(i = 0; i < numFaces; i++, efa= efa->next) {
00593                                                 if(!(efa->f & 1) && efa->h==0 && !((efa->v1->f&1)+(efa->v2->f&1)+(efa->v3->f&1)+(efa->v4?efa->v4->f&1:0))) {
00594                                                         float co[4][3];
00595                                                         VECCOPY(co[0], vert[ face[i].v1 ].co);
00596                                                         VECCOPY(co[1], vert[ face[i].v2 ].co);
00597                                                         VECCOPY(co[2], vert[ face[i].v3 ].co);
00598                                                         if(face[i].v4)
00599                                                                 VECCOPY(co[3], vert[ face[i].v4 ].co);
00600                                         
00601                                                         BLI_bvhtree_insert(tree, i, co[0], face[i].v4 ? 4 : 3);
00602                                                 }
00603                                         }
00604                                 }
00605                                 else {
00606                                         for(i = 0; i < numFaces; i++) {
00607                                                 float co[4][3];
00608                                                 VECCOPY(co[0], vert[ face[i].v1 ].co);
00609                                                 VECCOPY(co[1], vert[ face[i].v2 ].co);
00610                                                 VECCOPY(co[2], vert[ face[i].v3 ].co);
00611                                                 if(face[i].v4)
00612                                                         VECCOPY(co[3], vert[ face[i].v4 ].co);
00613                                 
00614                                                 BLI_bvhtree_insert(tree, i, co[0], face[i].v4 ? 4 : 3);
00615                                         }
00616                                 }
00617                                 BLI_bvhtree_balance(tree);
00618 
00619                                 //Save on cache for later use
00620 //                              printf("BVHTree built and saved on cache\n");
00621                                 bvhcache_insert(&mesh->bvhCache, tree, BVHTREE_FROM_FACES);
00622                         }
00623                 }
00624         }
00625         else
00626         {
00627 //              printf("BVHTree is already build, using cached tree\n");
00628         }
00629 
00630 
00631         //Setup BVHTreeFromMesh
00632         memset(data, 0, sizeof(*data));
00633         data->tree = tree;
00634 
00635         if(data->tree)
00636         {
00637                 data->cached = TRUE;
00638 
00639                 data->nearest_callback = mesh_faces_nearest_point;
00640                 data->raycast_callback = mesh_faces_spherecast;
00641 
00642                 data->mesh = mesh;
00643                 data->vert = mesh->getVertDataArray(mesh, CD_MVERT);
00644                 data->face = mesh->getFaceDataArray(mesh, CD_MFACE);
00645 
00646                 data->sphere_radius = epsilon;
00647         }
00648         return data->tree;
00649 
00650 }
00651 
00652 // Builds a bvh tree.. where nodes are the faces of the given mesh.
00653 BVHTree* bvhtree_from_mesh_edges(BVHTreeFromMesh *data, DerivedMesh *mesh, float epsilon, int tree_type, int axis)
00654 {
00655         BVHTree *tree = bvhcache_find(&mesh->bvhCache, BVHTREE_FROM_EDGES);
00656 
00657         //Not in cache
00658         if(tree == NULL)
00659         {
00660                 int i;
00661                 int numEdges= mesh->getNumEdges(mesh);
00662                 MVert *vert     = mesh->getVertDataArray(mesh, CD_MVERT);
00663                 MEdge *edge = mesh->getEdgeDataArray(mesh, CD_MEDGE);
00664 
00665                 if(vert != NULL && edge != NULL)
00666                 {
00667                         /* Create a bvh-tree of the given target */
00668                         tree = BLI_bvhtree_new(numEdges, epsilon, tree_type, axis);
00669                         if(tree != NULL)
00670                         {
00671                                 for(i = 0; i < numEdges; i++)
00672                                 {
00673                                         float co[4][3];
00674                                         VECCOPY(co[0], vert[ edge[i].v1 ].co);
00675                                         VECCOPY(co[1], vert[ edge[i].v2 ].co);
00676                         
00677                                         BLI_bvhtree_insert(tree, i, co[0], 2);
00678                                 }
00679                                 BLI_bvhtree_balance(tree);
00680 
00681                                 //Save on cache for later use
00682 //                              printf("BVHTree built and saved on cache\n");
00683                                 bvhcache_insert(&mesh->bvhCache, tree, BVHTREE_FROM_EDGES);
00684                         }
00685                 }
00686         }
00687         else
00688         {
00689 //              printf("BVHTree is already build, using cached tree\n");
00690         }
00691 
00692 
00693         //Setup BVHTreeFromMesh
00694         memset(data, 0, sizeof(*data));
00695         data->tree = tree;
00696 
00697         if(data->tree)
00698         {
00699                 data->cached = TRUE;
00700 
00701                 data->nearest_callback = mesh_edges_nearest_point;
00702                 data->raycast_callback = NULL;
00703 
00704                 data->mesh = mesh;
00705                 data->vert = mesh->getVertDataArray(mesh, CD_MVERT);
00706                 data->edge = mesh->getEdgeDataArray(mesh, CD_MEDGE);
00707 
00708                 data->sphere_radius = epsilon;
00709         }
00710         return data->tree;
00711 
00712 }
00713 
00714 // Frees data allocated by a call to bvhtree_from_mesh_*.
00715 void free_bvhtree_from_mesh(struct BVHTreeFromMesh *data)
00716 {
00717         if(data->tree)
00718         {
00719                 if(!data->cached)
00720                         BLI_bvhtree_free(data->tree);
00721 
00722                 memset( data, 0, sizeof(*data) );
00723         }
00724 }
00725 
00726 
00727 /* BVHCache */
00728 typedef struct BVHCacheItem
00729 {
00730         int type;
00731         BVHTree *tree;
00732 
00733 } BVHCacheItem;
00734 
00735 static void bvhcacheitem_set_if_match(void *_cached, void *_search)
00736 {
00737         BVHCacheItem * cached = (BVHCacheItem *)_cached;
00738         BVHCacheItem * search = (BVHCacheItem *)_search;
00739 
00740         if(search->type == cached->type)
00741         {
00742                 search->tree = cached->tree;            
00743         }
00744 } 
00745 
00746 BVHTree *bvhcache_find(BVHCache *cache, int type)
00747 {
00748         BVHCacheItem item;
00749         item.type = type;
00750         item.tree = NULL;
00751 
00752         BLI_linklist_apply(*cache, bvhcacheitem_set_if_match, &item);
00753         return item.tree;
00754 }
00755 
00756 void bvhcache_insert(BVHCache *cache, BVHTree *tree, int type)
00757 {
00758         BVHCacheItem *item = NULL;
00759 
00760         assert( tree != NULL );
00761         assert( bvhcache_find(cache, type) == NULL );
00762 
00763         item = MEM_mallocN(sizeof(BVHCacheItem), "BVHCacheItem");
00764         assert( item != NULL );
00765 
00766         item->type = type;
00767         item->tree = tree;
00768 
00769         BLI_linklist_prepend( cache, item );
00770 }
00771 
00772 
00773 void bvhcache_init(BVHCache *cache)
00774 {
00775         *cache = NULL;
00776 }
00777 
00778 static void bvhcacheitem_free(void *_item)
00779 {
00780         BVHCacheItem *item = (BVHCacheItem *)_item;
00781 
00782         BLI_bvhtree_free(item->tree);
00783         MEM_freeN(item);
00784 }
00785 
00786 
00787 void bvhcache_free(BVHCache *cache)
00788 {
00789         BLI_linklist_free(*cache, (LinkNodeFreeFP)bvhcacheitem_free);
00790         *cache = NULL;
00791 }
00792 
00793