Merge pull request #283 from achaikou/assure_4_samples

Assure enough samples for calculation

internal/core/core_attributes_test.go

@@ -1293,3 +1293,119 @@ func TestAttributesOnSamplesNearTraceBoundary(t *testing.T) {
 		require.Equalf(t, expected, *value, "[%v]", attr)
 	}
 }
+
+func TestAttributesInsufficientInterpolationSamples(t *testing.T) {
+	testcases := []struct {
+		name        string
+		top         [][]float32
+		bottom      [][]float32
+		expectedMin [][]float32
+		expectedMax [][]float32
+	}{
+		// top survey boundary at 4.0, bottom boundary at 40.0
+		{
+			name:        "No samples between top and bottom border in the middle",
+			top:         [][]float32{{22.00}},
+			bottom:      [][]float32{{23.00}},
+			expectedMin: [][]float32{{0.00}},
+			expectedMax: [][]float32{{0.25}},
+		},
+		{
+			name:        "No samples between top and bottom border, two to the top survey boundary",
+			top:         [][]float32{{10.00}},
+			bottom:      [][]float32{{11.00}},
+			expectedMin: [][]float32{{-3.00}},
+			expectedMax: [][]float32{{-2.75}},
+		},
+		{
+			name:        "No samples between top and bottom border, one to the top survey boundary",
+			top:         [][]float32{{6.00}},
+			bottom:      [][]float32{{7.00}},
+			expectedMin: [][]float32{{-4.00}},
+			expectedMax: [][]float32{{-3.75}},
+		},
+		{
+			name:        "One sample between top and bottom border, one to the top survey boundary",
+			top:         [][]float32{{6.00}},
+			bottom:      [][]float32{{11.00}},
+			expectedMin: [][]float32{{-4.00}},
+			expectedMax: [][]float32{{-2.75}},
+		},
+		{
+			name:        "Top boundary on the top survey boundary, no additional samples between borders",
+			top:         [][]float32{{4.00}},
+			bottom:      [][]float32{{7.00}},
+			expectedMin: [][]float32{{-4.50}},
+			expectedMax: [][]float32{{-3.75}},
+		},
+		{
+			name:        "Top boundary on the top survey boundary + additional sample between borders",
+			top:         [][]float32{{4.00}},
+			bottom:      [][]float32{{11.00}},
+			expectedMin: [][]float32{{-4.50}},
+			expectedMax: [][]float32{{-2.75}},
+		},
+		{
+			name:        "No samples between top and bottom border, one to the bottom survey boundary",
+			top:         [][]float32{{36.50}},
+			bottom:      [][]float32{{39.00}},
+			expectedMin: [][]float32{{3.625}},
+			expectedMax: [][]float32{{4.25}},
+		},
+		{
+			name:        "One sample between top and bottom border, one to the bottom survey boundary",
+			top:         [][]float32{{35.50}},
+			bottom:      [][]float32{{39.00}},
+			expectedMin: [][]float32{{3.375}},
+			expectedMax: [][]float32{{4.25}},
+		},
+		{
+			name:        "Bottom boundary on the bottom survey boundary, no additional samples between borders",
+			top:         [][]float32{{39.50}},
+			bottom:      [][]float32{{40.00}},
+			expectedMin: [][]float32{{4.375}},
+			expectedMax: [][]float32{{4.5}},
+		},
+	}
+
+	targetAttributes := []string{"min_at", "min", "max_at", "max"}
+	const stepsize = float32(0.5)
+	interpolationMethod, _ := GetInterpolationMethod("linear")
+
+	for _, testcase := range testcases {
+		topSurface := samples10Surface(testcase.top)
+		bottomSurface := samples10Surface(testcase.bottom)
+
+		handle, _ := NewDSHandle(samples10)
+		defer handle.Close()
+		buf, boundsErr := handle.GetAttributesBetweenSurfaces(
+			topSurface,
+			bottomSurface,
+			stepsize,
+			targetAttributes,
+			interpolationMethod,
+		)
+
+		require.NoErrorf(t, boundsErr,
+			"[%s] Expected enough samples between top %v and bottom %v",
+			testcase.name,
+			testcase.top,
+			testcase.bottom,
+		)
+
+		require.Len(t, buf, len(targetAttributes), "Wrong number of attributes")
+
+		checkAttributes := func(buf [][]byte, attrn int, expected [][]float32, label string) {
+			for i, attr := range buf[len(buf)/4*attrn : len(buf)/4*(attrn+1)] {
+				result, err := toFloat32(attr)
+				require.NoErrorf(t, err, "Couldn't convert to float32")
+
+				require.Equalf(t, expected[i], *result, "[%v]: Wrong %s values", testcase.name, label)
+			}
+		}
+		expected := [][][]float32{testcase.top, testcase.expectedMin, testcase.bottom, testcase.expectedMax}
+		for i, attr := range targetAttributes {
+			checkAttributes(buf, i, expected[i], attr)
+		}
+	}
+}

internal/core/subvolume.cpp

@@ -1,3 +1,4 @@
+#include <cassert>
 #include <cmath>
 #include <stdexcept>
 
@@ -133,13 +134,39 @@ SurfaceBoundedSubVolume* make_subvolume(
         };
 
         std::int8_t top_margin = calculate_margin(Border::Top);
