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

SCM Repository

[diderot] Diff of /branches/lamont/test/iso2d.diderot
ViewVC logotype

Diff of /branches/lamont/test/iso2d.diderot

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 2306, Thu Mar 21 08:03:03 2013 UTC revision 2307, Thu Mar 21 09:06:39 2013 UTC
# Line 5  Line 5 
5  // step of Newton-Raphson.  // step of Newton-Raphson.
6  //  //
7  // Process output with:  // Process output with:
8  // unu jhisto -i iso2d.txt -b 512 512 -min 0 0 -max 1 1 | unu 2op neq - 0 | unu quantize -b 8  -o iso2d.png  // unu jhisto -i pos.nrrd -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;  int gridSize = 300;
11  field#1(2)[] F = ctmr ⊛ image("../data/ddro-80.nrrd");  field#1(2)[] F = ctmr ⊛ image("../data/ddro-80.nrrd");
# Line 19  Line 19 
19      real isoval = 50.0 if F(pos) >= 40.0      real isoval = 50.0 if F(pos) >= 40.0
20                         else 30.0 if F(pos) >= 20.0                         else 30.0 if F(pos) >= 20.0
21                                   else 10.0;                                   else 10.0;
22    // the following won't work, because the compiler doesn't know how to deal with state vars that
23    // are fields (really just a local binding)!
24    //    field#1(2)[] G = F - isoval;
25      int steps = 0;      int steps = 0;
26      update {      update {
27            field#1(2)[] G = F - isoval;
28          // We bail if we're no longer inside or taken too many steps.          // We bail if we're no longer inside or taken too many steps.
29          if (!inside(pos, F) || steps > stepsMax) {          if (!inside(pos, G) || steps > stepsMax) {
30              die;              die;
31          }          }
32          // GLKs recent changes (revision 1329) were made to make the code          // GLKs recent changes (revision 1329) were made to make the code
33          // a more obvious implementation of Newton-Raphson, and also to create          // a more obvious implementation of Newton-Raphson, and also to create
34          // some obvious and not-so-obvious opportunties for optimization          // some obvious and not-so-obvious opportunties for optimization
35          vec2 grad = ∇F(pos);          vec2 grad = ∇G(pos);
36          if (|grad| == 0.0) {    // can't compute step if |∇F|, so have to bail          if (|grad| == 0.0) {    // can't compute step if |∇F|, so have to bail
37              die;              die;
38          }          }
39          vec2 norm = normalize(grad);          vec2 norm = normalize(grad);
40          vec2 delta = -((F(pos) - isoval)/|grad|)*norm;  // Newton-Raphson step          vec2 delta = -(G(pos)/|grad|)*norm;  // Newton-Raphson step
41          if (|delta| < epsilon) {    // we've converged if step is small enough          if (|delta| < epsilon) {    // we've converged if step is small enough
42              stabilize;              stabilize;
43          }          }

Legend:
Removed from v.2306  
changed lines
  Added in v.2307

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