likelihood_field_model.hpp

Go to the documentation of this file.

// Copyright 2022-2023 Ekumen, Inc.
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
//     http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
 
#ifndef BELUGA_SENSOR_LIKELIHOOD_FIELD_MODEL_HPP
#define BELUGA_SENSOR_LIKELIHOOD_FIELD_MODEL_HPP
 
#include <algorithm>
#include <cmath>
#include <random>
#include <vector>
 
#include <beluga/algorithm/distance_map.hpp>
#include <beluga/sensor/data/occupancy_grid.hpp>
#include <beluga/sensor/data/value_grid.hpp>
#include <range/v3/range/conversion.hpp>
#include <range/v3/view/all.hpp>
#include <range/v3/view/transform.hpp>
#include <sophus/se2.hpp>
#include <sophus/so2.hpp>
 
namespace beluga {
 
 
struct LikelihoodFieldModelParam {
 
  double max_obstacle_distance = 100.0;
  double max_laser_distance = 2.0;
  double z_hit = 0.5;
  double z_random = 0.5;
 
  double sigma_hit = 0.2;
};
 
 
template <class OccupancyGrid>
class LikelihoodFieldModel {
 public:
  using state_type = Sophus::SE2d;
  using weight_type = double;
  using measurement_type = std::vector<std::pair<double, double>>;
  using map_type = OccupancyGrid;
  using param_type = LikelihoodFieldModelParam;
 
 
  explicit LikelihoodFieldModel(const param_type& params, const map_type& grid)
      : params_{params},
        likelihood_field_{make_likelihood_field(params, grid)},
        world_to_likelihood_field_transform_{grid.origin().inverse()} {}
 
  [[nodiscard]] const auto& likelihood_field() const { return likelihood_field_; }
 
 
  [[nodiscard]] auto operator()(measurement_type&& points) const {
    return [this, points = std::move(points)](const state_type& state) -> weight_type {
      const auto transform = world_to_likelihood_field_transform_ * state;
      const auto x_offset = transform.translation().x();
      const auto y_offset = transform.translation().y();
      const auto cos_theta = transform.so2().unit_complex().x();
      const auto sin_theta = transform.so2().unit_complex().y();
      const auto unknown_space_occupancy_prob = 1. / params_.max_laser_distance;
      // TODO(glpuga): Investigate why AMCL and QuickMCL both use this formula for the weight.
      // See https://github.com/Ekumen-OS/beluga/issues/153
      const auto unknown_space_occupancy_likelihood_cubed =
          unknown_space_occupancy_prob * unknown_space_occupancy_prob * unknown_space_occupancy_prob;
      return std::transform_reduce(
          points.cbegin(), points.cend(), 1.0, std::plus{},
          [this, x_offset, y_offset, cos_theta, sin_theta,
           unknown_space_occupancy_likelihood_cubed](const auto& point) {
            // Transform the end point of the laser to the grid local coordinate system.
            // Not using Eigen/Sophus because they make the routine x10 slower.
            // See `benchmark_likelihood_field_model.cpp` for reference.
            const auto x = point.first * cos_theta - point.second * sin_theta + x_offset;
            const auto y = point.first * sin_theta + point.second * cos_theta + y_offset;
            // for performance, we store the likelihood already elevated to the cube
            return likelihood_field_.data_near(x, y).value_or(unknown_space_occupancy_likelihood_cubed);
          });
    };
  }
 
 
  void update_map(const map_type& grid) {
    likelihood_field_ = make_likelihood_field(params_, grid);
    world_to_likelihood_field_transform_ = grid.origin().inverse();
  }
 
 private:
  param_type params_;
  ValueGrid2<double> likelihood_field_;
  Sophus::SE2d world_to_likelihood_field_transform_;
 
  static ValueGrid2<double> make_likelihood_field(const LikelihoodFieldModelParam& params, const OccupancyGrid& grid) {
    const auto squared_distance = [&grid,
                                   squared_max_distance = params.max_obstacle_distance * params.max_obstacle_distance](
                                      std::size_t first, std::size_t second) {
      return std::min((grid.coordinates_at(first) - grid.coordinates_at(second)).squaredNorm(), squared_max_distance);
    };
 
    const auto neighborhood = [&grid](std::size_t index) { return grid.neighborhood4(index); };
 
    const auto distance_map = nearest_obstacle_distance_map(grid.obstacle_data(), squared_distance, neighborhood);
 
    const auto to_likelihood = [amplitude =
                                    params.z_hit / (params.sigma_hit * std::sqrt(2 * Sophus::Constants<double>::pi())),
                                two_squared_sigma = 2 * params.sigma_hit * params.sigma_hit,
                                offset = params.z_random / params.max_laser_distance](double squared_distance) {
      assert(two_squared_sigma > 0.0);
      assert(amplitude > 0.0);
      return amplitude * std::exp(-squared_distance / two_squared_sigma) + offset;
    };
 
    // we store the likelihood elevated to the cube to save a few runtime computations
    // when calculating the importance weight
    const auto to_the_cube = [](auto likelihood) { return likelihood * likelihood * likelihood; };
 
    auto likelihood_data = distance_map | ranges::views::transform(to_likelihood) |
                           ranges::views::transform(to_the_cube) | ranges::to<std::vector>;
    return ValueGrid2<double>{std::move(likelihood_data), grid.width(), grid.resolution()};
  }
};
 
}  // namespace beluga
 
#endif