Subversion Repositories gelsvn

Rev

Rev 667 | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 667 Rev 677
Line 53... Line 53...
53 
                swap(vec,new_vec);
 53 
                swap(vec,new_vec);
54 
            }
54 
            }
55 
        }
 55 
        }
56 
    }
 56 
    }
57 
 
 57 
 
58 
    double voronoi_area(const Manifold& m, VertexID v)
 58 
    double mixed_area(const Manifold& m, VertexID v)
59 
    {
 59 
    {
60 
        double area_mixed = 0;
 60 
        double area_mixed = 0;
61 
        //For each triangle T from the 1-ring neighborhood of x
61 
        //For each triangle T from the 1-ring neighborhood of x
62 
        for(Walker w = m.walker(v); !w.full_circle(); w = w.circulate_vertex_cw()){
 62 
        for(Walker w = m.walker(v); !w.full_circle(); w = w.circulate_vertex_cw()){
63 
            double f_area = area(m, w.face());
 63 
            double f_area = area(m, w.face());
Line 78... Line 78...
78 
                    (1.0/tan(a2)) * sqr_length(v1-v0));
 78 
                    (1.0/tan(a2)) * sqr_length(v1-v0));
79 
            }
 79 
            }
80 
            else // Voronoi inappropriate
80 
            else // Voronoi inappropriate
81 
            {
 81 
            {
82 
                // Add either area(f)/4 or area(f)/2
82 
                // Add either area(f)/4 or area(f)/2
83 
                if(a0>M_PI/2.0)// the angle of f at x is obtuse
- 
 
84 
                    area_mixed += f_area/2;
 - 
 
85 
                else
 - 
 
86 
                    area_mixed += f_area/4;
 83 
                area_mixed += f_area/3;
87 
            }
 84 
            }
88 
        }
 85 
        }
89 
        return area_mixed;
 86 
        return area_mixed;
90 
    }
 87 
    }
91 
 
 88 
 
Line 115... Line 112...
115 
            double d_left = dot(cond_normalize(nbr-left),cond_normalize(vertex-left));
 112 
            double d_left = dot(cond_normalize(nbr-left),cond_normalize(vertex-left));
116 
            double d_right = dot(cond_normalize(nbr-right),cond_normalize(vertex-right));
 113 
            double d_right = dot(cond_normalize(nbr-right),cond_normalize(vertex-right));
117 
            double a_left  = acos(min(1.0, max(-1.0, d_left)));
 114 
            double a_left  = acos(min(1.0, max(-1.0, d_left)));
118 
            double a_right = acos(min(1.0, max(-1.0, d_right)));
 115 
            double a_right = acos(min(1.0, max(-1.0, d_right)));
119 
 
 116 
 
120 
            double w = 1.0/(1e-300+tan(a_left));
 117 
//            double w = 1.0/(1e-300+tan(a_left));
121 
            w += 1.0/(1e-300+tan(a_right));
 118 
//            w += 1.0/(1e-300+tan(a_right));
122 
//            double w = sin(a_left + a_right) / (1e-300 + sin(a_left)*sin(a_right));
119 
            double w = sin(a_left + a_right) / (1e-300 + sin(a_left)*sin(a_right));
123 
            curv_normal += w * (nbr-vertex);
 120 
            curv_normal += w * (nbr-vertex);
124 
            w_sum += w;
 121 
            w_sum += w;
125 
        }
 122 
        }
126 
 
 123 
 
127 
    }
 124 
    }
Line 130... Line 127...
130 
    {
 127 
    {
131 
        Vec3d curv_normal;
 128 
        Vec3d curv_normal;
132 
        double w_sum;
 129 
        double w_sum;
133 
        unnormalized_mean_curvature_normal(m, v, curv_normal, w_sum);
 130 
        unnormalized_mean_curvature_normal(m, v, curv_normal, w_sum);
134 
 
 131 
 
135 
        return curv_normal / (4*voronoi_area(m, v));
 132 
        return curv_normal / (4*mixed_area(m, v));
136 
    }
 133 
    }
137 
 
 134 
 
