Point Cloud Library (PCL) 1.14.0
Loading...
Searching...
No Matches
projection_matrix.hpp
1/*
2 * Software License Agreement (BSD License)
3 *
4 * Point Cloud Library (PCL) - www.pointclouds.org
5 * Copyright (c) 2012-, Open Perception, Inc.
6 *
7 * All rights reserved.
8 *
9 * Redistribution and use in source and binary forms, with or without
10 * modification, are permitted provided that the following conditions
11 * are met:
12 *
13 * * Redistributions of source code must retain the above copyright
14 * notice, this list of conditions and the following disclaimer.
15 * * Redistributions in binary form must reproduce the above
16 * copyright notice, this list of conditions and the following
17 * disclaimer in the documentation and/or other materials provided
18 * with the distribution.
19 * * Neither the name of the copyright holder(s) nor the names of its
20 * contributors may be used to endorse or promote products derived
21 * from this software without specific prior written permission.
22 *
23 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
24 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
25 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
26 * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
27 * COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
28 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
29 * BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
30 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
31 * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
32 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
33 * ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
34 * POSSIBILITY OF SUCH DAMAGE.
35 *
36 */
37
38#pragma once
39
40#include <pcl/common/projection_matrix.h>
41#include <pcl/console/print.h> // for PCL_ERROR
42#include <pcl/cloud_iterator.h>
43
44#include <Eigen/Eigenvalues> // for SelfAdjointEigenSolver
45
46namespace pcl
47{
48
49namespace common
50{
51
52namespace internal
53{
54
55template <typename MatrixT> void
56makeSymmetric (MatrixT& matrix, bool use_upper_triangular = true)
57{
58 if (use_upper_triangular && (MatrixT::Flags & Eigen::RowMajorBit))
59 {
60 matrix.coeffRef (4) = matrix.coeff (1);
61 matrix.coeffRef (8) = matrix.coeff (2);
62 matrix.coeffRef (9) = matrix.coeff (6);
63 matrix.coeffRef (12) = matrix.coeff (3);
64 matrix.coeffRef (13) = matrix.coeff (7);
65 matrix.coeffRef (14) = matrix.coeff (11);
66 }
67 else
68 {
69 matrix.coeffRef (1) = matrix.coeff (4);
70 matrix.coeffRef (2) = matrix.coeff (8);
71 matrix.coeffRef (6) = matrix.coeff (9);
72 matrix.coeffRef (3) = matrix.coeff (12);
73 matrix.coeffRef (7) = matrix.coeff (13);
74 matrix.coeffRef (11) = matrix.coeff (14);
75 }
76}
77
78} // namespace internal
79} // namespace common
80
81
82template <typename PointT> double
85 Eigen::Matrix<float, 3, 4, Eigen::RowMajor>& projection_matrix,
86 const Indices& indices)
87{
88 // internally we calculate with double but store the result into float matrices.
89 using Scalar = double;
90 projection_matrix.setZero ();
91 if (cloud->height == 1 || cloud->width == 1)
92 {
93 PCL_ERROR ("[pcl::estimateProjectionMatrix] Input dataset is not organized!\n");
94 return (-1.0);
95 }
96
97 Eigen::Matrix<Scalar, 4, 4, Eigen::RowMajor> A = Eigen::Matrix<Scalar, 4, 4, Eigen::RowMajor>::Zero ();
98 Eigen::Matrix<Scalar, 4, 4, Eigen::RowMajor> B = Eigen::Matrix<Scalar, 4, 4, Eigen::RowMajor>::Zero ();
99 Eigen::Matrix<Scalar, 4, 4, Eigen::RowMajor> C = Eigen::Matrix<Scalar, 4, 4, Eigen::RowMajor>::Zero ();
100 Eigen::Matrix<Scalar, 4, 4, Eigen::RowMajor> D = Eigen::Matrix<Scalar, 4, 4, Eigen::RowMajor>::Zero ();
101
102 pcl::ConstCloudIterator <PointT> pointIt (*cloud, indices);
103
104 while (pointIt)
105 {
106 unsigned yIdx = pointIt.getCurrentPointIndex () / cloud->width;
107 unsigned xIdx = pointIt.getCurrentPointIndex () % cloud->width;
108
109 const PointT& point = *pointIt;
110 if (std::isfinite (point.x))
111 {
112 Scalar xx = point.x * point.x;
113 Scalar xy = point.x * point.y;
114 Scalar xz = point.x * point.z;
115 Scalar yy = point.y * point.y;
116 Scalar yz = point.y * point.z;
117 Scalar zz = point.z * point.z;
118 Scalar xx_yy = xIdx * xIdx + yIdx * yIdx;
119
120 A.coeffRef (0) += xx;
121 A.coeffRef (1) += xy;
122 A.coeffRef (2) += xz;
123 A.coeffRef (3) += point.x;
124
125 A.coeffRef (5) += yy;
126 A.coeffRef (6) += yz;
127 A.coeffRef (7) += point.y;
128
129 A.coeffRef (10) += zz;
130 A.coeffRef (11) += point.z;
131 A.coeffRef (15) += 1.0;
132
133 B.coeffRef (0) -= xx * xIdx;
134 B.coeffRef (1) -= xy * xIdx;
135 B.coeffRef (2) -= xz * xIdx;
136 B.coeffRef (3) -= point.x * static_cast<double>(xIdx);
137
138 B.coeffRef (5) -= yy * xIdx;
139 B.coeffRef (6) -= yz * xIdx;
140 B.coeffRef (7) -= point.y * static_cast<double>(xIdx);
141
142 B.coeffRef (10) -= zz * xIdx;
143 B.coeffRef (11) -= point.z * static_cast<double>(xIdx);
144
145 B.coeffRef (15) -= xIdx;
146
147 C.coeffRef (0) -= xx * yIdx;
148 C.coeffRef (1) -= xy * yIdx;
149 C.coeffRef (2) -= xz * yIdx;
150 C.coeffRef (3) -= point.x * static_cast<double>(yIdx);
151
152 C.coeffRef (5) -= yy * yIdx;
153 C.coeffRef (6) -= yz * yIdx;
154 C.coeffRef (7) -= point.y * static_cast<double>(yIdx);
155
156 C.coeffRef (10) -= zz * yIdx;
157 C.coeffRef (11) -= point.z * static_cast<double>(yIdx);
158
159 C.coeffRef (15) -= yIdx;
160
161 D.coeffRef (0) += xx * xx_yy;
162 D.coeffRef (1) += xy * xx_yy;
163 D.coeffRef (2) += xz * xx_yy;
164 D.coeffRef (3) += point.x * xx_yy;
165
166 D.coeffRef (5) += yy * xx_yy;
167 D.coeffRef (6) += yz * xx_yy;
168 D.coeffRef (7) += point.y * xx_yy;
169
170 D.coeffRef (10) += zz * xx_yy;
171 D.coeffRef (11) += point.z * xx_yy;
172
173 D.coeffRef (15) += xx_yy;
174 }
175
176 ++pointIt;
177 } // while
178
183
184 Eigen::Matrix<Scalar, 12, 12, Eigen::RowMajor> X = Eigen::Matrix<Scalar, 12, 12, Eigen::RowMajor>::Zero ();
185 X.topLeftCorner<4,4> ().matrix () = A;
186 X.block<4,4> (0, 8).matrix () = B;
187 X.block<4,4> (8, 0).matrix () = B;
188 X.block<4,4> (4, 4).matrix () = A;
189 X.block<4,4> (4, 8).matrix () = C;
190 X.block<4,4> (8, 4).matrix () = C;
191 X.block<4,4> (8, 8).matrix () = D;
192
193 Eigen::SelfAdjointEigenSolver<Eigen::Matrix<Scalar, 12, 12, Eigen::RowMajor> > ei_symm (X);
194 Eigen::Matrix<Scalar, 12, 12, Eigen::RowMajor> eigen_vectors = ei_symm.eigenvectors ();
195
196 // check whether the residual MSE is low. If its high, the cloud was not captured from a projective device.
197 Eigen::Matrix<Scalar, 1, 1> residual_sqr = eigen_vectors.col (0).transpose () * X * eigen_vectors.col (0);
198
199 double residual = residual_sqr.coeff (0);
200
201 projection_matrix.coeffRef (0) = static_cast <float> (eigen_vectors.coeff (0));
202 projection_matrix.coeffRef (1) = static_cast <float> (eigen_vectors.coeff (12));
203 projection_matrix.coeffRef (2) = static_cast <float> (eigen_vectors.coeff (24));
204 projection_matrix.coeffRef (3) = static_cast <float> (eigen_vectors.coeff (36));
205 projection_matrix.coeffRef (4) = static_cast <float> (eigen_vectors.coeff (48));
206 projection_matrix.coeffRef (5) = static_cast <float> (eigen_vectors.coeff (60));
207 projection_matrix.coeffRef (6) = static_cast <float> (eigen_vectors.coeff (72));
208 projection_matrix.coeffRef (7) = static_cast <float> (eigen_vectors.coeff (84));
209 projection_matrix.coeffRef (8) = static_cast <float> (eigen_vectors.coeff (96));
210 projection_matrix.coeffRef (9) = static_cast <float> (eigen_vectors.coeff (108));
211 projection_matrix.coeffRef (10) = static_cast <float> (eigen_vectors.coeff (120));
212 projection_matrix.coeffRef (11) = static_cast <float> (eigen_vectors.coeff (132));
213
214 if (projection_matrix.coeff (0) < 0)
215 projection_matrix *= -1.0;
216
217 return (residual);
218}
219
220} // namespace pcl
221
unsigned getCurrentPointIndex() const
std::uint32_t width
The point cloud width (if organized as an image-structure).
std::uint32_t height
The point cloud height (if organized as an image-structure).
shared_ptr< const PointCloud< PointT > > ConstPtr
@ B
Definition norms.h:54
void makeSymmetric(MatrixT &matrix, bool use_upper_triangular=true)
double estimateProjectionMatrix(typename pcl::PointCloud< PointT >::ConstPtr cloud, Eigen::Matrix< float, 3, 4, Eigen::RowMajor > &projection_matrix, const Indices &indices)
Estimates the projection matrix P = K * (R|-R*t) from organized point clouds, with K = [[fx,...
IndicesAllocator<> Indices
Type used for indices in PCL.
Definition types.h:133
A point structure representing Euclidean xyz coordinates, and the RGB color.