POV-Ray : Newsgroups : povray.newusers : isosurface problem Server Time
5 Sep 2024 04:18:51 EDT (-0400)
  isosurface problem (Message 1 to 5 of 5)  
From: pete
Subject: isosurface problem
Date: 3 Feb 2002 18:51:48
Message: <3c5dcd14@news.povray.org>
Hi,

I'm trying to use the isosurface functionality of POV-Ray to
depict (in a loose, not physically accurate, way ;) an isosurface
of magnetic field strength around a moving point charge.

I'm using the "Biot-Savart law" where, if the particle has
velocity v and you put the particle at the origin, the field
strength at some point p is proportional to | ((v x p) / |p|^3) |
(where 'x' is the vector cross-product, and '|z|' means "vector
length of z"). 

I've worked around the difficulty of not being able to access
vector functions inside POV-Ray function definitions by defining
vector length and "cross product magnitude" functions, but I still
can't get it to work. Instead of a nice donut facing in the
direction of v, I just get whatever bounding box I choose for
the isosurface, completely filled up.

Here's my povray source

8<---C-U-T--H-E-R-E----
#version 3.5;
#include "colors.inc"
#include "functions.inc"
#declare CamLoc = < -28000000, 0,0>;
camera { location CamLoc look_at <0,0,0> angle 45 }
light_source {CamLoc color White*0.35}
light_source {<-50, 150,-75> color White}
#declare IsoFinish = finish { ambient 1 diffuse 1 specular 1
                              roughness 0.02 brilliance 2 }
#declare vel = <35,25,45>;
#declare velx = vel.x;
#declare vely = vel.y;
#declare velz = vel.z;
#declare f_vlen = function(x,y,z) { sqrt(x^2 + y^2 + z^2) }
#declare f_vcrossmag = function(v1x,v1y,v1z,v2x,v2y,v2z) {
        sqrt( (v1y*v2z - v1z*v2y)^2 + (v1z*v2x - v1x*v2z)^2
              + (v1x*v2y - v1y*v2x)^2 )
}
#declare f_bfield = function(x,y,z,vx,vy,vz) {
        f_vcrossmag(vx,vy,vz,x,y,z) / pow(f_vlen(x,y,z),3)
}
isosurface { function {f_bfield(x,y,z,velx,vely,velz) - 6}
        max_gradient 35000 threshold 0 contained_by {box{-8000000, 8000000}}
        texture { pigment {color Red} finish {IsoFinish} }
        scale 1/vlength(1)
}
8<---C-U-T--H-E-R-E----

Being convinced that the equation should yield more sensible
results, I wrote a perl script that positions spheres
along vectors chosen at random, attempting to find the point
along the vector for which the equation gives a magnitude
of 6 for the field, so it does a very rough approximation of
what I expected the POV-Ray isosurface to produce. In order
to make this a good debugging test, I implemented the functions
in perl in exactly the same way I've done them in povray.

This seems to work as expected, giving the plump donut shape,
lined up with the velocity vector.

As I'm by no means fully conversant with either the maths or the
POV-Ray isosurfaces and functions, it's quite possible I've
made a mistake somewhere - perhaps in my basic understanding of
what isosurfaces do and their limits, which is why I'm posting
this in the newusers group - but I *have* checked over the
equivalence between my perl and POV-Ray definitions of the
functions quite carefully, so I'd be interested if anyone can
tell me where my POV-Ray script is going wrong!

Just for reference, here's the perl script I used - if you want
to run it, you'll need to send the output to a file, like
this:
   ./bf.pl > bfield.pov

