diff --git a/include/merian-nodes/nodes/accumulate/accumulate.hpp b/include/merian-nodes/nodes/accumulate/accumulate.hpp
index bf9e51f..0c57103 100644
--- a/include/merian-nodes/nodes/accumulate/accumulate.hpp
+++ b/include/merian-nodes/nodes/accumulate/accumulate.hpp
@@ -37,6 +37,8 @@ class Accumulate : public Node {
 
         float adaptive_alpha_reduction = 0.0;
         float adaptive_alpha_ipr_factor = 1.5;
+
+        uint32_t iteration = 0;
     };
 
   public:
diff --git a/include/merian-shaders/reprojection.glsl b/include/merian-shaders/reprojection.glsl
index 06cd449..3c00c1d 100644
--- a/include/merian-shaders/reprojection.glsl
+++ b/include/merian-shaders/reprojection.glsl
@@ -43,4 +43,29 @@ bool reprojection_intersect_border(inout vec2 prev_pos, const vec2 mv, const vec
     return false;
 }
 
+ivec2 reproject_pixel_nearest(const vec2 pixel) {
+    return ivec2(round(pixel));
+}
+
+// Performs stochastic bilinear interpolation when reprojecting
+ivec2 reproject_pixel_stochastic(const vec2 pixel, const float random) {
+    const vec2 relative_pos = fract(pixel);
+
+    // (0, 0)
+    float bary_sum = relative_pos.x * relative_pos.y;
+    if (random <= bary_sum)
+        return ivec2(ceil(pixel));
+
+    bary_sum += relative_pos.x * (1. - relative_pos.y);
+    if (random <= bary_sum)
+        return ivec2(0, 1) * ivec2(floor(pixel)) + (1 - ivec2(0, 1)) * ivec2(ceil(pixel));
+
+    bary_sum += (1. - relative_pos.x) * relative_pos.y;
+    if (random <= bary_sum)
+        return ivec2(1, 0) * ivec2(floor(pixel)) + (1 - ivec2(1, 0)) * ivec2(ceil(pixel));
+
+    // (1, 1)
+    return ivec2(floor(pixel));
+}
+
 #endif
diff --git a/src/merian-nodes/nodes/accumulate/accumulate.comp b/src/merian-nodes/nodes/accumulate/accumulate.comp
index b27e043..75eac77 100644
--- a/src/merian-nodes/nodes/accumulate/accumulate.comp
+++ b/src/merian-nodes/nodes/accumulate/accumulate.comp
@@ -7,6 +7,7 @@
 
 #define FILTER_MODE_NEAREST 0
 #define FILTER_MODE_BILINEAR 1
+#define FILTER_MODE_STOCHASTIC_BILINEAR 2
 
 layout(constant_id = 2) const uint WG_ROUNDED_IRR_SIZE_X = 1;
 layout(constant_id = 3) const uint WG_ROUNDED_IRR_SIZE_Y = 1;
@@ -31,20 +32,23 @@ layout(push_constant, std140) uniform pc_t {
 
     float adaptive_alpha_reduction;
     float adaptive_alpha_ipr_factor;
+
+    uint iteration;
 } pc;
 
 #include "merian-shaders/normal_encode.glsl"
 #include "merian-shaders/reprojection.glsl"
+#include "merian-shaders/hash.glsl"
+#include "merian-shaders/random.glsl"
 #include "merian-shaders/color/colors_yuv.glsl"
 
-void get_prev_nearest(const vec2 prev_pos,
+void get_prev_nearest(const ivec2 prev_ipos,
                       const float linear_z,
                       const float vel_z,
                       const vec3 normal,
                       inout vec4 prev_irr_histlen,
                       inout vec2 prev_moments,
                       inout float sum_w) {
-  const ivec2 prev_ipos = ivec2(round(prev_pos));
   if(any(greaterThanEqual(prev_ipos, textureSize(img_prev_accum, 0))) || any(lessThan(prev_ipos, ivec2(0))))
     return;
 
@@ -180,15 +184,15 @@ void main() {
     const GBuffer gbuf = gbuffer[gbuffer_index(ipos, imageSize(img_accum))];
     const vec3 normal = geo_decode_normal(gbuf.enc_normal);
 
-    if (FILTER_MODE == FILTER_MODE_NEAREST)
-      get_prev_nearest(prev_pos,
+    if (FILTER_MODE == FILTER_MODE_NEAREST) {
+      get_prev_nearest(ivec2(round(prev_pos)),
                        gbuf.linear_z,
                        gbuf.vel_z,
                        normal,
                        prev_irr_histlen,
                        prev_moments,
                        sum_w);
-    else if (FILTER_MODE == FILTER_MODE_BILINEAR)
+    } else if (FILTER_MODE == FILTER_MODE_BILINEAR) {
       get_prev_bilinear(prev_pos,
                         gbuf.linear_z,
                         gbuf.vel_z,
@@ -196,6 +200,16 @@ void main() {
                         prev_irr_histlen,
                         prev_moments,
                         sum_w);
+    } else if (FILTER_MODE == FILTER_MODE_STOCHASTIC_BILINEAR) {
+      uint rng_state = pcg3d16(uvec3(ipos, pc.iteration));
+      get_prev_nearest(reproject_pixel_stochastic(prev_pos, XorShift32(rng_state)),
+                       gbuf.linear_z,
+                       gbuf.vel_z,
+                       normal,
+                       prev_irr_histlen,
+                       prev_moments,
+                       sum_w);
+    }
 
     if (EXTENDED_SEARCH) {
       if (sum_w <= 0.01) {
diff --git a/src/merian-nodes/nodes/accumulate/accumulate.cpp b/src/merian-nodes/nodes/accumulate/accumulate.cpp
index e5faa63..f3df816 100644
--- a/src/merian-nodes/nodes/accumulate/accumulate.cpp
+++ b/src/merian-nodes/nodes/accumulate/accumulate.cpp
@@ -126,6 +126,7 @@ void Accumulate::process(GraphRun& run,
                          const vk::CommandBuffer& cmd,
                          const DescriptorSetHandle& descriptor_set,
                          [[maybe_unused]] const NodeIO& io) {
+    accumulate_pc.iteration = run.get_total_iteration();
 
     if (accumulate_pc.firefly_filter_enable || accumulate_pc.adaptive_alpha_reduction > 0.0f) {
         MERIAN_PROFILE_SCOPE_GPU(run.get_profiler(), cmd, "compute percentiles");
@@ -183,7 +184,8 @@ Accumulate::NodeStatusFlags Accumulate::properties(Properties& config) {
     accumulate_pc.normal_reject_cos = glm::cos(angle);
     config.config_percent("depth threshold", accumulate_pc.depth_reject_percent,
                           "Reject points with depths farther apart (relative to the max)");
-    needs_rebuild |= config.config_options("filter mode", filter_mode, {"nearest", "linear"});
+    needs_rebuild |= config.config_options("filter mode", filter_mode,
+                                           {"nearest", "bilinear", "stochastic bilinear"});
     needs_rebuild |=
         config.config_bool("extended search", extended_search,
                            "search in a 3x3 radius with weakened rejection thresholds for valid "