Skip to content

Commit

Permalink
lloyds alg start
Browse files Browse the repository at this point in the history
Signed-off-by: Christian Henkel <[email protected]>
  • Loading branch information
ct2034 committed Nov 1, 2024
1 parent bdc35ce commit d02ddcd
Show file tree
Hide file tree
Showing 7 changed files with 303 additions and 19 deletions.
128 changes: 128 additions & 0 deletions learn/lloyds_algorithm/alg.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,128 @@
#!/usr/bin/env python3


from typing import Dict, Tuple, List
import numpy as np
import matplotlib.pyplot as plt
import cv2
from scipy.spatial import Voronoi


def load_map(filename: str) -> np.ndarray:
"""
Load a png file as map.
"""
img = cv2.imread(filename, cv2.IMREAD_COLOR)
# Convert to binary image
img = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
img = cv2.threshold(img, 200, 255, cv2.THRESH_BINARY)[1]
return img


def initialize_points(img: np.ndarray, k: int) -> np.ndarray:
"""
Initialize the k points.
"""
points = []
for i in range(k):
x = np.random.randint(0, img.shape[0])
y = np.random.randint(0, img.shape[1])
while img[x, y] == 0:
x = np.random.randint(0, img.shape[0])
y = np.random.randint(0, img.shape[1])
points.append((x, y))
return np.array(points)


def compute_centroids(vor: Voronoi, img: np.ndarray) -> np.ndarray:
"""
Compute the centroids of the Voronoi cells.
"""
centroids = []
for region in vor.regions:
if not region or -1 in region:
continue
polygon = [vor.vertices[i] for i in region]
polygon = np.array(polygon)
# Compute the centroid of the polygon
centroid = np.mean(polygon, axis=0)
x, y = int(centroid[0]), int(centroid[1])
# Ensure the centroid is within the image bounds and not in an obstacle
if 0 <= x < img.shape[1] and 0 <= y < img.shape[0] and img[y, x] != 0:
centroids.append([x, y])
return np.array(centroids)


def one_iter(points: np.ndarray, img: np.ndarray, k: int) -> np.ndarray:
"""
Perform Lloyd's algorithm to refine the points.
"""
vor = Voronoi(points)
centroids = compute_centroids(vor, img)
if len(centroids) == k:
points = centroids
return points


def lloyds_algorithm(points: np.ndarray, img: np.ndarray, k: int) -> np.ndarray:
"""
Perform Lloyd's algorithm.
"""
for _ in range(100):
points = one_iter(points, img, k)
return points


def make_colors(n: int) -> List[Tuple[float, float, float]]:
"""
Make n colors.
"""
colors = []
for i in range(n):
r = np.random.rand()
g = np.random.rand()
b = np.random.rand()
colors.append((r, g, b))
return colors


def show_img(
img: np.ndarray,
colors: List[Tuple[float, float, float]],
points: np.ndarray,
_cells: Dict[int, List[Tuple[int, int]]],
):
"""
Show the image. With points.
"""
# image
plt.imshow(img, cmap="gray")

# points
plt.scatter(points[:, 1], points[:, 0], c="r", marker="+")

# cells
print(f"{len(_cells)=}")
if _cells is not None:
for i_p, cell in enumerate(_cells.values()):
print(cell)
if len(cell) == 0:
continue
npcell = np.array(cell)
plt.fill(npcell[:, 1], npcell[:, 0], color=colors[i_p], alpha=0.5)

# simplify plot
plt.axis("off")
plt.subplots_adjust(left=0, right=1, top=1, bottom=0)
plt.show()


if __name__ == "__main__":
np.random.seed(0)
img = load_map("map.png")
n_nodes = 10
n_iterations = 100
points = initialize_points(img, n_nodes)
points = lloyds_algorithm(points, img, n_nodes)
colors = make_colors(n_nodes)
show_img(img, colors, points, {})
Binary file added learn/lloyds_algorithm/map.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
1 change: 1 addition & 0 deletions roadmaps/cvt/lodepng
Submodule lodepng added at d398e0
1 change: 1 addition & 0 deletions roadmaps/gsorm/lodepng
Submodule lodepng added at d398e0
3 changes: 0 additions & 3 deletions roadmaps/ompl/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -19,8 +19,6 @@ FIND_PACKAGE(Boost 1.71.0 REQUIRED COMPONENTS
serialization
python${PYTHON_VERSION}
)
FIND_PACKAGE(Eigen3 REQUIRED)
FIND_PACKAGE(yaml-cpp REQUIRED)
FIND_PACKAGE(ompl REQUIRED)
find_package( OpenCV REQUIRED )

