Accuracy problems with -evaluate, -fx, -function

[alosca]

Calling

convert mock.pgm -depth 16 -evaluate pow 2 mock-square.pgm

to compute the square value of each pixel in the input image mock.pgm is producing odd results, i.e., the resulting numbers are not actually the square values and they are off by quite a margin of error.

I made sure square values of the input pixels in image mock.pgm would not exceed 16 bits depth - indeed, the max value in mock.pgm is 89, see below. I also tried -fx "u*u" and -function polynomial 1,0,0 but to no avail (they probably use the same floating point engine), the resulting numbers have the same error. If I manually change the bit depth in the input file from 8 (255) to 16 (65535) and keep pixel values the same results are worse: square values are all 0's!

A few values that were reported that are not correct:

Input - Output - Correct
9 - 82 - 81
10 - 101 - 100
89 - 7983 - 7921
84 - 7111 - 7056

Any idea of what I am missing or this is a known bug? I moved '-depth 16' to different places in the command line but that also didn't help.

Ultimately, what I want to do is to compute the sum of the squared difference between two images. Is there a built in routine already in Magick for that?

Below are the input and output .pgm files in plain text format (they are only 9x9 images) in case someone wants to try with the same data.

Thanks,
- Alex

--- mock.pgm

P2
9 9
255
0 12 6 6 3 10 36 89 64
3 3 1 4 9 12 58 78 61
20 3 3 14 11 19 74 73 27
6 9 9 7 15 40 82 84 69
1 18 6 0 10 53 81 59 47
5 6 6 33 19 64 89 83 35
4 10 11 29 46 89 83 52 35
12 5 4 41 80 80 50 41 13
1 6 14 40 67 71 61 31 9

-- mock-square.pgm

P2
9 9
65535
0 145 36 36 9 101 1306 7983 4128
9 9 1 16 82 145 3390 6132 3750
403 9 9 198 122 364 5519 5371 735
36 82 82 49 227 1613 6777 7111 4798
1 327 36 0 101 2831 6612 3508 2226
25 36 36 1098 364 4128 7983 6943 1235
16 101 122 848 2133 7983 6943 2725 1235
145 25 16 1694 6450 6450 2520 1694 170
1 36 198 1613 4524 5081 3750 969 82

[fmw42]

I get similar results in IM 6.7.0.10 Q16 Mac OSX Tiger from the following commands:

convert test.pgm -evaluate pow 2 -depth 16 txt:
convert test.pgm -gamma 0.5 -depth 16 txt:
convert test.pgm fx "u*u" -depth 16 txt:
convert test.pgm test.pgm -compose multiply -composite -depth 16 txt:


For example the two items in red come from (12^2 and 10^2)

