← Back to branch summary

~ubuntu-branches/ubuntu/trusty/blender/trusty

« back to all changes in this revision

Viewing changes to extern/eltopo/common/makelevelset2.cpp

  • Committer: Package Import Robot
  • Author(s): Jeremy Bicha
  • Date: 2013-03-06 12:08:47 UTC
  • mfrom: (1.5.1) (14.1.8 experimental)
  • Revision ID: package-import@ubuntu.com-20130306120847-frjfaryb2zrotwcg
Tags: 2.66a-1ubuntu1
https://launchpad.net/bugs/1076930
https://launchpad.net/bugs/1089256
https://launchpad.net/bugs/1052743
https://launchpad.net/bugs/999024
* Resynchronize with Debian (LP: #1076930, #1089256, #1052743, #999024,
  #1122888, #1147084)
* debian/control:
  - Lower build-depends on libavcodec-dev since we're not
    doing the libav9 transition in Ubuntu yet

Show diffs side-by-side

added added

removed removed

Lines of Context:
extern/eltopo/common/makelevelset2.cpp
1
 
#include "makelevelset2.h"
2
 
 
3
 
// find distance x0 is from segment x1-x2
4
 
static double point_segment_distance(const Vec2d &x0, const Vec2d &x1, const Vec2d &x2)
5
 
{
6
 
   Vec2d dx(x2-x1);
7
 
   double m2=mag2(dx);
8
 
   // find parameter value of closest point on segment
9
 
   double s12=(double)( dot(x2-x0, dx)/m2 );
10
 
   if(s12<0){
11
 
      s12=0;
12
 
   }else if(s12>1){
13
 
      s12=1;
14
 
   }
15
 
   // and find the distance
16
 
   return mag(s12*x1+(1-s12)*x2 - x0);
17
 
}
18
 
 
19
 
static void check_neighbour(const std::vector<Vec2ui> &edge, const std::vector<Vec2d> &x,
20
 
                            Array2d &phi, Array2i &closest_edge,
21
 
                            const Vec2d &gx, int i0, int j0, int i1, int j1)
22
 
{
23
 
   if(closest_edge(i0,j0)==closest_edge(i1,j1))
24
 
      return;
25
 
   if(closest_edge(i1,j1)>=0){
26
 
      unsigned int p, q; assign(edge[closest_edge(i1,j1)], p, q);
27
 
      double d=point_segment_distance(gx, x[p], x[q]);
28
 
      if(d<phi(i0,j0)){
29
 
         phi(i0,j0)=d;
30
 
         closest_edge(i0,j0)=closest_edge(i1,j1);
31
 
      }
32
 
   }
33
 
}
34
 
 
35
 
void make_level_set2(const std::vector<Vec2ui> &edge, const std::vector<Vec2d> &x,
36
 
                     const Vec2d &origin, double dx, int nx, int ny,
37
 
                     Array2d &phi)
38
 
{
39
 
   phi.resize(nx, ny);
40
 
   phi.assign((nx+ny)*dx); // upper bound on distance
41
 
   Array2i closest_edge(nx, ny, -1);
42
 
   Array2i intersection_count(nx, ny, 0); // intersection_count(i,j) is # of edge intersections in (i-1,i]x{j}
43
 
   // we begin by initializing distances near the mesh, and figuring out intersection counts
44
 
   Vec2d ijmin, ijmax;
45
 
   for(unsigned int e=0; e<edge.size(); ++e){
46
 
      unsigned int p, q; assign(edge[e], p, q);
47
 
      minmax((x[p]-origin)/dx, (x[q]-origin)/dx, ijmin, ijmax);
48
 
      // do distances
49
 
      int i0=clamp(int(ijmin[0])-1, 0, nx-1), i1=clamp(int(ijmax[0])+2, 0, nx-1);
50
 
      int j0=clamp(int(ijmin[1])-1, 0, ny-1), j1=clamp(int(ijmax[1])+2, 0, ny-1);
51
 
      for(int j=j0; j<=j1; ++j) for(int i=i0; i<=i1; ++i){
52
 
         Vec2d gx(i*dx+origin[0], j*dx+origin[1]);
53
 
         double d=point_segment_distance(gx, x[p], x[q]);
54
 
         if(d<phi(i,j)){
55
 
            phi(i,j)=d;
56
 
            closest_edge(i,j)=e;
57
 
         }
58
 
      }
59
 
      // and do intersection counts
60
 
      if(x[p][1]!=x[q][1]){ // if it's not a horizontal edge
61
 
         double fi0, fj0, fi1, fj1;
62
 
         if(x[p][1]<x[q][1]){
63
 
            fi0=((double)x[p][0]-origin[0])/dx;
64
 
            fj0=((double)x[p][1]-origin[1])/dx;
65
 
            fi1=((double)x[q][0]-origin[0])/dx;
66
 
            fj1=((double)x[q][1]-origin[1])/dx;
67
 
         }else{
68
 
            fi0=((double)x[q][0]-origin[0])/dx;
69
 
            fj0=((double)x[q][1]-origin[1])/dx;
70
 
            fi1=((double)x[p][0]-origin[0])/dx;
71
 
            fj1=((double)x[p][1]-origin[1])/dx;
72
 
         }
73
 
         j0=clamp(int(std::floor(fj0)), 0, ny-1)+1;
74
 
         j1=clamp(int(std::floor(fj1)), 0, ny-1);
75
 
         for(int j=j0; j<=j1; ++j){
76
 
            double alpha=(j-fj0)/(fj1-fj0);
77
 
            double fi=(1-alpha)*fi0 + alpha*fi1; // where the edge intersects the j'th grid line
78
 
            int i_interval=int(std::ceil(fi)); // intersection is in (i_interval-1,i_interval]
79
 
            if(i_interval<0) i_interval=0; // we enlarge the first interval to include everything to the left
80
 
            if(i_interval<nx){
81
 
               ++intersection_count(i_interval,j);
82
 
            } // we ignore intersections that are beyond the right edge of the grid
83
 
         }
84
 
      }
85
 
   }
86
 
   // and now we fill in the rest of the distances with fast sweeping
87
 
   for(unsigned int pass=0; pass<2; ++pass){
88
 
      for(int j=1; j<ny; ++j) for(int i=1; i<nx; ++i){
89
 
         Vec2d gx(i*dx+origin[0], j*dx+origin[1]);
90
 
         check_neighbour(edge, x, phi, closest_edge, gx, i, j, i-1, j-1);
91
 
         check_neighbour(edge, x, phi, closest_edge, gx, i, j, i-1, j);
92
 
         check_neighbour(edge, x, phi, closest_edge, gx, i, j, i, j-1);
93
 
      }
94
 
      for(int j=ny-2; j>=0; --j) for(int i=nx-2; i>=0; --i){
95
 
         Vec2d gx(i*dx+origin[0], j*dx+origin[1]);
96
 
         check_neighbour(edge, x, phi, closest_edge, gx, i, j, i+1, j+1);
97
 
         check_neighbour(edge, x, phi, closest_edge, gx, i, j, i+1, j);
98
 
         check_neighbour(edge, x, phi, closest_edge, gx, i, j, i, j+1);
99
 
      }
100
 
      for(int j=ny-2; j>=0; --j) for(int i=1; i<nx; ++i){
101
 
         Vec2d gx(i*dx+origin[0], j*dx+origin[1]);
102
 
         check_neighbour(edge, x, phi, closest_edge, gx, i, j, i-1, j+1);
103
 
         check_neighbour(edge, x, phi, closest_edge, gx, i, j, i-1, j);
104
 
         check_neighbour(edge, x, phi, closest_edge, gx, i, j, i, j+1);
105
 
      }
106
 
      for(int j=1; j<ny; ++j) for(int i=nx-2; i>=0; --i){
107
 
         Vec2d gx(i*dx+origin[0], j*dx+origin[1]);
108
 
         check_neighbour(edge, x, phi, closest_edge, gx, i, j, i+1, j-1);
109
 
         check_neighbour(edge, x, phi, closest_edge, gx, i, j, i+1, j);
110
 
         check_neighbour(edge, x, phi, closest_edge, gx, i, j, i, j-1);
111
 
      }
112
 
   }
113
 
   // then figure out signs (inside/outside) from intersection counts
114
 
   for(int j=0; j<ny; ++j){
115
 
      int total_count=0;
116
 
      for(int i=0; i<nx; ++i){
117
 
         total_count+=intersection_count(i,j);
118
 
         if(total_count%2==1){ // if parity of intersections so far is odd,
119
 
            phi(i,j)=-phi(i,j); // we are inside the mesh
120
 
         }
121
 
      }
122
 
   }
123
 
}
124