|
 |
Mike Horvath <mik### [at] gmail com> wrote:
> On 7/29/2018 3:57 PM, Tor Olav Kristensen wrote:
> > fn_Gradient() is not a function. It is a macro within math.inc. When you pass a
> > function to this macro it will create and return a new function that estimates
> > the magnitude of the gradient of the function that you passed to it.
> >
> > The function that you pass to this macro must take 3 arguments; usually x, y and
> > z (but their names doesn't really matter). And the function that it returns
> > takes 3 arguments.
> >
> > --
> > Tor Olav
> > http://subcube.com
> >
>
> I tried this code:
>
> #declare f_test = function(var1,var2,var3)
> {pow(var1,2)+pow(var2,2)+pow(var3,2)-1-0.5*f_noise3d(var1*3,var2*3,var3*3)}
> #declare f_normalized = function(var1,var2,var3,varA,varB,varC)
>
{f_test(var1,var2,var3)/sqrt(4*pow(var1,2)/pow(varA,4)+4*pow(var2,2)/pow(varB,4)+4*pow(var3,2)/pow(varC,4))}
>
> and this code:
>
> SetGradientAccuracy(0.0001)
> #declare f_test= function {x*x+y*y+z*z-1-0.5*f_noise3d(x*3,y*3,z*3)}
> #declare f_input = function {f_test(x,y,z)}
> #declare GradientFn = fn_Gradient(f_input)
> #declare f_normalized = function {f_input(x,y,z)/GradientFn(x,y,z)}
>
> And the latter was much slower. Is there a way to speed it up?
Those two code snippets does not calculate the same.
I do not know what is calculated in the denominator in
f_normalized() in the the first expression. (I even suspect that
it is not a useful calculation to do.)
But in the denominator in f_normalized() in the second expression,
the magnitude of the gradient for the f_test() function is
calculated - and this is a somewhat expensive calculation.
Note that there is no need to embed the f_test() function into
another function; f_input(), that just calls f_test(). The
fn_Gradient() macro can deal with f_test() directly, like this:
#include "math.inc"
#declare f_test =
function { x*x+y*y+z*z-1-0.5*f_noise3d(x*3,y*3,z*3) }
#declare f_gradient = fn_Gradient(f_test)
#declare f_normalized =
function { f_test(x,y,z)/f_gradient(x,y,z) }
/*
Note that pow(x, 2) may be faster to evaluate than x*x.
Also note that if there is no argument list after function,
then function(x, y, z) { } is assumed.
Your code can be made easier for others to read if you
present it with some consistent indentation and extra white
space. (It will then also be easier for others to help you.)
E.g. like this:
*/
#declare D = 3;
#declare E = 3;
#declare F = 3;
#undef f_test // To allow for a new version
#declare f_test =
function {
pow(x, 2) + pow(y, 2) + pow(z, 2) - 1
-0.5*f_noise3d(D*x, E*y, F*z)
}
#undef f_normalized // To allow for a new version
#declare f_normalized =
function(x, y, z, a_, b_, c_) {
f_test(x, y, z)
/
sqrt(
4*pow(x, 2)/pow(a_, 4)
+4*pow(y, 2)/pow(b_, 4)
+4*pow(z, 2)/pow(c_, 4)
)
}
/*
If your varA, varB and varC values are constants at render-time,
then they do not need to be in the argument list of
f_normalized(). This function can then be written like this:
*/
// You'll have to fill in the numbers you want here:
#declare A = 1.0;
#declare B = 0.7;
#declare C = 1.0;
#declare AA = pow(A, 2);
#declare BB = pow(B, 2);
#declare CC = pow(C, 2);
#undef f_normalized
#declare f_normalized =
function {
f_test(x, y, z)
/
sqrt(
4*pow(x, 2)/pow(AA, 2)
+4*pow(y, 2)/pow(BB, 2)
+4*pow(z, 2)/pow(CC, 2)
)
}
/*
- or like this:
*/
#undef f_normalized
#declare f_normalized =
function {
f_test(x, y, z)
/
(
2*
sqrt(
pow(x/AA, 2)
+pow(y/BB, 2)
+pow(z/CC, 2)
)
)
}
/*
The f_r() function can be found in functions.inc.
It can be used instead of the sqrt() and pow() expressions above.
Documentation for f_r():
http://www.povray.org/documentation/view/3.7.0/448/
*/
#undef f_normalized
#declare f_normalized =
function {
f_test(x, y, z)/(2*f_r(x/AA, y/BB, z/CC))
}
/*
Now onto the macros for estimating the magnitude of the gradient
for 3D-functions.
Here's some reading that is relevant:
Numerical differentiation of functions:
https://en.wikipedia.org/wiki/Numerical_differentiation
The gradient of functions:
https://math.oregonstate.edu/home/programs/undergrad/CalculusQuestStudyGuides/vcalc/grad/grad.html
As you can read in the Wikipedia article, there are several ways
to estimate the derivative of a function.
Below are some example macros for doing this with a passed
function.
The 3b macro below, which is very simular to the fn_Gradient()
macro in math.inc, is the one to use for most accurate results.
But it does return a function that takes longer time to
evaluate than the functions returned from the 1b or 1c macros.
*/
#macro MakeGradientMagnitudeFunction_1a(Fn, h)
function {
sqrt(
pow((Fn(x + h, y, z) - Fn(x, y, z))/h, 2)
+pow((Fn(x, y + h, z) - Fn(x, y, z))/h, 2)
+pow((Fn(x, y, z + h) - Fn(x, y, z))/h, 2)
)
}
#end // MakeGradientMagnitudeFunction_1a
#macro MakeGradientMagnitudeFunction_2a(Fn, h)
function {
sqrt(
pow((Fn(x, y, z) - Fn(x - h, y, z))/h, 2)
+pow((Fn(x, y, z) - Fn(x, y - h, z))/h, 2)
+pow((Fn(x, y, z) - Fn(x, y, z - h))/h, 2)
)
}
#end // MakeGradientMagnitudeFunction_2a
#macro MakeGradientMagnitudeFunction_3a(Fn, h)
function {
sqrt(
pow((Fn(x + h, y, z) - Fn(x - h, y, z))/(2*h), 2)
+pow((Fn(x, y + h, z) - Fn(x, y - h, z))/(2*h), 2)
+pow((Fn(x, y, z + h) - Fn(x, y, z - h))/(2*h), 2)
)
}
#end // MakeGradientMagnitudeFunction_3a
/*
Now replace sqrt() and pow() with f_r().
As long as h (or h2) is positive we can divide by h (or h2)
outside the f_r() expression.
*/
#macro MakeGradientMagnitudeFunction_1b(Fn, h)
function {
f_r(
Fn(x + h, y, z) - Fn(x, y, z),
Fn(x, y + h, z) - Fn(x, y, z),
Fn(x, y, z + h) - Fn(x, y, z)
)/h
}
#end // MakeGradientMagnitudeFunction_1b
#macro MakeGradientMagnitudeFunction_2b(Fn, h)
function {
f_r(
Fn(x, y, z) - Fn(x - h, y, z),
Fn(x, y, z) - Fn(x, y - h, z),
Fn(x, y, z) - Fn(x, y, z - h)
)/h
}
#end // MakeGradientMagnitudeFunction_2b
#macro MakeGradientMagnitudeFunction_3b(Fn, h)
#local h2 = h*2;
function {
f_r(
Fn(x + h, y, z) - Fn(x - h, y, z),
Fn(x, y + h, z) - Fn(x, y - h, z),
Fn(x, y, z + h) - Fn(x, y, z - h)
)/h2
}
#end // MakeGradientMagnitudeFunction_3b
/*
In the 1b macro there are 6 calls to Fn(), where 3 of these
calls return the same value. If we rewrite this macro we can
reduce this to 4 calls to Fn().
But then we must add a call to a helper function, so we may
not gain much...
*/
#macro MakeGradientMagnitudeFunction_1c(Fn, h)
#local G_Fn_ =
function(f_, fx_, fy_, fz_) {
f_r(fx_ - f_, fy_ - f_, fz_ - f_)/h
}
function {
G_Fn_(
Fn(x , y , z ),
Fn(x + h, y , z ),
Fn(x , y + h, z ),
Fn(x , y , z + h)
)
}
#end // MakeGradientMagnitudeFunction_1c
// Here's some tests of the macros above:
#include "functions.inc"
#declare A = 3;
#declare B = 3;
#declare C = 3;
#undef f_test
#declare f_test =
function {
pow(x, 2) + pow(y, 2) + pow(z, 2) - 1
- 0.5*f_noise3d(A*x, B*y, C*z)
}
#declare H = 1e-5;
#declare Fn1a = MakeGradientMagnitudeFunction_1a(f_test, H)
#declare Fn1b = MakeGradientMagnitudeFunction_1b(f_test, H)
#declare Fn1c = MakeGradientMagnitudeFunction_1c(f_test, H)
#declare Fn2a = MakeGradientMagnitudeFunction_2a(f_test, H)
#declare Fn2b = MakeGradientMagnitudeFunction_2b(f_test, H)
#declare Fn3a = MakeGradientMagnitudeFunction_3a(f_test, H)
#declare Fn3b = MakeGradientMagnitudeFunction_3b(f_test, H)
#declare P = <12, -4, 3>;
#declare G1a = Fn1a(P.x, P.y, P.z);
#declare G1b = Fn1b(P.x, P.y, P.z);
#declare G1c = Fn1c(P.x, P.y, P.z);
#declare G2a = Fn2a(P.x, P.y, P.z);
#declare G2b = Fn2b(P.x, P.y, P.z);
#declare G3a = Fn3a(P.x, P.y, P.z);
#declare G3b = Fn3b(P.x, P.y, P.z);
#debug "\n\n"
#debug concat("G1a: ", str(G1a, 0, -1), "\n")
#debug concat("G1b: ", str(G1b, 0, -1), "\n")
#debug concat("G1c: ", str(G1c, 0, -1), "\n")
#debug "\n"
#debug concat("G2a: ", str(G2a, 0, -1), "\n")
#debug concat("G2b: ", str(G2b, 0, -1), "\n")
#debug "\n"
#debug concat("G3a: ", str(G3a, 0, -1), "\n")
#debug concat("G3b: ", str(G3b, 0, -1), "\n")
#debug concat("(G1b+G2b)/2: ", str((G1b + G2b)/2, 0, -1), "\n")
#debug "\n\n"
/*
And finally...
If you really need a, b and c to be in the argument list for
the passed function, then the macro can be written like this:
*/
#macro MakeGradientMagnitudeFunction(Fn, h)
function(x, y, z, a_, b_, c_) {
f_r(
Fn(x + h, y, z, a_, b_, c_) - Fn(x, y, z, a_, b_, c_),
Fn(x, y + h, z, a_, b_, c_) - Fn(x, y, z, a_, b_, c_),
Fn(x, y, z + h, a_, b_, c_) - Fn(x, y, z, a_, b_, c_)
)/h
}
#end // MakeGradientMagnitudeFunction
/*
- But then it does not make sense, and it will not work, if the
f_test() function does not also have them in its argument list...
#undef f_test
#declare f_test =
function(x, y, z, a_, b_, c_) { ... }
#declare f_gradient_magnitude =
MakeGradientMagnitudeFunction(f_test, H)
#undef f_normalized
#declare f_normalized =
function(x, y, z, a_, b_, c_) {
f_test(x, y, z, a_, b_, c_)
/
f_gradient_magnitude(x, y, z, a_, b_, c_)
}
*/
--
Tor Olav
http://subcube.com
Post a reply to this message
|
|
