| 1 : |
jhr |
1115 |
// iso2d
|
| 2 : |
|
|
//
|
| 3 : |
|
|
// Demo of finding isocontours via Newton-Raphson method.
|
| 4 : |
|
|
// Initializes positions on a grid, and each update applies one
|
| 5 : |
|
|
// step of Newton-Raphson.
|
| 6 : |
|
|
//
|
| 7 : |
|
|
// Process output with:
|
| 8 : |
|
|
// unu jhisto -i mip.txt -b 512 512 -min 0 0 -max 1 1 | unu 2op neq - 0 | unu quantize -b 8 -o iso2d.png
|
| 9 : |
|
|
|
| 10 : |
|
|
int gridSize = 300;
|
| 11 : |
|
|
field#1(2)[] F = load("../data/ddro-80.nrrd") ⊛ ctmr;
|
| 12 : |
|
|
int stepsMax = 10;
|
| 13 : |
|
|
real epsilon = 0.0001;
|
| 14 : |
|
|
|
| 15 : |
|
|
strand sample (int ui, int vi) {
|
| 16 : |
|
|
output vec2 pos = [lerp(0.0, 1.0, -0.5, real(ui), real(gridSize)-0.5),
|
| 17 : |
|
|
lerp(0.0, 1.0, -0.5, real(vi), real(gridSize)-0.5)];
|
| 18 : |
|
|
// set the isvalue to 50, 30, or 10, depending on whichever we're closest to
|
| 19 : |
|
|
real isoval = 50.0 if F(pos) >= 40.0 else 30.0 if F(pos) >= 20.0 else 10.0;
|
| 20 : |
|
|
int steps = 0;
|
| 21 : |
|
|
update {
|
| 22 : |
|
|
// We bail if we're no longer inside or taken too many steps.
|
| 23 : |
|
|
// We tag this strand result as garbage by setting it to a
|
| 24 : |
|
|
// non-finite position; these will be filtered later.
|
| 25 : |
|
|
// This is really what "die" is for, though (see below).
|
| 26 : |
|
|
if (!inside(pos, F) || steps > stepsMax) {
|
| 27 : |
|
|
pos = [nan,nan];
|
| 28 : |
|
|
// HEY (BUG) we shouldn't even have to use non-finite values;
|
| 29 : |
|
|
// we should be able to just say "die;". Currently, using
|
| 30 : |
|
|
// "die" leads to the program never finishing...
|
| 31 : |
|
|
stabilize;
|
| 32 : |
|
|
}
|
| 33 : |
|
|
vec2 grad = ∇F(pos);
|
| 34 : |
|
|
real gmsq = grad • grad;
|
| 35 : |
|
|
// If |∇F|=0 we can't compute update; we have to bail
|
| 36 : |
|
|
if (gmsq == 0.0) {
|
| 37 : |
|
|
pos = [nan,nan];
|
| 38 : |
|
|
stabilize;
|
| 39 : |
|
|
}
|
| 40 : |
|
|
// delta = Newton-Raphson step
|
| 41 : |
|
|
vec2 delta = -(F(pos) - isoval)*grad/gmsq;
|
| 42 : |
|
|
// We've converged if the update is small enough
|
| 43 : |
|
|
if (|delta| < epsilon) {
|
| 44 : |
|
|
stabilize;
|
| 45 : |
|
|
}
|
| 46 : |
|
|
pos = pos + delta;
|
| 47 : |
|
|
steps = steps + 1;
|
| 48 : |
|
|
}
|
| 49 : |
|
|
}
|
| 50 : |
|
|
|
| 51 : |
|
|
initially { sample(ui, vi) | vi in 0..(gridSize-1), ui in 0..(gridSize-1) };
|
Parent Directory
| 
Revision Log
