| 688 |
khor |
1 |
//
|
|
|
2 |
// Implicit.cpp
|
|
|
3 |
// PointReconstruction
|
|
|
4 |
//
|
|
|
5 |
// Created by J. Andreas Bærentzen on 16/03/13.
|
|
|
6 |
// Copyright (c) 2013 J. Andreas Bærentzen. All rights reserved.
|
|
|
7 |
//
|
|
|
8 |
|
|
|
9 |
#include "../CGLA/CGLA.h"
|
|
|
10 |
#include "../Geometry/Neighbours.h"
|
|
|
11 |
#include "../Geometry/GridAlgorithm.h"
|
|
|
12 |
#include "XForm.h"
|
|
|
13 |
#include "Implicit.h"
|
|
|
14 |
|
|
|
15 |
using namespace std;
|
|
|
16 |
using namespace CGLA;
|
|
|
17 |
|
|
|
18 |
namespace Geometry
|
|
|
19 |
{
|
|
|
20 |
void Implicit::push_to_surface(Vec3d& p, double tau, double max_dist) const
|
|
|
21 |
{
|
|
|
22 |
Vec3d g = grad(p);
|
|
|
23 |
|
|
|
24 |
double sl = sqr_length(g);
|
|
|
25 |
double d = (eval(p)-tau)/sl;
|
|
|
26 |
Vec3d disp = g*d;
|
|
|
27 |
double disp_len = length(disp);
|
|
|
28 |
double clamped_disp_len = min(max_dist, disp_len);
|
|
|
29 |
p = p - clamped_disp_len*disp/disp_len;
|
|
|
30 |
}
|
|
|
31 |
|
|
|
32 |
XForm grid_sample(const Implicit& imp, const CGLA::Vec3d& llf, const CGLA::Vec3d& urt,
|
|
|
33 |
Geometry::RGrid<float>& grid)
|
|
|
34 |
{
|
|
|
35 |
XForm xform(llf,urt, grid.get_dims(), 0.0);
|
|
|
36 |
for(int i=0;i<xform.get_dims()[0];++i)
|
|
|
37 |
for(int j=0;j<xform.get_dims()[1];++j)
|
|
|
38 |
for(int k=0;k<xform.get_dims()[2];++k)
|
|
|
39 |
{
|
|
|
40 |
Vec3d p = xform.inverse(Vec3d(i,j,k));
|
|
|
41 |
float f = imp.eval(p);
|
|
|
42 |
grid[Vec3i(i,j,k)] = f;
|
|
|
43 |
}
|
|
|
44 |
return xform;
|
|
|
45 |
}
|
|
|
46 |
|
|
|
47 |
|
|
|
48 |
float interpolate(const RGrid<float>& grid, const CGLA::Vec3d& _v)
|
|
|
49 |
{
|
|
|
50 |
Vec3d v = v_min(Vec3d(grid.get_dims()-Vec3i(1)), v_max(_v,Vec3d(0)));
|
|
|
51 |
|
|
|
52 |
Vec3i c0i(v);
|
|
|
53 |
|
|
|
54 |
const float alpha = v[0] - float(c0i[0]);
|
|
|
55 |
const float beta = v[1] - float(c0i[1]);
|
|
|
56 |
const float gamm = v[2] - float(c0i[2]);
|
|
|
57 |
float m_alpha = 1.0 - alpha;
|
|
|
58 |
float m_beta = 1.0 - beta;
|
|
|
59 |
float m_gamm = 1.0 - gamm;
|
|
|
60 |
float weights[8];
|
|
|
61 |
weights[0] = (m_alpha*m_beta*m_gamm);
|
|
|
62 |
weights[1] = (alpha*m_beta*m_gamm);
|
|
|
63 |
weights[2] = (m_alpha*beta*m_gamm);
|
|
|
64 |
weights[3] = (alpha*beta*m_gamm);
|
|
|
65 |
weights[4] = (m_alpha*m_beta*gamm);
|
|
|
66 |
weights[5] = (alpha*m_beta*gamm);
|
|
|
67 |
weights[6] = (m_alpha*beta*gamm);
|
|
|
68 |
weights[7] = (alpha*beta*gamm);
|
|
|
69 |
|
|
|
70 |
float f = 0;
|
|
|
71 |
for(int i=0;i<8;++i)
|
|
|
72 |
if(weights[i]>1e-10)
|
|
|
73 |
f += weights[i]*grid[c0i+Geometry::CubeCorners8i[i]];
|
|
|
74 |
|
|
|
75 |
return f;
|
|
|
76 |
}
|
|
|
77 |
|
|
|
78 |
|
|
|
79 |
double VolumetricImplicit::eval(const CGLA::Vec3d& p) const
|
|
|
80 |
{
|
|
|
81 |
return interpolate(grid,xform.apply(p));
|
|
|
82 |
}
|
|
|
83 |
|
|
|
84 |
CGLA::Vec3d VolumetricImplicit::grad(const CGLA::Vec3d& _p) const
|
|
|
85 |
{
|
|
|
86 |
Vec3d p = xform.apply(_p);
|
|
|
87 |
Vec3d g(interpolate(grid,p+Vec3d(1,0,0))-interpolate(grid,p-Vec3d(1,0,0)),
|
|
|
88 |
interpolate(grid,p+Vec3d(0,1,0))-interpolate(grid,p-Vec3d(0,1,0)),
|
|
|
89 |
interpolate(grid,p+Vec3d(0,0,1))-interpolate(grid,p-Vec3d(0,0,1)));
|
|
|
90 |
g *= 0.5 / xform.inv_scale();
|
|
|
91 |
return g;
|
|
|
92 |
}
|
|
|
93 |
|
|
|
94 |
}
|