Expand All @@ -45,7 +43,6 @@ add_library(${OMPLPY} SHARED
TARGET_LINK_LIBRARIES(
${OMPLPY}
${Boost_LIBRARIES}
${YAML_CPP_LIBRARIES}
${OpenCV_LIBS}
ompl
)
Expand Down
159 changes: 151 additions & 8 deletions roadmaps/ompl/src/omplpy.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -8,12 +8,6 @@
#include <boost/program_options.hpp>
#include <boost/foreach.hpp>

// Eigen
// #include <Eigen/Core>

// Yaml
#include <yaml-cpp/yaml.h>

// OMPL
#include <ompl/base/spaces/RealVectorStateSpace.h>
#include <ompl/geometric/planners/prm/SPARS.h>
Expand Down Expand Up @@ -160,6 +154,12 @@ class Ompl
}

std::cout << "Loaded environment: " << width << "x" << height << std::endl;
// for (
// std::vector<unsigned char>::const_iterator i = image.begin();
// i != image.end(); ++i)
// {
// std::cout << *i << ' ';
// }

base::StateSpacePtr space(new base::RealVectorStateSpace(dimension));
base::RealVectorBounds bounds(dimension);
Expand All @@ -184,7 +184,7 @@ class Ompl
base::ScopedState<base::RealVectorStateSpace> start(space);
base::ScopedState<base::RealVectorStateSpace> goal(space);

geometric::SPARStwo p(si);
geometric::SPARS p(si);
p.setProblemDefinition(pdef);

p.setDenseDeltaFraction(denseDelta / si->getMaximumExtent());
Expand Down Expand Up @@ -212,7 +212,7 @@ class Ompl
std::cout << "#edges " << boost::num_edges(roadmapOMPL) << std::endl;
// // const auto& roadmapOMPL = p.getDenseGraph();
auto stateProperty = boost::get(
geometric::SPARStwo::vertex_state_t(), roadmapOMPL);
geometric::SPARS::vertex_state_t(), roadmapOMPL);

// output
boost::python::list edges;
Expand Down Expand Up @@ -241,6 +241,138 @@ class Ompl
}
}

return boost::python::make_tuple(edges, duration_ms);

}

