floating-point, pixel-wise image division

[user001]

I am trying to port a simple image division function from Matlab to ImageMagick, but the results are different in the two languages. I believe the discrepancy is too great to be simply a consequence of floating point error and instead insinuates some mistake in my ImageMagick syntax.

Assume two grayscale images generated thus:

Code: Select all

magick -size 64x64 gradient:            -blur 1x1 gradient_u.tif
magick -size 64x64 gradient: -rotate 90 -blur 1x1 gradient_v.tif

These images are saved as 16-bit tiff files.

I calculate the ratio as follows, saving the output with double (FP64) precision:

Code: Select all

magick gradient_u.tif -define quantum:format=floating-point -depth 64 \
       gradient_v.tif -define quantum:format=floating-point -depth 64 \
       -fx 'u/v' quotient.tif

The same operation in Matlab would be carried out as follows:

Code: Select all

u = im2double(imread('gradient_u.tif'));
v = im2double(imread('gradient_v.tif'));
q = u ./ v;

However, the two images are unequal

Code: Select all

Q = imread('quotient.tif');
assert(isequal(Q,q)); % assertion fails

The difference is small but larger than I would expect for floating point error:

Code: Select all

d = q-Q;
err = d./q;
max(err(:)); %  1.6550e-07
min(err(:)); % -1.4829e-07

I have seen errors as high as 1E-4 with non-synthetic images, which seems almost certainly outside the realm of floating point error. Have I made some mistake in my ImageMagick implementation? Thanks.

[fmw42]

Using IM 7.0.7.1 Q16 HDRI Mac OS X Sierra, I get the following average over the image divisions two different ways that match. Using "%[fx:mean]" gives the mean in the range 0 to 1. (I have not tried Q32 HDRI)

Code: Select all

magick -size 64x64 gradient:            -blur 1x1 -define quantum:format=floating-point -depth 64 gradient_u.tif
magick -size 64x64 gradient: -rotate 90 -blur 1x1 -define quantum:format=floating-point -depth 64 gradient_v.tif

Code: Select all

magick gradient_u.tif gradient_v.tif -fx 'u/v' -define quantum:format=floating-point -depth 64 -format "%[fx:mean]\n" info:

4.09436

Code: Select all

magick gradient_u.tif gradient_v.tif +swap -define compose:clamp=off -compose divide -composite -define quantum:format=floating-point -depth 64 -format "%[fx:mean]\n" info:

4.09436

[snibgo]

What version IM are you running? In particular, what is the Q number and is it HDRI? Commonly, installations are Q16 integer, so precision is limited to one part in 65536.

EDIT: I cross-posted with Fred. As he implies in his answer, many results of the division will be outside the normal range, ie whiter than white. If integer IM is used, these will be clamped to the normal range.

When asking questions, please always say which version of IM you are using.

[user001]

Good point. I tried the command on two different systems:

One had the following version:
Version: ImageMagick 6.7.2-7 2017-03-22 Q16 http://www.imagemagick.org
Copyright: Copyright (C) 1999-2011 ImageMagick Studio LLC
Features: OpenMP

The second was:
Version: ImageMagick 7.0.4-4 Q16 x86_64 2017-01-20 http://www.imagemagick.org
Copyright: © 1999-2017 ImageMagick Studio LLC
License: http://www.imagemagick.org/script/license.php
Features: Cipher DPC HDRI Modules
Delegates (built-in): bzlib fontconfig freetype jng jpeg ltdl lzma png tiff x xml zlib

Is there a particular version that you would recommend? I see only Q8 and Q16 versions on this page: https://www.imagemagick.org/script/binary-releases.php

[user001]

Note that since the differences between the two quotients are both positive and negative, the mean is fairly similar:

>> format long g
>> mean2(Q(:)) % 4.09287688023606 [the ImageMagick result]
>> mean2(q(:)) % 4.09287688706628 [the Matlab result]

I would be interested to know if the operation could be done with the full precision of IEEE64. Thanks again.

[user001]

