Fix severe boundary artifacts in Box filter caused by floating-point precision

[DavidDiao]

Description

This PR fixes a severe accuracy bug specifically affecting the Box filter during downsampling.

Users migrating from or comparing outputs with PIL/Pillow might notice that at certain downsampling ratios (e.g., 1000px -> 288px), pica produces periodic visual glitches.

The Root Cause

Unlike continuous filters (Lanczos, Bicubic) which gracefully decay to 0.0 at their window edges, the Box filter is a discrete step function. It relies on a hard cliff: a pixel is either fully included (weight 1.0) or completely excluded (weight 0.0).

When calculating the normalized distance of a boundary pixel to the center diff = (pxl + 0.5) - center:

1. Mathematically, a boundary pixel might have an exact distance of 0.5 or -0.5.
2. However, due to IEEE 754 double-precision floating-point arithmetic in JS, the distance is often computed as something like -0.49999999999999994.
3. The original Box filter logic Math.abs(x) < 0.5 evaluates this as true.

As a result, an out-of-bounds pixel that should have been discarded (weight 0.0) incorrectly receives a full weight of 1.0. This accidentally widens the averaging window by 1 extra pixel for certain periodic target pixels, drastically changing the denominator and introducing a neighboring color, causing the massive deviation.

The Fix