8<----C-U-T--H-E-R-E----
#! /usr/bin/perl
$velx = 35; $vely = 25; $velz = 45; $cmag = 6; $points = 10000;
print STDERR "velocity <$velx,$vely,$velz>, magnitude $cmag.\n";
sub vlen($$$) {
        my($x,$y,$z) = @_;
        sqrt($x*$x + $y*$y + $z*$z);
}
sub vcrossmag($$$$$$) {
        my ($ux,$uy,$uz,$vx,$vy,$vz) = @_;
        sqrt(($uy*$vz - $uz*$vy)**2 + ($uz*$vx - $ux*$vz)**2
             + ($ux*$vy - $uy*$vx)**2);
}
sub f_bfield($$$$$$) {
        my ($x,$y,$z,$vx,$vy,$vz) = @_; 
        vcrossmag($vx,$vy,$vz,$x,$y,$z) / ((vlen($x,$y,$z))**3);
}
sub findmag($$$$) {
        my ($x,$y,$z,$mag) = @_;
        my $scale = 1; my $tried = 0; my $bestdiff = $mag;
        my $bestscale = $scale;
        while($tried++ <= 1000) {
                $m = f_bfield($x*$scale,$y*$scale,$z*$scale,$velx,$vely,$velz);
                if (abs($m - $mag) < $bestdiff) {
                        $bestdiff = abs($m - $mag);
                        $bestscale = $scale;
                }
                if ($m == 0) {return 0;}
                if (abs($m-$mag) <= $mag/100) {return $scale;}
                $scale *= sqrt($m/$mag);
        }
        print STDERR "kludge not working, returning $bestscale with bestdiff
$bestdiff\n";
        return $bestscale;
}
$vvscale = vlen($velx,$vely,$velz);
$vvx = $velx/$vvscale;
$vvy = $vely/$vvscale;
$vvz = $velz/$vvscale;
print "#declare velvect = <$vvx,$vvy,$vvz>;\n";
print "#declare velscale = $vvscale;\n";
print "#declare bobbins = union {\n";
for (1..$points) {
        $x = rand(0.2) - 0.1;
        $y = rand(0.2) - 0.1;
        $z = rand(0.2) - 0.1;
        $scale = findmag($x,$y,$z,$cmag);
        ($x,$y,$z) = (($x*$scale), ($y*$scale), ($z*$scale));
        print "\tsphere { <$x,$y,$z>, 0.05 }\n";
}
print "}\n";
print <<END_POV;
#include "colors.inc"
global_settings {
        assumed_gamma 2.2
}
cylinder { 0, velvect*velscale 0.1 pigment{rgb<1,1,1>}}
#declare CamLoc = vnormalize(vcross(y,vcross(x,velvect))) * velscale/5;
camera {
        location CamLoc
        look_at < 0, 0, 0>
//        direction <0,0,1.3333>
        angle 45
}
light_source {CamLoc color White*0.35}
light_source {<-50, 150,-75> color White}
#declare IsoFinish =
finish {
        ambient 0.31 diffuse 1
        specular 1 roughness 0.02
        brilliance 2
}
object {
        bobbins
        rotate y*clock*180
        pigment { rgb<0.9,0.1,0.1>}
        finish {IsoFinish}
}
END_POV
8<----C-U-T--H-E-R-E----

Well, I hope that's not been *too* long a post, and thanks for having
a look!

- Pete.


Post a reply to this message

From: R  Suzuki
Subject: Re: isosurface problem
Date: 3 Feb 2002 22:21:24
Message: <3c5dfe34@news.povray.org>
<pet### [at] somewhereinuk> wrote 
> Here's my povray source

Try the below code.  
Is this what you want?

R. Suzuki
-----------------------------------------------------
#version 3.5;
#include "colors.inc"
#include "functions.inc"
#declare CamLoc = <-8, 0,0>;
camera { location CamLoc look_at <0,0,0> angle 45 }
light_source {CamLoc color White*0.35}
light_source {<-50, 150,-75> color White}
#declare IsoFinish = finish { ambient 0.35 diffuse 1 specular 1
                              roughness 0.02 brilliance 2 }
#declare vel = <35,25,45>;
#declare velx = vel.x;
#declare vely = vel.y;
#declare velz = vel.z;
#declare f_vcrossmag = function(v1x,v1y,v1z,v2x,v2y,v2z) {
        sqrt( (v1y*v2z - v1z*v2y)^2 + (v1z*v2x - v1x*v2z)^2
              + (v1x*v2y - v1y*v2x)^2 )
}
#declare f_bfield = function(x,y,z,vx,vy,vz) {
        f_vcrossmag(vx,vy,vz,x,y,z) / pow(f_r(x,y,z),3)
}
isosurface { function {
               max(-1, 
                   -(f_bfield(x,y,z,velx,vely,velz) - 6)
               )
                  }
        max_gradient 20 threshold 0 contained_by {box{-3.2, 3.2}}  
        accuracy 0.01      
        evaluate 10,1.1,0.99
        texture { pigment {color Red} finish {IsoFinish} }
        scale 1/vlength(1)
}


Post a reply to this message

From: pete
Subject: Re: isosurface problem
Date: 4 Feb 2002 00:25:24
Message: <3c5e1b44$1@news.povray.org>
>   From: "R. Suzuki" <r-s### [at] aistgojp>

> Try the below code.
> Is this what you want?

Absolutely, it does just what I was trying to do! Thanks!

Still, I would like to understand why your version
works and mine doesn't.

[...]
>        f_vcrossmag(vx,vy,vz,x,y,z) / pow(f_r(x,y,z),3)

Aha, so there is a provided vlength function, f_r(). Good!