It seems that I'll need to make a custom ImageMagick build. I just found a couple pages on this (https://www.imagemagick.org/script/install-source.php, https://www.imagemagick.org/script/adva ... lation.php), which I will read. From a cursory glance, it seems that I will need something like `--enable-hdri` and: `--with-quantum-depth=64` (the second of the linked web pages states that 32-bit is the highest supported version, but I have found others commenting on Q64, so I am assuming the information on that page is dated). Please advise if there are any other options you would recommend for maximum precision.

[user001]

When I try to install (`./configure --prefix=/home/<username>/imq64 --with-quantum-depth=64 --enable-hdri && make && make install`), I get the following warning; this is being done on a CentOS v6.9 system.

Code: Select all

CDPATH="${ZSH_VERSION+.}:" && cd . && /bin/sh /home/<username>/imq64/ImageMagick/config/missing autoconf
aclocal.m4:17: warning: this file was generated for autoconf 2.69.
You have another version of autoconf.  It may work, but is not guaranteed to.
If you have problems, you may need to regenerate the build system entirely.
To do so, use the procedure documented by the package, typically 'autoreconf'.
configure.ac:22: error: Autoconf version 2.69 or higher is required
configure.ac:22: the top level
autom4te: /usr/bin/m4 failed with exit status: 63
WARNING: 'autoconf' is probably too old.
         You should only need it if you modified 'configure.ac',
         or m4 files included by it.
         The 'autoconf' program is part of the GNU Autoconf package:
         <http://www.gnu.org/software/autoconf/>
         It also requires GNU m4 and Perl in order to run:
         <http://www.gnu.org/software/m4/>
         <http://www.perl.org/>
make: *** [configure] Error 63

$ autoconf --version
autoconf (GNU Autoconf) 2.63

I don't have privileges to modify the system on which I am trying to install. Do you have any suggestions?

Perhaps I'll try IM v6 instead of v7.

[magick]

Use the default quantum (Q16) and just add --enable-hdri. Since the pixels are floating point, this configure command-line option is sufficient. Also, quantum 64 has never been officially blessed. It may work but we have not certified it yet.

[user001]

Same problem with IM v6.9.9. Any ideas?

[user001]

Thanks, magick. I did try using the default quantum but I still have the same autoconf error. Would you know how to find the most recent version of IM that is compatible with my version of autoconf? Perhaps I should just try v6.7.2-7, since that is the version already available on the system on which I am trying to install.

Also, I am a bit confused by the difference between hdri (high dynamic range) and quantum depth. You said that Q16 is sufficient for 64-bit calculations as long as HDRI is enabled. Does HDRI simply allow use of floating point pixel as opposed to integer values? If the depth is set to Q16, which I assume means the maximum bit depth for a pixel component is 16, how can one take advantage of 64-bit precision?

In fact, I did try the calculation on a system with Q16 IM with HDRI apparently enabled ("ImageMagick 7.0.4-4 Q16 ... Features: Cipher DPC HDRI Modules"), but the result is still different from what I get in matlab.

Thanks again.

[fmw42]

I believe that computations are done internally in double floating point precision if HDRI is on. The issue is with saving to an image. But with the tiff defines you can save to double precision as I did above. I will defer to the IM developers or others who know the code better to correct me.

[user001]

Thanks for the info, fmw42. It seems that I am still able to write 64-bit tiff output using Q16+HDRI IM:
$ identify quotient.tif
quotient.tif TIFF 64x64 64x64+0+0 64-bit Grayscale Gray 33.1KB 0.000u 0:00.000

[fmw42]

On Q16 HDRI I get:

Code: Select all

magick -size 64x64 gradient:            -blur 1x1 -define quantum:format=floating-point -depth 64 gradient_u.tif
magick -size 64x64 gradient: -rotate 90 -blur 1x1 -define quantum:format=floating-point -depth 64 gradient_v.tif

magick gradient_u.tif gradient_v.tif -fx 'u/v' -define quantum:format=floating-point -depth 64 divide.tif
magick identify -verbose -precision 15 divide.tif
      min: 291.016 (0.00444061)
      max: 1.47581e+07 (225.194)
      mean: 268324 (4.09437)

[fmw42]

See my bug report about fx calculations at viewtopic.php?f=3&t=32689&p=149672#p149672

[user001]

Thanks, what is divide.pfm by the way?