Barron Gillon wrote:
> 
> Does anyone know (preferably have a macro they could send me) how to find
> the surface normal of an isosurface at an arbitrary <x,y,z>?  If I remember
> correctly, a macro that returns <g'(x),h'(y),i'(z)> where f(x,y,z) =
> <g(x),h(y),i(z)> and x, y, and z are specified will suffice.  The catch is
> of course that f() would be an arbitrary function, and I don't know how to
> find numeric derivatives.  Did I miss something in the standard includes
> that would do this?  Has anyone else done this?  Rune?  Thanks
>...

Warp has shown in another post how to find the gradient
of a function f(x, y, z) for any point in 3D space by
analytical calculation of the 3 partial derivatives of
the function.

But if you want to have POV to estimate this gradient
by applying some numerical methods, (in order to find
the partial derivatives), then you may have a look at
the code I provided in my answer to Jan Walzer's thread
12. Nov. 2001:

"Help on isosurface - distances - other stuff..."

news://news.povray.org/3BF0541A.62BC8C1F%40hotmail.com
http://news.povray.org/povray.advanced-users/19992/


Tor Olav

Post a reply to this message

Barron Gillon wrote:
> 
> Does anyone know (preferably have a macro they could send me) how to find
> the surface normal of an isosurface at an arbitrary <x,y,z>?  If I remember
> correctly, a macro that returns <g'(x),h'(y),i'(z)> where f(x,y,z) =
> <g(x),h(y),i(z)> and x, y, and z are specified will suffice.  The catch is
> of course that f() would be an arbitrary function, and I don't know how to
> find numeric derivatives.
>...

At the end of my code below you'll find such a macro.

There's also an analytical solution provided, so that
you can compare the results from both the methods.
(I have chosen a torus function for this example.)

POV-Ray v3.5 are needed for the code below.
(Because it allows passing of functions as parameters
to macros.)

Note that there are other and more accurate (and
slower) numerical methods to find derivatives for
such functions, but the one used below will suffice
for many applications.


Tor Olav


// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7
// Copyright 2001 by Tor Olav Kristensen
// Email: tor### [at] hotmailcom
// http://www.crosswinds.net/~tok
// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7

#version 3.5;

// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7

#macro PrintVector(v0)

 #debug concat("<", vstr(3, v0, ", ", 0, -1), ">")

#end // macro PrintVector

// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7
// Define a torus function

#declare TorusFn =
function(x, y, z, Rmaj, Rmin) {
  sqrt((sqrt(x^2 + z^2) - Rmaj)^2 + y^2) - Rmin
}

// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7
// Define constants for the torus function that we are investigating.
// Also define a point in space that we are going to find the gradient
// vector at.

#declare Rmajor = 10;
#declare Rminor = 2;
#declare pA = <-3, 2, -5>*3;

// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7
// Find the analytical solution
// Note that the partial derivatives are not dependant of R minor.)

#declare DerxTorusFn =
function(x, y, z, Rmaj) {
  x/sqrt((sqrt(x^2 + z^2) - Rmaj)^2 + y^2)
  *(sqrt(x^2 + z^2) - Rmaj)/sqrt(x^2 + z^2)
}

#declare DeryTorusFn =
function(x, y, z, Rmaj) {
  y/sqrt((sqrt(x^2 + z^2) - Rmaj)^2 + y^2)
}

#declare DerzTorusFn =
function(x, y, z, Rmaj) {
  z/sqrt((sqrt(x^2 + z^2) - Rmaj)^2 + y^2)
  *(sqrt(x^2 + z^2) - Rmaj)/sqrt(x^2 + z^2)
}

#declare vGradAnalytical =
<
  DerxTorusFn(pA.x, pA.y, pA.z, Rmajor),
  DeryTorusFn(pA.x, pA.y, pA.z, Rmajor),
  DerzTorusFn(pA.x, pA.y, pA.z, Rmajor)
>;

#debug "\n"
PrintVector(vGradAnalytical)
#debug "\n"

// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7
// Find a numerical solution

#macro vGradient(Fn, p0)
  
  #local H = 0.00001;

  <Fn(p0.x + H, p0.y, p0.z) - Fn(p0.x - H, p0.y, p0.z),
   Fn(p0.x, p0.y + H, p0.z) - Fn(p0.x, p0.y - H, p0.z),
   Fn(p0.x, p0.y, p0.z + H) - Fn(p0.x, p0.y, p0.z - H)>/2/H;

#end // macro vGradient

#declare vGradNumerical =
vGradient(function { TorusFn(x, y, z, Rmajor, Rminor) }, pA);

#debug "\n"
PrintVector(vGradNumerical)
#debug "\n\n"

// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7

Post a reply to this message

Wasn't it Barron Gillon who wrote:
>not trace!  I need to find the normal at an arbitrary point that
>_may_or_may_not be on the surface.  Trace finds an intersection with a
>surface and then returns the normal at the intersection, but I need the
>normal from an arbitrary point given an arbitrary function.  Trace is good,
>but not quite what I'm looking for.  Thanks anyway.

Here's an example


global_settings {assumed_gamma 1.0}

camera { location  <-1, 1, -4.5> look_at <0, 0, 0>}

sky_sphere { pigment { rgb 1} }

light_source {<-100,200,-500> colour rgb 1}

// ----------------------------------------

// An arbitrary function
#declare F=function { sin(1.5*x)*y }

// An isosurface using this function
isosurface {function {F(x,y,z)}
        threshold (F(0.2,-0.2,0))
        max_gradient 5
        contained_by{box{-2,2}} open
        pigment {rgb 0.9}
        finish {phong 0.5 phong_size 10}
}

// Scan for "normals"
#declare Norm=<0,0,0>;
#declare X=-2;
#while (X<2)
 #declare Y=-2;
 #while (Y<2)
  #declare P=<X,Y,-0>;//Current point of interest
  // Declare a temporary isosurface 
  #declare Temp = 
    isosurface {function {F(x,y,z)}
      threshold F(X,Y,0)
      max_gradient 5000
      contained_by{box{-2.1,2.1}} open
      }
  // Trace it
  #declare Q=trace(Temp, P-y*0.1, y, Norm);
  // Just in case P-y*0.1 happens to be in the surface
  #if (vlength(Norm)=0)
    #declare Q=trace(Temp, P-x*0.1, x, Norm);
  #end
  #if (vlength(Norm)=0)
    #declare Q=trace(Temp, P-z*0.1, z, Norm);
  #end

  // A cylinder pointing in the direction of the "normal"        
  cylinder {P,P+Norm*0.2,0.02 pigment {rgb x} no_shadow}
 #declare Y=Y+0.2;
 #end
#declare X=X+0.2;
#end

Post a reply to this message

Tor Olav Kristensen wrote:
>...
> 
> #macro vGradient(Fn, p0)
> 
>   #local H = 0.00001;
> 
>   <Fn(p0.x + H, p0.y, p0.z) - Fn(p0.x - H, p0.y, p0.z),
>    Fn(p0.x, p0.y + H, p0.z) - Fn(p0.x, p0.y - H, p0.z),
>    Fn(p0.x, p0.y, p0.z + H) - Fn(p0.x, p0.y, p0.z - H)>/2/H;
> 
> #end // macro vGradient
>...

I just noticed that this macro had a little error.
Here's a better version of it:


#macro vGradient(Fn, p0)

  #local H = 0.00001;

  (<Fn(p0.x + H, p0.y, p0.z) - Fn(p0.x - H, p0.y, p0.z),
    Fn(p0.x, p0.y + H, p0.z) - Fn(p0.x, p0.y - H, p0.z),
    Fn(p0.x, p0.y, p0.z + H) - Fn(p0.x, p0.y, p0.z - H)>/2/H)

#end // macro vGradient


Tor Olav

Post a reply to this message