"virtualmeet" <tah### [at] yahoofr> wrote:

> Example: Cos(x) + cos(y) + cos(z)  (Schwartz). Well, beleave it or not but a
> grid of 1000x1000x1000 will require only 3000 cosinus fct calculations (I'm
> already able to use this kind of simplifications), rather than the normal 10^9.
Correction: rather than the normal 3x10^9.
PS: This works not only for "cos" but for all kind of unary fcts, ie:sin,
atan, log, "pow" fct with a constant...with the same ratio of calculs
simplifications (1/10^6 for a grid of 10^3x10^3x10^3).
Taha

Post a reply to this message

virtualmeet wrote:
> grid of 1000x1000x1000 

You do realize that such a grid would also require eight gigabyte (!!!) of
memory? Using the fastest personal computer memory available, using a 128
bit wide bus, you still need more than a second to even fill that grid with
any single identical value.

Assuming you render a 1280 * 1028 pixel image, in the time it takes to fill
the memory of your grid, you can actually compute the intersections of an
isosurface at those 1.25 million points for a cosine function!

Your enthusiasm is certainly commendable, but unfortunately also
misdirected. What you suggest simply won't be faster for anything but very,
very, very special case in very, very, very limited circumstance.

	Thorsten

Post a reply to this message

Thorsten Froehlich <tho### [at] trfde> wrote:

> You do realize that such a grid would also require eight gigabyte (!!!) of
> memory? Using the fastest personal computer memory available, using a 128
> bit wide bus, you still need more than a second to even fill that grid with
> any single identical value.
Hello,
I'm aware of that since my "workstation" is only a P3 550MH/256M :). This
parser is useful also for smaller grid and "big" isosurfaces fcts: I'm
usually using grids of 120x120x120 and some isosurfaces can take from 1 to
60 minutes to compute. Isosurfaces aren't usually as simples as Schwartz
example and can take hours to compute even for smaller grids.
Perhaps I didn't  explain the kind of simplification we can do :
An isosurface formulas with handreds of "cos(x)", "cos(y)" and "cos(z)"
terms inside it will also require as much cosinus fct calculations as
Schwartz!
A normal parser can make (very difficult to generalise but possible) some
simplifications like : cos(x) + cos(x)*y*log(abs(cos(x)) +1)) =
cos(x)*(1+y*log(abs(cos(x)) +1)). However, it still need to compute two
cosinus for every point of the grid (for a grid of 100x100x100, the number
of cosinus is 2x10^6).
Now for the geometrical parser:
1) we dont have to apply any factorisation.
2) The number of required cosinus calculation is 100 for the grid
100x100x100.
3) you can add as much as you want termes of "cos(x)" inside the isosurface
formula, the required number of cosinus calculations is always 100...
In fact, "cos(x)" termes are treated like constantes and this is possible
only because we use a geometrical parser (the parser is able to retreave
the value of "cos(x)" as easy as retreaving the value of the terme "x"
itself).
This example : cos(x)*(1+y*log(abs(cos(x)) +1))+ log(x*cos(x))+
pow(tan(abs(cos(x)+1))) + log(y*cos(x)*cos(x)) require also 100 cosinus
calculations without any hard work of factorisation.
3) This works for all unary fcts defined inside the parser (log, exp, ...)

I'm working on expanding this to binary defined fcts like "+" and "*" : this
is a little bit triky but there is a way to make termes like "cos(x*y)"
requiring as little as 10^4 cosinus calculations for the grid 100x100x100,
and this for any number of termes "cos(x*y)" inside the Isosurface formula.
Also, this way of factorisation has nothing to do with the inclusion of
mutiple instructions inside the evaluator fct of the parser, it's just
another way of sampling the calculations by using the cube symetrys.

> Your enthusiasm is certainly commendable, but unfortunately also
> misdirected.
We will try to find more convincing results :)
Cheers,
Taha

Post a reply to this message

It's unclear to me what is it that you are doing with the isosurfaces.
Are you rendering 3D images of them, or are you doing something completely
different?

  If you are rendering an image of the isosurface, what is it that you
need a grid of values for? Are you tesselating the isosurface (using
something like the marching cubes algorithm)?

  POV-Ray doesn't tesselate isosurfaces. It traces them on the fly.
A random ray would basically never hit exactly the node of a regular
3D grid of values. The only way to take advantage of such a grid would
be to interpolate from nearby values (which is, in fact, what tesselation
does basically; the "interpolation" is just done using triangles).

-- 
                                                          - Warp

Post a reply to this message