# ImageMagick pixel enumeration: 9,9,65535,gray
0,0: ( 0, 0, 0) #000000000000 black
1,0: ( 145, 145, 145) #009100910091 gray(0.221256%,0.221256%,0.221256%)
2,0: ( 36, 36, 36) #002400240024 gray(0.0549325%,0.0549325%,0.0549325%)
3,0: ( 36, 36, 36) #002400240024 gray(0.0549325%,0.0549325%,0.0549325%)
4,0: ( 9, 9, 9) #000900090009 gray(0.0137331%,0.0137331%,0.0137331%)
5,0: ( 101, 101, 101) #006500650065 gray(0.154116%,0.154116%,0.154116%)
6,0: ( 1306, 1306, 1306) #051A051A051A gray(1.99283%,1.99283%,1.99283%)
7,0: ( 7983, 7983, 7983) #1F2F1F2F1F2F gray(12.1813%,12.1813%,12.1813%)
8,0: ( 4128, 4128, 4128) #102010201020 gray(6.29892%,6.29892%,6.29892%)
0,1: ( 9, 9, 9) #000900090009 gray(0.0137331%,0.0137331%,0.0137331%)
1,1: ( 9, 9, 9) #000900090009 gray(0.0137331%,0.0137331%,0.0137331%)
2,1: ( 1, 1, 1) #000100010001 gray(0.0015259%,0.0015259%,0.0015259%)
3,1: ( 16, 16, 16) #001000100010 gray(0.0244144%,0.0244144%,0.0244144%)
4,1: ( 82, 82, 82) #005200520052 gray(0.125124%,0.125124%,0.125124%)
5,1: ( 145, 145, 145) #009100910091 gray(0.221256%,0.221256%,0.221256%)
6,1: ( 3390, 3390, 3390) #0D3E0D3E0D3E gray(5.17281%,5.17281%,5.17281%)
7,1: ( 6132, 6132, 6132) #17F417F417F4 gray(9.35683%,9.35683%,9.35683%)
8,1: ( 3750, 3750, 3750) #0EA60EA60EA6 gray(5.72213%,5.72213%,5.72213%)
0,2: ( 403, 403, 403) #019301930193 gray(0.614939%,0.614939%,0.614939%)
1,2: ( 9, 9, 9) #000900090009 gray(0.0137331%,0.0137331%,0.0137331%)
2,2: ( 9, 9, 9) #000900090009 gray(0.0137331%,0.0137331%,0.0137331%)
3,2: ( 198, 198, 198) #00C600C600C6 gray(0.302129%,0.302129%,0.302129%)
4,2: ( 122, 122, 122) #007A007A007A gray(0.18616%,0.18616%,0.18616%)
5,2: ( 364, 364, 364) #016C016C016C gray(0.555428%,0.555428%,0.555428%)
6,2: ( 5519, 5519, 5519) #158F158F158F gray(8.42145%,8.42145%,8.42145%)
7,2: ( 5371, 5371, 5371) #14FB14FB14FB gray(8.19562%,8.19562%,8.19562%)
8,2: ( 735, 735, 735) #02DF02DF02DF gray(1.12154%,1.12154%,1.12154%)
0,3: ( 36, 36, 36) #002400240024 gray(0.0549325%,0.0549325%,0.0549325%)
1,3: ( 82, 82, 82) #005200520052 gray(0.125124%,0.125124%,0.125124%)
2,3: ( 82, 82, 82) #005200520052 gray(0.125124%,0.125124%,0.125124%)
3,3: ( 49, 49, 49) #003100310031 gray(0.0747692%,0.0747692%,0.0747692%)
4,3: ( 227, 227, 227) #00E300E300E3 gray(0.34638%,0.34638%,0.34638%)
5,3: ( 1613, 1613, 1613) #064D064D064D gray(2.46128%,2.46128%,2.46128%)
6,3: ( 6777, 6777, 6777) #1A791A791A79 gray(10.341%,10.341%,10.341%)
7,3: ( 7111, 7111, 7111) #1BC71BC71BC7 gray(10.8507%,10.8507%,10.8507%)
8,3: ( 4798, 4798, 4798) #12BE12BE12BE gray(7.32128%,7.32128%,7.32128%)
0,4: ( 1, 1, 1) #000100010001 gray(0.0015259%,0.0015259%,0.0015259%)
1,4: ( 327, 327, 327) #014701470147 gray(0.49897%,0.49897%,0.49897%)
2,4: ( 36, 36, 36) #002400240024 gray(0.0549325%,0.0549325%,0.0549325%)
3,4: ( 0, 0, 0) #000000000000 black
4,4: ( 101, 101, 101) #006500650065 gray(0.154116%,0.154116%,0.154116%)
5,4: ( 2831, 2831, 2831) #0B0F0B0F0B0F gray(4.31983%,4.31983%,4.31983%)
6,4: ( 6612, 6612, 6612) #19D419D419D4 gray(10.0893%,10.0893%,10.0893%)
7,4: ( 3508, 3508, 3508) #0DB40DB40DB4 gray(5.35286%,5.35286%,5.35286%)
8,4: ( 2226, 2226, 2226) #08B208B208B2 gray(3.39666%,3.39666%,3.39666%)
0,5: ( 25, 25, 25) #001900190019 gray(0.0381476%,0.0381476%,0.0381476%)
1,5: ( 36, 36, 36) #002400240024 gray(0.0549325%,0.0549325%,0.0549325%)
2,5: ( 36, 36, 36) #002400240024 gray(0.0549325%,0.0549325%,0.0549325%)
3,5: ( 1098, 1098, 1098) #044A044A044A gray(1.67544%,1.67544%,1.67544%)
4,5: ( 364, 364, 364) #016C016C016C gray(0.555428%,0.555428%,0.555428%)
5,5: ( 4128, 4128, 4128) #102010201020 gray(6.29892%,6.29892%,6.29892%)
6,5: ( 7983, 7983, 7983) #1F2F1F2F1F2F gray(12.1813%,12.1813%,12.1813%)
7,5: ( 6943, 6943, 6943) #1B1F1B1F1B1F gray(10.5943%,10.5943%,10.5943%)
8,5: ( 1235, 1235, 1235) #04D304D304D3 gray(1.88449%,1.88449%,1.88449%)
0,6: ( 16, 16, 16) #001000100010 gray(0.0244144%,0.0244144%,0.0244144%)
1,6: ( 101, 101, 101) #006500650065 gray(0.154116%,0.154116%,0.154116%)
2,6: ( 122, 122, 122) #007A007A007A gray(0.18616%,0.18616%,0.18616%)
3,6: ( 848, 848, 848) #035003500350 gray(1.29397%,1.29397%,1.29397%)
4,6: ( 2133, 2133, 2133) #085508550855 gray(3.25475%,3.25475%,3.25475%)
5,6: ( 7983, 7983, 7983) #1F2F1F2F1F2F gray(12.1813%,12.1813%,12.1813%)
6,6: ( 6943, 6943, 6943) #1B1F1B1F1B1F gray(10.5943%,10.5943%,10.5943%)
7,6: ( 2725, 2725, 2725) #0AA50AA50AA5 gray(4.15808%,4.15808%,4.15808%)
8,6: ( 1235, 1235, 1235) #04D304D304D3 gray(1.88449%,1.88449%,1.88449%)
0,7: ( 145, 145, 145) #009100910091 gray(0.221256%,0.221256%,0.221256%)
1,7: ( 25, 25, 25) #001900190019 gray(0.0381476%,0.0381476%,0.0381476%)
2,7: ( 16, 16, 16) #001000100010 gray(0.0244144%,0.0244144%,0.0244144%)
3,7: ( 1694, 1694, 1694) #069E069E069E gray(2.58488%,2.58488%,2.58488%)
4,7: ( 6450, 6450, 6450) #193219321932 gray(9.84207%,9.84207%,9.84207%)
5,7: ( 6450, 6450, 6450) #193219321932 gray(9.84207%,9.84207%,9.84207%)
6,7: ( 2520, 2520, 2520) #09D809D809D8 gray(3.84527%,3.84527%,3.84527%)
7,7: ( 1694, 1694, 1694) #069E069E069E gray(2.58488%,2.58488%,2.58488%)
8,7: ( 170, 170, 170) #00AA00AA00AA gray(0.259403%,0.259403%,0.259403%)
0,8: ( 1, 1, 1) #000100010001 gray(0.0015259%,0.0015259%,0.0015259%)
1,8: ( 36, 36, 36) #002400240024 gray(0.0549325%,0.0549325%,0.0549325%)
2,8: ( 198, 198, 198) #00C600C600C6 gray(0.302129%,0.302129%,0.302129%)
3,8: ( 1613, 1613, 1613) #064D064D064D gray(2.46128%,2.46128%,2.46128%)
4,8: ( 4524, 4524, 4524) #11AC11AC11AC gray(6.90318%,6.90318%,6.90318%)
5,8: ( 5081, 5081, 5081) #13D913D913D9 gray(7.75311%,7.75311%,7.75311%)
6,8: ( 3750, 3750, 3750) #0EA60EA60EA6 gray(5.72213%,5.72213%,5.72213%)
7,8: ( 969, 969, 969) #03C903C903C9 gray(1.4786%,1.4786%,1.4786%)
8,8: ( 82, 82, 82) #005200520052 gray(0.125124%,0.125124%,0.125124%)

