diff --git a/Resample/Slice/KzComputation.cpp b/Resample/Slice/KzComputation.cpp
index ef5d804395632c2204f1323b2b72b77d88f3781c..20ed4589a8dc0a354bb23b042a584e4994db515a 100644
--- a/Resample/Slice/KzComputation.cpp
+++ b/Resample/Slice/KzComputation.cpp
@@ -44,15 +44,17 @@ complex_t checkForUnderflow(complex_t val)
 std::vector<complex_t> Compute::Kz::computeReducedKz(const SliceStack& slices, kvector_t k)
 {
     const size_t N = slices.size();
-    const double n_ref = slices[0].material().refractiveIndex(M_TWOPI / k.mag()).real();
+
+    const size_t i_ref = k.z() > 0. ? N-1 : 0;
+    const double n_ref = slices[i_ref].material().refractiveIndex(M_TWOPI / k.mag()).real();
     const double k_base = k.mag() * (k.z() > 0.0 ? -1 : 1);
 
     std::vector<complex_t> result(N);
     // Calculate refraction angle, expressed as k_z, for each layer.
-    complex_t rad = slices[0].scalarReducedPotential(k, n_ref);
-    result[0] = k_base * std::sqrt(rad);
-    for (size_t i = 1; i < N; ++i) {
-        rad = checkForUnderflow(slices[i].scalarReducedPotential(k, n_ref));
+    for (size_t i = 0; i < N; ++i) {
+        complex_t rad = slices[i].scalarReducedPotential(k, n_ref);
+        if (i!=i_ref)
+            rad = checkForUnderflow(rad);
         result[i] = k_base * std::sqrt(rad);
     }
     return result;
diff --git a/Tests/ReferenceData/Std/RectDetWithRoi.int.gz b/Tests/ReferenceData/Std/RectDetWithRoi.int.gz
index b699596c86db5ffb1657a09cddd7a2c4b6493b52..0bea92ef6de9726f96038b4b2666b39828e36889 100644
Binary files a/Tests/ReferenceData/Std/RectDetWithRoi.int.gz and b/Tests/ReferenceData/Std/RectDetWithRoi.int.gz differ
diff --git a/Tests/UnitTests/Core/Fresnel/KzComputationTest.cpp b/Tests/UnitTests/Core/Fresnel/KzComputationTest.cpp
index dd398ed5c531b1104a2bad3cffbc83daedc7be5a..315c6b1d203ba6ffcd2baf19cd17d8fc4d4a15fe 100644
--- a/Tests/UnitTests/Core/Fresnel/KzComputationTest.cpp
+++ b/Tests/UnitTests/Core/Fresnel/KzComputationTest.cpp
@@ -33,7 +33,7 @@ TEST_F(KzComputationTest, initial)
     mLayer.addLayer(layer3);
     mLayer.addLayer(layer4);
 
-    kvector_t k = vecOfLambdaAlphaPhi(1.0, 1.0 * Units::deg, 0.0);
+    kvector_t k = vecOfLambdaAlphaPhi(1.0, -1.0 * Units::deg, 0.0);
 
     const auto re_sample = ProcessedSample::make(mLayer, {});