2 //////////////////////////////////////////////////////////////////////////////////
3 // This program is free software: you can redistribute it and/or modify //
4 // it under the terms of the version 3 of the GNU General Public License //
5 // as published by the Free Software Foundation. //
7 // This program is distributed in the hope that it will be useful, but //
8 // WITHOUT ANY WARRANTY; without even the implied warranty of //
9 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU //
10 // General Public License for more details. //
12 // You should have received a copy of the GNU General Public License //
13 // along with this program. If not, see <http://www.gnu.org/licenses/>. //
15 // Written by Francois Fleuret //
16 // (C) Ecole Polytechnique Federale de Lausanne //
17 // Contact <pom@epfl.ch> for comments & bug reports //
18 //////////////////////////////////////////////////////////////////////////////////
26 #include "pom_solver.h"
29 //////////////////////////////////////////////////////////////////////
31 POMSolver::POMSolver(Room *room) : neg(room->view_width(), room->view_height()),
32 neg_view(room->view_width(), room->view_height()),
33 ii_neg(room->view_width(), room->view_height()),
34 ii_neg_view(room->view_width(), room->view_height()) {
35 global_difference.set(global_mu_image_density, global_sigma_image_density);
38 //////////////////////////////////////////////////////////////////////
40 void POMSolver::compute_average_images(int camera,
42 Vector<scalar_t> *proba_absence) {
45 for(int n = 0; n < room->nb_positions(); n++) if((*proba_absence)[n] <= global_proba_ignored) {
46 Rectangle *r = room->avatar(camera, n);
48 neg.multiply_subarray(r->xmin, r->ymin, r->xmax + 1, r->ymax + 1, (*proba_absence)[n]);
52 //////////////////////////////////////////////////////////////////////
54 void POMSolver::add_log_ratio(int camera,
57 Vector<scalar_t> *proba_absence,
58 Vector<scalar_t> *sum) {
60 // Computes the average on the complete picture
62 compute_average_images(camera, room, proba_absence);
64 double s = ii_neg.compute_sum(&neg);
65 double sv = ii_neg_view.compute_sum(&neg, view);
67 scalar_t noise_proba = 0.01; // 1% of the scene can remain unexplained
68 scalar_t average_surface = room->view_width() * room->view_height() * (1 + noise_proba) - s;
69 scalar_t average_diff = average_surface + sv;
71 // Cycles throw all positions and adds the log likelihood ratio to
72 // the total sum for each
74 for(int i = 0; i < room->nb_positions(); i++) {
75 Rectangle *r = room->avatar(camera, i);
77 scalar_t lambda = 1 - 1/(*proba_absence)[i];
79 scalar_t integral_neg = ii_neg.integral(r->xmin, r->ymin, r->xmax + 1, r->ymax + 1);
80 scalar_t average_surface_givpre = average_surface + integral_neg;
81 scalar_t average_surface_givabs = average_surface + lambda * integral_neg;
83 scalar_t integral_neg_view = ii_neg_view.integral(r->xmin, r->ymin, r->xmax + 1, r->ymax + 1);
84 scalar_t average_diff_givpre = average_diff + integral_neg - 2 * integral_neg_view;
85 scalar_t average_diff_givabs = average_diff + lambda * (integral_neg - 2 * integral_neg_view);
87 scalar_t log_mu0 = global_difference.log_proba(average_diff_givabs / average_surface_givabs);
88 scalar_t log_mu1 = global_difference.log_proba(average_diff_givpre / average_surface_givpre);
90 (*sum)[i] += log_mu1 - log_mu0;
96 void POMSolver::solve(Room *room,
97 Vector<scalar_t> *prior,
98 Vector<ProbaView *> *views,
99 Vector<scalar_t> *proba_presence,
101 char *convergence_file_format) {
103 Vector<scalar_t> log_prior_ratio(prior->length());
105 Vector<scalar_t> sum(room->nb_positions());
106 Vector<scalar_t> proba_absence(room->nb_positions());
108 for(int i = 0; i < room->nb_positions(); i++) {
109 log_prior_ratio[i] = log((*prior)[i]/(1 - (*prior)[i]));
110 proba_absence[i] = 1 - (*prior)[i];
115 for(int it = 0; (nb_stab < global_nb_stable_error_for_convergence) && (it < global_max_nb_solver_iterations); it++) {
118 for(int c = 0; c < room->nb_cameras(); c++)
119 add_log_ratio(c, room, (*views)[c], &proba_absence, &sum);
122 for(int i = 0; i < room->nb_positions(); i++) {
123 scalar_t np = global_smoothing_coefficient * proba_absence[i] + (1 - global_smoothing_coefficient) / (1 + exp(log_prior_ratio[i] + sum[i]));
124 if(abs(proba_absence[i] - np) > e) e = abs(proba_absence[i] - np);
125 proba_absence[i] = np;
128 if(e < global_error_max) nb_stab++; else nb_stab = 0;
130 if(convergence_file_format) {
131 char buffer[buffer_size];
132 for(int p = 0; p < room->nb_positions(); p++) (*proba_presence)[p] = 1 - proba_absence[p];
133 for(int c = 0; c < room->nb_cameras(); c++) {
134 pomsprintf(buffer, buffer_size, convergence_file_format, c, nb_frame, it);
135 cout << "Saving " << buffer << "\n"; cout.flush();
136 room->save_stochastic_view(buffer, c, (*views)[c], proba_presence);
142 for(int p = 0; p < room->nb_positions(); p++) (*proba_presence)[p] = 1 - proba_absence[p];