[alosca]

By the way, for completeness, the platform I ran the tests above is Linux Ubuntu 11.04 natty and IMagick 6.6.3-8 (6.6.2-6 2011-03-16 gives the same problem)

Linux XXX 2.6.38-8-generic #42-Ubuntu SMP Mon Apr 11 03:31:24 UTC 2011 x86_64 x86_64 x86_64 GNU/Linux

Version: ImageMagick 6.6.3-8 2010-08-21 Q16 http://www.imagemagick.org
Copyright: Copyright (C) 1999-2010 ImageMagick Studio LLC
Features: OpenMP


- Alex.

[fmw42]

Using normalized values

(10/255)*(10/255)*65535=100.7843137025 which rounds to 101.

This appears to be the basic issue of using normalized values rather than true floating point computations (which can be done in IM only using HDRI), if I am not mistaken. That is the trade one has to make by using different Q level compiles for efficiency.

[magick]

The -fx and -evaluate options works on normalized pixel values [0..1.0] and then multiplies by QuantumRange (65535 for Q16). For Q16 the results is 81.63280284419488. For Q8 we get 0.317647. If you have a better way that is invariant to the quantum-depth (Q8, Q16, Q32, Q64), let us know. Keep in mind that we multiply by 257 to upscale Q8 to Q16 (9*257 = 2313).

