diff --git a/CMakeLists.txt b/CMakeLists.txt index 59258b56..88e3fc32 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -131,6 +131,7 @@ set(CMAKE_CXX_FLAGS_RELEASE "-O2") #== GTest option(RUN_GTEST "Downloads google unit test API and runs google test scripts to test Nyxus" OFF) +option(RUN_BENCHMARKS "Build Nyxus benchmark executables" OFF) #==== Source files set(SOURCE @@ -432,6 +433,15 @@ if(BUILD_CLI) target_link_libraries(nyxus PRIVATE ${Nyxus_LIBRARIES}) endif() +if(RUN_BENCHMARKS) + add_executable(bench_convex_hull + benchmarks/bench_convex_hull.cpp + ${SOURCE} + ) + target_include_directories(bench_convex_hull PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}) + target_link_libraries(bench_convex_hull PRIVATE ${Nyxus_LIBRARIES}) +endif() + message(STATUS "Nyxus Dependencies Summary") message(STATUS "================================") message(STATUS "Support | Requested | Found") diff --git a/benchmarks/README.md b/benchmarks/README.md new file mode 100644 index 00000000..ecbade38 --- /dev/null +++ b/benchmarks/README.md @@ -0,0 +1,36 @@ +# Nyxus Benchmarks + +This directory contains small, dependency-light microbenchmarks for performance-sensitive feature code. + +## Convex Hull + +`bench_convex_hull` measures the 2D convex hull path used by `ConvexHullFeature::calculate()`. +It intentionally reports both timing and output shape (`hull_size`, `hull_area`, `solidity`) so performance comparisons can also catch ordering-sensitive correctness regressions. + +Build: + +```bash +cmake -S . -B build-bench -DRUN_BENCHMARKS=ON -DNOEXTRAS=ON +cmake --build build-bench --target bench_convex_hull -j +``` + +Run the default workload: + +```bash +./build-bench/bench_convex_hull +``` + +Run a focused workload: + +```bash +./build-bench/bench_convex_hull --cycles 5 --case y_major_raster --side 1024 +``` + +Input-order cases: + +- `x_major_sorted`: already sorted by the convex hull comparator (`x`, then `y`). +- `y_major_raster`: row-major scan order, matching the usual trivial ROI pixel ingestion path. +- `reversed`: reverse of the sorted comparator order. +- `shuffled`: deterministic pseudorandom order. + +For before/after comparisons, build the same target in two clean worktrees and compare median times for the same case/side/cycle settings. diff --git a/benchmarks/bench_convex_hull.cpp b/benchmarks/bench_convex_hull.cpp new file mode 100644 index 00000000..508e320e --- /dev/null +++ b/benchmarks/bench_convex_hull.cpp @@ -0,0 +1,285 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "../src/nyx/features/convex_hull.h" + +namespace +{ +using Clock = std::chrono::steady_clock; + +struct Workload +{ + std::string case_name; + int side = 0; +}; + +struct Options +{ + int cycles = 3; + std::vector workloads; +}; + +std::vector make_x_major_sorted(int side) +{ + std::vector pixels; + pixels.reserve(static_cast(side) * static_cast(side)); + + for (int x = 0; x < side; ++x) + for (int y = 0; y < side; ++y) + pixels.emplace_back(x, y, 1u); + + return pixels; +} + +std::vector make_y_major_raster(int side) +{ + std::vector pixels; + pixels.reserve(static_cast(side) * static_cast(side)); + + for (int y = 0; y < side; ++y) + for (int x = 0; x < side; ++x) + pixels.emplace_back(x, y, 1u); + + return pixels; +} + +std::vector make_pixels(const std::string& case_name, int side) +{ + if (case_name == "x_major_sorted") + return make_x_major_sorted(side); + + if (case_name == "y_major_raster") + return make_y_major_raster(side); + + auto pixels = make_x_major_sorted(side); + if (case_name == "reversed") + { + std::reverse(pixels.begin(), pixels.end()); + return pixels; + } + + if (case_name == "shuffled") + { + std::mt19937 rng(1234567u + static_cast(side)); + std::shuffle(pixels.begin(), pixels.end(), rng); + return pixels; + } + + throw std::runtime_error("unknown benchmark case: " + case_name); +} + +int reps_for_side(int side) +{ + if (side <= 64) + return 160; + if (side <= 128) + return 80; + if (side <= 256) + return 30; + if (side <= 512) + return 10; + return 4; +} + +double median(std::vector values) +{ + std::sort(values.begin(), values.end()); + const size_t n = values.size(); + if (n % 2 == 1) + return values[n / 2]; + return 0.5 * (values[n / 2 - 1] + values[n / 2]); +} + +void initialize_convex_roi(LR& roi, const std::vector& pixels, int side) +{ + roi.raw_pixels = pixels; + roi.convHull_CH.clear(); + roi.fvals.assign(static_cast(Nyxus::FeatureIMQ::_COUNT_), std::vector(1, 0.0)); + + // ConvexHullFeature depends on PERIMETER only for circularity; this keeps + // the benchmark focused on convex hull construction, not contour tracing. + roi.fvals[static_cast(Nyxus::Feature2D::PERIMETER)][0] = 4.0 * static_cast(side - 1); +} + +void calculate_convex_hull(LR& roi, ConvexHullFeature& feature, const Fsettings& settings) +{ + feature.calculate(roi, settings); + feature.save_value(roi.fvals); +} + +void run_workload(int cycle, const Workload& workload) +{ + const auto pixels = make_pixels(workload.case_name, workload.side); + const int reps = reps_for_side(workload.side); + + LR roi; + initialize_convex_roi(roi, pixels, workload.side); + + Fsettings settings; + ConvexHullFeature feature; + + uint64_t checksum = 0; + for (int i = 0; i < 2; ++i) + { + calculate_convex_hull(roi, feature, settings); + checksum += static_cast(roi.convHull_CH.size()); + checksum += static_cast(roi.fvals[static_cast(Nyxus::Feature2D::CONVEX_HULL_AREA)][0]); + } + + std::vector elapsed_ms; + elapsed_ms.reserve(static_cast(reps)); + + for (int i = 0; i < reps; ++i) + { + const auto start = Clock::now(); + calculate_convex_hull(roi, feature, settings); + const auto stop = Clock::now(); + + checksum += static_cast(roi.convHull_CH.size()); + checksum += static_cast(roi.fvals[static_cast(Nyxus::Feature2D::CONVEX_HULL_AREA)][0]); + elapsed_ms.push_back(std::chrono::duration(stop - start).count()); + } + + const double sum = std::accumulate(elapsed_ms.begin(), elapsed_ms.end(), 0.0); + const auto [min_it, max_it] = std::minmax_element(elapsed_ms.begin(), elapsed_ms.end()); + const double hull_area = roi.fvals[static_cast(Nyxus::Feature2D::CONVEX_HULL_AREA)][0]; + const double solidity = roi.fvals[static_cast(Nyxus::Feature2D::SOLIDITY)][0]; + + std::cout + << cycle << ',' + << workload.case_name << ',' + << workload.side << ',' + << pixels.size() << ',' + << reps << ',' + << std::fixed << std::setprecision(6) + << median(elapsed_ms) << ',' + << (sum / static_cast(elapsed_ms.size())) << ',' + << *min_it << ',' + << *max_it << ',' + << roi.convHull_CH.size() << ',' + << hull_area << ',' + << solidity << ',' + << checksum + << '\n'; +} + +void add_default_workloads(std::vector& workloads) +{ + for (int side : {128, 256, 512, 1024}) + workloads.push_back({"x_major_sorted", side}); + + for (int side : {64, 128, 256, 512, 1024}) + workloads.push_back({"y_major_raster", side}); + + // Keep pathological unsorted cases capped so this can also be run on older + // revisions where shuffled input can create a huge invalid hull. + for (int side : {64, 128, 256}) + workloads.push_back({"reversed", side}); + + for (int side : {64, 128}) + workloads.push_back({"shuffled", side}); +} + +void print_usage(const char* program) +{ + std::cerr + << "Usage: " << program << " [--cycles N] [--case NAME --side N ...]\n" + << "\n" + << "Cases: x_major_sorted, y_major_raster, reversed, shuffled\n" + << "If no --case/--side pairs are provided, a default workload is used.\n"; +} + +int parse_int(std::string_view value, std::string_view flag) +{ + const std::string text(value); + char* end = nullptr; + const long parsed = std::strtol(text.c_str(), &end, 10); + if (end == nullptr || *end != '\0' || parsed <= 0) + throw std::runtime_error("invalid positive integer for " + std::string(flag)); + return static_cast(parsed); +} + +Options parse_options(int argc, char** argv) +{ + Options options; + std::string pending_case; + + for (int i = 1; i < argc; ++i) + { + const std::string_view arg(argv[i]); + if (arg == "--help" || arg == "-h") + { + print_usage(argv[0]); + std::exit(0); + } + + if (arg == "--cycles") + { + if (++i >= argc) + throw std::runtime_error("--cycles requires a value"); + options.cycles = parse_int(argv[i], "--cycles"); + continue; + } + + if (arg == "--case") + { + if (++i >= argc) + throw std::runtime_error("--case requires a value"); + pending_case = argv[i]; + continue; + } + + if (arg == "--side") + { + if (pending_case.empty()) + throw std::runtime_error("--side must follow --case"); + if (++i >= argc) + throw std::runtime_error("--side requires a value"); + options.workloads.push_back({pending_case, parse_int(argv[i], "--side")}); + pending_case.clear(); + continue; + } + + throw std::runtime_error("unknown argument: " + std::string(arg)); + } + + if (!pending_case.empty()) + throw std::runtime_error("--case must be paired with --side"); + + if (options.workloads.empty()) + add_default_workloads(options.workloads); + + return options; +} +} + +int main(int argc, char** argv) +{ + try + { + const Options options = parse_options(argc, argv); + std::cout << "cycle,case,side,points,reps,median_ms,mean_ms,min_ms,max_ms,hull_size,hull_area,solidity,checksum\n"; + + for (int cycle = 1; cycle <= options.cycles; ++cycle) + for (const Workload& workload : options.workloads) + run_workload(cycle, workload); + } + catch (const std::exception& e) + { + std::cerr << "bench_convex_hull: " << e.what() << "\n"; + print_usage(argv[0]); + return 1; + } + + return 0; +} diff --git a/src/nyx/features/convex_hull_nontriv.cpp b/src/nyx/features/convex_hull_nontriv.cpp index c6ab0356..483d367a 100644 --- a/src/nyx/features/convex_hull_nontriv.cpp +++ b/src/nyx/features/convex_hull_nontriv.cpp @@ -1,54 +1,54 @@ #define _USE_MATH_DEFINES // For M_PI, etc. #include #include -#include "../dataset.h" -#include "../feature_method.h" -#include "../feature_settings.h" -#include "image_matrix_nontriv.h" -#include "convex_hull.h" - -using namespace Nyxus; - -ConvexHullFeature::ConvexHullFeature() : FeatureMethod("ConvexHullFeature") -{ - provide_features ({ Feature2D::CONVEX_HULL_AREA, Feature2D::SOLIDITY, Feature2D::CIRCULARITY }); - add_dependencies ({ Feature2D::PERIMETER }); -} - -void ConvexHullFeature::save_value(std::vector>& fvals) -{ - fvals [(int)Feature2D::CONVEX_HULL_AREA][0] = area; - fvals [(int)Feature2D::SOLIDITY][0] = solidity; - fvals [(int)Feature2D::CIRCULARITY][0] = circularity; -} - -void ConvexHullFeature::cleanup_instance() -{ - area = solidity = circularity = 0.0; -} - -void ConvexHullFeature::calculate (LR& r, const Fsettings& s) -{ - // Build the convex hull - build_convex_hull(r.raw_pixels, r.convHull_CH); - - // Calculate related features - double s_hull = polygon_area(r.convHull_CH), - s_roi = r.raw_pixels.size(), - p = r.fvals[(int)Feature2D::PERIMETER][0]; - area = s_hull; - solidity = s_roi / s_hull; - circularity = sqrt(4.0 * M_PI * s_roi / (p*p)); -} - +#include "../dataset.h" +#include "../feature_method.h" +#include "../feature_settings.h" +#include "image_matrix_nontriv.h" +#include "convex_hull.h" + +using namespace Nyxus; + +ConvexHullFeature::ConvexHullFeature() : FeatureMethod("ConvexHullFeature") +{ + provide_features ({ Feature2D::CONVEX_HULL_AREA, Feature2D::SOLIDITY, Feature2D::CIRCULARITY }); + add_dependencies ({ Feature2D::PERIMETER }); +} + +void ConvexHullFeature::save_value(std::vector>& fvals) +{ + fvals [(int)Feature2D::CONVEX_HULL_AREA][0] = area; + fvals [(int)Feature2D::SOLIDITY][0] = solidity; + fvals [(int)Feature2D::CIRCULARITY][0] = circularity; +} + +void ConvexHullFeature::cleanup_instance() +{ + area = solidity = circularity = 0.0; +} + +void ConvexHullFeature::calculate (LR& r, const Fsettings& s) +{ + // Build the convex hull + build_convex_hull(r.raw_pixels, r.convHull_CH); + + // Calculate related features + double s_hull = polygon_area(r.convHull_CH), + s_roi = r.raw_pixels.size(), + p = r.fvals[(int)Feature2D::PERIMETER][0]; + area = s_hull; + solidity = s_roi / s_hull; + circularity = sqrt(4.0 * M_PI * s_roi / (p*p)); +} + void ConvexHullFeature::build_convex_hull (const std::vector& cloud, std::vector& convhull) { convhull.clear(); - - // Skip calculation if the ROI is too small - if (cloud.size() < 2) - return; - + + // Skip calculation if the ROI is too small + if (cloud.size() < 2) + return; + std::vector& upperCH = convhull; std::vector lowerCH; @@ -86,23 +86,23 @@ void ConvexHullFeature::build_convex_hull (const std::vector& cloud, std lowerCH.pop_back(); lowerCH.push_back(ordered[n - i - 1]); } - - // We could use - // upperCH.insert (upperCH.end(), lowerCH.begin(), lowerCH.end()); - // but we don't need duplicate points in the result contour - for (auto& p : lowerCH) - if (std::find(upperCH.begin(), upperCH.end(), p) == upperCH.end()) - upperCH.push_back(p); -} - -void ConvexHullFeature::build_convex_hull (const OutOfRamPixelCloud& cloud, std::vector& convhull) -{ - convhull.clear(); - - // Skip calculation if the ROI is too small - if (cloud.size() < 2) - return; - + + // We could use + // upperCH.insert (upperCH.end(), lowerCH.begin(), lowerCH.end()); + // but we don't need duplicate points in the result contour + for (auto& p : lowerCH) + if (std::find(upperCH.begin(), upperCH.end(), p) == upperCH.end()) + upperCH.push_back(p); +} + +void ConvexHullFeature::build_convex_hull (const OutOfRamPixelCloud& cloud, std::vector& convhull) +{ + convhull.clear(); + + // Skip calculation if the ROI is too small + if (cloud.size() < 2) + return; + std::vector& upperCH = convhull; std::vector lowerCH; @@ -144,87 +144,87 @@ void ConvexHullFeature::build_convex_hull (const OutOfRamPixelCloud& cloud, std: lowerCH.pop_back(); lowerCH.push_back(orderedCloud[n - i - 1]); } - - // We could use - // upperCH.insert (upperCH.end(), lowerCH.begin(), lowerCH.end()); - // but we don't need duplicate points in the result contour - for (auto& p : lowerCH) - if (std::find(upperCH.begin(), upperCH.end(), p) == upperCH.end()) - upperCH.push_back(p); -} - -// Sort criterion: points are sorted with respect to their x-coordinate. -// If two points have the same x-coordinate then we compare their y-coordinates -bool ConvexHullFeature::compare_locations (const Pixel2& lhs, const Pixel2& rhs) -{ - return (lhs.x < rhs.x) || (lhs.x == rhs.x && lhs.y < rhs.y); -} - -// Check if three points make a right turn using cross product -bool ConvexHullFeature::right_turn (const Pixel2& P1, const Pixel2& P2, const Pixel2& P3) -{ - return ((P3.x - P1.x) * (P2.y - P1.y) - (P3.y - P1.y) * (P2.x - P1.x)) > 0; -} - -double ConvexHullFeature::polygon_area (const std::vector& vertices) -{ - // Blank polygon? - if (vertices.size() == 0) - return 0.0; - - // Normal polygon - double area = 0.0; - size_t n = vertices.size(); - for (size_t i = 0; i < n-1; i++) - { - const Pixel2 &p1 = vertices[i], - &p2 = vertices[i+1]; - area += p1.x * p2.y - p1.y * p2.x; - } - const Pixel2& p1 = vertices[0], - & p2 = vertices[n-1]; - area += p1.x * p2.y - p1.y * p2.x; - area = std::abs(area) / 2.0; - return area; -} - -void ConvexHullFeature::osized_calculate (LR& r, const Fsettings& s, ImageLoader& imloader) -{ - // Build the convex hull - build_convex_hull (r.raw_pixels_NT, r.convHull_CH); - - // Calculate related features - double s_hull = polygon_area(r.convHull_CH), - s_roi = r.raw_pixels_NT.size(), - p = r.fvals[(int)Feature2D::PERIMETER][0]; - area = s_hull; - solidity = s_roi / s_hull; - circularity = sqrt(4.0 * M_PI * s_roi / (p * p)); -} - -void ConvexHullFeature::extract (LR& r, const Fsettings& s) -{ - ConvexHullFeature f; - f.calculate (r, s); - f.save_value (r.fvals); -} - -namespace Nyxus -{ - void parallelReduceConvHull (size_t start, size_t end, std::vector* ptrLabels, std::unordered_map * ptrLabelData, const Fsettings & s, const Dataset & _) - { - for (auto i = start; i < end; i++) - { - int lab = (*ptrLabels)[i]; - LR& r = (*ptrLabelData)[lab]; - - if (r.has_bad_data()) - continue; - - ConvexHullFeature f; - f.calculate (r, s); - f.save_value (r.fvals); - } - } -} - + + // We could use + // upperCH.insert (upperCH.end(), lowerCH.begin(), lowerCH.end()); + // but we don't need duplicate points in the result contour + for (auto& p : lowerCH) + if (std::find(upperCH.begin(), upperCH.end(), p) == upperCH.end()) + upperCH.push_back(p); +} + +// Sort criterion: points are sorted with respect to their x-coordinate. +// If two points have the same x-coordinate then we compare their y-coordinates +bool ConvexHullFeature::compare_locations (const Pixel2& lhs, const Pixel2& rhs) +{ + return (lhs.x < rhs.x) || (lhs.x == rhs.x && lhs.y < rhs.y); +} + +// Check if three points make a right turn using cross product +bool ConvexHullFeature::right_turn (const Pixel2& P1, const Pixel2& P2, const Pixel2& P3) +{ + return ((P3.x - P1.x) * (P2.y - P1.y) - (P3.y - P1.y) * (P2.x - P1.x)) > 0; +} + +double ConvexHullFeature::polygon_area (const std::vector& vertices) +{ + // Blank polygon? + if (vertices.size() == 0) + return 0.0; + + // Normal polygon + double area = 0.0; + size_t n = vertices.size(); + for (size_t i = 0; i < n-1; i++) + { + const Pixel2 &p1 = vertices[i], + &p2 = vertices[i+1]; + area += p1.x * p2.y - p1.y * p2.x; + } + const Pixel2& p1 = vertices[n-1], + & p2 = vertices[0]; + area += p1.x * p2.y - p1.y * p2.x; + area = std::abs(area) / 2.0; + return area; +} + +void ConvexHullFeature::osized_calculate (LR& r, const Fsettings& s, ImageLoader& imloader) +{ + // Build the convex hull + build_convex_hull (r.raw_pixels_NT, r.convHull_CH); + + // Calculate related features + double s_hull = polygon_area(r.convHull_CH), + s_roi = r.raw_pixels_NT.size(), + p = r.fvals[(int)Feature2D::PERIMETER][0]; + area = s_hull; + solidity = s_roi / s_hull; + circularity = sqrt(4.0 * M_PI * s_roi / (p * p)); +} + +void ConvexHullFeature::extract (LR& r, const Fsettings& s) +{ + ConvexHullFeature f; + f.calculate (r, s); + f.save_value (r.fvals); +} + +namespace Nyxus +{ + void parallelReduceConvHull (size_t start, size_t end, std::vector* ptrLabels, std::unordered_map * ptrLabelData, const Fsettings & s, const Dataset & _) + { + for (auto i = start; i < end; i++) + { + int lab = (*ptrLabels)[i]; + LR& r = (*ptrLabelData)[lab]; + + if (r.has_bad_data()) + continue; + + ConvexHullFeature f; + f.calculate (r, s); + f.save_value (r.fvals); + } + } +} + diff --git a/tests/test_gldm.h b/tests/test_gldm.h index 0d2d1dd2..66b00c13 100644 --- a/tests/test_gldm.h +++ b/tests/test_gldm.h @@ -75,7 +75,7 @@ void test_gldm_lde() void test_gldm_lgle() { - test_gldm_feature(Nyxus::Feature2D::GLDM_SDE, "GLDM_SDE"); + test_gldm_feature(Nyxus::Feature2D::GLDM_LGLE, "GLDM_LGLE"); } void test_gldm_hgle() @@ -131,4 +131,4 @@ void test_gldm_dv() void test_gldm_de() { test_gldm_feature(Nyxus::Feature2D::GLDM_DE, "GLDM_DE"); -} \ No newline at end of file +} diff --git a/tests/test_ibsi_gldm.h b/tests/test_ibsi_gldm.h index 901ea779..0f332ca5 100644 --- a/tests/test_ibsi_gldm.h +++ b/tests/test_ibsi_gldm.h @@ -131,7 +131,7 @@ void test_ibsi_gldm_lde() void test_ibsi_gldm_lgle() { - test_ibsi_gldm_feature(Nyxus::Feature2D::GLDM_SDE, "GLDM_SDE"); + test_ibsi_gldm_feature(Nyxus::Feature2D::GLDM_LGLE, "GLDM_LGLE"); } void test_ibsi_gldm_hgle() @@ -187,4 +187,4 @@ void test_ibsi_gldm_dv() void test_ibsi_gldm_de() { test_ibsi_gldm_feature(Nyxus::Feature2D::GLDM_DE, "GLDM_DE"); -} \ No newline at end of file +} diff --git a/tests/test_shape_morphology_2d.h b/tests/test_shape_morphology_2d.h index 6e294339..561fd37c 100644 --- a/tests/test_shape_morphology_2d.h +++ b/tests/test_shape_morphology_2d.h @@ -49,8 +49,8 @@ static std::unordered_map shape2d_truth{ {"EDGE_MIN_INTENSITY", 12.0}, {"EDGE_INTEGRATED_INTENSITY", 753.0}, {"CIRCULARITY", 0.671081973229055}, - {"CONVEX_HULL_AREA", 24.0}, - {"SOLIDITY", 1.08333333333333}, + {"CONVEX_HULL_AREA", 20.0}, + {"SOLIDITY", 1.3}, {"EULER_NUMBER", -2.0}, {"FRACT_DIM_BOXCOUNT", -0.830074998557687}, {"FRACT_DIM_PERIMETER", -1.97227924244155},