Skip to content

Commit

Permalink
Merge pull request #361 from clEsperanto/fix-otsu
Browse files Browse the repository at this point in the history 
Fix-otsu
  • Loading branch information
@StRigaud
StRigaud authored Oct 14, 2024
2 parents 3a357b4 + 5bb5c90 commit 2326bfe
Show file tree
Hide file tree
Showing 7 changed files with 189 additions and 63 deletions.
2 changes: 1 addition & 1 deletion CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@ cmake_minimum_required(VERSION 3.20)

project(CLIc VERSION 0.13.3)

set(kernel_version_tag "3.0.0" CACHE STRING "clEsperanto kernel version tag")
set(kernel_version_tag "3.0.2" CACHE STRING "clEsperanto kernel version tag")
set(eigen_lib_version_tag "3.4.0" CACHE STRING "Eigen library version tag")

# if not set, set the default build type to Release
Expand Down
12 changes: 6 additions & 6 deletions clic/src/statistics.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -55,8 +55,8 @@ _std_per_label(const Device::Pointer & device,
const Array::Pointer & intensity,
int nb_labels) -> Array::Pointer
{
const size_t height = label->height();
const size_t depth = label->depth();
const auto height = label->height();
const auto depth = label->depth();

auto label_statistics_stack = Array::create(nb_labels, height, 6, 3, dType::FLOAT, mType::BUFFER, device);
label_statistics_stack->fill(0);
Expand All @@ -82,9 +82,9 @@ compute_statistics_per_labels(const Device::Pointer & device,


// initialize variables, output, and constants
const size_t offset = 1;
const size_t nb_labels = static_cast<size_t>(tier2::maximum_of_all_pixels_func(device, label)) + offset;
const size_t nb_measurements = nb_labels - offset;
const auto offset = 1;
const auto nb_labels = static_cast<size_t>(tier2::maximum_of_all_pixels_func(device, label)) + offset;
const auto nb_measurements = nb_labels - offset;
const RangeArray origin = { 0, 0, 0 };
const RangeArray region = { nb_measurements, 1, 1 };

Expand All @@ -105,7 +105,7 @@ compute_statistics_per_labels(const Device::Pointer & device,

// 0 labels
std::vector<float> labels_list(nb_measurements);
std::iota(labels_list.begin(), labels_list.end(), offset);
std::iota(labels_list.begin(), labels_list.end(), static_cast<float>(offset));
region_props["label"] = std::move(labels_list);

// 1-6 bbox x, y, z, width, height, depth
Expand Down
152 changes: 112 additions & 40 deletions clic/src/tier4/threshold_otsu.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,62 +5,134 @@
#include "tier4.hpp"

#include "utils.hpp"
#include <functional>
#include <numeric>

namespace cle::tier4
{

//
// auto
// threshold_otsu_func(const Device::Pointer & device, const Array::Pointer & src, Array::Pointer dst) -> Array::Pointer
// {
// constexpr int bin = 256;
// const float min_intensity = tier2::minimum_of_all_pixels_func(device, src);
// const float max_intensity = tier2::maximum_of_all_pixels_func(device, src);
// auto hist_array = Array::create(bin, 1, 1, 1, dType::FLOAT, mType::BUFFER, src->device());
// tier3::histogram_func(device, src, hist_array, bin, min_intensity, max_intensity);
// std::vector<float> histogram_array(hist_array->size());
// hist_array->readTo(histogram_array.data());
// double threshold = -1;
// double max_variance = -1;
// double variance = 0;
// double sum_1 = 0;
// double sum_2 = 0;
// double weight_1 = 0;
// double weight_2 = 0;
// double mean_1 = 0;
// double mean_2 = 0;
// const double nb_pixels = src->size();
// const double intensity_factor = static_cast<double>(max_intensity - min_intensity);
// // const double intensity_factor = static_cast<double>(max_intensity - min_intensity) / static_cast<double>(bin -
// 1);
// // old implementation - to be removed
// // std::vector<float> range(histogram_array.size());
// // std::iota(range.begin(), range.end(), 0.0);
// // std::transform(range.begin(), range.end(), range.begin(), [intensity_factor, min_intensity](float intensity) {
// // return intensity * intensity_factor + static_cast<float>(min_intensity);
// // });
// //
// std::vector<double> range(bin);
// std::iota(range.begin(), range.end(), 0.0f);
// std::transform(range.begin(), range.end(), range.begin(), [intensity_factor, min_intensity](double value) {
// return (value * intensity_factor) / (bin - 1) + static_cast<double>(min_intensity);
// });
// sum_1 = std::transform_reduce(
// range.begin(), range.end(), histogram_array.begin(), 0.0f, std::plus<>(), [](double intensity, float hist_value)
// {
// return intensity * static_cast<double>(hist_value);
// });
// for (size_t index = 0; index < range.size(); ++index)
// {
// if (histogram_array[index] == 0)
// {
// continue;
// }
// weight_1 += histogram_array[index];
// weight_2 = nb_pixels - weight_1;
// sum_2 += range[index] * static_cast<double>(histogram_array[index]);
// mean_1 = sum_2 / weight_1;
// mean_2 = (sum_1 - sum_2) / weight_2;
// variance = weight_1 * weight_2 * ((mean_1 - mean_2) * (mean_1 - mean_2));
// if (variance > max_variance)
// {
// threshold = range[index];
// max_variance = variance;
// }
// }
// std::cout << "Threshold: " << threshold << std::endl;
// tier0::create_like(src, dst, dType::BINARY);
// return tier1::greater_constant_func(device, src, dst, static_cast<float>(threshold));
// }


auto
threshold_otsu_func(const Device::Pointer & device, const Array::Pointer & src, Array::Pointer dst) -> Array::Pointer
{
constexpr int bin = 256;
const float min_intensity = tier2::minimum_of_all_pixels_func(device, src);
const float max_intensity = tier2::maximum_of_all_pixels_func(device, src);
auto hist_array = Array::create(bin, 1, 1, 1, dType::FLOAT, mType::BUFFER, src->device());

auto hist_array = Array::create(bin, 1, 1, 1, dType::FLOAT, mType::BUFFER, src->device());
tier3::histogram_func(device, src, hist_array, bin, min_intensity, max_intensity);
std::vector<float> histogram_array(hist_array->size());
hist_array->readTo(histogram_array.data());
float threshold = -1;
float max_variance = -1;
float variance = 0;
float sum_1 = 0;
float sum_2 = 0;
float weight_1 = 0;
float weight_2 = 0;
float mean_1 = 0;
float mean_2 = 0;
const float nb_pixels = src->size();
const float intensity_factor = (max_intensity - min_intensity) / (bin - 1);

std::vector<float> range(histogram_array.size());
std::iota(range.begin(), range.end(), 0.0f);
std::transform(range.begin(), range.end(), range.begin(), [intensity_factor, min_intensity](float intensity) {
return intensity * intensity_factor + min_intensity;
});
sum_1 = std::transform_reduce(
range.begin(), range.end(), histogram_array.begin(), 0.0f, std::plus<>(), [](float intensity, float hist_value) {
return intensity * hist_value;
std::vector<float> counts(hist_array->size());
hist_array->readTo(counts.data());


double range = max_intensity - min_intensity;

std::vector<double> bin_centers(bin);
std::iota(bin_centers.begin(), bin_centers.end(), 0.0);
std::transform(
bin_centers.begin(), bin_centers.end(), bin_centers.begin(), [range, min_intensity, bin](double value) {
return (value * range) / (bin - 1) + static_cast<double>(min_intensity);
});
for (size_t index = 0; index < range.size(); ++index)

std::vector<double> weight1(bin), weight2(bin), mean1(bin), mean2(bin), variance12(bin - 1);

// Compute weight1
std::partial_sum(counts.begin(), counts.end(), weight1.begin());

// Compute weight2
std::vector<double> reversed_counts(counts.rbegin(), counts.rend());
std::partial_sum(reversed_counts.begin(), reversed_counts.end(), weight2.rbegin());

// Compute mean1
std::vector<double> counts_bin_centers(bin);
std::transform(counts.begin(), counts.end(), bin_centers.begin(), counts_bin_centers.begin(), std::multiplies<>());
std::partial_sum(counts_bin_centers.begin(), counts_bin_centers.end(), mean1.begin());
std::transform(mean1.begin(), mean1.end(), weight1.begin(), mean1.begin(), std::divides<>());

// Compute mean2
std::partial_sum(counts_bin_centers.rbegin(), counts_bin_centers.rend(), mean2.rbegin());
std::transform(mean2.begin(), mean2.end(), weight2.begin(), mean2.begin(), std::divides<>());

// Compute variance12
for (size_t i = 0; i < bin - 1; ++i)
{
if (histogram_array[index] == 0)
{
continue;
}
weight_1 += histogram_array[index];
weight_2 = nb_pixels - weight_1;
sum_2 += range[index] * histogram_array[index];
mean_1 = sum_2 / weight_1;
mean_2 = (sum_1 - sum_2) / weight_2;
variance = weight_1 * weight_2 * ((mean_1 - mean_2) * (mean_1 - mean_2));
if (variance > max_variance)
{
threshold = range[index];
max_variance = variance;
}
variance12[i] = weight1[i] * weight2[i + 1] * (mean1[i] - mean2[i + 1]) * (mean1[i] - mean2[i + 1]);
}

auto max_it = std::max_element(variance12.begin(), variance12.end());
size_t idx = std::distance(variance12.begin(), max_it);

double threshold = bin_centers[idx];

std::cout << "Threshold: " << threshold << std::endl;

tier0::create_like(src, dst, dType::BINARY);
return tier1::greater_constant_func(device, src, dst, threshold);
return tier1::greater_constant_func(device, src, dst, static_cast<float>(threshold));
}


} // namespace cle::tier4
2 changes: 1 addition & 1 deletion clic/src/tier7/eroded_otsu_labeling.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,7 @@ eroded_otsu_labeling_func(const Device::Pointer & device,
binary->copyTo(eroded1);
for (int i = 0; i < number_of_erosions; i++)
{
tier1::erode_box_func(device, eroded1, eroded2);
tier1::erode_func(device, eroded1, eroded2, "box");
std::swap(eroded1, eroded2);
}
return tier6::masked_voronoi_labeling_func(device, eroded1, binary, dst);
Expand Down
2 changes: 1 addition & 1 deletion clic/src/tier7/voronoi_otsu_labeling.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,7 @@ voronoi_otsu_labeling_func(const Device::Pointer & device,
{
tier0::create_like(src, dst, dType::LABEL);
auto temp = tier1::gaussian_blur_func(device, src, nullptr, spot_sigma, spot_sigma, spot_sigma);
auto spot = tier2::detect_maxima_box_func(device, temp, nullptr, 0, 0, 0);
auto spot = tier2::detect_maxima_func(device, temp, nullptr, 0, 0, 0, "box");
temp = tier1::gaussian_blur_func(device, src, nullptr, outline_sigma, outline_sigma, outline_sigma);
auto segmentation = tier4::threshold_otsu_func(device, temp, nullptr);
auto binary = tier1::binary_and_func(device, spot, segmentation, nullptr);
Expand Down
57 changes: 57 additions & 0 deletions tests/tier4/test_threshold_otsu.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,22 @@ class TestThresholdOtsu : public ::testing::TestWithParam<std::string>
const std::array<float, 3 * 2 * 2> input = { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12 }; // 6.0058594 skimage
const std::array<uint8_t, 3 * 2 * 2> valid = { 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1 }; // 6.00392 clic
std::array<uint8_t, 3 * 2 * 2> output;


const std::array<float, 3 * 3 * 1> input_low = { 0, 0, 0, 0, 0.003, 0, 0, 0, 0 };
const std::array<uint8_t, 3 * 3 * 1> valid_low = { 0, 0, 0, 0, 1, 0, 0, 0, 0 };
std::array<uint8_t, 3 * 3 * 1> output_low;


const std::vector<float> test = { 0.2295937, 0.92521435, 1.8537747, 1.8731794, 0.9393289, 0.22589178, 0.3740187,
1.3714234, 2.6094873, 2.590501, 1.3086381, 0.32493836, 0.2948968, 0.9629241,
1.6921631, 1.6221197, 0.826933, 0.21854304, 0.1166786, 0.36253795, 0.68537563,
0.8279946, 0.6108529, 0.26835373, 0.04627159, 0.25866497, 0.8760503, 1.6040901,
1.5428433, 0.8239908, 0.0561676, 0.3889102, 1.3680687, 2.4913776, 2.392372,
1.2908775, 0.06311502, 0.3839612, 1.1596216, 1.9070312, 1.7505941, 0.9107878 };

const std::vector<float> valid_test = { 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0 };
};

TEST_P(TestThresholdOtsu, execute)
Expand All @@ -31,6 +47,47 @@ TEST_P(TestThresholdOtsu, execute)
}
}

TEST_P(TestThresholdOtsu, lowIntensity)
{
std::string param = GetParam();
cle::BackendManager::getInstance().setBackend(param);
auto device = cle::BackendManager::getInstance().getBackend().getDevice("", "all");
device->setWaitToFinish(true);

auto gpu_input = cle::Array::create(3, 3, 1, 2, cle::dType::FLOAT, cle::mType::BUFFER, device);
gpu_input->writeFrom(input_low.data());

auto gpu_output = cle::tier4::threshold_otsu_func(device, gpu_input, nullptr);

gpu_output->readTo(output_low.data());
for (int i = 0; i < output_low.size(); i++)
{
EXPECT_EQ(output_low[i], valid_low[i]);
}
}


TEST_P(TestThresholdOtsu, gauss)
{
std::string param = GetParam();
cle::BackendManager::getInstance().setBackend(param);
auto device = cle::BackendManager::getInstance().getBackend().getDevice("", "all");
device->setWaitToFinish(true);

auto gpu_input = cle::Array::create(6, 7, 1, 2, cle::dType::FLOAT, cle::mType::BUFFER, device);
gpu_input->writeFrom(test.data());

auto gpu_output = cle::tier4::threshold_otsu_func(device, gpu_input, nullptr);

std::vector<uint8_t> output_test(gpu_output->size());
gpu_output->readTo(output_test.data());
for (int i = 0; i < output_test.size(); i++)
{
EXPECT_EQ(output_test[i], valid_test[i]);
}
}


std::vector<std::string>
getParameters()
{
Expand Down
25 changes: 11 additions & 14 deletions tests/tier7/test_voronoi_otsu_labeling.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,19 +6,16 @@
class TestVoronoiOtsuLabeling : public ::testing::TestWithParam<std::string>
{
protected:
const std::array<float, 7 * 7 * 3> input = {
0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 1, 1,
1, 1, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 3, 1, 0, 0, 0, 1, 5, 8, 2, 0, 0, 0, 1, 1, 5, 1, 1, 1, 1, 1,
1, 1, 1, 1, 5, 1, 1, 0, 0, 0, 2, 8, 5, 1, 0, 0, 0, 1, 3, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1,
1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0
};
const std::array<uint32_t, 7 * 7 * 3> valid = {
0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 1, 1, 1, 2, 2, 2, 2, 0, 0, 0, 2, 2,
2, 2, 0, 0, 0, 2, 2, 2, 2, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 1, 0, 0, 2, 1, 1, 1, 1, 2, 2, 2, 2,
1, 1, 1, 2, 2, 2, 2, 2, 0, 0, 2, 2, 2, 2, 0, 0, 0, 2, 2, 2, 2, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1,
1, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 1, 1, 1, 2, 2, 2, 2, 0, 0, 0, 2, 2, 2, 2, 0, 0, 0, 2, 2, 2, 2, 0, 0, 0
const std::array<float, 6 * 7 * 1> input = {

0, 0, 1, 1, 0, 0, 0, 1, 8, 9, 1, 0, 0, 1, 7, 6, 1, 0, 0, 0, 1,
1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 1, 8, 7, 1, 0, 0, 1, 1, 1, 0


};
std::array<uint32_t, 7 * 7 * 3> output;
const std::array<uint32_t, 6 * 7 * 1> valid = { 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 0, 0, 0, 1,
2, 2, 0, 0, 0, 0, 2, 2, 0, 0, 0, 2, 2, 2, 0, 0, 0, 0, 2, 2, 0 };
std::array<uint32_t, 6 * 7 * 1> output;
};

TEST_P(TestVoronoiOtsuLabeling, execute)
Expand All @@ -28,10 +25,10 @@ TEST_P(TestVoronoiOtsuLabeling, execute)
auto device = cle::BackendManager::getInstance().getBackend().getDevice("", "all");
device->setWaitToFinish(true);

auto gpu_input = cle::Array::create(7, 7, 3, 3, cle::dType::FLOAT, cle::mType::BUFFER, device);
auto gpu_input = cle::Array::create(6, 7, 1, 1, cle::dType::FLOAT, cle::mType::BUFFER, device);
gpu_input->writeFrom(input.data());

auto gpu_output = cle::tier7::voronoi_otsu_labeling_func(device, gpu_input, nullptr, 0, 1);
auto gpu_output = cle::tier7::voronoi_otsu_labeling_func(device, gpu_input, nullptr, 1, 1);

gpu_output->readTo(output.data());
for (int i = 0; i < output.size(); i++)
Expand Down

0 comments on commit 2326bfe

Please sign in to comment.