[...]
>isosurface { function {
>               max(-1,
>                   -(f_bfield(x,y,z,velx,vely,velz) - 6)
>               )
>                  }

I admit I don't get this at all.

I notice I can take the max() part out (is it some kind of
error trap?), leaving the function definition as just

       -(f_bfield(x,y,z,velx,vely,velz) - 6)

and I can change this to

        6 - f_bfield(x,y,z,velx,vely,velz)

and it still works. But

        f_bfield(x,y,z,velx,vely,velz) - 6

(as in my original) doesn't work.

I thought that the surface would be drawn where the function becomes
zero, but then if 6 - f_bfield() is zero then f_bfield() - 6 should be
zero too! So I don't understand why this change makes such a
huge difference.
 
[...]
>        evaluate 10,1.1,0.99

I'll have to look this one up in the docs, too, since it sounds
fundamental, and I've no idea what it's doing!

Many thanks for the speedy and very helpful reply.

- Pete.


Post a reply to this message

From: R  Suzuki
Subject: Re: isosurface problem
Date: 4 Feb 2002 05:34:57
Message: <3c5e63d1@news.povray.org>
<pet### [at] somewhereinuk> wrote
>>isosurface { function {
>>               max(-1,
>>                   -(f_bfield(x,y,z,velx,vely,velz) - 6)
>>               )
>>                  }
>
>I admit I don't get this at all.
>
>I notice I can take the max() part out (is it some kind of
>error trap?), leaving the function definition as just
>
>       -(f_bfield(x,y,z,velx,vely,velz) - 6)
>
>and I can change this to
>
>        6 - f_bfield(x,y,z,velx,vely,velz)
>
>and it still works.

The max() and evaluate are tips for faster rendering.
You don't need them if the rendering speed is acceptable for you.

Without the max(), the maximum gradient of the function 
is extremely high.  It slows down rendering.


> But
>
>        f_bfield(x,y,z,velx,vely,velz) - 6
>
>(as in my original) doesn't work.
>
>I thought that the surface would be drawn where the function becomes
>zero, but then if 6 - f_bfield() is zero then f_bfield() - 6 should be
>zero too! So I don't understand why this change makes such a
>huge difference.

Because the isosurface is not a surface object---it is a solid object
by default.  If you want to see only the surface of your original 
function, you should add 'open' keyword. (see POV-Help 6.5.4.1)

>>        evaluate 10,1.1,0.99
>
>I'll have to look this one up in the docs, too, since it sounds
>fundamental, and I've no idea what it's doing!

Well, this is the dynamic max_gradient estimation technique and
it has not been documented yet.  See my messages
http://news.povray.org/3bd7f6e8$1@news.povray.org
http://news.povray.org/3be10cf9@news.povray.org

This technique is effective in this case because large gradient 
values are located in a small region and the gradient changes smoothly.

>Many thanks for the speedy and very helpful reply.

You are welcome.

R. Suzuki


Post a reply to this message

From: pete
Subject: Re: isosurface problem
Date: 4 Feb 2002 08:51:37
Message: <3c5e91e9@news.povray.org>
>   From: "R. Suzuki"

>The max() and evaluate are tips for faster rendering.
>You don't need them if the rendering speed is acceptable for you.

>Without the max(), the maximum gradient of the function
>is extremely high.  It slows down rendering.

Ok.

>>I thought that the surface would be drawn where the function becomes
>>zero, but then if 6 - f_bfield() is zero then f_bfield() - 6 should be
>>zero too! So I don't understand why this change makes such a
>>huge difference.
>
>Because the isosurface is not a surface object---it is a solid object
>by default.  If you want to see only the surface of your original
>function, you should add 'open' keyword. (see POV-Help 6.5.4.1)

Ah, yes! my problem was that the bfield function gets *less* as you
go away from the origin - the opposite of "simpler" functions
what I'm used to, where  the function gets bigger as x, y and z
increase.

So my function was, essentially, inside-out!

POV-Ray is counting all the values which are less than the threshold
as "inside", and those greater as "outside" this solid object.

With "6 - f_bfield()" the negative bits (ie less than threshold, with
threshold == 0) are all where they should be (between the origin and
the surface) and with "f_bfield - 6" they are all where they
*shouldn't* be - everywhere "outside" the surface, using "outside"
to mean the side of the surface that the camera is on, rather than
what POV-Ray means by it. So the change serves to align POV-Ray's
notion of "inside" and "outside" (in this case) with my simple
intuitive one.

I guess when I was rendering with "f_bfield() - 6" there was a
donut-shaped hole right in the middle of my huge red cube, which I
would have seen if I'd chopped it open ;-)  

This makes sense. As someone once said: "Now I can go on!"

>>>        evaluate 10,1.1,0.99

>Well, this is the dynamic max_gradient estimation technique and
>it has not been documented yet.

Right. Having re-read the isosurface part of the beta documentation,
it seems to be saying "don't worry about it, just put these magic
numbers in your script" ;-)

> See my messages
>http://news.povray.org/3bd7f6e8$1@news.povray.org
>http://news.povray.org/3be10cf9@news.povray.org

I'll have a look. Thanks again.

- Pete.


Post a reply to this message

Copyright 2003-2023 Persistence of Vision Raytracer Pty. Ltd.