-        if (top_margin != segment_blueprint.preferred_margin()) {
-            subvolume->m_segment_top_margins.emplace(i, top_margin);
-        }
+        bool is_top_margin_atypical = (top_margin != segment_blueprint.preferred_margin());
 
         std::int8_t bottom_margin = calculate_margin(Border::Bottom);
-        if (bottom_margin != segment_blueprint.preferred_margin()) {
-            subvolume->m_segment_bottom_margins.emplace(i, bottom_margin);
+        bool is_bottom_margin_atypical = (bottom_margin != segment_blueprint.preferred_margin());
+
+        // limitation from makima samples interpolation algorithm
+        const int min_samples = 4;
+        assert(
+            (void("Current logic relies on relationship between min_samples and preferred_margin"),
+             min_samples == 2 * segment_blueprint.preferred_margin())
+        );
+        auto size = segment_blueprint.size(top_depth, bottom_depth, top_margin, bottom_margin);
+        if (size < min_samples) {
+            if (is_top_margin_atypical && is_bottom_margin_atypical) {
+                throw std::runtime_error(
+                    "Segment size is too small. Top margin: " +
+                    std::to_string(top_margin) + ", bottom margin: " +
+                    std::to_string(bottom_margin)
+                );
+            }
+
+            int diff = min_samples - size;
+            if (is_top_margin_atypical) {
+                bottom_margin += diff;
+                is_bottom_margin_atypical = true;
+            } else {
+                top_margin += diff;
+                is_top_margin_atypical = true;
+            }
+        }
+
+        if (is_top_margin_atypical) {
+            subvolume->m_segment_top_margins.emplace(i, top_margin);
         }
 
         subvolume->m_segment_offsets[i + 1] =
@@ -156,7 +183,7 @@ void SurfaceBoundedSubVolume::reinitialize(
 ) const {
     segment.reinitialize(
         m_ref[index], m_top[index], m_bottom[index],
-        top_margin(index), bottom_margin(index),
+        top_margin(index),
         m_data.begin() + m_segment_offsets[index], m_data.begin() + m_segment_offsets[index + 1]
     );
 }

internal/core/subvolume.hpp

@@ -323,7 +323,6 @@ class RawSegment : public Segment {
         const float top_boundary,
         const float bottom_boundary,
         std::uint8_t top_margin,
-        std::uint8_t bottom_margin,
         std::vector<float>::const_iterator data_begin,
         std::vector<float>::const_iterator data_end,
         RawSegmentBlueprint const* blueprint
@@ -332,22 +331,19 @@ class RawSegment : public Segment {
           m_blueprint(blueprint),
           m_data_begin(data_begin),
           m_data_end(data_end),
-          m_top_margin(top_margin),
-          m_bottom_margin(bottom_margin)
+          m_top_margin(top_margin)
           {}
 
     void reinitialize(
         float reference,
         float top_boundary,
         float bottom_boundary,
         std::uint8_t top_margin,
-        std::uint8_t bottom_margin,
         std::vector<float>::const_iterator data_begin,
         std::vector<float>::const_iterator data_end
     ) noexcept {
         Segment::reinitialize(reference, top_boundary, bottom_boundary);
         this->m_top_margin = top_margin;
-        this->m_bottom_margin = bottom_margin;
         this->m_data_begin = data_begin;
         this->m_data_end = data_end;
     }
@@ -380,7 +376,6 @@ class RawSegment : public Segment {
     std::vector<float>::const_iterator m_data_begin;
     std::vector<float>::const_iterator m_data_end;
     std::uint8_t m_top_margin;
-    std::uint8_t m_bottom_margin;
 };
 
 /**
@@ -476,19 +471,12 @@ class SurfaceBoundedSubVolume {
                    : this->m_segment_blueprint.preferred_margin();
     }
 
-    std::uint8_t bottom_margin(std::size_t index) const {
-        return this->m_segment_bottom_margins.count(index)
-                   ? this->m_segment_bottom_margins.at(index)
-                   : this->m_segment_blueprint.preferred_margin();
-    }
-
     RawSegment vertical_segment(std::size_t index) const noexcept {
         return RawSegment(
             this->m_ref[index],
             this->m_top[index],
             this->m_bottom[index],
             this->top_margin(index),
-            this->bottom_margin(index),
             m_data.begin() + m_segment_offsets[index],
             m_data.begin() + m_segment_offsets[index + 1],
             &this->m_segment_blueprint
@@ -551,7 +539,6 @@ class SurfaceBoundedSubVolume {
      * that are different from preferred blueprint margin.
      */
     std::unordered_map<std::size_t, std::uint8_t> m_segment_top_margins;
-    std::unordered_map<std::size_t, std::uint8_t> m_segment_bottom_margins;
 
     RegularSurface const& m_ref;
     RegularSurface const& m_top;

tests/gtest/test_utils.cpp

@@ -258,7 +258,7 @@ void Datahandle10SamplesTest::check_attribute(
     // make exclusive
     high_index[0] += 1;
     high_index[1] += 1;
-    high_index[2] += subvolume.bottom_margin(0) + 1;
+    high_index[2] = low_index[2] + subvolume.vertical_segment(0).size();
 
     std::size_t nr_of_values = subvolume.nsamples(
         0, (high_index[0] - low_index[0]) * (high_index[1] - low_index[1])