size()-1; ++i )
{
int j = i + ((unsigned)rand()) % (samples->size()-i);
Coord_Diff tmp = (*samples)[j];
(*samples)[i] = (*samples)[j];
(*samples)[j] = tmp;
}
}
// Allocate the matrix and vector
mtx_l.resize(p+3, m+3);
mtx_v.resize(p+3, 2);
mtx_orig_k.resize(p, m);
// Fill K (p x m, upper left of L)
for ( unsigned i=0; i& pts )
{
assert( samples.get() && "Morpher no longer owns 'samples'");
const unsigned m = mtx_orig_k.size2();
for ( std::vector::iterator ite=pts.begin(), end=pts.end();
ite != end;
++ite )
{
double x=ite->x, y=ite->y;
double dx = mtx_v(m+0, 0) + mtx_v(m+1, 0)*x + mtx_v(m+2, 0)*y;
double dy = mtx_v(m+0, 1) + mtx_v(m+1, 1)*x + mtx_v(m+2, 1)*y;
std::vector::const_iterator diff_ite = samples->begin();
Matrix_Col cv0(mtx_v,0), cv1(mtx_v,1);
Matrix_Col::const_iterator cv0_ite(cv0.begin()), cv1_ite(cv1.begin());
for ( unsigned i=0; ix - x) + SQR(diff_ite->y - y) );
dx += (*cv0_ite) * d;
dy += (*cv1_ite) * d;
}
ite->x += dx;
ite->y += dy;
}
}
/// Calculate bending energy, or if subsampling_factor
/// for constructor was < 1.0, approximate it.
double calc_bending_energy()
{
assert( samples.get() && "Morpher no longer owns 'samples'");
// bending energy = trace( W^T * A * W ),
// where A = upper left m x m block of mtx_orig_k
const unsigned m = mtx_orig_k.size2();
ublas::matrix_range mtx_w(mtx_v,
ublas::range(0, m),
ublas::range(0, 2));
ublas::matrix_range mtx_a(mtx_orig_k,
ublas::range(0, m),
ublas::range(0, m));
Matrix bm = prod( Matrix(prod(trans(mtx_w), mtx_a)), mtx_w);
assert( bm.size1() == bm.size2() && bm.size1() == 2 );
return bm(0,0) + bm(1,1);
}
/// Takes away the ownership of 'samples' from the morpher.
/// After this, calling other functions becomes illegal.
std::auto_ptr > grab_samples()
{
return samples;
}
private:
inline double base_func(double r2)
{
// same as r*r * log(r), but for r^2:
return ( r2==0 )
? 0.0 // function limit at 0
: r2 * log(r2) * 0.217147241; // = 1/(2*log(10))
}
std::auto_ptr > samples;
Matrix mtx_l, mtx_v, mtx_orig_k;
};