22
33#include < console_bridge/console.h>
44
5- MAPGMM::MAPGMM (int n_components, const std::vector<geo::Vec3>& points, const GMMParams & params) : K_(n_components), inlier_component_(0 ), params_(params) {
5+ MAPGMM::MAPGMM (int n_components, const std::vector<geo::Vec3>& points, const GMMParams & params) : K_(n_components), inlier_component_(0 ), params_(params)
6+ {
67 // Initialize all containers to proper sizes
78 means_.resize (K_);
89 covs_.resize (K_);
@@ -16,7 +17,8 @@ MAPGMM::MAPGMM(int n_components, const std::vector<geo::Vec3>& points, const GMM
1617 nu0_.resize (K_);
1718
1819 // Initialize with identity matrices to avoid uninitialized memory
19- for (int k = 0 ; k < K_; k++) {
20+ for (int k = 0 ; k < K_; k++)
21+ {
2022 means_[k] = Eigen::Vector3d::Zero ();
2123 covs_[k] = Eigen::Matrix3d::Identity ();
2224 mu0_[k] = Eigen::Vector3d::Zero ();
@@ -27,13 +29,15 @@ MAPGMM::MAPGMM(int n_components, const std::vector<geo::Vec3>& points, const GMM
2729 setupPriors (points);
2830}
2931
30- void MAPGMM::fit (const std::vector<geo::Vec3>& points, const geo::Pose3D& sensor_pose) {
32+ void MAPGMM::fit (const std::vector<geo::Vec3>& points, const geo::Pose3D& sensor_pose)
33+ {
3134 // Convert points to Eigen matrix
3235 int N = points.size ();
3336 Eigen::MatrixXd data (N, 3 );
3437
3538 // Transform points to sensor frame
36- for (int i = 0 ; i < N; i++) {
39+ for (int i = 0 ; i < N; i++)
40+ {
3741 geo::Vec3 p_map = sensor_pose * points[i];
3842 data (i, 0 ) = p_map.x ;
3943 data (i, 1 ) = p_map.y ;
@@ -49,7 +53,8 @@ void MAPGMM::fit(const std::vector<geo::Vec3>& points, const geo::Pose3D& sensor
4953 covs_.resize (K_);
5054
5155 // Simple initialization: random points as centers for means and small covariances
52- for (int k = 0 ; k < K_; k++) {
56+ for (int k = 0 ; k < K_; k++)
57+ {
5358 int idx = rand () % N;
5459 means_[k] = data.row (idx).transpose ();
5560 covs_[k] = Eigen::Matrix3d::Identity () * 0.01 ; // Covariance here is not the prior but the initial guess of the cluster.
@@ -66,13 +71,15 @@ void MAPGMM::fit(const std::vector<geo::Vec3>& points, const geo::Pose3D& sensor
6671
6772 resp_ = Eigen::MatrixXd::Zero (N, K_);
6873
69- for (int iter = 0 ; iter < max_iter; iter++) {
74+ for (int iter = 0 ; iter < max_iter; iter++)
75+ {
7076 // E-step: Calculate responsibilities
7177 double log_likelihood = eStep (data, resp_);
7278
7379 // Check convergence
7480 double change = std::abs ((log_likelihood - prev_log_likelihood) / (std::abs (prev_log_likelihood) + 1e-10 ));
75- if (iter > 0 && change < likelihood_change) {
81+ if (iter > 0 && change < likelihood_change)
82+ {
7683 CONSOLE_BRIDGE_logInform (" MAP-GMM converged after %d iterations" , iter);
7784 break ;
7885 }
@@ -84,7 +91,8 @@ void MAPGMM::fit(const std::vector<geo::Vec3>& points, const geo::Pose3D& sensor
8491
8592 // Assign labels based on highest responsibility
8693 labels_.resize (N);
87- for (int i = 0 ; i < N; i++) {
94+ for (int i = 0 ; i < N; i++)
95+ {
8896 Eigen::VectorXd r = resp_.row (i);
8997 int max_idx = 0 ;
9098 r.maxCoeff (&max_idx);
@@ -95,31 +103,37 @@ void MAPGMM::fit(const std::vector<geo::Vec3>& points, const geo::Pose3D& sensor
95103 determineInlierComponent ();
96104}
97105
98- std::vector<int > MAPGMM::get_labels () const {
106+ std::vector<int > MAPGMM::get_labels () const
107+ {
99108 return labels_;
100109}
101- int MAPGMM::get_inlier_component () const {
110+ int MAPGMM::get_inlier_component () const
111+ {
102112 return inlier_component_;
103113}
104114
105115
106116
107- void MAPGMM::setupPriors (const std::vector<geo::Vec3>& points) {
117+ void MAPGMM::setupPriors (const std::vector<geo::Vec3>& points)
118+ {
108119 if (mu0_.size () != static_cast <size_t >(K_) ||
109120 kappa0_.size () != static_cast <size_t >(K_) ||
110121 Psi0_.size () != static_cast <size_t >(K_) ||
111- nu0_.size () != static_cast <size_t >(K_)) {
122+ nu0_.size () != static_cast <size_t >(K_))
123+ {
112124 CONSOLE_BRIDGE_logError (" Prior vectors not properly initialized in MAPGMM constructor" );
113125 return ;
114126 }
115127 // Dirichlet prior for weights (alpha>1 favors more uniform weights)
116128 alpha_ = params_.alpha ;
117129 Eigen::Vector3d means = Eigen::Vector3d::Zero ();
118- for (const auto & p : points) {
130+ for (const auto & p : points)
131+ {
119132 means += Eigen::Vector3d (p.x , p.y , p.z );
120133 }
121134 means /= points.size ();
122- for (int k = 0 ; k < K_; k++) {
135+ for (int k = 0 ; k < K_; k++)
136+ {
123137 // Same prior for all components initially
124138 mu0_[k] = means;
125139
@@ -132,35 +146,42 @@ void MAPGMM::setupPriors(const std::vector<geo::Vec3>& points) {
132146 }
133147}
134148
135- void MAPGMM::computeBoundingVolume (const Eigen::MatrixXd& data) {
149+ void MAPGMM::computeBoundingVolume (const Eigen::MatrixXd& data)
150+ {
136151 Eigen::Vector3d min_vals = data.colwise ().minCoeff ();
137152 Eigen::Vector3d max_vals = data.colwise ().maxCoeff ();
138153 volume_ = (max_vals - min_vals).prod (); // volume = (xmax - xmin) * (ymax - ymin) * (zmax - zmin)
139- if (volume_ <= 1e-12 ) {
154+ if (volume_ <= 1e-12 )
155+ {
140156 volume_ = 1e-6 ; // avoid division by zero in uniform component
141157 CONSOLE_BRIDGE_logWarn (" MAP-GMM: bounding volume too small, clamping to %g" , volume_);
142158 }
143159}
144160
145- double MAPGMM::eStep (const Eigen::MatrixXd& data, Eigen::MatrixXd& resp_) {
161+ double MAPGMM::eStep (const Eigen::MatrixXd& data, Eigen::MatrixXd& resp_)
162+ {
146163 // Number of data points
147164 int N = data.rows ();
148165 double log_likelihood = 0.0 ;
149166
150167 // Calculate log probabilities for each data point and component (number of components is K_ which is uniform outlier + actual clusters)
151168 Eigen::MatrixXd log_probs (N, K_);
152169
153- for (int k = 0 ; k < K_; k++) {
170+ for (int k = 0 ; k < K_; k++)
171+ {
154172 // Compute log probability of each point under component k
155- for (int i = 0 ; i < N; i++) {
173+ for (int i = 0 ; i < N; i++)
174+ {
156175 Eigen::Vector3d x = data.row (i).transpose ();
157- if (k == 0 ) {
176+ if (k == 0 )
177+ {
158178 // For the outlier component, use a uniform distribution over the bounding volume
159179 double uniform_prob = 1.0 / volume_;
160180 log_probs (i, k) = std::log (weights_[k]) + std::log (uniform_prob);
161181 continue ;
162182 }
163- else {
183+ else
184+ {
164185 Eigen::Vector3d diff = x - means_[k];
165186
166187 // Calculate log probability of x under Gaussian component k (Math is here)
@@ -173,7 +194,8 @@ double MAPGMM::eStep(const Eigen::MatrixXd& data, Eigen::MatrixXd& resp_) {
173194 // Regularize covariance for numerical stability
174195 Eigen::Matrix3d Sigma = covs_[k] + Eigen::Matrix3d::Identity () * 1e-9 ;
175196 double det = Sigma.determinant ();
176- if (det <= 1e-18 || !std::isfinite (det)) {
197+ if (det <= 1e-18 || !std::isfinite (det))
198+ {
177199 Sigma += Eigen::Matrix3d::Identity () * 1e-6 ;
178200 det = Sigma.determinant ();
179201 }
@@ -193,7 +215,8 @@ double MAPGMM::eStep(const Eigen::MatrixXd& data, Eigen::MatrixXd& resp_) {
193215 }
194216
195217 // Calculate responsibilities (and normalize)
196- for (int i = 0 ; i < N; i++) {
218+ for (int i = 0 ; i < N; i++)
219+ {
197220 // Numerical stability: subtract max value
198221 double max_log_prob = log_probs.row (i).maxCoeff ();
199222 Eigen::VectorXd exp_probs = (log_probs.row (i).array () - max_log_prob).exp ();
@@ -209,11 +232,13 @@ double MAPGMM::eStep(const Eigen::MatrixXd& data, Eigen::MatrixXd& resp_) {
209232 return log_likelihood;
210233}
211234
212- void MAPGMM::mStep (const Eigen::MatrixXd& data, const Eigen::MatrixXd& resp_) {
235+ void MAPGMM::mStep (const Eigen::MatrixXd& data, const Eigen::MatrixXd& resp_)
236+ {
213237 int N = data.rows ();
214238
215239 // For each component
216- for (int k = 0 ; k < K_; k++) {
240+ for (int k = 0 ; k < K_; k++)
241+ {
217242 // Component responsibility sum
218243 double Nk = resp_.col (k).sum ();
219244
@@ -224,7 +249,8 @@ void MAPGMM::mStep(const Eigen::MatrixXd& data, const Eigen::MatrixXd& resp_) {
224249
225250 // Calculate mean of data
226251 Eigen::Vector3d mean_data = Eigen::Vector3d::Zero ();
227- for (int i = 0 ; i < N; i++) {
252+ for (int i = 0 ; i < N; i++)
253+ {
228254 mean_data += resp_ (i, k) * data.row (i).transpose ();
229255 }
230256 mean_data /= Nk;
@@ -234,7 +260,8 @@ void MAPGMM::mStep(const Eigen::MatrixXd& data, const Eigen::MatrixXd& resp_) {
234260
235261 // Calculate weighted covariance of data
236262 Eigen::Matrix3d cov_data = Eigen::Matrix3d::Zero ();
237- for (int i = 0 ; i < N; i++) {
263+ for (int i = 0 ; i < N; i++)
264+ {
238265 Eigen::Vector3d diff = data.row (i).transpose () - mean_data;
239266 cov_data += resp_ (i, k) * diff * diff.transpose ();
240267 }
@@ -256,22 +283,27 @@ void MAPGMM::mStep(const Eigen::MatrixXd& data, const Eigen::MatrixXd& resp_) {
256283 weights_ /= weights_.sum ();
257284}
258285
259- void MAPGMM::determineInlierComponent () {
286+ void MAPGMM::determineInlierComponent ()
287+ {
260288 // Compute N_k: total responsibility mass for each component
261289 std::vector<double > Nk (K_, 0.0 );
262290 int N = labels_.size (); // or you can use resp.rows()
263291
264- for (int i = 0 ; i < N; i++) {
265- for (int k = 0 ; k < K_; k++) {
292+ for (int i = 0 ; i < N; i++)
293+ {
294+ for (int k = 0 ; k < K_; k++)
295+ {
266296 Nk[k] += resp_ (i, k); // You'll need to save `resp_` as a class member
267297 }
268298 }
269299
270300 // Find the component with max Nk
271301 // Skip the outlier uniform distribution component (k=0)
272302 inlier_component_ = 1 ;
273- for (int k = 1 ; k < K_; k++) {
274- if (Nk[k] > Nk[inlier_component_]) {
303+ for (int k = 1 ; k < K_; k++)
304+ {
305+ if (Nk[k] > Nk[inlier_component_])
306+ {
275307 inlier_component_ = k;
276308 }
277309 }
0 commit comments