Hi Warp,
In K3DSurf, I'm using mainly three steps to draw Isosurfaces:
1) Grid initialisation:
Preparing the Grid: the grid is a 3D matrix that holds values obtained from
parsing the Isosurface math formula. This process can be simplified to that
code :
++++++++++++++++++++++++++++++++++++++++++++++++++++
void InitGrig(){
     FunctionParser Isosurface;
     double Var[3];
  double X[100], y[100], Z[100];
  // init the parser:
     Isosurface.Parse("cos(x)+cos(y)+cos(z)+log(abs(cos(x))+1)", "x,y,z");
     // Init the axes to a Cube [-1,1][-1,1][-1,1]
  for(int i=0; i<100; i++) X[i] = y[i] = Z[i] = -1+i*(2/100);
  // Compute the Grid:
     for(i=0; i<limitX; i++) {
         Var[0] = X[i];  // Axe x
      for(j=0; j<limitY; j++) {
         Var[2] = Y[j];  // Axe y
       for(k=0; k<limitZ; k++) {
         Var[3] = Z[k];  // Axe z
         Grid[i][j][k] = Isosurface.Eval(Var);
       }}}
}
+++++++++++++++++++++++++++++++++++++++++++++++++++++
That's it. This simple example show how to compute the grid for the
isosurface : F(x, y, z) = cos(x)+cos(y)+cos(z)+log(abs(cos(x))+1)
2) Generate Triangles from the grid:
Use the marching cube algorithm to generate triangles from the matrix Grid.
3) Draw the image in the window:
Use OpenGl to draw thoses triangles.

The steps 2 and 3 are now extremely fast in K3DSurf: This part of code is
able to treat
(as Thorsten said) millions of triangles in less than a second.
Infortunatly, the first step is very slow compared the the others.
PovRay is not using the marching cube algorithm but it's also using a
mathematical parser in the process of generating the 3D images: even if
PovRay holds the fastest implementation of the raytracing algo, it's going
to be very slow because of it's math parser and this is especially true
when computing big grids and complicated Isosurfaces formulas.
I know that claiming to make a new parser "technologie" is quite crazy but
I'm here to show you that it's really what it looks like. The only
requirement to understand easily what I'm doing is just to have some
understanding of the 3D drawing technics because this solution is inspired
mainly from some basic geometry thoughts.
Now I'll try to explain what's the relationship between the parser and
Isosurfaces:
The original parser you created was calculating a single point at a time.
This is abvious when you look at the code shown in the top and specialy the
line:
>> Grid[i][j][k] = Isosurface.Eval(Var); // Var holds p(x,y,z) of a single point
Imagine now that we want to define a new Eval() fct that can take as
parameter sets of N points defined like this: P(X[N], Y[N], Z[N]) where X, Y
and Z and vectors.
For N =64,the code of the newly defined fct NewEval() will looks like this:

The original code:
-------------------------------------------
    for(IP=0; IP<ByteCodeSize; ++IP){
        switch(ByteCode[IP])
        {
         case   cCos: Stack[SP] = cos(Stack[SP]);
                      break;
........

// Variables:
          default:
              Stack[++SP] = Vars[ByteCode[IP]-VarBegin];
        }
    }

Where Stack = new double[StackSize];
and   Vars  = new double[3];
---------------------------------------------

