← Back to branch summary

~siretart/ubuntu/utopic/blender/libav10

« back to all changes in this revision

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

  • Committer: Package Import Robot
  • Author(s): Matteo F. Vescovi
  • Date: 2012-07-23 08:54:18 UTC
  • mfrom: (14.2.16 sid)
  • mto: (14.2.19 sid)
  • mto: This revision was merged to the branch mainline in revision 42.
  • Revision ID: package-import@ubuntu.com-20120723085418-9foz30v6afaf5ffs
Tags: 2.63a-2
http://bugs.debian.org/658075
* debian/: Cycles support added (Closes: #658075)
  For now, this top feature has been enabled only
  on [any-amd64 any-i386] architectures because
  of OpenImageIO failing on all others
* debian/: scripts installation path changed
  from /usr/lib to /usr/share:
  + debian/patches/: patchset re-worked for path changing
  + debian/control: "Breaks" field added on yafaray-exporter

Show diffs side-by-side

added added

removed removed

Lines of Context:
extern/eltopo/common/levelset.cpp
 
1
#if 0
 
2
#include "levelset.h"
 
3
#include "util.h"
 
4
#include "array2_utils.h"
 
5
#include "array3_utils.h"
 
6
 
 
7
float interpolate_phi(const Vec2f& point, const Array2f& grid, const Vec2f& origin, const float dx) {
 
8
   float inv_dx = 1/dx;
 
9
   Vec2f temp = (point-origin)*inv_dx;
 
10
   return interpolate_value(temp, grid);
 
11
}
 
12
 
 
13
float interpolate_phi(const Vec3f& point, const Array3f& grid, const Vec3f& origin, const float dx) {
 
14
   float inv_dx = 1/dx;
 
15
   Vec3f temp = (point-origin)*inv_dx;
 
16
   return interpolate_value(temp, grid);
 
17
}
 
18
 
 
19
float interpolate_normal(Vec2f& normal, const Vec2f& point, const Array2f& grid, const Vec2f& origin, const float dx) {
 
20
   float inv_dx = 1/dx;
 
21
   Vec2f temp = (point-origin)*inv_dx;
 
22
   float value = interpolate_gradient(normal, temp, grid);
 
23
   if(mag(normal) != 0)
 
24
      normalize(normal);
 
25
   return value;
 
26
}
 
27
 
 
28
float interpolate_normal(Vec3f& normal, const Vec3f& point, const Array3f& grid, const Vec3f& origin, const float dx) {
 
29
   float inv_dx = 1/dx;
 
30
   Vec3f temp = (point-origin)*inv_dx;
 
31
   float value = interpolate_gradient(normal, temp, grid);
 
32
   if(mag(normal) != 0)
 
33
      normalize(normal);
 
34
   return value;
 
35
}
 
36
 
 
37
void project_to_isosurface(Vec2f& point, const float target_value, const Array2f& grid, const Vec2f& origin, const float dx) {
 
38
   float tol = 0.01f*dx; //some fraction of a grid cell;
 
39
   int max_iter = 5;
 
40
   
 
41
   int iter = 0;
 
42
   Vec2f normal;
 
43
   float phi = interpolate_normal(normal, point, grid, origin, dx);
 
44
   while(fabs(phi - target_value) > tol && iter++ < max_iter) {
 
45
      point -= (phi - target_value) * normal;
 
46
      phi = interpolate_normal(normal, point, grid, origin, dx);
 
47
   }
 
48
}
 
49
 
 
50
void project_to_isosurface(Vec3f& point, const float target_value, const Array3f& grid, const Vec3f& origin, const float dx) {
 
51
   float tol = 0.01f*dx; //some fraction of a grid cell;
 
52
   int max_iter = 5;
 
53
   
 
54
   int iter = 0;
 
55
   Vec3f normal;
 
56
   float phi = interpolate_normal(normal, point, grid, origin, dx);
 
57
   while(fabs(phi - target_value) > tol && iter++ < max_iter) {
 
58
      point -= (phi - target_value) * normal;
 
59
      phi = interpolate_normal(normal, point, grid, origin, dx);
 
60
   }
 
61
}
 
62
 
 
63
 
 
64
//On a signed distance field, compute the fraction inside the (negative) isosurface
 
65
//along the 1D line segment joining phi_left and phi_right sample points.
 
66
float fraction_inside(float phi_left, float phi_right)
 
67
{
 
68
   //Note: should not generate divide by zero, because if
 
69
   //signs are different, and both values are != 0,
 
70
   //abs(left - right) is necessarily >= left, or right, alone, ie. not 0
 
71
   
 
72
   return 
 
73
      (phi_left >= 0 && phi_right >= 0)? //all empty
 
74
         0.0f : 
 
75
         (  (phi_left < 0 && phi_right < 0)? //all full
 
76
            1.0f:
 
77
            (
 
78
               (phi_left >= 0)?
 
79
               (1 - phi_left / (phi_left - phi_right)): //right inside 
 
80
               (phi_left / (phi_left - phi_right)) //left inside
 
81
            )
 
82
         );
 
83
}
 
84
 
 
85
//On a signed distance field, compute the fraction inside the (negative) isosurface
 
86
//along the 1D line segment joining phi_left and phi_right sample points.
 
87
//Except now there are two level sets, and we want the fraction that is the union
 
88
//of their interiors
 
89
float fraction_inside_either(float phi_left0, float phi_right0, float phi_left1, float phi_right1)
 
90
{
 
91
   if(phi_left0 <= 0) {
 
92
      if(phi_right0 <= 0) {
 
93
         //entirely inside solid0 [-/-][?]
 
94
         return 1;
 
95
      }
 
96
      else { //phi_right0 > 0
 
97
         if(phi_left1 <= 0) {
 
98
            if(phi_right1 <= 0) {
 
99
               //entirely inside solid1 -> [-/+][-/-]
 
100
               return 1;
 
101
            }
 
102
            else {//both left sides are inside, neither right side [-/+][-/+]
 
103
               return max( fraction_inside(phi_left0, phi_right0), 
 
104
                           fraction_inside(phi_left1, phi_right1) );
 
105
            }
 
106
         }
 
107
         else { //phi_left1 > 0 
 
108
            if(phi_right1 <= 0) { //opposite sides are interior [-/+][+/-]
 
109
               float frac0 = fraction_inside(phi_left0, phi_right0);
 
110
               float frac1 = fraction_inside(phi_left1, phi_right1);
 
111
               float total =  frac0+frac1;
 
112
               if(total <= 1)
 
113
                  return total;
 
114
               else
 
115
                  return 1;
 
116
 
 
117
            }
 
118
            else {//phi_right1 > 0
 
119
               //phi1 plays no role, both outside [-/+][+/+]
 
120
               return fraction_inside(phi_left0, phi_right0);
 
121
            }
 
122
         }
 
123
      }
 
124
   }
 
125
   else {
 
126
      if(phi_right0 <= 0) {
 
127
         if(phi_left1 <= 0) {
 
128
            if(phi_right1 <= 0) {
 
129
               //entirely inside solid1[+/-][-/-]
 
130
               return 1;
 
131
            }
 
132
            else {
 
133
               //coming in from opposing sides [+/-][-/+]
 
134
               float frac0 = fraction_inside(phi_left0, phi_right0);
 
135
               float frac1 = fraction_inside(phi_left1, phi_right1);
 
136
               float total =  frac0+frac1;
 
137
               if(total <= 1)
 
138
                  return total;
 
139
               else
 
140
                  return 1;
 
141
            }
 
142
         }
 
143
         else { //phi_left1 > 0 
 
144
            if(phi_right1 <= 0) {
 
145
               //coming from the same side, take the larger one [+/-][+/-]
 
146
               return max( fraction_inside(phi_left0, phi_right0), 
 
147
                           fraction_inside(phi_left1, phi_right1) );
 
148
            }
 
149
            else { //phi_right > 0
 
150
               //Only have to worry about phi_0 [+/-][+/+]
 
151
               return fraction_inside(phi_left0, phi_right0);
 
152
            }
 
153
             
 
154
         }
 
155
      }
 
156
      else {
 
157
         //Only have to worry about phi_1 [+/+][?]
 
158
         return fraction_inside(phi_left1, phi_right1);
 
159
      }
 
160
   }
 
161
}
 
162
 
 
163
void compute_volume_fractions(const Array2f& levelset, const Vec2f& ls_origin, float ls_dx, Array2f& volumes, const Vec2f& v_origin, float v_dx, int subdivisions) {
 
164
   
 
165
   float sub_dx = v_dx / (float)(subdivisions+1);
 
166
   
 
167
   for(int j = 0; j < volumes.nj; ++j) for(int i = 0; i < volumes.ni; ++i) {
 
168
      //centre of the volume cells
 
169
      Vec2f bottom_left = v_origin + Vec2f(i*v_dx, j*v_dx);
 
170
      int inside_samples = 0;
 
171
      for(int subj = 0; subj < subdivisions+1; ++subj) for(int subi = 0; subi < subdivisions+1; ++subi) {
 
172
         Vec2f point = bottom_left + Vec2f( (subi+0.5f)*sub_dx, (subj+0.5f)*sub_dx);
 
173
         float data = interpolate_phi(point, levelset, ls_origin, ls_dx);
 
174
         inside_samples += (data < 0)?1:0;
 
175
      }
 
176
      volumes(i,j) = (float)inside_samples / (float)sqr(subdivisions+1);
 
177
   }
 
178
}
 
179
 
 
180
void compute_volume_fractions(const Array3f& levelset, const Vec3f& ls_origin, float ls_dx, Array3f& volumes, const Vec3f& v_origin, float v_dx, int subdivisions) {
 
181
   
 
182
   float sub_dx = v_dx / (float)(subdivisions+1);
 
183
   
 
184
   for(int k = 0; k < volumes.nk; ++k) for(int j = 0; j < volumes.nj; ++j) for(int i = 0; i < volumes.ni; ++i) {
 
185
      //centre of the volume cells
 
186
      Vec3f bottom_left = v_origin + Vec3f((i-0.5f)*v_dx, (j-0.5f)*v_dx, (k-0.5f)*v_dx);
 
187
      int inside_samples = 0;
 
188
      
 
189
      //Speedup! Test the centre point, and if it's more than a grid cell away from the interface, we can assume 
 
190
      //the cell is either totally full or totally empty
 
191
      float estimate = interpolate_phi(bottom_left + 0.5f*v_dx*Vec3f(1,1,1), levelset, ls_origin, ls_dx);
 
192
      if(estimate > v_dx) {
 
193
         volumes(i,j,k) = 0;
 
194
         continue;
 
195
      }
 
196
      else if(estimate < -v_dx) {
 
197
         volumes(i,j,k) = 1;
 
198
         continue;
 
199
      }
 
200
 
 
201
      for(int subk = 0; subk < subdivisions+1; ++subk) for(int subj = 0; subj < subdivisions+1; ++subj) for(int subi = 0; subi < subdivisions+1; ++subi) {
 
202
         Vec3f point = bottom_left + Vec3f( (subi+0.5f)*sub_dx, (subj+0.5f)*sub_dx, (subk+0.5f)*sub_dx);
 
203
         float data = interpolate_phi(point, levelset, ls_origin, ls_dx);
 
204
         inside_samples += (data < 0)?1:0;
 
205
      }
 
206
      volumes(i,j,k) = (float)inside_samples / (float)cube(subdivisions+1);
 
207
   }
 
208
}
 
209
 
 
210
#endif