boost::python::tuple runPRMStar(
std::string mapFile, int seed,
float denseDelta, float sparseDelta, float stretchFactor,
int maxFailures, double maxTime, int maxIter)
{
ompl::RNG::setSeed(seed);

const size_t dimension = 2;

std::vector<unsigned char> image; // the raw pixels
unsigned width, height;

// decode
unsigned error =
lodepng::decode(image, width, height, mapFile, LCT_GREY, 8);

// if there's an error, display it
if (error) {
std::cerr << "decoder error " << error << ": " <<
lodepng_error_text(error) << std::endl;
return {};
}

std::cout << "Loaded environment: " << width << "x" << height << std::endl;
// for (
// std::vector<unsigned char>::const_iterator i = image.begin();
// i != image.end(); ++i)
// {
// std::cout << *i << ' ';
// }

base::StateSpacePtr space(new base::RealVectorStateSpace(dimension));
base::RealVectorBounds bounds(dimension);
bounds.setLow(0, 0);
bounds.setHigh(0, width);
bounds.setLow(1, 0);
bounds.setHigh(1, height);
space->as<base::RealVectorStateSpace>()->setBounds(bounds);
base::SpaceInformationPtr si(new base::SpaceInformation(space));
base::StateValidityCheckerPtr stateValidityChecker(
new ImageStateValidityChecker(
si, image,
width, height));
si->setStateValidityChecker(stateValidityChecker);
si->setStateValidityCheckingResolution(
1. / (float)std::max(
width,
height) / 5.);
si->setup();

base::ProblemDefinitionPtr pdef(new base::ProblemDefinition(si));
base::ScopedState<base::RealVectorStateSpace> start(space);
base::ScopedState<base::RealVectorStateSpace> goal(space);

geometric::PRMstar p(si);
p.setProblemDefinition(pdef);

// p.setDenseDeltaFraction(denseDelta / si->getMaximumExtent());
// p.setSparseDeltaFraction(sparseDelta / si->getMaximumExtent());
// p.setStretchFactor(stretchFactor);
// p.setMaxFailures(maxFailures);

std::cout << "a" << std::endl;
auto t_start = high_resolution_clock::now();
std::cout << "b" << std::endl;
p.constructRoadmap(
base::IterationTerminationCondition(maxIter));
std::cout << "c" << std::endl;
auto t_stop = high_resolution_clock::now();
std::cout << "d" << std::endl;
auto duration = duration_cast<microseconds>(t_stop - t_start);
std::cout << "e" << std::endl;
float duration_ms = static_cast<float>(duration.count()) / 1000.;
std::cout << "f" << std::endl;

// p.printDebug();

const auto & roadmapOMPL = p.getRoadmap();
int n_vertices = boost::num_vertices(roadmapOMPL);
std::cout << "n_vertices " << n_vertices << std::endl;
int n_edges = boost::num_edges(roadmapOMPL);
std::cout << "n_edges " << n_edges << std::endl;
// // const auto& roadmapOMPL = p.getDenseGraph();
auto stateProperty = boost::get(
geometric::PRMstar::vertex_state_t(), roadmapOMPL);

// output
boost::python::list edges;

// read the edges
BOOST_FOREACH(
const geometric::PRMstar::Edge e, boost::edges(
roadmapOMPL))
{
int i = e.m_source;
int j = e.m_target;

if(i >= n_vertices || j >= n_vertices) {
continue;
}
if(boost::python::len(edges) >= 100) {
break;
}

std::cout << "i: " << i << " j: " << j << std::endl;

base::State * state_i = stateProperty[i];
base::State * state_j = stateProperty[j];
if (state_i && state_j) {
const auto typedState_i =
state_i->as<base::RealVectorStateSpace::StateType>();
float x_i = (*typedState_i)[0];
float y_i = (*typedState_i)[1];

const auto typedState_j =
state_j->as<base::RealVectorStateSpace::StateType>();
float x_j = (*typedState_j)[0];
float y_j = (*typedState_j)[1];

std::cout << "x_i: " << x_i << " y_i: " << y_i << " x_j: " << x_j << " y_j: " << y_j << std::endl;

edges.append(boost::python::make_tuple(i, x_i, y_i, j, x_j, y_j));
} else {
std::cout << state_i << " " << state_j << std::endl;
}
}


return boost::python::make_tuple(edges, duration_ms);

}
Expand All @@ -262,5 +394,16 @@ BOOST_PYTHON_MODULE(libomplpy)
bp::arg("maxFailures"),
bp::arg("maxTime"),
bp::arg("maxIter"))
).def(
"runPRMStar", &Ompl::runPRMStar,
(
bp::arg("mapFile"),
bp::arg("seed"),
bp::arg("denseDelta"),
bp::arg("sparseDelta"),
bp::arg("stretchFactor"),
bp::arg("maxFailures"),
bp::arg("maxTime"),
bp::arg("maxIter"))
);
}
30 changes: 22 additions & 8 deletions roadmaps/plots.py
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,7 @@
MAP_INFLATED_PATH = './roadmaps/benchmark_examples/Berlin_1_256_inflated_1000.png'
OUTPUT_PATH = 'roadmaps/plots/'

N_SAMPLES = 1000
N_SAMPLES = 100


def make_plot(graph, map_fname, path, _):
Expand Down Expand Up @@ -72,19 +72,33 @@ def make_plot(graph, map_fname, path, _):
# )
# make_plot(prm.g, MAP_PATH, os.path.join(OUTPUT_PATH, 'PRM.png'), 'PRM')

# SPARS2
# PRM Star
ompl = Ompl()
edges, duration_ms = ompl.runSparsTwo(
edges, duration_ms = ompl.runPRMStar(
MAP_INFLATED_PATH, #mapFile
1, #seed
2.5, #denseDelta
2034, #seed
3., #denseDelta
10., #sparseDelta
5., #stretchFactor
3., #stretchFactor
1000, #maxFailures
4., #maxTime
40000, #maxIter
2., #maxTime
400, #maxIter
)
print(edges)

# SPARS2
# ompl = Ompl()
# edges, duration_ms = ompl.runSparsTwo(
# MAP_INFLATED_PATH, #mapFile
# 2034, #seed
# 3., #denseDelta
# 10., #sparseDelta
# 3., #stretchFactor
# 1000, #maxFailures
# 2., #maxTime
# 400, #maxIter
# )
# print(edges)

# cleanup
os.remove(MAP_INFLATED_PATH)

0 comments on commit d02ddcd

Please sign in to comment.