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

SCM Repository

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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2950 - (download) (annotate)
Wed Mar 4 20:52:35 2015 UTC (4 years, 3 months ago) by lamonts
File size: 5583 byte(s)
Adding the iso2d demo with particle interaction
// 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 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; 

    bool foundContour = false; 
    int contourSteps = 0; 
    int pSteps = 0; 
    real energy = 0;     // HEY: can do convergence test based on variance of energies
    vec2 force = [0,0];   // or can test convergence based on sum of |force|
    output vec2 outPos = pos; 

    vec2 posOld1 = pos;   // remember last TWO positions
    vec2 posOld2 = pos;
    real hh = hhInit;

    update {

        if(foundContour) 
        {
          if( contourSteps > iterMax)  {
            stabilize;
          } 

          // compute energy and forces on us
          energy = 0;
          force = [0,0];
          
          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) {
                energy += (1 - |r_ij|)^4;
                force += - (-4*(1 - |r_ij|)^3) * d_ij;
              }
            }
          }
          force /= rr;     // smaller particles make larger forces
          real stepMax = rr;  // iosVal50 if isoval == 5.0 
                      // else isoVal10 if isoval == 30.0
                             //    else isoVal10;  // limit on distance to travel per iter

          // update position based on force
          posOld2 = posOld1;   // shuffle saved positions down
          posOld1 = pos;
          if (energy > 0.0) {  // we have neighbors
              tensor[2,2] pten = identity[2] - normalize(pos)⊗normalize(pos);
              force = pten•force;  // project force onto tangent surface
              vec2 step = hh*force;
              if (|step| > stepMax) {
                  // decrease hh by factor by which step was too big
                  hh *= stepMax/|step|;
                  // and find smaller step
                  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; 
        }else {
            // We bail if we're no longer inside or taken too many steps.
            if (!inside(pos, F) || pSteps > stepsMax) {
    	         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;
                 }
            }
            pos += stepScale * delta;
            pSteps += 1;
        }
        outPos = pos; 
    }
}
global{
   real energyMean = mean{P.energy | P in Particle.all};
   real var = 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);
}

initially { Particle(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