POV-Ray

Analysis shows that this is due to a precision-related issue in a part of some input image handling function (namely map_pos()). More specifically, the program code is based on the erroneous assumption that X+image_width with negative X always evaluates to a value smaller than image_width, ignoring the fact that if X is sufficiently close to 0 it may instead evaluate to image_width due to rounding. This triggers a sanity check leading to the "Picture index out of range" error. Ironically, all the code further downstream is designed in such a way that it is either robust against out-of-range conditions, or must expect such conditions anyway, so the sanity check can safely be removed.

However, code analysis shows that the same wrong assumption appears to have been made further downstream in another function as well (namely no_interpolation()), so the code there needs some careful inspection along the way.

Should be fixed with change #5678 (including the no_interpolation() overhaul).

Confirmed, change #5678 does solves the issue, at least on the original scene.

After installing RC6 I am still getting the same error.
Below is a link to the scene rendered with one thread. You can see that the render fails before the block reaches the image border. For comparison I have included a link to the scene rendered with eight threads.

1 thread
http://imgur.com/8BG0X

8 threads
http://imgur.com/HTldu

Code:

#version 3.7;

global_settings {
adc_bailout 0.0039
ambient_light rgb <1.000,1.000,1.000>
assumed_gamma 1.00
irid_wavelength rgb <0.250,0.180,0.140>
max_trace_level 5
number_of_waves 10
noise_generator 3
charset ascii
}

background { colour rgb <0.000,0.000,0.000> }

#declare Ring_Texture1 =
texture {
uv_mapping
pigment {

image_map{
 png "sat_ring_color.png"
 interpolate 2 
 map_type 0
}
rotate    <90.000,90.000,0.000>

}

finish {

ambient     rgb <0.100,0.100,0.100> *5
brilliance  1.000
crand       0.000
diffuse     0.600
metallic    0.000
phong       0.000
phong_size  40.000
specular    0.000
roughness   0.050

}

}

#declare Camera0 =
camera {
perspective
location <3843.816,38.892,-2660.667>
up y
right 1.333*x
angle 33.000
sky ←0.004,1.000,0.002>
look_at < 0.449, 18.943, 0.102 >
} end Camera0 disc { Disc0
0,y,4.400000,2.500000
texture{ Ring_Texture1 }

  scale     <750.000000,750.000000,750.000000>     

rotate <0.000000,0.000000,-20.000000>
translate < -3500.000000,900.000000,900.000000>
} // end Disc0

camera{ Camera0 }

No surprise about the RC6 behaviour, as it doesn't include change #5678.

I just pulled #5678 compiled fine on my linux, but when I run the sample scene above I get:

povray: support/imageutil.cpp:994: int pov::map_pos(const double*, const TPATTERN*, double*, double*): Assertion `(*xcoor >= 0.0) && (*xcoor < (double)image→iwidth)' failed.

Jim, can you please replace

    assert ((*xcoor >= 0.0) && (*xcoor < (DBL)image->iwidth));
    assert ((*ycoor >= 0.0) && (*ycoor < (DBL)image->iheight));

with

    assert (*xcoor >= 0.0);
    assert (*xcoor < (DBL)image->iwidth);
    assert (*ycoor >= 0.0);
    assert (*ycoor < (DBL)image->iheight);

in pov::map_pos() to get a clearer picture of what exactly goes wrong there?

I'd also like to invite anyone to review the code of the warp() function (see types.h), whether there's any possible source of error that I might have missed. At present I can only imagine that fmod() on linux fails to adhere to the C++ standard, but I definitely don't like that thought.

Same error as before:
povray: support/imageutil.cpp:997: int pov::map_pos(const double*, const TPATTERN*, double*, double*): Assertion `*xcoor < (double)image→iwidth' failed.

which corosponds to the following:
assert (*xcoor < (DBL)image→iwidth);

Did some research and I think I've found something that may apply:
glibc bug

I don't think any of the two known glibc fmod() bugs are relevant here. One applies to very small values of both dividend and divisor (we know our divisor is 1 or greater) and causes a NaN result (which can't be the case here because the first assert would fail already), the other applies only to ia64 (Itanium architecture).

I guess there is yet another fmod() bug, and we'll need to somehow work around it; Jim, can you please try again with change #5684 applied?

Nope ... not yet. Basically the same failure message.
I made sure I'd removed the changes you had me make in source/backend/support/imageutil.cpp. I just removed file from my workspace and did sync to head in perforce. I did NOT re-run configure as this was a changed file rather than a new one. As always I did a make clean before starting the re-build.

Jim, I'd like you to add a few more assert statements:

	// verify that the result is not NaN
	assert ((*xcoor < 0.0) || (*xcoor >= 0.0));
	assert ((*ycoor < 0.0) || (*ycoor >= 0.0));

	// verify that the result is finite
	assert (abs(*xcoor) <= numeric_limits<DBL>::max());
	assert (abs(*ycoor) <= numeric_limits<DBL>::max());

	// verify that the result complies with the lower bound
	assert (*xcoor >= 0.0);
	assert (*ycoor >= 0.0);

	// verify that the result complies with the upper bound
	assert (*xcoor <= (DBL)image->iwidth);
	assert (*ycoor <= (DBL)image->iheight);

	// verify that the result strictly complies with the upper bound
	assert (*xcoor < (DBL)image->iwidth);
	assert (*ycoor < (DBL)image->iheight);