[alosca]

I don't know the internals of IM but wouldn't using double and, if needed, quad precision floating point numbers to perform the operations be sufficient? Only after the operations are done results would be mapped to whatever quantization is necessary to hold the numbers. No need for normalization here.

If double precision floating numbers are used internally to hold pixel values then 'most' operations would be safe. Of course one cannot completely avoid under/overflow unless you use arbitrary precision numbers but this comes at a cost. In scientific computing this is the standard to avoid round off errors as much as possible.

Using normalized values in the multiplication will necessarily incur in a error, approximately

(P/255)*(P/255)*65,535 = P^2 * (65535/(255*255)) = P^2 * 1.00784313725490196078

That's why

9*9 = 81 * 1.00784313725490196078 = 81.63529411764705882318 ~= 82
89*89 = 7921 * 1.00784313725490196078 = 7983.12549019607843133838 ~= 7983

That's definitely bad!

I am probably also misunderstanding '-depth' since what I wanted in my initial post was to guarantee that I have enough bytes to store the result of the multiplied numbers and not actually linearly map the result to the [0, 65535] range. In other words, I want 10*10 to always give me 100 independent of the quantization - sure, it has to be large enough to hold that number. Maybe there is another operator that I can use?

[anthony]

NOTE that to most mathematical functions, a puixel value of '9' is a numberical value of 9/(2^depth-1)
that color values are really a numerical value from 0.0 to 1.0.

If the value is 'stored' into the image, the value is rounded into an integer, UNLESS you use HDRI.
However even in HDRI, the values are stored with a fixed multiple of the QUANTURM RANGE (EG: 2^depth-1), so that routines that normally are only designed to deal with integers, still work as expected.

If you are directly using say -fx expressions to generate a string (percent escape) then the value will be processed as a double and converted to a ASCII string directly from that doulbe without 'storing' into an integer and the associated rounding. EG:

Code: Select all

   convert null: -format '%[fx:9.0*9.0]' info:

will print '81.0'

BUT a value of 9 in an image is NOT is numerical value whcih depends on the image 'depth' when reading in, and out. On the other hand it in-memory its 'depth' depends on the IM's compile time Quality, and if HDRI is enabled.

In summery IM preserves the values, but always treats them as a 'fixed point integer in a 0.0 to 1.0 numerical value range, which corresponds to the appropriate color range (black to white, transparent to opaque, no-red to full-red, etc).

[alosca]

Unfortunately it is the [0.0, 1.0] representation of color/intensity values that is troublesome. It is a smart decision to adopt a canonical representation, to have a uniform way to seamlessly handle all possible different quantizations, but it seems that has its drawbacks as it was shown by the multiplication example. In my view, it would be hard to convince someone that it is OK to return something different of 81 for 9*9, or 7921 for 89*89, and so on.

9.pgm is an image with all pixel values equal to 9. The following

convert 9.pgm -evaluate pow 2

results in an image with 0's for all pixel values. It just doesn't feel right! Something, I don't know what, would need to be changed in IM to fix the problem. Saying that a pixel of value 9 is not actually 9, it is 0.0352941176471, and thus the 0's one obtains as a result of pow(.,2) because 0.0352941176471 * 0.0352941176471 * 255 = 0.317647058824 which is 0 when truncated, is not convincing.

Maybe there is a way to still use [0,1] and make the calculations be correct ?!

"convert 9.pgm -evaluate pow 0.5" returns 48, not 3, for pixel values. That's not what one would expect from the 'pow' function.

[fmw42]

I think you have to account for the scaling between 0-1 and quantum range (twice - once manually and once because of the scaling that IM does automatically on output) in your convert command.

