45 double y_max,
double resolution)
48 m_rfgm_run_update_upon_clear(true),
49 m_insertOptions_common(
nullptr),
52 m_hasToRecoverMeanAndCov(true),
54 m_average_normreadings_mean(0),
55 m_average_normreadings_var(0),
56 m_average_normreadings_count(0)
69 const double x_min,
const double x_max,
const double y_min,
70 const double y_max,
const double resolution,
74 x_min, x_max, y_min, y_max, resolution, fill_value);
81 m_gmrf_connectivity = new_connectivity_descriptor;
90 m_insertOptions_common =
91 getCommonInsertOptions();
93 m_average_normreadings_mean = 0;
94 m_average_normreadings_var = 0;
95 m_average_normreadings_count = 0;
111 "[clear] Setting covariance matrix to %ux%u\n",
112 (
unsigned int)(m_size_y * m_size_x),
113 (
unsigned int)(m_size_y * m_size_x));
116 m_insertOptions_common->KF_defaultCellMeanValue,
117 m_insertOptions_common->KF_initialCellStd
123 m_cov.setSize(m_size_y * m_size_x, m_size_y * m_size_x);
126 const double KF_covSigma2 =
127 square(m_insertOptions_common->KF_covSigma);
128 const double res2 =
square(m_resolution);
129 const double std0sqr =
130 square(m_insertOptions_common->KF_initialCellStd);
132 for (
int i = 0; i < m_cov.rows(); i++)
134 int cx1 = (i % m_size_x);
135 int cy1 = (i / m_size_x);
137 for (
int j = i; j < m_cov.cols(); j++)
139 int cx2 = (j % m_size_x);
140 int cy2 = (j / m_size_x);
144 m_cov(i, j) = std0sqr;
150 exp(-0.5 * (res2 *
static_cast<double>(
154 m_cov(j, i) = m_cov(i, j);
163 case mrKalmanApproximate:
165 m_hasToRecoverMeanAndCov =
true;
171 "[CRandomFieldGridMap2D::clear] Resetting compressed cov. "
172 "matrix and cells\n");
174 m_insertOptions_common->KF_defaultCellMeanValue,
175 m_insertOptions_common->KF_initialCellStd
182 const signed W = m_insertOptions_common->KF_W_size;
183 const size_t N = m_map.size();
184 const size_t K = 2 * W * (W + 1) + 1;
186 const double KF_covSigma2 =
187 square(m_insertOptions_common->KF_covSigma);
188 const double std0sqr =
189 square(m_insertOptions_common->KF_initialCellStd);
190 const double res2 =
square(m_resolution);
192 m_stackedCov.setSize(N, K);
197 const double* ptr_first_row = m_stackedCov.get_unsafe_row(0);
199 for (
size_t i = 0; i < N; i++)
201 double* ptr = m_stackedCov.get_unsafe_row(i);
210 for (Acx = 1; Acx <= W; Acx++)
213 exp(-0.5 * (res2 *
static_cast<double>(
218 for (Acy = 1; Acy <= W; Acy++)
219 for (Acx = -W; Acx <= W; Acx++)
223 (res2 *
static_cast<double>(
230 memcpy(ptr, ptr_first_row,
sizeof(
double) * K);
235 "[clear] %ux%u cells done in %.03fms\n",
unsigned(m_size_x),
236 unsigned(m_size_y), 1000 * tictac.
Tac());
246 "[clear] Generating Prior based on 'Squared Differences'\n");
248 "[clear] Initial map dimension: %u cells, x=(%.2f,%.2f) "
250 static_cast<unsigned int>(m_map.size()), getXMin(), getXMax(),
251 getYMin(), getYMax());
258 float res_coef = 1.f;
260 if (this->m_insertOptions_common->GMRF_use_occupancy_information)
262 if (!m_insertOptions_common->GMRF_simplemap_file.empty())
266 this->m_insertOptions_common->GMRF_simplemap_file);
274 !m_insertOptions_common->GMRF_gridmap_image_file.empty())
278 this->m_insertOptions_common->GMRF_gridmap_image_file,
279 this->m_insertOptions_common->GMRF_gridmap_image_res,
280 this->m_insertOptions_common->GMRF_gridmap_image_cx,
281 this->m_insertOptions_common->GMRF_gridmap_image_cy);
284 this->getResolution() /
285 this->m_insertOptions_common->GMRF_gridmap_image_res;
290 "Neither `simplemap_file` nor `gridmap_image_file` "
291 "found in config/mission file. Quitting.");
296 "Resizing m_map to match Occupancy Gridmap dimensions");
303 "Occupancy Gridmap dimensions: x=(%.2f,%.2f)m "
308 "Occupancy Gridmap dimensions: %u cells, cx=%i cy=%i\n\n",
309 static_cast<unsigned int>(
313 "New map dimensions: %u cells, x=(%.2f,%.2f) "
315 static_cast<unsigned int>(m_map.size()), getXMin(),
316 getXMax(), getYMin(), getYMax());
318 "New map dimensions: %u cells, cx=%u cy=%u\n",
319 static_cast<unsigned int>(m_map.size()),
320 static_cast<unsigned int>(getSizeX()),
321 static_cast<unsigned int>(getSizeY()));
324 "./obstacle_map_from_MRPT_for_GMRF.png");
328 const size_t nodeCount = m_map.size();
331 const size_t nPriorFactors =
332 (this->getSizeX() - 1) * this->getSizeY() +
338 "[clear] Generating "
340 <<
" prior factors for a map size of N=" << nodeCount << endl);
343 m_gmrf.initialize(nodeCount);
345 m_mrf_factors_activeObs.clear();
346 m_mrf_factors_activeObs.resize(
349 m_mrf_factors_priors.clear();
354 if (!m_gmrf_connectivity &&
355 this->m_insertOptions_common->GMRF_use_occupancy_information)
359 "MRF Map Dimmensions: %u x %u cells \n",
360 static_cast<unsigned int>(m_size_x),
361 static_cast<unsigned int>(m_size_y));
363 "Occupancy map Dimmensions: %i x %i cells \n",
372 size_t cxoj_min, cxoj_max, cyoj_min, cyoj_max, seed_cxo,
374 size_t cxoi_min, cxoi_max, cyoi_min, cyoi_max, objective_cxo,
376 size_t cxo_min, cxo_max, cyo_min,
380 std::multimap<size_t, size_t>
381 cell_interconnections;
385 for (
size_t j = 0; j < nodeCount;
389 cxoj_min = floor(cx * res_coef);
390 cxoj_max = cxoj_min + ceil(res_coef - 1);
391 cyoj_min = floor(cy * res_coef);
392 cyoj_max = cyoj_min + ceil(res_coef - 1);
394 seed_cxo = cxoj_min + ceil(res_coef / 2 - 1);
395 seed_cyo = cyoj_min + ceil(res_coef / 2 - 1);
399 if (m_Ocgridmap.
getCell(seed_cxo, seed_cyo) < 0.5)
407 m_mrf_factors_activeObs[j].push_back(new_obs);
408 m_gmrf.addConstraint(
409 *m_mrf_factors_activeObs[j]
416 for (
int neighbor = 0; neighbor < 2; neighbor++)
422 if (cx >= (m_size_x - 1))
continue;
430 if (cy >= (m_size_y - 1))
continue;
437 cxoi_min = floor(cxi * res_coef);
438 cxoi_max = cxoi_min + ceil(res_coef - 1);
439 cyoi_min = floor(cyi * res_coef);
440 cyoi_max = cyoi_min + ceil(res_coef - 1);
442 objective_cxo = cxoi_min + ceil(res_coef / 2 - 1);
443 objective_cyo = cyoi_min + ceil(res_coef / 2 - 1);
446 cxo_min =
min(cxoj_min, cxoi_min);
447 cxo_max = max(cxoj_max, cxoi_max);
448 cyo_min =
min(cyoj_min, cyoi_min);
449 cyo_max = max(cyoj_max, cyoi_max);
453 if (exist_relation_between2cells(
454 &m_Ocgridmap, cxo_min, cxo_max, cyo_min,
455 cyo_max, seed_cxo, seed_cyo, objective_cxo,
462 m_insertOptions_common->GMRF_lambdaPrior;
464 m_mrf_factors_priors.push_back(new_prior);
465 m_gmrf.addConstraint(
466 *m_mrf_factors_priors
470 cell_interconnections.insert(
471 std::pair<size_t, size_t>(j, i));
472 cell_interconnections.insert(
473 std::pair<size_t, size_t>(i, j));
479 if (++cx >= m_size_x)
489 m_gmrf_connectivity.get();
491 if (custom_connectivity !=
nullptr)
493 "[CRandomFieldGridMap2D::clear] Initiating prior "
494 "(using user-supplied connectivity pattern)");
497 "[CRandomFieldGridMap2D::clear] Initiating prior "
498 "(fully connected)");
504 size_t cx = 0, cy = 0;
505 for (
size_t j = 0; j < nodeCount; j++)
511 const size_t c_idx_to_check[2] = {cx, cy};
512 const size_t c_idx_to_check_limits[2] = {m_size_x - 1,
514 const size_t c_neighbor_idx_incr[2] = {1, m_size_x};
516 for (
int dir = 0; dir < 2; dir++)
518 if (c_idx_to_check[dir] >= c_idx_to_check_limits[dir])
521 const size_t i = j + c_neighbor_idx_incr[dir];
524 double edge_lamdba = .0;
525 if (custom_connectivity !=
nullptr)
527 const bool is_connected =
529 this, cx, cy, cx + (dir == 0 ? 1 : 0),
530 cy + (dir == 1 ? 1 : 0), edge_lamdba);
531 if (!is_connected)
continue;
536 m_insertOptions_common->GMRF_lambdaPrior;
541 new_prior.
Lambda = edge_lamdba;
543 m_mrf_factors_priors.push_back(new_prior);
544 m_gmrf.addConstraint(
545 *m_mrf_factors_priors.rbegin());
549 if (++cx >= m_size_x)
558 "[clear] Prior built in " << tictac.
Tac() <<
" s ----------");
560 if (m_rfgm_run_update_upon_clear)
563 updateMapEstimation_GMRF();
569 cerr <<
"MAP TYPE NOT RECOGNIZED... QUITTING" << endl;
595 point.
x - m_insertOptions_common->cutoffRadius * 2,
596 point.
x + m_insertOptions_common->cutoffRadius * 2,
597 point.
y - m_insertOptions_common->cutoffRadius * 2,
598 point.
y + m_insertOptions_common->cutoffRadius * 2, defCell);
602 int Ac_cutoff =
round(m_insertOptions_common->cutoffRadius / m_resolution);
603 unsigned Ac_all = 1 + 2 * Ac_cutoff;
604 double minWinValueAtCutOff =
606 m_insertOptions_common->cutoffRadius /
607 m_insertOptions_common->sigma));
609 if (m_DM_lastCutOff != m_insertOptions_common->cutoffRadius ||
610 m_DM_gaussWindow.size() !=
square(Ac_all))
613 "[CRandomFieldGridMap2D::insertObservation_KernelDM_DMV] "
614 "Precomputing window %ux%u\n",
618 double std = m_insertOptions_common->sigma;
621 m_DM_gaussWindow.resize(Ac_all * Ac_all);
622 m_DM_lastCutOff = m_insertOptions_common->cutoffRadius;
627 for (
unsigned cx = 0; cx < Ac_all; cx++)
629 for (
unsigned cy = 0; cy < Ac_all; cy++)
631 dist = m_resolution * sqrt(
633 square(Ac_cutoff + 1 - cx) +
634 square(Ac_cutoff + 1 - cy)));
635 *(it++) = std::exp(-
square(dist / std));
640 "[CRandomFieldGridMap2D::insertObservation_KernelDM_DMV] Done!");
645 const int sensor_cx = x2idx(point.
x);
646 const int sensor_cy = y2idx(point.
y);
650 for (
int Acx = -Ac_cutoff; Acx <= Ac_cutoff; Acx++)
652 for (
int Acy = -Ac_cutoff; Acy <= Ac_cutoff; ++Acy, ++windowIt)
654 const double windowValue = *windowIt;
656 if (windowValue > minWinValueAtCutOff)
658 cell = cellByIndex(sensor_cx + Acx, sensor_cy + Acy);
661 cell->
dm_mean += windowValue * normReading;
664 const double cell_var =
665 square(normReading - computeMeanCellValue_DM_DMV(cell));
680 cutoffRadius(sigma * 3.0),
687 KF_initialCellStd(1.0),
688 KF_observationModelNoise(0),
689 KF_defaultCellMeanValue(0),
692 GMRF_lambdaPrior(0.01f),
694 GMRF_lambdaObs(10.0f),
698 GMRF_lambdaObsLoss(0.0f),
700 GMRF_use_occupancy_information(false),
701 GMRF_simplemap_file(
""),
702 GMRF_gridmap_image_file(
""),
703 GMRF_gridmap_image_res(0.01f),
704 GMRF_gridmap_image_cx(0),
705 GMRF_gridmap_image_cy(0),
707 GMRF_saturate_min(-std::numeric_limits<double>::max()),
708 GMRF_saturate_max(std::numeric_limits<double>::max()),
709 GMRF_skip_variance(false)
720 "sigma = %f\n", sigma);
722 "cutoffRadius = %f\n", cutoffRadius);
724 "R_min = %f\n", R_min);
726 "R_max = %f\n", R_max);
728 "dm_sigma_omega = %f\n", dm_sigma_omega);
731 "KF_covSigma = %f\n", KF_covSigma);
733 "KF_initialCellStd = %f\n", KF_initialCellStd);
735 "KF_observationModelNoise = %f\n",
736 KF_observationModelNoise);
738 "KF_defaultCellMeanValue = %f\n",
739 KF_defaultCellMeanValue);
741 "KF_W_size = %u\n", (
unsigned)KF_W_size);
744 "GMRF_lambdaPrior = %f\n", GMRF_lambdaPrior);
746 "GMRF_lambdaObs = %f\n", GMRF_lambdaObs);
748 "GMRF_lambdaObsLoss = %f\n", GMRF_lambdaObs);
751 "GMRF_use_occupancy_information = %s\n",
752 GMRF_use_occupancy_information ?
"YES" :
"NO");
754 "GMRF_simplemap_file = %s\n",
755 GMRF_simplemap_file.c_str());
757 "GMRF_gridmap_image_file = %s\n",
758 GMRF_gridmap_image_file.c_str());
760 "GMRF_gridmap_image_res = %f\n",
761 GMRF_gridmap_image_res);
763 "GMRF_gridmap_image_cx = %u\n",
764 static_cast<unsigned int>(GMRF_gridmap_image_cx));
766 "GMRF_gridmap_image_cy = %u\n",
767 static_cast<unsigned int>(GMRF_gridmap_image_cy));
778 sigma =
iniFile.read_float(section.c_str(),
"sigma", sigma);
780 iniFile.read_float(section.c_str(),
"cutoffRadius", cutoffRadius);
781 R_min =
iniFile.read_float(section.c_str(),
"R_min", R_min);
782 R_max =
iniFile.read_float(section.c_str(),
"R_max", R_max);
786 iniFile.read_float(section.c_str(),
"KF_covSigma", KF_covSigma);
787 KF_initialCellStd =
iniFile.read_float(
788 section.c_str(),
"KF_initialCellStd", KF_initialCellStd);
789 KF_observationModelNoise =
iniFile.read_float(
790 section.c_str(),
"KF_observationModelNoise", KF_observationModelNoise);
791 KF_defaultCellMeanValue =
iniFile.read_float(
792 section.c_str(),
"KF_defaultCellMeanValue", KF_defaultCellMeanValue);
795 GMRF_lambdaPrior =
iniFile.read_float(
796 section.c_str(),
"GMRF_lambdaPrior", GMRF_lambdaPrior);
798 iniFile.read_float(section.c_str(),
"GMRF_lambdaObs", GMRF_lambdaObs);
799 GMRF_lambdaObsLoss =
iniFile.read_float(
800 section.c_str(),
"GMRF_lambdaObsLoss", GMRF_lambdaObsLoss);
802 GMRF_use_occupancy_information =
iniFile.read_bool(
803 section.c_str(),
"GMRF_use_occupancy_information",
false,
false);
804 GMRF_simplemap_file =
805 iniFile.read_string(section.c_str(),
"simplemap_file",
"",
false);
806 GMRF_gridmap_image_file =
807 iniFile.read_string(section.c_str(),
"gridmap_image_file",
"",
false);
808 GMRF_gridmap_image_res =
809 iniFile.read_float(section.c_str(),
"gridmap_image_res", 0.01f,
false);
810 GMRF_gridmap_image_cx =
811 iniFile.read_int(section.c_str(),
"gridmap_image_cx", 0,
false);
812 GMRF_gridmap_image_cy =
813 iniFile.read_int(section.c_str(),
"gridmap_image_cy", 0,
false);
825 img.saveToFile(filName);
890 double new_x_min,
double new_x_max,
double new_y_min,
double new_y_max,
891 const TRandomFieldCell& defaultValueNewCells,
double additionalMarginMeters)
902 new_x_min, new_x_max, new_y_min, new_y_max, defaultValueNewCells,
903 additionalMarginMeters);
922 "[CRandomFieldGridMap2D::resize] Resizing from %ux%u to %ux%u "
924 static_cast<unsigned>(old_sizeX),
925 static_cast<unsigned>(old_sizeY),
939 for (i = 0; i < N; i++)
944 bool C1_isFromOldMap =
945 Acx_left <= cx1 && cx1 < (Acx_left + old_sizeX) &&
946 Acy_bottom <= cy1 && cy1 < (Acy_bottom + old_sizeY);
950 for (j = i; j < N; j++)
955 bool C2_isFromOldMap =
956 Acx_left <= cx2 && cx2 < (Acx_left + old_sizeX) &&
957 Acy_bottom <= cy2 && cy2 < (Acy_bottom + old_sizeY);
960 if (C1_isFromOldMap && C2_isFromOldMap)
963 unsigned int idx_c1 =
965 old_sizeX * (cy1 - Acy_bottom));
966 unsigned int idx_c2 =
968 old_sizeX * (cy2 - Acy_bottom));
974 ASSERT_((cx1 - Acx_left) < old_sizeX);
975 ASSERT_((cy1 - Acy_bottom) < old_sizeY);
979 ASSERT_((cx2 - Acx_left) < old_sizeX);
980 ASSERT_((cy2 - Acy_bottom) < old_sizeY);
982 ASSERT_(idx_c1 < old_sizeX * old_sizeY);
983 ASSERT_(idx_c2 < old_sizeX * old_sizeY);
986 printf(
"cx1=%u\n",
static_cast<unsigned>(cx1));
987 printf(
"cy1=%u\n",
static_cast<unsigned>(cy1));
988 printf(
"cx2=%u\n",
static_cast<unsigned>(cx2));
989 printf(
"cy2=%u\n",
static_cast<unsigned>(cy2));
992 static_cast<unsigned>(Acx_left));
995 static_cast<unsigned>(Acy_bottom));
998 static_cast<unsigned>(idx_c1));
1001 static_cast<unsigned>(idx_c2)););
1003 m_cov(i, j) = oldCov(idx_c1, idx_c2);
1006 if (i == j)
ASSERT_(idx_c1 == idx_c2);
1008 if (i == j &&
m_cov(i, i) < 0)
1011 "\ni=%u \nj=%i \ncx1,cy1 = %u,%u \n "
1012 "cx2,cy2=%u,%u \nidx_c1=%u \nidx_c2=%u "
1013 "\nAcx_left=%u \nAcy_bottom=%u "
1014 "\nold_sizeX=%u \n",
1015 static_cast<unsigned>(i),
1016 static_cast<unsigned>(j),
1017 static_cast<unsigned>(cx1),
1018 static_cast<unsigned>(cy1),
1019 static_cast<unsigned>(cx2),
1020 static_cast<unsigned>(cy2),
1021 static_cast<unsigned>(idx_c1),
1022 static_cast<unsigned>(idx_c2),
1023 static_cast<unsigned>(Acx_left),
1024 static_cast<unsigned>(Acy_bottom),
1025 static_cast<unsigned>(old_sizeX));
1038 for (i = 0; i < N; i++)
1043 bool C1_isFromOldMap =
1044 Acx_left <= cx1 && cx1 < (Acx_left + old_sizeX) &&
1045 Acy_bottom <= cy1 && cy1 < (Acy_bottom + old_sizeY);
1046 for (j = i; j < N; j++)
1051 bool C2_isFromOldMap =
1052 Acx_left <= cx2 && cx2 < (Acx_left + old_sizeX) &&
1053 Acy_bottom <= cy2 && cy2 < (Acy_bottom + old_sizeY);
1058 if (!C1_isFromOldMap || !C2_isFromOldMap)
1071 static_cast<double>(
1078 "c(i,i)=%e c(j,j)=%e\n",
m_cov(i, i),
1099 printf(
"[CRandomFieldGridMap2D::resize] Done\n");
1114 "[resize] Resizing from %ux%u to %ux%u (%u cells)\n",
1115 static_cast<unsigned>(old_sizeX),
1116 static_cast<unsigned>(old_sizeY),
1123 const size_t N =
m_map.size();
1124 const size_t K = 2 * W * (W + 1) + 1;
1137 const double std0sqr =
1139 double* ptr = &template_row[0];
1141 const double KF_covSigma2 =
1149 for (Acx = 1; Acx <= W; Acx++)
1152 exp(-0.5 * (res2 *
static_cast<double>(
1157 for (Acy = 1; Acy <= W; Acy++)
1158 for (Acx = -W; Acx <= W; Acx++)
1161 exp(-0.5 * (res2 *
static_cast<double>(
1169 for (
size_t i = N - 1; i < N; i--)
1174 const int old_idx_of_i =
1175 (cx - Acx_left) + (cy - Acy_bottom) * old_sizeX;
1185 memcpy(new_row, &template_row[0],
sizeof(
double) * K);
1192 if (old_idx_of_i !=
int(i))
1194 const double* ptr_old =
1197 memcpy(ptr_new, ptr_old,
sizeof(
double) * K);
1221 resize(point.
x - 1, point.
x + 1, point.
y - 1, point.
y + 1, defCell);
1235 int cellIdx =
xy2idx(point.
x, point.
y);
1240 double yk = normReading - cell->
kf_mean;
1242 m_cov(cellIdx, cellIdx) +
1245 ->KF_observationModelNoise);
1246 double sk_1 = 1.0 / sk;
1252 MRPT_LOG_DEBUG(
"[insertObservation_KF] Updating mean values...");
1257 for (i = 0, it =
m_map.begin(); it !=
m_map.end(); ++it, ++i)
1259 it->kf_mean += yk * sk_1 *
m_cov(i, cellIdx);
1267 MRPT_LOG_DEBUG(
"[insertObservation_KF] Updating covariance matrix...");
1271 double* oldCov = (
double*) malloc(
sizeof(
double) * N * N);
1272 double* oldCov_ptr = oldCov;
1273 for (i = 0; i < N; i++)
1275 memcpy(oldCov_ptr,
m_cov.get_unsafe_row(i),
sizeof(
double) * N);
1280 "Copy matrix %ux%u: %.06fms\n", (
unsigned)
m_cov.rows(),
1281 (
unsigned)
m_cov.cols(), tictac.
Tac() * 1000);
1286 const double* oldCov_row_c = oldCov + cellIdx * N;
1287 for (i = 0; i < N; i++)
1289 const double oldCov_i_c = oldCov[i * N + cellIdx];
1290 const double sk_1_oldCov_i_c = sk_1 * oldCov_i_c;
1292 const double* oldCov_row_i = oldCov + i * N;
1293 for (j = i; j < N; j++)
1296 oldCov_row_i[j] - sk_1_oldCov_i_c * oldCov_row_c[j];
1299 m_cov.set_unsafe(i, j, new_cov_ij);
1300 m_cov.set_unsafe(j, i, new_cov_ij);
1305 if (
m_cov(i, i) < 0)
1308 "Wrong insertion in KF! m_cov(%u,%u) = %.5f",
1309 static_cast<unsigned int>(i),
1310 static_cast<unsigned int>(i),
m_cov(i, i));
1314 m_map[i].kf_std = sqrt(new_cov_ij);
1348 DIMs.saveToTextFile(
1351 "% Grid limits: [x_min x_max y_min y_max]\n");
1371 all_means.saveToTextFile(
1375 all_vars.saveToTextFile(
1378 all_confs.saveToTextFile(
1394 for (
size_t i = 0; i <
m_size_y; i++)
1395 for (
size_t j = 0; j <
m_size_x; j++)
1401 MEAN.saveToTextFile(
1403 STDs.saveToTextFile(
1410 filNamePrefix +
std::string(
"_mean_compressed_cov.txt"),
1416 m_cov.saveToTextFile(
1422 CImage img_cov(STDs,
true);
1440 for (
size_t i = 0; i <
m_size_y; ++i)
1442 for (
size_t j = 0; j <
m_size_x; ++j, ++idx)
1447 XYZ(idx, 0) =
idx2x(j);
1448 XYZ(idx, 1) =
idx2y(i);
1454 MEAN.saveToTextFile(
1456 STDs.saveToTextFile(
1462 "% Columns: GRID_X GRID_Y ESTIMATED_Z "
1463 "STD_DEV_OF_ESTIMATED_Z \n");
1480 const double std_times = 3;
1488 FILE* f =
os::fopen(filName.c_str(),
"wt");
1492 f,
"%%-------------------------------------------------------\n");
1493 os::fprintf(f,
"%% File automatically generated using the MRPT method:\n");
1494 os::fprintf(f,
"%%'CRandomFieldGridMap2D::saveAsMatlab3DGraph'\n");
1498 f,
"%% Jose Luis Blanco Claraco, University of Malaga @ 2006-2007\n");
1501 f,
"%%-------------------------------------------------------\n\n");
1503 unsigned int cx, cy;
1504 vector<double> xs, ys;
1586 os::fprintf(f,
"set(S1,'FaceAlpha',0.9,'EdgeAlpha',0.9);\n");
1587 os::fprintf(f,
"set(S2,'FaceAlpha',0.4,'EdgeAlpha',0.4);\n");
1588 os::fprintf(f,
"set(S3,'FaceAlpha',0.2,'EdgeAlpha',0.2);\n");
1591 f,
"set(gca,'PlotBoxAspectRatio',[%f %f %f]);\n",
m_x_max -
m_x_min,
1599 "\nfigure; imagesc(xs,ys,z_mean);axis equal;title('Mean "
1600 "value');colorbar;");
1603 "\nfigure; imagesc(xs,ys,(z_upper-z_mean)/%f);axis equal;title('Std "
1604 "dev of estimated value');colorbar;",
1622 mrpt::make_aligned_shared<mrpt::opengl::CSetOfObjects>();
1639 unsigned int cx, cy;
1640 vector<double> xs, ys;
1660 mrpt::make_aligned_shared<opengl::CSetOfTriangles>();
1661 obj_m->enableTransparency(
true);
1663 mrpt::make_aligned_shared<opengl::CSetOfTriangles>();
1664 obj_v->enableTransparency(
true);
1668 double maxVal_m = 0, minVal_m = 1, AMaxMin_m, maxVal_v = 0,
1669 minVal_v = 1, AMaxMin_v;
1679 minVal_m =
min(minVal_m,
c);
1680 maxVal_m = max(maxVal_m,
c);
1683 minVal_v =
min(minVal_v,
v);
1684 maxVal_v = max(maxVal_v,
v);
1688 AMaxMin_m = maxVal_m - minVal_m;
1689 if (AMaxMin_m == 0) AMaxMin_m = 1;
1690 AMaxMin_v = maxVal_v - minVal_v;
1691 if (AMaxMin_v == 0) AMaxMin_v = 1;
1696 triag.
a[0] = triag.
a[1] = triag.
a[2] =
1706 ASSERT_(cell_x_1y !=
nullptr);
1708 ASSERT_(cell_xy_1 !=
nullptr);
1711 ASSERT_(cell_x_1y_1 !=
nullptr);
1732 double col_xy = c_xy;
1733 double col_x_1y = c_x_1y;
1734 double col_xy_1 = c_xy_1;
1735 double col_x_1y_1 = c_x_1y_1;
1738 triag.
x[0] = xs[cx];
1739 triag.
y[0] = ys[cy];
1741 triag.
x[1] = xs[cx];
1742 triag.
y[1] = ys[cy - 1];
1743 triag.
z[1] = c_xy_1;
1744 triag.
x[2] = xs[cx - 1];
1745 triag.
y[2] = ys[cy - 1];
1746 triag.
z[2] = c_x_1y_1;
1747 jet2rgb(col_xy, triag.
r[0], triag.
g[0], triag.
b[0]);
1748 jet2rgb(col_xy_1, triag.
r[1], triag.
g[1], triag.
b[1]);
1749 jet2rgb(col_x_1y_1, triag.
r[2], triag.
g[2], triag.
b[2]);
1750 obj_m->insertTriangle(triag);
1753 triag.
x[0] = xs[cx];
1754 triag.
y[0] = ys[cy];
1756 triag.
x[1] = xs[cx - 1];
1757 triag.
y[1] = ys[cy - 1];
1758 triag.
z[1] = c_x_1y_1;
1759 triag.
x[2] = xs[cx - 1];
1760 triag.
y[2] = ys[cy];
1761 triag.
z[2] = c_x_1y;
1762 jet2rgb(col_xy, triag.
r[0], triag.
g[0], triag.
b[0]);
1763 jet2rgb(col_x_1y_1, triag.
r[1], triag.
g[1], triag.
b[1]);
1764 jet2rgb(col_x_1y, triag.
r[2], triag.
g[2], triag.
b[2]);
1765 obj_m->insertTriangle(triag);
1780 col_x_1y = v_x_1y / 10;
1782 col_xy_1 = v_xy_1 / 10;
1784 col_x_1y_1 = v_x_1y_1 / 10;
1788 triag.
x[0] = xs[cx];
1789 triag.
y[0] = ys[cy];
1790 triag.
z[0] = c_xy + v_xy;
1791 triag.
x[1] = xs[cx];
1792 triag.
y[1] = ys[cy - 1];
1793 triag.
z[1] = c_xy_1 + v_xy_1;
1794 triag.
x[2] = xs[cx - 1];
1795 triag.
y[2] = ys[cy - 1];
1796 triag.
z[2] = c_x_1y_1 + v_x_1y_1;
1797 jet2rgb(col_xy, triag.
r[0], triag.
g[0], triag.
b[0]);
1798 jet2rgb(col_xy_1, triag.
r[1], triag.
g[1], triag.
b[1]);
1799 jet2rgb(col_x_1y_1, triag.
r[2], triag.
g[2], triag.
b[2]);
1800 obj_v->insertTriangle(triag);
1803 triag.
x[0] = xs[cx];
1804 triag.
y[0] = ys[cy];
1805 triag.
z[0] = c_xy + v_xy;
1806 triag.
x[1] = xs[cx - 1];
1807 triag.
y[1] = ys[cy - 1];
1808 triag.
z[1] = c_x_1y_1 + v_x_1y_1;
1809 triag.
x[2] = xs[cx - 1];
1810 triag.
y[2] = ys[cy];
1811 triag.
z[2] = c_x_1y + v_x_1y;
1812 jet2rgb(col_xy, triag.
r[0], triag.
g[0], triag.
b[0]);
1813 jet2rgb(col_x_1y_1, triag.
r[1], triag.
g[1], triag.
b[1]);
1814 jet2rgb(col_x_1y, triag.
r[2], triag.
g[2], triag.
b[2]);
1815 obj_v->insertTriangle(triag);
1819 meanObj->insert(obj_m);
1820 varObj->insert(obj_v);
1830 mrpt::make_aligned_shared<opengl::CSetOfTriangles>();
1831 obj_m->enableTransparency(
true);
1833 mrpt::make_aligned_shared<opengl::CSetOfTriangles>();
1834 obj_v->enableTransparency(
true);
1838 double maxVal_m = 0, minVal_m = 1, AMaxMin_m, maxVal_v = 0,
1839 minVal_v = 1, AMaxMin_v;
1849 minVal_m =
min(minVal_m,
c);
1850 maxVal_m = max(maxVal_m,
c);
1853 minVal_v =
min(minVal_v,
v);
1854 maxVal_v = max(maxVal_v,
v);
1858 AMaxMin_m = maxVal_m - minVal_m;
1859 if (AMaxMin_m == 0) AMaxMin_m = 1;
1860 AMaxMin_v = maxVal_v - minVal_v;
1861 if (AMaxMin_v == 0) AMaxMin_v = 1;
1866 triag.
a[0] = triag.
a[1] = triag.
a[2] =
1876 ASSERT_(cell_x_1y !=
nullptr);
1878 ASSERT_(cell_xy_1 !=
nullptr);
1881 ASSERT_(cell_x_1y_1 !=
nullptr);
1902 double col_xy = c_xy;
1903 double col_x_1y = c_x_1y;
1904 double col_xy_1 = c_xy_1;
1905 double col_x_1y_1 = c_x_1y_1;
1908 triag.
x[0] = xs[cx];
1909 triag.
y[0] = ys[cy];
1911 triag.
x[1] = xs[cx];
1912 triag.
y[1] = ys[cy - 1];
1913 triag.
z[1] = c_xy_1;
1914 triag.
x[2] = xs[cx - 1];
1915 triag.
y[2] = ys[cy - 1];
1916 triag.
z[2] = c_x_1y_1;
1917 jet2rgb(col_xy, triag.
r[0], triag.
g[0], triag.
b[0]);
1918 jet2rgb(col_xy_1, triag.
r[1], triag.
g[1], triag.
b[1]);
1919 jet2rgb(col_x_1y_1, triag.
r[2], triag.
g[2], triag.
b[2]);
1921 obj_m->insertTriangle(triag);
1924 triag.
x[0] = xs[cx];
1925 triag.
y[0] = ys[cy];
1927 triag.
x[1] = xs[cx - 1];
1928 triag.
y[1] = ys[cy - 1];
1929 triag.
z[1] = c_x_1y_1;
1930 triag.
x[2] = xs[cx - 1];
1931 triag.
y[2] = ys[cy];
1932 triag.
z[2] = c_x_1y;
1933 jet2rgb(col_xy, triag.
r[0], triag.
g[0], triag.
b[0]);
1934 jet2rgb(col_x_1y_1, triag.
r[1], triag.
g[1], triag.
b[1]);
1935 jet2rgb(col_x_1y, triag.
r[2], triag.
g[2], triag.
b[2]);
1936 obj_m->insertTriangle(triag);
1942 double v_x_1y =
min(
1944 double v_xy_1 =
min(
1946 double v_x_1y_1 =
min(
1952 col_x_1y_1 = v_x_1y_1;
1955 triag.
x[0] = xs[cx];
1956 triag.
y[0] = ys[cy];
1958 triag.
x[1] = xs[cx];
1959 triag.
y[1] = ys[cy - 1];
1960 triag.
z[1] = v_xy_1;
1961 triag.
x[2] = xs[cx - 1];
1962 triag.
y[2] = ys[cy - 1];
1963 triag.
z[2] = v_x_1y_1;
1964 jet2rgb(col_xy, triag.
r[0], triag.
g[0], triag.
b[0]);
1965 jet2rgb(col_xy_1, triag.
r[1], triag.
g[1], triag.
b[1]);
1966 jet2rgb(col_x_1y_1, triag.
r[2], triag.
g[2], triag.
b[2]);
1968 obj_v->insertTriangle(triag);
1971 triag.
x[0] = xs[cx];
1972 triag.
y[0] = ys[cy];
1974 triag.
x[1] = xs[cx - 1];
1975 triag.
y[1] = ys[cy - 1];
1976 triag.
z[1] = v_x_1y_1;
1977 triag.
x[2] = xs[cx - 1];
1978 triag.
y[2] = ys[cy];
1979 triag.
z[2] = v_x_1y;
1980 jet2rgb(col_xy, triag.
r[0], triag.
g[0], triag.
b[0]);
1981 jet2rgb(col_x_1y_1, triag.
r[1], triag.
g[1], triag.
b[1]);
1982 jet2rgb(col_x_1y, triag.
r[2], triag.
g[2], triag.
b[2]);
1983 obj_v->insertTriangle(triag);
1987 meanObj->insert(obj_m);
1988 varObj->insert(obj_v);
2014 const double alpha =
2018 const double r_val =
2030 const double alpha =
2034 const double r_val =
2050 const double x,
const double y,
double& out_predict_response,
2051 double& out_predict_response_variance,
bool do_sensor_normalization,
2056 vector<TInterpQuery> queries;
2057 switch (interp_method)
2061 queries[0].cx =
x2idx(
x);
2062 queries[0].cy =
y2idx(
y);
2063 queries[0].coef = 1.0;
2074 queries[0].cx =
x2idx(
x);
2075 queries[0].cy =
y2idx(
y);
2076 queries[0].coef = 1.0;
2095 queries[0].coef = K_1 * (
idx2x(queries[3].cx) -
x) *
2096 (
idx2y(queries[3].cy) -
y);
2097 queries[1].coef = K_1 * (
idx2x(queries[3].cx) -
x) *
2098 (
y -
idx2y(queries[0].cy));
2099 queries[2].coef = K_1 * (
x -
idx2x(queries[0].cx)) *
2100 (
idx2y(queries[3].cy) -
y);
2101 queries[3].coef = K_1 * (
x -
idx2x(queries[0].cx)) *
2102 (
y -
idx2y(queries[0].cy));
2110 for (
size_t i = 0; i < queries.size(); i++)
2179 out_predict_response = 0;
2180 out_predict_response_variance = 0;
2181 for (
size_t i = 0; i < queries.size(); i++)
2183 out_predict_response += queries[i].val * queries[i].coef;
2184 out_predict_response_variance += queries[i].var * queries[i].coef;
2188 if (do_sensor_normalization)
2189 out_predict_response =
2191 out_predict_response *
2206 "Inserting KF2: (" << normReading <<
") at Position" << point << endl);
2209 const size_t K = 2 * W * (W + 1) + 1;
2210 const size_t W21 = 2 * W + 1;
2211 const size_t W21sqr = W21 * W21;
2227 point.
x - Aspace, point.
x + Aspace, point.
y - Aspace, point.
y + Aspace,
2233 const size_t N =
m_map.size();
2250 const int cellIdx =
xy2idx(point.
x, point.
y);
2255 double yk = normReading - cell->
kf_mean;
2261 double sk_1 = 1.0 / sk;
2264 MRPT_LOG_DEBUG(
"[insertObservation_KF2] Updating mean values...");
2273 const int cx_c =
x2idx(point.
x);
2274 const int cy_c =
y2idx(point.
y);
2276 const int Acx0 = max(-W, -cx_c);
2277 const int Acy0 = max(-W, -cy_c);
2278 const int Acx1 =
min(W,
int(
m_size_x) - 1 - cx_c);
2279 const int Acy1 =
min(W,
int(
m_size_y) - 1 - cy_c);
2286 std::vector<int> window_idxs(W21sqr, -1);
2292 for (
int Acy = Acy0; Acy <= 0; Acy++)
2294 int limit_cx = Acy < 0 ? Acx1 : -1;
2296 size_t idx = cx_c + Acx0 +
m_size_x * (cy_c + Acy);
2297 int idx_c_in_idx = -Acy * W21 - Acx0;
2299 int window_idx = Acx0 + W + (Acy + W) * W21;
2301 for (
int Acx = Acx0; Acx <= limit_cx; Acx++)
2304 const double cov_i_c =
m_stackedCov(idx, idx_c_in_idx);
2307 m_map[idx].kf_mean += yk * sk_1 * cov_i_c;
2310 cross_covs_c_i[window_idx] = cov_i_c;
2311 window_idxs[window_idx++] = idx;
2319 for (
int Acy = 0; Acy <= Acy1; Acy++)
2321 int start_cx = Acy > 0 ? Acx0 : 0;
2323 size_t idx = cx_c + start_cx +
m_size_x * (cy_c + Acy);
2327 (W + 1) + (Acy - 1) * W21 + (start_cx + W);
2331 int window_idx = start_cx + W + (Acy + W) * W21;
2333 for (
int Acx = start_cx; Acx <= Acx1; Acx++)
2335 ASSERT_(idx_i_in_c >= 0 && idx_i_in_c <
int(K));
2338 m_map[idx].kf_mean += yk * sk_1 * cov_i_c;
2341 cross_covs_c_i[window_idx] = cov_i_c;
2342 window_idxs[window_idx++] = idx;
2357 for (
size_t i = 0; i < W21sqr; i++)
2359 const int idx_i = window_idxs[i];
2360 if (idx_i < 0)
continue;
2366 const double cov_c_i = cross_covs_c_i[i];
2368 for (
size_t j = i; j < W21sqr; j++)
2370 const int idx_j = window_idxs[j];
2371 if (idx_j < 0)
continue;
2377 const int Ax = cx_j - cx_i;
2382 const int Ay = cy_j - cy_i;
2385 const double cov_c_j = cross_covs_c_i[j];
2389 idx_j_in_i = Ax + W + (Ay - 1) * W21 + W + 1;
2394 double& cov_to_change =
m_stackedCov(idx_i, idx_j_in_i);
2395 double Delta_cov = cov_c_j * cov_c_i * sk_1;
2396 if (i == j && cov_to_change < Delta_cov)
2398 "Negative variance value appeared! Please increase the "
2399 "size of the window "
2400 "(W).\n(m_insertOptions_common->KF_covSigma=%f)",
2403 cov_to_change -= Delta_cov;
2421 const size_t N =
m_map.size();
2422 for (
size_t i = 0; i < N; i++)
2433 out_means.resize(N);
2434 for (
size_t i = 0; i < N; ++i) out_means[i] =
BASE::m_map[i].kf_mean;
2447 out_means.resize(N);
2450 for (
size_t i = 0; i < N; ++i)
2465 ASSERT_(N ==
size_t(in_means.size()));
2466 ASSERT_(N ==
size_t(in_std.size()));
2469 for (
size_t i = 0; i < N; ++i)
2471 m_map[i].kf_mean = in_means[i];
2498 "insertObservation() isn't implemented for selected 'mapType'")
2504 const bool update_map,
const bool time_invariant,
2505 const double reading_stddev)
2523 sensorReading, point, update_map, time_invariant,
2524 reading_stddev == .0
2531 "insertObservation() isn't implemented for selected 'mapType'")
2540 const bool update_map,
const bool time_invariant,
2541 const double reading_information)
2546 const int cellIdx =
xy2idx(point.
x, point.
y);
2554 new_obs.
Lambda = reading_information;
2561 catch (std::exception e)
2563 cerr <<
"Exception while Inserting new Observation: " << e.what()
2576 Eigen::VectorXd x_incr, x_var;
2580 ASSERT_(
size_t(
m_map.size()) ==
size_t(x_incr.size()));
2583 size_t(
m_map.size()) ==
size_t(x_var.size()));
2586 for (
size_t j = 0; j <
m_map.size(); j++)
2590 : std::sqrt(x_var[j]);
2591 m_map[j].gmrf_mean += x_incr[j];
2604 for (
auto ito = obs.begin(); ito != obs.end();)
2606 if (!ito->time_invariant)
2613 if (ito->Lambda < 0)
2616 ito = obs.erase(ito);
2627 size_t cxo_max,
size_t cyo_min,
size_t cyo_max,
const size_t seed_cxo,
2628 const size_t seed_cyo,
const size_t objective_cxo,
2629 const size_t objective_cyo)
2635 cxo_min = max(cxo_min, (
size_t)0);
2636 cxo_max =
min(cxo_max, (
size_t)m_Ocgridmap->
getSizeX() - 1);
2637 cyo_min = max(cyo_min, (
size_t)0);
2638 cyo_max =
min(cyo_max, (
size_t)m_Ocgridmap->
getSizeY() - 1);
2644 if ((seed_cxo < cxo_min) || (seed_cxo >= cxo_max) || (seed_cyo < cyo_min) ||
2645 (seed_cyo >= cyo_max))
2650 if ((objective_cxo < cxo_min) || (objective_cxo >= cxo_max) ||
2651 (objective_cyo < cyo_min) || (objective_cyo >= cyo_max))
2658 if ((m_Ocgridmap->
getCell(seed_cxo, seed_cyo) < 0.5) !=
2659 (m_Ocgridmap->
getCell(objective_cxo, objective_cyo) < 0.5))
2667 cxo_max - cxo_min + 1, cyo_max - cyo_min + 1);
2674 matExp(seed_cxo - cxo_min, seed_cyo - cyo_min) = 1;
2679 while (seedsOld < seedsNew)
2681 seedsOld = seedsNew;
2683 for (
int col = 0; col < matExp.cols(); col++)
2685 for (
int row = 0;
row < matExp.rows();
row++)
2688 if (matExp(
row, col) == 1)
2690 matExp(
row, col) = 2;
2692 for (
int i = -1; i <= 1; i++)
2694 for (
int j = -1; j <= 1; j++)
2697 if ((
int(
row) + j >= 0) &&
2698 (
int(
row) + j <=
int(matExp.rows() - 1)) &&
2699 (
int(col) + i >= 0) &&
2700 (
int(col) + i <=
int(matExp.cols()) - 1))
2702 if (!((i == 0 && j == 0) ||
2703 !(matExp(
row + j, col + i) == 0)))
2707 row + cxo_min, col + cyo_min) <
2711 col + i + cyo_min) < 0.5))
2713 if ((
row + j + cxo_min ==
2715 (col + i + cyo_min ==
2722 matExp(
row + j, col + i) = 1;
2752 return m_parent->m_map[this->node_id].gmrf_mean - this->obsValue;
2756 return this->Lambda;
2765 return m_parent->m_map[this->node_id_i].gmrf_mean -
2766 m_parent->m_map[this->node_id_j].gmrf_mean;
2770 return this->Lambda;
2773 double& dr_dx_i,
double& dr_dx_j)
const