Home My Page Projects Code Snippets Project Openings diderot
Summary Activity Tracker Tasks SCM

SCM Repository

[diderot] View of /branches/pure-cfg/test/iso2d.diderot
ViewVC logotype

View of /branches/pure-cfg/test/iso2d.diderot

Parent Directory Parent Directory | Revision Log Revision Log


Revision 957 - (download) (annotate)
Sat Apr 23 12:51:20 2011 UTC (8 years, 4 months ago) by glk
File size: 2337 byte(s)
better documenting of nan bug
// iso2d
//
// Demo of finding isocontours via Newton-Raphson method.
// Initializes positions on a grid, and each update applies one
// step of Newton-Raphson.
//
// Process output with:
// unu jhisto -i mip.txt -b 512 512 -min 0 0 -max 63 63 | unu 2op neq - 0 | unu quantize -b 8  -o iso2d.png

int gridSize = 100;
real isoval = 30.0;
image(2)[] img = load("../data/ddro-64-unit.nrrd");
field#1(2)[] F = img ⊛ ctmr;
int stepsMax = 10;
real epsilon = 0.001;

strand sample (int ui, int vi) {
    output vec2 pos = [lerp(1.0, 62.0, -0.5, real(ui), real(gridSize)-0.5),
                       lerp(1.0, 62.0, -0.5, real(vi), real(gridSize)-0.5)];
    int steps = 0;
    update {
        // We bail if we're no longer inside or taken too many steps.
        // We tag this strand result as garbage by setting it to a
        // non-finite position; these will be filtered later.  
        // This is really what "die" is for, though (see below).
        if (!inside(pos, F) || steps > stepsMax) {
          // HEY (BUG) not able to use NaN as in "pos = [nan,nan]"
          // generated C code "vec2f(nanf, nanf)" causes problems:
          // iso2d.c: In function ‘sample_update’:
          // iso2d.c:86: error: incompatible type for argument 1 of ‘vec2f’
          // iso2d.c:86: error: incompatible type for argument 2 of ‘vec2f’
          // uncaught exception Fail [Fail: error compiling/linking]
          //   raised at c-target/c-target.sml:323.14-323.44
          pos = [∞,∞];  // should be: pos = [nan,nan];
          // HEY (BUG) we shouldn't even have to use non-finite values;
          // we should be able to just say "die;".  Currently, using
          // "die" leads to the program never finishing...
          stabilize;
        }
        vec2 grad = ∇F(pos);
        real gmsq = grad • grad;
        // If |∇F|=0 we can't compute update; we have to bail
        if (gmsq == 0.0) {    
          pos = [∞,∞];
          stabilize;
        }
        // delta = Newton-Raphson step
        vec2 delta = -(F(pos) - isoval)*grad/gmsq;
        // We've converged if the update is small enough
        if (|delta| < epsilon) { 
          stabilize;
        }
        pos = pos + delta;
        steps = steps + 1;
    }
}

initially [ sample(ui, vi) | vi in 0..(gridSize-1), ui in 0..(gridSize-1) ];

root@smlnj-gforge.cs.uchicago.edu
ViewVC Help
Powered by ViewVC 1.0.0