138 
    double sum_curvatures(const Manifold& m, VertexAttributeVector<double>& curvature)
 135 
    double sum_curvatures(const Manifold& m, VertexAttributeVector<double>& curvature)
139 
    {
 136 
    {
140 
        double sum = 0;
 137 
        double sum = 0;
141 
        for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){
 138 
        for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){
142 
            if(boundary(m, *v))
 139 
            if(boundary(m, *v))
143 
                continue;
140 
                continue;
144 
            sum += curvature[*v] * voronoi_area(m, *v);
 141 
            sum += curvature[*v] * mixed_area(m, *v);
145 
        }
 142 
        }
146 
        return sum;
 143 
        return sum;
147 
    }
 144 
    }
148 
 
 145 
 
149 
 
 146 
 
Line 164... Line 161...
164 
        {
 161 
        {
165 
            double dot_prod =
162 
            double dot_prod =
166 
                std::max(-1.0, std::min(1.0, dot(edges[i],edges[(i+1)%N])));
 163 
                std::max(-1.0, std::min(1.0, dot(edges[i],edges[(i+1)%N])));
167 
            angle_sum += acos(dot_prod);
 164 
            angle_sum += acos(dot_prod);
168 
        }
 165 
        }
169 
        return (2*M_PI - angle_sum)/voronoi_area(m, v);
 166 
        return (2*M_PI - angle_sum)/mixed_area(m, v);
170 
 
 167 
 
171 
    }
 168 
    }
172 
 
 169 
 
173 
    Mat3x3d curvature_tensor(const Manifold& m, HalfEdgeID h)
 170 
    Mat3x3d curvature_tensor(const Manifold& m, HalfEdgeID h)
174 
    {
 171 
    {
Line 201... Line 198...
201 
            return curv_tensor;
 198 
            return curv_tensor;
202 
 
 199 
 
203 
        for(Walker w = m.walker(v); !w.full_circle(); w = w.circulate_vertex_cw())
 200 
        for(Walker w = m.walker(v); !w.full_circle(); w = w.circulate_vertex_cw())
204 
            curv_tensor += 0.5*curvature_tensor(m, w.halfedge());
 201 
            curv_tensor += 0.5*curvature_tensor(m, w.halfedge());
205 
 
 202 
 
206 
        curv_tensor /= voronoi_area(m, v);
 203 
        curv_tensor /= mixed_area(m, v);
207 
 
 204 
 
208 
        return curv_tensor;
 205 
        return curv_tensor;
209 
    }
 206 
    }
210 
 
 207 
 
211 
 
 208 
 
Line 263... Line 260...
263 
        double tmp_area;
 260 
        double tmp_area;
264 
 
 261 
 
265 
        for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){
 262 
        for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){
266 
            if(boundary(m, *v))
 263 
            if(boundary(m, *v))
267 
                continue;
 264 
                continue;
268 
            double a = voronoi_area(m, *v);
 265 
            double a = mixed_area(m, *v);
269 
            tmp_curvature_tensors[*v] = curvature_tensors[*v] * a;
 266 
            tmp_curvature_tensors[*v] = curvature_tensors[*v] * a;
270 
            tmp_area = a;
 267 
            tmp_area = a;
271 
            for(Walker w = m.walker(*v); !w.full_circle(); w = w.circulate_vertex_cw()){
 268 
            for(Walker w = m.walker(*v); !w.full_circle(); w = w.circulate_vertex_cw()){
272 
                if(!boundary(m, w.vertex())){
 269 
                if(!boundary(m, w.vertex())){
273 
                    double a = voronoi_area(m, w.vertex());
 270 
                    double a = mixed_area(m, w.vertex());
274 
                    tmp_curvature_tensors[*v] += curvature_tensors[w.vertex()]*a;
 271 
                    tmp_curvature_tensors[*v] += curvature_tensors[w.vertex()]*a;
275 
                    tmp_area += a;
 272 
                    tmp_area += a;
276 
                }
 273 
                }
277 
                tmp_curvature_tensors[*v] /= tmp_area;
 274 
                tmp_curvature_tensors[*v] /= tmp_area;
278 
            }
 275 
            }