73 template <
typename VectorT >
76 center_of_gravity(
const std::vector< VectorT >& _points ) {
79 for (uint i = 0 ; i < _points.size() ; ++i ) {
80 cog = cog + _points[i];
83 cog = cog / (
typename VectorT::value_type )_points.size();
89 template <
typename VectorT ,
typename QuaternionT >
91 icp(
const std::vector< VectorT >& _points1 ,
const std::vector< VectorT >& _points2 ,
VectorT& _cog1 ,
VectorT& _cog2 ,
double& _scale1 ,
double& _scale2 , QuaternionT& _rotation )
94 _cog1 = center_of_gravity(_points1);
95 _cog2 = center_of_gravity(_points2);
101 for ( uint i = 0 ; i < _points1.size() ; ++i ) {
102 diff = _points1[i] - _cog1;
103 _scale1 = std::max( _scale1 , sqrt( diff.norm() ) );
104 diff = _points2[i] - _cog2;
105 _scale2 = std::max( _scale2 , sqrt( diff.norm() ) );
120 for ( uint i = 0 ; i < _points1.size() ; ++i ) {
121 Sxx += ( _points1[i][0] - _cog1[0] ) * ( _points2[i][0] - _cog2[0] );
122 Sxy += ( _points1[i][0] - _cog1[0] ) * ( _points2[i][1] - _cog2[1] );
123 Sxz += ( _points1[i][0] - _cog1[0] ) * ( _points2[i][2] - _cog2[2] );
125 Syx += ( _points1[i][1] - _cog1[1] ) * ( _points2[i][0] - _cog2[0] );
126 Syy += ( _points1[i][1] - _cog1[1] ) * ( _points2[i][1] - _cog2[1] );
127 Syz += ( _points1[i][1] - _cog1[1] ) * ( _points2[i][2] - _cog2[2] );
129 Szx += ( _points1[i][2] - _cog1[2] ) * ( _points2[i][0] - _cog2[0] );
130 Szy += ( _points1[i][2] - _cog1[2] ) * ( _points2[i][1] - _cog2[1] );
131 Szz += ( _points1[i][2] - _cog1[2] ) * ( _points2[i][2] - _cog2[2] );
134 const double scale = _scale1 * _scale2;
136 Sxx = Sxx * 1.0 / scale;
137 Sxy = Sxy * 1.0 / scale;
138 Sxz = Sxz * 1.0 / scale;
139 Syx = Syx * 1.0 / scale;
140 Syy = Syy * 1.0 / scale;
141 Syz = Syz * 1.0 / scale;
142 Szx = Szx * 1.0 / scale;
143 Szy = Szy * 1.0 / scale;
144 Szz = Szz * 1.0 / scale;
146 gmm::dense_matrix<double> N; gmm::resize(N,4,4); gmm::clear(N);
147 N(0,0) = Sxx + Syy + Szz;
153 N(1,1) = Sxx - Syy - Szz;
159 N(2,2) = -Sxx + Syy - Szz;
165 N(3,3) = -Sxx - Syy + Szz;
167 std::vector< double > eigval; eigval.resize(4);
168 gmm::dense_matrix< double > eigvect; gmm::resize(eigvect, 4,4); gmm::clear(eigvect);
169 gmm::symmetric_qr_algorithm(N, eigval, eigvect);
172 double max = -FLT_MAX;
173 for (
int i = 0; i < 4; i++)
180 _rotation[0] = eigvect(0,gr);
181 _rotation[1] = eigvect(1,gr);
182 _rotation[2] = eigvect(2,gr);
183 _rotation[3] = eigvect(3,gr);
189 template <
typename VectorT ,
typename QuaternionT >
191 icp(
const std::vector< VectorT >& _points1 ,
const std::vector< VectorT >& _points2 ,
VectorT& _cog1 ,
VectorT& _cog2 , QuaternionT& _rotation )
194 _cog1 = center_of_gravity(_points1);
195 _cog2 = center_of_gravity(_points2);
207 for ( uint i = 0 ; i < _points1.size() ; ++i ) {
208 Sxx += ( _points1[i][0] - _cog1[0] ) * ( _points2[i][0] - _cog2[0] );
209 Sxy += ( _points1[i][0] - _cog1[0] ) * ( _points2[i][1] - _cog2[1] );
210 Sxz += ( _points1[i][0] - _cog1[0] ) * ( _points2[i][2] - _cog2[2] );
212 Syx += ( _points1[i][1] - _cog1[1] ) * ( _points2[i][0] - _cog2[0] );
213 Syy += ( _points1[i][1] - _cog1[1] ) * ( _points2[i][1] - _cog2[1] );
214 Syz += ( _points1[i][1] - _cog1[1] ) * ( _points2[i][2] - _cog2[2] );
216 Szx += ( _points1[i][2] - _cog1[2] ) * ( _points2[i][0] - _cog2[0] );
217 Szy += ( _points1[i][2] - _cog1[2] ) * ( _points2[i][1] - _cog2[1] );
218 Szz += ( _points1[i][2] - _cog1[2] ) * ( _points2[i][2] - _cog2[2] );
221 gmm::dense_matrix<double> N; gmm::resize(N,4,4); gmm::clear(N);
222 N(0,0) = Sxx + Syy + Szz;
228 N(1,1) = Sxx - Syy - Szz;
234 N(2,2) = -Sxx + Syy - Szz;
240 N(3,3) = -Sxx - Syy + Szz;
242 std::vector< double > eigval; eigval.resize(4);
243 gmm::dense_matrix< double > eigvect; gmm::resize(eigvect, 4,4); gmm::clear(eigvect);
244 gmm::symmetric_qr_algorithm(N, eigval, eigvect);
247 double max = -FLT_MAX;
248 for (
int i = 0; i < 4; i++)
255 _rotation[0] = eigvect(0,gr);
256 _rotation[1] = eigvect(1,gr);
257 _rotation[2] = eigvect(2,gr);
258 _rotation[3] = eigvect(3,gr);
Classes for ICP Algorithm.
void icp(const std::vector< VectorT > &_points1, const std::vector< VectorT > &_points2, VectorT &_cog1, VectorT &_cog2, double &_scale1, double &_scale2, QuaternionT &_rotation)
Compute rigid transformation from first point set to second point set.