That is IM is multiplying 9/255 * 9/255 then multiply by 255 on output. So you need to provide another multiply by 255

convert -size 1x1 xc:"gray(9)" -depth 8 txt:
# ImageMagick pixel enumeration: 1,1,255,rgb
0,0: ( 9, 9, 9) #090909 rgb(9,9,9)



convert -size 1x1 xc:"gray(9)" -evaluate pow 2 -depth 8 txt:
# ImageMagick pixel enumeration: 1,1,255,rgb
0,0: ( 0, 0, 0) #000000 black

IM 6.7.0.10 Q16 Mac OSX Tiger.


whereas:

convert -size 1x1 xc:"gray(9)" -evaluate pow 2 -evaluate multiply 255 -depth 8 txt:
# ImageMagick pixel enumeration: 1,1,255,rgb
0,0: ( 81, 81, 81) #515151 rgb(81,81,81)

convert -size 1x1 xc:"gray(3)" -evaluate pow 2 -evaluate multiply 255 -depth 8 txt:
# ImageMagick pixel enumeration: 1,1,255,rgb
0,0: ( 9, 9, 9) #090909 rgb(9,9,9)

convert -size 1x1 xc:"gray(15)" -evaluate pow 2 -evaluate multiply 255 -depth 8 txt:
# ImageMagick pixel enumeration: 1,1,255,rgb
0,0: (225,225,225) #E1E1E1 rgb(225,225,225)

convert -size 1x1 xc:"gray(15)" -evaluate pow 2 -evaluate multiply 0.99221789883 -depth 16 txt:
# ImageMagick pixel enumeration: 1,1,65535,rgb
0,0: ( 225, 225, 225) #00E100E100E1 rgb(0.343328%,0.343328%,0.343328%)

convert -size 1x1 xc:"gray(25)" -evaluate pow 2 -evaluate multiply 0.99221789883 -depth 16 txt:
# ImageMagick pixel enumeration: 1,1,65535,rgb
0,0: ( 625, 625, 625) #027102710271 rgb(0.953689%,0.953689%,0.953689%)

where .00221789883 = 255*255/65535



For -fx, you need to do the scaling twice so that fx sees 9x9.


convert -size 1x1 xc:"gray(9)" -depth 8 -format "%[fx:u*u]" info:
0.00124567

whereas

convert -size 1x1 xc:"gray(9)" -depth 8 -format "%[fx:255*255*u*u]" info:
81

convert xc: -format "%[fx:9*9]" info:
81

[alosca]

Thanks for the tip. Unfortunately this is a cumbersome solution and probably not right when using high order polynomials (my primary motivation when I started this thread) specially when the polynomial order is only known at run time. Also, for numbers larger than 15 all gets mapped either to 255 (sure, we need more than 8 bits to store the result)

convert -size 1x1 xc:"gray(16)" -evaluate pow 2 -evaluate multiply 255 -depth 8 txt:
# ImageMagick pixel enumeration: 1,1,255,rgb
0,0: (255,255,255) #FFFFFF white

or, when using -depth 16, to something not exactly the right answer.

I still don't know what would be the right solution. Is there a way to ask IM to map values back to the original quantization *before* they are used e.g. with fx ? Something like a '-quant' operator:

convert image-8bits.pgm -quant -format "%[fx:u*u]"

or even

convert image-8bits.pgm -quant -evaluate pow 2 -depth 16 output.pgm

which gets pixel values mapped back to their original, non normalized values before evaluating pow? This would create a 16 bits image output.pgm with the right values. Or even '-quant Q' which would tell IM that we need as many bits as in Q necessary to store the processed values. Just a thought!

[fmw42]

See my edits above for -depth 16. Yes, it is cumbersome and this solution is only for doing a square, though the same principle can be used for any computation, though, it is not simple.

I know of no way for IM to handle this invisibly, nor any other command to help.

What IM needs is a true floating point configuration -- i.e. one that does not do the scaling from image values to the range 0-1. It should process as double precision (64-bits) so that it can handle any format up to simple 32-bit float. I think that is probably the only way that might work. The scaling to the range 0 to 1 is ingrained into IM v6. So something like this would have to be considered for IM 7.

