Test failure with vips_scale on Windows x86

[kleisauke]

It looks like the compiler flags -ffast-math -ftree-vectorize (added since 8.7.0) causes issues on Windows 32-bit (x86) on the TestConversion.test_scaleimage testcase.

================================== FAILURES ===================================
_______________________ TestConversion.test_scaleimage ________________________

self = <test_conversion.TestConversion object at 0x04FBA2B0>

    def test_scaleimage(self):
        for fmt in noncomplex_formats:
            test = self.colour.cast(fmt)

            im = test.scaleimage()
            assert im.max() == 255
            assert im.min() == 0

            im = test.scaleimage(log=True)
>           assert im.max() == 255
E           assert 254.0 == 255
E            +  where 254.0 = <function Image.__getattr__.<locals>.call_function at 0x04B6B1E0>()
E            +    where <function Image.__getattr__.<locals>.call_function at 0x04B6B1E0> = <pyvips.Image 100x100 uchar, 3 bands, srgb>.max

test\test-suite\test_conversion.py:723: AssertionError
==================== 1 failed, 92 passed in 28.28 seconds =====================

All noncomplex_formats fails except BandFormat.DOUBLE. This will work fine, for example:

    def test_scaleimage(self):
        test = self.colour.cast(pyvips.BandFormat.DOUBLE)

        im = test.scaleimage()
        assert im.max() == 255
        assert im.min() == 0

        im = test.scaleimage(log=True)
        assert im.max() == 255

These vips-dev-w32-* binaries will currently fail on the TestConversion.test_scaleimage testcase: https://github.com/libvips/build-win64-mxe/releases/tag/v8.7.0 https://github.com/libvips/build-win64-mxe/releases/tag/v8.7.1 https://github.com/libvips/build-win64-mxe/releases/tag/v8.7.2

Note that the 64-bit, v8.7.0-rc3 release and all vips-dev-w32-* binaries from here will work as expected.

/cc @jcupitt

[kleisauke]

@jcupitt Could this be a cause of rounding errors? With this patch it's working correctly (tested on the Windows x86/x64 binaries):

diff --git a/libvips/conversion/scale.c b/libvips/conversion/scale.c
index 1111111..2222222 100644
--- a/libvips/conversion/scale.c
+++ b/libvips/conversion/scale.c
@@ -106,7 +106,19 @@ vips_scale_build( VipsObject *object )
            return( -1 );
    }
    else if( scale->log ) { 
-       double f = 255.0 / log10( 1.0 + pow( mx, scale->exp ) );
+       double f;
+
+       // Try to convert x ** (1/4) into square roots.
+       if( scale->exp == 0.25 ) {
+           // pow(mx, 0.25) --> sqrt( sqrt(mx) )
+           f = 255.0 / log10( 1.0 + sqrt( sqrt(mx) ) );
+       } else {
+           f = 255.0 / log10( 1.0 + pow( mx, scale->exp ) );
+       }
+
+       /* Add .5 to get round-to-nearest.
+        */
+       f += 0.5;

        if( vips_pow_const1( scale->in, &t[2], scale->exp, NULL ) ||
            vips_linear1( t[2], &t[3], 1.0, 1.0, NULL ) ||

I've also tried to optimize it a bit, it's often faster to do the calculations using sqrt than calling the pow function (the exponent for log scale is 0.25 by default).

[kleisauke]

Perhaps the sqrt trick should be part of POW: https://github.com/libvips/libvips/blob/10c4831a709063189d67b7d6076b1b26bf62b2b2/libvips/arithmetic/math2.c#L136-L145

FWIW, GCC also expand pow calls with exponents 0.25, 0.50, 0.75, 1./3., and 1./6. constants into a series of sqrt and cbrt calls under -ffast-math, see: https://github.com/gcc-mirror/gcc/commit/a02df28b537920a1156d483d012f9c47e51c3a9d and https://gcc.godbolt.org/z/ZKFfOZ.

[jcupitt]

Could it just be the +0.5 that's fixing the problem? I'd be surprised if pow() -> sqrt() changed accuracy much. This code is out of the critical path (not in a pixel loop), so speed is less important than clarity here.

Yes, I think pow() is not used much in libvips now because of the performance problems. I agree, that POW() macro could have some more cases.

[kleisauke]

I can confirm that with +0.5 fixes this problem (tested on the Windows x86/x64 binaries):

diff --git a/libvips/conversion/scale.c b/libvips/conversion/scale.c
index 1111111..2222222 100644
--- a/libvips/conversion/scale.c
+++ b/libvips/conversion/scale.c
@@ -106,7 +106,9 @@ vips_scale_build( VipsObject *object )
            return( -1 );
    }
    else if( scale->log ) { 
-       double f = 255.0 / log10( 1.0 + pow( mx, scale->exp ) );
+       /* Add .5 to get round-to-nearest.
+        */
+       double f = 255.0 / log10( 1.0 + pow( mx, scale->exp ) ) + 0.5;

        if( vips_pow_const1( scale->in, &t[2], scale->exp, NULL ) ||
            vips_linear1( t[2], &t[3], 1.0, 1.0, NULL ) ||

This formula with +0.5 is also mentioned in the comment here.

[jcupitt]

Great! Would you like to make this change, or shall I? You have commit rights to libvips, I think.

[kleisauke]

For changes in libvips I prefer to open a pull request rather than commit it directly in the repo itself.

If you don't mind, you can add this change yourself (because I'm not sure if this should be incorporated in the 8.7 branch as well).

[jcupitt]

Oh duh I'm so dumb, should that be:

                        vips_linear1( t[4], &t[5], f, 0.5,

ie. put the 0.5 in the add part of linear1 on the log branch.

I guess this is a bug-fix, so you're right, it should go into 8.7.

[jcupitt]

I pushed +0.5 in log mode to the head of 8.7. Does that look OK?

[kleisauke]

I can confirm that https://github.com/libvips/libvips/commit/9d66420ad5997510389abc8b9d5315fb63867da7 solves this issue (tested on Windows x86 and x64).

I'll close this issue, thanks for fixing this!

[jcupitt]

Great, I'll merge to master as well. And well done for finding this annoying thing!