diff --git a/.gitattributes b/.gitattributes new file mode 100644 index 0000000..68bc67b --- /dev/null +++ b/.gitattributes @@ -0,0 +1,20 @@ +# Set the default behavior, in case people don't have core.autocrlf set. +* text=auto + +# Explicitly declare text files you want to always be normalized and converted +# to native line endings on checkout. +*.tex text +*.bib text +*.svg text +*.py text +*.vbs text +*.cpp text +*.hpp text +Makefile text + +# Declare files that will always have CRLF line endings on checkout. +*.sln text eol=crlf + +# Denote all files that are truly binary and should not be modified. +*.png binary +*.jpg binary diff --git a/.gitignore b/.gitignore index 259148f..6ab33f6 100644 --- a/.gitignore +++ b/.gitignore @@ -30,3 +30,10 @@ *.exe *.out *.app + +# Extra directories +.idea/ +.vscode/ +/cmake-build-* +/build* +/external* diff --git a/.travis.yml b/.travis.yml new file mode 100644 index 0000000..d479abf --- /dev/null +++ b/.travis.yml @@ -0,0 +1,34 @@ +dist: trusty +sudo: false +language: cpp +cache: ccache +matrix: + include: + - os: linux + compiler: gcc-7 + addons: + apt: + sources: + - ubuntu-toolchain-r-test + packages: + - gcc-7 + - g++-7 + - xorg-dev + - libglu1-mesa-dev + env: + - MATRIX_EVAL="export CC=gcc-7 && CXX=g++-7 && CONFIG=Debug && NPROC=2" + - os: osx + compiler: clang + env: + - MATRIX_EVAL="export CONFIG=Debug && NPROC=2" + +install: +- if [[ "$TRAVIS_OS_NAME" == "osx" ]]; then brew install ccache; fi +- if [[ "$TRAVIS_OS_NAME" == "osx" ]]; then export PATH="/usr/local/opt/ccache/libexec:$PATH"; fi +- eval "${MATRIX_EVAL}" + +script: +- mkdir build +- cd build +- cmake -DCMAKE_BUILD_TYPE=$CONFIG .. +- make -j ${NPROC} diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..9dfc544 --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,73 @@ +################################################################################ +# General Information +################################################################################ + +cmake_minimum_required(VERSION 3.8) +project(voroffset) + +################################################################################ + +set(VOROFFSET_EXTERNAL ${CMAKE_CURRENT_SOURCE_DIR}/external/) +list(APPEND CMAKE_MODULE_PATH ${CMAKE_CURRENT_SOURCE_DIR}/cmake) + +if(NOT CMAKE_BUILD_TYPE) + message(STATUS "No build type selected, default to Release") + set(CMAKE_BUILD_TYPE "Release") +endif() + +################################################################################ + +# Use folder in Visual Studio +set_property(GLOBAL PROPERTY USE_FOLDERS ON) + +# Color output +include(UseColors) + +# Export compile flags(used for autocompletion of the C++ code) +set(CMAKE_EXPORT_COMPILE_COMMANDS 1) + +# CMake plugin for vscode +include(CMakeToolsHelpers OPTIONAL) + +# Generate position independent code +set(CMAKE_POSITION_INDEPENDENT_CODE ON) + +# Build shared or static? +# set(BUILD_SHARED_LIBRARIES OFF) + +################################################################################ + +# CMake options +option(VOROFFSET_WITH_SANITIZERS "Enable sanitizers" OFF) +option(VOROFFSET_WITH_TBB "Enable TBB" ON) + +option(SANITIZE_ADDRESS "Sanitize Address" OFF) +option(SANITIZE_MEMORY "Sanitize Memory" OFF) +option(SANITIZE_THREAD "Sanitize Thread" OFF) +option(SANITIZE_UNDEFINED "Sanitize Undefined" OFF) + +# Override cached options +# set(SANITIZE_ADDRESS OFF CACHE BOOL "" FORCE) +# set(VOROFFSET_WITH_TBB ON CACHE BOOL "" FORCE) + +################################################################################ +# Dependencies +################################################################################ + +# Sanitizers +if(VOROFFSET_WITH_SANITIZERS) + list(APPEND CMAKE_MODULE_PATH ${VOROFFSET_EXTERNAL}/sanitizers-cmake/cmake) +endif() + +include(VoroffsetDependencies) + +################################################################################ + +# 2D version prototype +add_subdirectory(src/vor2d) + +# 3D version prototype +add_subdirectory(src/vor3d) + +# Binary executables +add_subdirectory(app) diff --git a/README.md b/README.md new file mode 100644 index 0000000..5694c6b --- /dev/null +++ b/README.md @@ -0,0 +1,61 @@ +# Voroffset + +Discrete mesh offsetting based on half-space Voronoi diagrams and Power diagrams. + +### Compilation + +Dependencies are downloaded automatically by CMake at configuration time. To compile the code, just type: + +```bash +mkdir build +cd build +cmake -j8 +``` + +### Running the code + +See possible options with: + +```bash +./offset3d -h +``` + +Example usage: + +``` +./offset3d filigree.ply -n 512 -p 10 -r 5 -x dilation +``` + +##### offset2d + +Takes a .svg file as input, and saves the result of the dilation as a quad-mesh (.obj). + +##### offset3d + +Takes a triangle mesh as input (.stl, .obj, .off), and saves the dilated output as a mesh (quad-mesh with .obj, hex-mesh with .mesh, etc.) + +``` +Offset3D +Usage: ./offset3d [OPTIONS] input [output] + +Positionals: + input TEXT Input model + output TEXT=output.obj Output model + +Options: + -h,--help Print this help message and exit + -i,--input TEXT Input model + -o,--output TEXT=output.obj Output model + -j,--json TEXT Output json file + -d,--dexels_size FLOAT=1 Size of a dexel (in mm) + -n,--num_dexels INT=256 Number of dexels (-1 to use dexel size instead) + -p,--padding INT Padding (in #dexels) + -t,--num_thread UINT=6 Number of threads + -r,--radius FLOAT=8 Dilation/erosion radius (in #dexels) + -m,--method TEXT in {brute_force,ours} + The method to use + -x,--apply TEXT in {closing,dilation,erosion,noop,opening}=dilation + Morphological operation to apply + -f,--force Overwrite output file + -u,--radius_in_mm Radius is given in mm instead +``` diff --git a/app/CMakeLists.txt b/app/CMakeLists.txt new file mode 100644 index 0000000..cac3ec8 --- /dev/null +++ b/app/CMakeLists.txt @@ -0,0 +1,2 @@ +add_subdirectory(cli2d) +add_subdirectory(cli3d) diff --git a/app/cli2d/CMakeLists.txt b/app/cli2d/CMakeLists.txt new file mode 100644 index 0000000..be264c3 --- /dev/null +++ b/app/cli2d/CMakeLists.txt @@ -0,0 +1,27 @@ +################################################################################ + +cmake_minimum_required(VERSION 3.3) +project(offset2d) + +################################################################################ + +# Add executable +add_executable(${PROJECT_NAME} offset2d.cpp) + +# Let's get a little bit paranoid +include(SetWarnings) +target_compile_options(${PROJECT_NAME} PRIVATE ${ALL_WARNINGS}) + +# Sanitizers +if(VOROFFSET_WITH_SANITIZERS) + add_sanitizers(${PROJECT_NAME}) +endif() + +# Dependencies +target_link_libraries(${PROJECT_NAME} PRIVATE vor2d yimg::yimg CLI11::CLI11) + +# Output directory for binaries +set_target_properties(${PROJECT_NAME} + PROPERTIES + RUNTIME_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}" +) diff --git a/app/cli2d/offset2d.cpp b/app/cli2d/offset2d.cpp new file mode 100644 index 0000000..f2c60e8 --- /dev/null +++ b/app/cli2d/offset2d.cpp @@ -0,0 +1,83 @@ +//////////////////////////////////////////////////////////////////////////////// +#include +#include +#include +#include +#include +#include +//////////////////////////////////////////////////////////////////////////////// + +#define WHITE YImage::YPixel({255, 255, 255, 255}) +#define BLACK YImage::YPixel({0, 0, 0, 255}) + +bool pixel_is_white(YImage::YPixel p) { + return p.r > 128 || p.g > 128 || p.b > 128 || p.a > 128; +} + +int main(int argc, char *argv[]) { + // Default arguments + struct { + std::string input; + std::string output = "out.png"; + double radius = 0; + bool erode = false; + bool force = false; + bool transpose = false; + bool negate = false; + } args; + + // Parse arguments + CLI::App app("Offset2D"); + app.add_option("input,-i,--input", args.input, "Input image.") + ->required() + ->check(CLI::ExistingFile); + app.add_option("output,-o,--output", args.output, "Output image."); + app.add_option("-r,--radius", args.radius, "Dilation/erosion radius."); + app.add_flag("-e,--erode", args.erode, "Erode instead of dilate."); + app.add_flag("-f,--force", args.force, "Overwrite output file."); + app.add_flag("-t,--transpose", args.transpose, "Transpose input image."); + app.add_flag("-n,--negate", args.negate, "Negate input image."); + try { + app.parse(argc, argv); + } catch (const CLI::ParseError &e) { + return app.exit(e); + } + + // Load SVG + vor::DoubleCompressedImage dexels = vor::create_dexels(args.input.c_str()); + + // Pre-processing operations + if (args.transpose) { + dexels.transposeInPlace(); + } + if (args.negate) { + dexels.negate(); + } + + // Offset + if (args.radius > 0) { + std::cout << "-- Performing offset by radius r = " << args.radius + << std::endl; + if (args.erode) { + dexels.erode(args.radius); + } else { + dexels.dilate(args.radius); + } + } + + // Save output image + if (std::ifstream(args.output)) { + // Output file exists! + if (args.force) { + std::cout << "-- Overwriting output file: " << args.output << std::endl; + vor::dexel_dump(args.output.c_str(), dexels); + } else { + std::cerr << "-- Output file already exists. Please use -f to force overwriting." + << std::endl; + } + } else { + std::cout << "-- Saving" << std::endl; + vor::dexel_dump(args.output.c_str(), dexels); + } + return 0; +} diff --git a/app/cli3d/CMakeLists.txt b/app/cli3d/CMakeLists.txt new file mode 100644 index 0000000..bdcc948 --- /dev/null +++ b/app/cli3d/CMakeLists.txt @@ -0,0 +1,27 @@ +################################################################################ + +cmake_minimum_required(VERSION 3.3) +project(offset3d) + +################################################################################ + +# Add executable +add_executable(${PROJECT_NAME} offset3d.cpp) + +# Let's get a little bit paranoid +include(SetWarnings) +target_compile_options(${PROJECT_NAME} PRIVATE ${ALL_WARNINGS}) + +# Sanitizers +if(VOROFFSET_WITH_SANITIZERS) + add_sanitizers(${PROJECT_NAME}) +endif() + +# Dependencies +target_link_libraries(${PROJECT_NAME} PRIVATE vor3d CLI11::CLI11 json::json) + +# Output directory for binaries +set_target_properties(${PROJECT_NAME} + PROPERTIES + RUNTIME_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}" +) diff --git a/app/cli3d/offset3d.cpp b/app/cli3d/offset3d.cpp new file mode 100644 index 0000000..a4a9a17 --- /dev/null +++ b/app/cli3d/offset3d.cpp @@ -0,0 +1,179 @@ +//////////////////////////////////////////////////////////////////////////////// +#include +#include +#include +#include +#include +#ifdef USE_TBB +#include +#endif +#include +#include +#include +#include +#include +#include +#include +#include +//////////////////////////////////////////////////////////////////////////////// + +int main(int argc, char * argv[]) { + // Default arguments + struct { + std::string input; + std::string output_mesh = "output.obj"; + std::string output_json = ""; + std::string method = "ours"; + std::string operation = "dilation"; + double radius = 8; + double dexels_size = 1; + int padding = 0; + int num_dexels = 256; + unsigned int num_thread = std::max(1u, std::thread::hardware_concurrency()); + bool force = false; + bool radius_in_mm = false; + } args; + + // Parse arguments + CLI::App app("Offset3D"); + app.add_option("input,-i,--input", args.input, "Input model")->required()->check(CLI::ExistingFile); + app.add_option("output,-o,--output", args.output_mesh, "Output model", true); + app.add_option("-j,--json", args.output_json, "Output json file"); + app.add_option("-d,--dexels_size", args.dexels_size, "Size of a dexel (in mm)", true); + app.add_option("-n,--num_dexels", args.num_dexels, "Number of dexels (-1 to use dexel size instead)", true); + app.add_option("-p,--padding", args.padding, "Padding (in #dexels)"); + app.add_option("-t,--num_thread", args.num_thread, "Number of threads", true); + app.add_option("-r,--radius", args.radius, "Dilation/erosion radius (in #dexels)", true); + app.add_set("-m,--method", args.method, {"ours","brute_force"}, "The method to use"); + app.add_set("-x,--apply", args.operation, {"noop","dilation","erosion","closing","opening"}, + "Morphological operation to apply", true); + app.add_flag("-f,--force", args.force, "Overwrite output file"); + app.add_flag("-u,--radius_in_mm", args.radius_in_mm, "Radius is given in mm instead"); + try { + app.parse(argc, argv); + } catch (const CLI::ParseError &e) { + return app.exit(e); + } + + // Declare stuff + vor3d::CompressedVolume input; + vor3d::CompressedVolume output; + double time = 0; + double time_1, time_2; + nlohmann::json json_data; + + // Load input model and dexelize + Timer t; + { + GEO::Stopwatch W("Dexelize"); + input = vor3d::create_dexels(args.input, args.dexels_size, args.padding, args.num_dexels); + } + + // Compute radius in #dexels if needed, and display some stats + if (args.radius_in_mm) { + args.radius /= input.spacing(); + } + GEO::Logger::out("Stats") << "Grid size: " << input.gridSize().transpose() << std::endl; + GEO::Logger::out("Stats") << "Spacing (in mm): " << input.spacing() << std::endl; + GEO::Logger::out("Stats") << "Origin (in mm): " << input.origin().transpose() << std::endl; + GEO::Logger::out("Stats") << "Extent (in mm): " << input.extent().transpose() << std::endl; + GEO::Logger::out("Stats") << "Radius (in #dexels): " << args.radius << std::endl; + GEO::Logger::out("Stats") << "Number of threads: " << args.num_thread << std::endl; + + if (!args.output_json.empty()) { + json_data = { + { "method", args.method }, + { "num_threads", args.num_thread }, + { "model_name", args.input }, + { "voxel_size", input.spacing() }, + { "padding", input.padding() }, + { "num_dexels", args.num_dexels }, + { "radius", args.radius }, + { "grid_size", { input.gridSize()(0),input.gridSize()(1) } }, + { "num_segments", input.numSegments() }, + { "operation", args.operation } + }; + } + + #ifdef USE_TBB + tbb::task_scheduler_init init(args.num_thread); + #endif + + // Create offset operator + GEO::Logger::div("Offseting"); + std::unique_ptr op; + if (args.method == "ours") { + op = std::make_unique(); + } else if (args.method == "brute_force") { + op = std::make_unique(); + } else { + GEO::Logger::err("Offset") << "Invalid method: " << args.method << std::endl; + return 1; + } + assert(op); + + // Apply operation + if (args.operation == "noop") { + output = input; + } else if (args.operation == "erosion") { + GEO::Stopwatch W("Erosion"); + op->erosion(input, output, args.radius, time_1, time_2); + } else if (args.operation == "dilation") { + GEO::Stopwatch W("Dilation"); + op->dilation(input, output, args.radius, time_1, time_2); + } else if (args.operation == "closing") { + GEO::Stopwatch W("Closing"); + vor3d::CompressedVolume tmp; + op->dilation(input, tmp, args.radius, time_1, time_2); + op->erosion(tmp, output, args.radius, time_1, time_2); + } else if (args.operation == "opening") { + GEO::Stopwatch W("Opening"); + vor3d::CompressedVolume tmp; + op->erosion(input, tmp, args.radius, time_1, time_2); + op->dilation(tmp, output, args.radius, time_1, time_2); + } else { + throw std::invalid_argument("Operation"); + } + + GEO::Logger::div("Saving"); + if (!args.output_mesh.empty()) { + if (std::ifstream(args.output_mesh)) { + // Output file exists! + if (args.force) { + GEO::Logger::out("Save") << "Overwriting output file: " << args.output_mesh << std::endl; + GEO::Stopwatch W("Save"); + vor3d::dexel_dump(args.output_mesh, output); + } else { + GEO::Logger::out("Save") << "Output mesh already exists. Please use -f to force overwriting." << std::endl; + } + } else { + GEO::Stopwatch W("Save"); + std::ofstream out(args.output_mesh.c_str()); + vor3d::dexel_dump(args.output_mesh, output); + } + } + + time = t.get(); + if (!args.output_json.empty()) { + json_data["time"] = time; + json_data["time_first_pass"] = time_1; + json_data["time_second_pass"] = time_2; + } + + if (!args.output_json.empty()) { + if (std::ifstream(args.output_json)) { + // Output file exists! + if (args.force) { + GEO::Logger::out("Save") << "Overwriting output file: " << args.output_json << std::endl; + std::ofstream o(args.output_json); + o << std::setw(4) << json_data << std::endl; + } else { + GEO::Logger::out("Save") << "Output json already exists. Please use -f to force overwriting." << std::endl; + } + } else { + std::ofstream o(args.output_json); + o << std::setw(4) << json_data << std::endl; + } + } + return 0; +} diff --git a/cmake/DownloadProject.CMakeLists.cmake.in b/cmake/DownloadProject.CMakeLists.cmake.in new file mode 100644 index 0000000..89be4fd --- /dev/null +++ b/cmake/DownloadProject.CMakeLists.cmake.in @@ -0,0 +1,17 @@ +# Distributed under the OSI-approved MIT License. See accompanying +# file LICENSE or https://github.com/Crascit/DownloadProject for details. + +cmake_minimum_required(VERSION 2.8.2) + +project(${DL_ARGS_PROJ}-download NONE) + +include(ExternalProject) +ExternalProject_Add(${DL_ARGS_PROJ}-download + ${DL_ARGS_UNPARSED_ARGUMENTS} + SOURCE_DIR "${DL_ARGS_SOURCE_DIR}" + BINARY_DIR "${DL_ARGS_BINARY_DIR}" + CONFIGURE_COMMAND "" + BUILD_COMMAND "" + INSTALL_COMMAND "" + TEST_COMMAND "" +) diff --git a/cmake/DownloadProject.cmake b/cmake/DownloadProject.cmake new file mode 100644 index 0000000..e300f42 --- /dev/null +++ b/cmake/DownloadProject.cmake @@ -0,0 +1,182 @@ +# Distributed under the OSI-approved MIT License. See accompanying +# file LICENSE or https://github.com/Crascit/DownloadProject for details. +# +# MODULE: DownloadProject +# +# PROVIDES: +# download_project( PROJ projectName +# [PREFIX prefixDir] +# [DOWNLOAD_DIR downloadDir] +# [SOURCE_DIR srcDir] +# [BINARY_DIR binDir] +# [QUIET] +# ... +# ) +# +# Provides the ability to download and unpack a tarball, zip file, git repository, +# etc. at configure time (i.e. when the cmake command is run). How the downloaded +# and unpacked contents are used is up to the caller, but the motivating case is +# to download source code which can then be included directly in the build with +# add_subdirectory() after the call to download_project(). Source and build +# directories are set up with this in mind. +# +# The PROJ argument is required. The projectName value will be used to construct +# the following variables upon exit (obviously replace projectName with its actual +# value): +# +# projectName_SOURCE_DIR +# projectName_BINARY_DIR +# +# The SOURCE_DIR and BINARY_DIR arguments are optional and would not typically +# need to be provided. They can be specified if you want the downloaded source +# and build directories to be located in a specific place. The contents of +# projectName_SOURCE_DIR and projectName_BINARY_DIR will be populated with the +# locations used whether you provide SOURCE_DIR/BINARY_DIR or not. +# +# The DOWNLOAD_DIR argument does not normally need to be set. It controls the +# location of the temporary CMake build used to perform the download. +# +# The PREFIX argument can be provided to change the base location of the default +# values of DOWNLOAD_DIR, SOURCE_DIR and BINARY_DIR. If all of those three arguments +# are provided, then PREFIX will have no effect. The default value for PREFIX is +# CMAKE_BINARY_DIR. +# +# The QUIET option can be given if you do not want to show the output associated +# with downloading the specified project. +# +# In addition to the above, any other options are passed through unmodified to +# ExternalProject_Add() to perform the actual download, patch and update steps. +# The following ExternalProject_Add() options are explicitly prohibited (they +# are reserved for use by the download_project() command): +# +# CONFIGURE_COMMAND +# BUILD_COMMAND +# INSTALL_COMMAND +# TEST_COMMAND +# +# Only those ExternalProject_Add() arguments which relate to downloading, patching +# and updating of the project sources are intended to be used. Also note that at +# least one set of download-related arguments are required. +# +# If using CMake 3.2 or later, the UPDATE_DISCONNECTED option can be used to +# prevent a check at the remote end for changes every time CMake is run +# after the first successful download. See the documentation of the ExternalProject +# module for more information. It is likely you will want to use this option if it +# is available to you. Note, however, that the ExternalProject implementation contains +# bugs which result in incorrect handling of the UPDATE_DISCONNECTED option when +# using the URL download method or when specifying a SOURCE_DIR with no download +# method. Fixes for these have been created, the last of which is scheduled for +# inclusion in CMake 3.8.0. Details can be found here: +# +# https://gitlab.kitware.com/cmake/cmake/commit/bdca68388bd57f8302d3c1d83d691034b7ffa70c +# https://gitlab.kitware.com/cmake/cmake/issues/16428 +# +# If you experience build errors related to the update step, consider avoiding +# the use of UPDATE_DISCONNECTED. +# +# EXAMPLE USAGE: +# +# include(DownloadProject) +# download_project(PROJ googletest +# GIT_REPOSITORY https://github.com/google/googletest.git +# GIT_TAG master +# UPDATE_DISCONNECTED 1 +# QUIET +# ) +# +# add_subdirectory(${googletest_SOURCE_DIR} ${googletest_BINARY_DIR}) +# +#======================================================================================== + + +set(_DownloadProjectDir "${CMAKE_CURRENT_LIST_DIR}") + +include(CMakeParseArguments) + +function(download_project) + + set(options QUIET) + set(oneValueArgs + PROJ + PREFIX + DOWNLOAD_DIR + SOURCE_DIR + BINARY_DIR + # Prevent the following from being passed through + CONFIGURE_COMMAND + BUILD_COMMAND + INSTALL_COMMAND + TEST_COMMAND + ) + set(multiValueArgs "") + + cmake_parse_arguments(DL_ARGS "${options}" "${oneValueArgs}" "${multiValueArgs}" ${ARGN}) + + # Hide output if requested + if (DL_ARGS_QUIET) + set(OUTPUT_QUIET "OUTPUT_QUIET") + else() + unset(OUTPUT_QUIET) + message(STATUS "Downloading/updating ${DL_ARGS_PROJ}") + endif() + + # Set up where we will put our temporary CMakeLists.txt file and also + # the base point below which the default source and binary dirs will be. + # The prefix must always be an absolute path. + if (NOT DL_ARGS_PREFIX) + set(DL_ARGS_PREFIX "${CMAKE_BINARY_DIR}") + else() + get_filename_component(DL_ARGS_PREFIX "${DL_ARGS_PREFIX}" ABSOLUTE + BASE_DIR "${CMAKE_CURRENT_BINARY_DIR}") + endif() + if (NOT DL_ARGS_DOWNLOAD_DIR) + set(DL_ARGS_DOWNLOAD_DIR "${DL_ARGS_PREFIX}/${DL_ARGS_PROJ}-download") + endif() + + # Ensure the caller can know where to find the source and build directories + if (NOT DL_ARGS_SOURCE_DIR) + set(DL_ARGS_SOURCE_DIR "${DL_ARGS_PREFIX}/${DL_ARGS_PROJ}-src") + endif() + if (NOT DL_ARGS_BINARY_DIR) + set(DL_ARGS_BINARY_DIR "${DL_ARGS_PREFIX}/${DL_ARGS_PROJ}-build") + endif() + set(${DL_ARGS_PROJ}_SOURCE_DIR "${DL_ARGS_SOURCE_DIR}" PARENT_SCOPE) + set(${DL_ARGS_PROJ}_BINARY_DIR "${DL_ARGS_BINARY_DIR}" PARENT_SCOPE) + + # The way that CLion manages multiple configurations, it causes a copy of + # the CMakeCache.txt to be copied across due to it not expecting there to + # be a project within a project. This causes the hard-coded paths in the + # cache to be copied and builds to fail. To mitigate this, we simply + # remove the cache if it exists before we configure the new project. It + # is safe to do so because it will be re-generated. Since this is only + # executed at the configure step, it should not cause additional builds or + # downloads. + file(REMOVE "${DL_ARGS_DOWNLOAD_DIR}/CMakeCache.txt") + + # Create and build a separate CMake project to carry out the download. + # If we've already previously done these steps, they will not cause + # anything to be updated, so extra rebuilds of the project won't occur. + # Make sure to pass through CMAKE_MAKE_PROGRAM in case the main project + # has this set to something not findable on the PATH. + configure_file("${_DownloadProjectDir}/DownloadProject.CMakeLists.cmake.in" + "${DL_ARGS_DOWNLOAD_DIR}/CMakeLists.txt") + execute_process(COMMAND ${CMAKE_COMMAND} -G "${CMAKE_GENERATOR}" + -D "CMAKE_MAKE_PROGRAM:FILE=${CMAKE_MAKE_PROGRAM}" + . + RESULT_VARIABLE result + ${OUTPUT_QUIET} + WORKING_DIRECTORY "${DL_ARGS_DOWNLOAD_DIR}" + ) + if(result) + message(FATAL_ERROR "CMake step for ${DL_ARGS_PROJ} failed: ${result}") + endif() + execute_process(COMMAND ${CMAKE_COMMAND} --build . + RESULT_VARIABLE result + ${OUTPUT_QUIET} + WORKING_DIRECTORY "${DL_ARGS_DOWNLOAD_DIR}" + ) + if(result) + message(FATAL_ERROR "Build step for ${DL_ARGS_PROJ} failed: ${result}") + endif() + +endfunction() diff --git a/cmake/SetWarnings.cmake b/cmake/SetWarnings.cmake new file mode 100644 index 0000000..38ad09a --- /dev/null +++ b/cmake/SetWarnings.cmake @@ -0,0 +1,145 @@ +################################################################################ +cmake_minimum_required(VERSION 2.6.3) +################################################################################ +# See comments and discussions here: +# http://stackoverflow.com/questions/5088460/flags-to-enable-thorough-and-verbose-g-warnings +################################################################################ + +set(MY_FLAGS + -Wall + -Wextra + -pedantic + + # -Wconversion + #-Wunsafe-loop-optimizations # broken with C++11 loops + -Wunused + + -Wno-long-long + -Wpointer-arith + -Wformat=2 + -Wuninitialized + -Wcast-qual + -Wmissing-noreturn + -Wmissing-format-attribute + -Wredundant-decls + + -Werror=implicit + -Werror=nonnull + -Werror=init-self + -Werror=main + -Werror=missing-braces + -Werror=sequence-point + -Werror=return-type + -Werror=trigraphs + -Werror=array-bounds + -Werror=write-strings + -Werror=address + -Werror=int-to-pointer-cast + -Werror=pointer-to-int-cast + + -Wno-unused-variable + -Wunused-but-set-variable + -Wno-unused-parameter + + #-Weffc++ + -Wno-old-style-cast + -Wno-sign-conversion + #-Wsign-conversion + -Wno-sign-compare + + -Wshadow + + -Wstrict-null-sentinel + -Woverloaded-virtual + -Wsign-promo + -Wstack-protector + -Wstrict-aliasing + -Wstrict-aliasing=2 + -Wswitch-default + -Wswitch-enum + -Wswitch-unreachable + + -Wcast-align + -Wdisabled-optimization + #-Winline # produces warning on default implicit destructor + -Winvalid-pch + -Wmissing-include-dirs + -Wpacked + -Wno-padded + # -Wstrict-overflow + # -Wstrict-overflow=2 + + -Wctor-dtor-privacy + -Wlogical-op + -Wnoexcept + -Woverloaded-virtual + -Wundef + + -Wnon-virtual-dtor + -Wdelete-non-virtual-dtor + + ########### + # GCC 6.1 # + ########### + + -Wnull-dereference + -fdelete-null-pointer-checks + -Wduplicated-cond + -Wmisleading-indentation + + #-Weverything + + ########################### + # Enabled by -Weverything # + ########################### + + #-Wdocumentation + #-Wdocumentation-unknown-command + #-Wfloat-equal + #-Wcovered-switch-default + + #-Wglobal-constructors + #-Wexit-time-destructors + #-Wmissing-variable-declarations + #-Wextra-semi + #-Wweak-vtables + #-Wno-source-uses-openmp + #-Wdeprecated + #-Wnewline-eof + #-Wmissing-prototypes + + #-Wno-c++98-compat + #-Wno-c++98-compat-pedantic + + ########################### + # Need to check if those are still valid today + ########################### + + #-Wimplicit-atomic-properties + #-Wmissing-declarations + #-Wmissing-prototypes + #-Wstrict-selector-match + #-Wundeclared-selector + #-Wunreachable-code + + # Not a warning, but enable link-time-optimization + #-flto + + # Gives meaningful stack traces + -fno-omit-frame-pointer +) + +# Flags above don't make sense for MSVC +if(${CMAKE_CXX_COMPILER_ID} STREQUAL "MSVC") + set(MY_FLAGS) +endif() + +include(CheckCXXCompilerFlag) + +set(ALL_WARNINGS) +foreach(FLAG ${MY_FLAGS}) + check_cxx_compiler_flag("${FLAG}" IS_SUPPORTED_${FLAG}) + if(IS_SUPPORTED_${FLAG}) + set(ALL_WARNINGS ${ALL_WARNINGS} ${FLAG}) + endif() +endforeach() diff --git a/cmake/UseColors.cmake b/cmake/UseColors.cmake new file mode 100644 index 0000000..32c3c7e --- /dev/null +++ b/cmake/UseColors.cmake @@ -0,0 +1,45 @@ +################################################################################ +# When using Clang, there is nothing to do: colors are enabled by default +# When using GCC >= 4.9, colored diagnostics can be enabled natively +# When using an older version, one can use gccfilter (a perl script) +# +# I do not recommend using gccfilter as of now (May 2014), because it seems to +# be bugged. But if you still want to try, here is how to install it on Ubuntu: +# +# +# 1) Download the perl script and add it to you $PATH +# mkdir -p ~/.local/bin +# wget -P ~/.local/bin http://www.mixtion.org/gccfilter/gccfilter +# chmod +x ~/local/bin/gccfilter +# echo 'PATH="$HOME/.local/bin:$PATH"' >> ~/.bashrc +# +# 2) Install the dependencies +# * Term::ANSIColor +# sudo cpan +# cpan> install Term::ANSIColor +# * The module "Getopt::Long" is included in "perl-base" +# * For Getopt::ArgvFile and Regexp::Common ... +# sudo apt-get install libgetopt-argvfile-perl libregexp-common-perl +# +################################################################################ + +if(CMAKE_COMPILER_IS_GNUCXX) + # If GCC >= 4.9, just activate the right option + if(NOT CMAKE_CXX_COMPILER_VERSION VERSION_LESS 4.9) + message(STATUS "GCC >= 4.9 detected, enabling colored diagnostics") + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fdiagnostics-color=auto") + set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -fdiagnostics-color=auto") + set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -fdiagnostics-color=auto") + return() + endif() + # If GCC < 4.9, maybe we can use gccfilter + find_program(GCC_FILTER gccfilter) + if(GCC_FILTER) + option(COLOR_GCC "Use GCCFilter to color compiler output messages" OFF) + set(COLOR_GCC_OPTIONS "-c -r -w" CACHE STRING "Arguments that are passed to gccfilter when output coloring is switchend on. Defaults to -c -r -w.") + if(COLOR_GCC) + set_property(GLOBAL PROPERTY RULE_LAUNCH_COMPILE "${GCC_FILTER} ${COLOR_GCC_OPTIONS}") + message(STATUS "Using gccfilter for colored diagnostics") + endif() + endif() +endif() diff --git a/cmake/VoroffsetDependencies.cmake b/cmake/VoroffsetDependencies.cmake new file mode 100644 index 0000000..cf4212f --- /dev/null +++ b/cmake/VoroffsetDependencies.cmake @@ -0,0 +1,83 @@ +################################################################################ +# CMake download helpers +################################################################################ + +# Download external dependencies +include(VoroffsetDownloadExternal) + +# Sanitizers +if(VOROFFSET_WITH_SANITIZERS) + voroffset_download_sanitizers() + find_package(Sanitizers) +endif() + +################################################################################ +# Required libraries +################################################################################ + +# Eigen +if(NOT TARGET Eigen3::Eigen) + add_library(voroffset_eigen INTERFACE) + voroffset_download_eigen() + target_include_directories(voroffset_eigen SYSTEM INTERFACE ${VOROFFSET_EXTERNAL}/eigen) + add_library(Eigen3::Eigen ALIAS voroffset_eigen) +endif() + +# CL11 +if(NOT TARGET CLI11::CLI11) + voroffset_download_cli11() + add_subdirectory(${VOROFFSET_EXTERNAL}/cli11) + target_compile_definitions(CLI11 INTERFACE -DCLI11_STD_OPTIONAL=0) + target_compile_definitions(CLI11 INTERFACE -DCLI11_EXPERIMENTAL_OPTIONAL=0) +endif() + +# Geogram +if(NOT TARGET geogram::geogram) + voroffset_download_geogram() + set(GEOGRAM_SEARCH_PATHS "${VOROFFSET_EXTERNAL}/geogram") + include(geogram) +endif() + +# json +if(NOT TARGET json::json) + add_library(voroffset_json INTERFACE) + voroffset_download_json() + target_include_directories(voroffset_json SYSTEM INTERFACE ${VOROFFSET_EXTERNAL}/json) + target_include_directories(voroffset_json SYSTEM INTERFACE ${VOROFFSET_EXTERNAL}/json/nlohmann) + add_library(json::json ALIAS voroffset_json) +endif() + +# Nanosvg +if(NOT TARGET nanosvg::nanosvg) + voroffset_download_nanosvg() + add_library(nanosvg_nanosvg INTERFACE) + add_library(nanosvg::nanosvg ALIAS nanosvg_nanosvg) + target_include_directories(nanosvg_nanosvg SYSTEM INTERFACE ${VOROFFSET_EXTERNAL}/nanosvg/src) +endif() + +# TBB +if(VOROFFSET_WITH_TBB AND NOT TARGET tbb::tbb) + set(TBB_BUILD_STATIC ON CACHE BOOL " " FORCE) + set(TBB_BUILD_SHARED OFF CACHE BOOL " " FORCE) + set(TBB_BUILD_TBBMALLOC OFF CACHE BOOL " " FORCE) + set(TBB_BUILD_TBBMALLOC_PROXY OFF CACHE BOOL " " FORCE) + set(TBB_BUILD_TESTS OFF CACHE BOOL " " FORCE) + + voroffset_download_tbb() + add_subdirectory(${VOROFFSET_EXTERNAL}/tbb tbb) + set_property(TARGET tbb_static tbb_def_files PROPERTY FOLDER "dependencies") + set_target_properties(tbb_static PROPERTIES COMPILE_FLAGS "-Wno-implicit-fallthrough -Wno-missing-field-initializers -Wno-unused-parameter -Wno-keyword-macro") + + add_library(voroffset_tbb INTERFACE) + target_include_directories(voroffset_tbb SYSTEM INTERFACE ${VOROFFSET_EXTERNAL}/tbb/include) + target_link_libraries(voroffset_tbb INTERFACE tbb_static) + add_library(tbb::tbb ALIAS voroffset_tbb) +endif() + +# yimg +if(NOT TARGET ymg::ymg) + voroffset_download_yimg() + add_library(yimg_yimg STATIC ${VOROFFSET_EXTERNAL}/yimg/YImage.cpp) + target_include_directories(yimg_yimg PUBLIC ${VOROFFSET_EXTERNAL}/yimg) + add_library(yimg::yimg ALIAS yimg_yimg) +endif() diff --git a/cmake/VoroffsetDownloadExternal.cmake b/cmake/VoroffsetDownloadExternal.cmake new file mode 100644 index 0000000..135b4cf --- /dev/null +++ b/cmake/VoroffsetDownloadExternal.cmake @@ -0,0 +1,95 @@ +################################################################################ +include(DownloadProject) + +# With CMake 3.8 and above, we can hide warnings about git being in a +# detached head by passing an extra GIT_CONFIG option +if(NOT (${CMAKE_VERSION} VERSION_LESS "3.8.0")) + set(VOROFFSET_EXTRA_OPTIONS "GIT_CONFIG advice.detachedHead=false") +else() + set(VOROFFSET_EXTRA_OPTIONS "") +endif() + +# Shortcut function +function(voroffset_download_project name) + download_project( + PROJ ${name} + SOURCE_DIR ${VOROFFSET_EXTERNAL}/${name} + DOWNLOAD_DIR ${VOROFFSET_EXTERNAL}/.cache/${name} + QUIET + ${VOROFFSET_EXTRA_OPTIONS} + ${ARGN} + ) +endfunction() + +################################################################################ + +## Eigen +function(voroffset_download_eigen) + voroffset_download_project(eigen + URL http://bitbucket.org/eigen/eigen/get/3.3.7.tar.gz + URL_MD5 f2a417d083fe8ca4b8ed2bc613d20f07 + ) +endfunction() + +## CppNumericalSolvers +function(voroffset_download_cppoptlib) + voroffset_download_project(CppNumericalSolvers + GIT_REPOSITORY https://github.com/PatWie/CppNumericalSolvers.git + GIT_TAG 7eddf28fa5a8872a956d3c8666055cac2f5a535d + ) +endfunction() + +## CLI11 +function(voroffset_download_cli11) + voroffset_download_project(cli11 + URL https://github.com/CLIUtils/CLI11/archive/v1.7.1.tar.gz + ) +endfunction() + +## nanosvg +function(voroffset_download_nanosvg) + voroffset_download_project(nanosvg + GIT_REPOSITORY https://github.com/memononen/nanosvg.git + GIT_TAG 2b08deeb553c723d151f908d786c64136d26d576 + ) +endfunction() + +## Sanitizers +function(voroffset_download_sanitizers) + voroffset_download_project(sanitizers-cmake + GIT_REPOSITORY https://github.com/arsenm/sanitizers-cmake.git + GIT_TAG 99e159ec9bc8dd362b08d18436bd40ff0648417b + ) +endfunction() + +## Json +function(voroffset_download_json) + voroffset_download_project(json + URL https://github.com/nlohmann/json/releases/download/v3.1.2/include.zip + URL_HASH SHA256=495362ee1b9d03d9526ba9ccf1b4a9c37691abe3a642ddbced13e5778c16660c + ) +endfunction() + +## geogram +function(voroffset_download_geogram) + voroffset_download_project(geogram + GIT_REPOSITORY https://github.com/alicevision/geogram.git + GIT_TAG v1.6.9 + ) +endfunction() + +# yimg +function(voroffset_download_yimg) + voroffset_download_project(yimg + GIT_REPOSITORY https://github.com/jdumas/yimg.git + GIT_TAG 6d924fdb038c86242d272d69d5e4525086fe60cc + ) +endfunction() + +## TBB +function(voroffset_download_tbb) + voroffset_download_project(tbb + GIT_REPOSITORY https://github.com/wjakob/tbb.git + GIT_TAG 4c3ffe5a5f37addef0dd6283c74c4402a3b4ebc9 + ) +endfunction() diff --git a/cmake/geogram.cmake b/cmake/geogram.cmake new file mode 100644 index 0000000..f350287 --- /dev/null +++ b/cmake/geogram.cmake @@ -0,0 +1,77 @@ +################################################################################ +# Find Geogram and build it as part of the current build +################################################################################ + +if(TARGET geogram::geogram) + return() +endif() + +################################################################################ + +find_path(GEOGRAM_SOURCE_INCLUDE_DIR + geogram/basic/common.h + PATHS ${GEOGRAM_SEARCH_PATHS} + PATH_SUFFIXES src/lib + NO_DEFAULT_PATH +) + +set(GEOGRAM_ROOT ${GEOGRAM_SOURCE_INCLUDE_DIR}/../..) + +################################################################################ + +if(${CMAKE_SYSTEM_NAME} MATCHES "Windows") + set(VORPALINE_ARCH_64 TRUE CACHE BOOL "" FORCE) + set(VORPALINE_PLATFORM Win-vs-generic CACHE STRING "" FORCE) +elseif(${CMAKE_SYSTEM_NAME} MATCHES "Linux") + set(VORPALINE_PLATFORM Linux64-gcc-dynamic CACHE STRING "" FORCE) +elseif(${CMAKE_SYSTEM_NAME} MATCHES "Darwin") + set(VORPALINE_PLATFORM Darwin-clang-dynamic CACHE STRING "" FORCE) +endif() + +option(GEOGRAM_WITH_GRAPHICS "Viewers and geogram_gfx library" OFF) +option(GEOGRAM_WITH_LEGACY_NUMERICS "Legacy numerical libraries" OFF) +option(GEOGRAM_WITH_HLBFGS "Non-linear solver (Yang Liu's HLBFGS)" ON) +option(GEOGRAM_WITH_TETGEN "Tetrahedral mesher (Hang Si's TetGen)" OFF) +option(GEOGRAM_WITH_TRIANGLE "Triangle mesher (Jonathan Shewchuk's triangle)" OFF) +option(GEOGRAM_WITH_EXPLORAGRAM "Experimental code (hexahedral meshing vpipeline and optimal transport)" OFF) +option(GEOGRAM_WITH_LUA "Built-in LUA interpreter" OFF) +option(GEOGRAM_LIB_ONLY "Libraries only (no example programs/no viewer)" ON) +option(GEOGRAM_WITH_FPG "Predicate generator (Sylvain Pion's FPG)" OFF) +option(GEOGRAM_USE_SYSTEM_GLFW3 "Use the version of GLFW3 installed in the system if found" OFF) + +################################################################################ + +add_subdirectory(${GEOGRAM_ROOT} geogram) +target_include_directories(geogram SYSTEM PUBLIC ${GEOGRAM_SOURCE_INCLUDE_DIR}) + +################################################################################ + +if(${CMAKE_SYSTEM_NAME} MATCHES "Windows") + # remove warning for multiply defined symbols (caused by multiple + # instanciations of STL templates) + #target_compile_options(geogram INTERFACE /wd4251) + + # remove all unused stuff from windows.h + # target_compile_definitions(geogram INTERFACE -DWIN32_LEAN_AND_MEAN) + # target_compile_definitions(geogram INTERFACE -DVC_EXTRALEAN) + + # do not define a min() and a max() macro, breaks + # std::min() and std::max() !! + target_compile_definitions(geogram INTERFACE -DNOMINMAX) + + # we want M_PI etc... + target_compile_definitions(geogram INTERFACE -D_USE_MATH_DEFINES) + + if(NOT VORPALINE_BUILD_DYNAMIC) + # If we use static library, we link with the static C++ runtime. + foreach(config ${CMAKE_CONFIGURATION_TYPES}) + string(TOUPPER ${config} config) + string(REPLACE /MD /MT CMAKE_C_FLAGS_${config} "${CMAKE_C_FLAGS_${config}}") + string(REPLACE /MD /MT CMAKE_CXX_FLAGS_${config} "${CMAKE_CXX_FLAGS_${config}}") + endforeach() + endif() +endif() + +################################################################################ + +add_library(geogram::geogram ALIAS geogram) diff --git a/src/vor2d/CMakeLists.txt b/src/vor2d/CMakeLists.txt new file mode 100644 index 0000000..7d91b22 --- /dev/null +++ b/src/vor2d/CMakeLists.txt @@ -0,0 +1,53 @@ +################################################################################ + +cmake_minimum_required(VERSION 3.3) +project(vor2d) + +################################################################################ + +add_library(${PROJECT_NAME} + Common.cpp + Common.h + CompressedImage.cpp + CompressedImage.h + DistanceTransform.cpp + DistanceTransform.h + Dexelize.cpp + Dexelize.h + DoubleCompressedImage.cpp + DoubleCompressedImage.h + DoubleVoronoi.cpp + DoubleVoronoi.h + Image.h + MorphologyOperators.cpp + MorphologyOperators.h + Voronoi.cpp + Voronoi.h +) + +target_include_directories(${PROJECT_NAME} PUBLIC ..) + +################################################################################ + +# Let's get a little bit paranoid +include(SetWarnings) +target_compile_options(${PROJECT_NAME} PRIVATE ${ALL_WARNINGS}) + +# Use C++14 +target_compile_features(${PROJECT_NAME} PUBLIC cxx_std_14) + +# Generate position independent code +set_target_properties(${PROJECT_NAME} PROPERTIES POSITION_INDEPENDENT_CODE ON) + +# Sanitizers +if(VOROFFSET_WITH_SANITIZERS) + add_sanitizers(${PROJECT_NAME}) +endif() + +# Dependencies +target_link_libraries(${PROJECT_NAME} + PUBLIC + Eigen3::Eigen + geogram::geogram + nanosvg::nanosvg +) diff --git a/src/vor2d/Common.cpp b/src/vor2d/Common.cpp new file mode 100644 index 0000000..4012d32 --- /dev/null +++ b/src/vor2d/Common.cpp @@ -0,0 +1,17 @@ +//////////////////////////////////////////////////////////////////////////////// +#include "vor2d/Common.h" +#include +#include +//////////////////////////////////////////////////////////////////////////////// + +[[noreturn]] void voroffset::assertion_failed( + const std::string& condition_string, + const std::string& file, int line) +{ + std::ostringstream os; + os << "Assertion failed: " << condition_string << ".\n"; + os << "File: " << file << ",\n"; + os << "Line: " << line; + + throw std::runtime_error(os.str()); +} diff --git a/src/vor2d/Common.h b/src/vor2d/Common.h new file mode 100644 index 0000000..2fc5841 --- /dev/null +++ b/src/vor2d/Common.h @@ -0,0 +1,81 @@ +#pragma once + +//////////////////////////////////////////////////////////////////////////////// +#include +#include +#include +//////////////////////////////////////////////////////////////////////////////// + +namespace voroffset +{ + +// 2D index coordinates +typedef Eigen::Vector2i Vector2i; +typedef double Scalar; +typedef std::complex PointF; +typedef std::vector Curve; +struct m_Segment +{ + Scalar _begin_pt, _end_pt, _radius; + m_Segment() {}; + m_Segment(Scalar begin_pt, Scalar end_pt, Scalar radius) + :_begin_pt(begin_pt), _end_pt(end_pt), _radius(radius) + {} + bool operator <(const m_Segment &s) const + { + return (_begin_pt + _end_pt < s._begin_pt + s._end_pt) || (_begin_pt + _end_pt == s._begin_pt + s._end_pt && _radius < s._radius); + } +}; + +// 2D array of pairs of integers +typedef Eigen::Matrix MatrixXci; + +// Custom assert function +[[noreturn]] void assertion_failed( + const std::string& condition_string, + const std::string& file, int line); + +} // namespace voroffset + +// Shortcut alias +namespace vor = voroffset; + +//////////////////////////////////////////////////////////////////////////////// + +#define vor_assert_msg(x, s) \ +do \ +{ \ + if (!(x)) \ + { \ + voroffset::assertion_failed((s), __FILE__, __LINE__); \ + } \ +} while(0) + +// ----------------------------------------------------------------------------- + +#define vor_assert(x) \ +do \ +{ \ + if(!(x)) \ + { \ + voroffset::assertion_failed(#x, __FILE__, __LINE__); \ + } \ +} while (0) + + +#ifdef voroffset_NO_DEBUG + +#define vor_debug(x) + +#else + +#define vor_debug(x) \ +do \ +{ \ + if(!(x)) \ + { \ + voroffset::assertion_failed(#x, __FILE__, __LINE__); \ + } \ +} while (0) + +#endif diff --git a/src/vor2d/CompressedImage.cpp b/src/vor2d/CompressedImage.cpp new file mode 100644 index 0000000..ae9c5dc --- /dev/null +++ b/src/vor2d/CompressedImage.cpp @@ -0,0 +1,647 @@ +//////////////////////////////////////////////////////////////////////////////// +#include "vor2d/CompressedImage.h" +#include "vor2d/Voronoi.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +//////////////////////////////////////////////////////////////////////////////// + +using namespace voroffset; + +//////////////////////////////////////////////////////////////////////////////// +// Input/Output +//////////////////////////////////////////////////////////////////////////////// + +// Import from 2D image +void voroffset::CompressedImage::fromImage(std::function read_pixel_func) +{ + for (int j = 0; j < (int) m_Rays.size(); ++j) + { + m_Rays[j].clear(); + bool prev = false; + for (int i = 0; i < m_XSize; ++i) + { + bool curr = read_pixel_func(i, j); + if (curr ^ prev) + { + m_Rays[j].push_back(i); + } + prev = curr; + } + if (prev != false) + { + m_Rays[j].push_back(m_XSize); + } + } +} + +// Export to 2D image +void voroffset::CompressedImage::toImage(std::function write_pixel_func) const +{ + for (int j = 0; j < (int) m_Rays.size(); ++j) + { + int i = 0; + bool status = false; + for (const auto & val : m_Rays[j]) + { + for (; i < val; ++i) + { + write_pixel_func(i, j, status); + } + status = !status; + } + for (; i < m_XSize; ++i) + { + write_pixel_func(i, j, status); + } + } +} + +#include + +// Marshalling +void voroffset::CompressedImage::save(std::ostream &out) const +{ + //out << std::setprecision(std::numeric_limits::max_digits10); + out << m_XSize << " " << m_Rays.size() << "\n"; + for (const auto & row : m_Rays) + { + out << row.size(); + for (const auto & val : row) + out << ' ' << val; + { + } + out << "\n"; + } +} + +void voroffset::CompressedImage::load(std::istream &in) +{ + unsigned int num_cols; + in >> m_XSize >> num_cols; + m_Rays.assign(num_cols, {}); + for (auto & row : m_Rays) + { + size_t size; in >> size; + row.resize(size); + for (auto & val : row) + { + in >> val; + } + } +} + +//////////////////////////////////////////////////////////////////////////////// +// Access +//////////////////////////////////////////////////////////////////////////////// + +// Resize +void voroffset::CompressedImage::resize(int w, int h) +{ + m_XSize = w; + m_Rays.assign(h, {}); +} + +// Sanity check +bool voroffset::CompressedImage::isValid() const +{ + for (const auto & row : m_Rays) + { + if (row.size() % 2 != 0) + { + std::cerr << "row size % 2 != 0" << std::endl; + return false; + } + Scalar lastEvent = -1; + for (const auto &val : row) + { + if (val <= lastEvent) + { + std::cerr << "val > lastEvent" << std::endl; + return false; + } + lastEvent = val; + } + if (lastEvent > m_XSize) + { + std::cerr << "lastEvent <= xsize" << std::endl; + return false; + } + } + return true; +} + +// Empty? +bool voroffset::CompressedImage::empty() const +{ + for (const auto & row : m_Rays) + { + if (!row.empty()) { return false; } + } + return true; +} + +// Apply a function to each segment in the structure +void voroffset::CompressedImage::iterate(std::function func) const +{ + for (int i = 0; i < (int) m_Rays.size(); ++i) + { + if (m_Rays.size() % 2 != 0) + { + throw std::runtime_error("Invalid number of events along a ray."); + } + for (size_t j = 0; 2*j < m_Rays[i].size(); ++j) + { + func(i, m_Rays[i][2*j], m_Rays[i][2*j+1]); + } + } +} + +// Writes a pixel at a given location +// return true if the value was changed +bool voroffset::CompressedImage::write(int i, int j, bool newVal) +{ + if (j < 0 || j >= (int) m_Rays.size()) + { + return false; + } + std::vector &row = m_Rays[j]; + if (row.empty()) + { + // poke in an empty line + if (!newVal) { return false; } + row = {i, i+1}; + return true; + } + if (row[0] > i) + { + // poke to the left of the segments + if (!newVal) { return false; } + if (row[0] > i+1) + { + row.insert(row.begin(), 2, i+1); + } + row[0] = i; + return true; + } + size_t a = 0; + size_t b = row.size() - 1; + if (row[b] <= i) + { + // poke to the right of the segments + if (!newVal) { return false; } + if (row[b] == i) + { + row[b] = i+1; + } + else + { + row.push_back(i); + row.push_back(i+1); + } + return true; + } + // Assumes that row[a] ≤ i < row[b] + while (a + 1 < b) + { + size_t c = (a + b) / 2; + if (row[c] <= i) + { + a = c; + } + else + { + b = c; + } + } + bool isInside = ((a % 2) == 0); + if (isInside == newVal) { return false; } + int len = row[a+1] - row[a]; + if (len == 1) + { // erase all + row.erase(row.begin()+a, row.begin()+a+2); + } + else if (i == row[a]) + { // move left end to the right + row[a] = i+1; + } + else if (i+1 == row[a+1]) + { // move right end to the left + row[a+1] = i; + } + else + { + row.insert(row.begin()+a+1, 2, i+1); + row[a+1] = i; + } + return true; +} + +// Dichotomy search +bool voroffset::CompressedImage::at(int i, int j) const +{ + if (j < 0 || j >= (int) m_Rays.size()) + { + return false; + } + const std::vector &row = m_Rays[j]; + size_t a = 0; + size_t b = row.size() - 1; + if (row.empty() || row[a] > i) { return false; } + if (row[b] <= i) { return (b % 2) == 0; } + // Assumes that row[a] ≤ i < row[b] + while (a + 1 < b) + { + size_t c = (a + b) / 2; + if (row[c] <= i) + { + a = c; + } + else + { + b = c; + } + } + return (a % 2) == 0; +} + +// Returns the min k such that row[k, i] = 0 (or i+1 if not possible) +int voroffset::CompressedImage::lowerBound(int i, int j) const +{ + const int kmin = std::numeric_limits::min(); + if (j < 0 || j >= (int) m_Rays.size()) + { + return kmin; + } + const std::vector &row = m_Rays[j]; + size_t a = 0; + size_t b = row.size() - 1; + if (row.empty() || row[a] > i) { return kmin; } + if (row[b] <= i) { return (b % 2 == 0 ? i + 1 : row[b]); } + // Assumes that row[a] ≤ i < row[b] + while (a + 1 < b) + { + size_t c = (a + b) / 2; + if (row[c] <= i) + { + a = c; + } + else + { + b = c; + } + } + return (a % 2 == 0 ? i + 1 : row[a]); +} + +// Returns the max k such that row[i, k-1] = 0 (or i if not possible) +int voroffset::CompressedImage::upperBound(int i, int j) const +{ + const int kmax = std::numeric_limits::max(); + if (j < 0 || j >= (int) m_Rays.size()) + { + return kmax; + } + const std::vector &row = m_Rays[j]; + size_t a = 0; + size_t b = row.size() - 1; + if (row.empty() || row[b] <= i) { return kmax; } + if (row[a] > i) { return (a % 2 == 1 ? i : row[a]); } + // Assumes that row[a] ≤ i < row[b] + while (a + 1 < b) + { + size_t c = (a + b) / 2; + if (row[c] <= i) + { + a = c; + } + else + { + b = c; + } + } + return (b % 2 == 1 ? i : row[b]); +} + +// Returns 1 if row j has a non-zero element in [u,v] +bool voroffset::CompressedImage::between(int u, int v, int j) const +{ + if (j < 0 || j >= (int) m_Rays.size()) + { + return false; + } + const std::vector &row = m_Rays[j]; + size_t a = 0; + size_t b = row.size() - 1; + if (row.empty()) { return false; } + if (row[a] > v) { return false; } + if (row[b] <= v) { return (b % 2) == 0 || (row[b] > u); } + // Assumes that row[a] ≤ b < row[b] + while (a + 1 < b) + { + size_t c = (a + b) / 2; + if (row[c] <= v) + { + a = c; + } + else + { + b = c; + } + } + return (a % 2) == 0 || (row[a] > u); +} + +//////////////////////////////////////////////////////////////////////////////// +// Global operations +//////////////////////////////////////////////////////////////////////////////// + +// Negates the map +void voroffset::CompressedImage::negate() +{ + for (auto & row : m_Rays) + { + if (row.empty()) + { + row = {0, m_XSize}; + } + else + { + // First event + if (row[0] == 0) + { + row.erase(row.begin()); + } + else + { + row.insert(row.begin(), 0); + } + // Last event + if (row.back() == m_XSize) + { + row.pop_back(); + } + else + { + row.push_back(m_XSize); + } + } + } +} + +// Transpose operation to get a column-compressed map +voroffset::CompressedImage voroffset::CompressedImage::transposed() const +{ + CompressedImage tr(*this); + tr.transposeInPlace(); + return tr; +} + +void voroffset::CompressedImage::transposeInPlace() +{ + // Store all events and sort them transposed + std::set > allEvents; + for (int j = 0; j < (int) m_Rays.size(); ++j) + { + for (const auto & val : m_Rays[j]) + { + allEvents.emplace(val, j); + } + } + // Resize internal array and swap dims + { + int tmp = m_XSize; + m_XSize = (int) m_Rays.size(); + m_Rays.assign(tmp, {}); + } + + // Helper function: toggle an event in a line + auto toggle = [](std::set &s, int x) + { + if (s.count(x)) + { + s.erase(x); + } + else + { + s.insert(x); + } + }; + +#if 0 + int prevLine = 0; + std::set currentLine; + for (auto it = allEvents.begin(); it != allEvents.end(); ++it) + { + int i = it->first; + int j = it->second; + + // Register current line + if (prevLine != i) + { + for (int ii = prevLine; ii < i; ++ii) + { + for (auto jt = currentLine.begin(); jt != currentLine.end(); ++jt) + { + m_Rays[ii].push_back(*jt); + } + } + } + + // Toggle segment (i, j) -- (i, j + 1) + toggle(currentLine, j); + toggle(currentLine, j + 1); + prevLine = i; + } + // Register last line + m_Rays[prevLine].insert(m_Rays[prevLine].end(), currentLine.begin(), currentLine.end()); +#else + int startCol = -1; + int endCol = -1; + int prevLine = -1; + std::set currentLine; + for (auto it = allEvents.begin(); it != allEvents.end(); ++it) + { + int i = it->first; + int j = it->second; + + // Register current line + if (prevLine != i) + { + // Toggle segment (i, startCol) -- (i, endCol) + toggle(currentLine, startCol); + toggle(currentLine, endCol); + startCol = endCol = j; + + if (prevLine != -1) + { + for (int ii = prevLine; ii < i; ++ii) + { + for (auto jt = currentLine.begin(); jt != currentLine.end(); ++jt) + { + m_Rays[ii].push_back(*jt); + } + } + } + } + + if (j != endCol) + { + // Toggle segment (i, startCol) -- (i, endCol) + toggle(currentLine, startCol); + toggle(currentLine, endCol); + startCol = j; + } + endCol = j + 1; + prevLine = i; + } + + // Toggle segment (i, startCol) -- (i, endCol) + toggle(currentLine, startCol); + toggle(currentLine, endCol); + + // Register last line + m_Rays[prevLine].insert(m_Rays[prevLine].end(), currentLine.begin(), currentLine.end()); +#endif +} + +//////////////////////////////////////////////////////////////////////////////// +// Morphology operators +//////////////////////////////////////////////////////////////////////////////// + +namespace +{ + typedef std::vector IntVector; + + template + void addSegmentAtTheEnd(Container & dst, int start, int end) + { + if (dst.empty() || (dst.back() < start)) + { + dst.push_back(start); + dst.push_back(end); + } + else if (end > dst.back()) + { + dst.back() = end; + } + } + + void getNextSegment(const IntVector & segs, IntVector::const_iterator & it, + int & left, int & right, bool & full) + { + if (it == segs.end()) + { + full = false; + } + else + { + left = *it++; + right = *it++; + } + } + + // Computes the union of two sorted lists of segments. + // Uses a parallel sweep of both lists, akin to merge-sort. + void unionSegs(const IntVector & a, const IntVector & b, IntVector & result) + { + if (a.empty()) + { + result.insert(result.end(), b.begin(), b.end()); + return; + } + if (b.empty()) + { + result.insert(result.end(), a.begin(), a.end()); + return; + } + IntVector::const_iterator ia(a.begin()), ib(b.begin()); + int al, ar, bl, br; + al = *ia++; ar = *ia++; + bl = *ib++; br = *ib++; + bool bfull(true), afull(true); + while (bfull || afull) + { + if (bfull && bl <= al) + { + addSegmentAtTheEnd(result, bl, br); + getNextSegment(b, ib, bl, br, bfull); + } + else if (afull) + { + addSegmentAtTheEnd(result, al, ar); + if (ia == a.end()) + { + afull = false; + al = std::numeric_limits::max(); + } + else + { + al = *ia++; + ar = *ia++; + } + } + } + } + + template + void unionMap(Iterator1 a, Iterator2 b, std::vector > &res) + { + for (size_t i = 0; i < res.size(); ++i) + { + IntVector temp; + unionSegs(a[i], b[i], temp); + res[i].swap(temp); + } + } + +} // anonymous namespace + +//////////////////////////////////////////////////////////////////////////////// + +void CompressedImage::dilate(double radius) +{ + typedef std::reverse_iterator rev_t; + + std::vector > r1, r2; + voronoi_half_dilate(height(), width(), radius, m_Rays.begin(), r1); + voronoi_half_dilate(height(), width(), radius, rev_t(m_Rays.end()), r2); + unionMap(r1.cbegin(), r2.crbegin(), m_Rays); + vor_assert(isValid() == true); +} + +// ----------------------------------------------------------------------------- + +void CompressedImage::erode(double radius) +{ + typedef std::reverse_iterator rev_t; + + std::vector > r1, r2; + voronoi_half_erode(height(), width(), radius, m_Rays.begin(), r1); + voronoi_half_erode(height(), width(), radius, rev_t(m_Rays.end()), r2); + unionMap(r1.cbegin(), r2.crbegin(), m_Rays); + negate(); + vor_assert(isValid()); +} + +// ----------------------------------------------------------------------------- + +void CompressedImage::close(double r) +{ + dilate(r); + erode(r); +} + +// ----------------------------------------------------------------------------- + +void CompressedImage::open(double r) +{ + erode(r); + dilate(r); +} diff --git a/src/vor2d/CompressedImage.h b/src/vor2d/CompressedImage.h new file mode 100644 index 0000000..6e45796 --- /dev/null +++ b/src/vor2d/CompressedImage.h @@ -0,0 +1,105 @@ +#pragma once + +//////////////////////////////////////////////////////////////////////////////// +#include "vor2d/Common.h" +#include "vor2d/Image.h" +#include +#include +#include +#include +//////////////////////////////////////////////////////////////////////////////// + +namespace voroffset +{ + +class CompressedImage +{ + typedef int Scalar; + +private: + ////////////////// + // Data members // + ////////////////// + + int m_XSize; + std::vector> m_Rays; + +public: + ////////////////// + // Input/Output // + ////////////////// + + CompressedImage(int w = 0, int h = 0) : m_XSize(w), m_Rays(h) { } + CompressedImage(std::istream &in) { this->load(in); } + + // Import/export + void fromImage(std::function read_pixel_func); + void toImage(std::function write_pixel_func) const; + + // Marshalling + void save(std::ostream &out) const; + void load(std::istream &in); + +public: + //////////// + // Access // + //////////// + + // Accessors + int width() const { return m_XSize; } + int height() const { return (int) m_Rays.size(); } + + // Resize + void resize(int w, int h); + + // Sanity check + bool isValid() const; + + // Empty? + bool empty() const; + + // Apply a function to each segment in the structure + void iterate(std::function func) const; + + // Writes a pixel at a given location + // return true if the value was changed + bool write(int i, int j, bool newVal); + + // Dichotomy search + bool at(int i, int j) const; + + // Returns the min k such that x[k, i] = 0 (or i+1 if not possible) + int lowerBound(int i, int j) const; + + // Returns the max k such that x[i, k-1] = 0 (or i if not possible) + int upperBound(int i, int j) const; + + // Returns 1 if column j has a non-zero element between u and v (inclusive) + bool between(int u, int v, int j) const; + +public: + /////////////////////// + // Global operations // + /////////////////////// + + // Negates the map + void negate(); + + // Transpose operation to get a column-compressed map + CompressedImage transposed() const; + + // Transpose operation to get a column-compressed map + void transposeInPlace(); + +public: + ////////////////////////// + // Morphology operators // + ////////////////////////// + + void dilate(double r); + void erode(double r); + void close(double r); + void open(double r); +}; + +} // namespace voroffset diff --git a/src/vor2d/Dexelize.cpp b/src/vor2d/Dexelize.cpp new file mode 100644 index 0000000..043234e --- /dev/null +++ b/src/vor2d/Dexelize.cpp @@ -0,0 +1,89 @@ +#include "Dexelize.h" +#define NANOSVG_IMPLEMENTATION +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#ifndef DPI +#define DPI 90 +#endif // ! DPI + + +voroffset::DoubleCompressedImage voroffset::create_dexels(const std::string &file) +{ + std::vector contours; + std::cerr << "[loading] " << file << " ... "; + double tol = 1.5; + NSVGimage* image = nsvgParseFromFile(file.c_str(), "mm", DPI); + // For all curves in the file + for (NSVGshape *shape = image->shapes; shape != NULL; shape = shape->next) + { + // Flatten the curve in a series of segments + Curve poly; + for (NSVGpath *path = shape->paths; path != NULL; path = path->next) + { + for (int i = 0; i < path->npts - 1; i++) + poly.push_back(PointF(path->pts[2 * i], path->pts[2 * i + 1])); + } + contours.push_back(poly); + } + std::cerr << "Read a SVG of size : " << image->width << " x " << image->height + << " (" << contours.size() << " polygones)" << std::endl; + DoubleCompressedImage dexels(std::ceil(image->height*25.4 / DPI), std::ceil(image->width*25.4 / DPI)); + dexels.fromImage(contours); + nsvgDelete(image); + return dexels; +} + +void voroffset::dexel_dump(const std::string &filename, const DoubleCompressedImage &dexels) +{ + // Initialize the Geogram library + GEO::initialize(); + + // Import standard command line arguments, and custom ones + GEO::CmdLine::import_arg_group("standard"); + + GEO::Mesh mesh; + + for (int x = 0; x < dexels.height(); ++x) + { + for (int i = 0; 2 * i < dexels.at(x).size(); ++i) + { + GEO::vec3 xyz_min, xyz_max; + xyz_min[0] = x; + xyz_min[1] = dexels.at(x)[2 * i + 0]; + xyz_min[2] = 1; + xyz_max[0] = x + 1; + xyz_max[1] = dexels.at(x)[2 * i + 1]; + xyz_max[2] = 1; + GEO::vec3 diff[4] = + { + GEO::vec3(0,0,0), GEO::vec3(1,0,0), GEO::vec3(0,1,0), GEO::vec3(1,1,0), + }; + int v = mesh.vertices.nb(); + for (int lv = 0; lv < 4; ++lv) + { + for (int d = 0; d < 3; ++d) + { + diff[lv][d] = xyz_min[d] + diff[lv][d] * (xyz_max[d] - xyz_min[d]); + } + mesh.vertices.create_vertex(&diff[lv][0]); + } + mesh.facets.create_quad(v, v + 1, v + 2, v + 3); + } + } + mesh.facets.compute_borders(); + mesh.facets.connect(); + mesh.vertices.remove_isolated(); + GEO::mesh_save(mesh, filename); +} diff --git a/src/vor2d/Dexelize.h b/src/vor2d/Dexelize.h new file mode 100644 index 0000000..45289d9 --- /dev/null +++ b/src/vor2d/Dexelize.h @@ -0,0 +1,18 @@ +#pragma once + +//////////////////////////////////////////////////////////////////////////////// +#include "vor2d/Common.h" +#include "vor2d/DoubleCompressedImage.h" +#include +#include +#include +#include +//////////////////////////////////////////////////////////////////////////////// + +namespace voroffset +{ + DoubleCompressedImage create_dexels(const std::string &filename); + + void dexel_dump(const std::string &filename, const DoubleCompressedImage &dexels); + +} // namespace voroffset3d diff --git a/src/vor2d/DistanceTransform.cpp b/src/vor2d/DistanceTransform.cpp new file mode 100644 index 0000000..675ffbb --- /dev/null +++ b/src/vor2d/DistanceTransform.cpp @@ -0,0 +1,177 @@ +//////////////////////////////////////////////////////////////////////////////// +#include "vor2d/DistanceTransform.h" +#include +#include +#include +//////////////////////////////////////////////////////////////////////////////// + +namespace +{ + + constexpr int DIST_MAX = std::numeric_limits::max(); + + // v++ avoiding overflows + int safe_incr(int v) { return (v == DIST_MAX ? v : v + 1); } + +} // anonymous namespace + +//////////////////////////////////////////////////////////////////////////////// + +// Meijster, A., Roerdink, J. B., & Hesselink, W. H. (2002). +// A general algorithm for computing distance transforms in linear time. +void voroffset::computeDistanceTransformExact(const Image &img, Image &dt) +{ + const int width = (int) img.width(); + const int height = (int) img.height(); + + // Initialization + dt.resize(width, height); + + ///////////// + // PHASE 1 // + ///////////// + + #pragma omp parallel for + for (int x = 0; x < width; ++x) + { + // Scan 1 + dt(x, 0) = (img(x, 0) == 0 ? 0 : DIST_MAX); + for (int y = 1; y < height; ++y) + { + dt(x, y) = (img(x, y) == 0 ? 0 : safe_incr(dt(x, y - 1))); + } + // Scan 2 + for (int y = height - 2; y >= 0; --y) + { + if (dt(x, y + 1) < dt(x,y)) + { + dt(x, y) = 1 + dt(x, y + 1); + } + } + } + + ///////////// + // PHASE 2 // + ///////////// + + #pragma omp parallel for + for (int y = 0; y < height; ++y) + { + Eigen::VectorXi s(width); + Eigen::VectorXi t(width); + Eigen::VectorXi r(width); // Temporary row + + // See [Meijster et al. 2002] for explanations of 'f' and 'sep' + auto f = [y, &dt] (int x, int i) + { + return (dt(i, y) == DIST_MAX ? DIST_MAX : (x - i)*(x - i) + dt(i, y)); + }; + auto sep = [y, &dt] (int i, int u) + { + if (dt(u, y) == DIST_MAX || dt(i, y) == DIST_MAX) + { + return DIST_MAX; + } + else + { + return (u*u - i*i + dt(u, y) - dt(i, y)) / (2 * (u - i)); + } + }; + + // Scan 3 + int q = 0; s[0] = 0; t[0] = 0; + for (int u = 1; u < width; ++u) + { + vor_assert(q < width); + while (q >= 0 && f(t[q], s[q]) > f(t[q], u)) + { + --q; + } + if (q < 0) + { + q = 0; s[0] = u; + } + else + { + int w = safe_incr(sep(s[q], u)); + if (w < width) + { + q = q + 1; s[q] = u; t[q] = w; + } + } + } + + // Scan 4 + for (int u = width - 1; u >= 0; --u) + { + r[u] = f(u, s[q]); + if (u == t[q]) { --q; } + } + dt.data().row(y) = r; + } +} + +//////////////////////////////////////////////////////////////////////////////// + +// Danielsson, P. E. (1980). Euclidean distance mapping. +void voroffset::computeDistanceTransformApprox(const Image &img, Image &dist) +{ + const int xmax = (int) dist.width(); + const int ymax = (int) dist.height(); + + // Initialization + dist.resize(xmax, ymax); + for (int y = 0; y < ymax; ++y) + { + for (int x = 0; x < xmax; ++x) + { + dist(x, y) = (img(x, y) == 0 ? Vector2i(0, 0) : Vector2i(xmax, ymax)); + } + } + + //////////// + // SCAN 1 // + //////////// + + for (int y = 1; y < ymax; ++y) + { + for (int x = 0; x < xmax; ++x) + { + Vector2i vt = dist(x, y-1) + Vector2i(0, -1); + if (vt.squaredNorm() < dist(x, y).squaredNorm()) { dist(x, y) = vt; } + } + for (int x = 1; x < xmax; ++x) + { + Vector2i vt = dist(x-1, y) + Vector2i(-1, 0); + if (vt.squaredNorm() < dist(x, y).squaredNorm()) { dist(x, y) = vt; } + } + for (int x = (int) xmax - 2; x >= 0; --x) + { + Vector2i vt = dist(x+1, y) + Vector2i(1, 0); + if (vt.squaredNorm() < dist(x, y).squaredNorm()) { dist(x, y) = vt; } + } + } + + //////////// + // SCAN 2 // + //////////// + + for (int y = ymax - 2; y >= 0; --y) + { + for (int x = 0; x < xmax; ++x) + { + Vector2i vt = dist(x, y+1) + Vector2i(0, 1); + if (vt.squaredNorm() < dist(x, y).squaredNorm()) { dist(x, y) = vt; } + } + for (int x = 1; x < xmax; ++x) + { + Vector2i vt = dist(x-1, y) + Vector2i(-1, 0); + if (vt.squaredNorm() < dist(x, y).squaredNorm()) { dist(x, y) = vt; } + } + for (int x = xmax - 2; x >= 0; --x) + { + Vector2i vt = dist(x+1, y) + Vector2i(1, 0); + if (vt.squaredNorm() < dist(x, y).squaredNorm()) { dist(x, y) = vt; } + } + } +} diff --git a/src/vor2d/DistanceTransform.h b/src/vor2d/DistanceTransform.h new file mode 100644 index 0000000..362d768 --- /dev/null +++ b/src/vor2d/DistanceTransform.h @@ -0,0 +1,44 @@ +#pragma once + +//////////////////////////////////////////////////////////////////////////////// +#include "vor2d/Common.h" +#include "vor2d/Image.h" +#include +#include +//////////////////////////////////////////////////////////////////////////////// + +namespace voroffset +{ + +// +// @brief Compute the distance transform of an image. +// +// @param[in] img { Input image, where 0 represents a feature. } +// @param[out] dt { Distance to the nearest feature in the image. } +// +// The input image must contain 0 where the pixels are on the object, +// and high values (int max works) for pixels on the outside. +// +// Reference: +// +// Meijster, A., Roerdink, J. and Hesselink, W. 2000. +// A general algorithm for computing distance transforms in linear time. +// In Proceedings of the 5th International Conference on Mathematical +// Morphology and its Applications to Image and Signal Processing. +// +void computeDistanceTransformExact(const Image &img, Image &dt); + +// +// @brief Compute Euclidean distance map. +// +// @param[in] img { Input image, where 0 represents a feature. } +// @param[out] dist { 2D array of vectors corresponding to the distance to the +// nearest feature. } +// Reference: +// +// Danielsson, Per-Erik. Euclidean distance mapping. +// Computer Graphics and image processing 14.3 (1980): 227-248. +// +void computeDistanceTransformApprox(const Image &img, Image &dist); + +} // namespace voroffset diff --git a/src/vor2d/DoubleCompressedImage.cpp b/src/vor2d/DoubleCompressedImage.cpp new file mode 100644 index 0000000..ff8c18f --- /dev/null +++ b/src/vor2d/DoubleCompressedImage.cpp @@ -0,0 +1,772 @@ +//////////////////////////////////////////////////////////////////////////////// +#include "vor2d/DoubleCompressedImage.h" +#include "vor2d/DoubleVoronoi.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +//////////////////////////////////////////////////////////////////////////////// + +using namespace voroffset; +#ifndef MIN_SEG_SIZE +#define MIN_SEG_SIZE 1e-10 +#endif + +//////////////////////////////////////////////////////////////////////////////// +// Input/Output +//////////////////////////////////////////////////////////////////////////////// + +// Import from 2D image +void voroffset::DoubleCompressedImage::fromImage(std::vector input_curves) +{ + int size_ = input_curves.size(); + std::vector intersections_; + for (int j = 0; j < (int)m_Rays.size(); ++j) + { + m_Rays[j].clear(); + for (int k = 0; k < size_; k++) + { + scanLine(intersections_,j,input_curves[k]); + unionIntersections(m_Rays[j], intersections_); + intersections_.clear(); + } + } +} + + +void voroffset::DoubleCompressedImage::scanLine(std::vector &intersections, int line_x, Curve curve) +// Scanline algorithm for the line[i] +{ + intersections.clear(); + int i, j, k; + Scalar s; + Scalar y_in_; + bool status; + int size_ = curve.size(); + for (i = 0; i line_x) + { + s = (curve[(i + 1) % size_].real() - line_x) / ( curve[(i + 1) % size_].real() - curve[i % size_].real() ); + y_in_ = s*curve[i].imag() + (1 - s)*curve[(i + 1) % size_].imag(); + intersections.push_back(y_in_); + } + else + { + j = 1; + while (curve[(i + j) % size_].real() == line_x) + j++; + if (curve[(i + j) % size_].real() > line_x) + intersections.push_back(curve[(i + j - 1) % size_].imag()); + } + } + else if (curve[i].real() > line_x) + { + if (curve[(i + 1) % size_].real() < line_x) + { + s = (curve[(i + 1) % size_].real() - line_x) / (curve[(i + 1) % size_].real() - curve[i].real()); + y_in_ = s*curve[i].imag() + (1 - s)*curve[(i + 1) % size_].imag(); + intersections.push_back(y_in_); + } + else + { + j = 1; + while (curve[(i + j) % size_].real() == line_x) + j++; + if (curve[(i + j) % size_].real() < line_x) + intersections.push_back(curve[(i + j - 1) % size_].imag()); + } + } + + sort(intersections.begin(), intersections.end()); + +} +void voroffset::DoubleCompressedImage::unionIntersections(std::vector &intersections_1, std::vector intersections_2) +// Due to we assume two closed curves have no intersections, so the imageinary part of the this two intersections have no overlaps +{ + if (intersections_2.size() == 0) + return; + if (intersections_1.size()==0) + { + for (int i = 0; i < intersections_2.size(); i++) + intersections_1.push_back(intersections_2[i]); + } + else if(intersections_1[intersections_1.size()-1]<=intersections_2[0]) // intersections_2 is totally on the right side of the intersections_1 + { + for (int i = 0; i < intersections_2.size(); i++) + intersections_1.push_back(intersections_2[i]); + } + else // intersections_2 is totally on the left side of the intersections_1 + { + std::vector::iterator it_ = intersections_1.begin(); + for (int i = intersections_2.size() - 1; i >= 0; i--) + it_ = intersections_1.insert(it_, intersections_2[i]); + } +} + +void voroffset::DoubleCompressedImage::copyFrom(const voroffset::DoubleCompressedImage source_dexel) +{ + m_XSize = source_dexel.m_XSize; + int size_ = source_dexel.m_Rays.size(); + m_Rays.assign(size_, std::vector()); + for (int i = 0; i < size_; i++) + { + int inter_num = source_dexel.m_Rays[i].size(); + vor_assert(inter_num % 2 == 0); + m_Rays[i].resize(inter_num); + std::copy_n(source_dexel.m_Rays[i].begin(), inter_num, m_Rays[i].begin()); + } +} + +// Export to 2D image +/*void voroffset::CompressedImageFloat::toImage(std::function write_pixel_func) const { + for (int j = 0; j < (int)m_Rays.size(); ++j) { + int i = 0; + bool status = false; + for (const auto & val : m_Rays[j]) { + for (; i < val; ++i) { + write_pixel_func(i, j, status); + } + status = !status; + } + for (; i < m_XSize; ++i) { + write_pixel_func(i, j, status); + } + } +} +*/ +// Marshalling +void voroffset::DoubleCompressedImage::save(std::ostream &out) const +{ + out << m_XSize << " " << m_Rays.size() << "\n"; + size_t num_rows = m_Rays.size(); + for (int i = 0; i < num_rows; i++) + { + size_t size=m_Rays[i].size(); + out << size; + vor_assert(size % 2 == 0); + for (int k = 0; k < size; k+=2) + { + out << ' ' << m_Rays[i][k]; + out << ' ' << m_Rays[i][k + 1]; + } + out << '\n'; + + } +} + +void voroffset::DoubleCompressedImage::load(std::istream &in) +{ + unsigned int num_cols; + in >> m_XSize >> num_cols; + m_Rays.assign(num_cols, {}); + int flag; + for (int i = 0; i < num_cols; i++) + { + size_t size; + in >> size; + vor_assert(size % 2 == 0); + m_Rays[i].resize(size); + for (int k = 0; k < size; k += 2) + { + in >> m_Rays[i][k]>>m_Rays[i][k+1]; + } + + + } +} + +//////////////////////////////////////////////////////////////////////////////// +// Access +//////////////////////////////////////////////////////////////////////////////// + +// Resize +void voroffset::DoubleCompressedImage::resize(int w, int h) +{ + m_XSize = w; + m_Rays.assign(h, {}); +} + +// Sanity check +bool voroffset::DoubleCompressedImage::isValid() const +{ + for (const auto & row : m_Rays) + { + if (row.size() % 2 != 0) + { + std::cerr << "row size % 2 != 0" << std::endl; + return false; + } + Scalar lastEvent = -1; + for (const auto &val : row) + { + if (val < lastEvent) + { + std::cerr << "val > lastEvent" << std::endl; + return false; + } + lastEvent = val; + } + if (lastEvent > m_XSize) + { + std::cerr << "lastEvent <= xsize" << std::endl; + return false; + } + } + return true; +} + +// Empty? +bool voroffset::DoubleCompressedImage::empty() const +{ + for (const auto & row : m_Rays) + { + if (!row.empty()) { return false; } + } + return true; +} + +// Apply a function to each segment in the structure +void voroffset::DoubleCompressedImage::iterate(std::function func) const +{ + for (int i = 0; i < (int)m_Rays.size(); ++i) + { + if (m_Rays[i].size() % 2 != 0) + { + throw std::runtime_error("Invalid number of events along a ray."); + } + for (size_t j = 0; 2 * j < m_Rays[i].size(); ++j) + { + func(i, m_Rays[i][2 * j], m_Rays[i][2 * j + 1]); + } + } +} + +// Writes a pixel at a given location +// return true if the value was changed +bool voroffset::DoubleCompressedImage::write(int i, int j, bool newVal) +{ +/* if (j < 0 || j >= (int)m_Rays.size()) { + return false; + } + std::vector &row = m_Rays[j]; + if (row.empty()) { + // poke in an empty line + if (!newVal) { return false; } + row = { i, i + 1 }; + return true; + } + if (row[0] > i) { + // poke to the left of the segments + if (!newVal) { return false; } + if (row[0] > i + 1) { + row.insert(row.begin(), 2, i + 1); + } + row[0] = i; + return true; + } + size_t a = 0; + size_t b = row.size() - 1; + if (row[b] <= i) { + // poke to the right of the segments + if (!newVal) { return false; } + if (row[b] == i) { + row[b] = i + 1; + } + else { + row.push_back(i); + row.push_back(i + 1); + } + return true; + } + // Assumes that row[a] �� i < row[b] + while (a + 1 < b) { + size_t c = (a + b) / 2; + if (row[c] <= i) { + a = c; + } + else { + b = c; + } + } + bool isInside = ((a % 2) == 0); + if (isInside == newVal) { return false; } + int len = row[a + 1] - row[a]; + if (len == 1) { // erase all + row.erase(row.begin() + a, row.begin() + a + 2); + } + else if (i == row[a]) { // move left end to the right + row[a] = i + 1; + } + else if (i + 1 == row[a + 1]) { // move right end to the left + row[a + 1] = i; + } + else { + row.insert(row.begin() + a + 1, 2, i + 1); + row[a + 1] = i; + } + return true;*/ + return true; +} + +// Dichotomy search +bool voroffset::DoubleCompressedImage::at(int i, int j) const +{ + if (j < 0 || j >= (int)m_Rays.size()) + { + return false; + } + const std::vector &row = m_Rays[j]; + size_t a = 0; + size_t b = row.size() - 1; + if (row.empty() || row[a] > i) { return false; } + if (row[b] <= i) { return (b % 2) == 0; } + // Assumes that row[a] �� i < row[b] + while (a + 1 < b) + { + size_t c = (a + b) / 2; + if (row[c] <= i) + { + a = c; + } + else + { + b = c; + } + } + return (a % 2) == 0; +} + +// Returns the min k such that row[k, i] = 0 (or i+1 if not possible) +int voroffset::DoubleCompressedImage::lowerBound(int i, int j) const +{ + const int kmin = std::numeric_limits::min(); + if (j < 0 || j >= (int)m_Rays.size()) + { + return kmin; + } + const std::vector &row = m_Rays[j]; + size_t a = 0; + size_t b = row.size() - 1; + if (row.empty() || row[a] > i) { return kmin; } + if (row[b] <= i) { return (b % 2 == 0 ? i + 1 : row[b]); } + // Assumes that row[a] �� i < row[b] + while (a + 1 < b) + { + size_t c = (a + b) / 2; + if (row[c] <= i) + { + a = c; + } + else + { + b = c; + } + } + return (a % 2 == 0 ? i + 1 : row[a]); +} + +// Returns the max k such that row[i, k-1] = 0 (or i if not possible) +int voroffset::DoubleCompressedImage::upperBound(int i, int j) const +{ + const int kmax = std::numeric_limits::max(); + if (j < 0 || j >= (int)m_Rays.size()) + { + return kmax; + } + const std::vector &row = m_Rays[j]; + size_t a = 0; + size_t b = row.size() - 1; + if (row.empty() || row[b] <= i) { return kmax; } + if (row[a] > i) { return (a % 2 == 1 ? i : row[a]); } + // Assumes that row[a] �� i < row[b] + while (a + 1 < b) + { + size_t c = (a + b) / 2; + if (row[c] <= i) + { + a = c; + } + else + { + b = c; + } + } + return (b % 2 == 1 ? i : row[b]); +} + +// Returns 1 if row j has a non-zero element in [u,v] +bool voroffset::DoubleCompressedImage::between(int u, int v, int j) const +{ + if (j < 0 || j >= (int)m_Rays.size()) + { + return false; + } + const std::vector &row = m_Rays[j]; + size_t a = 0; + size_t b = row.size() - 1; + if (row.empty()) { return false; } + if (row[a] > v) { return false; } + if (row[b] <= v) { return (b % 2) == 0 || (row[b] > u); } + // Assumes that row[a] �� b < row[b] + while (a + 1 < b) + { + size_t c = (a + b) / 2; + if (row[c] <= v) + { + a = c; + } + else + { + b = c; + } + } + return (a % 2) == 0 || (row[a] > u); +} + +//////////////////////////////////////////////////////////////////////////////// +// Global operations +//////////////////////////////////////////////////////////////////////////////// + +// Negates the map +void voroffset::DoubleCompressedImage::negate() +{ + for (auto & row : m_Rays) + { + if (row.empty()) + { + row = { 0, m_XSize*1.0f }; + } + else + { + // First event + if (row[0] == 0) + { + row.erase(row.begin()); + } + else + { + row.insert(row.begin(), 0); + } + // Last event + if (row.back() == m_XSize) + { + row.pop_back(); + } + else + { + row.push_back(m_XSize); + } + } + } +} + +// Transpose operation to get a column-compressed map +voroffset::DoubleCompressedImage voroffset::DoubleCompressedImage::transposed() const +{ + DoubleCompressedImage tr(*this); + tr.transposeInPlace(); + return tr; +} + +void voroffset::DoubleCompressedImage::transposeInPlace() +{ + // Store all events and sort them transposed + std::set > allEvents; + for (int j = 0; j < (int)m_Rays.size(); ++j) + { + for (const auto & val : m_Rays[j]) + { + allEvents.emplace(val, j); + } + } + // Resize internal array and swap dims + { + int tmp = m_XSize; + m_XSize = (int)m_Rays.size(); + m_Rays.assign(tmp, {}); + } + + // Helper function: toggle an event in a line + auto toggle = [](std::set &s, int x) + { + if (s.count(x)) + { + s.erase(x); + } + else + { + s.insert(x); + } + }; + +#if 0 + int prevLine = 0; + std::set currentLine; + for (auto it = allEvents.begin(); it != allEvents.end(); ++it) + { + int i = it->first; + int j = it->second; + + // Register current line + if (prevLine != i) + { + for (int ii = prevLine; ii < i; ++ii) + { + for (auto jt = currentLine.begin(); jt != currentLine.end(); ++jt) + { + m_Rays[ii].push_back(*jt); + } + } + } + + // Toggle segment (i, j) -- (i, j + 1) + toggle(currentLine, j); + toggle(currentLine, j + 1); + prevLine = i; + } + // Register last line + m_Rays[prevLine].insert(m_Rays[prevLine].end(), currentLine.begin(), currentLine.end()); +#else + int startCol = -1; + int endCol = -1; + int prevLine = -1; + std::set currentLine; + for (auto it = allEvents.begin(); it != allEvents.end(); ++it) + { + int i = it->first; + int j = it->second; + + // Register current line + if (prevLine != i) + { + // Toggle segment (i, startCol) -- (i, endCol) + toggle(currentLine, startCol); + toggle(currentLine, endCol); + startCol = endCol = j; + + if (prevLine != -1) + { + for (int ii = prevLine; ii < i; ++ii) + { + for (auto jt = currentLine.begin(); jt != currentLine.end(); ++jt) + { + m_Rays[ii].push_back(*jt); + } + } + } + } + + if (j != endCol) + { + // Toggle segment (i, startCol) -- (i, endCol) + toggle(currentLine, startCol); + toggle(currentLine, endCol); + startCol = j; + } + endCol = j + 1; + prevLine = i; + } + + // Toggle segment (i, startCol) -- (i, endCol) + toggle(currentLine, startCol); + toggle(currentLine, endCol); + + // Register last line + m_Rays[prevLine].insert(m_Rays[prevLine].end(), currentLine.begin(), currentLine.end()); +#endif +} + +//////////////////////////////////////////////////////////////////////////////// +// Morphology operators +//////////////////////////////////////////////////////////////////////////////// + +namespace +{ + typedef std::vector DoubleVector; + + template + void addSegmentAtTheEnd(Container & dst, double start, double end) + { + if (dst.empty() || (dst.back() < start)) + { + dst.push_back(start); + dst.push_back(end); + } + else if (end > dst.back()) + { + dst.back() = end; + } + } + + void getNextSegment(const DoubleVector & segs, DoubleVector::const_iterator & it, + double & left, double & right, bool & full) + { + if (it == segs.end()) + { + full = false; + } + else + { + left = *it++; + right = *it++; + } + } + + // Computes the union of two sorted lists of segments. + // Uses a parallel sweep of both lists, akin to merge-sort. + void unionSegs(const DoubleVector & a, const DoubleVector & b, DoubleVector & result) + { + if (a.empty()) + { + result.insert(result.end(), b.begin(), b.end()); + return; + } + if (b.empty()) + { + result.insert(result.end(), a.begin(), a.end()); + return; + } + DoubleVector::const_iterator ia(a.begin()), ib(b.begin()); + double al, ar, bl, br; + al = *ia++; ar = *ia++; + bl = *ib++; br = *ib++; + bool bfull(true), afull(true); + while (bfull || afull) + { + if (bfull && bl <= al) + { + addSegmentAtTheEnd(result, bl, br); + getNextSegment(b, ib, bl, br, bfull); + } + else if (afull) + { + addSegmentAtTheEnd(result, al, ar); + if (ia == a.end()) + { + afull = false; + al = std::numeric_limits::max(); + } + else + { + al = *ia++; + ar = *ia++; + } + } + } + } + + template + void unionMap(Iterator1 a, Iterator2 b, std::vector > &res) + { + for (size_t i = 0; i < res.size(); ++i) + { + DoubleVector temp; + unionSegs(a[i], b[i], temp); + res[i].swap(temp); + } + } + +} // anonymous namespace + + //////////////////////////////////////////////////////////////////////////////// + +void DoubleCompressedImage::dilate(double radius) +{ + typedef std::reverse_iterator rev_t; + + std::vector > r1, r2; + voronoiF_half_dilate(height(), width(), radius * m_Rays.size(), m_Rays.begin(), r1); + voronoiF_half_dilate(height(), width(), radius * m_Rays.size(), rev_t(m_Rays.end()), r2); + unionMap(r1.cbegin(), r2.crbegin(), m_Rays); + vor_assert(isValid() == true); +} + +// ----------------------------------------------------------------------------- + +void DoubleCompressedImage::erode(double radius) +{ + typedef std::reverse_iterator rev_t; + + std::vector > r1, r2; + voronoiF_half_erode(height(), width(), radius, m_Rays.begin(), r1); + voronoiF_half_erode(height(), width(), radius, rev_t(m_Rays.end()), r2); + unionMap(r1.cbegin(), r2.crbegin(), m_Rays); + negate(); + vor_assert(isValid()); +} + +// ----------------------------------------------------------------------------- + +void DoubleCompressedImage::close(double r) +{ + dilate(r); + erode(r); +} + +// ----------------------------------------------------------------------------- + +void DoubleCompressedImage::open(double r) +{ + erode(r); + dilate(r); +} + +void DoubleCompressedImage::negateRay(std::vector &result, double y_min, double y_max) +{ + int size_ = result.size(); + + if (result.empty()) + { + result = { y_min, y_max }; + } + else + { + // First event + if (result[0] == y_min) + { + result.erase(result.begin()); + } + else + { + result.insert(result.begin(), y_min); + } + // Last event + if (result.back() == y_max) + { + result.pop_back(); + } + else + { + result.push_back(y_max); + } + } +} + +void DoubleCompressedImage::removePoint(const std::vector segs, std::vector &res) +{ + int size_ = segs.size(); + for (int i = 0; i < size_; i += 2) + { + if (segs[i + 1] - segs[i] < MIN_SEG_SIZE) + continue; + res.push_back(segs[i]); + res.push_back(segs[i + 1]); + + } +} + +double DoubleCompressedImage::area() +{ + double area = 0.f; + for (int i = 0; i < m_Rays.size(); i++) + for (int k = 0; k < m_Rays[i].size(); k+=2) + area += 1.0*(m_Rays[i][k + 1] - m_Rays[i][k]); + return area; +} \ No newline at end of file diff --git a/src/vor2d/DoubleCompressedImage.h b/src/vor2d/DoubleCompressedImage.h new file mode 100644 index 0000000..183aad3 --- /dev/null +++ b/src/vor2d/DoubleCompressedImage.h @@ -0,0 +1,117 @@ +#pragma once + +//////////////////////////////////////////////////////////////////////////////// +#include "vor2d/Common.h" +#include "vor2d/Image.h" +#include +#include +#include +#include +//////////////////////////////////////////////////////////////////////////////// + +//////////////////////////////////////////////////////////////////////////////// +namespace voroffset +{ + +class DoubleCompressedImage +{ + + +private: + ////////////////// + // Data members // + ////////////////// + int m_XSize; +public: + std::vector> m_Rays; + std::vector> m_Segs; + +public: + ////////////////// + // Input/Output // + ////////////////// + + DoubleCompressedImage(int w = 0, int h = 0) : m_XSize(w), m_Rays(h) { } + DoubleCompressedImage(std::istream &in) { this->load(in); } + + + // Import/export + void fromImage(std::vectorinput_curves); + // Marshalling + void save(std::ostream &out) const; + void load(std::istream &in); + +private: + void scanLine(std::vector &intersections, int i, Curve curve); // Scanline algorithm for the line[i] + void unionIntersections(std::vector &intersections_1, std::vector intersections_2); + // Union intersections, the results store in the first one +public: + //////////// + // Access // + //////////// + // Accessors + int width() const { return m_XSize; } + int height() const { return (int) m_Rays.size(); } + + // Resize + void resize(int w, int h); + + // Sanity check + bool isValid() const; + + // Empty? + bool empty() const; + + // Apply a function to each segment in the structure + void iterate(std::function func) const; + + const std::vector & at(int x) const { return m_Rays[x]; } + + // Writes a pixel at a given location + // return true if the value was changed + bool write(int i, int j, bool newVal); + + // Dichotomy search + bool at(int i, int j) const; + + // Returns the min k such that x[k, i] = 0 (or i+1 if not possible) + int lowerBound(int i, int j) const; + + // Returns the max k such that x[i, k-1] = 0 (or i if not possible) + int upperBound(int i, int j) const; + + // Returns 1 if column j has a non-zero element between u and v (inclusive) + bool between(int u, int v, int j) const; + +public: + /////////////////////// + // Global operations // + /////////////////////// + + // Negates the map + void negate(); + + // Copy from another dexel + void copyFrom(const DoubleCompressedImage source_dexel); + + // Transpose operation to get a column-compressed map + DoubleCompressedImage transposed() const; + + // Transpose operation to get a column-compressed map + void transposeInPlace(); + +public: + ////////////////////////// + // Morphology operators // + ////////////////////////// + + void dilate(double r); + void erode(double r); + void close(double r); + void open(double r); + void negateRay(std::vector &result, double y_min, double y_max); + double area(); + void removePoint(const std::vector segs, std::vector &res); +}; + +} // namespace voroffset diff --git a/src/vor2d/DoubleVoronoi.cpp b/src/vor2d/DoubleVoronoi.cpp new file mode 100644 index 0000000..dc2de01 --- /dev/null +++ b/src/vor2d/DoubleVoronoi.cpp @@ -0,0 +1,731 @@ +//////////////////////////////////////////////////////////////////////////////// +#include "vor2d/DoubleVoronoi.h" +#include +#include +#include +#include +#include +#include +#include +//////////////////////////////////////////////////////////////////////////////// + +using namespace voroffset; + +//////////////////////////////////////////////////////////////////////////////// + +namespace +{ + + // ----------------------------------------------------------------------------- + + template + T sign(T a) { return a > 0 ? T(1) : T(-1); } + + // Works only for non-negative numbers (≥ 0) + /*template + IntType ceilDiv(IntType x, IntType y) { + return (x / y + (x % y != 0)); + }*/ + + // Auxiliary methods + template + T det(T a, T b, T c, T d) + { + return a * d - b * c; + } + + // Returns true iff a ≤ sqrt(x) + template + bool infSqrt(T a, T x) + { + return (a <= 0 || a*a <= x); + } + + // Returns true iff sqrt(x) ≤ b + template + bool supSqrt(T b, T x) + { + return (b >= 0 && x <= b*b); + } + + static inline bool ceilIfGreateq(double x, double xmin, double &xdst) + { + if (x >= xmin) + { + xdst = x; + return true; + } + else + { + return false; + } + } + + // Typedefs + typedef VoronoiMorphoF::SegmentF Segment; + //typedef double real; + + double X(PointF p) { return p.real(); } + double Y(PointF p) { return p.imag(); } + + //////////////////////////////////////////////////////////////////////////////// + + //template + struct Ray + { + typedef double Scalar; + + Scalar a, b, c; + + Ray(Scalar _a, Scalar _b, Scalar _c) + : a(_a), b(_b), c(_c) + { } + + // Perpendicular bisector of p and q, equation is a*x + b*y = c + Ray(PointF p, PointF q) + : a(2 * (X(p) - X(q))) + , b(2 * (Y(p) - Y(q))) + , c(X(p) * X(p) - X(q) * X(q) + Y(p) * Y(p) - Y(q) * Y(q)) + {} + + inline double inter(Scalar x) const + { + return (c - a*x)/b; + } + + inline bool interBefore(int x, double yr) const + { + // Precondition: b==0 + return ((c - a*x) * sign(b) <= yr * b * sign(b)); + } + + inline bool interAfter(int x, double yl) const + { + // Precondition: b==0 + return (yl * b * sign(b) <= (c - a*x) * sign(b)); + } + }; + + //////////////////////////////////////////////////////////////////////////////// + + /* + * Assumes that y_p < y_q < y_r + * Returns true if the intersection is between y1 and y2 + */ + bool ray_intersect_range(PointF p, PointF q, PointF r, double y1, double y2, double &xinter) + { + if (X(p) == X(q) && X(q) == X(r)) + { + return false; + } + else + { + const Ray r1(p, q); + const Ray r2(q, r); + const double denom = det(r1.a, r2.a, r1.b, r2.b); + if (denom == 0) + { + return false; + } + else + { + const double ynumer = det(r1.a, r2.a, r1.c, r2.c); + if (y1 * denom * sign(denom) <= ynumer * sign(denom) + && ynumer * sign(denom) <= y2 * denom * sign(denom)) + { + const double xnumer = det(r1.c, r2.c, r1.b, r2.b); + if (xnumer * sign(denom) > X(q) * denom * sign(denom)) + { + xinter = ((double)xnumer / denom); + return true; + } + } + return false; + } + } + } + + bool ray_intersect_after(PointF p, PointF q, PointF r, double y, double &xinter) + { + if (X(p) == X(q) && X(q) == X(r)) + { + return false; + } + else + { + const Ray r1(p, q); + const Ray r2(q, r); + const double denom = det(r1.a, r2.a, r1.b, r2.b); + if (denom == 0) + { + return false; + } + else + { + const double ynumer = det(r1.a, r2.a, r1.c, r2.c); + if (y * denom * sign(denom) <= ynumer * sign(denom)) + { + const double xnumer = det(r1.c, r2.c, r1.b, r2.b); + if (xnumer * sign(denom) > X(q) * denom * sign(denom)) + { + xinter = ((double)xnumer / denom); + return true; + } + } + return false; + } + } + } + + bool ray_intersect_before(PointF p, PointF q, PointF r, double y, double &xinter) + { + if (X(p) == X(q) && X(q) == X(r)) + { + return false; + } + else + { + const Ray r1(p, q); + const Ray r2(q, r); + const double denom = det(r1.a, r2.a, r1.b, r2.b); + if (denom == 0) + { + return false; + } + else + { + const double ynumer = det(r1.a, r2.a, r1.c, r2.c); + if (ynumer * sign(denom) <= y * denom * sign(denom)) + { + const double xnumer = det(r1.c, r2.c, r1.b, r2.b); + if (xnumer * sign(denom) > X(q) * denom * sign(denom)) + { + xinter = ((double)xnumer / denom); + return true; + } + } + return false; + } + } + } + + /* + * Finds a point that is at equal distance of p, q and s + * Assumes that s.x < p.x && s.x < q.x + */ + bool voronoi_vertex(Segment s, PointF p, PointF q, double &xinter) + { + /* + * We seek the next voronoi point within the range [s.y1, s.y2]. + * Let (xm, ym) be its coordinates. Per the previous sentence we must have: + * (a) xm >= s.x + * (b) s.y1 ≤ ym ≤ s.y2 + * + * Let r be the bisector of p and q. We suppose that the point we seek is on r. + * As the squared distance from (xm,ym) to s is given by (xm - s.x)^2, we must have: + * + * (xm - s.x)^2 = (xm - xp)^2 + (ym - yp)^2 (1) + * = (xm - xq)^2 + (ym - yq)^2 (2) + * + * The following code solves (1) for ym, and returns true iff (a) and (b) + */ + const Ray r(p, q); + const double u = 2.0 * (X(p) - s.x); + const double w = s.x * s.x - X(p) * X(p); + double A = r.a; + double B = r.b * u - 2.0 * r.a * Y(p); + double C = r.a * Y(p) * Y(p) - w * r.a - r.c * u; + //reduce(A, B, C); + const double delta = B * B - 4.0 * A * C; + //const ll alpha = 2 * A * A * s.y1 + 2 * A * B; + //const ll beta = 2 * A * A * s.y2 + 2 * A * B; + const double alpha = 2.0 * std::abs(A) * s.y1 + sign(A) * B; + const double beta = 2.0 * std::abs(A) * s.y2 + sign(A) * B; + const bool pos = (A > 0); + // We must have alpha ≤ sign(a) * ±sqrt(delta) ≤ beta + if (A == 0) + { + if (2.0 * s.y1 <= Y(p) + Y(q) && Y(p) + Y(q) <= 2.0 * s.y2) + { + vor_assert(X(p) != s.x); + const double xm = ((Y(p) - Y(q)) * (Y(p) - Y(q)) - 4.0 * w) / (4.0 * u); + return ceilIfGreateq(xm, s.x, xinter); + } + else + { + return false; + } + } + else if (delta < 0) + { + return false; + } + else if (delta == 0 && alpha <= 0 && 0 <= beta) + { + // root y0 = -B / 2A is within range + double xm = (2.0 * A * r.c + r.b * B) / (2.0 * r.a * A); + return ceilIfGreateq(xm, s.x, xinter); + } + else if ( + pos ? infSqrt(alpha, delta) && supSqrt(beta, delta) + : supSqrt(-alpha, delta) && infSqrt(-beta, delta)) + { + // root y1 = (-B+sqrt(delta)) / 2A is within range + const double xm = (2.0 * A * r.c + r.b * B - r.b * std::sqrt(delta)) / (2.0 * r.a * A); + return ceilIfGreateq(xm, s.x, xinter); + } + else if ( + !pos ? infSqrt(alpha, delta) && supSqrt(beta, delta) + : supSqrt(-alpha, delta) && infSqrt(-beta, delta)) + { + // root y2 = (-B-sqrt(delta)) / 2A is within range + const double xm = (2.0 * A * r.c + r.b * B + r.b * std::sqrt(delta)) / (2.0 * r.a * A); + return ceilIfGreateq(xm, s.x, xinter); + } + else + { + return false; + } + } + + // ----------------------------------------------------------------------------- + +} // anonymous namespace + + //////////////////////////////////////////////////////////////////////////////// + // VoronoiMorpho + //////////////////////////////////////////////////////////////////////////////// + + /* + * Assumes that y_p < y_q < y_r + */ +double VoronoiMorphoF::rayIntersect(PointF p, PointF q, PointF r) const +{ + if (X(p) == X(q) && X(q) == X(r)) + { + return m_XMax*1.0f; + } + else + { + const Ray r1(p, q); + const Ray r2(q, r); + const double denom = det(r1.a, r2.a, r1.b, r2.b); + if (denom == 0) + { + // Intersection at x >= m_XMax are discarded; + return m_XMax*1.0f; + } + else + { + const double numer = det(r1.c, r2.c, r1.b, r2.b); + if (numer * sign(denom) <= X(q) * denom * sign(denom)) + { + return m_XMax*1.0f; + } + else + { + return ((double)numer / denom); + } + } + } +} + +/* +* Assumes that y_lp < y_ab < y_qr +*/ +double VoronoiMorphoF::treatSegments(SegmentF lp, SegmentF ab, SegmentF qr) const +{ + if (ab.x >= lp.x && ab.x >= qr.x) + { + // Nothing to do + return m_XMax; + } + else if (ab.isPoint()) + { + if (lp.isPoint() && qr.isPoint()) + { + + // Case: lp = point, ab = point, qr = point + return rayIntersect(lp.any(), ab.any(), qr.any()); + } + else if (lp.isPoint()) + { + // Case: lp = point, ab = point, qr = segment + double xinter = m_XMax*1.0f; + if (ray_intersect_before(lp.any(), ab.any(), qr.left(), qr.y1, xinter)) + { + return xinter; + } + else if (ray_intersect_after(lp.any(), ab.any(), qr.right(), qr.y2, xinter)) + { + return xinter; + } + else if (voronoi_vertex(qr, lp.any(), ab.any(), xinter)) + { + return xinter; + } + else + { + return xinter; + } + } + else if (qr.isPoint()) + { + // Case: lp = segment, ab = point, qr = point + double xinter = m_XMax*1.0f; + if (ray_intersect_before(lp.left(), ab.any(), qr.any(), lp.y1, xinter)) + { + return xinter; + } + else if (ray_intersect_after(lp.right(), ab.any(), qr.any(), lp.y2, xinter)) + { + return xinter; + } + else if (voronoi_vertex(lp, ab.any(), qr.any(), xinter)) + { + return xinter; + } + else + { + return xinter; + } + } + else + { + // Case: lp = segment, ab = point, qr = segment + double xinter = m_XMax*1.0f; + if (ray_intersect_before(lp.left(), ab.any(), qr.left(), lp.y1, xinter)) + { + return xinter; + } + else if (ray_intersect_after(lp.right(), ab.any(), qr.right(), qr.y2, xinter)) + { + return xinter; + } + else if (ray_intersect_range(lp.right(), ab.any(), qr.left(), lp.y2, qr.y1, xinter)) + { + return xinter; + } + else if (voronoi_vertex(lp, ab.any(), qr.left(), xinter)) + { + return xinter; + } + else if (voronoi_vertex(qr, lp.right(), ab.any(), xinter)) + { + return xinter; + } + else + { + return xinter; + } + } + } + else + { + if (lp.isPoint() && qr.isPoint()) + { + // Case: lp = point, ab = segment, qr = point + double xinter = m_XMax*1.0f; + if (ray_intersect_before(lp.any(), ab.left(), qr.any(), ab.y1, xinter)) + { + return xinter; + } + else if (ray_intersect_after(lp.any(), ab.right(), qr.any(), ab.y2, xinter)) + { + return xinter; + } + else if (voronoi_vertex(ab, lp.any(), qr.any(), xinter)) + { + return xinter; + } + else + { + return xinter; + } + } + else if (lp.isPoint()) + { + // Case: lp = point, ab = segment, qr = segment + double xinter = m_XMax*1.0f; + if (ray_intersect_before(lp.any(), ab.left(), qr.left(), ab.y1, xinter)) + { + return xinter; + } + else if (ray_intersect_after(lp.any(), ab.right(), qr.right(), qr.y2, xinter)) + { + return xinter; + } + else if (ray_intersect_range(lp.any(), ab.right(), qr.left(), ab.y2, qr.y1, xinter)) + { + return xinter; + } + else if (voronoi_vertex(ab, lp.any(), qr.left(), xinter)) + { + return xinter; + } + else if (voronoi_vertex(qr, lp.any(), ab.right(), xinter)) + { + return xinter; + } + else + { + return xinter; + } + } + else if (qr.isPoint()) + { + // Case: lp = segment, ab = segment, qr = point + double xinter = m_XMax*1.0f; + + if (ray_intersect_before(lp.left(), ab.left(), qr.any(), lp.y1, xinter)) + { + return xinter; + } + else if (ray_intersect_after(lp.right(), ab.right(), qr.any(), ab.y2, xinter)) + { + return xinter; + } + else if (ray_intersect_range(lp.right(), ab.left(), qr.any(), lp.y2, ab.y1, xinter)) + { + return xinter; + } + else if (voronoi_vertex(ab, lp.right(), qr.any(), xinter)) + { + return xinter; + } + else if (voronoi_vertex(lp, ab.left(), qr.any(), xinter)) { + + return xinter; + } + else + { + return xinter; + } + } + else + { + // Case: lp = segment, ab = segment, qr = segment + double xinter = m_XMax*1.0f; + if (ray_intersect_before(lp.left(), ab.left(), qr.left(), lp.y1, xinter)) + { + return xinter; + } + else if (ray_intersect_after(lp.right(), ab.right(), qr.right(), qr.y2, xinter)) + { + return xinter; + } + else if (ray_intersect_range(lp.right(), ab.left(), qr.left(), lp.y2, ab.y1, xinter)) + { + return xinter; + } + else if (ray_intersect_range(lp.right(), ab.right(), qr.left(), ab.y2, qr.y1, xinter)) + { + return xinter; + } + else if (voronoi_vertex(ab, lp.right(), qr.left(), xinter)) + { + return xinter; + } + else if (voronoi_vertex(lp, ab.left(), qr.left(), xinter)) + { + return xinter; + } + else if (voronoi_vertex(qr, lp.right(), ab.right(), xinter)) + { + return xinter; + } + else + { + return xinter; + } + } + } +} + +// ----------------------------------------------------------------------------- + +// Assumes that it != m_S.end() +void VoronoiMorphoF::exploreLeft(S_const_iter it, SegmentF qr, int i) +{ + while (it != m_S.begin()) + { + double xinter = treatSegments(*std::prev(it), *it, qr); + if (xinter <= i) + { + it = std::prev(m_S.erase(it)); + } + else + { + if (std::ceil(xinter) < m_XMax) + { + m_Q[std::ceil(xinter)].push_back(*it); + } + return; + } + } +} + +// Assumes that it != m_S.end() +void VoronoiMorphoF::exploreRight(S_const_iter it, SegmentF lp, int i) +{ + while (std::next(it) != m_S.end()) + { + double xinter = treatSegments(lp, *it, *std::next(it)); + if (xinter <= i) + { + it = m_S.erase(it); + } + else + { + if (std::ceil(xinter) < m_XMax) + { + m_Q[std::ceil(xinter)].push_back(*it); + } + return; + } + } +} + +// ----------------------------------------------------------------------------- + +// Insert a new seed segment (i,j1) -- (i,j2) +void VoronoiMorphoF::insertSegment(int i, double j1, double j2) +{ + const SegmentF s(i, j1, j2); + // 1st step: remove all segments occluded by s + S_const_iter it = m_S.lower_bound(s); + while (it != m_S.end() && s.contains(*it)) + { + it = m_S.erase(it); + } + while (it != m_S.begin() && s.contains(*std::prev(it))) + { + m_S.erase(std::prev(it)); + } + // 2nd step: See if we have to split the neighboring segments on the right + if (it != m_S.end()) + { + if (it->y1 <= s.y2) + { + int xsup = it->x + (int)std::floor(m_Radius) + 1; + // Special case: we may have to split the overlapping Segment + if (it->y1 < s.y1) + { + SegmentF lp(it->x, it->y1, s.y1); + m_S.insert(it, lp); + if (xsup < m_XMax) { m_Q[xsup].push_back(lp); } + } + SegmentF ab(it->x, s.y2, it->y2); + it = m_S.erase(it); + it = m_S.insert(it, ab); + if (xsup < m_XMax) { m_Q[xsup].push_back(ab); } + } + } + // And now look at the neighboring segments on the left + if (it != m_S.begin()) + { + S_const_iter jt = std::prev(it); + if (jt->y2 >= s.y1) + { + int xsup = jt->x + (int)std::floor(m_Radius) + 1; + // Special case: we may have to split the overlapping Segment + if (jt->y2 > s.y2) + { + SegmentF lp(jt->x, s.y2 , jt->y2); + it = m_S.insert(it, lp); + if (xsup < m_XMax) { m_Q[xsup].push_back(lp); } + } + SegmentF ab(jt->x, jt->y1, s.y1); + jt = m_S.erase(jt); + m_S.insert(jt, ab); + if (xsup < m_XMax) { m_Q[xsup].push_back(ab); } + } + } + // 3rd step: Insert the new segment in the list of seeds + S_const_iter jt = m_S.insert(it, s); + int xsup = s.x + (int)std::floor(m_Radius) + 1; + if (xsup < m_XMax) { m_Q[xsup].push_back(s); } + // 4th step: Inspect what's on the right and on the left + if (it != m_S.end()) + { + exploreRight(it, s, i); + } + if (jt != m_S.begin()) + { + exploreLeft(std::prev(jt), s, i); + } +} + +// ----------------------------------------------------------------------------- + +// Remove seeds that are not contributing anymore to the current sweep line (x==i) +void VoronoiMorphoF::removeInactiveSegments(int i) +{ + std::sort(m_Q[i].begin(), m_Q[i].end()); + while (!m_Q[i].empty()) + { + S_const_iter it = m_S.find(m_Q[i].back()); + m_Q[i].pop_back(); + if (it != m_S.end()) + { + it = m_S.erase(it); + if (it != m_S.begin() && it != m_S.end()) + { + exploreLeft(std::prev(it), *it, i); + exploreRight(it, *std::prev(it), i); + } + } + } + m_Q[i].clear(); +} + +//////////////////////////////////////////////////////////////////////////////// + +namespace { + + // Appends segment [a,b[ to the given line + void appendSegment(std::vector &nl, double a, double b) + { + // Ensures nl[-2] < a + while (!nl.empty() && nl.rbegin()[1] >= a) + { + b = std::max(b, nl.back()); + nl.pop_back(); + nl.pop_back(); + } + if (nl.empty() || nl.back() < a) + { + nl.push_back(a); + nl.push_back(b); + } + else if (nl.back() < b) + { + nl.back() = b; + } + vor_assert(nl.size() % 2 == 0); + } + +} // anonymous namespace + + // ----------------------------------------------------------------------------- + + // Extract the result for the current line +void VoronoiMorphoF::getLine(int i, std::vector &nl, int ysize) +{ + for (Segment s : m_S) + { + vor_assert((i - s.x) <= m_Radius); + const double dy = (std::sqrt(m_Radius * m_Radius - (i - s.x) * (i - s.x))); + const double a = s.y1 - dy > 0 ? s.y1-dy : 0; // Clamp at 0 no matter what + const double b = s.y2 + dy < ysize ? s.y2 + dy : ysize; // ysize is the actual line size + if (a > ysize || b < 0) + continue; // The segment is totally beyond the range + appendSegment(nl, a, b); + } +} + +//////////////////////////////////////////////////////////////////////////////// + +#ifndef WIN32 +// static_assert(std::is_standard_layout::value, "VoronoiMorpho is not standard layout!"); +#endif diff --git a/src/vor2d/DoubleVoronoi.h b/src/vor2d/DoubleVoronoi.h new file mode 100644 index 0000000..806140a --- /dev/null +++ b/src/vor2d/DoubleVoronoi.h @@ -0,0 +1,153 @@ +#pragma once + +//////////////////////////////////////////////////////////////////////////////// +#include "vor2d/Common.h" +#include "vor2d/Image.h" +#include +#include +#include +#include +//////////////////////////////////////////////////////////////////////////////// + +namespace voroffset +{ + + struct VoronoiMorphoF + { + //////////////// + // Data types // + // A horizontal segment + struct SegmentF + { + double x, y1, y2; + + SegmentF(double _x, double _y1, double _y2) + : x(_x), y1(_y1), y2(_y2) + { } + + inline PointF left() const { return PointF(x, y1); } + inline PointF right() const { return PointF(x, y2); } + inline PointF any() const { return PointF(x, y1); } + + bool operator <(const SegmentF &s) const + { + return (y1 + y2 < s.y1 + s.y2) || (y1 + y2 == s.y1 + s.y2 && x < s.x); + } + + inline bool contains(const SegmentF &o) const + { + return y1 <= o.y1 && o.y2 <= y2; + } + + inline bool isPoint() const + { + return y1 == y2; + } + }; + //////////////// + ////////////////// + // Data members // + ////////////////// + + // Constant parameters for the current dilation + const int m_XMax, m_YMax; + const double m_Radius; + + // Seeds above sorted by ascending y of their midpoint + std::set m_S; + + // Candidate Voronoi vertices sorted by ascending x + std::vector > m_Q; + + // Typedefs + typedef std::set::const_iterator S_const_iter; + + ///////////// + // Methods // + ///////////// + + // Assumes that y_p < y_q < y_r + double rayIntersect(PointF p, PointF q, PointF r) const; + + // Assumes that y_lp < y_ab < y_qr + double treatSegments(SegmentF lp, SegmentF ab, SegmentF qr) const; + + // Assumes that it != m_S.end() + void exploreLeft(S_const_iter it, SegmentF qr, int i); + + // Assumes that it != m_S.end() + void exploreRight(S_const_iter it, SegmentF lp, int i); + + // Insert a new seed segment [(i,j1), (i,j2)] + void insertSegment(int i, double j1, double j2); + + // Remove seeds that are not contributing anymore to the current sweep line (y==i) + void removeInactiveSegments(int i); + + // Extract the result for the current line + void getLine(int i, std::vector &nl, int ysize); + + VoronoiMorphoF(int _xmax, int _ymax, double _radius) + : m_XMax(_xmax) + , m_YMax(_ymax) + , m_Radius(_radius) + , m_Q(m_XMax) + { } + }; + + //////////////////////////////////////////////////////////////////////////////// + + // Forward sweep for dilation operation + template + void voronoiF_half_dilate( + int xsize, int ysize, double radius, + Iterator it, std::vector > &result) + { + VoronoiMorphoF voronoi(xsize, ysize, radius); + result.assign(xsize, std::vector()); + for (int i = 0; i < xsize; ++i) + { + voronoi.removeInactiveSegments(i); + vor_assert(it[i].size() % 2 == 0); + for (int k = (int)it[i].size() - 2; k >= 0; k -= 2) + { + voronoi.insertSegment(i, it[i][k], it[i][k + 1]); + } + voronoi.getLine(i, result[i], ysize); + } + } + + // ----------------------------------------------------------------------------- + + // Forward sweep for erosion operation + template + void voronoiF_half_erode( + int xsize, int ysize, double radius, + Iterator it, std::vector > &result) + { + VoronoiMorphoF voronoi(xsize, ysize + 1, radius); + result.assign(xsize, std::vector()); + if (xsize == 0) { return; } + // Extra segment before first row + voronoi.insertSegment(-1, -1, ysize); + for (int i = 0; i < xsize; ++i) + { + voronoi.removeInactiveSegments(i); + vor_assert(it[i].size() % 2 == 0); + // Flip segments + insert an extra one at the extremities + double j2 = ysize*1.0f; + for (int k = (int)it[i].size() - 2; k >= 0; k -= 2) + { + voronoi.insertSegment(i, it[i][k + 1], j2); + j2 = it[i][k]; + } + voronoi.insertSegment(i, -1, j2); + // Retrieve result + voronoi.getLine(i, result[i], ysize); + } + } + + //////////////////////////////////////////////////////////////////////////////// + +} // namespace voroffset + diff --git a/src/vor2d/Image.h b/src/vor2d/Image.h new file mode 100644 index 0000000..29faf54 --- /dev/null +++ b/src/vor2d/Image.h @@ -0,0 +1,45 @@ +#pragma once + +//////////////////////////////////////////////////////////////////////////////// +#include "vor2d/Common.h" +#include +#include +#include +//////////////////////////////////////////////////////////////////////////////// + +namespace voroffset +{ + +// Thin wrapper around Eigen::Matrix for representing images. +// In a matrix M, elements are indexed by [row,col], while we are used to index +// elements of an image by [x,y] coordinates. So, we will use a small proxy so +// that we don't have to swap the [x,y] coefficients when accessing an element. +// +// Also, we store images as row-major matrices, as we expect pixels to be layout +// by increasing y*w+x index. + +template +class Image +{ +private: + typedef Eigen::Matrix MatrixXs; + + MatrixXs _data; + +public: + Image(int w = 0, int h = 0) : _data(h, w) { } + + void resize(int w, int h) { _data.resize(h, w); } + + auto width() const { return _data.cols(); } + auto height() const { return _data.rows(); } + + auto operator() (int x, int y) const { return _data(y, x); } + auto & operator() (int x, int y) { return _data(y, x); } + + // low-level data access + const MatrixXs & data() const { return _data; } + MatrixXs & data() { return _data; } +}; + +} // namespace voroffset diff --git a/src/vor2d/MorphologyOperators.cpp b/src/vor2d/MorphologyOperators.cpp new file mode 100644 index 0000000..5d705ab --- /dev/null +++ b/src/vor2d/MorphologyOperators.cpp @@ -0,0 +1,307 @@ +#include"vor2d/MorphologyOperators.h" + +using namespace voroffset; + +#ifndef MIN_SEG_SIZE +#define MIN_SEG_SIZE 1e-10 +#endif + +// Appends segment [a,b] to the given sorted line +void voroffset::appendSegment(std::vector &nl, double a, double b) +{ + // Ensures nl[-2] < a + while (!nl.empty() && nl.rbegin()[1] >= a) + { + b = std::max(b, nl.back()); + nl.pop_back(); + nl.pop_back(); + } + if (nl.empty() || nl.back() < a) + { + nl.push_back(a); + nl.push_back(b); + } + else if (nl.back() < b) + { + nl.back() = b; + } + vor_assert(nl.size() % 2 == 0); +} + +void voroffset::appendSegment(std::vector &n1, m_Segment seg) +{ + while (!n1.empty() && n1.back()._begin_pt >= seg._begin_pt) + { + seg._end_pt = std::max(seg._end_pt, n1.back()._end_pt); + n1.pop_back(); + } + if (n1.empty() || n1.back()._end_pt < seg._begin_pt) + n1.push_back(seg); + else if (n1.back()._end_pt < seg._end_pt) + n1.back()._end_pt = seg._end_pt; +} +////////////////////////////////////////////////////////////////////////////// + +// add segment at the end of the sorted segment vector +template +void voroffset::addSegmentAtTheEnd(Container & dst, double start, double end) +{ + if (dst.empty() || (dst.back() < start)) + { + dst.push_back(start); + dst.push_back(end); + } + else if (end > dst.back()) + { + dst.back() = end; + } +} + +template +void voroffset::addSegmentAtTheEnd(Container & dst, m_Segment segment) +{ + if (dst.empty() || (dst.back()._end_pt < segment._begin_pt)) + dst.push_back(segment); + else if (segment._end_pt > dst.back()._end_pt) + { + dst.back()._end_pt = segment._end_pt; + } +} +//////////////////////////////////////////////////////////////////// + +// get the next segment of the given position +void voroffset::getNextSegment(const std::vector & segs, std::vector::const_iterator & it, + double & left, double & right, bool & full) +{ + if (it == segs.end()) + { + full = false; + } + else + { + left = *it++; + right = *it++; + } +} + +void voroffset::getNextSegment(const std::vector & segs, std::vector::const_iterator & it, + m_Segment &next_seg, bool & full) +{ + if (it == segs.end()) + { + full = false; + } + else + { + next_seg = *it++; + } +} +//////////////////////////////////////////////////////////////////////////////////// + +// Computes the union of two sorted lists of segments. +// Uses a parallel sweep of both lists, akin to merge-sort. +void voroffset::unionSegs(const std::vector & a, const std::vector & b, std::vector & result) +{ + if (a.empty()) + { + result.insert(result.end(), b.begin(), b.end()); + return; + } + if (b.empty()) + { + result.insert(result.end(), a.begin(), a.end()); + return; + } + std::vector::const_iterator ia(a.begin()), ib(b.begin()); + double al, ar, bl, br; + al = *ia++; ar = *ia++; + bl = *ib++; br = *ib++; + bool bfull(true), afull(true); + while (bfull || afull) + { + if (bfull && bl <= al) + { + addSegmentAtTheEnd(result, bl, br); + getNextSegment(b, ib, bl, br, bfull); + } + else if (afull) + { + addSegmentAtTheEnd(result, al, ar); + if (ia == a.end()) + { + afull = false; + al = std::numeric_limits::max(); + } + else + { + al = *ia++; + ar = *ia++; + } + } + } +} + +void voroffset::unionSegs(const std::vector & a, const std::vector & b, std::vector & result) +{ + if (a.empty()) + { + result.insert(result.end(), b.begin(), b.end()); + return; + } + if (b.empty()) + { + result.insert(result.end(), a.begin(), a.end()); + return; + } + std::vector::const_iterator ia(a.begin()), ib(b.begin()); + m_Segment a_1, b_1; + a_1 = *ia++; + b_1 = *ib++; + bool bfull(true), afull(true); + while (bfull || afull) + { + if (bfull && b_1._begin_pt <= a_1._begin_pt) + { + addSegmentAtTheEnd(result, b_1); + getNextSegment(b, ib, b_1, bfull); + } + else if (afull) + { + addSegmentAtTheEnd(result, a_1); + if (ia == a.end()) + { + afull = false; + a_1._begin_pt = std::numeric_limits::max(); + } + else + a_1 = *ia++; + } + } +} +////////////////////////////////////////////////////////////////////////////// + +// Unoin map +/* + +template +void voroffset::unionMap(Iterator1 a, Iterator2 b, std::vector &res) +{ + for (size_t i = 0; i < res.size(); ++i) + { + SegVector temp; + unionSegs(a[i], b[i], temp); + res[i].swap(temp); + } +}*/ +////////////////////////////////////////////////////////////////////////////////// + +// negate ray +void voroffset::negate_ray(std::vector &result, double y_min, double y_max) +{ + int size_ = result.size(); + + if (result.empty()) + { + result = { y_min, y_max }; + } + else + { + // First event + if (result[0] == y_min) + { + result.erase(result.begin()); + } + else + { + result.insert(result.begin(), y_min); + } + // Last event + if (result.back() == y_max) + { + result.pop_back(); + } + else + { + result.push_back(y_max); + } + } +} + +void voroffset::negate_ray(std::vector input, std::vector &result, double y_min, double y_max) +{ + int size_ = input.size(); + result.clear(); + if (input.empty()) + { + result = { m_Segment(y_min, y_max ,1) }; + } + else + { + // First event + if (input[0]._begin_pt > y_min) + result.push_back(m_Segment(y_min, input[0]._begin_pt, 1)); + for (int i = 0; i < size_ - 1; i++) + result.push_back(m_Segment(input[i]._end_pt, input[i + 1]._begin_pt, 1)); + // End event + if (input[size_ - 1]._end_pt < y_max) + result.push_back(m_Segment(input[size_ - 1]._end_pt, y_max, 1)); + } +} +//////////////////////////////////////////////////////////////////////////////////////// + +// calculate the xor between two vectors of segments +void voroffset::calculate_ray_xor(std::vector ray_1, std::vector ray_2, std::vector &ray_xor, double y_min, double y_max) +{ + ray_xor.clear(); + std::vector record_1, record_2, tmpt_union_1, tmpt_union_2, before_remove_pts; + negate_ray(ray_1, record_1, y_min, y_max); + negate_ray(ray_2, record_2, y_min, y_max); + unionSegs(record_1, ray_2, tmpt_union_1); + unionSegs(record_2, ray_1, tmpt_union_2); + negate_ray(tmpt_union_1, ray_1, y_min, y_max); + negate_ray(tmpt_union_2, ray_2, y_min, y_max); + unionSegs(ray_1, ray_2, before_remove_pts); + removepoint(before_remove_pts, ray_xor); +} + +void voroffset::calculate_ray_xor(std::vector ray_1, std::vector ray_2, std::vector &ray_xor, double y_min, double y_max) +{ + ray_xor.clear(); + std::vector record_1, record_2, tmpt_union_1, tmpt_union_2, before_remove_pts; + record_1.resize(ray_1.size()); + record_2.resize(ray_2.size()); + std::copy_n(ray_1.begin(), ray_1.size(), record_1.begin()); + std::copy_n(ray_2.begin(), ray_2.size(), record_2.begin()); + negate_ray(record_1, y_min, y_max); + negate_ray(record_2, y_min, y_max); + unionSegs(record_1, ray_2, tmpt_union_1); + unionSegs(record_2, ray_1, tmpt_union_2); + negate_ray(tmpt_union_1, y_min, y_max); + negate_ray(tmpt_union_2, y_min, y_max); + unionSegs(tmpt_union_1, tmpt_union_2, before_remove_pts); + removepoint(before_remove_pts, ray_xor); +} +///////////////////////////////////////////////////////////////////////////// + +// remove the segment when its size is too small to be seen +void voroffset::removepoint(const std::vector segs, std::vector &res) +{ + int size_ = segs.size(); + for (int i = 0; i < size_; i++) + { + if (segs[i]._end_pt - segs[i]._begin_pt < MIN_SEG_SIZE) + continue; + res.push_back(segs[i]); + } +} +void voroffset::removepoint(const std::vector segs, std::vector &res) +{ + int size_ = segs.size(); + for (int i = 0; i < size_; i += 2) + { + if (segs[i + 1] - segs[i] < MIN_SEG_SIZE) + continue; + res.push_back(segs[i]); + res.push_back(segs[i + 1]); + } +} \ No newline at end of file diff --git a/src/vor2d/MorphologyOperators.h b/src/vor2d/MorphologyOperators.h new file mode 100644 index 0000000..4dcf772 --- /dev/null +++ b/src/vor2d/MorphologyOperators.h @@ -0,0 +1,41 @@ +#pragma once +#include "vor2d/Common.h" +// Dealing with the operation on the same sweepline, such as union segments, negate ray, append segment to the line and so on + +namespace voroffset +{ + // Appends segment [a,b] to the given sorted line + void appendSegment(std::vector &nl, double a, double b); + void appendSegment(std::vector &n1, m_Segment seg); + + // add segment at the end of the sorted segment vector + template + void addSegmentAtTheEnd(Container & dst, double start, double end); + template + void addSegmentAtTheEnd(Container & dst, m_Segment segment); + + + // get the next segment of the given position + void getNextSegment(const std::vector & segs, std::vector::const_iterator & it, + double & left, double & right, bool & full); + void getNextSegment(const std::vector & segs, std::vector::const_iterator & it, + m_Segment & next_seg, bool & full); + + // Computes the union of two sorted lists of segments. + // Uses a parallel sweep of both lists, akin to merge-sort. + void unionSegs(const std::vector & a, const std::vector & b, std::vector & result); + void unionSegs(const std::vector & a, const std::vector & b, std::vector & result); + + // negate ray + void negate_ray(std::vector &result, double y_min, double y_max); + void negate_ray(std::vector input, std::vector &result, double y_min, double y_max); + + // calculate the xor of two rays + void calculate_ray_xor(std::vector ray_1, std::vector ray_2, std::vector &ray_xor, double y_min, double y_max); + void calculate_ray_xor(std::vector ray_1, std::vector ray_2, std::vector &ray_xor, double y_min, double y_max); + + // remove the segment when its size is too small to be seen, which can modify the viwer of visualization + void removepoint(const std::vector segs, std::vector &res); + void removepoint(const std::vector segs, std::vector &res); + +}//namespace voroffset \ No newline at end of file diff --git a/src/vor2d/Voronoi.cpp b/src/vor2d/Voronoi.cpp new file mode 100644 index 0000000..e730b28 --- /dev/null +++ b/src/vor2d/Voronoi.cpp @@ -0,0 +1,765 @@ +//////////////////////////////////////////////////////////////////////////////// +#include "vor2d/Voronoi.h" +#include +#include +#include +#include +#include +#include +#include +//////////////////////////////////////////////////////////////////////////////// + +using namespace voroffset; + +//////////////////////////////////////////////////////////////////////////////// + +namespace +{ + +// ----------------------------------------------------------------------------- + +template +IntType sign(IntType a) { return a > 0 ? IntType(1) : IntType(-1); } + +// Works only for non-negative numbers (≥ 0) +template +IntType ceilDiv(IntType x, IntType y) +{ + return (x / y + (x % y != 0)); +} + +// Auxiliary methods +template +IntType det(IntType a, IntType b, IntType c, IntType d) +{ + return a * d - b * c; +} + +template +IntType gcd(IntType a, IntType b) +{ + IntType aa = (a < 0 ? -a : a); + IntType bb = (b < 0 ? -b : b); + if (aa < bb) { std::swap(aa, bb); } + while (bb) + { + const IntType rest = (aa % bb); + aa = bb; + bb = rest; + } + return aa; +} + +template +double divideDouble(IntType a, IntType b) +{ + IntType d = gcd(a, b); + return (double) (a/d) / (double) (b/d); +} + +template +void reduce(IntType &a, IntType &b, IntType &c) +{ + IntType d = gcd(gcd(a, b), c); + if (d != 0) + { + a /= d; + b /= d; + c /= d; + } +} + +// Returns true iff a ≤ sqrt(x) +template +bool infSqrt(IntType a, IntType x) +{ + return (a <= 0 || a*a <= x); +} + +// Returns true iff sqrt(x) ≤ b +template +bool supSqrt(IntType b, IntType x) +{ + return (b >= 0 && x <= b*b); +} + +static inline bool ceilIfGreateq(double x, int xmin, int &xdst) +{ + if (x >= xmin) + { + xdst = (int) std::ceil(x); + return true; + } + else + { + return false; + } +} + +// Typedefs +typedef VoronoiMorpho::Segment Segment; +typedef VoronoiMorpho::Point Point; +typedef long long ll; + +ll X(Point p) { return p.real(); } +ll Y(Point p) { return p.imag(); } + +//////////////////////////////////////////////////////////////////////////////// + +//template +struct Ray +{ + typedef ll Scalar; + + Scalar a, b, c; + + Ray(Scalar _a, Scalar _b, Scalar _c) + : a(_a), b(_b), c(_c) + { } + + // Perpendicular bisector of p and q, equation is a*x + b*y = c + Ray(Point p, Point q) + : a(2 * (X(p) - X(q))) + , b(2 * (Y(p) - Y(q))) + , c(X(p) * X(p) - X(q) * X(q) + Y(p) * Y(p) - Y(q) * Y(q)) + { + reduce(a, b, c); + } + + inline double inter(Scalar x) const + { + return divideDouble(c - a*x, b); + } + + inline bool interBefore(int x, int yr) const + { + // Precondition: b==0 + return ((c - a*x) * sign(b) <= yr * b * sign(b)); + } + + inline bool interAfter(int x, int yl) const + { + // Precondition: b==0 + return (yl * b * sign(b) <= (c - a*x) * sign(b)); + } +}; + +//////////////////////////////////////////////////////////////////////////////// + +/* + * Assumes that y_p < y_q < y_r + * Returns true if the intersection is between y1 and y2 + */ +bool ray_intersect_range(Point p, Point q, Point r, int y1, int y2, int &xinter) +{ + if (X(p) == X(q) && X(q) == X(r)) + { + return false; + } + else + { + const Ray r1(p, q); + const Ray r2(q, r); + const ll denom = det(r1.a, r2.a, r1.b, r2.b); + if (denom == 0) + { + return false; + } + else + { + const ll ynumer = det(r1.a, r2.a, r1.c, r2.c); + if (y1 * denom * sign(denom) <= ynumer * sign(denom) + && ynumer * sign(denom) <= y2 * denom * sign(denom)) + { + const ll xnumer = det(r1.c, r2.c, r1.b, r2.b); + if (xnumer * sign(denom) > X(q) * denom * sign(denom)) + { + xinter = (int) std::ceil((double)xnumer / denom); + return true; + } + } + return false; + } + } +} + +bool ray_intersect_after(Point p, Point q, Point r, int y, int &xinter) +{ + if (X(p) == X(q) && X(q) == X(r)) + { + return false; + } + else + { + const Ray r1(p, q); + const Ray r2(q, r); + const ll denom = det(r1.a, r2.a, r1.b, r2.b); + if (denom == 0) + { + return false; + } else + { + const ll ynumer = det(r1.a, r2.a, r1.c, r2.c); + if (y * denom * sign(denom) <= ynumer * sign(denom)) + { + const ll xnumer = det(r1.c, r2.c, r1.b, r2.b); + if (xnumer * sign(denom) > X(q) * denom * sign(denom)) + { + xinter = (int) std::ceil((double) xnumer / denom); + return true; + } + } + return false; + } + } +} + +bool ray_intersect_before(Point p, Point q, Point r, int y, int &xinter) +{ + if (X(p) == X(q) && X(q) == X(r)) + { + return false; + } + else + { + const Ray r1(p, q); + const Ray r2(q, r); + const ll denom = det(r1.a, r2.a, r1.b, r2.b); + if (denom == 0) + { + return false; + } + else + { + const ll ynumer = det(r1.a, r2.a, r1.c, r2.c); + if (ynumer * sign(denom) <= y * denom * sign(denom)) + { + const ll xnumer = det(r1.c, r2.c, r1.b, r2.b); + if (xnumer * sign(denom) > X(q) * denom * sign(denom)) + { + xinter = (int) std::ceil((double) xnumer / denom); + return true; + } + } + return false; + } + } +} + +/* + * Finds a point that is at equal distance of p, q and s + * Assumes that s.x < p.x && s.x < q.x + */ +bool voronoi_vertex(Segment s, Point p, Point q, int &xinter) +{ + /* + * We seek the next voronoi point within the range [s.y1, s.y2]. + * Let (xm, ym) be its coordinates. Per the previous sentence we must have: + * (a) xm >= s.x + * (b) s.y1 ≤ ym ≤ s.y2 + * + * Let r be the bisector of p and q. We suppose that the point we seek is on r. + * As the squared distance from (xm,ym) to s is given by (xm - s.x)^2, we must have: + * + * (xm - s.x)^2 = (xm - xp)^2 + (ym - yp)^2 (1) + * = (xm - xq)^2 + (ym - yq)^2 (2) + * + * The following code solves (1) for ym, and returns true iff (a) and (b) + */ + const Ray r(p, q); + const ll u = 2ll * (X(p) - s.x); + const ll w = s.x * s.x - X(p) * X(p); + ll A = r.a; + ll B = r.b * u - 2ll * r.a * Y(p); + ll C = r.a * Y(p) * Y(p) - w * r.a - r.c * u; + reduce(A, B, C); + const ll delta = B * B - 4ll * A * C; + //const ll alpha = 2 * A * A * s.y1 + 2 * A * B; + //const ll beta = 2 * A * A * s.y2 + 2 * A * B; + const ll alpha = 2ll * std::abs(A) * s.y1 + sign(A) * B; + const ll beta = 2ll * std::abs(A) * s.y2 + sign(A) * B; + const bool pos = (A > 0); + // We must have alpha ≤ sign(a) * ±sqrt(delta) ≤ beta + if (A == 0) + { + if (2ll * s.y1 <= Y(p) + Y(q) && Y(p) + Y(q) <= 2ll * s.y2) { + vor_assert(X(p) != s.x); + const double xm = ((Y(p) - Y(q)) * (Y(p) - Y(q)) - 4.0 * w) / (4.0 * u); + return ceilIfGreateq(xm, s.x, xinter); + } + else + { + return false; + } + } + else if (delta < 0) + { + return false; + } + else if (delta == 0 && alpha <= 0 && 0 <= beta) + { + // root y0 = -B / 2A is within range + double xm = (2.0 * A * r.c + r.b * B) / (2.0 * r.a * A); + return ceilIfGreateq(xm, s.x, xinter); + } + else if ( + pos ? infSqrt( alpha, delta) && supSqrt( beta, delta) + : supSqrt(-alpha, delta) && infSqrt(-beta, delta)) + { + // root y1 = (-B+sqrt(delta)) / 2A is within range + const double xm = (2.0 * A * r.c + r.b * B - r.b * std::sqrt(delta)) / (2.0 * r.a * A); + return ceilIfGreateq(xm, s.x, xinter); + } + else if ( + !pos ? infSqrt( alpha, delta) && supSqrt( beta, delta) + : supSqrt(-alpha, delta) && infSqrt(-beta, delta)) + { + // root y2 = (-B-sqrt(delta)) / 2A is within range + const double xm = (2.0 * A * r.c + r.b * B + r.b * std::sqrt(delta)) / (2.0 * r.a * A); + return ceilIfGreateq(xm, s.x, xinter); + } + else + { + return false; + } +} + +// ----------------------------------------------------------------------------- + +} // anonymous namespace + +//////////////////////////////////////////////////////////////////////////////// +// VoronoiMorpho +//////////////////////////////////////////////////////////////////////////////// + +/* + * Assumes that y_p < y_q < y_r + */ +int VoronoiMorpho::rayIntersect(Point p, Point q, Point r) const +{ + if (X(p) == X(q) && X(q) == X(r)) + { + return m_XMax; + } + else + { + const Ray r1(p, q); + const Ray r2(q, r); + const ll denom = det(r1.a, r2.a, r1.b, r2.b); + if (denom == 0) + { + // Intersection at x >= m_XMax are discarded; + return m_XMax; + } + else + { + const ll numer = det(r1.c, r2.c, r1.b, r2.b); + if (numer * sign(denom) <= X(q) * denom * sign(denom)) + { + return m_XMax; + } + else + { + return (int) std::ceil((double)numer / denom); + } + } + } +} + +/* + * Assumes that y_lp < y_ab < y_qr + */ +int VoronoiMorpho::treatSegments(Segment lp, Segment ab, Segment qr) const +{ + if (ab.x >= lp.x && ab.x >= qr.x) + { + // Nothing to do + return m_XMax; + } + else if (ab.isPoint()) + { + if (lp.isPoint() && qr.isPoint()) + { + // Case: lp = point, ab = point, qr = point + return rayIntersect(lp.any(), ab.any(), qr.any()); + } + else if (lp.isPoint()) + { + // Case: lp = point, ab = point, qr = segment + int xinter = m_XMax; + if (ray_intersect_before(lp.any(), ab.any(), qr.left(), qr.y1, xinter)) + { + return xinter; + } + else if (ray_intersect_after(lp.any(), ab.any(), qr.right(), qr.y2, xinter)) + { + return xinter; + } + else if (voronoi_vertex(qr, lp.any(), ab.any(), xinter)) + { + return xinter; + } + else + { + return xinter; + } + } + else if (qr.isPoint()) + { + // Case: lp = segment, ab = point, qr = point + int xinter = m_XMax; + if (ray_intersect_before(lp.left(), ab.any(), qr.any(), lp.y1, xinter)) + { + return xinter; + } + else if (ray_intersect_after(lp.right(), ab.any(), qr.any(), lp.y2, xinter)) + { + return xinter; + } + else if (voronoi_vertex(lp, ab.any(), qr.any(), xinter)) + { + return xinter; + } + else + { + return xinter; + } + } + else + { + // Case: lp = segment, ab = point, qr = segment + int xinter = m_XMax; + if (ray_intersect_before(lp.left(), ab.any(), qr.left(), lp.y1, xinter)) + { + return xinter; + } + else if (ray_intersect_after(lp.right(), ab.any(), qr.right(), qr.y2, xinter)) + { + return xinter; + } + else if (ray_intersect_range(lp.right(), ab.any(), qr.left(), lp.y2, qr.y1, xinter)) + { + return xinter; + } + else if (voronoi_vertex(lp, ab.any(), qr.left(), xinter)) + { + return xinter; + } + else if (voronoi_vertex(qr, lp.right(), ab.any(), xinter)) + { + return xinter; + } + else + { + return xinter; + } + } + } + else + { + if (lp.isPoint() && qr.isPoint()) + { + // Case: lp = point, ab = segment, qr = point + int xinter = m_XMax; + if (ray_intersect_before(lp.any(), ab.left(), qr.any(), ab.y1, xinter)) + { + return xinter; + } + else if (ray_intersect_after(lp.any(), ab.right(), qr.any(), ab.y2, xinter)) + { + return xinter; + } + else if (voronoi_vertex(ab, lp.any(), qr.any(), xinter)) + { + return xinter; + } + else + { + return xinter; + } + } + else if (lp.isPoint()) + { + // Case: lp = point, ab = segment, qr = segment + int xinter = m_XMax; + if (ray_intersect_before(lp.any(), ab.left(), qr.left(), ab.y1, xinter)) + { + return xinter; + } + else if (ray_intersect_after(lp.any(), ab.right(), qr.right(), qr.y2, xinter)) + { + return xinter; + } + else if (ray_intersect_range(lp.any(), ab.right(), qr.left(), ab.y2, qr.y1, xinter)) + { + return xinter; + } + else if (voronoi_vertex(ab, lp.any(), qr.left(), xinter)) + { + return xinter; + } + else if (voronoi_vertex(qr, lp.any(), ab.right(), xinter)) + { + return xinter; + } + else + { + return xinter; + } + } + else if (qr.isPoint()) + { + // Case: lp = segment, ab = segment, qr = point + int xinter = m_XMax; + if (ray_intersect_before(lp.left(), ab.left(), qr.any(), lp.y1, xinter)) + { + return xinter; + } + else if (ray_intersect_after(lp.right(), ab.right(), qr.any(), ab.y2, xinter)) + { + return xinter; + } + else if (ray_intersect_range(lp.right(), ab.left(), qr.any(), lp.y2, ab.y1, xinter)) + { + return xinter; + } + else if (voronoi_vertex(ab, lp.right(), qr.any(), xinter)) + { + return xinter; + } + else if (voronoi_vertex(lp, ab.left(), qr.any(), xinter)) + { + return xinter; + } + else + { + return xinter; + } + } + + else + { + // Case: lp = segment, ab = segment, qr = segment + int xinter = m_XMax; + if (ray_intersect_before(lp.left(), ab.left(), qr.left(), lp.y1, xinter)) + { + return xinter; + } + else if (ray_intersect_after(lp.right(), ab.right(), qr.right(), qr.y2, xinter)) + { + return xinter; + } + else if (ray_intersect_range(lp.right(), ab.left(), qr.left(), lp.y2, ab.y1, xinter)) + { + return xinter; + } + else if (ray_intersect_range(lp.right(), ab.right(), qr.left(), ab.y2, qr.y1, xinter)) + { + return xinter; + } + else if (voronoi_vertex(ab, lp.right(), qr.left(), xinter)) + { + return xinter; + } + else if (voronoi_vertex(lp, ab.left(), qr.left(), xinter)) + { + return xinter; + } + else if (voronoi_vertex(qr, lp.right(), ab.right(), xinter)) + { + return xinter; + } + else + { + return xinter; + } + } + } +} + +// ----------------------------------------------------------------------------- + +// Assumes that it != m_S.end() +void VoronoiMorpho::exploreLeft(S_const_iter it, Segment qr, int i) +{ + while (it != m_S.begin()) + { + int xinter = treatSegments(*std::prev(it), *it, qr); + if (xinter <= i) + { + it = std::prev(m_S.erase(it)); + } + else + { + if (xinter < m_XMax) + { + m_Q[xinter].push_back(*it); + } + return; + } + } +} + +// Assumes that it != m_S.end() +void VoronoiMorpho::exploreRight(S_const_iter it, Segment lp, int i) +{ + while (std::next(it) != m_S.end()) + { + int xinter = treatSegments(lp, *it, *std::next(it)); + if (xinter <= i) + { + it = m_S.erase(it); + } + else + { + if (xinter < m_XMax) + { + m_Q[xinter].push_back(*it); + } + return; + } + } +} + +// ----------------------------------------------------------------------------- + +// Insert a new seed segment (i,j1) -- (i,j2) +void VoronoiMorpho::insertSegment(int i, int j1, int j2) +{ + const Segment s(i, j1, j2); + // 1st step: remove all segments occluded by s + S_const_iter it = m_S.lower_bound(s); + while (it != m_S.end() && s.contains(*it)) + { + it = m_S.erase(it); + } + while (it != m_S.begin() && s.contains(*std::prev(it))) + { + m_S.erase(std::prev(it)); + } + // 2nd step: See if we have to split the neighboring segments on the right + if (it != m_S.end()) + { + if (it->y1 <= s.y2) + { + int xsup = it->x + (int) std::floor(m_Radius) + 1; + // Special case: we may have to split the overlapping Segment + if (it->y1 < s.y1) + { + Segment lp(it->x, it->y1, s.y1 - 1); + m_S.insert(it, lp); + if (xsup < m_XMax) { m_Q[xsup].push_back(lp); } + } + Segment ab(it->x, s.y2 + 1, it->y2); + it = m_S.erase(it); + it = m_S.insert(it, ab); + if (xsup < m_XMax) { m_Q[xsup].push_back(ab); } + } + } + // And now look at the neighboring segments on the left + if (it != m_S.begin()) + { + S_const_iter jt = std::prev(it); + if (jt->y2 >= s.y1) + { + int xsup = jt->x + (int) std::floor(m_Radius) + 1; + // Special case: we may have to split the overlapping Segment + if (jt->y2 > s.y2) + { + Segment lp(jt->x, s.y2 + 1, jt->y2); + it = m_S.insert(it, lp); + if (xsup < m_XMax) { m_Q[xsup].push_back(lp); } + } + Segment ab(jt->x, jt->y1, s.y1 - 1); + jt = m_S.erase(jt); + m_S.insert(jt, ab); + if (xsup < m_XMax) { m_Q[xsup].push_back(ab); } + } + } + // 3rd step: Insert the new segment in the list of seeds + S_const_iter jt = m_S.insert(it, s); + int xsup = s.x + (int) std::floor(m_Radius) + 1; + if (xsup < m_XMax) { m_Q[xsup].push_back(s); } + // 4th step: Inspect what's on the right and on the left + if (it != m_S.end()) + { + exploreRight(it, s, i); + } + if (jt != m_S.begin()) + { + exploreLeft(std::prev(jt), s, i); + } +} + +// ----------------------------------------------------------------------------- + +// Remove seeds that are not contributing anymore to the current sweep line (y==i) +void VoronoiMorpho::removeInactiveSegments(int i) +{ + std::sort(m_Q[i].begin(), m_Q[i].end()); + while (!m_Q[i].empty()) + { + S_const_iter it = m_S.find(m_Q[i].back()); + m_Q[i].pop_back(); + if (it != m_S.end()) + { + it = m_S.erase(it); + if (it != m_S.begin() && it != m_S.end()) + { + exploreLeft(std::prev(it), *it, i); + exploreRight(it, *std::prev(it), i); + } + } + } + m_Q[i].clear(); +} + +//////////////////////////////////////////////////////////////////////////////// + +namespace +{ + + // Appends segment [a,b[ to the given line + void appendSegment(std::vector &nl, int a, int b) + { + // Ensures nl[-2] < a + while (!nl.empty() && nl.rbegin()[1] >= a) + { + b = std::max(b, nl.back()); + nl.pop_back(); + nl.pop_back(); + } + if (nl.empty() || nl.back() < a) + { + nl.push_back(a); + nl.push_back(b); + } + else if (nl.back() < b) + { + nl.back() = b; + } + vor_assert(nl.size() % 2 == 0); + } + +} // anonymous namespace + +// ----------------------------------------------------------------------------- + +// Extract the result for the current line +void VoronoiMorpho::getLine(int i, std::vector &nl, int ysize) +{ + for (Segment s : m_S) + { + vor_assert((i - s.x) <= m_Radius); + const int dy = std::floor(std::sqrt(m_Radius * m_Radius - (i - s.x) * (i - s.x))); + const int a = std::max(s.y1 - dy, 0); // Clamp at 0 no matter what + const int b = std::min(s.y2 + dy + 1, ysize); // ysize is the actual line size + appendSegment(nl, a, b); + } +} + +//////////////////////////////////////////////////////////////////////////////// + +#ifndef WIN32 +// static_assert(std::is_standard_layout::value, "VoronoiMorpho is not standard layout!"); +#endif diff --git a/src/vor2d/Voronoi.h b/src/vor2d/Voronoi.h new file mode 100644 index 0000000..1269c4a --- /dev/null +++ b/src/vor2d/Voronoi.h @@ -0,0 +1,156 @@ +#pragma once + +//////////////////////////////////////////////////////////////////////////////// +#include "vor2d/Common.h" +#include "vor2d/Image.h" +#include +#include +#include +//////////////////////////////////////////////////////////////////////////////// + +namespace voroffset +{ + +struct VoronoiMorpho +{ + //////////////// + // Data types // + //////////////// + + // A 2D point + typedef std::complex Point; + + // A horizontal segment + struct Segment + { + int x, y1, y2; + + Segment(int _x, int _y1, int _y2) + : x(_x), y1(_y1), y2(_y2) + { } + + inline Point left() const { return Point(x, y1); } + inline Point right() const { return Point(x, y2); } + inline Point any() const { return Point(x, y1); } + + bool operator <(const Segment &s) const + { + return (y1 + y2 < s.y1 + s.y2) || (y1 + y2 == s.y1 + s.y2 && x < s.x); + } + + inline bool contains(const Segment &o) const + { + return y1 <= o.y1 && o.y2 <= y2; + } + + inline bool isPoint() const + { + return y1 == y2; + } + }; + + ////////////////// + // Data members // + ////////////////// + + // Constant parameters for the current dilation + const int m_XMax, m_YMax; + const double m_Radius; + + // Seeds above sorted by ascending y of their midpoint + std::set m_S; + + // Candidate Voronoi vertices sorted by ascending x + std::vector > m_Q; + + // Typedefs + typedef std::set::const_iterator S_const_iter; + + ///////////// + // Methods // + ///////////// + + // Assumes that y_p < y_q < y_r + int rayIntersect(Point p, Point q, Point r) const; + + // Assumes that y_lp < y_ab < y_qr + int treatSegments(Segment lp, Segment ab, Segment qr) const; + + // Assumes that it != m_S.end() + void exploreLeft(S_const_iter it, Segment qr, int i); + + // Assumes that it != m_S.end() + void exploreRight(S_const_iter it, Segment lp, int i); + + // Insert a new seed segment [(i,j1), (i,j2)] + void insertSegment(int i, int j1, int j2); + + // Remove seeds that are not contributing anymore to the current sweep line (y==i) + void removeInactiveSegments(int i); + + // Extract the result for the current line + void getLine(int i, std::vector &nl, int ysize); + + VoronoiMorpho(int _xmax, int _ymax, double _radius) + : m_XMax(_xmax) + , m_YMax(_ymax) + , m_Radius(_radius) + , m_Q(m_XMax) + { } +}; + +//////////////////////////////////////////////////////////////////////////////// + +// Forward sweep for dilation operation +template +void voronoi_half_dilate( + int xsize, int ysize, double radius, + Iterator it, std::vector > &result) +{ + VoronoiMorpho voronoi(xsize, ysize, radius); + result.assign(xsize, std::vector()); + for (int i = 0; i < xsize; ++i) + { + voronoi.removeInactiveSegments(i); + vor_assert(it[i].size() % 2 == 0); + for (int k = (int) it[i].size() - 2; k >= 0; k -= 2) + { + voronoi.insertSegment(i, it[i][k], it[i][k + 1] - 1); + } + voronoi.getLine(i, result[i], ysize); + } +} + +// ----------------------------------------------------------------------------- + +// Forward sweep for erosion operation +template +void voronoi_half_erode( + int xsize, int ysize, double radius, + Iterator it, std::vector > &result) +{ + VoronoiMorpho voronoi(xsize, ysize + 1, radius); + result.assign(xsize, std::vector()); + if (xsize == 0) { return; } + // Extra segment before first row + voronoi.insertSegment(-1, -1, ysize); + for (int i = 0; i < xsize; ++i) + { + voronoi.removeInactiveSegments(i); + vor_assert(it[i].size() % 2 == 0); + // Flip segments + insert an extra one at the extremities + int j2 = ysize; + for (int k = (int) it[i].size() - 2; k >= 0; k -= 2) + { + voronoi.insertSegment(i, it[i][k + 1], j2); + j2 = it[i][k] - 1; + } + voronoi.insertSegment(i, -1, j2); + // Retrieve result + voronoi.getLine(i, result[i], ysize); + } +} + +//////////////////////////////////////////////////////////////////////////////// + +} // namespace voroffset diff --git a/src/vor3d/CMakeLists.txt b/src/vor3d/CMakeLists.txt new file mode 100644 index 0000000..6f685d4 --- /dev/null +++ b/src/vor3d/CMakeLists.txt @@ -0,0 +1,66 @@ +################################################################################ + +cmake_minimum_required(VERSION 3.3) +project(vor3d) + +################################################################################ + +add_library(${PROJECT_NAME} + Common.cpp + Common.h + CompressedVolume.cpp + CompressedVolume.h + CompressedVolumeBase.cpp + CompressedVolumeBase.h + CompressedVolumeWithRadii.cpp + CompressedVolumeWithRadii.h + Dexelize.cpp + Dexelize.h + HalfDilationOperator.cpp + HalfDilationOperator.h + HalfDilationOperator.hpp + Morpho2D.cpp + Morpho2D.h + MorphologyOperators.cpp + MorphologyOperators.h + MorphologyOperators.hpp + SeparatePower2D.cpp + SeparatePower2D.h + Timer.cpp + Timer.h + Voronoi.cpp + Voronoi.h + Voronoi2D.cpp + Voronoi2D.h + VoronoiBruteForce.cpp + VoronoiBruteForce.h + VoronoiVorPower.cpp + VoronoiVorPower.h +) + +target_include_directories(${PROJECT_NAME} PUBLIC ..) + +################################################################################ + +# Let's get a little bit paranoid +include(SetWarnings) +target_compile_options(${PROJECT_NAME} PRIVATE ${ALL_WARNINGS}) + +# Use C++14 +target_compile_features(${PROJECT_NAME} PUBLIC cxx_std_14) + +# Generate position independent code +set_target_properties(${PROJECT_NAME} PROPERTIES POSITION_INDEPENDENT_CODE ON) + +# Sanitizers +if(VOROFFSET_WITH_SANITIZERS) + add_sanitizers(${PROJECT_NAME}) +endif() + +# Dependencies +target_link_libraries(${PROJECT_NAME} PUBLIC Eigen3::Eigen geogram::geogram) + +# TBB library +if(VOROFFSET_WITH_TBB) + target_link_libraries(${PROJECT_NAME} PUBLIC tbb::tbb) +endif() diff --git a/src/vor3d/Common.cpp b/src/vor3d/Common.cpp new file mode 100644 index 0000000..815c5a1 --- /dev/null +++ b/src/vor3d/Common.cpp @@ -0,0 +1,18 @@ +//////////////////////////////////////////////////////////////////////////////// +#include "vor3d/Common.h" +#include +#include +#include +//////////////////////////////////////////////////////////////////////////////// + +[[noreturn]] void voroffset3d::assertion_failed( + const std::string& condition_string, + const std::string& file, int line) +{ + std::ostringstream os; + os << "Assertion failed: " << condition_string << ".\n"; + os << "File: " << file << ",\n"; + os << "Line: " << line; + + throw std::runtime_error(os.str()); +} \ No newline at end of file diff --git a/src/vor3d/Common.h b/src/vor3d/Common.h new file mode 100644 index 0000000..343d239 --- /dev/null +++ b/src/vor3d/Common.h @@ -0,0 +1,253 @@ +#pragma once + +//////////////////////////////////////////////////////////////////////////////// +#include +#include +#include +//////////////////////////////////////////////////////////////////////////////// +#ifndef EPS +#define EPS 1e-8 +#endif + +namespace voroffset3d +{ + // typedef// + typedef double Scalar; + typedef std::complex PointF; + typedef std::vector::iterator data_iterator; + /////////////////////////////////////////////////////////////////////////////////////////////////// + //Data Struct + ////////////////////////////////////////////////////////////////////////////////////////////////// + //Point struct + ////////////////////////////////////////////////////////////////////////////////////////////////// + struct Point // for voronoi diagram + { + double x, y; + Point(double _x, double _y) + :x(_x), y(_y) + {} + Point() = default; + bool operator <(const Point &p) const + { + return y < p.y || (y == p.y && x < p.x); + } + inline bool contains(const Point &o) const + { + return std::sqrt((x - o.x)*(x - o.x) + (y - o.y)*(y - o.y)) < EPS; + } + }; + + struct PointWithRadius // for initial data CompressedVolumeWithRadii + { + Scalar _pt, _radius; + PointWithRadius() = default; + PointWithRadius(Scalar pt, Scalar radius) + :_pt(pt), _radius(radius) + {} + bool operator<(const PointWithRadius &p) const + { + return _pt < p._pt || (_pt == p._pt && _radius < p._radius); + } + inline bool contains(const PointWithRadius &o) const + { + return _radius - o._radius > std::abs(_pt - o._pt) - EPS; + } + }; + + struct PointWithRadiusX // for power diagram + { + double x, y, r; + PointWithRadiusX(double _x, double _y, double _r) + :x(_x), y(_y), r(_r) + {} + PointWithRadiusX() = default; + bool operator <(const PointWithRadiusX &p) const + { + return y < p.y || (y == p.y && x < p.x) || (y == p.y && x == p.x && r < p.r); + } + + inline bool contains(const PointWithRadiusX &o) const + { + //return (y1 < o.y1 || std::abs(y1 - o.y1) < EPS) && (o.y2 < y2 || std::abs(o.y2 - y2)); + return r - o.r > std::sqrt(std::pow(x - o.x, 2.0) + std::pow(y - o.y, 2.0)) - EPS; + } + }; + + ///////////////////////////////////////////////////////////////////////////////////////////////////// + //Segment struct + //////////////////////////////////////////////////////////////////////////////////////////////////// + struct Segment // for voronoi diagram + { + double x, y1, y2; + + Segment(double _x, double _y1, double _y2) + : x(_x), y1(_y1), y2(_y2) + { } + Segment() = default; + inline Point left() const { return Point(x, y1); } + inline Point right() const { return Point(x, y2); } + inline Point any() const { return Point(x, y1); } + + bool operator <(const Segment &s) const + { + return (y1 + y2 < s.y1 + s.y2) || (y1 + y2 == s.y1 + s.y2 && x < s.x); + } + + inline bool contains(const Segment &o) const + { + return (y1 <= o.y1) && (o.y2 <= y2); + } + + inline bool isPoint() const + { + return std::abs(y1 - y2) < EPS; + } + inline bool isSplit(const Segment &o) const + /* return true if we need to split the current segment when we append o, which overlaps with current segment + Notice: the current segment and o have this below postion relationship: + z1 <= o.z1 + */ + { + return x < o.x && !isPoint(); + } + inline bool isOcclude(const Segment &o) const + { + return (y1 < o.y1 + EPS) && (o.y2 < y2 + EPS); + } + }; + + struct SegmentWithRadius // for initial data CompressedVolumeWithRadii + { + Scalar y1, y2, r; + SegmentWithRadius() = default; + SegmentWithRadius(Scalar begin_pt, Scalar end_pt, Scalar radius) + :y1(begin_pt), y2(end_pt), r(radius) + {} + inline PointWithRadius left() const { return PointWithRadius(y1, r); } + inline PointWithRadius right() const { return PointWithRadius(y2, r); } + inline PointWithRadius any() const { return PointWithRadius(y1, r); } + bool operator <(const SegmentWithRadius &s) const + { + return (y1 + y2 < s.y1 + s.y2) || (y1 + y2 == s.y1 + s.y2 && r < s.r); + } + inline bool contains(const SegmentWithRadius &o) const + { + //return (y1 < o.y1 || std::abs(y1 - o.y1) < EPS) && (o.y2 < y2 || std::abs(o.y2 - y2)); + return (y1 < o.y1 + EPS) && (o.y2 < y2 + EPS); + } + + inline bool isPoint() const + { + return std::abs(y1 - y2) < EPS; + } + inline bool isSplit(const SegmentWithRadius &o) const + /* return true if we need to split the current segment when we append o, which overlaps with current segment + Notice: the current segment and o have this below postion relationship: + z1 <= o.z1 + */ + { + return r < o.r && !isPoint(); + } + inline bool isOcclude(const SegmentWithRadius &o) const + { + return (y1 < o.y1 + EPS) && (o.y2 < y2 + EPS) && (r > o.r - EPS); + } + }; + + + struct SegmentWithRadiusX // for power diagram + { + double x, y1, y2, r; + + SegmentWithRadiusX(double _x, double _y1, double _y2, double _r) + : x(_x), y1(_y1), y2(_y2), r(_r) + { } + + SegmentWithRadiusX() = default; + + inline PointWithRadiusX left() const { return PointWithRadiusX(x, y1, r); } + inline PointWithRadiusX right() const { return PointWithRadiusX(x, y2, r); } + inline PointWithRadiusX any() const { return PointWithRadiusX(x, y1, r); } + + bool operator <(const SegmentWithRadiusX &s) const + { + return (y1 + y2 < s.y1 + s.y2) || (y1 + y2 == s.y1 + s.y2 && x < s.x) || (y1 + y2 == s.y1 + s.y2 && x == s.x && r < s.r); + } + + inline bool contains(const SegmentWithRadiusX &o) const + { + //return (y1 < o.y1 || std::abs(y1 - o.y1) < EPS) && (o.y2 < y2 || std::abs(o.y2 - y2)); + return (y1 < o.y1 + EPS) && (o.y2 < y2 + EPS); + } + + inline bool isPoint() const + { + return std::abs(y1 - y2) < EPS; + } + + inline bool isSplit(const SegmentWithRadiusX &o) const + /* return true if we need to split the current segment when we append o, which overlaps with current segment + Notice: the current segment and o have this below postion relationship: + z1 <= o.z1 + */ + { + return (r + x < o.r + o.x) && !isPoint(); + } + inline bool isOcclude(const SegmentWithRadiusX &o) const + { + return (y1 < o.y1 + EPS) && (o.y2 < y2 + EPS) && (r + x > o.r + o.x - EPS); + } + }; + + +//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +// Custom assert function +[[noreturn]] void assertion_failed( + const std::string& condition_string, + const std::string& file, int line); + +} // namespace voroffset + +// Shortcut alias +namespace vor3d = voroffset3d; + +//////////////////////////////////////////////////////////////////////////////// + +#define vor_assert_msg(x, s) \ +do \ +{ \ + if (!(x)) \ + { \ + voroffset3d::assertion_failed((s), __FILE__, __LINE__); \ + } \ +} while(0) + +// ----------------------------------------------------------------------------- + +#define vor_assert(x) \ +do \ +{ \ + if(!(x)) \ + { \ + voroffset3d::assertion_failed(#x, __FILE__, __LINE__); \ + } \ +} while (0) + +#ifdef voroffset_NO_DEBUG + +#define vor_debug(x) + +#else + +#define vor_debug(x) \ +do \ +{ \ + if(!(x)) \ + { \ + voroffset3d::assertion_failed(#x, __FILE__, __LINE__); \ + } \ +} while (0) + +#endif + diff --git a/src/vor3d/CompressedVolume.cpp b/src/vor3d/CompressedVolume.cpp new file mode 100644 index 0000000..947f0bd --- /dev/null +++ b/src/vor3d/CompressedVolume.cpp @@ -0,0 +1,152 @@ +//////////////////////////////////////////////////////////////////////////////// +#include "vor3d/CompressedVolume.h" +#include "vor3d/MorphologyOperators.h" +#include +//////////////////////////////////////////////////////////////////////////////// + +using namespace voroffset3d; + +//////////////////////////////////////////////////////////////////////////////// + +CompressedVolume::CompressedVolume( + Eigen::Vector3d origin, Eigen::Vector3d extent, Scalar voxel_size, int padding) +{ + m_Origin = origin; + m_Extent = extent; + m_Spacing = voxel_size; + m_Padding = padding; + m_Origin.array() -= padding * voxel_size; + m_GridSize[0] = (int)(std::ceil(extent[0] / m_Spacing) + 2 * padding); + m_GridSize[1] = (int)(std::ceil(extent[1] / m_Spacing) + 2 * padding); + m_Data.assign(m_GridSize[0] * m_GridSize[1], {}); + +} + +void CompressedVolume::iterate(int i, int j, std::function func) +{ + const auto &ray = at(i, j); + for (size_t k = 0; k+1 < ray.size(); k+=2) + { + func(ray[k], ray[k + 1]); + } +} + +void CompressedVolume::iterate(int i, int j, std::function func) +{ + const auto &ray = at(i, j); + for (size_t k = 0; k + 1 < ray.size(); k += 2) + { + func(ray[k], ray[k + 1], 0.f); + } +} + +void CompressedVolume::copy_volume_from(const CompressedVolume voxel) +{ + m_Origin = voxel.origin(); + m_Extent = voxel.extent(); + m_GridSize = voxel.gridSize(); + m_Padding = voxel.padding(); + m_Spacing = voxel.spacing(); + m_Data.assign(m_GridSize[0] * m_GridSize[1], {}); + //m_Data_Points.assign(m_GridSize[0] * m_GridSize[1], {}); + for(int x=0;x> m_Origin(0) >> m_Origin(1) >> m_Origin(2); + in >> m_Extent(0) >> m_Extent(1) >> m_Extent(2); + in >> m_GridSize(0) >> m_GridSize(1); + in >> m_Padding; + in >> m_Spacing; + m_Data.assign(m_GridSize(0)*m_GridSize(1), {}); + for (auto & row : m_Data) + { + size_t size; + in >> size; + row.resize(size); + for (auto & val : row) + { + in >> val; + } + } +} + +void CompressedVolume::save(std::ostream &out) const +{ + out << m_Origin(0) << ' ' << m_Origin(1) << ' ' << m_Origin(2) << "\n"; + out << m_Extent(0) << ' ' << m_Extent(1) << ' ' << m_Extent(2) << "\n"; + out << m_GridSize(0) << ' ' << m_GridSize(1) << "\n"; + out << m_Padding << "\n"; + out << m_Spacing << "\n"; + for (const auto & row : m_Data) + { + out << row.size(); + for (const auto & val : row) + { + out << ' ' << val; + } + out << "\n"; + } +} diff --git a/src/vor3d/CompressedVolume.h b/src/vor3d/CompressedVolume.h new file mode 100644 index 0000000..3469e5e --- /dev/null +++ b/src/vor3d/CompressedVolume.h @@ -0,0 +1,43 @@ +#pragma once + +//////////////////////////////////////////////////////////////////////////////// +#include "vor3d/CompressedVolumeBase.h" +//////////////////////////////////////////////////////////////////////////////// + +namespace voroffset3d +{ + + class CompressedVolume : public CompressedVolumeBase + { + private: + // Member data + std::vector > m_Data; + + public: + // Interface + CompressedVolume(Eigen::Vector3d origin, Eigen::Vector3d extent, Scalar voxel_size, int padding); + CompressedVolume() = default; + + // Marshalling + void save(std::ostream &out) const; + void load(std::istream &in); + + void copy_volume_from(const CompressedVolume voxel); + double get_volume(); + int numSegments(); + const std::vector & at(int x, int y) const { return m_Data[x + m_GridSize[0] * y]; } + std::vector & at(int x, int y) { return m_Data[x + m_GridSize[0] * y]; } + + + // Apply a function to each segment in the structure + virtual void iterate(int i, int j, std::function func) override; + virtual void iterate(int i, int j, std::function func1) override; + + virtual void appendSegment(int i, int j, Scalar begin_pt, Scalar end_pt, Scalar radius) override; + + virtual void reshape(int xsize, int ysize) override; + virtual void resize(int xsize, int ysize) override; + virtual void clear() override; + }; + +} // namespace voroffset3d diff --git a/src/vor3d/CompressedVolumeBase.cpp b/src/vor3d/CompressedVolumeBase.cpp new file mode 100644 index 0000000..b84c1a8 --- /dev/null +++ b/src/vor3d/CompressedVolumeBase.cpp @@ -0,0 +1,21 @@ +#include "vor3d/CompressedVolumeBase.h" + +using namespace voroffset3d; + +Eigen::Vector2d CompressedVolumeBase::dexelCenter(int x, int y) const +{ + Eigen::Vector2d pos; + pos[0] = (x + 0.5) * m_Spacing; + pos[1] = (y + 0.5) * m_Spacing; + return pos + m_Origin.head<2>(); +} + +void CompressedVolumeBase::reset(Eigen::Vector3d origin, Eigen::Vector3d extent, Scalar voxel_size, int padding, + int xsize, int ysize) +{ + m_Origin = origin; + m_Extent = extent; + m_Spacing = voxel_size; + m_Padding = padding; + reshape(xsize, ysize); +} diff --git a/src/vor3d/CompressedVolumeBase.h b/src/vor3d/CompressedVolumeBase.h new file mode 100644 index 0000000..5c8e629 --- /dev/null +++ b/src/vor3d/CompressedVolumeBase.h @@ -0,0 +1,51 @@ +#pragma once + +//////////////////////////////////////////////////////////////////////////////// +#include "vor3d/Common.h" +#include +#include +#include +#include +//////////////////////////////////////////////////////////////////////////////// + +namespace voroffset3d +{ + + class CompressedVolumeBase + { + protected: + Eigen::Vector3d m_Origin; + Eigen::Vector3d m_Extent; + Scalar m_Spacing; // voxel size (in mm) + int m_Padding; + Eigen::Vector2i m_GridSize; + + public: + // Interface + virtual ~CompressedVolumeBase() = default; + //CompressedVolumeBase() {} + + // Apply a function to each segment in the structure + virtual void iterate(int i, int j, std::function func) = 0; + virtual void iterate(int i, int j, std::function func) = 0; + + // Append a segment to the given segments vector + virtual void appendSegment(int i, int j, Scalar begin_pt, Scalar end_pt, Scalar radius) = 0; + + virtual void reshape(int xsize, int ysize) = 0; + void reset(Eigen::Vector3d origin, Eigen::Vector3d extent, Scalar voxel_size, int padding, + int xsize, int ysize); + virtual void resize(int xsize, int ysize) = 0; + virtual void clear() = 0; + + int numDexels() const { return m_GridSize[0] * m_GridSize[1]; } + Eigen::Vector2d dexelCenter(int x, int y) const; + Eigen::Vector2i gridSize() const { return m_GridSize; } + Eigen::Vector3d origin() const { return m_Origin; } + double spacing() const { return m_Spacing; } + Eigen::Vector3d extent() const { return m_Extent; } + int padding() const { return m_Padding; } + + }; + +} // namespace voroffset3d \ No newline at end of file diff --git a/src/vor3d/CompressedVolumeWithRadii.cpp b/src/vor3d/CompressedVolumeWithRadii.cpp new file mode 100644 index 0000000..9ee6403 --- /dev/null +++ b/src/vor3d/CompressedVolumeWithRadii.cpp @@ -0,0 +1,112 @@ +//////////////////////////////////////////////////////////////////////////////// +#include "vor3d/CompressedVolumeWithRadii.h" +#include "vor3d/MorphologyOperators.h" +//////////////////////////////////////////////////////////////////////////////// + +using namespace voroffset3d; + +//////////////////////////////////////////////////////////////////////////////// + +CompressedVolumeWithRadii::CompressedVolumeWithRadii( + Eigen::Vector3d origin, Eigen::Vector3d extent, double voxel_size, int padding) +{ + m_Origin = origin; + m_Extent = extent; + m_Spacing = voxel_size; + m_Padding = padding; + m_Origin.array() -= padding * voxel_size; + m_GridSize[0] = (int)(std::ceil(extent[0] / m_Spacing) + 2 * padding); + m_GridSize[1] = (int)(std::ceil(extent[1] / m_Spacing) + 2 * padding); + m_Data.assign(m_GridSize[0] * m_GridSize[1], {}); + +} +// ---------------------------------------------------------------------------- +void CompressedVolumeWithRadii::iterate(int i, int j, std::function func) +{ + auto &m_ray = at(i, j); + for (const SegmentWithRadius s : m_ray) + func(s.y1, s.y2, s.r); + +} + +void CompressedVolumeWithRadii::iterate(int i, int j, std::function func) +{ + auto &m_ray = at(i, j); + for (const SegmentWithRadius s : m_ray) + func(s.y1, s.y2); + +} + +void CompressedVolumeWithRadii::copy_volume_from(const CompressedVolumeWithRadii voxel) +{ + m_Origin = voxel.origin(); + m_Extent = voxel.extent(); + m_GridSize = voxel.gridSize(); + m_Padding = voxel.padding(); + m_Spacing = voxel.spacing(); + m_Data.assign(m_GridSize[0] * m_GridSize[1], {}); + for(int x=0;x> m_Data; + std::vector m_segs_vec; + + public: + // Interface + CompressedVolumeWithRadii(Eigen::Vector3d origin, Eigen::Vector3d extent, double voxel_size, int padding); + CompressedVolumeWithRadii() = default; + + void copy_volume_from(const CompressedVolumeWithRadii voxel); + double get_volume(); + int numSegments(); + const std::vector & at(int x, int y) const { return m_Data[x + m_GridSize[0] * y]; } + std::vector & at(int x, int y) { return m_Data[x + m_GridSize[0] * y]; } + + // Apply a function to each segment in the structure + //void iterate(std::function func) const; + virtual void iterate(int i, int j, std::function func) override; + virtual void iterate(int i, int j, std::function func1) override; + + virtual void appendSegment(int i, int j, Scalar begin_pt, Scalar end_pt, Scalar radius) override; + + virtual void reshape(int xsize, int ysize) override; + virtual void resize(int xsize, int ysize) override; + virtual void clear() override; + + }; + +} // namespace voroffset3d diff --git a/src/vor3d/Dexelize.cpp b/src/vor3d/Dexelize.cpp new file mode 100644 index 0000000..300abb5 --- /dev/null +++ b/src/vor3d/Dexelize.cpp @@ -0,0 +1,352 @@ +#include "vor3d/Dexelize.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +//#include "tbb/blocked_range.h" +//#include "tbb/parallel_for.h" +//#include "tbb/task_scheduler_init.h" + +//////////////////////////////////////////////////////////////////////////////// +#ifndef TEST_WID +#define TEST_WID 25 +#endif +#ifndef TEST_LEN +#define TEST_LEN 30 +#endif +#ifndef TEST_HEI +#define TEST_HEI 40 +#endif + +namespace +{ + + template + inline std::ostream& operator<<(std::ostream &out, std::array v) + { + out << "{"; + if (!v.empty()) + { + std::copy(v.begin(), v.end() - 1, std::ostream_iterator(out, "; ")); + out << v.back(); + } + out << "}"; + return out; + } + + //////////////////////////////////////////////////////////////////////////////// + // NOTE: Function `point_in_triangle_2d` comes from SDFGen by Christopher Batty. + // https://github.com/christopherbatty/SDFGen/blob/master/makelevelset3.cpp + //////////////////////////////////////////////////////////////////////////////// + + // calculate twice signed area of triangle (0,0)-(x1,y1)-(x2,y2) + // return an SOS-determined sign (-1, +1, or 0 only if it's a truly degenerate triangle) + int orientation( + double x1, double y1, double x2, double y2, double &twice_signed_area) + { + twice_signed_area = y1 * x2 - x1 * y2; + if (twice_signed_area > 0) return 1; + else if (twice_signed_area < 0) return -1; + else if (y2 > y1) return 1; + else if (y2 < y1) return -1; + else if (x1 > x2) return 1; + else if (x1 < x2) return -1; + else return 0; // only true when x1==x2 and y1==y2 + } + + // ----------------------------------------------------------------------------- + + // robust test of (x0,y0) in the triangle (x1,y1)-(x2,y2)-(x3,y3) + // if true is returned, the barycentric coordinates are set in a,b,c. + bool point_in_triangle_2d( + double x0, double y0, double x1, double y1, + double x2, double y2, double x3, double y3, + double &a, double &b, double &c) + { + x1 -= x0; x2 -= x0; x3 -= x0; + y1 -= y0; y2 -= y0; y3 -= y0; + int signa = orientation(x2, y2, x3, y3, a); + if (signa == 0) return false; + int signb = orientation(x3, y3, x1, y1, b); + if (signb != signa) return false; + int signc = orientation(x1, y1, x2, y2, c); + if (signc != signa) return false; + double sum = a + b + c; + geo_assert(sum != 0); // if the SOS signs match and are nonzero, there's no way all of a, b, and c are zero. + a /= sum; + b /= sum; + c /= sum; + return true; + } + + // ----------------------------------------------------------------------------- + + // \brief Computes the (approximate) orientation predicate in 2d. + // \details Computes the sign of the (approximate) signed volume of + // the triangle p0, p1, p2 + // \param[in] p0 first vertex of the triangle + // \param[in] p1 second vertex of the triangle + // \param[in] p2 third vertex of the triangle + // \retval POSITIVE if the triangle is oriented positively + // \retval ZERO if the triangle is flat + // \retval NEGATIVE if the triangle is oriented negatively + // \todo check whether orientation is inverted as compared to + // Shewchuk's version. + inline GEO::Sign orient_2d_inexact(GEO::vec2 p0, GEO::vec2 p1, GEO::vec2 p2) + { + double a11 = p1[0] - p0[0]; + double a12 = p1[1] - p0[1]; + + double a21 = p2[0] - p0[0]; + double a22 = p2[1] - p0[1]; + + double Delta = GEO::det2x2( + a11, a12, + a21, a22 + ); + + return GEO::geo_sgn(Delta); + } + + //////////////////////////////////////////////////////////////////////////////// + + + /** + * @brief { Intersect a vertical ray with a triangle } + * + * @param[in] M { Mesh containing the triangle to intersect } + * @param[in] f { Index of the facet to intersect } + * @param[in] q { Query point (only XY coordinates are used) } + * @param[out] z { Intersection } + * + * @return { description_of_the_return_value } + */ + template + int intersect_ray_z(const GEO::Mesh &M, GEO::index_t f, const GEO::vec3 &q, double &z) + { + using namespace GEO; + + index_t c = M.facets.corners_begin(f); + const vec3& p1 = Geom::mesh_vertex(M, M.facet_corners.vertex(c++)); + const vec3& p2 = Geom::mesh_vertex(M, M.facet_corners.vertex(c++)); + const vec3& p3 = Geom::mesh_vertex(M, M.facet_corners.vertex(c)); + + double u, v, w; + if (point_in_triangle_2d( + q[X], q[Y], p1[X], p1[Y], p2[X], p2[Y], p3[X], p3[Y], u, v, w)) + { + z = u * p1[Z] + v * p2[Z] + w * p3[Z]; + auto sign = orient_2d_inexact(vec2(p1[X], p1[Y]), vec2(p2[X], p2[Y]), vec2(p3[X], p3[Y])); + switch (sign) + { + case GEO::POSITIVE: return 1; + case GEO::NEGATIVE: return -1; + case GEO::ZERO: + default: return 0; + } + } + + return 0; + } + + // ----------------------------------------------------------------------------- + + void compute_sign(const GEO::Mesh &M, + const GEO::MeshFacetsAABB &aabb_tree, vor3d::CompressedVolume &dexels) + { + auto size = dexels.gridSize(); + + try + { + GEO::ProgressTask task("Ray marching", 100); + + GEO::vec3 min_corner, max_corner; + GEO::get_bbox(M, &min_corner[0], &max_corner[0]); + + const GEO::vec3 origin(dexels.origin().data()); + const double origin_z = origin[2]; + const double spacing = dexels.spacing(); + //GEO::parallel_for([&](int y) + for (int y = 0; y < size[1]; ++y) + { + //if (GEO::Thread::current()->id() == 0) { + // task.progress((int) (100.0 * y / size[1] * GEO::Process::number_of_cores())); + //} + for (int x = 0; x < size[0]; ++x) + { + GEO::vec2 center_xy(dexels.dexelCenter(x, y).data()); + GEO::vec3 center(center_xy[0], center_xy[1], 0); + + GEO::Box box; + box.xyz_min[0] = box.xyz_max[0] = center[0]; + box.xyz_min[1] = box.xyz_max[1] = center[1]; + box.xyz_min[2] = min_corner[2] - spacing; + box.xyz_max[2] = max_corner[2] + spacing; + + std::vector inter; + auto action = [&M, &inter, ¢er, &spacing](GEO::index_t f) + { + double z; + if (intersect_ray_z(M, f, center, z)) + { + inter.push_back(z / spacing); + } + }; + aabb_tree.compute_bbox_facet_bbox_intersections(box, action); + std::sort(inter.begin(), inter.end()); + + dexels.at(x, y).resize(inter.size()); + //dexels.at(x, y).resize(inter.size() / 2); + size_t n = dexels.at(x, y).size(); + std::copy_n(inter.begin(), inter.size(), dexels.at(x, y).begin()); + /*for (int i = 0; i < n; i++) + { + dexels.at(x, y)[i] = vor3d::m_Segment(inter[2 * i], inter[2 * i + 1], radius); + }*/ + } + }//, 0, size[1]); + } + catch (const GEO::TaskCanceled&) + { + // Do early cleanup + } + } + +} + +//////////////////////////////////////////////////////////////////////////////// + +voroffset3d::CompressedVolume voroffset3d::create_dexels( + const std::string &filename, double &voxel_size, int padding, int num_voxels) +{ + // Initialize the Geogram library + GEO::initialize(); + + // Import standard command line arguments, and custom ones + GEO::CmdLine::import_arg_group("standard"); + + // Declare a mesh + GEO::Mesh M; + + // Load the mesh and display timings + GEO::Logger::div("Loading"); + { + GEO::Stopwatch W("Load"); + if (!GEO::mesh_load(filename, M)) + { + throw std::runtime_error("Invalid input mesh."); + } + geo_assert(M.vertices.dimension() == 3); + } + + // Initialize voxel grid and AABB tree + GEO::vec3 min_corner, max_corner; + GEO::get_bbox(M, &min_corner[0], &max_corner[0]); + GEO::vec3 extent = max_corner - min_corner; + if (num_voxels > 0) + { + // Force number of voxels along longest axis + double max_extent = std::max(extent[0], std::max(extent[1], extent[2])); + voxel_size = max_extent / num_voxels; + } + GEO::MeshFacetsAABB aabb_tree(M); + + // Dexelize the input mesh + GEO::Logger::div("Dexelizing"); + CompressedVolume dexels( + Eigen::Vector3d(min_corner[0], min_corner[1], min_corner[2]), + Eigen::Vector3d(extent[0], extent[1], extent[2]), + voxel_size, padding); + compute_sign(M, aabb_tree, dexels); + return dexels; +} + +//////////////////////////////////////////////////////////////////////////////// +// NOTE: Function `dexel_dump` comes from SDFGen by Christopher Batty. +// https://github.com/jdumas/geotools/blob/master/voxmesh +//////////////////////////////////////////////////////////////////////////////// + +bool endswith(const std::string &str, const std::string &suffix) { + if (str.length() >= suffix.length()) { + return (0 == str.compare(str.length() - suffix.length(), suffix.length(), suffix)); + } else { + return false; + } +} + +void points_dump(const std::string &filename, const voroffset3d::CompressedVolume &dexels) +{ + GEO::Mesh mesh; + + for (int y = 0; y < dexels.gridSize()[1]; ++y) { + for (int x = 0; x < dexels.gridSize()[0]; ++x) { + for (int i = 0; 2 * i < dexels.at(x, y).size(); ++i) { + GEO::vec3 xyz_min, xyz_max; + xyz_min[0] = dexels.origin()[0] + (x + 0.5) * dexels.spacing(); + xyz_min[1] = dexels.origin()[1] + (y + 0.5) * dexels.spacing(); + xyz_min[2] = dexels.at(x, y)[2 * i + 0] * dexels.spacing(); + xyz_max[0] = dexels.origin()[0] + (x + 0.5) * dexels.spacing(); + xyz_max[1] = dexels.origin()[1] + (y + 0.5) * dexels.spacing(); + xyz_max[2] = dexels.at(x, y)[2 * i + 1] * dexels.spacing(); + mesh.vertices.create_vertex(&xyz_min[0]); + mesh.vertices.create_vertex(&xyz_max[0]); + } + } + } + + GEO::mesh_save(mesh, filename); +} + +void voroffset3d::dexel_dump(const std::string &filename, const CompressedVolume &dexels) +{ + if (endswith(filename, ".xyz")) { + points_dump(filename, dexels); + return; + } + GEO::Mesh mesh; + + for (int y = 0; y < dexels.gridSize()[1]; ++y) { + for (int x = 0; x < dexels.gridSize()[0]; ++x) { + for (int i = 0; 2 * i < dexels.at(x, y).size(); ++i) { + GEO::vec3 xyz_min, xyz_max; + xyz_min[0] = dexels.origin()[0] + x * dexels.spacing(); + xyz_min[1] = dexels.origin()[1] + y * dexels.spacing(); + xyz_min[2] = dexels.at(x, y)[2 * i + 0] * dexels.spacing(); + xyz_max[0] = dexels.origin()[0] + (x + 1) * dexels.spacing(); + xyz_max[1] = dexels.origin()[1] + (y + 1) * dexels.spacing(); + xyz_max[2] = dexels.at(x, y)[2 * i + 1] * dexels.spacing(); + GEO::vec3 diff[8] = { + GEO::vec3(0,0,1), GEO::vec3(1,0,1), GEO::vec3(0,1,1), GEO::vec3(1,1,1), + GEO::vec3(0,0,0), GEO::vec3(1,0,0), GEO::vec3(0,1,0), GEO::vec3(1,1,0), + }; + int v = mesh.vertices.nb(); + for (int lv = 0; lv < 8; ++lv) { + for (int d = 0; d < 3; ++d) { + diff[lv][d] = xyz_min[d] + diff[lv][d] * (xyz_max[d] - xyz_min[d]); + } + mesh.vertices.create_vertex(&diff[lv][0]); + } + mesh.cells.create_hex(v, v + 1, v + 2, v + 3, v + 4, v + 5, v + 6, v + 7); + } + } + } + + mesh.cells.compute_borders(); + mesh.cells.connect(); + mesh.vertices.remove_isolated(); + + GEO::mesh_save(mesh, filename); +} + diff --git a/src/vor3d/Dexelize.h b/src/vor3d/Dexelize.h new file mode 100644 index 0000000..9ea6dc8 --- /dev/null +++ b/src/vor3d/Dexelize.h @@ -0,0 +1,30 @@ +#pragma once + +//////////////////////////////////////////////////////////////////////////////// +#include "vor3d/Common.h" +#include "vor3d/CompressedVolume.h" +#include +#include +#include +#include +//////////////////////////////////////////////////////////////////////////////// + +namespace voroffset3d +{ + +/** + * @brief Creates dexels from a triangle mesh. + * + * @param[in] filename { Filename of the mesh to load. } + * @param[in,out] voxel_size { Voxel size for the 2D dexel grid. } + * @param[in] padding { Additional padding on each side. } + * @param[in] num_voxels { Explicitly set the 2D grid size (max length), before padding. } + * + * @return { The dexelized volume. } + */ +CompressedVolume create_dexels(const std::string &filename, + double &voxel_size, int padding = 0, int num_voxels = -1); + +void dexel_dump(const std::string &filename, const CompressedVolume &voxels); + +} // namespace voroffset3d diff --git a/src/vor3d/HalfDilationOperator.cpp b/src/vor3d/HalfDilationOperator.cpp new file mode 100644 index 0000000..b37510e --- /dev/null +++ b/src/vor3d/HalfDilationOperator.cpp @@ -0,0 +1,31 @@ +#include "vor3d/HalfDilationOperator.h" +#include "vor3d/MorphologyOperators.h" +namespace voroffset3d +{ + // Forward sweep for dilation operation + void halfDilate( + Morpho2D &vor, + bool is_apply_power_alg, + CompressedVolumeBase &input, + CompressedVolumeBase &output, + int x0, int y0, int deltaX, int deltaY) + { + //output.clear(); + for (int i = x0, j = y0, k = 0; i < input.gridSize()(0) && i >= 0 && j < input.gridSize()(1) && j >= 0; i += deltaX, j += deltaY, ++k) + { + vor.removeInactiveSegments(k); + vor.removeInactivePoints(k); + input.iterate(i, j, [&vor, k](double z1, double z2, double r) + { + vor.insertSegment(k, z1, z2, r); + }); + vor.flushLine(i); + vor.getLine(is_apply_power_alg, x0, y0, deltaX, deltaY, k, [&output](int posX, int posY, double z1, double z2, double r) + { + output.appendSegment(posX, posY, z1, z2, r); + }); + } + + } + +} \ No newline at end of file diff --git a/src/vor3d/HalfDilationOperator.h b/src/vor3d/HalfDilationOperator.h new file mode 100644 index 0000000..cd1bbb7 --- /dev/null +++ b/src/vor3d/HalfDilationOperator.h @@ -0,0 +1,32 @@ +#pragma once +#include "vor3d/CompressedVolumeBase.h" +#include "vor3d/CompressedVolume.h" +#include "vor3d/CompressedVolumeWithRadii.h" +#include "vor3d/Morpho2D.h" + +namespace voroffset3d +{ + + + /** + * @brief Forward sweep for dilation operation. + * + * @param[in] vor { The algorithm for computing the dilation. } + * @param[in] is_apply_power_alg { Whether we apply power algorithm, which means that we need to change the dilation radius. } + * @param[in] input { Input data. } + * @param[in] output { Output data. } + * @param[in] (x0,y0,deltaX, deltaY) + { Arguments for locate the plane as well as decide the sweepline direction. } + **/ + void halfDilate( + Morpho2D &vor, + bool is_apply_power_alg, + CompressedVolumeBase &input, + CompressedVolumeBase &output, + int x0, int y0, int deltaX, int deltaY); + + template + void unionMap(CompressedVolumeType voxel_1, CompressedVolumeType voxel_2, CompressedVolumeType &result); +} + +#include "vor3d/HalfDilationOperator.hpp" \ No newline at end of file diff --git a/src/vor3d/HalfDilationOperator.hpp b/src/vor3d/HalfDilationOperator.hpp new file mode 100644 index 0000000..9464a12 --- /dev/null +++ b/src/vor3d/HalfDilationOperator.hpp @@ -0,0 +1,16 @@ +#include "vor3d/HalfDilationOperator.h" + +namespace voroffset3d +{ + template + void unionMap(CompressedVolumeType voxel_1, CompressedVolumeType voxel_2, CompressedVolumeType &result) + { + int x_size = voxel_1.gridSize()(0); + int y_size = voxel_1.gridSize()(1); + for (int i = 0; i < x_size; i++) + for (int j = 0; j < y_size; j++) + { + unionSegs(voxel_1.at(i, j), voxel_2.at(i, j), result.at(i, j)); + } + } +} \ No newline at end of file diff --git a/src/vor3d/Morpho2D.cpp b/src/vor3d/Morpho2D.cpp new file mode 100644 index 0000000..4fd6929 --- /dev/null +++ b/src/vor3d/Morpho2D.cpp @@ -0,0 +1 @@ +#include "vor3d/Morpho2D.h" \ No newline at end of file diff --git a/src/vor3d/Morpho2D.h b/src/vor3d/Morpho2D.h new file mode 100644 index 0000000..cf2cc2d --- /dev/null +++ b/src/vor3d/Morpho2D.h @@ -0,0 +1,52 @@ +#pragma once + +//////////////////////////////////////////////////////////////////////////////// +#include "vor3d/Common.h" +#include +#include +#include +#include +//////////////////////////////////////////////////////////////////////////////// + + +namespace voroffset3d +{ + class Morpho2D + { + public: + // Constant parameters for the current dilation + int m_XMax; + double m_YMax, m_YMin; + double m_Radius; + double m_DexelSize; + + public: + virtual ~Morpho2D() = default; + + ///////////// + // Methods // + ///////////// + public: + // Insert a new seed segment + virtual void insertSegment(int i, double j1, double j2, double r) = 0; + + // Remove seeds that are not contributing anymore to the current sweep line + virtual void removeInactiveSegments(int i) = 0; + + // Remove seeds that are not contributing anymore to the current sweep line + virtual void removeInactivePoints(int i) {} + + virtual void flushLine(int i) {} + + // Extract the result for the current line + virtual void getLine(bool is_set_radii, int posX, int posY, int deltaX, int deltaY, int i, + std::function _appendSegment)=0; + + // Reset data + virtual void resetData() = 0; + }; + + + //////////////////////////////////////////////////////////////////////////////// + +} // namespace voroffset3d diff --git a/src/vor3d/MorphologyOperators.cpp b/src/vor3d/MorphologyOperators.cpp new file mode 100644 index 0000000..a539672 --- /dev/null +++ b/src/vor3d/MorphologyOperators.cpp @@ -0,0 +1,390 @@ +#include"vor3d/MorphologyOperators.h" +using namespace voroffset3d; + +#ifndef MIN_SEG_SIZE +#define MIN_SEG_SIZE 1e-10 +#endif + + + +namespace voroffset3d +{ + + // Appends segment [a,b] to the given sorted line + // Assume that inserted segment is larger than all of the segments in the segments vector + void appendSegment(std::vector &nl, double a, double b) + { + // Ensures nl[-2] < a + while (!nl.empty() && nl.rbegin()[1] >= a) + { + b = std::max(b, nl.back()); + nl.pop_back(); + nl.pop_back(); + } + if (nl.empty() || nl.back() < a) + { + nl.push_back(a); + nl.push_back(b); + } + else if (nl.back() < b) + { + nl.back() = b; + } + vor_assert(nl.size() % 2 == 0); + } + + // Appends segment [a,b] to the given sorted line + void appendSegment(std::vector &n1, SegmentWithRadius seg) + { + std::vector _segs_vec; + while (!n1.empty() && n1.back().y1 >= seg.y1) + { + if (n1.back().y1 > seg.y2) + { + _segs_vec.push_back(n1.back()); + } + else + { + if (n1.back().y2 >= seg.y2) + { + if (n1.back().r > seg.r) + { + seg.y2 = n1.back().y1; + _segs_vec.push_back(n1.back()); + } + else if (n1.back().r == seg.r) + { + seg.y2 = n1.back().y2; + } + else + { + n1.back().y1 = seg.y2; + _segs_vec.push_back(n1.back()); + } + } + else + { + if (n1.back().r > seg.r) + { + _segs_vec.push_back(vor3d::SegmentWithRadius(n1.back().y2, seg.y2, seg.r)); + _segs_vec.push_back(n1.back()); + seg.y2 = n1.back().y1; + } + } + } + n1.pop_back(); + } + if (n1.empty() || n1.back().y2 < seg.y1) + n1.push_back(seg); + else if (n1.back().y2 < seg.y2) + { + if (n1.back().r == seg.r) + n1.back().y2 = seg.y2; + else if (n1.back().r > seg.r) + { + seg.y1 = n1.back().y2; + _segs_vec.push_back(seg); + } + else + { + n1.back().y2 = seg.y1; + _segs_vec.push_back(seg); + } + } + n1.insert(n1.end(), _segs_vec.rbegin(), _segs_vec.rend()); + } + + // get the next segment of the given position + void getNextSegment(const std::vector & segs, std::vector::const_iterator & it, + double & left, double & right, bool & full) + { + if (it == segs.end()) + { + full = false; + } + else + { + left = *it++; + right = *it++; + } + } + + //////////////////////////////////////////////////////////////////////////////////// + + // Computes the union of two sorted and nonoverlapping lists of segments. + // output is also sorted and nonoverlapping segments + // Uses a parallel sweep of both lists, akin to merge-sort. + void unionSegs(const std::vector & a, const std::vector & b, std::vector & result) + { + if (result.size()) + result.clear(); + if (a.empty()) + { + result.insert(result.end(), b.begin(), b.end()); + return; + } + if (b.empty()) + { + result.insert(result.end(), a.begin(), a.end()); + return; + } + std::vector::const_iterator ia(a.begin()), ib(b.begin()); + double al, ar, bl, br; + al = *ia++; ar = *ia++; + bl = *ib++; br = *ib++; + bool bfull(true), afull(true); + while (bfull || afull) + { + if (bfull && bl <= al) + { + addSegmentAtTheEnd(result, bl, br); + getNextSegment(b, ib, bl, br, bfull); + } + else if (afull) + { + addSegmentAtTheEnd(result, al, ar); + if (ia == a.end()) + { + afull = false; + al = std::numeric_limits::max(); + } + else + { + al = *ia++; + ar = *ia++; + } + } + } + } + + + void unionPoints(const std::vector & a, const std::vector & b + , std::vector & result) + { + if (result.size()) + result.clear(); + if (a.empty()) + { + result.insert(result.end(), b.begin(), b.end()); + return; + } + if (b.empty()) + { + result.insert(result.end(), a.begin(), a.end()); + return; + } + std::vector::const_iterator ia(a.begin()), ib(b.begin()); + PointWithRadius a_1, b_1; + a_1 = *ia++; + b_1 = *ib++; + bool bfull(true), afull(true); + while (bfull || afull) + { + if (bfull && b_1._pt <= a_1._pt) + { + result.push_back(b_1); + if (ib == b.end()) + { + bfull = false; + b_1._pt = std::numeric_limits::max(); + } + else + { + b_1 = *ib++; + } + } + else if (afull) + { + result.push_back(a_1); + if (ia == a.end()) + { + afull = false; + a_1._pt = std::numeric_limits::max(); + } + else + { + a_1 = *ia++; + } + } + } + } + + ////////////////////////////////////////////////////////////////////////////// + + // Unoin map + /* + + template + void voroffset::unionMap(Iterator1 a, Iterator2 b, std::vector &res) + { + for (size_t i = 0; i < res.size(); ++i) + { + SegVector temp; + unionSegs(a[i], b[i], temp); + res[i].swap(temp); + } + }*/ + ////////////////////////////////////////////////////////////////////////////////// + + // negate ray + void negate_ray(std::vector &result, double y_min, double y_max) + { + size_t size_ = result.size(); + + if (result.empty()) + { + result = { y_min, y_max }; + } + else + { + // First event + if (result[0] == y_min) + { + result.erase(result.begin()); + } + else + { + result.insert(result.begin(), y_min); + } + // Last event + if (result.back() == y_max) + { + result.pop_back(); + } + else + { + result.push_back(y_max); + } + } + } + + void negate_ray(std::vector input, std::vector &result, double y_min, double y_max) + { + size_t size_ = input.size(); + result.clear(); + if (input.empty()) + { + result = { SegmentWithRadius(y_min, y_max , 0) }; + } + else + { + // First event + if (input[0].y1 > y_min) + result.push_back(SegmentWithRadius(y_min, input[0].y1, 0)); + for (int i = 0; i < size_ - 1; i++) + result.push_back(SegmentWithRadius(input[i].y2, input[i + 1].y1, 0)); + // End event + if (input[size_ - 1].y2 < y_max) + result.push_back(SegmentWithRadius(input[size_ - 1].y2, y_max, 0)); + } + } + + void negate_ray_range(std::vector &result, double y_min, double y_max) + { + size_t size_ = result.size(); + + if (result.empty()) + { + result = { y_min, y_max }; + } + else + { + int count_first = 0; + int count_last = 0; + while (result[0] <= y_min) + { + result.erase(result.begin()); + count_first++; + if(result.empty()) + break; + } + if (count_first % 2 == 0) + result.insert(result.begin(), y_min); + + while (result.back() >= y_max) + { + result.pop_back(); + count_last++; + if(result.empty()) + break; + } + if (count_last % 2 == 0) + result.push_back(y_max); + + } + } + //////////////////////////////////////////////////////////////////////////////////////// + + // calculate the xor between two vectors of segments + void calculate_ray_xor(std::vector ray_1, std::vector ray_2 + , std::vector &ray_xor, double y_min, double y_max) + { + ray_xor.clear(); + std::vector record_1, record_2, tmp_union_1, tmp_union_2, before_remove_pts; + negate_ray(ray_1, record_1, y_min, y_max); + negate_ray(ray_2, record_2, y_min, y_max); + unionSegs(record_1, ray_2, tmp_union_1); + unionSegs(record_2, ray_1, tmp_union_2); + negate_ray(tmp_union_1, ray_1, y_min, y_max); + negate_ray(tmp_union_2, ray_2, y_min, y_max); + unionSegs(ray_1, ray_2, before_remove_pts); + removepoint(before_remove_pts, ray_xor); + } + + void calculate_ray_xor(std::vector ray_1, std::vector ray_2, std::vector &ray_xor, double y_min, double y_max) + { + ray_xor.clear(); + std::vector record_1, record_2, tmp_union_1, tmp_union_2, before_remove_pts; + record_1.resize(ray_1.size()); + record_2.resize(ray_2.size()); + std::copy_n(ray_1.begin(), ray_1.size(), record_1.begin()); + std::copy_n(ray_2.begin(), ray_2.size(), record_2.begin()); + negate_ray(record_1, y_min, y_max); + negate_ray(record_2, y_min, y_max); + unionSegs(record_1, ray_2, tmp_union_1); + unionSegs(record_2, ray_1, tmp_union_2); + negate_ray(tmp_union_1, y_min, y_max); + negate_ray(tmp_union_2, y_min, y_max); + unionSegs(tmp_union_1, tmp_union_2, before_remove_pts); + removepoint(before_remove_pts, ray_xor); + } + ///////////////////////////////////////////////////////////////////////////// + + // remove the segment when its size is too small to be seen + void removepoint(const std::vector segs, std::vector &res) + { + size_t size_ = segs.size(); + for (int i = 0; i < size_; i++) + { + if (segs[i].y2 - segs[i].y1 < MIN_SEG_SIZE) + continue; + res.push_back(segs[i]); + } + } + void removepoint(const std::vector segs, std::vector &res) + { + size_t size_ = segs.size(); + for (int i = 0; i < size_; i += 2) + { + if (segs[i + 1] - segs[i] < MIN_SEG_SIZE) + continue; + res.push_back(segs[i]); + res.push_back(segs[i + 1]); + } + } + + + int index3dToIndex2d(int i, int j, int xsize, int ysize) + { + vor_assert(i < xsize && j < ysize && i >= 0 && j >= 0); + return i + j*xsize; + } + + Eigen::Vector2i index2dToIndex3d(int index, int xsize, int ysize) + { + vor_assert(index < xsize*ysize && index >= 0); + Eigen::Vector2i res(index%xsize, index / xsize); + return res; + } + +} \ No newline at end of file diff --git a/src/vor3d/MorphologyOperators.h b/src/vor3d/MorphologyOperators.h new file mode 100644 index 0000000..db3c9a9 --- /dev/null +++ b/src/vor3d/MorphologyOperators.h @@ -0,0 +1,53 @@ +#pragma once +#include "vor3d/Common.h" +// Dealing with the operation on the same sweepline, such as union segments, negate ray, append segment to the line and so on + +namespace voroffset3d +{ + // Appends segment [a,b] to the given sorted line + void appendSegment(std::vector &nl, double a, double b); + void appendSegment(std::vector &n1, SegmentWithRadius seg); + + // add segment at the end of the sorted segment vector + template + void addSegmentAtTheEnd(Container & dst, double start, double end); + template + void addSegmentAtTheEnd(Container & dst, SegmentType segment); + + + // get the next segment of the given position + void getNextSegment(const std::vector & segs, std::vector::const_iterator & it, + double & left, double & right, bool & full); + + // Computes the union of two sorted lists of segments. + // Uses a parallel sweep of both lists, akin to merge-sort. + void unionSegs(const std::vector & a, const std::vector & b, std::vector & result); + template + void unionSegs(const std::vector & a, const std::vector & b, std::vector & result); + + void unionPoints(const std::vector & a, const std::vector & b, std::vector & result); + + // negate ray which is occluded by [y_min, y_max] + void negate_ray(std::vector &result, double y_min, double y_max); + void negate_ray(std::vector input, std::vector &result, double y_min, double y_max); + + // negate ray based on the range [y_min, y_max], i.e. delete the segments which don't overlap with [y_min, y_max] + // as well as negate the segments which are occluded by [y_min, y_max] + void negate_ray_range(std::vector &result, double y_min, double y_max); + + // calculate the xor of two rays + void calculate_ray_xor(std::vector ray_1, std::vector ray_2, + std::vector &ray_xor, double y_min, double y_max); + void calculate_ray_xor(std::vector ray_1, std::vector ray_2, std::vector &ray_xor, double y_min, double y_max); + + // remove the segment when its size is too small to be seen, which can modify the viwer of visualization + void removepoint(const std::vector segs, std::vector &res); + void removepoint(const std::vector segs, std::vector &res); + + // coordinate convertion + int index3dToIndex2d(int i, int j, int xsize, int ysize); + Eigen::Vector2i index2dToIndex3d(int index2d, int xsize, int ysize); + +}//namespace voroffset + +#include "MorphologyOperators.hpp" diff --git a/src/vor3d/MorphologyOperators.hpp b/src/vor3d/MorphologyOperators.hpp new file mode 100644 index 0000000..e0d7204 --- /dev/null +++ b/src/vor3d/MorphologyOperators.hpp @@ -0,0 +1,130 @@ +#include"vor3d/MorphologyOperators.h" + + +namespace voroffset3d +{ + // add segment at the end of the sorted segment vector + template + void addSegmentAtTheEnd(Container & dst, double start, double end) + { + if (dst.empty() || (dst.back() < start)) + { + dst.push_back(start); + dst.push_back(end); + } + else if (end > dst.back()) + { + dst.back() = end; + } + } + + template + void addSegmentAtTheEnd(Container & dst, SegmentType segment) + { + if (dst.empty()) + { + dst.push_back(segment); + return; + } + else if (dst.back().isOcclude(segment)) // Attention: we have dst.back().y1<=segment.y1 + { + return; + } + else if (segment.isOcclude(dst.back())) // Attention: we have dst.back().y1<=segment.y1 + { + dst.pop_back(); + dst.push_back(segment); + return; + } + else if (dst.back().y2 < segment.y1) + { + dst.push_back(segment); + return; + } + else if (segment.y2 > dst.back().y2) + { + if (dst.back().isSplit(segment)) + { + dst.back().y2 = segment.y1; + dst.push_back(segment); + return; + } + else + { + segment.y1 = dst.back().y2; + dst.push_back(segment); + return; + } + } + else // segment is occluded by dst.back(), for we have assumed that dst.back()._begin_pt <= segment._begin_pt + { + if (dst.back().isSplit(segment)) + { + SegmentType new_segment(dst.back()); + new_segment.y1 = segment.y2; + dst.back().y2 = segment.y1; + dst.push_back(segment); + dst.push_back(new_segment); + return; + } + } + } + + // Computes the union of two sorted and nonoverlapping lists of segments. + // output is also sorted and nonoverlapping segments + template + void unionSegs(const std::vector & a, const std::vector & b + , std::vector & result) + { + if (result.size()) + result.clear(); + if (a.empty()) + { + result.insert(result.end(), b.begin(), b.end()); + return; + } + if (b.empty()) + { + result.insert(result.end(), a.begin(), a.end()); + return; + } + auto ia = a.begin(); + auto ib = b.begin(); + SegmentType a_1, b_1; + a_1 = *ia++; + b_1 = *ib++; + bool bfull(true), afull(true); + while (bfull || afull) + { + if (bfull && b_1.y1 <= a_1.y1) + { + addSegmentAtTheEnd(result, b_1); + if (ib == b.end()) + { + bfull = false; + b_1.y1 = std::numeric_limits::max(); + } + else + b_1 = *ib++; + } + else if (afull) + { + addSegmentAtTheEnd(result, a_1); + if (ia == a.end()) + { + afull = false; + a_1.y1 = std::numeric_limits::max(); + } + else + a_1 = *ia++; + } + } + // clear the capacity + /*a.clear(); + a.shrink_to_fit(); + b.clear(); + b.shrink_to_fit();*/ + } +} +/////////////////////////////////////////////////////////////////////////////// + diff --git a/src/vor3d/SeparatePower2D.cpp b/src/vor3d/SeparatePower2D.cpp new file mode 100644 index 0000000..ad9be16 --- /dev/null +++ b/src/vor3d/SeparatePower2D.cpp @@ -0,0 +1,546 @@ +//////////////////////////////////////////////////////////////////////////////// +#include "vor3d/SeparatePower2D.h" +#include "vor3d/MorphologyOperators.h" +#include +#include +#include +#include +#include +#include +#include +//////////////////////////////////////////////////////////////////////////////// + +using namespace voroffset3d; + +//////////////////////////////////////////////////////////////////////////////// + +namespace +{ + + // ----------------------------------------------------------------------------- + + template + T sign(T a) + { + return a > 0 ? T(1) : T(-1); + } + + + // Auxiliary methods + template + T det(T a, T b, T c, T d) + { + return a * d - b * c; + } + // Returns true iff a ¡Ü sqrt(x) + template + bool infSqrt(T a, T x) + { + return (a <= 0 || a*a <= x); + } + + // Returns true iff sqrt(x) ¡Ü b + template + bool supSqrt(T b, T x) + { + return (b >= 0 && x <= b*b); + } + + static inline bool ceilIfGreateq(double x, double xmin, double &xdst) + { + if (x >= xmin) + { + xdst = x; + return true; + } + else + { + return false; + } + } + + // Typedefs + typedef SegmentWithRadiusX Segment; + typedef PointWithRadiusX Point; + + //////////////////////////////////////////////////////////////////////////////// + + //////////////////////////////////////////////////////////////////////////////// + + //template + struct Ray + { + typedef double Scalar; + + Scalar a, b, c; + + Ray(Scalar _a, Scalar _b, Scalar _c) + : a(_a), b(_b), c(_c) + {} + + // Perpendicular bisector of p and q, equation is a*x + b*y = c + Ray(Point p, Point q) + : a(2 * (p.x - q.x)) + , b(2 * (p.y - q.y)) + , c(p.x * p.x - q.x * q.x + p.y * p.y - q.y * q.y + q.r * q.r - p.r * p.r) + { + //reduce(a, b, c); + } + + inline double inter(Scalar x) const + { + //return divideDouble(c - a*x, b); + return (c - a*x) / b; + } + + inline bool interBefore(int x, double yr) const + { + // Precondition: b==0 + return ((c - a*x) * sign(b) <= yr * b * sign(b)); + } + + inline bool interAfter(int x, double yl) const + { + // Precondition: b==0 + return (yl * b * sign(b) <= (c - a*x) * sign(b)); + } + }; + +} // anonymous namespace + + //////////////////////////////////////////////////////////////////////////////// + // VoronoiMorpho + //////////////////////////////////////////////////////////////////////////////// + + /* + * Assumes that y_p < y_q < y_r + */ +double SeparatePowerMorpho2D::rayIntersect(PointWithRadiusX p, PointWithRadiusX q, PointWithRadiusX r) const +{ + if (std::abs(p.x - q.x) < EPS && std::abs(q.x - r.x) < EPS) + { + return m_XMax*1.0f; + } + else + { + Ray r1(p, q); + Ray r2(q, r); + if (r1.b*r2.b < 0) + { + r1.a = -r1.a; + r1.b = -r1.b; + r1.c = -r1.c; + } + const double denom = det(r1.a, r2.a, r1.b, r2.b); + //std::cout << r1.a << " " << r1.b << " " << r1.c << " , " << r2.a << " " << r2.b << " " << r2.c << std::endl; + if (std::abs(denom) < EPS) + { + // Intersection at x >= m_XMax are discarded; + return m_XMax; + } + else if (denom > 0) + { + return m_XMax; + } + else + { + const double numer = det(r1.c, r2.c, r1.b, r2.b); + if (numer / denom < q.x + EPS || numer / denom >(q.x + q.r)) + { + return m_XMax; + } + else + { + return ((double)numer / denom); + } + } + } +} + +// ----------------------------------------------------------------------------- + +// Assumes that it != m_S.end() +void SeparatePowerMorpho2D::exploreLeft(P_const_iter it, PointWithRadiusX p, int i) +{ + while (it != m_P.begin()) + { + double xinter = rayIntersect(*std::prev(it), *it, p); + if (xinter < i + EPS) + { + it = std::prev(m_P.erase(it)); + } + else + { + int x_sup = (int)std::floor(xinter) + 1; + if (xinter < m_XMax) + { + if (x_sup < m_XMax) + { + m_Q_P[x_sup].push_back(*it); + } + } + return; + } + } +} + +// Assumes that it != m_S.end() +void SeparatePowerMorpho2D::exploreRight(P_const_iter it, PointWithRadiusX p, int i) +{ + while (std::next(it) != m_P.end()) + { + double xinter = rayIntersect(p, *it, *std::next(it)); + if (xinter < i + EPS) + { + it = m_P.erase(it); + } + else + { + int x_sup = (int)std::floor(xinter) + 1; + if (xinter < m_XMax) + { + if (x_sup < m_XMax) + { + m_Q_P[x_sup].push_back(*it); + } + } + return; + } + } +} + +// ----------------------------------------------------------------------------- + +// Insert a new seed segment (i,j1) -- (i,j2) +void SeparatePowerMorpho2D::insertSegment(int i, double j1, double j2, double r) +{ + if (std::abs(j1 - j2) < EPS) + { + m_P_Tmp.push_back(PointWithRadiusX(i, j1, r)); + return; + } + SegmentWithRadiusX new_seg(i, j1, j2, r); + m_S_Tmp.push_back(new_seg); + m_P_Tmp.push_back(new_seg.left()); + m_P_Tmp.push_back(new_seg.right()); +} + +void SeparatePowerMorpho2D::insertPoint(PointWithRadiusX p) +{ + //PointWithRadiusX p(i, j1, r); + P_const_iter it = m_P.lower_bound(p); + while (it != m_P.end() && p.contains(*it)) + { + it = m_P.erase(it); + } + while (it != m_P.begin() && p.contains(*std::prev(it))) + { + m_P.erase(std::prev(it)); + } + P_const_iter jt = m_P.insert(it, p); + int xsup = (int)std::floor(p.x + p.r) + 1; + if (xsup < m_XMax) { m_Q_P[xsup].push_back(p); } + it = std::next(jt); + if (it != m_P.end()) + { + exploreRight(it, p, (int)p.x); + } + if (jt != m_P.begin()) + { + exploreLeft(std::prev(jt), p, (int)p.x); + } + +} + +// ----------------------------------------------------------------------------- + +// Remove seeds that are not contributing anymore to the current sweep line (ith) +void SeparatePowerMorpho2D::removeInactiveSegments(int i) +{ + // remove inactive segments + int k, j; + k = 0; j = 0; + size_t size = m_S.size(); + for (; k < size; k++) + { + if ((int)std::floor(m_S[k].x + m_S[k].r) + 1 > i) + { + m_S[j] = m_S[k]; + j++; + } + } + m_S.resize(j); +} + +void SeparatePowerMorpho2D::removeInactivePoints(int i) +{ + std::sort(m_Q_P[i].begin(), m_Q_P[i].end()); + while (!m_Q_P[i].empty()) + { + P_const_iter it = m_P.find(m_Q_P[i].back()); + m_Q_P[i].pop_back(); + if (it != m_P.end()) + { + it = m_P.erase(it); + if (it != m_P.begin() && it != m_P.end()) + { + exploreLeft(std::prev(it), *it, i); + exploreRight(it, *std::prev(it), i); + } + } + } + m_Q_P[i].clear(); +} + +void SeparatePowerMorpho2D::getLine(bool is_set_radii, int posX, int posY, int deltaX, int deltaY, int i, + std::function _appendSegment) +{ + int pos_x = posX + i*deltaX; + int pos_y = posY + i*deltaY; + ray_P.clear(); + ray_S.clear(); + ray_U.clear(); + for (SegmentWithRadiusX s : m_S) + { + vor_assert((i - s.x) <= s.r + EPS); + double a, b; + if (i - s.x > s.r) + { + a = s.y1; + b = s.y2; + } + else + { + const double dy = std::sqrt(s.r * s.r - (i - s.x) * (i - s.x)); + a = s.y1 - dy; + b = s.y2 + dy; + } + + appendSegment(ray_S, a, b); //all the radii are same for the first phase + } + for (PointWithRadiusX p : m_P) + { + vor_assert((i - p.x) <= p.r + EPS); + if (i - p.x > p.r) + { + continue; + } + else + { + const double dy = std::sqrt(p.r * p.r - (i - p.x) * (i - p.x)); + double a = p.y - dy; + double b = p.y + dy; + appendSegment(ray_P, a, b); + } + } + vor3d::unionSegs(ray_P, ray_S, ray_U); + for (int k = 0; k < ray_U.size(); k += 2) + { + _appendSegment(pos_x, pos_y, ray_U[k], ray_U[k + 1], 0.f); + } +} + +void SeparatePowerMorpho2D::resetData() +{ + m_S.clear(); + m_S_New.clear(); + m_S_Tmp.clear(); + m_P.clear(); + current_line = -1; +} +void SeparatePowerMorpho2D::unionSegs() +{ + if (m_S_New.size()) + m_S_New.clear(); + if (m_S.empty()) + { + m_S_New.insert(m_S_New.end(), m_S_Tmp.begin(), m_S_Tmp.end()); + return; + } + if (m_S_Tmp.empty()) + { + m_S_New.insert(m_S_New.end(), m_S.begin(), m_S.end()); + return; + } + //std::vector::const_iterator ia(m_S.begin()), ib(m_S_Tmp.begin()); + size_t size_a, size_b; + size_a = m_S.size(); + size_b = m_S_Tmp.size(); + int i, j; + i = 0; j = 0; + SegmentWithRadiusX a_1(m_S[0]), b_1(m_S_Tmp[0]); + i++; + j++; + bool bfull(true), afull(true); + while (bfull || afull) + { + if (afull && a_1.y1 <= b_1.y1) + { + addSegment(a_1); + if (i == size_a) + { + afull = false; + a_1.y1 = std::numeric_limits::max(); + } + else + a_1 = SegmentWithRadiusX(m_S[i++]); + } + else + { + // appendsegment + addSegment(b_1); + if (j == size_b) + { + bfull = false; + b_1.y1 = std::numeric_limits::max(); + } + else + b_1 = SegmentWithRadiusX(m_S_Tmp[j++]); + } + } +} + +void SeparatePowerMorpho2D::addSegment(SegmentWithRadiusX segment) +{ + if (segment.isPoint()) + return; + if (m_S_New.empty()) + { + m_S_New.push_back(segment); + } + else if (m_S_New.back().isOcclude(segment)) // Attention: we have dst.back().y1<=segment.y1 + { + // do nothing + } + else if (m_S_New.back().y2 < segment.y1) + { + m_S_New.push_back(segment); + } + else if (segment.y2 > m_S_New.back().y2) + { + if (m_S_New.back().isSplit(segment)) + { + if (std::abs(m_S_New.back().y1 - segment.y1) < EPS) + m_S_New.pop_back(); + else + { + m_S_New.back().y2 = segment.y1; + } + m_S_New.push_back(segment); + } + else + { + segment.y1 = m_S_New.back().y2; + m_S_New.push_back(segment); + } + } + else // segment is contained by m_S_New.back(), for we have assumed that dst.back()._begin_pt <= segment._begin_pt + { + if (m_S_New.back().isSplit(segment)) + { + + SegmentWithRadiusX new_segment(m_S_New.back()); + new_segment.y1 = segment.y2; + m_S_New.back().y2 = segment.y1; + if (m_S_New.back().isPoint()) + m_S_New.pop_back(); + m_S_New.push_back(segment); + if (!new_segment.isPoint()) + { + m_S_New.push_back(new_segment); + } + } + } +} + +void SeparatePowerMorpho2D::flushLine(int i) +{ + // removing inactive points when a new segment comes + /*for (int k = 0; k < m_S_Tmp.size(); k++) + { + auto it = m_P.lower_bound(m_S_Tmp[k].left()); + if (it != m_P.begin() && std::abs(std::prev(it)->y - m_S_Tmp[k].y1) < EPS) + { + it = std::prev(it); + } + while (it != m_P.end()) + { + if (it->y > m_S_Tmp[k].y2) + break; + double delta_x = m_S_Tmp[k].x - it->x; + double delta_y = std::min(it->y - m_S_Tmp[k].y1, m_S_Tmp[k].y2 - it->y); + if (it->r < std::sqrt(delta_x*delta_x + delta_y*delta_y)) + { + it = m_P.erase(it); + } + else + { + it = std::next(it); + } + } + }*/ + // removing the points lying on the current sweepline + for (int k = 0; k < m_P_Tmp.size(); k++) + { + int pos = locatePoint(m_P_Tmp[k]); + if (pos != -1) + { + double r_m = std::min(std::min(m_S[pos].r - m_P_Tmp[k].x + m_S[pos].x, m_P_Tmp[k].y - m_S[pos].y1), m_S[pos].y2 - m_P_Tmp[k].y); + if (m_P_Tmp[k].r <= r_m) + continue; + else + insertPoint(m_P_Tmp[k]); + } + else + { + insertPoint(m_P_Tmp[k]); + } + } + unionSegs(); + m_S.swap(m_S_New); + m_S_New.clear(); + m_S_Tmp.clear(); + m_P_Tmp.clear(); + +} +int SeparatePowerMorpho2D::locatePoint(PointWithRadiusX p) +/* using binary search to find the first segment in m_S which occuludes the point p + if succeed, return the index of that segment, + otherwise, return -1 +*/ +{ + int left = 0; + int right = (int)m_S.size() - 1; + if (right < 0) + return -1; + int middle; + if (p.y > m_S[right].y2 || p.y < m_S[left].y1) + return -1; + while (left < right - 1) + { + middle = left + ((right - left) >> 1); + + if (m_S[middle].y1 > p.y) + { + right = middle; + } + else if (m_S[middle].y1 < p.y) + { + left = middle; + } + else + { + left = middle; + break; + } + } + if (m_S[left].y2 > p.y) + return left; + else + return -1; +} + + +#ifndef WIN32 +// static_assert(std::is_standard_layout::value, "VoronoiMorpho is not standard layout!"); +#endif \ No newline at end of file diff --git a/src/vor3d/SeparatePower2D.h b/src/vor3d/SeparatePower2D.h new file mode 100644 index 0000000..b239276 --- /dev/null +++ b/src/vor3d/SeparatePower2D.h @@ -0,0 +1,104 @@ +#pragma once + +//////////////////////////////////////////////////////////////////////////////// +#include "vor3d/Common.h" +#include "vor3d/Morpho2D.h" +#include +#include +#include +#include +//////////////////////////////////////////////////////////////////////////////// + +namespace voroffset3d +{ + // the coordinate of our sweep line are range from dexelcenter(0,0) to dexelcenter(gridsize(0),0) if scanning from direction y + // else it should be from dexelcenter(0,0) to dexelcenter(0, gridsize(1)) + struct SeparatePowerMorpho2D : public Morpho2D + { + //////////////// + ////////////////// + // Data members // + ////////////////// + + // Constant parameters for the current dilation + int m_XMax; + double m_YMax, m_YMin; + double m_Radius; + double m_DexelSize; + int current_line; + // Seeds above sorted by ascending y of their midpoint + std::vector m_S; + std::vector m_S_New; + std::vector m_S_Tmp; + std::set m_P; + std::vector m_P_Tmp; + + // Candidate Voronoi vertices sorted by ascending x + std::vector> m_Q_P; + + // Typedefs + typedef std::vector::const_iterator S_const_iter; + typedef std::set::const_iterator P_const_iter; + + ///////////// + // Methods // + ///////////// + + // Assumes that y_p < y_q < y_r + double rayIntersect(PointWithRadiusX p, PointWithRadiusX q, PointWithRadiusX r) const; + + // Assumes that it != m_S.end() + void exploreLeft(P_const_iter it, PointWithRadiusX p, int i); + + // Assumes that it != m_S.end() + void exploreRight(P_const_iter it, PointWithRadiusX p, int i); + + // Insert a new seed segment + virtual void insertSegment(int i, double j1, double j2, double r) override; + void insertPoint(PointWithRadiusX p); + + // Remove seeds that are not contributing anymore to the current sweep line + virtual void removeInactiveSegments(int i) override; + void removeInactivePoints(int i); + + virtual void resetData() override; + + void flushLine(int i); + + int locatePoint(PointWithRadiusX p); + /* using binary search to find the first segment in m_S which occuludes the point p + if succeed, return the index of that segment, + otherwise, return -1 + */ + + // Extract the result for the current line + void getLine(bool is_set_radii, int posX, int posY, int deltaX, int deltaY, int i, + std::function _appendSegment); + + SeparatePowerMorpho2D(int _xmax, double _ymin, double _ymax, double _dexel_size) + : m_XMax(_xmax) + , m_YMax(_ymax) + , m_YMin(_ymin) + , m_DexelSize(_dexel_size) + , m_Q_P(m_XMax) + { + current_line = -1; + } + + SeparatePowerMorpho2D() = default; + + private: + std::vector ray_S; + std::vector ray_P; + std::vector ray_U; + + private: + void unionSegs(); // union segments in m_S and m_S_Tmp, at mean time, we upate m_S + void addSegment(SegmentWithRadiusX segment); // append the segment to the end + + }; + + + //////////////////////////////////////////////////////////////////////////////// + +} // namespace voroffset3d diff --git a/src/vor3d/Timer.cpp b/src/vor3d/Timer.cpp new file mode 100644 index 0000000..b6a7199 --- /dev/null +++ b/src/vor3d/Timer.cpp @@ -0,0 +1,21 @@ +#include "vor3d/Timer.h" + +double Timer::show() +{ + typedef std::chrono::duration double_milli; + auto t = elapsed(); + auto t_s = std::chrono::duration_cast(t); + auto t_ms = std::chrono::duration_cast(t - t_s); + auto total_ms = std::chrono::duration_cast(t); + return std::chrono::duration_cast(t).count(); +} + + + +double Timer::get() +{ + typedef std::chrono::duration double_milli; + auto t = elapsed(); + auto t_ms = std::chrono::duration_cast(t); + return t_ms.count(); +} \ No newline at end of file diff --git a/src/vor3d/Timer.h b/src/vor3d/Timer.h new file mode 100644 index 0000000..43e2345 --- /dev/null +++ b/src/vor3d/Timer.h @@ -0,0 +1,44 @@ +#pragma once +#include +#include + +struct Timer +{ +public: + // Public time point typedef + typedef std::chrono::steady_clock::time_point TimePoint; + + // Compute time difference + template + static auto elapsed(std::chrono::time_point &t1) + { + using namespace std::chrono; + auto t2 = clock_t::now(); + auto time_span = t2 - t1; + t1 = t2; + return time_span; + } + +private: + // Member variables + TimePoint m_CurrentTime; + +public: + // Get current time + static TimePoint now() { return std::chrono::steady_clock::now(); } + + // Constructor + Timer() + : m_CurrentTime(Timer::now()) + { } + + // Compute elapsed time + auto elapsed() { return Timer::elapsed(m_CurrentTime); } + + // Print elapsed time, return elapsed time (in ms) + double show(); + + // Retrieve elapsed time (in ms) + double get(); + +}; \ No newline at end of file diff --git a/src/vor3d/Voronoi.cpp b/src/vor3d/Voronoi.cpp new file mode 100644 index 0000000..537db70 --- /dev/null +++ b/src/vor3d/Voronoi.cpp @@ -0,0 +1,112 @@ +#include"vor3d/Voronoi.h" +#include"vor3d/MorphologyOperators.h" +using namespace voroffset3d; + +#ifndef MY_TYPE_EPS +#define MY_TYPE_EPS 1e-10 +#endif +void VoronoiMorpho::erosion(CompressedVolume input, CompressedVolume &result, double radius, double &time_1, double &time_2) +{ + double z_min = input.origin()(2) / input.spacing() - 1; + double z_max = input.origin()(2) / input.spacing() + 2 * input.padding() + input.extent()(2) / input.spacing() + 1; + negate(input, z_min, z_max); + //result.copy_volume_from(input); + result.reshape(input.gridSize()(0), input.gridSize()(1)); + dilation(input,result,radius, time_1, time_2); + negateInv(result, z_min + 1, z_max - 1); +} +void VoronoiMorpho::negate(CompressedVolume &result, double z_min, double z_max) +{ + int xsize = result.gridSize()(0); + int ysize = result.gridSize()(1); + // First step, add a border to the current grid + result.resize(xsize + 2, ysize + 2); + int read = xsize*ysize - 1; + int write = index3dToIndex2d(xsize, ysize, xsize + 2, ysize + 2); + int count = 0; + for (; read >= 0; read--) + { + count++; + Eigen::Vector2i pre_coor = index2dToIndex3d(read, xsize + 2, ysize + 2); + Eigen::Vector2i cur_coor = index2dToIndex3d(write, xsize + 2, ysize + 2); + result.at(cur_coor(0), cur_coor(1)).resize(result.at(pre_coor(0), pre_coor(1)).size()); + std::copy_n(result.at(pre_coor(0), pre_coor(1)).begin(), result.at(pre_coor(0), pre_coor(1)).size(), + result.at(cur_coor(0), cur_coor(1)).begin()); + if (count == xsize) + { + count = 0; + write -= 3; + } + else + write--; + } + for (int i = 0; i < xsize + 2; i++) + result.at(i, 0).clear(); + for (int j = 0; j < ysize + 2; j++) + result.at(0, j).clear(); + for (int j = 0; j < ysize + 2; j++) + result.at(xsize + 1, j).clear(); + // Second step, negate ray + for (int x = 0; x < result.gridSize()(0); x++) + for (int y = 0; y < result.gridSize()(1); y++) + { + negate_ray(result.at(x, y), z_min, z_max); + } +} + +void VoronoiMorpho::negateInv(CompressedVolume &result, double z_min, double z_max) +{ + int xsize = result.gridSize()(0); + int ysize = result.gridSize()(1); + // First step, add a border to the current grid + int read = index3dToIndex2d(1, 1, xsize, ysize); + int write = 0; + int count = 0; + for (; read <= xsize*ysize - xsize - 2; write++) + { + count++; + Eigen::Vector2i pre_coor = index2dToIndex3d(read, xsize, ysize); + Eigen::Vector2i cur_coor = index2dToIndex3d(write, xsize, ysize); + result.at(cur_coor(0), cur_coor(1)).resize(result.at(pre_coor(0), pre_coor(1)).size()); + std::copy_n(result.at(pre_coor(0), pre_coor(1)).begin(), result.at(pre_coor(0), pre_coor(1)).size(), + result.at(cur_coor(0), cur_coor(1)).begin()); + if (count == xsize - 2) + { + count = 0; + read += 3; + } + else + read++; + + } + result.resize(xsize - 2, ysize - 2); + // Second step, negate ray + for (int x = 0; x < result.gridSize()(0); x++) + for (int y = 0; y < result.gridSize()(1); y++) + { + negate_ray_range(result.at(x, y), z_min, z_max); + } +} + +double VoronoiMorpho::calculateXor(CompressedVolume voxel_1, CompressedVolume voxel_2, CompressedVolume &result) +/** +* @brief calculate the xor between two voxels, with the assumption that these two voxels have the same gridesize +* +* @param[in] voxel_1 { One of voxel. } +* @param[in] voxel_2 { Another voxel. } +* @param[in] result { the result of xor operation. } +* +* @return { The volume of the xor voxel. } +*/ +{ + int x_size = voxel_1.gridSize()(0); + int y_size = voxel_1.gridSize()(1); + double z_min = voxel_1.origin()(2) / voxel_1.spacing(); + double z_max = voxel_1.origin()(2) / voxel_1.spacing() + 2 * voxel_1.padding() + voxel_1.extent()(2) / voxel_1.spacing(); + result.reset(voxel_1.origin(),voxel_1.extent(),voxel_1.spacing(),voxel_1.padding(), x_size, y_size); + for (int x = 0; x < x_size; x++) + for (int y = 0; y < y_size; y++) + calculate_ray_xor(voxel_1.at(x, y), voxel_2.at(x, y), result.at(x, y),z_min,z_max); + return result.get_volume(); +} + diff --git a/src/vor3d/Voronoi.h b/src/vor3d/Voronoi.h new file mode 100644 index 0000000..019c717 --- /dev/null +++ b/src/vor3d/Voronoi.h @@ -0,0 +1,46 @@ +#pragma once +#include "vor3d/Common.h" +#include "vor3d/CompressedVolumeBase.h" +#include "vor3d/CompressedVolume.h" +#include "vor3d/CompressedVolumeWithRadii.h" +#include +#include +#include +#include + +namespace voroffset3d +{ + class VoronoiMorpho + { + + public: + virtual ~VoronoiMorpho() = default; + virtual void dilation(CompressedVolume input, CompressedVolume &result, double radius, double &time_1, double &time_2) = 0; + /** + * @brief dilate the input with given radius + * + * @param[in] input { Input data. } + * @param[in] result { Output data. } + * @param[in] radius { Dilated radius. } + * @param[in] time_1 { time cost by first pass} + * @param[in] time_2 { time cost by second pass} + * + **/ + virtual void erosion(CompressedVolume input, CompressedVolume &result, double radius, double &time_1, double &time_2); + double calculateXor(CompressedVolume voxel_1, CompressedVolume voxel_2, CompressedVolume &result); + /** + * @brief calculate the xor between two voxels, with the assumption that these two voxels have the same gridesize + * + * @param[in] voxel_1 { One of voxel. } + * @param[in] voxel_2 { Another voxel. } + * @param[in] result { the result of xor operation. } + * + * @return { The volume of the xor voxel. } + */ + + protected: + void negate(CompressedVolume &result, double z_min, double z_max); + void negateInv(CompressedVolume &result, double z_min, double z_max); + }; +} + \ No newline at end of file diff --git a/src/vor3d/Voronoi2D.cpp b/src/vor3d/Voronoi2D.cpp new file mode 100644 index 0000000..e1f7508 --- /dev/null +++ b/src/vor3d/Voronoi2D.cpp @@ -0,0 +1,749 @@ +//////////////////////////////////////////////////////////////////////////////// +#include "vor3d/Voronoi2D.h" +#include "vor3d/MorphologyOperators.h" +#include +#include +#include +#include +#include +#include +#include +//////////////////////////////////////////////////////////////////////////////// + +using namespace voroffset3d; +#ifndef TOLENRANCE +#define TOLENRANCE 1e-6 +#endif +//////////////////////////////////////////////////////////////////////////////// + +namespace +{ + + // ----------------------------------------------------------------------------- + + template + T sign(T a) + { + return a > 0 ? T(1) : T(-1); + } + + + // Auxiliary methods + template + T det(T a, T b, T c, T d) + { + return a * d - b * c; + } + // Returns true iff a ¡Ü sqrt(x) + template + bool infSqrt(T a, T x) + { + return (a <= 0 || a*a <= x); + } + + // Returns true iff sqrt(x) ¡Ü b + template + bool supSqrt(T b, T x) + { + return (b >= 0 && x <= b*b); + } + + static inline bool ceilIfGreateq(double x, double xmin, double &xdst) + { + if (x >= xmin) + { + xdst = x; + return true; + } + else + { + return false; + } + } + + //////////////////////////////////////////////////////////////////////////////// + + //template + struct Ray + { + Scalar a, b, c; + + Ray(Scalar _a, Scalar _b, Scalar _c) + : a(_a), b(_b), c(_c) + { } + + // Perpendicular bisector of p and q, equation is a*x + b*y = c + Ray(Point p, Point q) + : a(2 * (p.x - q.x)) + , b(2 * (p.y - q.y)) + , c(p.x * p.x - q.x * q.x + p.y * p.y - q.y * q.y) + { + //reduce(a, b, c); + } + + inline double inter(Scalar x) const + { + //return divideDouble(c - a*x, b); + return (c - a*x) / b; + } + + inline bool interBefore(int x, double yr) const + { + // Precondition: b==0 + return ((c - a*x) * sign(b) <= yr * b * sign(b)); + } + + inline bool interAfter(int x, double yl) const + { + // Precondition: b==0 + return (yl * b * sign(b) <= (c - a*x) * sign(b)); + } + }; + + //////////////////////////////////////////////////////////////////////////////// + + /* + * Assumes that y_p < y_q < y_r + * Returns true if the intersection is between y1 and y2 + */ + bool ray_intersect_range(Point p, Point q, Point r, double y1, double y2, double &xinter) + { + if (std::abs(p.x - q.x) q.x) + { + xinter = ((double)xnumer / denom); + return true; + } + } + return false; + } + } + } + + bool ray_intersect_after(Point p, Point q, Point r, double y, double &xinter) + { + if (std::abs(p.x - q.x) q.x) + { + xinter = ((double)xnumer / denom); + return true; + } + } + return false; + } + } + } + + bool ray_intersect_before(Point p, Point q, Point r, double y, double &xinter) + { + if (std::abs(p.x - q.x) q.x) + { + xinter = ((double)xnumer / denom); + return true; + } + } + return false; + } + } + } + + /* + * Finds a point that is at equal distance of p, q and s + * Assumes that s.x < p.x && s.x < q.x + */ + bool voronoi_vertex(Segment s, Point p, Point q, double &xinter) + { + /* + * We seek the next voronoi point within the range [s.y1, s.y2]. + * Let (xm, ym) be its coordinates. Per the previous sentence we must have: + * (a) xm >= s.x + * (b) s.y1 ¡Ü ym ¡Ü s.y2 + * + * Let r be the bisector of p and q. We suppose that the point we seek is on r. + * As the squared distance from (xm,ym) to s is given by (xm - s.x)^2, we must have: + * + * (xm - s.x)^2 = (xm - xp)^2 + (ym - yp)^2 (1) + * = (xm - xq)^2 + (ym - yq)^2 (2) + * + * The following code solves (1) for ym, and returns true iff (a) and (b) + */ + const Ray r(p, q); + const double u = 2.0 * (p.x - s.x); + const double w = s.x * s.x - p.x * p.x; + double A = r.a; + double B = r.b * u - 2.0 * r.a * p.y; + double C = r.a * p.y * p.y - w * r.a - r.c * u; + //reduce(A, B, C); + const double delta = B * B - 4.0 * A * C; + //const ll alpha = 2 * A * A * s.y1 + 2 * A * B; + //const ll beta = 2 * A * A * s.y2 + 2 * A * B; + const double alpha = 2.0 * std::abs(A) * s.y1 + sign(A) * B; + const double beta = 2.0 * std::abs(A) * s.y2 + sign(A) * B; + const bool pos = (A > 0); + // We must have alpha ¡Ü sign(a) * ¡Àsqrt(delta) ¡Ü beta + if (std::abs(A)= m_XMax are discarded; + return m_XMax*1.0f; + } + else + { + const double numer = det(r1.c, r2.c, r1.b, r2.b); + if (numer /denom < q.x + EPS) + { + return m_XMax*1.0f; + } + else + { + return ((double)numer / denom); + } + } + } +} + +/* +* Assumes that y_lp < y_ab < y_qr +*/ +double VoronoiMorpho2D::treatSegments(Segment lp, Segment ab, Segment qr) const +{ + if (ab.x + EPS > lp.x && ab.x + EPS > qr.x) + { + // Nothing to do + return m_XMax; + } + else if (ab.isPoint()) + { + if (lp.isPoint() && qr.isPoint()) + { + // Case: lp = point, ab = point, qr = point + return rayIntersect(lp.any(), ab.any(), qr.any()); + } + else if (lp.isPoint()) + { + // Case: lp = point, ab = point, qr = segment + double xinter = m_XMax*1.0f; + if (ray_intersect_before(lp.any(), ab.any(), qr.left(), qr.y1, xinter)) + { + return xinter; + } + else if (ray_intersect_after(lp.any(), ab.any(), qr.right(), qr.y2, xinter)) + { + return xinter; + } + else if (voronoi_vertex(qr, lp.any(), ab.any(), xinter)) + { + return xinter; + } + else + { + return xinter; + } + } + else if (qr.isPoint()) + { + // Case: lp = segment, ab = point, qr = point + double xinter = m_XMax*1.0f; + if (ray_intersect_before(lp.left(), ab.any(), qr.any(), lp.y1, xinter)) + { + return xinter; + } + else if (ray_intersect_after(lp.right(), ab.any(), qr.any(), lp.y2, xinter)) + { + return xinter; + } + else if (voronoi_vertex(lp, ab.any(), qr.any(), xinter)) + { + return xinter; + } + else + { + return xinter; + } + } + else + { + // Case: lp = segment, ab = point, qr = segment + double xinter = m_XMax*1.0f; + if (ray_intersect_before(lp.left(), ab.any(), qr.left(), lp.y1, xinter)) + { + return xinter; + } + else if (ray_intersect_after(lp.right(), ab.any(), qr.right(), qr.y2, xinter)) + { + return xinter; + } + else if (ray_intersect_range(lp.right(), ab.any(), qr.left(), lp.y2, qr.y1, xinter)) + { + return xinter; + } + else if (voronoi_vertex(lp, ab.any(), qr.left(), xinter)) + { + return xinter; + } + else if (voronoi_vertex(qr, lp.right(), ab.any(), xinter)) + { + return xinter; + } + else + { + return xinter; + } + } + } + else + { + if (lp.isPoint() && qr.isPoint()) + { + // Case: lp = point, ab = segment, qr = point + double xinter = m_XMax*1.0f; + if (ray_intersect_before(lp.any(), ab.left(), qr.any(), ab.y1, xinter)) + { + return xinter; + } + else if (ray_intersect_after(lp.any(), ab.right(), qr.any(), ab.y2, xinter)) + { + return xinter; + } + else if (voronoi_vertex(ab, lp.any(), qr.any(), xinter)) + { + return xinter; + } + else + { + return xinter; + } + } + else if (lp.isPoint()) + { + // Case: lp = point, ab = segment, qr = segment + double xinter = m_XMax*1.0f; + if (ray_intersect_before(lp.any(), ab.left(), qr.left(), ab.y1, xinter)) + { + return xinter; + } + else if (ray_intersect_after(lp.any(), ab.right(), qr.right(), qr.y2, xinter)) + { + return xinter; + } + else if (ray_intersect_range(lp.any(), ab.right(), qr.left(), ab.y2, qr.y1, xinter)) + { + return xinter; + } + else if (voronoi_vertex(ab, lp.any(), qr.left(), xinter)) + { + return xinter; + } + else if (voronoi_vertex(qr, lp.any(), ab.right(), xinter)) + { + return xinter; + } + else + { + return xinter; + } + } + else if (qr.isPoint()) + { + // Case: lp = segment, ab = segment, qr = point + double xinter = m_XMax*1.0f; + if (ray_intersect_before(lp.left(), ab.left(), qr.any(), lp.y1, xinter)) + { + return xinter; + } + else if (ray_intersect_after(lp.right(), ab.right(), qr.any(), ab.y2, xinter)) + { + return xinter; + } + else if (ray_intersect_range(lp.right(), ab.left(), qr.any(), lp.y2, ab.y1, xinter)) + { + return xinter; + } + else if (voronoi_vertex(ab, lp.right(), qr.any(), xinter)) + { + return xinter; + } + else if (voronoi_vertex(lp, ab.left(), qr.any(), xinter)) + { + return xinter; + } + else + { + return xinter; + } + } + else + { + // Case: lp = segment, ab = segment, qr = segment + double xinter = m_XMax*1.0f; + if (ray_intersect_before(lp.left(), ab.left(), qr.left(), lp.y1, xinter)) + { + return xinter; + } + else if (ray_intersect_after(lp.right(), ab.right(), qr.right(), qr.y2, xinter)) + { + return xinter; + } + else if (ray_intersect_range(lp.right(), ab.left(), qr.left(), lp.y2, ab.y1, xinter)) + { + return xinter; + } + else if (ray_intersect_range(lp.right(), ab.right(), qr.left(), ab.y2, qr.y1, xinter)) + { + return xinter; + } + else if (voronoi_vertex(ab, lp.right(), qr.left(), xinter)) + { + return xinter; + } + else if (voronoi_vertex(lp, ab.left(), qr.left(), xinter)) + { + return xinter; + } + else if (voronoi_vertex(qr, lp.right(), ab.right(), xinter)) + { + return xinter; + } + else + { + return xinter; + } + } + } +} + +// ----------------------------------------------------------------------------- + +// Assumes that it != m_S.end() +void VoronoiMorpho2D::exploreLeft(S_const_iter it, Segment qr, int i) +{ + while (it != m_S.begin()) + { + double xinter = treatSegments(*std::prev(it), *it, qr); + if (xinter < i + EPS) + { + it = std::prev(m_S.erase(it)); + } + else + { + int x_sup = (int)std::floor(xinter) + 1; + if (xinter < m_XMax) + { + if (x_sup < m_XMax) + { + m_Q[x_sup].push_back(*it); + } + } + return; + } + } +} + +// Assumes that it != m_S.end() +void VoronoiMorpho2D::exploreRight(S_const_iter it, Segment lp, int i) +{ + while (std::next(it) != m_S.end()) + { + double xinter = treatSegments(lp, *it, *std::next(it)); + if (xinter < i + EPS) + { + it = m_S.erase(it); + } + else + { + int x_sup = (int)std::floor(xinter) + 1; + if (xinter < m_XMax) + { + if (x_sup < m_XMax) + { + m_Q[x_sup].push_back(*it); + } + } + return; + } + } +} + +// ----------------------------------------------------------------------------- + +// Insert a new seed segment (i,j1) -- (i,j2) +void VoronoiMorpho2D::insertSegment(int i, double j1, double j2, double r) +{ + m_new_segs.clear(); + const Segment s(i, j1, j2); + // 1st step: remove all segments occluded by s + S_const_iter it = m_S.lower_bound(s); + while (it != m_S.end() && s.contains(*it)) + { + it = m_S.erase(it); + } + while (it != m_S.begin() && s.contains(*std::prev(it))) + { + m_S.erase(std::prev(it)); + } + // 2nd step: See if we have to split the neighboring segments on the right + if (it != m_S.end()) + { + if (it->y1 <= s.y2) + { + int xsup = (int)std::floor(it->x + m_Radius)+1; + // Special case: we may have to split the overlapping Segment + if (it->y1 < s.y1) + { + Segment lp(it->x, it->y1, s.y1); + //m_S.insert(it, lp); + //if (xsup < m_XMax) { m_Q[xsup].push_back(lp); } + m_new_segs.push_back(lp); + } + Segment ab(it->x, s.y2, it->y2); + it = m_S.erase(it); + //it = m_S.insert(it, ab); + //if (xsup < m_XMax) { m_Q[xsup].push_back(ab); } + m_new_segs.push_back(ab); + } + } + // And now look at the neighboring segments on the left + if (it != m_S.begin()) + { + S_const_iter jt = std::prev(it); + if (jt->y2 >= s.y1) + { + int xsup = (int)std::floor(jt->x + m_Radius)+1; + // Special case: we may have to split the overlapping Segment + if (jt->y2 > s.y2) + { + Segment lp(jt->x, s.y2, jt->y2); + //it = m_S.insert(it, lp); + //if (xsup < m_XMax) { m_Q[xsup].push_back(lp); } + m_new_segs.push_back(lp); + } + Segment ab(jt->x, jt->y1, s.y1); + jt = m_S.erase(jt); + //m_S.insert(jt, ab); + //if (xsup < m_XMax) { m_Q[xsup].push_back(ab); } + m_new_segs.push_back(ab); + } + } + for (int k = 0; k < m_new_segs.size(); k++) + { + m_S.insert(it, m_new_segs[k]); + int xsup = (int)std::floor(m_new_segs[k].x + m_Radius) + 1; + if (xsup < m_XMax) { m_Q[xsup].push_back(m_new_segs[k]); } + } + // 3rd step: Insert the new segment in the list of seeds + //it = m_S.lower_bound(s); + S_const_iter jt = m_S.insert(it, s); + it = std::next(jt); + int xsup = (int)std::floor(s.x + m_Radius)+1; + if (xsup < m_XMax) { m_Q[xsup].push_back(s); } + // 4th step: Inspect what's on the right and on the left + if (it != m_S.end()) + { + exploreRight(it, s, i); + } + if (jt != m_S.begin()) + { + exploreLeft(std::prev(jt), s, i); + } +} + +// ----------------------------------------------------------------------------- + +// Remove seeds that are not contributing anymore to the current sweep line (ith) +void VoronoiMorpho2D::removeInactiveSegments(int i) +{ + std::sort(m_Q[i].begin(), m_Q[i].end()); + while (!m_Q[i].empty()) + { + S_const_iter it = m_S.find(m_Q[i].back()); + m_Q[i].pop_back(); + if (it != m_S.end()) + { + it = m_S.erase(it); + if (it != m_S.begin() && it != m_S.end()) + { + exploreLeft(std::prev(it), *it, i); + exploreRight(it, *std::prev(it), i); + } + } + } + m_Q[i].clear(); +} + +//////////////////////////////////////////////////////////////////////////////// + + + +// Extract the result for the current line +void VoronoiMorpho2D::getLine(bool is_set_radii, int posX, int posY, int deltaX, int deltaY, int i, + std::function _appendSegment) +{ + if (is_set_radii) + { + for (Segment s : m_S) + { + double a, b; + double r; + double delta_x = i - s.x < m_Radius ? i - s.x : m_Radius; + if (i == s.x) + { + r = m_Radius; + } + else + { + r = std::sqrt(m_Radius * m_Radius - delta_x*delta_x); + } + a = s.y1; + b = s.y2; + int pos_x = posX + i*deltaX; + int pos_y = posY + i*deltaY; + _appendSegment(pos_x, pos_y, a, b, r); + } + } + else + { + for (Segment s : m_S) + { + vor_assert((i - s.x) <= m_Radius + TOLENRANCE*m_DexelSize); + double delta_x = i - s.x < m_Radius ? i - s.x : m_Radius; + const double dy = (std::sqrt(m_Radius * m_Radius - delta_x * delta_x)); + const double a = s.y1 - dy; + const double b = s.y2 + dy; + int pos_x = posX + i*deltaX; + int pos_y = posY + i*deltaY; + _appendSegment(pos_x, pos_y, a, b, 0.f); + } + + } +} + +void VoronoiMorpho2D::resetData() +{ + m_S.clear(); +} + + +#ifndef WIN32 +// static_assert(std::is_standard_layout::value, "VoronoiMorpho is not standard layout!"); +#endif diff --git a/src/vor3d/Voronoi2D.h b/src/vor3d/Voronoi2D.h new file mode 100644 index 0000000..fe8ac3f --- /dev/null +++ b/src/vor3d/Voronoi2D.h @@ -0,0 +1,93 @@ +#pragma once + +//////////////////////////////////////////////////////////////////////////////// +#include "vor3d/Common.h" +#include "vor3d/Morpho2D.h" +#include +#include +#include +#include +//////////////////////////////////////////////////////////////////////////////// + +namespace voroffset3d +{ + // the coordinate of our sweep line are range from dexelcenter(0,0) to dexelcenter(gridsize(0),0) if scanning from direction y + // else it should be from dexelcenter(0,0) to dexelcenter(0, gridsize(1)) + struct VoronoiMorpho2D : public Morpho2D + { + //////////////// + // Data types // + // A horizontal segment + + + //////////////// + ////////////////// + // Data members // + ////////////////// + + //void normalDilation(std::vector> &results); + //void midphaseDilation(std::vector> input, std::vector> &result); + + // Constant parameters for the current dilation + int m_XMax; + double m_YMax, m_YMin; + double m_Radius; + double m_DexelSize; + // Seeds above sorted by ascending y of their midpoint + std::set m_S; + + // Candidate Voronoi vertices sorted by ascending x + std::vector > m_Q; + + // Typedefs + typedef std::set::const_iterator S_const_iter; + + ///////////// + // Methods // + ///////////// + + // Assumes that y_p < y_q < y_r + double rayIntersect(Point p, Point q, Point r) const; + + // Assumes that y_lp < y_ab < y_qr + double treatSegments(Segment lp, Segment ab, Segment qr) const; + + // Assumes that it != m_S.end() + void exploreLeft(S_const_iter it, Segment qr, int i); + + // Assumes that it != m_S.end() + void exploreRight(S_const_iter it, Segment lp, int i); + + // Insert a new seed segment [(i,j1), (i,j2)] + virtual void insertSegment(int i, double j1, double j2, double r) override; + + // Remove seeds that are not contributing anymore to the current sweep line (y==i) + virtual void removeInactiveSegments(int i) override; + + virtual void resetData() override; + + // Extract the result for the current line + void getLine(bool is_set_radii, int posX, int posY, int deltaX, int deltaY, int i, + std::function _appendSegment); + + + VoronoiMorpho2D(int _xmax, double _ymin, double _ymax, double _radius, double _dexel_size) + : m_XMax(_xmax) + , m_YMax(_ymax) + , m_YMin(_ymin) + , m_Radius(_radius) + , m_DexelSize(_dexel_size) + , m_Q(m_XMax) + {} + + VoronoiMorpho2D() + {} + + private: + std::vector m_new_segs; + }; + + + //////////////////////////////////////////////////////////////////////////////// + +} // namespace voroffset3d diff --git a/src/vor3d/VoronoiBruteForce.cpp b/src/vor3d/VoronoiBruteForce.cpp new file mode 100644 index 0000000..a7c9590 --- /dev/null +++ b/src/vor3d/VoronoiBruteForce.cpp @@ -0,0 +1,101 @@ +#include "vor3d/VoronoiBruteForce.h" +#include "vor3d/MorphologyOperators.h" +#ifdef USE_TBB +#include "tbb/blocked_range.h" +#include "tbb/parallel_for.h" +#include "tbb/task_scheduler_init.h" + +#ifndef GRAIN_SIZE +#define GRAIN_SIZE 10 +#endif +#endif + +using namespace voroffset3d; + + +void VoronoiMorphoBruteForce::dilation(CompressedVolume input, CompressedVolume &result, double radius, double &time_1, double &time_2) +{ + int x_size = input.gridSize()(0); + int y_size = input.gridSize()(1); + result.reset(input.origin(), input.extent(), input.spacing(), input.padding(), x_size, y_size); +#ifdef USE_TBB + auto operation = [&](const tbb::blocked_range &range) + { + std::vector tmp_ray; + for (uint32_t phaseIdx = range.begin(); phaseIdx < range.end(); ++phaseIdx) + { + + int x = phaseIdx % x_size; + int y = phaseIdx / x_size; +#else + std::vector tmp_ray; + for(int i = 0;i < x_size * y_size;i++) + { + int x = i % x_size; + int y = i / x_size; +#endif + + tmp_ray.clear(); + for (int range_y = (int) std::ceil(y - radius); + range_y <= (int) std::floor(y + radius) && range_y <= y_size - 1; range_y++) + { + for (int range_x = (int) std::ceil(x - radius); + range_x <= (int) std::floor(x + radius) && range_x <= x_size - 1; range_x++) + { + //overflow judgement + range_x = range_x > 0 ? range_x : 0; + range_y = range_y > 0 ? range_y : 0; + if (pow(range_x - x, 2.0) + pow(range_y - y, 2.0) <= radius * radius) + { + double dz = std::sqrt(radius * radius - pow(range_x - x, 2.0) - pow(range_y - y, 2.0)); + for (int k = 0; k < input.at(range_x, range_y).size(); k += 2) + { + double down_z = input.at(range_x, range_y)[k] - dz; + double up_z = input.at(range_x, range_y)[k + 1] + dz; + tmp_ray.push_back(MyPoint(down_z, -1)); + tmp_ray.push_back(MyPoint(up_z, 1)); + } + } + } + } + unionMap(tmp_ray, result.at(x, y)); + } +#ifdef USE_TBB + }; + tbb::blocked_range range(0u, (uint32_t)x_size * y_size, GRAIN_SIZE*x_size); + tbb::parallel_for(range, operation); +#endif +} + + +void VoronoiMorphoBruteForce::unionMap(std::vector vec, std::vector &result) +{ + result.clear(); + size_t inter_num = vec.size(); + if(inter_num == 0) + return; // do nothing + vor_assert(inter_num % 2 == 0); + std::sort(vec.begin(), vec.end(), [](auto &left, auto &right) + { + return left.first < right.first; + }); + int flag_ = vec[0].second; + result.push_back(vec[0].first); + for (int k = 1; k < inter_num; ) + { + while (flag_ != 0) + { + flag_ += vec[k].second; + k++; + } + result.push_back(vec[k - 1].first); + if (k < inter_num) + { + result.push_back(vec[k].first); + flag_ = vec[k].second; + k++; + } + } + +} + diff --git a/src/vor3d/VoronoiBruteForce.h b/src/vor3d/VoronoiBruteForce.h new file mode 100644 index 0000000..5d541e4 --- /dev/null +++ b/src/vor3d/VoronoiBruteForce.h @@ -0,0 +1,15 @@ +#pragma once + +#include"vor3d/Voronoi.h" + +namespace voroffset3d +{ + class VoronoiMorphoBruteForce : public VoronoiMorpho + { + typedef std::pair MyPoint; + public: + virtual void dilation(CompressedVolume input, CompressedVolume &result, double radius, double &time_1, double &time_2) override; + private: + void unionMap(std::vector vec, std::vector &result); + }; +} \ No newline at end of file diff --git a/src/vor3d/VoronoiVorPower.cpp b/src/vor3d/VoronoiVorPower.cpp new file mode 100644 index 0000000..ff10c7d --- /dev/null +++ b/src/vor3d/VoronoiVorPower.cpp @@ -0,0 +1,99 @@ +//////////////////////////////////////////////////////////////////////////////// +#include "vor3d/VoronoiVorPower.h" +#include "vor3d/Morpho2D.h" +#include "vor3d/SeparatePower2D.h" +#include "vor3d/Voronoi2D.h" +#include "vor3d/MorphologyOperators.h" +#include "vor3d/HalfDilationOperator.h" +#include "vor3d/Timer.h" +//////////////////////////////////////////////////////////////////////////////// +#ifdef USE_TBB +#include "tbb/blocked_range.h" +#include "tbb/parallel_for.h" +#include "tbb/task_scheduler_init.h" +//////////////////////////////////////////////////////////////////////////////// +#ifndef GRAIN_SIZE +#define GRAIN_SIZE 10 +#endif +#endif + +//////////////////////////////////////////////////////////////////////////////// + +using namespace voroffset3d; + +void VoronoiMorphoVorPower::dilation(CompressedVolume input, CompressedVolume &result, double radius, double &time_1, double &time_2) +{ + int xsize = input.gridSize()(0); + int ysize = input.gridSize()(1); + m_zmin = input.origin()(2) / input.spacing(); + m_zmax = input.origin()(2) / input.spacing() + 2 * input.padding() + input.extent()(2) / input.spacing(); + CompressedVolumeWithRadii output1, output2, mid_output; + CompressedVolume output3, output4; + output1.reshape(xsize, ysize); + output2.reshape(xsize, ysize); + output3.reshape(xsize, ysize); + output4.reshape(xsize, ysize); + mid_output.reshape(xsize, ysize); + result.reset(input.origin(), input.extent(), input.spacing(), input.padding(), xsize, ysize); + + // 1st pass + Timer time_pass_1; +#ifdef USE_TBB + auto firstpass = [&](const tbb::blocked_range &range) + { + VoronoiMorpho2D op_x(ysize, m_zmin, m_zmax, radius, input.spacing()); + for (uint32_t phaseIdx = range.begin(); phaseIdx < range.end(); ++phaseIdx) + { + const uint32_t x = phaseIdx; +#else + // x-direction + VoronoiMorpho2D op_x(ysize, m_zmin, m_zmax, radius, input.spacing()); + for (int x = 0; x < xsize; x++) { +#endif + halfDilate(op_x, true, input, output1, x, 0, 0, +1); + op_x.resetData(); + halfDilate(op_x, true, input, output2, x, ysize - 1, 0, -1); + op_x.resetData(); + } +#ifdef USE_TBB + }; + + tbb::blocked_range rangex(0u, (uint32_t)xsize, GRAIN_SIZE); + tbb::parallel_for(rangex, firstpass); +#endif + + unionMap(output1, output2, mid_output); + time_1 = time_pass_1.get(); + + // 2nd pass + Timer time_pass_2; +#ifdef USE_TBB + auto secondpass = [&](const tbb::blocked_range &range) + { + SeparatePowerMorpho2D op_y(xsize, m_zmin, m_zmax, input.spacing()); + for (uint32_t phaseIdy = range.begin(); phaseIdy < range.end(); ++phaseIdy) + { + const uint32_t y = phaseIdy; +#else + SeparatePowerMorpho2D op_y(xsize, m_zmin, m_zmax, input.spacing()); + for(int y=0;y rangey(0u, (uint32_t)ysize, GRAIN_SIZE); + tbb::parallel_for(rangey, secondpass); +#endif + + unionMap(output3, output4, result); + time_2 = time_pass_2.get(); +} + + + diff --git a/src/vor3d/VoronoiVorPower.h b/src/vor3d/VoronoiVorPower.h new file mode 100644 index 0000000..14112ed --- /dev/null +++ b/src/vor3d/VoronoiVorPower.h @@ -0,0 +1,15 @@ +#pragma once + +#include"vor3d/Voronoi.h" + +namespace voroffset3d +{ + class VoronoiMorphoVorPower : public VoronoiMorpho + { + public: + virtual void dilation(CompressedVolume input, CompressedVolume &result, double radius, double &time_1, double &time_2) override; + + private: + double m_zmin, m_zmax; + }; +} \ No newline at end of file