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

SCM Repository

[diderot] Diff of /examples/iso2d-spatial/iso2d.diderot
ViewVC logotype

Diff of /examples/iso2d-spatial/iso2d.diderot

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

revision 3012, Mon Mar 9 11:24:46 2015 UTC revision 3013, Mon Mar 9 13:06:51 2015 UTC
# Line 1  Line 1 
1  // iso2d  int gridSize = 200;
2  //  real isoval = 30;
3  // Demo of finding isocontours via Newton-Raphson method.  input real rad = 0.3;      // actual particle radius
4  // Initializes positions on a grid, and each update applies one  real hhInit = rad/4;
5  // step of Newton-Raphson.  real evariance = ∞;
 //  
 // Process output with:  
 // 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  
   
 int gridSize = 100;  
 image(2)[] diderot = image("data/ddro-80.nrrd");  
 field#1(2)[] F = ctmr ⊛ diderot;  
 input int stepsMax = 50;  
 input real stepScale = 1.0;  
 real epsilon = 0.0001;  
 //////////////// Global Variables for Spacing Particles //////////////////////  
 input real rr = 2.0;      // actual particle radius  
 real stdev = ∞;  
 real RR = rr+0.0;         // neighbor query radius (MUST be >= rr; should get same results for any RR >= rr)  
 real hhInit = 10.0;       // initial integration step size; can err too big, will be trimmed down during iterations  
 int iterSpacing = 1;             // which iteration we're on  
 int iterMax = 200;  
 real iosVal50 = 0.0;  
 real isoVal30 = 0.0;  
 real isoVal10 = 0.0;  
 /////////////////////////////////////////////////////////////////////////////  
   
 strand Particle (int ui, int vi) {  
   // world is 1x1 centered at (0.5, 0.5)  
     vec2 pos = [lerp(0.0, 1.0, -0.5, real(ui), real(gridSize)-0.5),  
                        lerp(0.0, 1.0, -0.5, real(vi), real(gridSize)-0.5)];  
   // set the isvalue to 50, 30, or 10, depending on whichever we're closest to  
     real isoval = 50.0 if F(pos) >= 40.0  
                        else 30.0 if F(pos) >= 20.0  
                                  else 10.0;  
   
     real pIsoVal50 = 0.0;  
     real pIsoVal30 = 0.0;  
     real pIsoVal10 = 0.0;  
6    
7      bool foundContour = false;  input int newtonIterMax = 10;
8      int contourSteps = 0;  input int forceIterMax = 20;
9      int pSteps = 0;  
10      real energy = 0;     // HEY: can do convergence test based on variance of energies  field#1(2)[] F = bspln3 ⊛ image("data/ddro-80.nrrd");
11      vec2 force = [0,0];   // or can test convergence based on sum of |force|  
12      output vec2 outPos = pos;  strand Particle (int ID, vec2 pos0) {
13      output vec2 pos = pos0;
14    
15      vec2 posOld1 = pos;   // remember last TWO positions    bool foundContour = false;
16      vec2 posOld2 = pos;    int iter = 0;
17      int forceIter = 0;  // want to get rid of this with global convergence criterion
18      real energy = 0;    // has to be a member variable so that global reduction can see it
19      vec2 force = [0,0];
20      real hh = hhInit;      real hh = hhInit;
21    
22      update {      update {
23        if (!foundContour) {
24          if(foundContour)        if (!inside(pos, F) || iter == newtonIterMax) {
25          {          die; // quit if outside field or took too many steps
26            if( contourSteps > iterMax)  {        }
27          // Newton-Raphson step
28          vec2 step = -normalize(∇F(pos))*(F(pos) - isoval)/|∇F(pos)|;
29          pos += step;
30          iter += 1;
31          if (|step| < 0.001) {
32            foundContour = true;
33          }
34        } else {
35          if (forceIter == forceIterMax) {
36              stabilize;              stabilize;
37            }            }
38          forceIter +=1;
39            // compute energy and forces on us        energy = 0; force = [0,0];
40            energy = 0;        foreach (Particle p_j in sphere(rad)) {
41            force = [0,0];          vec2 r_ij = (pos - p_j.pos)/rad;
   
           foreach (Particle p_j in sphere(RR)) {  
   
             if(p_j.foundContour)  
             {  
               vec2 r_ij = (pos - p_j.pos)/rr;  
               vec2 d_ij = normalize(r_ij);  
               if (|r_ij| < 1) {  
42                  energy += (1 - |r_ij|)^4;                  energy += (1 - |r_ij|)^4;
43                  force += - (-4*(1 - |r_ij|)^3) * d_ij;          force -= normalize(r_ij)*(-4*(1 - |r_ij|)^3)/rad;
44                }                }
45              }        force -= normalize(∇F(pos))⊗normalize(∇F(pos))•force;
46            }        if (|force| > 0) {
47            force /= rr;     // smaller particles make larger forces          // project force onto tangent plane
48            real stepMax = rr;  // iosVal50 if isoval == 5.0          // take gradient descent step
49                        // else isoVal10 if isoval == 30.0          vec2 down = normalize(force);
50                               //    else isoVal10;  // limit on distance to travel per iter          vec2 posLast = pos;
51            pos += hh*down;
52            // update position based on force          // take Newton-Raphson step back to surface
53            posOld2 = posOld1;   // shuffle saved positions down          pos += -((F(pos) - isoval)/|∇F(pos)|)*normalize(∇F(pos));
54            posOld1 = pos;          pos += -((F(pos) - isoval)/|∇F(pos)|)*normalize(∇F(pos));
55            if (energy > 0.0) {  // we have neighbors          real newenergy = 0;
56                tensor[2,2] pten = identity[2] - normalize(pos)⊗normalize(pos);          foreach (Particle p_j2 in sphere(rad)) {
57                force = pten•force;  // project force onto tangent surface            vec2 r_ij = (pos - p_j2.pos)/rad;
58                vec2 step = hh*force;            newenergy += (1 - |r_ij|)^4;
59                if (|step| > stepMax) {          }
60                    // decrease hh by factor by which step was too big          if (newenergy > energy - 0.5*hh*|force|) {
61                    hh *= stepMax/|step|;            hh *= 0.5; // backtrack
62                    // and find smaller step            pos = posLast;
                   step = hh*force;  
               }  
               // take step and re-find implicit surface  
               pos = normalize(pos + step);  
               real travel = |pos - posOld1| + |posOld1 - posOld2|;  
               if (travel > 0) {  
                   // if we've moved in the past two steps, but we've moved back  
                   // to where we were two steps ago, we're oscillating ==>  
                   // okay = 0. Two steps in the same direction ==> okay = 1.  
                   real okay = |pos - posOld2|/travel;  
                   // slow down if oscillating, speed up a little if not  
                   hh *= lerp(0.8, 1.001, 0, okay, 1);  
               }  
           }  
           contourSteps+=1;  
63          }else {          }else {
64              // We bail if we're no longer inside or taken too many steps.            // re-initialize for new particle neighborhood
65              if (!inside(pos, F) || pSteps > stepsMax) {            hh = hhInit;
                  die;  
             }  
             vec2 grad = ∇F(pos);  
             if (|grad| == 0.0) {    // can't compute step if |∇F|, so have to bail  
                  die;  
             }  
             vec2 norm = normalize(grad);  
             vec2 delta = -((F(pos) - isoval)/|grad|)*norm;  // Newton-Raphson step  
             if (|delta| < epsilon) {    // we've converged if step is small enough  
                    foundContour = true;  
                  print("converged\n");  
                  if(isoval == 50) {  
                     pIsoVal50 = 1.0;  
                  }else if (isoval == 30) {  
                     pIsoVal30 = 1.0;  
                  }else {  
                     pIsoVal10 = 1.0;  
66                   }                   }
67              }              }
68              pos += stepScale * delta;      } // if (!foundContour)
             pSteps += 1;  
         }  
         outPos = pos;  
69      }      }
70  }  }
71  global{  global{
72     real energyMean = mean{P.energy | P in Particle.all};     real energyMean = mean{P.energy | P in Particle.all};
73     real var = mean{ (P.energy - energyMean) * (P.energy - energyMean) | P in Particle.all};    evariance = mean{ (P.energy - energyMean) * (P.energy - energyMean) | P in Particle.all};
    iosVal50 = sum{P.pIsoVal50 | P in Particle.all};  
    isoVal30 = sum{P.pIsoVal30 | P in Particle.all};  
    isoVal10 = sum{P.pIsoVal10 | P in Particle.all};  
   
  // real var = variance{P.energy | P in Particle.all};  
   stdev = sqrt(var);  
74  }  }
75    
 initially { Particle(ui, vi) | vi in 0..(gridSize-1), ui in 0..(gridSize-1) };  
76    initially { Particle(ui + gridSize*vi,
77                         [lerp(-1.5, 1.5, -0.5, real(ui), real(gridSize)-0.5),
78                          lerp(-1.5, 1.5, -0.5, real(vi), real(gridSize)-0.5)])
79                 | vi in 0..(gridSize-1), ui in 0..(gridSize-1) };

Legend:
Removed from v.3012  
changed lines
  Added in v.3013

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