This PR implements a mathematically robust boundary extraction specifically to defend against floating-point jitter at exactly ±0.5 for the Box filter. By offsetting the Math.floor bounds by +0.5 (a standard technique used in PIL's C implementation), we physically exclude the floating-point-corrupted pixels from entering the accumulation loop.

Impact

• Box filter: Completely eliminates the periodic ~15 pixel value errors. Output now flawlessly matches Pillow's Box downsampling at the integer level.
• Other filters (Lanczos, Bicubic, etc.): Unaffected. Since continuous filters evaluate to 0.0 at x = support, floating-point boundary jitter does not impact their sum.

[puzrin]

I understand the problem somehow..., but don't understand your fix.

if (x > -0.5 && x <= 0.5) return 1.0;
      return 0.0;

Looks non-symmetric. TBH, the box filter is shitty, and it's really worth trashing it instead of fixing.

Could you explain the details of the fix and why the problem will not happen at a different scale ratio?

[DavidDiao]

Sorry for the delay and the long post. I wanted to reorganize my thoughts properly, so I used an AI to help structure and translate this explanation for better clarity.

1. The srcLast Index Bug

In the original code, srcLast is calculated as Math.ceil(srcPixel + srcWindow). However, this value actually represents the geometric right edge, not the maximum array index. The correct array index should be the right edge minus one: Math.ceil(...) - 1.
Normally, this off-by-one error goes unnoticed because the fn will return a weight of 0 for this extra pixel, but it means the loop is technically grabbing one pixel outside the intended window. Although it didn't cost much, it was somewhat related to the subsequent modifications.

2. Regarding the "Looks non-symmetric" concern

You are completely right that it is non-symmetric! But actually, in image processing, there is no single "correct" way to do Box/Area resampling—different libraries use different conventions.

Let's look at how three major libraries downsample [10, 20, 30, 40] from length 4 to 3 (Scale factor $S = 4/3 \approx 1.333$, Window $W \approx 0.666$). Here is exactly how they calculate the first destination pixel (index 0), whose geometric range in the source is [0, 1.333]:

• OpenCV (INTER_AREA): Calculates exact overlapping fractional areas.
  • Target 0 covers all of source 0 (value 10, area 1.0) and part of source 1 (value 20, area 0.333).
  • Calculation: (10 * 1.0 + 20 * 0.333) / 1.333 = 12.5. Output: [12.5, 25, 37.5].
• PyTorch (adaptive_avg_pool2d / area): Uses a discrete integer sliding window based on [floor(i * S), ceil((i+1) * S) - 1].
  • Target 0 maps to indices [floor(0), ceil(1.333) - 1] = [0, 1].
  • Calculation: Average of source 0 and 1 -> (10 + 20) / 2 = 15. Output: [15, 25, 35].
    (Fun fact: If we apply the ceil - 1 fix from Point 1, and just change pica's box fn to simply return 1.0, Pica becomes mathematically identical to PyTorch!)
• Pillow (Resampling.BOX): Uses distance to the mapped center. The Pillow docs state: "Each pixel of source image contributes to one pixel of the destination image". To clearly see its asymmetric nature, let's look at downsampling [10, 20, 30, 40, 50] from length 5 to 2 (Scale factor $S = 2.5$, Window $W = 1.25$):
  • Target 0: Its center maps to $0.5 \times 2.5 = 1.25$. The valid distance condition is $-W < (pxl + 0.5) - Center \le W$, which mathematically resolves to $-0.5 < pxl \le 2.0$. The only valid integer indices are 0, 1, 2. Calculation: (10 + 20 + 30) / 3 = 20.
  • Target 1: Its center maps to $1.5 \times 2.5 = 3.75$. The condition resolves to $2.0 < pxl \le 4.5$. The valid integer indices are 3, 4. Calculation: (40 + 50) / 2 = 45.
  • Output: [20, 45]. Notice how Target 0 averages 3 pixels while Target 1 only averages 2 pixels! This is exactly where the asymmetry comes from.

Pica's implementation and naming are inherently closest to Pillow. Therefore, strictly speaking, this PR doesn't fix a "wrong" algorithm, it just perfectly aligns Pica with Pillow's asymmetric convention.

3. Modifying floor/ceil to round

My PR changes the boundaries to:

srcFirst = Math.max(0, Math.floor(srcPixel - srcWindow + 0.5));
srcLast = Math.min(srcSize - 1, Math.floor(srcPixel + srcWindow + 0.5) - 1);

Mathematically, Math.floor(x + 0.5) is identical to Math.round(x). I essentially just changed floor/ceil to round(). (I used the +0.5 syntax in the PR because it's the standard way Pillow writes it in C).

From a geometric perspective, this makes total sense: srcPixel - srcWindow and srcPixel + srcWindow are the exact geometric left and right bounds. By using round(), we are simply snapping these floating-point geometric edges to the nearest discrete pixel indices.

From an algebraic perspective, this is exactly what the x > -win && x <= win condition inside your fn expects. Here is the simple proof:
Condition: -win < (pxl + 0.5) - srcPixel <= win
Isolate pxl: srcPixel - win - 0.5 < pxl <= srcPixel + win - 0.5
Since pxl must be an integer:

• The smallest integer strictly greater than A - 0.5 is floor(A - 0.5) + 1 = floor(A + 0.5) = round(A).
• The largest integer less than or equal to B - 0.5 is floor(B - 0.5) = floor(B + 0.5) - 1 = round(B) - 1.

This proves the valid integer range is mathematically guaranteed to be exactly [round(srcPixel - win), round(srcPixel + win) - 1].

Why pre-calculate this instead of relying on fn? IEEE-754 Precision.
Due to floating-point errors, relying on fn is dangerous for the Box filter. Sometimes an edge pixel evaluates to -0.49999994 instead of -0.5, tricking the fn into returning a full 1.0 weight for an out-of-bounds pixel. By calculating the boundaries upfront with round, we physically exclude these jittery pixels from the loop.

To be completely transparent, because of floating-point math, even round isn't 100% bulletproof. For example, scaling 13 → 6 (scale = 0.461538...):

• When dest=2: srcPixel + srcWindow = 6.499999999999999. (Rounds to 6, so srcLast is 5).
• When dest=3: srcPixel - srcWindow = 6.5. (Rounds to 7, so srcFirst is 7).
  As a result, source pixel 6 is skipped completely! Interestingly, Pillow has this exact same bug.

I wrote a quick AI-assisted script to test this. (This script is the only part I haven't checked myself.) Without this fix, the boundary failure rate is ~30% (and >80% for problematic scales). With this fix, the error rate drops to ~10%. Here is the test script:

Click to expand the test script (JS)

function runFloatingPointTest(maxSrc = 500) {
    let totalPairs = 0;
    let failedPairs = 0;

    const allScales = new Set();
    const failedScales = new Set();
    const failureSamples = [];

    for (let src = 2; src <= maxSrc; src++) {
        for (let dest = 1; dest < src; dest++) {
            // Since i + 1 < dest, we need dest > 1 to have at least one valid 'i'
            if (dest <= 1) continue;

            totalPairs++;
            const scale = dest / src;
            allScales.add(scale);

            let pairFailed = false;
            const scale_inv = 1 / scale;
            const half_over_scale = 0.5 / scale;

            for (let i = 0; i < dest - 1; i++) {
                const lhs = (i + 0.5) * scale_inv + half_over_scale;
                const rhs = (i + 1.5) * scale_inv - half_over_scale;

                const roundedLhs = Math.round(lhs);
                const roundedRhs = Math.round(rhs);

                if (roundedLhs !== roundedRhs) {
                    pairFailed = true;
                    failedScales.add(scale);

                    if (failureSamples.length < 5) {
                        failureSamples.push({
                            src,
                            dest,
                            i,
                            scale,
                            scale_inv,
                            lhs,
                            rhs,
                            roundedLhs,
                            roundedRhs
                        });
                    }
                }
            }

            if (pairFailed) {
                failedPairs++;
            }
        }
    }

    // Output Results
    console.log(`Test Configuration: src in [2, ${maxSrc}]`);
    console.log(`----------------------------------------`);
    console.log(`Total (src, dest) pairs tested : ${totalPairs}`);
    console.log(`Failed (src, dest) pairs       : ${failedPairs} (${((failedPairs / totalPairs) * 100).toFixed(4)}%)`);
    console.log(`----------------------------------------`);
    console.log(`Total unique scales tested     : ${allScales.size}`);
    console.log(`Failed unique scales           : ${failedScales.size} (${((failedScales.size / allScales.size) * 100).toFixed(4)}%)`);
    console.log(`----------------------------------------`);

    if (failureSamples.length > 0) {
        console.log("Failure Samples (showing up to 5):");
        failureSamples.forEach((sample, idx) => {
            console.log(`\nSample ${idx + 1}: src = ${sample.src}, dest = ${sample.dest}, i = ${sample.i}`);
            console.log(`  scale = ${sample.scale}`);
            console.log(`  scale_inv = ${sample.scale_inv}`);
            console.log(`  LHS = ${sample.lhs} (Math.round -> ${sample.roundedLhs})`);
            console.log(`  RHS = ${sample.rhs} (Math.round -> ${sample.roundedRhs})`);
            console.log(`  Difference (LHS - RHS) = ${sample.lhs - sample.rhs}`);
        });
    } else {
        console.log("No mismatches found within the tested range.");
    }
}

function runFloatingPointTest(maxSrc = 500) {
    let totalPairs = 0;
    let failedPairs = 0;

    const allScales = new Set();
    const failedScales = new Set();
    const failureSamples = [];

    for (let src = 2; src <= maxSrc; src++) {
        for (let dest = 1; dest < src; dest++) {
            // Every valid (src, dest) with 0 < dest < src is counted in the denominator
            totalPairs++;
            const scale = dest / src;
            allScales.add(scale);

            const scale_inv = 1 / scale;
            const half_over_scale = 0.5 / scale;

            let pairFailed = false;

            // Loop for i + 1 < dest (will not execute if dest <= 1)
            for (let i = 0; i < dest - 1; i++) {
                const mid = (i + 0.5) * scale_inv;
                const left = Math.floor(mid - half_over_scale);
                const right = Math.ceil(mid + half_over_scale);

                for (let j = left; j <= right; j++) {
                    // Mathematical condition: (2 * j + 1) * dest must be divisible by 2 * src
                    const mathIsHalf = ((2 * j + 1) * dest) % (2 * src) === 0;

                    if (mathIsHalf) {
                        // Actual floating-point calculation
                        const actual = (j + 0.5 - mid) * scale;
                        const actualIsHalf = Math.abs(actual % 1) === 0.5;

                        if (!actualIsHalf) {
                            pairFailed = true;
                            failedScales.add(scale);

                            if (failureSamples.length < 5) {
                                failureSamples.push({
                                    src,
                                    dest,
                                    i,
                                    j,
                                    scale,
                                    scale_inv,
                                    mid,
                                    actual,
                                    actualMod1: Math.abs(actual % 1)
                                });
                            }
                        }
                    }
                }
            }

            // If a pair fails the test, it is added to the numerator
            if (pairFailed) {
                failedPairs++;
            }
        }
    }

    // Output Results (Strictly in English)
    console.log(`Test Configuration: src in [2, ${maxSrc}]`);
    console.log(`----------------------------------------`);
    console.log(`Total (src, dest) pairs tested : ${totalPairs}`);
    console.log(`Failed (src, dest) pairs       : ${failedPairs} (${totalPairs > 0 ? ((failedPairs / totalPairs) * 100).toFixed(4) : 0}%)`);
    console.log(`----------------------------------------`);
    console.log(`Total unique scales tested     : ${allScales.size}`);
    console.log(`Failed unique scales           : ${failedScales.size} (${allScales.size > 0 ? ((failedScales.size / allScales.size) * 100).toFixed(4) : 0}%)`);
    console.log(`----------------------------------------`);

    if (failureSamples.length > 0) {
        console.log("Failure Samples (showing up to 5):");
        failureSamples.forEach((sample, idx) => {
            console.log(`\nSample ${idx + 1}: src = ${sample.src}, dest = ${sample.dest}, i = ${sample.i}, j = ${sample.j}`);
            console.log(`  scale = ${sample.scale}`);
            console.log(`  scale_inv = ${sample.scale_inv}`);
            console.log(`  mid = ${sample.mid}`);
            console.log(`  actual = (j + 0.5 - mid) * scale = ${sample.actual}`);
            console.log(`  |actual % 1| = ${sample.actualMod1} (Expected: 0.5)`);
        });
    } else {
        console.log("All qualified cases matched the mathematical expectation of 0.5.");
    }
}

(Note: This issue only affects the BOX filter because it's a discrete step function. Other continuous filters return 0.0 at the edges anyway, so boundary jitter doesn't matter).

Summary

To sum up, this PR strictly synchronizes Pica's Box filter with Pillow's Box filter. (Actually, Pica is slightly better because Pica preserves 15-bit precision between horizontal/vertical passes, whereas Pillow truncates to 8-bit). If you feel strict Pillow parity isn't necessary for Pica, feel free to close this.

As a side note on why I even dug into this: I have a use case where I need to downsample images to less than 1/3 of their original size. BILINEAR completely skips a lot of source pixels at that scale, while BOX is fast and guarantees all source pixels are averaged.

Let me know what you think!

[puzrin]

If you have no objection, I would keep this open for an uncertain time. Just because the reading is interesting :) . Lack of time to dive into every detail.

I have a use case where I need to downsample images to less than 1/3 of their original size. BILINEAR completely skips a lot of source pixels at that scale, while BOX is fast and guarantees all source pixels are averaged.

I suspect, here can be divergence in terminology

pica/src/mm_resize/resize_filter_info.ts

Line 16 in ce6a3b6

win: 0.5,

box filer (in this project) has 0.5 window (ignore 50% pixels in each direction, 75% total). This works acceptably only when you scale down by > 10x in multiple steps (make a small thumbnail of a huge source) and want the first steps to be much faster.

[DavidDiao]

I have no objection to keeping this PR open.
Additionally, win: 0.5 does not mean ignoring 50% of the pixels in each direction. Since the sampling range is [srcPixel - srcWindow, srcPixel + srcWindow], ±0.5 equates to a total range of exactly 1. You can also refer to the Pillow example I provided earlier: the end position of dest 0 and the start position of dest 1 are both exactly at src 2.

[puzrin]

Ah, you are right, 0.5 in both directions gives 1.0 in total. Shame on me :)