Changed to :
---------------------------------------------

    for(IP=0; IP<ByteCodeSize; ++IP)
    {
        switch(ByteCode[IP])
        {
         case   cCos:
NewStack = &StackArray[SP*64];
for(i=0; i<64; i++) NewStack[i] = cos(NewStack[i]);
break;
........

// Variables:
          default:
    ++SP;
    NewStack = &StackArray[SP*64];
     if((ByteCode[IP]-VarBegin) == 0) memcpy(NewStack, &Vars[0]  ,
64*sizeof(double));
else if((ByteCode[IP]-VarBegin) == 1) memcpy(NewStack, &Vars[64] ,
64*sizeof(double));
else if((ByteCode[IP]-VarBegin) == 2) memcpy(NewStack, &Vars[124],
64*sizeof(double));

Where StackArray = new double[64*StackSize];
and   Vars       = new double[3*64];
-------------------------------------------------
This code holds almost all the idea behind this new kind of parsers.
Now the best part is that I'm going to use my knowlege of the isosurfaces
geometry to make the newly created function NewEval() working for real:
Imgine we have an isosurface to evaluate for a grid of 100x100x100, your
original code was used like I showd in code on the top.
To use the second newly created fct, the code on the top should be changed
to something like that :

++++++++++++++++++++++++++++++++++++++++++++++++++++
void InitGrig(){
     FunctionParser Isosurface;
     double Var[3*64];
  double X[100], y[100], Z[100];
  // init the parser:
     Isosurface.Parse("cos(x)+cos(y)+cos(z)+log(abs(cos(x))+1)", "x,y,z");
     // Init the axes to a Cube [-1,1][-1,1][-1,1]
  for(int i=0; i<100; i++) X[i] = y[i] = Z[i] = -1+i*(2/100);
  // Compute the Grid:

for(StartI=0; StartI<100; StartI+=4){
    l=0;
    for(i=0; i<4; i++)
       for(j=0; j<4; j++)
          for(k=0; k<4; k++)
                Vars[l++]      = X[i+StartI];
for(StartJ=0; StartJ<100; StartJ+=4){
    l=0;
    for(i=0; i<4; i++)
       for(j=0; j<4; j++)
          for(k=0; k<4; k++)
                Vars[64+(l++)]  = Y[j+StartJ];
for(StartK=0; StartK<100; StartK+=4){
    l=0;
    for(i=0; i<4; i++)
       for(j=0; j<4; j++)
          for(k=0; k<4; k++)
                Vars[128+(l++)]  = Z[k+StartK];


   Grid[i][j][k] = Isosurface.NewEval(Var);
}}}
}
+++++++++++++++++++++++++++++++++++++++++++++++++++++

That's it that's all! This code can be resumed to this : I splited my
isosurfaces calculations to chunks of 64 operations and that amount of
operations is equivalent to a small grid of 4x4x4.
To calculate the hole grid of 100x100x100, I have to use NewEval() : n =
(100x100x100)/(4x4x4) times.
So, now we can use the parser to evaluate more than one operations in the
loop.
I didn't finish yet...I just discovered (and implemented) a very efficient
way of sampling the calculations by using cubic symetrys.
I'll talk about it later.
Cheers,
Taha

Post a reply to this message

Well, that's the problem: You are tesselating the isosurface, POV-Ray
doesn't. POV-Ray wouldn't benefit from your optimizations because they
can be used only when tesselating the isosurface.

-- 
                                                          - Warp

Post a reply to this message

Warp <war### [at] tagpovrayorg> wrote:
> Well, that's the problem: You are tesselating the isosurface, POV-Ray
> doesn't. POV-Ray wouldn't benefit from your optimizations because they
> can be used only when tesselating the isosurface.
Hi,
I know that the 3 steps in generating images by PovRay doesn't have clear
frontiers like in thoses used in K3DSurf. However, it seems to me that
everything can be splitted that way: If we don't see it, that's because we
aren't used to it.
How about changing the raytracing technique to make it able to treat 64 rays
at the time? Is it really impossible to imagine that process ? I don't think
so...
I think that the real problem in using this technic for a raytracer will be
that the internal parser will "hardly" take a full advantage of that
implementation: Sampling the calculations by using cubic symetrys.
I'll try to explain what's this simplification and why a geometrical parser
is so important.
I want to make my code able to simplify the calculations by calculating only
the first "cos(x)" in the formulas. All other occurence of that terme will
be copied from the first result (in the example below, the parser will have
to calculate only the first "cos(x)"):
cos(x)+cos(y)+cos(z)+log(abs(cos(x))+1)

------------------------------------------------
int cosfctx = -1; // defined outside NewEval()
int usedvariable = -1;
double cosfctxT[depth];

  for(IP=0; IP<ByteCodeSize; ++IP)
    {
        switch(ByteCode[IP])
        {
         case   cCos:
NewStack = &StackArray[SP*64];
// This part is executed when the vriable in the top of the stack is "X"
if(variable == 0){
  if(cosfctx != 1){
//This code is executed only once in the process of generating the grid:
  for(i=0; i<100; i++) cosfctxT[i] = cos(X[i]);
  cosfctx = 1;
  }
    l=0;
    for(i=0; i<8; i++)
       for(j=0; j<8; j++)
          for(k=0; k<8; k++)
                NewStack[l++]  = cosfctxT[i+StartI];
  }

variable = -1;
break;
}

//This part is not always reached:
for(i=0; i<64; i++) NewStack[i] = cos(NewStack[i]);
break;

.............


// Variables:
          default:
    ++SP;
    NewStack = &StackArray[SP*64];
     if((ByteCode[IP]-VarBegin) == 0) {
memcpy(NewStack, &Vars[0]  ,64*sizeof(double));
variable = 0; // tell the parser which variable it's using
}
else if((ByteCode[IP]-VarBegin) == 1) memcpy(NewStack, &Vars[64] ,
64*sizeof(double));
else if((ByteCode[IP]-VarBegin) == 2) memcpy(NewStack, &Vars[124],
64*sizeof(double));

----------------------------------------------------------------
That code is now able to calculate only one occurence of "cos(x)" and use it
for all other termes in the formulas. More over, the parser is calculating
only 100 cosinus fct in the hole process. We can do the same for all unary
fct defined by the parser and for all defined variables. Inlike Thorsten, I
think that those cases are very commun and can do huge time saving.
Moreover, we can expand this technique to heigher level.
Also, this kind of simplification is possible only because the parser
"knows" the geometry it's using: it's impossible to do such thing in
general. Multiple instructions inside the loop can work but not this kind
of optimisation and that can be the case of PovRay but I still didn't
investigate all the possibilitys.
Hope it's clear to you.
Cheers,
Taha

Post a reply to this message

"virtualmeet" <tah### [at] yahoofr> wrote:
>     l=0;
>     for(i=0; i<8; i++)
>        for(j=0; j<8; j++)
>           for(k=0; k<8; k++)
>                 NewStack[l++]  = cosfctxT[i+StartI];

Correction :
-----------------------------------------------------
     l=0;
     for(i=0; i<4; i++)
        for(j=0; j<4; j++)
           for(k=0; k<4; k++)
                 NewStack[l++]  = cosfctxT[i+StartI];
------------------------------------------------------

That's because I'm using an extraordinary 512=8x8x8 version :)...
Taha

Post a reply to this message

virtualmeet wrote:
> Warp <war### [at] tagpovrayorg> wrote:
>> Well, that's the problem: You are tesselating the isosurface, POV-Ray
>> doesn't. POV-Ray wouldn't benefit from your optimizations because they
>> can be used only when tesselating the isosurface.
> Hi,
> I know that the 3 steps in generating images by PovRay doesn't have clear
> frontiers like in thoses used in K3DSurf. However, it seems to me that
> everything can be splitted that way: If we don't see it, that's because we
> aren't used to it.
> How about changing the raytracing technique to make it able to treat 64 rays
> at the time? Is it really impossible to imagine that process ? I don't think
> so...

Your code does *NOT* render an isosurface!!! It is a simple voxel algorithm,
something not even closely related to an isosurface defined by a function.
As i told you very early on, you can already do exactly that in POV-ray with
the density file/function feature. The fact remains it is a *different*
algorithm that is *not* related to isosurfaces. As I have also explained
multiple times now, it does *not* make sense to implement isosurfaces based
on a voxel array. If you trace 64 or one ray is completely irrelevant, as
you voxel array will not be faster to traverse if you use one ray or 64 at a
time.

Warp and I have now multiple times explained things, and you keep repeating
the same incorrect statements. What you say just has neither head nor tail
when applied to ray tracing. In fact the algorithm you describe isn't even
the best scan-line rendering algorithm for voxels.

Please, first understand how ray tracing works in general, what isosurfaces
actually are, how intersection with isosurfaces are computed, and after all
that you will realise your misconceptions. Please!

	Thorsten

Post a reply to this message

virtualmeet <tah### [at] yahoofr> wrote:
> How about changing the raytracing technique to make it able to treat 64 rays
> at the time?

  I don't think it's possible to raytrace an isosurface 64 rays at a time
and get a speed advantage.

  There may be a speed advantage if the isosurface function is evaluated
*only once* for each ray. However, it isn't. Searching for the intersection
between the ray and the isosurface is an iterative process: the function
is evaluated many times for one intersection calculation, and each
consecutive evaluation depends on the previous evaluations (and thus
cannot be made in parallel). It's basically a refining process where
values closer and closer to the actual isosurface are obtained from
the previous calculations.

  Calculating 64 rays at a time would probably not give speed benefits:
The parameters for each function evaluation for each ray for each substep
are different, so they cannot be all calculated using the same variables.
Moreover, the amount of substeps to be done may change from ray to ray
(if I have understood the isosurface raytracing algorithm correctly).

  This would also be feasible only for the first recursion level (ie. the
rays starting from the camera). Reflected/refracted rays could not possibly
use this kind of parallelism in any reasonable way (because it's impossible
to tell which group of rays actually hit the isosurface when they have been
bouncing around the scene).

  And of course there would be problems inherent with tracing 64 rays
at a time: If part of those rays actually don't hit the isosurface (eg.
because there's another object in front), the function evaluations done
for them would go to waste.

-- 
                                                          - Warp

Post a reply to this message