POV-Ray : Newsgroups : povray.advanced-users : isosurface surface normals : Re: isosurface surface normals Server Time
28 Mar 2024 16:39:35 EDT (-0400)
  Re: isosurface surface normals  
From: Tor Olav Kristensen
Date: 25 Dec 2001 20:44:38
Message: <3C292BC5.5A6AB4A6@hotmail.com>
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

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