Mini ray tracer
Jon Harrop
I recently tried my hand at writing a ray tracer (something I haven't done
for 15 years!) and, in particular, I to write the shortest possible ray
tracer which was still comprehensible and interesting. The initial result
was a 222-line OCaml program which incrementally renders a sphere-flake via
OpenGL:
http://www.ffconsultancy.com/free/ray_tracer/index.html
I then cut this program down to a 66-line OCaml program and ported it into a
97-line C++ program. These programs are compared on this page:
http://www.ffconsultancy.com/free/ray_tracer/comparison.html
The cut-down versions don't do reflections, color or such a pretty
sphere-flake and output a greyscale PGM but they do still use hierarchical
spherical bounding volumes to accelerate rendering enough that they can
render scenes containing millions of spheres.
--
Dr Jon D Harrop, Flying Frog Consultancy
http://www.ffconsultancy.com
asama...@matf.bg.ac.yu
compiler? I got running time of about 36s for scene level 3 on
Pentium3 machine and running time of about 25s when all references to
"double" type in C++ code changed to "float" (and float is usually
"good enough" for ray tracing)...
Regards,
Alex
tbp...@gmail.com
> I then cut this program down to a 66-line OCaml program and ported it
into a
> 97-line C++ program. These programs are compared on this page:
>
> http://www.ffconsultancy.com/free/ray_tracer/comparison.html
Somewhere on an opteron 146 running debian64:
# g++ -g -O3 ray.cpp -o ray
# time ./ray >tt
real 0m13.959s
user 0m13.946s
sys 0m0.009s
# objdump -x ray
[2mn later]
# diff ray_orig.cpp ray.cpp
12,16c12,18
< Vec operator+(Vec a, Vec b) { return Vec(a.x + b.x, a.y + b.y, a.z +
b.z); }
< Vec operator-(Vec a, Vec b) { return Vec(a.x - b.x, a.y - b.y, a.z -
b.z); }
< Vec operator*(double a, Vec b) { return Vec(a * b.x, a * b.y, a *
b.z); }
< double dot(Vec a, Vec b) { return a.x*b.x + a.y*b.y + a.z*b.z; }
< Vec unitise(Vec a) { return (1 / sqrt(dot(a, a))) * a; }
---
>
> inline Vec operator+(const Vec &a, const Vec &b) { return Vec(a.x +
b.x, a.y + b.y, a.z + b.z); }
> inline Vec operator-(const Vec &a, const Vec &b) { return Vec(a.x -
b.x, a.y - b.y, a.z - b.z); }
> inline Vec operator*(const double a, const Vec &b) { return Vec(a *
b.x, a * b.y, a * b.z); }
> inline double dot(const Vec &a, const Vec &b) { return a.x*b.x +
a.y*b.y + a.z*b.z; }
> inline Vec unitise(const Vec &a) { return (1 / sqrt(dot(a, a))) * a;
}
>
# g++ -g -O3 ray.cpp -o ray
# time ./ray >tt
real 0m10.493s
user 0m10.483s
sys 0m0.005s
If your point was to show how to write unbelievably inneficient code in
C++, you've succeeded.
Can't wait for a C++ vs Java comparison.
Matt Pharr
tbp...@gmail.com writes:
> Jon Harrop wrote:
>> I then cut this program down to a 66-line OCaml program and ported it
> into a
>> 97-line C++ program. These programs are compared on this page:
>>
>> http://www.ffconsultancy.com/free/ray_tracer/comparison.html
> Is that comparison supposed to be fair in any way whatsoever?
>
> elided]
>
> If your point was to show how to write unbelievably inneficient code in
> C++, you've succeeded.
I think the point was to have some fun and to make a little demo that shows
how a functional language like OCaml can give rise to nice short expressive
programs while still delivering competitive performance. And for those
goals, the OP succeeded!
Note also that he mentions on his web page that he's doing this for "the
computer language shootout benchmarks", http://shootout.alioth.debian.org/.
That page clearly points out all the caveats around benchmarks like these
(to the point that they can be effectively meaningless in the real world.)
Specifically:
>> about The Language Shootout Benchmarks
>>
>>Our goals are to learn about programming languages, compare their
>>performance in various (possibly meaningless) ways and, most importantly,
>>have some fun!
I've at least had some fun looking at his implementation!
-matt
--
Matt Pharr ma...@pharr.org <URL:http://graphics.stanford.edu/~mmp>
=======================================================================
In a cruel and evil world, being cynical can allow you to get some
entertainment out of it. --Daniel Waters
tbp
> > [demonstration that "inline" vector math makes C++ programs faster
> > elided]
a) he's telling the OCaml compiler to inline everything it can while
not giving its c++ counterpart opportunities to do so; fairness?
b) more importantly he's passing large structures by _value_ and i've
turned them into _references_.
In fact if you do the same for other hotpath functions as in:
19c21
< virtual void intersect(double &, Vec &, Ray) = 0;
---
> virtual void intersect(double &, Vec &, const Ray &) = 0;
26c28
< double ray_sphere(Ray ray) { // Intersection of a ray with a sphere
---
> double ray_sphere(const Ray &ray) { // Intersection of a ray with a
sphere
35c37
< void intersect(double &lambda, Vec &normal, Ray ray) {
---
> void intersect(double &lambda, Vec &normal, const Ray &ray) {
50c52
< void intersect(double &lambda, Vec &normal, Ray ray) {
---
> void intersect(double &lambda, Vec &normal, const Ray &ray) {
then the c++ version becomes 50% faster than originally (and in fact,
unconditionally faster than the OCaml version), without adding a single
line.
You'll notice i haven't even fixed the gratitious use of virtual
functions, and other details that force standard compliant c++
compilers to pessimize a lot when facing such... hmm.. source.
I'm sure OCaml is a wonderful & expressive language, but crafting such
a fishy c++ equivalent to make it shine (because you have an agenda) is
a despicable practice.
> > If your point was to show how to write unbelievably inneficient
code in
> > C++, you've succeeded.
>
> I think the point was to have some fun and to make a little demo that
shows
> how a functional language like OCaml can give rise to nice short
expressive
> programs while still delivering competitive performance. And for
those
> goals, the OP succeeded!
>
> Note also that he mentions on his web page that he's doing this for
"the
> computer language shootout benchmarks",
http://shootout.alioth.debian.org/.
> That page clearly points out all the caveats around benchmarks like
these
> (to the point that they can be effectively meaningless in the real
world.)
I have nothing against such benchmarks when they are not clearly
rigged.
But i guess the main point of the whole operation wasn't fair
benchmarking, education, or fun but selling books: "This program has
demonstrated some of the ways that OCaml improves upon C++. For a
thorough introduction to OCaml, read our book "Objective CAML for
Scientists".
Or put another way, that's a sad PR stunt.
Have a nice day.
Jon Harrop
> Matt Pharr wrote:
>> > [demonstration that "inline" vector math makes C++ programs faster
>> > elided]
> You've missed the point.
My post is more about clarity and less about performance. Your proposed
alterations make the C++ implementation significantly longer and more
obfuscated and your littering of the source code with "inline" actually
slows the program down on my Athlon t-bird.
Moreover, your choices of initial optimisation are decidedly suboptimal. I'd
have gone for:
1. Terminating the intersection of shadow rays when the first intersection
is found (instead of finding the closest intersection).
2. Use an implicit scene, to avoid storing it explicitly.
3. Use single-precision storage (or no storage at all).
Both of these optimisations will give much bigger performance improvements
when averaged over different platforms and architectures.
As it happens, I had already implemented all of your optimisations and all
of these optimisations to both implementations before I made my original
post. None of them add anything new. Sometimes C++ is faster, sometimes
OCaml is faster. OCaml is always much more succinct. You can keep doing
this until you are blue in the face but I don't think you'll learn much of
interest.
> a) he's telling the OCaml compiler to inline everything it can while
> not giving its c++ counterpart opportunities to do so; fairness?
Both ocamlopt and g++ were allowed to inline code.
As OCaml is for symbolic use, it performs no inlining by default (unlike
g++) so it is common to ask OCaml to do some inlining on numerical code.
Thus, specifying "-inline 100" actually makes for a fairer comparison.
> b) more importantly he's passing large structures by _value_ and i've
> turned them into _references_.
Both implementations use pass by value as this is clearer, shorter and (as a
consequence) more common in real code.
> In fact if you do the same for other hotpath functions as in:
> ...
> then the c++ version becomes 50% faster than originally
If you want a fair comparison then you should also make equivalent changes
to the OCaml implementation.
> (and in fact, unconditionally faster than the OCaml version),
You are trying to compare optimised C++ against unoptimised OCaml, which
would be unfair. More worryingly, you seem to have skipped the part of the
experiment where you actually measure something.
> without adding a single line.
When restricted to 80 columns, your optimisations add several lines. Indeed,
they push the C++ program over 100 LOC limit imposed by the creators of the
shootout.
In fact, I'd already implemented your optimisations (and many more effective
optimisations) and had to throw them out because of this limit. Note that
the OCaml has 30 more lines with which to optimise. My optimised <100 LOC
OCaml program is much faster than my optimised <100 LOC C++ on all
platforms.
> You'll notice i haven't even fixed the gratitious use of virtual
> functions,
What would you recommend instead?
I chose an inheritance hierarchy because this is the closest C++ equivalent
to a variant type which I see in typical C++ programs.
Another optimisation that I made was to replace inheritance with a single
Scene struct which represented a group of child scenes or a sphere when it
had no children. I threw this out as well because it is not a fair
equivalent to OCaml's variant type. For example, you could not specify a
colour for each sphere without also specifying colours for all bounding
volumes.
> and other details
Can you elaborate on these other "details"?
> that force standard compliant c++
> compilers to pessimize a lot when facing such... hmm.. source.
Manually implementing pass constant by reference is both obfuscating and
error-prone whilst being theoretically unnecessary. The fact that the OCaml
compiler does this optimisation for you when the GNU C++ compiler does not
might interest some people.
> I'm sure OCaml is a wonderful & expressive language, but crafting such
> a fishy c++ equivalent to make it shine (because you have an agenda) is
> a despicable practice.
Everyone will be free to contribute to the shootout version of the ray
tracer. Perhaps you would like to contribute a less "fishy" C++
implementation?
Jon Harrop
> Great work. What is default floating point type used by OCaml
> compiler?
Double precision.
> I got running time of about 36s for scene level 3 on
> Pentium3 machine and running time of about 25s when all references to
> "double" type in C++ code changed to "float" (and float is usually
> "good enough" for ray tracing)...
I agree that single precision is usually good enough for ray tracing.
However, when I tried the same optimisation I found that it made the
program much faster on 32-bit x86 and much slower on 64-bit AMD64. I assume
this is a float alignment issue.
Also, it generates (numerically) different results, which makes it trickier
to validate against the OCaml.
Warp
> My post is more about clarity and less about performance. Your proposed
> alterations make the C++ implementation significantly longer and more
> obfuscated and your littering of the source code with "inline" actually
> slows the program down on my Athlon t-bird.
Then your compiler is stupid. Get a better one.
The keyword 'inline' is not a command for the compiler that says it must
inline the function or else it will be deleted immediately. And the other
way around, the lack of 'inline' keyword by no means stops the compiler
from inlining the function if it can and it estimates it would be useful.
In *theory* 'inline' is a hint to the compiler. In practice it has nothing
to do with inlining. Most of the modern compilers will simply ignore that
keyword when they are estimating whether a function should be inlined or
not. In this context 'inline' is a no-op. It's like the 'register'
keyword which is nowadays obsolete and a no-op.
What 'inline' *does* have some effect on is the linker. If a function
with an identical declaration is implemented in more than one place
then in normal circumstances the linker will give an error (because it
can't know which one of them it should use). However, if the functions
were declared 'inline', then the compiler will tell the linker that
all of those functions are actually one and the same, and then the linker
will just use one of them.
And by the way, if you are afraid of "litter" in C++ code, then you
should switch to another programming language.
That "litter", as you call it, is not useless junk as you seem to imply.
It has a reason and a significance.
> > b) more importantly he's passing large structures by _value_ and i've
> > turned them into _references_.
> Both implementations use pass by value as this is clearer, shorter and (as a
> consequence) more common in real code.
I hope you are not passing megabyte-sized vectors by value...
If your opinion of "clearer" and "shorter" is that the & symbol should
not be used, then you might want to either check your priorities or
change your programming language (because C++ is not for you).
> If you want a fair comparison then you should also make equivalent changes
> to the OCaml implementation.
You want to compare which programming language runs inefficient code
faster?
> > You'll notice i haven't even fixed the gratitious use of virtual
> > functions,
> What would you recommend instead?
> I chose an inheritance hierarchy because this is the closest C++ equivalent
> to a variant type which I see in typical C++ programs.
Inheritance has no overhead in C++. Dynamic binding has.
While virtual function calls are pretty fast they are still slightly
slower than regular function calls. If you want maximum speed you should
not make time-critical functions virtual unless you are really forced to.
> Another optimisation that I made was to replace inheritance with a single
> Scene struct which represented a group of child scenes or a sphere when it
> had no children.
Inheritance causes no overhead. That kind of optimization is useless.
--
- Warp
tbp
> tbp wrote:
> My post is more about clarity and less about performance. Your
proposed
> alterations make the C++ implementation significantly longer and more
> obfuscated and your littering of the source code with "inline"
actually
> slows the program down on my Athlon t-bird.
I've put explicit inlines because, ie for Vec operators, you've put
them outside of the class scope, hence screwing the common heuristic
saying that a member function declared & defined within the class body
should be inlined (and that's a much stronger hint than using an
'inline' directive).
Anyway, that's decoration compared to other issues in your source.
> Moreover, your choices of initial optimisation are decidedly
suboptimal. I'd
> have gone for:
Hmm? I tinker with a handful of lines, cut the runtime in half and
that's suboptimal?
Interesting.
> 1. Terminating the intersection of shadow rays when the first
intersection
> is found (instead of finding the closest intersection).
>
> 2. Use an implicit scene, to avoid storing it explicitly.
>
> 3. Use single-precision storage (or no storage at all).
>
> Both of these optimisations will give much bigger performance
improvements
> when averaged over different platforms and architectures.
Your aproach & space partition (as presented) will never give runtime
performance worth the time it takes to implement efficiently.
Such a sphere flake should render in a blink of the eye on modern
platforms, see the SPD for reference.
Or this: http://www.uni-koblenz.de/~cg/publikationen/cp_raytrace.pdf
I've merely fixed the most glaring flaws in your implementation.
Period.
> As it happens, I had already implemented all of your optimisations
and all
> of these optimisations to both implementations before I made my
original
> post. None of them add anything new. Sometimes C++ is faster,
sometimes
> OCaml is faster. OCaml is always much more succinct. You can keep
doing
> this until you are blue in the face but I don't think you'll learn
much of
> interest.
I've never contested that OCaml is more succint, expressive or elegant.
I'm just saying that, given the horrible c++ implementation you've
presented, you have absolutely no authority or qualification to say
anything about the relative merit of c++ & OCaml performance wise.
> Both ocamlopt and g++ were allowed to inline code.
It's all about exposing oportunities to do so.
> As OCaml is for symbolic use, it performs no inlining by default
(unlike
> g++) so it is common to ask OCaml to do some inlining on numerical
code.
> Thus, specifying "-inline 100" actually makes for a fairer
comparison.
No. You've used none of the c++ idiom that actually helps a c++
compiler to guess what to inline or not.
Know your compiler (or in that case, your language).
> Both implementations use pass by value as this is clearer, shorter
and (as a
> consequence) more common in real code.
Absolute non sense.
I can't beleive you really think that, but if that's the case, then you
should get back to your homework.
> If you want a fair comparison then you should also make equivalent
changes
> to the OCaml implementation.
You're the one trying to sell books about OCaml...
> You are trying to compare optimised C++ against unoptimised OCaml,
which
> would be unfair. More worryingly, you seem to have skipped the part
of the
> experiment where you actually measure something.
I've provided patches. See for yourself.
> When restricted to 80 columns, your optimisations add several lines.
Indeed,
> they push the C++ program over 100 LOC limit imposed by the creators
of the
> shootout.
You want to fit in those artificial constraints? Shorten names or
define some macros.
What was the fuss about already?
> In fact, I'd already implemented your optimisations (and many more
effective
> optimisations) and had to throw them out because of this limit. Note
that
> the OCaml has 30 more lines with which to optimise. My optimised <100
LOC
> OCaml program is much faster than my optimised <100 LOC C++ on all
> platforms.
Sure. I mean you've shown such a level of C++ wizardry that we're
supposed to take your word for it?
> > You'll notice i haven't even fixed the gratitious use of virtual
> > functions,
>
> What would you recommend instead?
Not using virtual functions in the hotpath.
Primo, that's uncalled for as your "hierarchy" doesn't mix types at a
given level, each object's type is known upfront.
Secundo, virtual functions are expensive because they forbid inlining
and are implemented as indirect branches , like:
401c83: callq *0x8(%rax)
It doesn't make sense to use an expensive feature in the hotpath if you
can avoid it.
And you definitely can.
> I chose an inheritance hierarchy because this is the closest C++
equivalent
> to a variant type which I see in typical C++ programs.
>
> Another optimisation that I made was to replace inheritance with a
single
> Scene struct which represented a group of child scenes or a sphere
when it
> had no children. I threw this out as well because it is not a fair
> equivalent to OCaml's variant type. For example, you could not
specify a
> colour for each sphere without also specifying colours for all
bounding
> volumes.
That's totally orthogonal.
Please do yourself a favor and read a good book about C++.
> > and other details
> Can you elaborate on these other "details"?
Strictly speaking of your implementation details (and not algorithm or
anything), you've also failed the const correctness test for example.
But i suppose you're going to say say that it's cruft or obfuscation.
> > that force standard compliant c++
> > compilers to pessimize a lot when facing such... hmm.. source.
>
> Manually implementing pass constant by reference is both obfuscating
and
> error-prone whilst being theoretically unnecessary. The fact that the
OCaml
> compiler does this optimisation for you when the GNU C++ compiler
does not
> might interest some people.
I think you've just discovered that OCaml & C++ are two different
language that require the programmer to express things in a totally
different way.
And before you propose ways to improve C++ you should learn the
language first.
> Everyone will be free to contribute to the shootout version of the
ray
> tracer. Perhaps you would like to contribute a less "fishy" C++
> implementation?
I don't have books to sell.
If/When i write a raytracer, i choose a path that at least has some
chances to perform decently.
See Wald & Havran.
Jon Harrop
> Jon Harrop <use...@jdh30.plus.com> wrote:
>> My post is more about clarity and less about performance. Your proposed
>> alterations make the C++ implementation significantly longer and more
>> obfuscated and your littering of the source code with "inline" actually
>> slows the program down on my Athlon t-bird.
>
> Then your compiler is stupid. Get a better one.
I agree that the compiler is stupid but I'm comparing g++ as this is a
common compiler. A better C++ compiler is likely to produce faster code, of
course.
> And by the way, if you are afraid of "litter" in C++ code, then you
> should switch to another programming language.
Absolutely. :-)
> That "litter", as you call it, is not useless junk as you seem to imply.
> It has a reason and a significance.
It also has unexpected platform dependent results, as "tbp" has discovered.
>> > b) more importantly he's passing large structures by _value_ and i've
>> > turned them into _references_.
>
>> Both implementations use pass by value as this is clearer, shorter and
>> (as a consequence) more common in real code.
>
> I hope you are not passing megabyte-sized vectors by value...
No, they are small, constant-sized vectors (3D), so pass by value does not
affect the asymptotic algorithmic complexity.
> If your opinion of "clearer" and "shorter" is that the & symbol should
> not be used, then you might want to either check your priorities or
> change your programming language (because C++ is not for you).
Unfortunately, you cannot just use "&" because it will take a reference to a
temporary, giving an error like:
ray.cpp: In member function `virtual void Sphere::intersect(double&, Vec&,
Ray)
':
ray.cpp:39: error: no match for 'operator+' in 'ray.Ray::orig +
operator*(double, Vec)(((&ray) + 24))'
ray.cpp:12: error: candidates are: Vec operator+(Vec&, Vec&)
/usr/include/c++/3.3/bits/stl_bvector.h:261: error:
std::_Bit_const_iterator std::operator+(int, const
std::_Bit_const_iterator&)
/usr/include/c++/3.3/bits/stl_bvector.h:202: error:
std::_Bit_iterator std::operator+(int, const std::_Bit_iterator&)
ray.cpp: In function `double ray_trace(double, Vec, Ray, Scene*)':
ray.cpp:62: error: no match for 'operator+' in 'ray.Ray::orig +
operator*(double, Vec)(((&ray) + 24))'
ray.cpp:12: error: candidates are: Vec operator+(Vec&, Vec&)
/usr/include/c++/3.3/bits/stl_bvector.h:261: error:
std::_Bit_const_iterator std::operator+(int, const
std::_Bit_const_iterator&)
/usr/include/c++/3.3/bits/stl_bvector.h:202: error:
std::_Bit_iterator std::operator+(int, const std::_Bit_iterator&)
You must either use "const ... &" instead, or make sure that all values are
created explicitly as local variables in the caller (rather than being
inline subexpressions in the function call).
It is precisely this obfuscation (and the fact that these are tiny
structures, and that this is not necessarily an optimisation, and that it
could be done by the compiler) that made me choose pass by value.
>> If you want a fair comparison then you should also make equivalent
>> changes to the OCaml implementation.
>
> You want to compare which programming language runs inefficient code
> faster?
If you completely optimise both implementations then you'll end up with the
same assembler. There is only value in comparing simpler implementations in
both languages. So yes, they will not be completely optimised but
subjectively "like" real code.
>> I chose an inheritance hierarchy because this is the closest C++
>> equivalent to a variant type which I see in typical C++ programs.
>
> Inheritance has no overhead in C++. Dynamic binding has.
> While virtual function calls are pretty fast they are still slightly
> slower than regular function calls. If you want maximum speed you should
> not make time-critical functions virtual unless you are really forced to.
I believe I am "really forced to" as virtual lookup is used instead of
switching on the tag of the variant type/union. In C++, virtual functions
are more elegant (and more common in real code) than C-style tagged unions.
>> Another optimisation that I made was to replace inheritance with a single
>> Scene struct which represented a group of child scenes or a sphere when
>> it had no children.
>
> Inheritance causes no overhead. That kind of optimization is useless.
Getting rid of inheritance gets rid of the virtual functions. However, it
necessarily incurs an equivalent test somewhere else.
Jon Harrop
>> When restricted to 80 columns, your optimisations add several lines.
> Indeed,
>> they push the C++ program over 100 LOC limit imposed by the creators
> of the
>> shootout.
> You want to fit in those artificial constraints?
Yes, it actually says that at the top of the web page.
> Shorten names or define some macros. What was the fuss about already?
If you recall, the fuss was about clarity. So shortening names and defining
macros doesn't really seem appropriate.
>> In fact, I'd already implemented your optimisations (and many more
> effective
>> optimisations) and had to throw them out because of this limit. Note
> that
>> the OCaml has 30 more lines with which to optimise. My optimised <100
> LOC
>> OCaml program is much faster than my optimised <100 LOC C++ on all
>> platforms.
> Sure. I mean you've shown such a level of C++ wizardry that we're
> supposed to take your word for it?
I'd be very happy to include a better C++ version, if you have one.
>> > You'll notice i haven't even fixed the gratitious use of virtual
>> > functions,
>>
>> What would you recommend instead?
> Not using virtual functions in the hotpath.
> Primo, that's uncalled for as your "hierarchy" doesn't mix types at a
> given level, each object's type is known upfront.
So you're talking about specialising the scene tree so that it can only
represent this scene, rather than scenes in general?
> Secundo, virtual functions are expensive because they forbid inlining
> and are implemented as indirect branches , like:
> 401c83: callq *0x8(%rax)
> It doesn't make sense to use an expensive feature in the hotpath if you
> can avoid it.
> And you definitely can.
Actually, this is another optimisation that I had already done. It makes no
significant difference to the performance because this isn't actually on
the "hotpath", as you call it. Most of the time is spent elsewhere.
>> > and other details
>> Can you elaborate on these other "details"?
> Strictly speaking of your implementation details (and not algorithm or
> anything), you've also failed the const correctness test for example.
> But i suppose you're going to say say that it's cruft or obfuscation.
Yes, that was also done deliberately to squeeze the C++ into 100 LOC.
>> Everyone will be free to contribute to the shootout version of the
> ray
>> tracer. Perhaps you would like to contribute a less "fishy" C++
>> implementation?
> I don't have books to sell.
> If/When i write a raytracer, i choose a path that at least has some
> chances to perform decently.
> See Wald & Havran.
I didn't think so. :-)
Jon Harrop
> Such a sphere flake should render in a blink of the eye on modern
> platforms, see the SPD for reference.
> Or this: http://www.uni-koblenz.de/~cg/publikationen/cp_raytrace.pdf
Although I'm sure that a production ray tracer will be a lot faster than my
tutorial ray tracer, the paper you cite describes a high-level program
which is not much faster than mine. The faster alternative is actually
written in assembler (using SIMD) and also fails to reach "interactive
frame rates" on such scenes.
Warp
> I agree that the compiler is stupid but I'm comparing g++ as this is a
> common compiler. A better C++ compiler is likely to produce faster code, of
> course.
Are you using the optimization flags for your platform?
> Unfortunately, you cannot just use "&" because it will take a reference to a
> temporary, giving an error like:
const &
You very seldom need non-const references.
> > Inheritance has no overhead in C++. Dynamic binding has.
> > While virtual function calls are pretty fast they are still slightly
> > slower than regular function calls. If you want maximum speed you should
> > not make time-critical functions virtual unless you are really forced to.
> I believe I am "really forced to" as virtual lookup is used instead of
> switching on the tag of the variant type/union. In C++, virtual functions
> are more elegant (and more common in real code) than C-style tagged unions.
I was not suggesting a "tagged union" as an alternative. That would
actually be slower than a dynamically bound class.
I was suggesting as an alternative a different design where you do not
need dynamic binding.
> > Inheritance causes no overhead. That kind of optimization is useless.
> Getting rid of inheritance gets rid of the virtual functions.
Well, yes, but you are using the wrong means. You don't need to drop
inheritance simply to stop using dynamic binding (for the latter you
simply need to drop the 'virtual' keywords; inheritance has nothing
to do with that).
--
- Warp
tbp
> tbp wrote:
> > Such a sphere flake should render in a blink of the eye on modern
> > platforms, see the SPD for reference.
> > Or this:
http://www.uni-koblenz.de/~cg/publikationen/cp_raytrace.pdf
>
> Although I'm sure that a production ray tracer will be a lot faster
than my
> tutorial ray tracer, the paper you cite describes a high-level
program
> which is not much faster than mine. The faster alternative is
actually
> written in assembler (using SIMD) and also fails to reach
"interactive
> frame rates" on such scenes.
Not assembler, intrinsics.
And the whole thing is not about SIMD but ray coherence & BVH which are
notions that have nothing to do with a specific platform.
It just happens to fit into the SIMD model quite well.
tbp
> > You want to fit in those artificial constraints?
> Yes, it actually says that at the top of the web page.
> > Shorten names or define some macros. What was the fuss about
already?
> If you recall, the fuss was about clarity. So shortening names and
defining
> macros doesn't really seem appropriate.
incidentally happens to be the worst way to implement it in C++.
> I'd be very happy to include a better C++ version, if you have one.
Better in what metric? I'm cautious because it seems we live in two
distinct reality.
> So you're talking about specialising the scene tree so that it can
only
> represent this scene, rather than scenes in general?
Scenes in general?
Please correct me if i'm wrong but it seems that your program is about
a sphere BVH enclosing... spheres.
Is that your idea of generality?
So, please focus on the topic at hand and you'll notice that using
dynamic binding is superfluous.
> > Secundo, virtual functions are expensive because they forbid
inlining
> > and are implemented as indirect branches , like:
> > 401c83: callq *0x8(%rax)
> > It doesn't make sense to use an expensive feature in the hotpath if
you
> > can avoid it.
> > And you definitely can.
>
> Actually, this is another optimisation that I had already done. It
makes no
> significant difference to the performance because this isn't actually
on
> the "hotpath", as you call it. Most of the time is spent elsewhere.
I see:
struct Scene { // Abstract base class representing a scene
virtual void intersect(double &, Vec &, Ray) = 0;
So you're saying intersections aren't part of the hotpath?
Well you maybe right, after all it's only called a few million times.
> >> Can you elaborate on these other "details"?
> > Strictly speaking of your implementation details (and not algorithm
or
> > anything), you've also failed the const correctness test for
example.
> > But i suppose you're going to say say that it's cruft or
obfuscation.
>
> Yes, that was also done deliberately to squeeze the C++ into 100 LOC.
Bingo!
And i suppose that the fact that this emasculated C++ version performed
slower than its OCaml counterpart was merely a side effect?
What's next? You're going to try to sell me a bridge?
Oops. A book first, then a bridge.
> I didn't think so. :-)
Oh! That was a test and i failed it?
Damn.
Jon Harrop
> Are you using the optimization flags for your platform?
Yes, I've played with various optimisation settings with both compilers on
both platforms. I didn't see anything particularly interesting when
tweaking them. The performance varies by about 20% but the ratio of the
performances stays the same (between OCaml and C++).
> You very seldom need non-const references.
In point of fact, the program uses non-const references to return multiple
results from functions. IMHO, this is often done in C++ (and pass by
non-const pointer in C) because it is too tedious to define a new type for
the return values of each function in C/C++.
As there are currently only two values being returned by the intersect
routine, I could have used STL pair<> instead. I might investigate that (I
think it will be simpler but slower). I have seen code which uses
pair<pair<...>, ...> to return three values but I don't much like it.
>> I believe I am "really forced to" as virtual lookup is used instead of
>> switching on the tag of the variant type/union. In C++, virtual functions
>> are more elegant (and more common in real code) than C-style tagged
>> unions.
>
> I was not suggesting a "tagged union" as an alternative. That would
> actually be slower than a dynamically bound class.
> I was suggesting as an alternative a different design where you do not
> need dynamic binding.
I have actually already done this "optimisation" because I also thought it
would speed things up. It turns out that the virtual lookup is taking an
insignificant amount of time. I ditched the resulting code because it is
less general and less clear whilst being no faster. It is significantly
shorter though.
Jon Harrop
> Jon Harrop wrote:
>> If you recall, the fuss was about clarity. So shortening names and
>> defining
>> macros doesn't really seem appropriate.
> If i sum it up, nothing's appropriate but what you did which
> incidentally happens to be the worst way to implement it in C++.
Perhaps there is a better way but I haven't seen it.
>> I'd be very happy to include a better C++ version, if you have one.
> Better in what metric? I'm cautious because it seems we live in two
> distinct reality.
Better => Faster whilst remaining clear and satisfying the requirements of
the shootout.
>> So you're talking about specialising the scene tree so that it can
> only
>> represent this scene, rather than scenes in general?
> Scenes in general?
> Please correct me if i'm wrong but it seems that your program is about
> a sphere BVH enclosing... spheres.
> Is that your idea of generality?
In both cases, the types of objects and bounding volumes are extensible.
To add a new kind of object to the C++, you create a new class derived from
Scene which implements its own "intersect" member function.
To add a new kind of object to the OCaml, you supplement the variant type
"scene" with another constructor and supplement the pattern match in the
"intersect" function with a case to handle your new object.
Similarly for new types of bounding volume.
Perhaps it would be interesting to add a new kind of object to both
implementations, to see how this is done in each language?
>> Actually, this is another optimisation that I had already done. It
> makes no
>> significant difference to the performance because this isn't actually
> on
>> the "hotpath", as you call it. Most of the time is spent elsewhere.
> I see:
> struct Scene { // Abstract base class representing a scene
> virtual void intersect(double &, Vec &, Ray) = 0;
>
> So you're saying intersections aren't part of the hotpath?
> Well you maybe right, after all it's only called a few million times.
Yes, exactly. It isn't really that surprising when you consider that
"ray_sphere" is called ~10x more than "intersect" and it contains several
branches, lots of floating point maths and function calls.
>> Yes, that was also done deliberately to squeeze the C++ into 100 LOC.
> Bingo!
>
> And i suppose that the fact that this emasculated C++ version performed
> slower than its OCaml counterpart was merely a side effect?
The C++ outperforms the OCaml on AMD64 (and IA64). I'd say that the
performance is comparable in all cases.
Jon Harrop
> Not assembler, intrinsics.
> And the whole thing is not about SIMD but ray coherence & BVH which are
> notions that have nothing to do with a specific platform.
> It just happens to fit into the SIMD model quite well.
Ok. I don't know anything about intrinsics but they don't seem to be
relevant to my comparison of a tutorial ray tracer in C++ and OCaml.
As I've said, if you want to optimise my ray tracer (or any other program)
then I would recommend algorithmic optimisations over low-level ones like
SIMD, cache coherency etc. There are plenty more algorithmic optimisations
to be done (although I've already done the main one - hierarchical bounding
volumes).
In this case, it looks as though you'd want to write specialised C code for
certain platforms. This could be called from C++ or OCaml and it would
speed both programs up by roughly the same amount. This might be a good
idea though, as OCaml is clearly preferable for the higher-level, not-
performance-critical bulk of the code. So you could get rid of C++
altogether. :-)
Warp
> > You very seldom need non-const references.
> In point of fact, the program uses non-const references to return multiple
> results from functions. IMHO, this is often done in C++ (and pass by
> non-const pointer in C) because it is too tedious to define a new type for
> the return values of each function in C/C++.
Often? It is usually rare that you have to return several values
from a function so that usually all of them are needed.
The problem of making one function which returns many values versus
making several functions, each one returning one value, is that when
you only need that one value returning many is useless overhead and
forces you to make clumsy code (you have to create dummy variables to
hold those return variables you don't need and then immediately discard
them).
When there's a logical reason to return many values at once, returning
a class or struct as the return value of the function may in fact be a
good idea. If you call the function like this:
ReturnValues values = theFunction(whatever);
where ReturnValues is a class or struct, C++ compilers (or at least gcc;
probably also most of the other compilers) will actually make that very
efficiently (they will perform the so-called return value optimization).
The return value optimization means that, in a code like the above,
the function being called actually builds its return value directly
into 'values' and will not create any temporaries nor will it call any
copy constructors. Perhaps a bit surprisingly this works even if the
function in question is not inline and its implementation is in a
completely separate compilation unit (ie. object file) which is just
linked later to the binary (and in fact it works even if the function
in question is in a dynamically loaded library).
It might even be that the return value optimization technique produces
slightly faster code than the function returning its values by writing
to references given to it as parameters.
> As there are currently only two values being returned by the intersect
> routine, I could have used STL pair<> instead. I might investigate that (I
> think it will be simpler but slower). I have seen code which uses
> pair<pair<...>, ...> to return three values but I don't much like it.
How about just creating your own struct or class?
> I have actually already done this "optimisation" because I also thought it
> would speed things up. It turns out that the virtual lookup is taking an
> insignificant amount of time.
Yes, virtual function calls are very fast in C++. However, what is fastest
depends on what you are comparing it to.
If you are comparing it to an equivalent functionality implemented with
a big switch-case block, then you should definitely use the virtual
functions instead.
However, if you can design the program so that the function to be
called can be resolved at compilation time (instead of at runtime as
happens with virtual functions or switch-case blocks) then you might
get a slight speedup, specially if the function is short and called
millions of times.
--
- Warp
Jon Harrop
>> multiple results from functions. IMHO, this is often done in C++ (and
>> pass by non-const pointer in C) because it is too tedious to define a new
>> type for the return values of each function in C/C++.
>
> Often? It is usually rare that you have to return several values
> from a function so that usually all of them are needed.
It may seem rare when you use a programming language that makes it
difficult, but lots of useful functions return multiple values.
> The problem of making one function which returns many values versus
> making several functions, each one returning one value, is that when
> you only need that one value returning many is useless overhead and
> forces you to make clumsy code (you have to create dummy variables to
> hold those return variables you don't need and then immediately discard
> them).
In the case of the "intersect" function, you would not split it into two
functions to compute the parameter and normal vector of the intersection,
for example.
> When there's a logical reason to return many values at once, returning
> a class or struct as the return value of the function may in fact be a
> good idea. If you call the function like this:
>
> ReturnValues values = theFunction(whatever);
>
> where ReturnValues is a class or struct, C++ compilers (or at least gcc;
> probably also most of the other compilers) will actually make that very
> efficiently (they will perform the so-called return value optimization).
Creating structs for the return values of every function that returns
multiple values may be a good idea, but the verbosity and clumsiness of
doing this in C++ makes it quite undesirable.
> The return value optimization means that...
If it is performed, yes.
>> As there are currently only two values being returned by the intersect
>> routine, I could have used STL pair<> instead. I might investigate that
>> (I think it will be simpler but slower). I have seen code which uses
>> pair<pair<...>, ...> to return three values but I don't much like it.
>
> How about just creating your own struct or class?
That would generate quite a bit more code and would probably break the 100
LOC limit.
> However, if you can design the program so that the function to be
> called can be resolved at compilation time (instead of at runtime as
> happens with virtual functions or switch-case blocks) then you might
> get a slight speedup, specially if the function is short and called
> millions of times.
I do not believe that can be done in this case (without breaking generality
and the equivalence with the OCaml implementation).
Warp
> Creating structs for the return values of every function that returns
> multiple values may be a good idea, but the verbosity and clumsiness of
> doing this in C++ makes it quite undesirable.
You mean modular?
Have you ever heard of abstraction?
--
- Warp
tbp
> tbp wrote:
> > incidentally happens to be the worst way to implement it in C++.
>
> Perhaps there is a better way but I haven't seen it.
>
> >> I'd be very happy to include a better C++ version, if you have
one.
> > Better in what metric? I'm cautious because it seems we live in two
> > distinct reality.
>
> Better => Faster whilst remaining clear and satisfying the
requirements of
> the shootout.
and, as you insist on discarding my previous patches & measures, i've
taken the time this morning to make this up:
http://ompf.org/ray/miniray.cpp
# grep -e virtual -e inline miniray.cpp
# wc -l miniray.cpp
97 miniray.cpp
# g++ -O2 -ffast-math miniray.cpp -o miniray
# time ./miniray >pix.ppm
real 0m5.862s
user 0m5.800s
sys 0m0.006s
# cmp pix.orig.ppm pix.ppm; echo $?
0
# composite -compose difference pix.orig.ppm pix.ppm diff.ppm
# display diff.ppm
[contemplate a black picture]
So, on my box, the original version took ~13s then ~10s after some
patching, ~8s after some more patching and now we're at ~6s with a
const correct non virtual rule compliant version.
Before you start nitpicking please consider that:
a) this version complies with the shootout rules AFAIK.
b) it outputs the exact same picture.
c) i had to raise the epsilon a bit to avoid acnea, it has to do with
the way i compute normals but i'm not going to fix that.
d) i removed the virtualness in the most straightforward way and it's
certainly not aesthetically pleasing
e) i'm not going to spend ages on this, hence c) & d).
f) it's algorithmically equivalent to your version (ie primary and
shadow rays aren't special cased etc...)
> In both cases, the types of objects and bounding volumes are
extensible.
That means nothing. Everything is extensible.
And in this *particular case*, using dynamic binding is, again,
superfluous.
I've resolved that in a really crude way, but that's just one
possibility.
> > So you're saying intersections aren't part of the hotpath?
> > Well you maybe right, after all it's only called a few million
times.
>
> Yes, exactly. It isn't really that surprising when you consider that
> "ray_sphere" is called ~10x more than "intersect" and it contains
several
> branches, lots of floating point maths and function calls.
I was sarcastic; "ray_sphere" calls that interection routine which
cannot be inlined because it's virtual etc...
Also, the fact is that this kind of sphere intersection can easily be
made branchless with something like:
if ((d < 0.) | (t2 < 0.)) return infinity;
else return t1 > 0. ? t1 : t2;
But you want to branch away from the square root as it is what's
wasting most cycles.
> >> Yes, that was also done deliberately to squeeze the C++ into 100
LOC.
> > Bingo!
> >
> > And i suppose that the fact that this emasculated C++ version
performed
> > slower than its OCaml counterpart was merely a side effect?
>
> The C++ outperforms the OCaml on AMD64 (and IA64). I'd say that the
> performance is comparable in all cases.
Again, your benchmark as presented it totally unfair.
tbp
> Ok. I don't know anything about intrinsics but they don't seem to be
> relevant to my comparison of a tutorial ray tracer in C++ and OCaml.
> As I've said, if you want to optimise my ray tracer (or any other
program)
> then I would recommend algorithmic optimisations over low-level ones
like
Huh. Really?
> SIMD, cache coherency etc. There are plenty more algorithmic
optimisations
> to be done (although I've already done the main one - hierarchical
bounding
> volumes).
So what? The fact is that implementation quality matters just as much.
A program has to run on a given machine at some point.
> In this case, it looks as though you'd want to write specialised C
code for
> certain platforms. This could be called from C++ or OCaml and it
would
> speed both programs up by roughly the same amount. This might be a
good
> idea though, as OCaml is clearly preferable for the higher-level,
not-
> performance-critical bulk of the code. So you could get rid of C++
> altogether. :-)
If you want to see it that way, that's fine.
But i wouldn't say that spending most cycles in c++ called from "higher
level" amounts to getting rid of it. Au contraire.
Jon Harrop
> Jon Harrop <use...@jdh30.plus.com> wrote:
>> Creating structs for the return values of every function that returns
>> multiple values may be a good idea, but the verbosity and clumsiness of
>> doing this in C++ makes it quite undesirable.
>
> You mean modular?
No.
> Have you ever heard of abstraction?
Yes.
For example, the following OCaml function returns the given number as well
as its square and cube:
let f n = n, n*n, n*n*n;;
The C++ equivalent of this using a struct is:
struct res { int n, n2, n3; };
res f(int n) {
res ans;
ans.n = n;
ans.n2 = n*n;
ans.n3 = n*n*n;
return ans;
}
The C++ equivalent of this using a STL pair is:
pair<int, pair<int, int> > f(int n) {
return pair<int, pair<int, int> >(n, n*n, n*n*n);
}
Do you see what I mean by verbose and clumsy?
Jon Harrop
>> SIMD, cache coherency etc. There are plenty more algorithmic
> optimisations
>> to be done (although I've already done the main one - hierarchical
> bounding
>> volumes).
> So what? The fact is that implementation quality matters just as much.
No, it doesn't. If you can reduce your algorithm from O(n) to O(log n) then
you should do so before you rewrite it in assembler.
> A program has to run on a given machine at some point.
True. I'm not saying that it is not worthwhile to optimise a ray tracer with
SIMD, just that it would be premature to do so on my ray tracer right now.
I would be very interested to see a highly optimised version of my ray
tracer. I think a lot of people could learn from a set of web pages which
show how to optimise the programs. If you're willing to help and I can find
the time then I'd like to do this sometime.
>> In this case, it looks as though you'd want to write specialised C
> code for
>> certain platforms. This could be called from C++ or OCaml and it
> would
>> speed both programs up by roughly the same amount. This might be a
> good
>> idea though, as OCaml is clearly preferable for the higher-level,
> not-
>> performance-critical bulk of the code. So you could get rid of C++
>> altogether. :-)
> If you want to see it that way, that's fine.
> But i wouldn't say that spending most cycles in c++ called from "higher
> level" amounts to getting rid of it. Au contraire.
I wouldn't call this kind of stuff C++:
a = _mm_load_ps(inp_sse1);
b = _mm_load_ps(inp_sse2);
c = _mm_add_ps(a, b);
_mm_store_ps(out_sse, c);
Jon Harrop
> I'll leave non objective criterions aside (how do you measure clarity?)
> and, as you insist on discarding my previous patches & measures, i've
> taken the time this morning to make this up:
I appreciate it, thank you.
On my 1.2GHz Athlon t-bird I get:
$ g++-3.4 -Wall -O3 -ffast-math -funroll-all-loops ray.cpp -o ray
$ time ./ray >image.ppm
real 0m55.256s
user 0m54.860s
sys 0m0.080s
$ g++-3.4 -Wall -O3 -ffast-math -funroll-all-loops miniray.cpp -o miniray
miniray.cpp:18: warning: division by zero in `1.0e+0 / 0.'
$ time ./miniray >image.ppm
real 0m28.587s
user 0m28.400s
sys 0m0.040s
So your optimised version is 1.9x faster.
My optimised version is basically the same speed as yours (although quite
different):
$ g++-3.4 -Wall -O3 -march=athlon-tbird -ffast-math -funroll-all-loops
ray4.cpp -o ray4
$ time ./ray4 >image.ppm
real 0m29.960s
user 0m27.150s
sys 0m0.040s
> Before you start nitpicking please consider that:
> a) this version complies with the shootout rules AFAIK.
Yes.
> b) it outputs the exact same picture.
Yes.
> c) i had to raise the epsilon a bit to avoid acnea, it has to do with
> the way i compute normals but i'm not going to fix that.
Your epsilon is actually smaller than my delta (which was sqrt epsilon).
> d) i removed the virtualness in the most straightforward way and it's
> certainly not aesthetically pleasing
I'll study this in more detail.
> e) i'm not going to spend ages on this, hence c) & d).
What better solution would you recommend for (d)?
My optimised version uses:
struct Tree {
vector<Tree> child;
Vec center;
double radius;
Tree(Vec c, double r) : child(), center(c), radius(r) {}
};
which I don't much like.
> f) it's algorithmically equivalent to your version (ie primary and
> shadow rays aren't special cased etc...)
My optimised version has a specialised intersect function for shadow rays
(but still fits in <100 LOC and is just a clear, IMHO).
>> In both cases, the types of objects and bounding volumes are
> extensible.
> That means nothing. Everything is extensible.
> And in this *particular case*, using dynamic binding is, again,
> superfluous.
> I've resolved that in a really crude way, but that's just one
> possibility.
Can you describe how a new kind of object/bound would be implemented?
>> Yes, exactly. It isn't really that surprising when you consider that
>> "ray_sphere" is called ~10x more than "intersect" and it contains
> several
>> branches, lots of floating point maths and function calls.
> I was sarcastic; "ray_sphere" calls that interection routine which
> cannot be inlined because it's virtual etc...
In my original implementation, "intersect" calls "ray_sphere" and not the
other way around. Do you mean "ray_sphere" could not have been inlined?
> Also, the fact is that this kind of sphere intersection can easily be
> made branchless with something like:
> if ((d < 0.) | (t2 < 0.)) return infinity;
> else return t1 > 0. ? t1 : t2;
> But you want to branch away from the square root as it is what's
> wasting most cycles.
I hadn't optimised "ray_sphere".
>> The C++ outperforms the OCaml on AMD64 (and IA64). I'd say that the
>> performance is comparable in all cases.
> Again, your benchmark as presented it totally unfair.
In what way is it unfair?
Here's the C++ I had:
#include <vector>
#include <iostream>
#include <limits>
#include <cmath>
using namespace std;
numeric_limits<double> dbl;
double delta = sqrt(dbl.epsilon()), infinity = dbl.infinity(), pi = M_PI;
struct Vec { // 3D vector
double x, y, z;
Vec(double x2, double y2, double z2) : x(x2), y(y2), z(z2) {}
};
inline Vec operator+(const Vec &a, const Vec &b)
{ return Vec(a.x + b.x, a.y + b.y, a.z + b.z); }
inline Vec operator-(const Vec &a, const Vec &b)
{ return Vec(a.x - b.x, a.y - b.y, a.z - b.z); }
inline Vec operator*(double a, const Vec &b)
{ return Vec(a * b.x, a * b.y, a * b.z); }
inline double dot(const Vec &a, const Vec &b) { return a.x*b.x + a.y*b.y +
a.z*b.z; }
inline Vec unitise(const Vec &a) { return (1 / sqrt(dot(a, a))) * a; }
struct Ray { Vec orig, dir; Ray(Vec o, Vec d) : orig(o), dir(d) {} };
struct Tree {
vector<Tree> child;
Vec center;
double radius;
Tree(Vec c, double r) : child(), center(c), radius(r) {}
};
inline double ray_sphere(const Ray &ray, const Tree &tree) {
Vec v = tree.center - ray.orig;
double b = dot(v, ray.dir), disc = b*b - dot(v, v) +
tree.radius*tree.radius;
if (disc < 0) return infinity;
double d = sqrt(disc), t2 = b + d;
if (t2 < 0) return infinity;
double t1 = b - d;
return (t1 > 0 ? t1 : t2);
}
void intersect(double &lambda, Vec &normal, const Ray &ray, const Tree
&tree) {
double l = ray_sphere(ray, tree);
if (l >= lambda) return;
if (tree.child.size() == 0) {
lambda = l;
normal = unitise(ray.orig + l * ray.dir - tree.center);
} else
for (vector<Tree>::const_iterator it=tree.child.begin();
it!=tree.child.end(); ++it)
intersect(lambda, normal, ray, *it);
}
bool intersect(const Ray &ray, const Tree &tree) {
if (ray_sphere(ray, tree) == infinity) return false;
if (tree.child.size() == 0) return true; else {
for (vector<Tree>::const_iterator it=tree.child.begin();
it != tree.child.end(); ++it)
if (intersect(ray, *it)) return true;
}
return false;
}
double ray_trace(const double weight, const Vec &light, const Ray &ray,
const Tree &tree) {
double lambda = infinity;
Vec normal(0, 0, 0);
intersect(lambda, normal, ray, tree);
if (lambda == infinity) return 0;
Vec o = ray.orig + lambda * ray.dir + delta * normal;
double g = -dot(normal, light);
if (g <= 0) return 0.;
return (intersect(Ray(o, Vec(0, 0, 0) - light), tree) ? 0 : g);
}
Tree create(int level, double r, double x, double y, double z) {
if (level == 1) return Tree(Vec(x, y, z), r);
Tree group = Tree(Vec(x, y, z), 3*r);
group.child.push_back(Tree(Vec(x, y, z), r));
double rn = 3*r/sqrt(12.);
for (int dz=-1; dz<=1; dz+=2)
for (int dx=-1; dx<=1; dx+=2)
group.child.push_back(create(level-1, r/2, x-dx*rn, y+rn, z-dz*rn));
return group;
}
int main(int argc, char *argv[]) {
int level = (argc==2 ? atoi(argv[1]) : 6), w = 512, h = 512, ss = 4;
Tree tree = create(level, 1, 0, -1, 0);
cout << "P2\n" << w << " " << h << "\n256\n";
for (int y=h-1; y>=0; --y) {
for (int x=0; x<w; ++x) {
double g=0;
for (int dx=0; dx<ss; ++dx)
for (int dy=0; dy<ss; ++dy) {
Vec d(x+double(dx)/ss-w/2, y+double(dy)/ss-h/2, max(w, h));
g += ray_trace(1, unitise(Vec(-1, -3, 2)), Ray(Vec(0, 0, -4),
unitise(d)), tree);
}
cout << min(255, int(256. * g / (ss*ss))) << " ";
}
cout << endl;
}
return 0;
tbp
> I appreciate it, thank you.
> On my 1.2GHz Athlon t-bird I get:
> $ g++-3.4 -Wall -O3 -ffast-math -funroll-all-loops miniray.cpp -o
miniray
> miniray.cpp:18: warning: division by zero in `1.0e+0 / 0.'
> $ time ./miniray >image.ppm
>
> real 0m28.587s
> user 0m28.400s
> sys 0m0.040s
>
> So your optimised version is 1.9x faster.
I've only quickly checked with g++-3.4, g++-4.x does a much better job
(at least on my opteron on x86-64, haven't checked on my t-bird).
And i wouldn't consider my version as optimized by any stretch, but
just one way to decently write a c++ implementation of that thing.
> My optimised version is basically the same speed as yours (although
quite
> different):
Hmm, need to check what it gives with a shadow ray shortcut, if i can
cram it up in there.
> > c) i had to raise the epsilon a bit to avoid acnea, it has to do
with
> > the way i compute normals but i'm not going to fix that.
>
> Your epsilon is actually smaller than my delta (which was sqrt
epsilon).
Ah! Missed that bit. Anyway that shouldn't make a difference as both
epsilons are too small to cull any part of the hierarchy away.
> > d) i removed the virtualness in the most straightforward way and
it's
> > certainly not aesthetically pleasing
>
> I'll study this in more detail.
>
> > e) i'm not going to spend ages on this, hence c) & d).
>
> What better solution would you recommend for (d)?
If that's better in the more efficient sense, none.
If that's better in the tigher memory foot print sense, you could
shrink it a bit by removing the use of a bounding sphere for 'final'
nodes (those sitting at the lowest lvl) and paying for that at runtime
with a branch or something.
Other than that it's just a matter of taste, and the virtual solution
is the most concise.
> My optimised version uses:
>
> struct Tree {
> vector<Tree> child;
> Vec center;
> double radius;
> Tree(Vec c, double r) : child(), center(c), radius(r) {}
> };
>
> which I don't much like.
It's a bit tighter than my solution, but both are variations on the
same theme.
Ideally i would use an explicit memory layout so things could be
properly aligned etc... And it would certainly end up being much
uglier.
> > f) it's algorithmically equivalent to your version (ie primary and
> > shadow rays aren't special cased etc...)
>
> My optimised version has a specialised intersect function for shadow
rays
> (but still fits in <100 LOC and is just a clear, IMHO).
I'll try to bolt something like that on top of my version.
> Can you describe how a new kind of object/bound would be implemented?
A BVH is fine for that kind of scenes, but it's all in the construction
of the hierarchy (or put another way, in its quality).
Now i seriously doubt i could cram a SAH compiler within 100 LOC :)
> In my original implementation, "intersect" calls "ray_sphere" and not
the
> other way around. Do you mean "ray_sphere" could not have been
inlined?
Sphere::intersect calls Sphere::ray_sphere, and that maybe inlined but
Sphere::intersect itself then cannot be inlined because it's virtual
and not statically resolved (and on top of that instead of a static
call, you pay for a indirect one).
Compare that to what happens in my version where the compiler has the
choice to inline to its heart content; that is placing (or not) a call
if that's cheaper at some point.
Like i've said earlier exposing things to the compiler is the name of
the game.
> > Also, the fact is that this kind of sphere intersection can easily
be
> > made branchless with something like:
> > if ((d < 0.) | (t2 < 0.)) return infinity;
> > else return t1 > 0. ? t1 : t2;
> > But you want to branch away from the square root as it is what's
> > wasting most cycles.
>
> I hadn't optimised "ray_sphere".
I was just saying that, in that function cycles are spent computing the
square root, not branching; so it's cheaper to avoid with a branch than
to unconditionnaly compute it (and then use some conditionnal moves).
> >> The C++ outperforms the OCaml on AMD64 (and IA64). I'd say that
the
> >> performance is comparable in all cases.
> > Again, your benchmark as presented it totally unfair.
>
> In what way is it unfair?
Because you obstructed the c++ compiler in many many ways, the most
prominent being passing structures by value instead of by reference (as
show in previous patches).
That's why i've also said the comparison was rigged as you didn't put
as much effort on the c++ side as on the OCaml; it's not surprising my
version is ~2x faster than the original.
Now if you present an amended version that at least fixes those issues,
i'll have no objections to whatever conclusion you can come up to.
Warp
> >> Creating structs for the return values of every function that returns
> >> multiple values may be a good idea, but the verbosity and clumsiness of
> >> doing this in C++ makes it quite undesirable.
> >
> > You mean modular?
> No.
Then IMO you are wrong.
> Do you see what I mean by verbose and clumsy?
No.
Granted, C++ does not offer some features other programming languages do,
but you gave a misleading example. You gave the impression that if you
want to return three values from a function you have to always write
the entire struct code plus all the separate assignments.
That is like saying that each time you want to use a string you have
to write a string class explicitly.
C++ offers tools to make that code of yours simpler and shorter.
The way you gave in your example is not the only nor by any means
the best way.
Besides, it usually comes down to design. It is rare that a function
just returns some values whose type and meaning are random. Usually
those returned values have a close relation to each other and form
a logical group of values.
For example, if you have a function which returns an intersection point
and a normal vector, does it make sense, from design point of view, to
return these values using generic nameless structures with no higher
logic in them?
Or does it make more sense to create a logical data container called
according to its usage, for example "Intersection", which contains the
two (or more) things which are closely related to each other?
Which do you consider more logical, a function which just returns
a pair of values (using a generic container), or a function which returns
a value of type "Intersection"?
--
- Warp
tbp
> > So what? The fact is that implementation quality matters just as
much.
>
> No, it doesn't. If you can reduce your algorithm from O(n) to O(log
n) then
> you should do so before you rewrite it in assembler.
it would pointless to argue further more.
> > A program has to run on a given machine at some point.
>
> True. I'm not saying that it is not worthwhile to optimise a ray
tracer with
> SIMD, just that it would be premature to do so on my ray tracer right
now.
I beg to differ. You can't simply plug in some SIMD scheme on top of
that as an afterthought, because it's tightly coupled with the way you
represent things.
That's one of the common pitfall of OO abstraction.
Wald discuss that point a bit in his thesis and i couldn't agree more.
But i know it's a not a very popular opinion atm.
> I would be very interested to see a highly optimised version of my
ray
> tracer. I think a lot of people could learn from a set of web pages
which
> show how to optimise the programs. If you're willing to help and I
can find
> the time then I'd like to do this sometime.
There's certainly lots of angles to attack that problem
(algorithmically or not), but i doubt any of them could fit within 100
LOC.
Beside C++ idiomatic issues, a better memory layout could certainly
help.
That's, perhaps, the lowest hanging fruit (that or 'cheating' with
lower numerical precision, ie for square roots).
Anyway as i've already invested some time in a c++ version, i'm open to
suggestions :)
> I wouldn't call this kind of stuff C++:
>
> a = _mm_load_ps(inp_sse1);
> b = _mm_load_ps(inp_sse2);
> c = _mm_add_ps(a, b);
> _mm_store_ps(out_sse, c);
But that's legit C++.
C++ is C on steroids while C itself is just a gloried assembler.
No matter how much OO or syntactic sugar you pour onto it.
I kinda apreciate to be able to code to-the-metal.
To each, its own.
Jon Harrop
> Besides, it usually comes down to design. It is rare that a function
> just returns some values whose type and meaning are random.
No, that is my point. It is common, not rare. Look at the prevelance of
tuple return types in ML programs and you'll see what I mean. You won't see
it in C++ code because C++ makes it too tedious for people to bother coding
it.
Many of the functions in the OCaml standard library return 2-tuples. None of
them need to be explicitly defined.
Jon Harrop
> You can't simply plug in some SIMD scheme on top of
> that as an afterthought, because it's tightly coupled with the way you
> represent things.
> That's one of the common pitfall of OO abstraction.
> Wald discuss that point a bit in his thesis and i couldn't agree more.
> But i know it's a not a very popular opinion atm.
I agree, actually. I know nothing of cache coherency in the context of ray
tracers though.
> There's certainly lots of angles to attack that problem
> (algorithmically or not), but i doubt any of them could fit within 100
> LOC.
Yes, it wouldn't make for a very good introductory tutorial either. Although
I'm sure a lot of people (me included) would be interested.
> Beside C++ idiomatic issues, a better memory layout could certainly
> help.
> That's, perhaps, the lowest hanging fruit (that or 'cheating' with
> lower numerical precision, ie for square roots).
It would be interesting to compare equivalently optimised C++ and OCaml
programs based upon float arrays. Handling alignment will be tricky in
OCaml.
>> I wouldn't call this kind of stuff C++:
>>
>> a = _mm_load_ps(inp_sse1);
>> b = _mm_load_ps(inp_sse2);
>> c = _mm_add_ps(a, b);
>> _mm_store_ps(out_sse, c);
> But that's legit C++.
Is that in the C++ standard or is it a common extension?
> C++ is C on steroids while C itself is just a gloried assembler.
> No matter how much OO or syntactic sugar you pour onto it.
> I kinda apreciate to be able to code to-the-metal.
> To each, its own.
True.
Jon Harrop
> I've only quickly checked with g++-3.4, g++-4.x does a much better job
> (at least on my opteron on x86-64, haven't checked on my t-bird).
I'll check it out when the Debian packages arrive. :-)
> And i wouldn't consider my version as optimized by any stretch, but
> just one way to decently write a c++ implementation of that thing.
"more optimised" then. :-)
>> My optimised version is basically the same speed as yours (although
> quite
>> different):
> Hmm, need to check what it gives with a shadow ray shortcut, if i can
> cram it up in there.
I don't think it actually made mine much faster so it probably isn't worth
it.
> It's a bit tighter than my solution, but both are variations on the
> same theme.
> Ideally i would use an explicit memory layout so things could be
> properly aligned etc... And it would certainly end up being much
> uglier.
The OCaml might be significantly easier to understand if suitable
higher-order functions can be factored out which make the memory layout
"look" like a tree.
>> Can you describe how a new kind of object/bound would be implemented?
> A BVH is fine for that kind of scenes, but it's all in the construction
> of the hierarchy (or put another way, in its quality).
> Now i seriously doubt i could cram a SAH compiler within 100 LOC :)
What is an SAH compiler?
>> In my original implementation, "intersect" calls "ray_sphere" and not
> the
>> other way around. Do you mean "ray_sphere" could not have been
> inlined?
> Sphere::intersect calls Sphere::ray_sphere, and that maybe inlined but
That's the important one though.
> Sphere::intersect itself then cannot be inlined because it's virtual
> and not statically resolved (and on top of that instead of a static
> call, you pay for a indirect one).
> Compare that to what happens in my version where the compiler has the
> choice to inline to its heart content; that is placing (or not) a call
> if that's cheaper at some point.
> Like i've said earlier exposing things to the compiler is the name of
> the game.
Yes. Although this is only useful to other people if it has a description of
what information is being provided to the compiler. It would be nice to
write this up.
>> In what way is it unfair?
> Because you obstructed the c++ compiler in many many ways,
I think the OCaml compiler was equivalently obstructed.
> the most
> prominent being passing structures by value instead of by reference (as
> show in previous patches).
I think the compiler could (should) have spotted that optimisation without
having to be told.
> That's why i've also said the comparison was rigged as you didn't put
> as much effort on the c++ side as on the OCaml; it's not surprising my
> version is ~2x faster than the original.
Well, I have a more optimised OCaml version which is just as fast on x86 and
20 lines shorter.
> Now if you present an amended version that at least fixes those issues,
> i'll have no objections to whatever conclusion you can come up to.
I'll try to write another page describing both optimised versions. I also
have a C version now, which is slower than both the C++ and the OCaml. :-)
tbp
> tbp wrote:
> > I've only quickly checked with g++-3.4, g++-4.x does a much better
job
> > (at least on my opteron on x86-64, haven't checked on my t-bird).
>
> I'll check it out when the Debian packages arrive. :-)
snapshot is straightforward, yet a bit off topic :)
> > Hmm, need to check what it gives with a shadow ray shortcut, if i
can
> > cram it up in there.
>
> I don't think it actually made mine much faster so it probably isn't
worth
> it.
A quick hack at that (within the rules) shaved a second of runtime
here.
# diff -u miniray.cpp miniray_shadow.cpp
--- miniray.cpp 2005-05-11 14:09:08.000000000 +0200
+++ miniray_shadow.cpp 2005-05-11 20:11:00.000000000 +0200
@@ -47,12 +47,13 @@
if (t >= hit.t) return;
else { hit.t = t; hit.n = leaf.get_normal(ray.o + ray.d*t); }
}
- void intersect(const ray_t &ray, hit_t &hit) const {
+ template<bool shadow> void intersect(const ray_t &ray, hit_t &hit)
const {
if (!objs.empty()) {
if (bound.intersect(ray) >= hit.t) return;
intersect_leaf(ray,hit);
+ if (shadow & (hit.t != infinity)) return; // early exit.
for (objs_t::const_iterator it=objs.begin(); it!=objs.end(); ++it)
- (*it)->intersect(ray, hit);
+ (*it)->intersect<shadow>(ray, hit);
}
else intersect_leaf(ray,hit);
}
@@ -68,11 +69,11 @@
else return new group_t(sphere_t(v, r),sphere_t(v, r));
}
static double ray_trace(const group_t * const scene, const ray_t &ray)
{
- hit_t hit; scene->intersect(ray, hit);
+ hit_t hit; scene->intersect<false>(ray, hit);
const double diffuse = (hit.t == infinity) ? 0. : -hit.n.dot(light);
if (diffuse <= 0.) return 0.;
const ray_t sray(ray.o + (ray.d*hit.t) + (hit.n*epsilon), -light);
- hit_t shit; scene->intersect(sray,shit);
+ hit_t shit; scene->intersect<true>(sray,shit);
return (shit.t == infinity) ? diffuse : 0.;
}
int main(int argc, char *argv[]) {
# time ./miniray >pix.ppm
real 0m4.856s
user 0m4.812s
sys 0m0.004s
> The OCaml might be significantly easier to understand if suitable
> higher-order functions can be factored out which make the memory
layout
> "look" like a tree.
I'm absolutely not qualified to comment on OCaml unfortunately :)
> >> Can you describe how a new kind of object/bound would be
implemented?
> > A BVH is fine for that kind of scenes, but it's all in the
construction
> > of the hierarchy (or put another way, in its quality).
> > Now i seriously doubt i could cram a SAH compiler within 100 LOC :)
>
> What is an SAH compiler?
Surface Area Heuristic, a way to evaluate costs when building a
hierarchy.
Does wonders on BVH, kd-tree... you name it. But it's tricky to get it
right and even more to make it not behave quadratically.
Havran's thesis is a must read, Wald builds upon it (but there's a
metric ton of litterature about it).
http://www.cgg.cvut.cz/~havran/phdthesis.html
http://www.mpi-sb.mpg.de/~wald/PhD/
> > Like i've said earlier exposing things to the compiler is the name
of
> > the game.
>
> Yes. Although this is only useful to other people if it has a
description of
> what information is being provided to the compiler. It would be nice
to
> write this up.
I'm pedagogically challenged :)
On top of peremptoric statements like "know the language" and "know
your compiler", it boils down to things like const correctness,
tediously constraining things and so on.
Or put another way, compilers are mental midgets and C++ is a language
made for compilers and not humans.
Ok, that wasn't really helpful.
But i'm sure someone already has put that in a nicer way somewhere, no?
:)
> >> In what way is it unfair?
> > Because you obstructed the c++ compiler in many many ways,
>
> I think the OCaml compiler was equivalently obstructed.
I can't say. But for sure, your C++ "style" was not typical.
> > the most
> > prominent being passing structures by value instead of by reference
(as
> > show in previous patches).
>
> I think the compiler could (should) have spotted that optimisation
without
> having to be told.
That's wishful thinking. And it really really doesn't work that way.
And don't blame g++ for that, as it's actually quite good at that game.
But that's another topic.
> > That's why i've also said the comparison was rigged as you didn't
put
> > as much effort on the c++ side as on the OCaml; it's not surprising
my
> > version is ~2x faster than the original.
>
> Well, I have a more optimised OCaml version which is just as fast on
x86 and
> 20 lines shorter.
Bring it on ;)
> > Now if you present an amended version that at least fixes those
issues,
> > i'll have no objections to whatever conclusion you can come up to.
>
> I'll try to write another page describing both optimised versions. I
also
> have a C version now, which is slower than both the C++ and the
OCaml. :-)
Hah.
tbp
> I agree, actually. I know nothing of cache coherency in the context
of ray
> tracers though.
about mem/cache coherency, but also _rays_ coherency (dunno if that was
clear).
Then cache coherency isn't an issue in your app, at least at the
default level.
Though, better alignement would certainly help wrt doubles on x86.
> > There's certainly lots of angles to attack that problem
> > (algorithmically or not), but i doubt any of them could fit within
100
> > LOC.
>
> Yes, it wouldn't make for a very good introductory tutorial either.
Although
> I'm sure a lot of people (me included) would be interested.
Jacco Bikker has produced some interesting articles on flipcode about
that.
But his articles are now lagging way behind his own implementation.
He's a lazy bastard :)
I've implemented, both a BVH & kdtree based coherent raytracer, but
even the BVH one (which is simpler) would require lots of background
description before addressing the meaty stuff.
That's quite infortunate (and the litterature on the topic is
absolutely horrible to read).
> > Beside C++ idiomatic issues, a better memory layout could certainly
> > help.
> > That's, perhaps, the lowest hanging fruit (that or 'cheating' with
> > lower numerical precision, ie for square roots).
>
> It would be interesting to compare equivalently optimised C++ and
OCaml
> programs based upon float arrays. Handling alignment will be tricky
in
> OCaml.
Float arrays? I fail to grok the context. Could you elaborate a bit?
> >> I wouldn't call this kind of stuff C++:
> >>
> >> a = _mm_load_ps(inp_sse1);
> >> b = _mm_load_ps(inp_sse2);
> >> c = _mm_add_ps(a, b);
> >> _mm_store_ps(out_sse, c);
> > But that's legit C++.
>
> Is that in the C++ standard or is it a common extension?
Common x86 extension yes (seen in gcc, msvc, icc etc).
An altivec equivalent is common on PPC platforms, i think.
Admitedly, it's in phase with the "C as a glorified assembler" i
alluded to earlier.
Warp
> No, that is my point. It is common, not rare.
So effectively crappy design is common, not rare. However, we all knew
that already, of course.
--
- Warp
tbp
tbp wrote:
> A quick hack at that (within the rules) shaved a second of runtime
> here.
> # time ./miniray >pix.ppm
>
> real 0m4.856s
> user 0m4.812s
> sys 0m0.004s
I'm replying to myself (don't you love that), because while tinkering
with it, again - that's quite addictive -, i've noticed 2 things:
a) the original deallocation routine, Group::del(), that i've
shamelessly re-used is bogus and only delete the first level.
b) more importantly, the BVH can be tighten a bit: at a given level, if
you place a sphere of radius r centered at v(x,y,z), then instead of
bounding the sub hierarchy with sphere(v, 3*r) you can bound it with
sphere(v+(0,r,0), 2.4330128*r).
Then,
# time ./miniray >pix.ppm
real 0m3.734s
user 0m3.721s
sys 0m0.003s
Voilà.
Jon Harrop
> tbp wrote:
>> A quick hack at that (within the rules) shaved a second of runtime
>> here.
> [snip]
>> # time ./miniray >pix.ppm
>>
>> real 0m4.856s
>> user 0m4.812s
>> sys 0m0.004s
>
> I'm replying to myself (don't you love that), because while tinkering
> with it, again - that's quite addictive -, i've noticed 2 things:
> a) the original deallocation routine, Group::del(), that i've
> shamelessly re-used is bogus and only delete the first level.
Oops, yes. How thoughtless of you to copy my mistake. :-)
> b) more importantly, the BVH can be tighten a bit: at a given level, if
> you place a sphere of radius r centered at v(x,y,z), then instead of
> bounding the sub hierarchy with sphere(v, 3*r) you can bound it with
> sphere(v+(0,r,0), 2.4330128*r).
Yes, I hadn't optimised that either (just copied the idea from the real
sphere-flake, for which the easy approach is more excusable). :-)
> Then,
> # time ./miniray >pix.ppm
> real 0m3.734s
> user 0m3.721s
> sys 0m0.003s
>
> Voilą.
I'm surprised it makes such a big difference. I shall check it out in more
detail. I really have spent too much time playing with this now. :-)
tbp
> Oops, yes. How thoughtless of you to copy my mistake. :-)
> I'm surprised it makes such a big difference. I shall check it out in
more
> detail. I really have spent too much time playing with this now. :-)
All the magic is in the hierarchy.
Now as i couldn't see how i could pack in the code for a better bvh and
this tree is complete i've addressed other points this morning:
. i turn the tree into an array with skip pointers - that can be done
in 1 pass, the tree is complete - while i build in place.
. then i can iterate over it when tracing.
[uglified 100 LOC source annexed]
[also here http://ompf.org/ray/bvh100.cpp]
So i was at
> > # time ./miniray >pix.ppm
> > real 0m3.734s
> > user 0m3.721s
> > sys 0m0.003s
and now
# g++ -O2 -ffast-math -o bvh bvh100.cpp
# time ./bvh >pix.ppm
real 0m2.985s
user 0m2.966s
sys 0m0.005s
but
# cmp pix.orig.ppm pix.ppm
pix.orig.ppm pix.ppm differ: byte 1040, line 4
even if a diff between the two results in a black picture for the eye.
It's not as fast as it could, but the source is stuffed enough as it is
(cough).
// sphere bvh raytracer v0.3, (c) thierry berger-perrin 2005.
// keywords: sphere, bvh, iterative, array form, skip pointers.
#include <cmath>
#include <iostream>
#define EXPECT(a) __builtin_expect((a),true) // (a)
struct vec_t {
double x,y,z; vec_t() {}
vec_t(double a,double b,double c): x(a),y(b),z(c) {}
vec_t operator+(const vec_t &v) const { return
vec_t(x+v.x,y+v.y,z+v.z); }
vec_t operator-(const vec_t &v) const { return
vec_t(x-v.x,y-v.y,z-v.z); }
vec_t operator-() const { return vec_t(-x,-y,-z); }
vec_t operator*(const double d) const { return vec_t(x*d,y*d,z*d); }
double dot(const vec_t &v) const { return x*v.x + y*v.y + z*v.z; }
double magsqr() const { return dot(*this); }
vec_t norm() const { return *this*(1./sqrt(magsqr())); } };
static const vec_t light(vec_t(-1, -3, 2).norm());
static const double infinity = 1./0, epsilon = 1e-12;
struct ray_t {
vec_t o, d; ray_t(const vec_t &v) : o(v) {}
ray_t(const vec_t &v, const vec_t &w) : o(v),d(w) {} };
struct hit_t { vec_t n; double t; hit_t() : n(vec_t(0,0,0)),
t(infinity) {} };
struct sphere_t { vec_t o; double r; sphere_t() {}
sphere_t(const vec_t &v, double d) : o(v),r(d) {} // 1/r, constant
vec_t get_normal(const vec_t &v) const { return (v-o)*(1./r); }
double intersect(const ray_t &ray) const {
const vec_t v(o-ray.o);
const double b = ray.d.dot(v), disc = b*b - v.magsqr() + r*r;
if (disc < 0.) return infinity; // branch away from the square root
const double d = sqrt(disc), t2 = b + d, t1 = b - d;
if (t2 < 0.) return infinity; else return (t1 > 0. ? t1 : t2); //
cond mv.
} };
struct node_t; static node_t *pool = 0, *end = 0;
struct node_t { // used in array form + skip for navigation.
sphere_t bound, leaf; long diff; // gah, 2*32 + some. diff should be
moved.
node_t() {}
node_t(const sphere_t &b,const sphere_t &l, const int jump)
: bound(b), leaf(l), diff(jump) {}
template<bool shadow> static void intersect(const ray_t &ray, hit_t
&hit) {
const node_t * __restrict__ p = pool; //const node_t * p = pool;
while (p < end) {
if (EXPECT(p->bound.intersect(ray) >= hit.t)) // missed bound
p = (p->diff > 0) ? p+p->diff : end; // skip
else {
const double t = p->leaf.intersect(ray);
if (t < hit.t) { // if hit, update, then break for shadows.
hit.t = t; if (shadow) break;
hit.n = p->leaf.get_normal(ray.o + ray.d*t);
}
++p; // next
} } } };
static double ray_trace(const node_t * const scene, const ray_t &ray) {
hit_t hit; scene->intersect<false>(ray, hit); // trace primary
const double diffuse = (hit.t == infinity) ? 0. : -hit.n.dot(light);
if (diffuse <= 0.) return 0.;
const ray_t sray(ray.o + (ray.d*hit.t) + (hit.n*epsilon), -light);
hit_t shit; scene->intersect<true>(sray,shit); // trace shadow
return (shit.t == infinity) ? diffuse : 0.; }
static node_t *create(node_t *node, const int level, int dist,
const double r, const vec_t &v)
{ if (level > 1) { // directly write the tree as an array via placement
new.
node = 1+new (node) node_t(
sphere_t(v+vec_t(0,r,0),2.4330128*r),sphere_t(v, r), dist+1);
dist = std::max((dist-4)/4, 0); // *(distance>1?4:0) + 4;
const double rn = (3/sqrt(12.))*r;
for (int dz=-1; dz<=1; dz+=2) { for (int dx=-1; dx<=1; dx+=2) { //
go down
node = create(node, level-1,
dist,0.5*r,v+vec_t(-dx*rn,+rn,-dz*rn));
} } return node;
} else return 1+new (node) node_t(sphere_t(v, r),sphere_t(v, r),1);
// final
}
template <int ss> static void trace_grid_aa(const int width, const int
height) {
ray_t ray(vec_t(0, 0, -4)); enum { ss_sqr = ss*ss };
const double w = width, h = height, scale = 256. / double(ss_sqr);
vec_t grid[ss*ss]; // pre compute a ss² regular grid perturbation
{ const double rcp_ss = 1./double(ss), halfw = w/2, halfh = h/2;
int idx = 0;
for (int dx=0; dx<ss; ++dx) for (int dy=0; dy<ss; ++dy) {
grid[idx] = vec_t(double(dx)*rcp_ss-halfw,
double(dy)*rcp_ss-halfh, 0);
++idx;
} } // ok, now let's get back to work
vec_t v(0,w-1, std::max(w,h)); // scan line
for (int i=height-1; i>=0; --i) { for (int j=width-1; j>=0; --j) {
double g=0; for (int idx = 0; idx < ss_sqr; ++idx) { //
apply grid
ray.d = (v+grid[idx]).norm(); g += ray_trace(pool, ray); }//
trace
std::cout << int(scale * g) << " "; v.x += 1; // next pixel
} v.x = 0; v.y -= 1; } // next line
std::cout << std::endl;
}
int main(int argc, char *argv[]) {
enum { w = 512, h = 512, ss = 4, alignment = 64 };
const int level = (argc==2 ? atoi(argv[1]) : 6);
int count = 4, dec = level; while (--dec > 1) count = (count*4) + 4;
pool = new node_t[count]; end = pool + count; // raw
create(pool, level, count, 1., vec_t(0, -1, 0)); // cooked
std::cout << "P2\n" << w << " " << h << "\n256\n";
trace_grid_aa<ss>(w,h); return 0; // served
}
Jon Harrop
> All the magic is in the hierarchy.
> Now as i couldn't see how i could pack in the code for a better bvh and
> this tree is complete i've addressed other points this morning:
What data structures and algorithms would you recommend? Perhaps they could
be implemented in the OCaml version...
> . i turn the tree into an array with skip pointers - that can be done
> in 1 pass, the tree is complete - while i build in place.
> . then i can iterate over it when tracing.
>
> > # time ./miniray >pix.ppm
> > real 0m3.734s
> > user 0m3.721s
> > sys 0m0.003s
> and now
> # g++ -O2 -ffast-math -o bvh bvh100.cpp
> # time ./bvh >pix.ppm
> real 0m2.985s
> user 0m2.966s
> sys 0m0.005s
On 1.2GHz Athlon t-bird, the unoptimised C++ on my webpage took:
$ g++-3.4 -Wall -O2 -ffast-math ray.cpp -o ray
$ time ./ray >image.pgm
real 1m51.071s
user 1m31.510s
sys 0m0.180s
and the unoptimised OCaml took:
$ ocamlopt -inline 100 -ffast-math ray11.ml -o ray11
$ time ./ray10 >image.pgm
real 0m53.916s
user 0m40.770s
sys 0m0.120s
Your latest version takes:
$ g++-3.4 -Wall -O2 -ffast-math bvh100.cpp -o bvh100
bvh100.cpp:17: warning: division by zero in `1.0e+0 / 0'
$ time ./bvh100 >image.pgm
real 0m18.492s
user 0m16.030s
sys 0m0.070s
> but
> # cmp pix.orig.ppm pix.ppm
> pix.orig.ppm pix.ppm differ: byte 1040, line 4
> even if a diff between the two results in a black picture for the eye.
A couple of points here:
1. I think the -ffast-math can cause gcc to violate IEEE-complaint FP
arithmetic, in which case it could theoretically alter the output and could
be banned as a consequence.
2. If we're happy to produce similar-looking results then switching to
single-precision takes over 40% off the time (on my T-bird, probably not on
AMD64):
$ g++-3.4 -Wall -O2 -ffast-math bvh100_2.cpp -o bvh100_2
bvh100_2.cpp:17: warning: division by zero in `1.0e+0 / 0'
$ time ./bvh100_2 >image2.pgm
real 0m11.282s
user 0m9.380s
sys 0m0.050s
For the shootout, the results will probably have to be identical. I'm more
than happy to generate correct-looking results (as this is the point of a
ray tracer!) for my web page though.
tbp
> > All the magic is in the hierarchy.
>
could
> be implemented in the OCaml version...
Forgetting about implementation complexity for a moment, there's 2
angles:
a) the nature of the hierarchy; could use a kd tree (tiny nodes, cheap
intersection, efficient culling) or a box bvh (maybe tighter); both can
be cheaply traversed iteratively.
b) regardless of the hierarchy type, use a cost driven scheme to build
it; ultimately the full SAH thing.
That would give us a better hierarchy, but then it will certainly take
more time to build it than to render... (ie my best SAH kdtree compiler
takes ~30s to compile a 250k triangle soup).
I think my iterative version could be ammended for better memory
alignment (at least on x86) just by splitting the node_t struct into 2
spheres on a side (2*32bytes) and skip pointers on the other. Would
lend to some uglyfication. Haven't checked what it would give with
floats.
On top of all that, you can also think of shooting more than 1 ray at
once and things like that ;)
> On 1.2GHz Athlon t-bird, the unoptimised C++ on my webpage took:
> $ g++-3.4 -Wall -O2 -ffast-math bvh100.cpp -o bvh100
> bvh100.cpp:17: warning: division by zero in `1.0e+0 / 0'
> $ time ./bvh100 >image.pgm
>
> real 0m18.492s
> user 0m16.030s
> sys 0m0.070s
Something's wrong. That's >6x slower than on my box, an opteron 146
(2ghz), which is supposedly only, say, 2x faster than yours.
But i'm using a gcc snapshot in 64bit mode:
a) there's more registers.
b) more importantly, SSE math is used (not an option with doubles on
your t-bird).
> A couple of points here:
>
> 1. I think the -ffast-math can cause gcc to violate IEEE-complaint FP
> arithmetic, in which case it could theoretically alter the output and
could
> be banned as a consequence.
There's much more "violation" when using the x86 fpu (because of the
internal 80bit representation etc) than SSE.
> 2. If we're happy to produce similar-looking results then switching
to
> single-precision takes over 40% off the time (on my T-bird, probably
not on
> AMD64):
It would help with the memory bandwidth. Other than that, only the sqrt
would be a bit faster.
Try changing the fpu precision on your t-bird.
> For the shootout, the results will probably have to be identical. I'm
more
Yes.
> than happy to generate correct-looking results (as this is the point
of a
> ray tracer!) for my web page though.
Yesterday i've done some experimentation to see if i could generate a
nicer scene (with that bvh in array yada yada).
Got:
http://ompf.org/ray/boo8.png
http://ompf.org/ray/ohh12.png
http://ompf.org/ray/blossom12.png
and finally, http://ompf.org/ray/holy_cabbage_of_antioch2.png
That final picture took 6.8s with a 2x2 RGSS.
I haven't tried to compact the source to 100 LOC yet, but i think it's
possible :)
Warp
> 1. I think the -ffast-math can cause gcc to violate IEEE-complaint FP
> arithmetic, in which case it could theoretically alter the output and could
> be banned as a consequence.
I have my doubts that it will cause any difference in the final image
(given that you round numbers correctly). Those optimizations, if they
cause any computation to give a different result, probably cause the
difference in the 15th decimal or so.
I don't know what kind of optimizations -ffast-math does, but one
that I could imagine is this kind:
double d = 0.1234;
double a = x/d, b = y/d, c = z/d;
That could be probably optimized internally to:
double d = 1/0.1234;
double a = x*d, b = y*d, c = z*d;
Due to how floating point works a, b and c might get just a slightly
different value (the difference would probably be in the least significant
decimal or the two least, which with doubles is something like the 15th
decimal or so).
--
- Warp
Francois Labreque
>
> and the unoptimised OCaml took:
>
> $ ocamlopt -inline 100 -ffast-math ray11.ml -o ray11
> $ time ./ray10 >image.pgm
>
> real 0m53.916s
> user 0m40.770s
> sys 0m0.120s
>
Is there a reason why you're compiling "ray11" and then running "ray10"
for the test?
--
Francois Labreque | The surest sign of the existence of extra-
flabreque | terrestrial intelligence is that they never
@ | bothered to come down here and visit us!
videotron.ca | - Calvin
Jon Harrop
> Jon Harrop wrote:
>>
>> and the unoptimised OCaml took:
>>
>> $ ocamlopt -inline 100 -ffast-math ray11.ml -o ray11
>> $ time ./ray10 >image.pgm
>>
>> real 0m53.916s
>> user 0m40.770s
>> sys 0m0.120s
>>
>
> Is there a reason why you're compiling "ray11" and then running "ray10"
> for the test?
Oops. :-)
Jon Harrop
> On top of all that, you can also think of shooting more than 1 ray at
> once and things like that ;)
I might try this if I can find the time.
>> On 1.2GHz Athlon t-bird, the unoptimised C++ on my webpage took:
>> $ g++-3.4 -Wall -O2 -ffast-math bvh100.cpp -o bvh100
>> bvh100.cpp:17: warning: division by zero in `1.0e+0 / 0'
>> $ time ./bvh100 >image.pgm
>>
>> real 0m18.492s
>> user 0m16.030s
>> sys 0m0.070s
> Something's wrong. That's >6x slower than on my box, an opteron 146
> (2ghz), which is supposedly only, say, 2x faster than yours.
I think gcc is suffering from poor FP code generation. Interestingly, it is
usually ocamlopt that suffers from this and not gcc. I've no idea why but
this ray tracer is the only program I've ever written which produces these
kind of results (faster OCaml on x86).
Someone has kindly ported it to Fortran on the shootout now. For the same
functionality, the Fortran is much faster than both the C++ and the OCaml
(and the SML compiled with Mlton, which is actually significantly faster
than C++).
> But i'm using a gcc snapshot in 64bit mode:
> a) there's more registers.
> b) more importantly, SSE math is used (not an option with doubles on
> your t-bird).
My AMD64 is also 6x faster (it is supposed to be 3x faster, according to
other benchmarks).
>> A couple of points here:
>>
>> 1. I think the -ffast-math can cause gcc to violate IEEE-complaint FP
>> arithmetic, in which case it could theoretically alter the output and
> could
>> be banned as a consequence.
> There's much more "violation" when using the x86 fpu (because of the
> internal 80bit representation etc) than SSE.
True.
>> 2. If we're happy to produce similar-looking results then switching
> to
>> single-precision takes over 40% off the time (on my T-bird, probably
> not on
>> AMD64):
> It would help with the memory bandwidth. Other than that, only the sqrt
> would be a bit faster.
> Try changing the fpu precision on your t-bird.
It doubles the performance. However, it almost halves the performance on
AMD64.
tbp
> I think gcc is suffering from poor FP code generation. Interestingly, it is
> usually ocamlopt that suffers from this and not gcc. I've no idea why but
> this ray tracer is the only program I've ever written which produces these
> kind of results (faster OCaml on x86).
I guess it's hitting some weakspot in gcc; that might be interesting to
report.
It's interesting to note that on x86-64 ICC 8.1 & 9.0beta can't beat
gcc (or more precisely a 4.1 snapshot).
> Someone has kindly ported it to Fortran on the shootout now. For the same
> functionality, the Fortran is much faster than both the C++ and the OCaml
> (and the SML compiled with Mlton, which is actually significantly faster
> than C++).
I guess you really mean faster on x86 (as opposed to x86-64), and then
i would put the blame on the 32bit fpu codegen thing too.
I mean C/C++ being as low level as it is, it should be faster given a
proper implementation; uglier, hackish, verbose etc but not slower :)
> > Try changing the fpu precision on your t-bird.
>
> It doubles the performance. However, it almost halves the performance on
> AMD64.
Huh?
On x86-64, SSE should be used and that op. set is impervious to FPU
precision modes.
If you're talking about swapping doubles for floats, i've checked that
and i get a marginal speed improvement (after fixing stuff & memory
layout).
So, can you elaborate a bit?
Last week i've cleaned things a bit and put stub here,
http://ompf.org/ray/sphereflake/, for a 100 LOC C++ version that
generates a sphere flake, or something looking closer to the real
thing.
I've settled for a 2x2 rotated grid to get decent runtime while
preserving image quality.
Regards,
tbp.
Darren New
>>Someone has kindly ported it to Fortran on the shootout now.
> I mean C/C++ being as low level as it is, it should be faster given a
> proper implementation; uglier, hackish, verbose etc but not slower :)
FORTRAN has various language rules (such as disallowing aliasing,
allowing better optimizations) that C and C++ don't have. (FORTRAN
originally was also very restricted in the forms that array indicies
could take, ensuring the compiler could make use of index registers, for
example.) So, no, FORTRAN is often faster than C for the same program.
--
Darren New / San Diego, CA, USA (PST)
The samba was clearly inspired
by the margarita.
tbp
> FORTRAN has various language rules (such as disallowing aliasing,
> allowing better optimizations) that C and C++ don't have. (FORTRAN
> originally was also very restricted in the forms that array indicies
> could take, ensuring the compiler could make use of index registers, for
> example.) So, no, FORTRAN is often faster than C for the same program.
of tightening aliasing rules as compiler options (with said compilers
being smarter & smarter about aliasing); no silver bullet but it shows
some progress on that side.
Anyway, if you said C/C++ was ugly, borked, byzantine or unproductive i
would have agreed. But when it comes to runtime performance, i simply
can't buy it: you can't get much lower than that without writing
straight asm (with the corrolary that you may have to express your
problem in compiler friendly way)... it all boils down to the
implementation's quality.
PS: i really didn't meant to start yet another sterile religious debate
about languages :)
Jon Harrop
> tbp wrote:
>>>Someone has kindly ported it to Fortran on the shootout now.
>
>> I mean C/C++ being as low level as it is, it should be faster given a
>> proper implementation; uglier, hackish, verbose etc but not slower :)
>
> FORTRAN has various language rules (such as disallowing aliasing,
> allowing better optimizations) that C and C++ don't have. (FORTRAN
> originally was also very restricted in the forms that array indicies
> could take, ensuring the compiler could make use of index registers, for
> example.) So, no, FORTRAN is often faster than C for the same program.
Not only is that untrue, it hasn't been true for many years. The main reason
that non-trivial Fortran programs are typically much slower than their
modern counterparts is that the verbosity and difficulty of implementing
more data structures and algorithms in Fortran makes non-trivial approaches
intractable.
Consider to Fortran port of my ray tracer. Not only is the resulting program
3x as long as the OCaml (and by far the longest implementation), not as
general (the scene tree is prohibitively difficult to implement in Fortran
without sacrificing its current genericity) but the resulting program is
not significantly faster than implementations written in modern languages
(slightly faster than OCaml, slightly slower than Mlton-compiled SML).
As TBP has said, this discussion doesn't belong here. If you want to know
more about this then read my book on OCaml for scientists:
http://www.ffconsultancy.com/products/ocaml_for_scientists/