diff --git a/src/pymodule/amplitudecorrection/functions/module.hpp b/src/pymodule/amplitudecorrection/functions/module.hpp index d013a84..cec70af 100644 --- a/src/pymodule/amplitudecorrection/functions/module.hpp +++ b/src/pymodule/amplitudecorrection/functions/module.hpp @@ -134,27 +134,6 @@ void init_functions(pybind11::module& m) py::arg("min_beam_index") = std::nullopt, py::arg("max_beam_index") = std::nullopt, py::arg("mp_cores") = 1); - m.def( - "inplace_beam_sample_correction2", - // &inplace_beam_sample_correction2, xt::pytensor> , - // use lambda here because for some reason xt::pytensor is seen as reference from python on - // xt::pytensor& is not seen as reference - [](xt::pytensor wci, - const xt::pytensor& per_beam_offset, - const xt::pytensor& per_sample_offset, - std::optional min_beam_index = std::nullopt, - std::optional max_beam_index = std::nullopt, - int mp_cores = 1) { - inplace_beam_sample_correction2( - wci, per_beam_offset, per_sample_offset, min_beam_index, max_beam_index, mp_cores); - }, - DOC_functions(inplace_beam_sample_correction2), - py::arg("wci"), - py::arg("per_beam_offset"), - py::arg("per_sample_offset"), - py::arg("min_beam_index") = std::nullopt, - py::arg("max_beam_index") = std::nullopt, - py::arg("mp_cores") = 1); m.def( "inplace_beam_correction", diff --git a/src/tests/algorithms/functions/absorption.test.cpp b/src/tests/algorithms/functions/absorption.test.cpp new file mode 100644 index 0000000..3abfcbc --- /dev/null +++ b/src/tests/algorithms/functions/absorption.test.cpp @@ -0,0 +1,61 @@ +// SPDX-FileCopyrightText: 2022 - 2023 Peter Urban, Ghent University +// +// SPDX-License-Identifier: MPL-2.0 + +#include +#include +#include + +#include "../../../themachinethatgoesping/algorithms/amplitudecorrection/functions/absorption.hpp" + +// using namespace testing; +using namespace std; +using namespace themachinethatgoesping::algorithms::amplitudecorrection::functions; + +#define TESTTAG "[location]" + +/** + * @brief These tests just make sure that the return values are consistent with previous versions + * of the functions. The values are not checked for correctness. + * In case the functions are changed and updated to produce 'better' results, the tests should be + * updated accordingly. However, these changes should be documented and justified und users should + * be informed about these changes as they will affect quantification results. + * + */ +TEST_CASE("Absorption functions should reproduce previously computed results", TESTTAG) +{ + using Catch::Approx; + + SECTION("its90_to_its68") + { + CHECK(its90_to_its68(0.0) == Approx(0.0)); + CHECK(its90_to_its68(0.0f) == Approx(0.0f)); + CHECK(its90_to_its68(1.0) == Approx(1.00024)); + CHECK(its90_to_its68(1.0f) == Approx(1.00024f)); + CHECK(its90_to_its68(10.0) == Approx(10.0024)); + CHECK(its90_to_its68(10.0f) == Approx(10.0024f)); + CHECK(its90_to_its68(-100.0) == Approx(-100.024)); + CHECK(its90_to_its68(-100.0f) == Approx(-100.024f)); + } + + SECTION("calc_sound_velocity") + { + CHECK(calc_sound_velocity(0.0, 0.0, 0.0) == Approx(1402.40099361970374048)); + CHECK(calc_sound_velocity(0.0, 0.0, 35.0) == Approx(1449.02536060550323782)); + CHECK(calc_sound_velocity(0.0, 10.0, 35.0) == Approx(1489.78894590650997998)); + CHECK(calc_sound_velocity(100.0, 10.0, 35.0) == Approx(1491.42729295043295679)); + CHECK(calc_sound_velocity(100.0, 10.0, 35.0, 45.0, 0.0) == Approx(1491.43205948529271154)); + CHECK(calc_sound_velocity(100.0, 10.0, 35.0, 45.0, 40.0) == Approx(1491.43029656977273589)); + } + + SECTION("calc_absorption_coefficient_db_m") + { + CHECK(calc_absorption_coefficient_db_m(1000.0, 0.0, 1500.0, 10.0, 35.0) == Approx(0.0000595565729641)); + CHECK(calc_absorption_coefficient_db_m(10000.0, 100.0, 1500.0, 3.0, 35.0) == Approx(0.00114356087217127)); + CHECK(calc_absorption_coefficient_db_m(100000.0, 0.0, 1400.0, 10.0, 35.0, 8.0) == Approx(0.03550226534746687)); + CHECK(calc_absorption_coefficient_db_m(1000000.0, 100.0, 1500.0, 4.0, 35.0, 8.0) == Approx(0.43118748483875202)); + CHECK(calc_absorption_coefficient_db_m(1000.0, 0.0, 1500.0, 10.0, 35.0, 12.0) == Approx(0.06725602852617596)); + CHECK(calc_absorption_coefficient_db_m(1000.0, 100.0, 1500.0, 7.0, 35.0, 12.0) == Approx(0.06766628550920389)); + CHECK(calc_absorption_coefficient_db_m(10000.0, 0.0, 1500.0, 10.0, 35.0, 12.0) == Approx(0.1500729200888182)); + } +} diff --git a/src/tests/algorithms/functions/rangecorrection.test.cpp b/src/tests/algorithms/functions/rangecorrection.test.cpp new file mode 100644 index 0000000..8d96919 --- /dev/null +++ b/src/tests/algorithms/functions/rangecorrection.test.cpp @@ -0,0 +1,112 @@ +// SPDX-FileCopyrightText: 2022 - 2023 Peter Urban, Ghent University +// +// SPDX-License-Identifier: MPL-2.0 + +#include +#include +#include + +#include +#include + +#include +#include + +// using namespace testing; +using namespace std; +using namespace themachinethatgoesping::algorithms::amplitudecorrection::functions; + +#define TESTTAG "[location]" + +/** + * @brief These tests just make sure that the return values are consistent with previous versions + * of the functions. The values are not checked for correctness. + * In case the functions are changed and updated to produce 'better' results, the tests should be + * updated accordingly. However, these changes should be documented and justified und users should + * be informed about these changes as they will affect quantification results. + * + */ +TEST_CASE("Rangecorrection functions should reproduce previously computed results", TESTTAG) +{ + using Catch::Approx; + + SECTION("get_sample_numbers_plus_half") + { + auto sample_numbers = get_sample_numbers_plus_half>(0, 10, 1); + REQUIRE(sample_numbers.size() == 11); + CHECK(sample_numbers(0) == Approx(0.5)); + CHECK(sample_numbers(1) == Approx(1.5)); + CHECK(sample_numbers(2) == Approx(2.5)); + CHECK(sample_numbers(3) == Approx(3.5)); + CHECK(sample_numbers(4) == Approx(4.5)); + CHECK(sample_numbers(5) == Approx(5.5)); + CHECK(sample_numbers(6) == Approx(6.5)); + CHECK(sample_numbers(7) == Approx(7.5)); + CHECK(sample_numbers(8) == Approx(8.5)); + CHECK(sample_numbers(9) == Approx(9.5)); + CHECK(sample_numbers(10) == Approx(10.5)); + } + + SECTION("approximate_range_factor") + { + CHECK(approximate_range_factor(0.0001, 1500.0) == Approx(0.075)); + CHECK(approximate_range_factor(0.0001f, 1500.0f) == Approx(0.075f)); + CHECK(approximate_range_factor(0.0001, 1480.0) == Approx(0.074)); + CHECK(approximate_range_factor(0.0001f, 1480.0f) == Approx(0.074f)); + CHECK(approximate_range_factor(0.0001, 1600.0) == Approx(0.08)); + CHECK(approximate_range_factor(0.0001f, 1600.0f) == Approx(0.08f)); + CHECK(approximate_range_factor(0.0001, 1400.0) == Approx(0.07)); + CHECK(approximate_range_factor(0.0001f, 1400.0f) == Approx(0.07f)); + } + + SECTION("approximate_ranges_1") + { + auto ranges = approximate_ranges>(0.0001, 1450.0, 0, 10, 1); + REQUIRE(ranges.size() == 11); + CHECK(ranges(0) == Approx(0.03625)); + CHECK(ranges(1) == Approx(0.10875000000000001)); + CHECK(ranges(2) == Approx(0.18125000000000002)); + CHECK(ranges(3) == Approx(0.25375000000000003)); + CHECK(ranges(4) == Approx(0.32625000000000004)); + CHECK(ranges(5) == Approx(0.39875000000000005)); + CHECK(ranges(6) == Approx(0.47125000000000006)); + CHECK(ranges(7) == Approx(0.54375000000000007)); + CHECK(ranges(8) == Approx(0.61625000000000008)); + CHECK(ranges(9) == Approx(0.68875000000000008)); + CHECK(ranges(10) == Approx(0.76125000000000009)); + } + + SECTION("approximate_ranges_2") + { + auto ranges2 = approximate_ranges>(0.001, 1450.0, 100, 300, 4); + auto ranges1 = get_sample_numbers_plus_half>(100, 300, 4) * + approximate_range_factor(0.001, 1450.0); + + REQUIRE(ranges1.size() == ranges2.size()); + for (size_t i = 0; i < ranges1.size(); ++i) + CHECK(ranges2(i) == Approx(ranges1(i))); + } + + SECTION("compute_cw_range_correction") + { + auto ranges = xt::eval(xt::linspace(0.5, 10.5, 10)); + + for (const std::optional& absorption : + std::vector>{ std::nullopt, 0.f, 0.0124f }) + { + for (const std::optional& tvg : + std::vector>{ std::nullopt, -13.f, 0.f, 14.f }) + { + auto correction = compute_cw_range_correction(ranges, absorption, tvg); + REQUIRE(correction.size() == 10); + + for (size_t i = 0; i < 10; ++i) + { + INFO(fmt::format("i = {}, ranges(i) = {}", i, ranges(i))); + REQUIRE(correction(i) == + Approx(ranges(i) * absorption.value_or(0) * 2 + tvg.value_or(0) * std::log10(ranges(i)))); + } + } + } + } +} diff --git a/src/tests/algorithms/functions/wcicorrection.test.cpp b/src/tests/algorithms/functions/wcicorrection.test.cpp new file mode 100644 index 0000000..bd0af04 --- /dev/null +++ b/src/tests/algorithms/functions/wcicorrection.test.cpp @@ -0,0 +1,346 @@ +// SPDX-FileCopyrightText: 2022 - 2023 Peter Urban, Ghent University +// +// SPDX-License-Identifier: MPL-2.0 + +#include +#include +#include + +#include +#include + +#include +#include +#include + +// using namespace testing; +using namespace std; +using namespace themachinethatgoesping::algorithms::amplitudecorrection::functions; + +#define TESTTAG "[location]" + +template +void check_wci(std::string name, + const t_xtensor_2d& result, + const t_xtensor_2d& wci, + const t_xtensor_1d& per_beam_offset, + const t_xtensor_1d& per_sample_offset, + std::optional min_beam_index = std::nullopt, + std::optional max_beam_index = std::nullopt) +{ + using Catch::Approx; + + size_t nsamples = 20; + size_t nbeams = 10; + + // max_beam_nr and min_beam_nr are allowed to be > wci.shape(0) and will be limited to the + // actual number of beams + size_t max_beam_nr = max_beam_index.value_or(wci.shape(0) - 1); + if (max_beam_nr >= wci.shape(0)) + max_beam_nr = wci.shape(0) - 1; + + REQUIRE(result.dimension() == 2); + REQUIRE(wci.dimension() == 2); + REQUIRE(per_beam_offset.dimension() == 1); + REQUIRE(per_sample_offset.dimension() == 1); + + REQUIRE(result.shape(0) == nbeams); + REQUIRE(result.shape(1) == nsamples); + REQUIRE(wci.shape(0) == nbeams); + REQUIRE(wci.shape(1) == nsamples); + REQUIRE(per_beam_offset.shape(0) == nbeams); + REQUIRE(per_sample_offset.shape(0) == nsamples); + + for (size_t bn = min_beam_index.value_or(0); bn <= max_beam_nr; ++bn) + { + for (size_t sn = 0; sn < nsamples; ++sn) + { + // std::string msg = fmt::format( + // "ERROR[{}]: per_beam_offset({}) = {}, per_sample_offset({}) = {}, wci({},{}) = {}", + // name, + // bn, + // per_beam_offset(bn), + // sn, + // per_sample_offset(sn), + // bn, + // sn, + // wci(bn, sn)); + + // if (min_beam_index.has_value()) + // msg += fmt::format(", min_beam_index = {}", min_beam_index.value()); + + // if (max_beam_index.has_value()) + // msg += fmt::format(", max_beam_index = {}", max_beam_index.value()); + + // INFO(msg); + + REQUIRE(result(bn, sn) == + Approx(wci(bn, sn) + per_beam_offset(bn) + per_sample_offset(sn))); + } + } + + // check that the rest of the array is not touched + if (min_beam_index.has_value()) + { + size_t min_beam_nr = min_beam_index.value_or(0); + if (min_beam_nr >= wci.shape(0)) + min_beam_nr = wci.shape(0); + + for (size_t bn = 0; bn < min_beam_nr; ++bn) + { + for (size_t sn = 0; sn < nsamples; ++sn) + { + // INFO( + // fmt::format("ERROR[{}]: Value should be untouchted (1)! bn = {}/min {}, sn = {}" + // "wci({},{}) = {}", + // name, + // bn, + // min_beam_index.value(), + // sn, + // bn, + // sn, + // wci(bn, sn))); + REQUIRE(result(bn, sn) == 1); + } + } + } + + if (max_beam_index.has_value()) + for (size_t bn = max_beam_index.value() + 1; bn < nbeams; ++bn) + { + for (size_t sn = 0; sn < nsamples; ++sn) + { + // INFO(fmt::format( + // "ERROR[{}]: Value should be untouchted (1)! bn = {}/max {}, sn = {}, " + // "wci({},{}) = {}", + // name, + // bn, + // max_beam_index.value(), + // sn, + // bn, + // sn, + // wci(bn, sn))); + REQUIRE(result(bn, sn) == 1); + } + } +} + +/** + * @brief These tests just make sure that the return values are consistent with previous versions + * of the functions. The values are not checked for correctness. + * In case the functions are changed and updated to produce 'better' results, the tests should be + * updated accordingly. However, these changes should be documented and justified und users should + * be informed about these changes as they will affect quantification results. + * + */ +TEST_CASE("WCICorrection functions should reproduce previously computed results", TESTTAG) +{ + using Catch::Approx; + + SECTION("wrong tensor shapes should throw") + { + xt::xtensor wci = xt::eval(xt::ones({ 10, 20 })); + xt::xtensor per_beam_offset = xt::eval(xt::linspace(-5.5, 10.5, 10)); + xt::xtensor per_beam_offset_wrong = xt::eval(xt::linspace(-5.5, 10.5, 11)); + xt::xtensor per_sample_offset = xt::eval(xt::linspace(-2.5, 35.2, 20)); + xt::xtensor per_sample_offset_wrong = + xt::eval(xt::linspace(-2.5, 35.2, 9)); + + // wrong beam/sample sizes + REQUIRE_THROWS_AS( + apply_beam_sample_correction(wci, per_beam_offset_wrong, per_sample_offset), + std::invalid_argument); + REQUIRE_THROWS_AS( + inplace_beam_sample_correction(wci, per_beam_offset_wrong, per_sample_offset), + std::invalid_argument); + REQUIRE_THROWS_AS( + apply_beam_sample_correction(wci, per_beam_offset, per_sample_offset_wrong), + std::invalid_argument); + REQUIRE_THROWS_AS( + inplace_beam_sample_correction(wci, per_beam_offset, per_sample_offset_wrong), + std::invalid_argument); + REQUIRE_THROWS_AS( + apply_beam_sample_correction(wci, per_beam_offset_wrong, per_sample_offset_wrong), + std::invalid_argument); + REQUIRE_THROWS_AS( + inplace_beam_sample_correction(wci, per_beam_offset_wrong, per_sample_offset_wrong), + std::invalid_argument); + + // beam sample arrays switched + REQUIRE_THROWS_AS(apply_beam_sample_correction(wci, per_sample_offset, per_beam_offset), + std::invalid_argument); + REQUIRE_THROWS_AS(inplace_beam_sample_correction(wci, per_sample_offset, per_beam_offset), + std::invalid_argument); + + REQUIRE_THROWS_AS(apply_beam_correction(wci, per_sample_offset), std::invalid_argument); + REQUIRE_THROWS_AS(inplace_beam_correction(wci, per_sample_offset), std::invalid_argument); + REQUIRE_THROWS_AS(apply_sample_correction(wci, per_beam_offset), std::invalid_argument); + REQUIRE_THROWS_AS(inplace_sample_correction(wci, per_beam_offset), std::invalid_argument); + } + + SECTION( + "apply_beam_sample_correction and inplace_beam_sample_correction (applied to full array)") + { + xt::xtensor wci = xt::eval(xt::ones({ 10, 20 })); + xt::xtensor per_beam_offset = xt::eval(xt::linspace(-5.5, 10.5, 10)); + auto per_beam_offset_0 = xt::zeros_like(per_beam_offset); + xt::xtensor per_sample_offset = xt::eval(xt::linspace(-2.5, 35.2, 20)); + auto per_sample_offset_0 = xt::zeros_like(per_sample_offset); + + for (size_t mp_cores : { 1, 0, 4 }) + { + // --- apply_beam_sample_correction --- + xt::xtensor result = + xt::eval(apply_beam_sample_correction(wci, per_beam_offset, per_sample_offset)); + check_wci("apply_beam_sample_correction" + fmt::format("({} cores)", mp_cores), + result, + wci, + per_beam_offset, + per_sample_offset); + + result = wci; + inplace_beam_sample_correction( + result, per_beam_offset, per_sample_offset, std::nullopt, std::nullopt, mp_cores); + check_wci("inplace_beam_sample_correction" + fmt::format("({} cores)", mp_cores), + result, + wci, + per_beam_offset, + per_sample_offset); + + // --- apply_sample_correction --- + result = xt::eval( + apply_beam_sample_correction(wci, per_beam_offset_0, per_sample_offset, mp_cores)); + check_wci("apply_sample_correction 1" + fmt::format("({} cores)", mp_cores), + result, + wci, + per_beam_offset_0, + per_sample_offset); + + result = xt::eval(apply_sample_correction(wci, per_sample_offset, mp_cores)); + check_wci("apply_sample_correction 2" + fmt::format("({} cores)", mp_cores), + result, + wci, + per_beam_offset_0, + per_sample_offset); + + result = wci; + inplace_sample_correction( + result, per_sample_offset, std::nullopt, std::nullopt, mp_cores); + check_wci("inplace_sample_correction" + fmt::format("({} cores)", mp_cores), + result, + wci, + per_beam_offset_0, + per_sample_offset); + + // --- apply_beam_correction --- + result = xt::eval( + apply_beam_sample_correction(wci, per_beam_offset, per_sample_offset_0, mp_cores)); + check_wci("apply_beam_correction 1" + fmt::format("({} cores)", mp_cores), + result, + wci, + per_beam_offset, + per_sample_offset_0); + + result = xt::eval(apply_beam_correction(wci, per_beam_offset, mp_cores)); + check_wci("apply_beam_correction 2" + fmt::format("({} cores)", mp_cores), + result, + wci, + per_beam_offset, + per_sample_offset_0); + + result = wci; + inplace_beam_correction(result, per_beam_offset, std::nullopt, std::nullopt, mp_cores); + check_wci("inplace_beam_correction" + fmt::format("({} cores)", mp_cores), + result, + wci, + per_beam_offset, + per_sample_offset_0); + + // --- apply system offset --- + result = xt::eval(apply_system_offset(wci, 12, mp_cores)); + check_wci("apply_system_offset" + fmt::format("({} cores)", mp_cores), + result, + xt::eval(wci + 12), + per_beam_offset_0, + per_sample_offset_0); + + result = wci; + inplace_system_offset(result, 12, std::nullopt, std::nullopt, mp_cores); + check_wci("inplace_system_offset" + fmt::format("({} cores)", mp_cores), + result, + xt::eval(wci + 12), + per_beam_offset_0, + per_sample_offset_0); + } + } + + SECTION("inplace_beam_sample_correction (applied to part of the array)") + { + + xt::xtensor wci = xt::eval(xt::ones({ 10, 20 })); + xt::xtensor per_beam_offset = xt::eval(xt::linspace(-5.5, 10.5, 10)); + auto per_beam_offset_0 = xt::zeros_like(per_beam_offset); + xt::xtensor per_sample_offset = xt::eval(xt::linspace(-2.5, 35.2, 20)); + auto per_sample_offset_0 = xt::zeros_like(per_sample_offset); + + for (size_t mp_cores : { 1, 0, 4 }) + for (std::optional min_bn : + std::vector>{ std::nullopt, 0, 3, 7, 9, 10, 11, 999999 }) + for (std::optional max_bn : + std::vector>{ std::nullopt, 0, 3, 7, 9, 10, 11, 999999 }) + { + + // --- inplace_beam_sample_correction --- + xt::xtensor result = wci; + inplace_beam_sample_correction( + result, per_beam_offset, per_sample_offset, min_bn, max_bn, mp_cores); + check_wci("PARTIAL inplace_beam_sample_correction" + + fmt::format("({} cores)", mp_cores), + result, + wci, + per_beam_offset, + per_sample_offset, + min_bn, + max_bn); + + // ---inplace_sample_correction--- + result = wci; + inplace_sample_correction(result, per_sample_offset, min_bn, max_bn, mp_cores); + check_wci("PARTIAL inplace_sample_correction" + + fmt::format("({} cores)", mp_cores), + result, + wci, + per_beam_offset_0, + per_sample_offset, + min_bn, + max_bn); + + // ---apply_beam_correction--- + result = wci; + inplace_beam_correction(result, per_beam_offset, min_bn, max_bn, mp_cores); + check_wci("PARTIAL inplace_beam_correction" + + fmt::format("({} cores)", mp_cores), + result, + wci, + per_beam_offset, + per_sample_offset_0, + min_bn, + max_bn); + + // --- apply system offset (inplace only) --- + result = wci; + auto modified_wci = wci; + xt::view(modified_wci, + xt::range(min_bn.value_or(0), + max_bn.value_or(modified_wci.shape(0) - 1) + 1), + xt::all()) += 12; + inplace_system_offset(result, 12, min_bn, max_bn, mp_cores); + check_wci("PARTIAL inplace_system_offset" + fmt::format("({} cores)", mp_cores), + result, + modified_wci, + per_beam_offset_0, + per_sample_offset_0); + } + } +} diff --git a/src/tests/meson.build b/src/tests/meson.build index b75fd60..70c56b9 100644 --- a/src/tests/meson.build +++ b/src/tests/meson.build @@ -11,6 +11,9 @@ sources = [ 'signalprocessing/datastructures/cwsignalparameters.test.cpp', 'signalprocessing/datastructures/fmsignalparameters.test.cpp', 'signalprocessing/datastructures/genericsignalparameters.test.cpp', + 'algorithms/functions/absorption.test.cpp', + 'algorithms/functions/rangecorrection.test.cpp', + 'algorithms/functions/wcicorrection.test.cpp', 'geoprocessing/datastructures/beamsampleparameters.test.cpp', 'geoprocessing/datastructures/raytraceresult.test.cpp', 'geoprocessing/datastructures/raytraceresults.test.cpp', diff --git a/src/themachinethatgoesping/algorithms/amplitudecorrection/functions/_shape_assertations.hpp b/src/themachinethatgoesping/algorithms/amplitudecorrection/functions/_shape_assertations.hpp new file mode 100644 index 0000000..e2d3dd4 --- /dev/null +++ b/src/themachinethatgoesping/algorithms/amplitudecorrection/functions/_shape_assertations.hpp @@ -0,0 +1,67 @@ +// SPDX-FileCopyrightText: 2024 Peter Urban, Ghent University +// SPDX-FileCopyrightText: 2022 GEOMAR Helmholtz Centre for Ocean Research Kiel +// +// SPDX-License-Identifier: MPL-2.0 + +/** + * @brief These functions are used for shape verifications in wcicorrection.hpp + * + */ + +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +namespace themachinethatgoesping { +namespace algorithms { +namespace amplitudecorrection { +namespace functions { + +// --- assertions --- + +template +inline void assert_wci_axis_shape(const t_xtensor_2d& wci, + const t_xtensor_1d& per_element_offset, + std::string_view axis_name) +{ + static_assert(tools::helper::c_xtensor_2d, + "Template parameter must be a 2D tensor"); + static_assert(tools::helper::c_xtensor_1d, + "Template parameter must be a 1D tensor"); + + if (axis >= wci.dimension()) + throw std::invalid_argument(fmt::format("ERROR[{}]: axis {} out of range", __func__, axis)); + + if (wci.shape(axis) != per_element_offset.shape(0)) + throw std::invalid_argument(fmt::format("ERROR[{}]: wci.shape({}) [{}] != {}.shape(0) [{}]", + __func__, + axis, + wci.shape(axis), + axis_name, + per_element_offset.shape(0))); +} + +template +inline void assert_wci_beam_sample_shape(const t_xtensor_2d& wci, + const t_xtensor_1d& per_beam_offset, + const t_xtensor_1d& per_sample_offset) +{ + assert_wci_axis_shape<0>(wci, per_beam_offset, "per_beam_offset"); + assert_wci_axis_shape<1>(wci, per_sample_offset, "per_sample_offset"); +} + +} +} +} +} diff --git a/src/themachinethatgoesping/algorithms/amplitudecorrection/functions/absorption.hpp b/src/themachinethatgoesping/algorithms/amplitudecorrection/functions/absorption.hpp index 34719ce..d784733 100644 --- a/src/themachinethatgoesping/algorithms/amplitudecorrection/functions/absorption.hpp +++ b/src/themachinethatgoesping/algorithms/amplitudecorrection/functions/absorption.hpp @@ -23,7 +23,8 @@ namespace functions { * @param T90 * @return double */ -inline double its90_to_its68(double T90) +template +inline t_float its90_to_its68(t_float T90) { // return (T90*T90)/gsw_t90_from_t68(T90); return T90 * 1.00024; diff --git a/src/themachinethatgoesping/algorithms/amplitudecorrection/functions/rangecorrection.hpp b/src/themachinethatgoesping/algorithms/amplitudecorrection/functions/rangecorrection.hpp index f690d09..51e183a 100644 --- a/src/themachinethatgoesping/algorithms/amplitudecorrection/functions/rangecorrection.hpp +++ b/src/themachinethatgoesping/algorithms/amplitudecorrection/functions/rangecorrection.hpp @@ -8,11 +8,12 @@ /* generated doc strings */ #include ".docstrings/rangecorrection.doc.hpp" -/* generated doc strings */ +#include #include #include #include +#include #include namespace themachinethatgoesping { @@ -22,6 +23,20 @@ namespace functions { // template +/** + * @brief Generates a 1D tensor of sample numbers incremented by half. (used for range compuation) + * + * This function calculates a range of sample numbers starting from + * `first_sample_nr + 0.5` to `last_sample_nr + 1.5` with a specified step. + * + * @tparam t_xtensor_1d A 1D tensor type that satisfies the `tools::helper::c_xtensor` concept. + * @tparam t_int An integer type for the sample numbers. + * @param first_sample_nr The starting sample number. + * @param last_sample_nr The ending sample number. + * @param step The step size for the range (default is 1). + * @return A 1D tensor of sample numbers incremented by half. + * @note The template parameter must be a 1D tensor. + */ template inline t_xtensor_1d get_sample_numbers_plus_half(t_int first_sample_nr, t_int last_sample_nr, @@ -36,6 +51,17 @@ inline t_xtensor_1d get_sample_numbers_plus_half(t_int first_sample_nr, return xt::arange(first_sample_nr + t_float(0.5), last_sample_nr + t_float(1.5), step); } +/** + * @brief Approximates the range factor based on the sample interval and a single sound velocity. + * + * This function calculates the range factor using the given sample interval and sound velocity. + * The formula used is: range_factor = sample_interval_s * sound_velocity_m_s * 0.5 + * + * @tparam t_float The floating-point type used for the calculations. + * @param sample_interval_s The sample interval in seconds. + * @param sound_velocity_m_s The sound velocity in meters per second. + * @return The approximated range factor. + */ template inline t_float approximate_range_factor(t_float sample_interval_s, t_float sound_velocity_m_s) { @@ -43,6 +69,21 @@ inline t_float approximate_range_factor(t_float sample_interval_s, t_float sound return sample_interval_s * sound_velocity_m_s * t_float(0.5); } +/** + * @brief Approximates the ranges based on the provided sample interval and a single sound velocity. + * + * This function calculates the approximate ranges for a given set of sample numbers + * by multiplying the sample numbers (plus half) with the approximate range factor. + * + * @tparam t_xtensor_1d A 1D tensor type that satisfies the c_xtensor_1d concept. + * @tparam t_int An integer type for sample numbers. + * @param sample_interval_s The sample interval in seconds. + * @param sound_velocity_m_s The sound velocity in meters per second. + * @param first_sample_nr The first sample number. + * @param last_sample_nr The last sample number. + * @param step The step size between sample numbers (default is 1). + * @return A 1D tensor containing the approximated ranges. + */ template inline t_xtensor_1d approximate_ranges( typename tools::helper::xtensor_datatype::type sample_interval_s, @@ -59,12 +100,27 @@ inline t_xtensor_1d approximate_ranges( approximate_range_factor(sample_interval_s, sound_velocity_m_s); } -template +/** + * @brief Approximates the ranges based on sample interval, a single sound velocity, and sample + * numbers. + * + * This function calculates the approximate ranges by using the provided sample interval, + * sound velocity, and sample numbers. The calculation is performed by adding 0.5 to each + * sample number and then multiplying by the approximate range factor. + * + * @tparam t_xtensor_1d A 1D tensor type that satisfies the tools::helper::c_xtensor concept. + * @tparam t_xtensor_1d_int A 1D tensor type for integers that satisfies the + * tools::helper::c_xtensor concept. + * @param sample_interval_s The interval between samples in seconds. + * @param sound_velocity_m_s The velocity of sound in meters per second. + * @param sample_numbers A 1D tensor containing the sample numbers. + * @return A 1D tensor containing the approximated ranges. + */ +template inline t_xtensor_1d approximate_ranges( typename tools::helper::xtensor_datatype::type sample_interval_s, typename tools::helper::xtensor_datatype::type sound_velocity_m_s, - const t_xtensor_1d_int& sample_numbers) + const t_xtensor_1d_int& sample_numbers) { static_assert(tools::helper::c_xtensor_1d, "Template parameter must be a 1D tensor"); @@ -76,28 +132,50 @@ inline t_xtensor_1d approximate_ranges( approximate_range_factor(sample_interval_s, sound_velocity_m_s); } +/** + * @brief Computes the continuous wave (CW) range correction. + * + * This function calculates the range correction based on the provided ranges, + * absorption coefficient, and time-varying gain (TVG) factor. The range correction + * is computed using the formula: + * + * \f[ + * \text{range correction} = 2 \cdot \text{absorption\_db\_m} \cdot \text{ranges\_m} + + * \text{tvg\_factor} \cdot \log_{10}(\text{ranges\_m}) \f] + * + * If the absorption coefficient is finite and non-zero, it is used in the calculation. + * If the TVG factor is finite and non-zero, it is also included in the calculation. + * If neither the absorption coefficient nor the TVG factor are finite and non-zero, + * the function returns a tensor of zeros with the same shape as the input ranges. + * + * @tparam t_xtensor_1d A 1D tensor type that satisfies the tools::helper::c_xtensor_1d concept. + * @param ranges_m A 1D tensor representing the ranges in meters. + * @param absorption_db_m The absorption coefficient in decibels per meter. + * @param tvg_factor The time-varying gain factor. + * @return A 1D tensor representing the computed range correction. + */ template inline t_xtensor_1d compute_cw_range_correction( - const t_xtensor_1d& ranges_m, - typename tools::helper::xtensor_datatype::type absorption_db_m, - typename tools::helper::xtensor_datatype::type tvg_factor) + const t_xtensor_1d& ranges_m, + std::optional::type> absorption_db_m, + std::optional::type> tvg_factor) { static_assert(tools::helper::c_xtensor_1d, "Template parameter must be a 1D tensor"); - //using t_float = typename tools::helper::xtensor_datatype::type; + // using t_float = typename tools::helper::xtensor_datatype::type; // range correction = absorption*R + tvg_factor*log10(R) - if (tools::helper::float_is_finite_and_not_zero(absorption_db_m)){ - if (tools::helper::float_is_finite_and_not_zero(tvg_factor)) - return (2*absorption_db_m) * ranges_m + tvg_factor * xt::log10(ranges_m); + if (absorption_db_m.has_value()) + { + if (tvg_factor.has_value()) + return (2 * absorption_db_m.value()) * ranges_m + tvg_factor.value() * xt::log10(ranges_m); - return (2*absorption_db_m) * ranges_m; + return (2 * absorption_db_m.value()) * ranges_m; } - if (tools::helper::float_is_finite_and_not_zero(tvg_factor)) - return tvg_factor * xt::log10(ranges_m); - + if (tvg_factor.has_value()) + return tvg_factor.value() * xt::log10(ranges_m); return xt::zeros_like(ranges_m); // range correction = absorption*R + tvg_factor*log10(R) diff --git a/src/themachinethatgoesping/algorithms/amplitudecorrection/functions/wcicorrection.hpp b/src/themachinethatgoesping/algorithms/amplitudecorrection/functions/wcicorrection.hpp index 94cd7b6..b4bdd19 100644 --- a/src/themachinethatgoesping/algorithms/amplitudecorrection/functions/wcicorrection.hpp +++ b/src/themachinethatgoesping/algorithms/amplitudecorrection/functions/wcicorrection.hpp @@ -21,6 +21,7 @@ #include +#include "_shape_assertations.hpp" #include "rangecorrection.hpp" namespace themachinethatgoesping { @@ -28,175 +29,22 @@ namespace algorithms { namespace amplitudecorrection { namespace functions { -// --- assertions --- - -template -inline void assert_wci_axis_shape_0(const t_xtensor_2d& wci, - std::string_view axis_name, - std::optional min_element_index = std::nullopt, - std::optional max_element_index = std::nullopt) -{ - static_assert(tools::helper::c_xtensor_2d, - "Template parameter must be a 2D tensor"); - - if (min_element_index.has_value()) - { - if (max_element_index.has_value()) - { - if (min_element_index.value() > max_element_index.value()) - throw std::invalid_argument( - fmt::format("ERROR[{}/{}]: min_element_index [{}] > max_element_index [{}]", - __func__, - axis_name, - min_element_index.value(), - max_element_index.value())); - - if (wci.shape(0) <= max_element_index.value()) - throw std::invalid_argument( - fmt::format("ERROR[{}/{}]: wci.shape(0) [{}] <= max_element_index [{}]", - __func__, - axis_name, - wci.shape(0), - max_element_index.value())); - } - else - { - if (wci.shape(0) <= min_element_index.value()) - throw std::invalid_argument( - fmt::format("ERROR[{}]: wci.shape(0) [{}] <= min_element_index [{}]", - __func__, - wci.shape(0), - min_element_index.value())); - } - } - else if (max_element_index.has_value()) - { - if (wci.shape(0) <= max_element_index.value()) - throw std::invalid_argument( - fmt::format("ERROR[{}]: wci.shape(0) [{}] <= max_element_index [{}]", - __func__, - wci.shape(0), - max_element_index.value())); - } -} - -template -inline void assert_wci_axis_shape(const t_xtensor_2d& wci, - const t_xtensor_1d& per_element_offset, - std::string_view axis_name, - std::optional min_element_index = std::nullopt, - std::optional max_element_index = std::nullopt) -{ - static_assert(tools::helper::c_xtensor_2d, - "Template parameter must be a 2D tensor"); - static_assert(tools::helper::c_xtensor_1d, - "Template parameter must be a 1D tensor"); - - if (axis >= wci.dimension()) - throw std::invalid_argument(fmt::format("ERROR[{}]: axis {} out of range", __func__, axis)); - - if (min_element_index.has_value()) - { - if (max_element_index.has_value()) - { - // --- min and max index is set --- - if (min_element_index.value() > max_element_index.value()) - throw std::invalid_argument( - fmt::format("ERROR[{}]: min_element_index [{}] > max_element_index [{}]", - __func__, - min_element_index.value(), - max_element_index.value())); - - if (wci.shape(axis) <= max_element_index.value()) - throw std::invalid_argument( - fmt::format("ERROR[{}]: wci.shape({}) [{}] <= max_element_index [{}]", - __func__, - axis, - wci.shape(axis), - max_element_index.value())); - - if (per_element_offset.size() != - max_element_index.value() - min_element_index.value() + 1) - throw std::invalid_argument( - fmt::format("ERROR[{}]: {}.size() [{}] != " - "max_element_index-min_element_index+1 [{}]", - __func__, - axis_name, - per_element_offset.size(), - max_element_index.value() - min_element_index.value() + 1)); - } - else - { - // --- only min index is set --- - if (wci.shape(axis) <= min_element_index.value()) - throw std::invalid_argument( - fmt::format("ERROR[{}]: wci.shape({}) [{}] <= min_element_index [{}]", - __func__, - axis, - wci.shape(axis), - min_element_index.value())); - - if (per_element_offset.size() != wci.shape(axis) - min_element_index.value()) - throw std::invalid_argument( - fmt::format("ERROR[{}]: {}.size() [{}] != " - "wci.shape(axis)-min_element_index [{}]", - __func__, - axis_name, - per_element_offset.size(), - wci.shape(axis) - min_element_index.value())); - } - } - else - { - if (max_element_index.has_value()) - { - // --- only max index is set --- - if (wci.shape(axis) <= max_element_index.value()) - throw std::invalid_argument( - fmt::format("ERROR[{}]: wci.shape({}) [{}] <= max_element_index [{}]", - __func__, - axis, - wci.shape(axis), - max_element_index.value())); - - if (per_element_offset.size() != max_element_index.value() + 1) - throw std::invalid_argument( - fmt::format("ERROR[{}]: {}.size() [{}] != max_element_index+1 [{}]", - __func__, - axis_name, - per_element_offset.size(), - max_element_index.value() + 1)); - } - else - { - // --- no index is set --- - if (wci.shape(axis) != per_element_offset.size()) - throw std::invalid_argument( - fmt::format("ERROR[{}/{}]: wci.shape({}) [{}] != {}.size() [{}]", - __func__, - axis_name, - axis, - wci.shape(axis), - axis_name, - per_element_offset.size())); - } - } -} - -template -inline void assert_wci_beam_sample_shape(const t_xtensor_2d& wci, - const t_xtensor_1d& per_beam_offset, - const t_xtensor_1d& per_sample_offset, - std::optional min_beam_index = std::nullopt, - std::optional max_beam_index = std::nullopt) -{ - assert_wci_axis_shape<0>( - wci, per_beam_offset, "per_beam_offset", min_beam_index, max_beam_index); - assert_wci_axis_shape<1>(wci, per_sample_offset, "per_sample_offset"); -} - // --- apply corrections --- +/** + * @brief Applies beam and sample corrections to the given 2D tensor. + * + * Apply beam and sample corrections to the input 2D tensor. per_beam correction is applied to each + * sample in a beam, per_sample correction to each sample nr for each beam + * + * @tparam t_xtensor_2d Type of the 2D tensor. + * @tparam t_xtensor_1d Type of the 1D tensor. + * @param wci The input 2D tensor to which corrections will be applied. + * @param per_beam_offset A 1D tensor containing the per-beam offsets. + * @param per_sample_offset A 1D tensor containing the per-sample offsets. + * @param mp_cores The number of cores to use for parallel processing (default is 1). + * @return A 2D tensor with the applied beam and sample corrections. + */ template inline t_xtensor_2d apply_beam_sample_correction(const t_xtensor_2d& wci, const t_xtensor_1d& per_beam_offset, @@ -224,36 +72,19 @@ inline void inplace_beam_sample_correction([[maybe_unused]] t_xtensor_2d& wci, std::optional max_beam_index = std::nullopt, int mp_cores = 1) { - assert_wci_beam_sample_shape( - wci, per_beam_offset, per_sample_offset, min_beam_index, max_beam_index); + assert_wci_beam_sample_shape(wci, per_beam_offset, per_sample_offset); - // Apply the range correction to each sample + size_t max_beam_nr = max_beam_index.value_or(wci.shape(0) - 1); + if (max_beam_nr >= wci.shape(0)) + max_beam_nr = wci.shape(0) - 1; + + // Apply the range correction to each sample #pragma omp parallel for num_threads(mp_cores) - for (unsigned int bi = min_beam_index.value_or(0); - bi <= max_beam_index.value_or(per_beam_offset.size() - 1); - ++bi) + for (unsigned int bi = min_beam_index.value_or(0); bi <= max_beam_nr; ++bi) xt::row(wci, bi) += (per_beam_offset.unchecked(bi) + per_sample_offset); // wci = xt::eval(wci); } -template -inline void inplace_beam_sample_correction2([[maybe_unused]] t_xtensor_2d& wci, - const t_xtensor_1d& per_beam_offset, - const t_xtensor_1d& per_sample_offset, - std::optional min_beam_index = std::nullopt, - std::optional max_beam_index = std::nullopt, - int mp_cores = 1) -{ - assert_wci_beam_sample_shape( - wci, per_beam_offset, per_sample_offset, min_beam_index, max_beam_index); - - // Apply the range correction to each sample -#pragma omp parallel for num_threads(mp_cores) - for (unsigned int bi = min_beam_index.value_or(0); - bi <= max_beam_index.value_or(per_beam_offset.size() - 1); - ++bi) - xt::row(wci, bi) += (per_beam_offset.unchecked(bi) + per_sample_offset); -} template inline t_xtensor_2d apply_beam_correction(const t_xtensor_2d& wci, @@ -278,13 +109,14 @@ inline void inplace_beam_correction([[maybe_unused]] t_xtensor_2d& wci, std::optional max_beam_index = std::nullopt, int mp_cores = 1) { - assert_wci_axis_shape<0>( - wci, per_beam_offset, "per_beam_offset", min_beam_index, max_beam_index); + assert_wci_axis_shape<0>(wci, per_beam_offset, "per_beam_offset"); + + size_t max_beam_nr = max_beam_index.value_or(wci.shape(0) - 1); + if (max_beam_nr >= wci.shape(0)) + max_beam_nr = wci.shape(0) - 1; #pragma omp parallel for num_threads(mp_cores) - for (unsigned int bi = min_beam_index.value_or(0); - bi <= max_beam_index.value_or(per_beam_offset.size() - 1); - ++bi) + for (unsigned int bi = min_beam_index.value_or(0); bi <= max_beam_nr; ++bi) xt::row(wci, bi) += per_beam_offset.unchecked(bi); // wci = xt::eval(wci); @@ -315,19 +147,27 @@ inline void inplace_sample_correction([[maybe_unused]] t_xtensor_2d& wci, std::optional max_beam_index = std::nullopt, [[maybe_unused]] int mp_cores = 1) { - assert_wci_axis_shape_0<0>(wci, "per_beam", min_beam_index, max_beam_index); assert_wci_axis_shape<1>(wci, per_sample_offset, "per_sample_offset"); + size_t max_beam_nr = max_beam_index.value_or(wci.shape(0) - 1); + if (max_beam_nr >= wci.shape(0)) + max_beam_nr = wci.shape(0) - 1; + if (mp_cores == 1) { - wci += xt::view(per_sample_offset, xt::newaxis(), xt::all()); + if (min_beam_index.has_value() || max_beam_index.has_value()) + { + xt::view(wci, xt::range(min_beam_index.value_or(0), max_beam_nr + 1), xt::all()) += + xt::view(per_sample_offset, xt::newaxis(), xt::all()); + } + else + wci += xt::view(per_sample_offset, xt::newaxis(), xt::all()); + return; } #pragma omp parallel for num_threads(mp_cores) - for (unsigned int bi = min_beam_index.value_or(0); - bi <= max_beam_index.value_or(wci.shape(0) - 1); - ++bi) + for (unsigned int bi = min_beam_index.value_or(0); bi <= max_beam_nr; ++bi) xt::row(wci, bi) += per_sample_offset; // wci = xt::eval(wci); @@ -358,18 +198,27 @@ inline void inplace_system_offset( std::optional max_beam_index = std::nullopt, int mp_cores = 1) { - assert_wci_axis_shape_0<0>(wci, "per_beam", min_beam_index, max_beam_index); + size_t max_beam_nr = max_beam_index.value_or(wci.shape(0) - 1); + if (max_beam_nr >= wci.shape(0)) + max_beam_nr = wci.shape(0) - 1; if (mp_cores == 1) { - wci += system_offset; - return; + if (min_beam_index.has_value() || max_beam_index.has_value()) + { + xt::view(wci, xt::range(min_beam_index.value_or(0), max_beam_nr + 1), xt::all()) += + system_offset; + return; + } + else + { + wci += system_offset; + return; + } } #pragma omp parallel for num_threads(mp_cores) - for (unsigned int bi = min_beam_index.value_or(0); - bi <= max_beam_index.value_or(wci.shape(0) - 1); - ++bi) + for (unsigned int bi = min_beam_index.value_or(0); bi <= max_beam_nr; ++bi) xt::row(wci, bi) += system_offset; // wci = xt::eval(wci); diff --git a/src/themachinethatgoesping/algorithms/meson.build b/src/themachinethatgoesping/algorithms/meson.build index 8cef8dc..eab13c1 100644 --- a/src/themachinethatgoesping/algorithms/meson.build +++ b/src/themachinethatgoesping/algorithms/meson.build @@ -14,6 +14,7 @@ headers = [ 'helloping.hpp', 'amplitudecorrection/functions.hpp', 'amplitudecorrection/.docstrings/functions.doc.hpp', + 'amplitudecorrection/functions/_shape_assertations.hpp', 'amplitudecorrection/functions/absorption.hpp', 'amplitudecorrection/functions/rangecorrection.hpp', 'amplitudecorrection/functions/wcicorrection.hpp',