But I am not an expert on IM configurations and don't know whether this kind of thing is even practical to implement if IM 7 keeps the same philosophy of scaling to 0-1 so that all data types are processed the same way. The current philosophy also provides a means of speeding up IM when you use Q8 rather than Q16 or Q32. But I would give up that speed if there were an option other than the current HDRI to avoid the scaling to the range 0-1.

This is something that the IM developers would need to consider now if they feel that it would be worthwile to implement such an option under IM 7. But as I said, it is a dramatic departure from the current philosphy, which I why I only suggest it as an option for consideration.

P.S. In the mean time, you might experiment with ImageJ


_____________________

Another odd example is consider a gradient that spans linearly from black to white. If one does

convert -size 100x100 gradient: -fx "log(u)" result

then one finds that the result is more linear than logarithmic as log(u) is done by adding 1 to u and then taking the log. So if u ranges from 0 to 1, we get log(1) to log(2) which does not show much curvature. That is why there is a scaled log function under -evaluate log scalevalue, so one can give it a larger scaling and get more curvature in the result.

_____________________

[alosca]

Let's hope for the best! Hopefully IM 7 will adopt the double precision representation, fix these problems, and make it more robust. Right now it makes me wonder if the errors in built in operations like convolution and such are not too detrimental and results reliable.

Thanks for the ImageJ tip. Besides Matlab, Netpbm has been handy and awk also helps to do math with images.

[anthony]

For general image processing it its far more common to regard color values as representing a 0 to 1 range. That will not change. even in IM v7. It is the way it should work!

The problem is when you want to deal with the values in the native quantum range rather than a 'normalized' range.

But then what is the quantum range you want to process in. The compiled quantum range, the depth the image was read or will right in, or a fixed depth of 8, or 16 bits, Perhaps even a percentage range (0 to100)?

All those ranges are equally valid, and completely depends on what you the user are trying to do! So what should IM do, that is fair to everyone, and not just the needs of a specific situation.



Note about FX espressions... log(x) should be the normal double Math library function, IM should not be adding a value of 1.0 to the number before it is processed. Just as atan(x) atan2(y,x) return values in +/- PI range, or that sin(), cos(), tan() request values in radians.

Yes it can make it awkward, but it mean in general the FX espression is not dependant on IM compile time quality, image depth, or other differences. The same goes for composition mathematical operators.

The -evaluate and -function operators however are meant to make them easy to use by users, and thus typically provide both input and output scaling and offset parameters.

[fmw42]

Note about FX espressions... log(x) should be the normal double Math library function, IM should not be adding a value of 1.0 to the number before it is processed. Just as atan(x) atan2(y,x) return values in +/- PI range, or that sin(), cos(), tan() request values in radians.

I was taking about image processing with -fx, not fx expressions.



This does not seem right to me:


convert -size 1x1 xc:black -fx "log(u)" txt:
# ImageMagick pixel enumeration: 1,1,65535,rgb
0,0: ( 0, 0, 0) #000000000000 black

convert -size 1x1 xc:white -fx "log(u)" txt:
# ImageMagick pixel enumeration: 1,1,65535,rgb
0,0: ( 0, 0, 0) #000000000000 black

Similarly with ln().

[anthony]

Considering that the log function is negative for any number in the 0 to 1 range, I am not suprised you have gotten zero! (clipped). Even a white value is 1.0, the log of which is also zero!

The %[fx:..] and -fx functions are exactly the same. The only deference is what values are being processed, and how many times it gets processed.

-fx is called once for each value in the first image only. %[fx:..] is called once for each image in the sequence.

Code: Select all

convert xc: -format '%[fx:log(2)]' info:
0.30103
convert xc: -format '%[fx:log(1)]' info:
0
convert xc: -format '%[fx:log(0.5)]' info:
-0.30103
convert xc: -format '%[fx:log(10)]' info:
1

which shows that it is base 10 without any offset, as it should be.

So and image value u,s,r,b,g,a, etc which are all 0-1 range will produce a negative, that is outside the 0-1 range so will be clipped if stored in a non-HDRI integer image.

A value of log(0) is of course -inf, and normal floating point will work with infinite values in a logical way, though if stored into an image integer is still a negate value so gets clipped to zero. A HDRI image file however will probably continue to store the infinite value.

Applying a function to a linear gradient and generating a profile graph will let you see what is happening more clearly.