OK narrowed it down a bit: Assertion `*xcoor < (double)image→iwidth' failed ... points to the next to the last assert in the code above.

This should bloody freaking be impossible. Jim, can you please make absolutely sure what you're running does include types.h revision #37?

Can you also please add the following assertions to the end of wrap() (types.h) right before the return statement:

	assert (!(tempVal >= upperLimit));
	assert (tempVal < upperLimit);

... on the proper version of types.h, I was pretty sure before, so just to make double sure I first did make distclean, prebuild.sh, configure, make ... etc to wipe the slate clean.

OK ... now I'm confused, shouldn't the addition of the asserts in wrap (types.h) should have caused it to loose it's lunch before returning to the caller in map_pos (imageutil.cpp)

Believe it or not ... here's the error:
povray: support/imageutil.cpp:994: int pov::map_pos(const double*, const TPATTERN*, double*, double*): Assertion `(*xcoor >= 0.0) && (*xcoor < (double)image→iwidth)' failed.

Yeah, we definitely have something weird going on there.

Does your build happen to be 32 bit, with x87 FPU mode (as opposed to SSE2)?

In that case there's still some potentially sane explanation left: The x87 internally uses 80 bits for computations, rather than 64 bits; if the 64-bit tempVal variable is optimized away and the interim result just stored in the 80-bit FPU register, we may have a value smaller than upperLimit that, when later stored in the 64-bit variable *xcoor, is rounded to upperLimit.

(64-bit builds should always use SSE2, which uses 64-bit precision internally.)

Jim, I've got another change for you to test: #5685. If it's indeed an "excessive precision" issue, I'd expect this to solve it.

My system is 32bit with x87 FPU ... I was in the process of looking for a configure option that might force one mode or another but found none, when I saw the email about #5685. Good news we're fat and happy now! I had a look at what you did in types.h. More than one way to skin a cat ... eh! Thanks for your persistence of vision ... and the clever fix!

For what it's worth, fmod is used in several files:

source/backend/support/imageutil.cpp
source/backend/render/tracetask.cpp
source/backend/pattern/pattern.cpp
source/backend/parser/express.cpp
source/backend/vm/fnpovfpu.cpp
source/backend/vm/fnintern.cpp

libraries/ilmbase/Imath/ImathMath.h
libraries/ilmbase/Imath/ImathEuler.h
source/base/types.h

I would suggest using

if((tempVal + EPSILON) >= upperLimit)

in line 139 of types.h.

Thorsten, I disagree; that would be a kludge introducing an artificial limit, rather than using the inherent limits of the double-precision floating point format.

I think POV-Ray already has more of such artificial limits than is good - most of them apparently poorly understood and quite arbitrary rather than based on some particular rationale; the errors/to-dos regarding numerical precision discrepancies between 3.6 and 3.7 are telling a story there.

Remember that we are talking about floating-point arithmetic here. You need an epsilon because of the accumulating error over many computations and the acceted limited precision of many compley mathematical operations (starting with division and multiplication, continuing to estimate square roots, etc). So there always is an epsilon to consider. This is one of those cases.

As for your other remark, these differences are the result of identifying all the file-local epsilons and small tolerances, which were frequently inhertited from theird party patches, many of which we integrated in 3.5. There was never motivation by anybody to clean those up, as it its tedious and involves testing on more than one platform, and ideally more than one architecture to find all the nasty assumtions some of that code makes. Many tolerances there are also way too big.

Apart from that we have various places where we are not tolerant and smart enough. This here is one such place!

Absolutely not, Thorsten; An epsilon does not belong in such a basic function - if any code further downstream processes the result in a way that would lead to additional precision issues, it needs to handle them of its own, in a manner that fits the issues to be expected there. Hard-coding such handling in wrap() would deny us any flexibility there.

Same in case we should have any precision issues upstream; but from the nature of the wrap() function I would not expect any such issues.

The function is well-defined with regards to precision limitations - that's all we need.

Pardon, but that is pure nonsense: The very reality of this "bug" already shows that your appraoch is doomed. Clipping is one of the core operations that exhibit precision-related problems. You cannot expect code further down the line to handle bugs in your clipping function, which clearly outputs out of range values. After all, the clipping function would be completely obsolete if code had to check afterwards if it actually performed the clipping correctly. The only possible solution for such code to deal with the detected problem then would be to clip the value itself. So what is the point of the clipping function then???

Thorsten, the very reality of this bug shows the necessity of such a function. The whole point of it is to provide a single point of implementation of this clipping (or rather, wrapping) operation, in order to once and for all fix any precision-related issues inherent to implementations of this particular operation.

The fact that it was a comparatively long way to identify and sort out all the precision-related issues does not negate the fact that a bug-free implementation (without an arbitrary epsilon) is possible and does make sense. Note that Jim has already confirmed by now that the current implementation does not output out of range values anymore.

Also remember that bugs are not fixed by randomly throwing epsilons (or other kludgy patches) at them. Bugs are fixed by identifying and understanding the underlying problem, and then address that problem. That has been done by now.

(In case you feel inclined to contest that the epsilon-based approach is kludgy, be reminded that epsilons are only good so long as the orders of magnitude involved are sufficiently "small"; beyond a certain limit the approach breaks down as well, which is not helpful at all if the code downstream relies on a robust implementation.)