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

SCM Repository

[diderot] View of /tests/lamont-tests/halftone.diderot
ViewVC logotype

View of /tests/lamont-tests/halftone.diderot

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4640 - (download) (annotate)
Tue Sep 27 20:54:47 2016 UTC (2 years, 11 months ago) by glk
File size: 3827 byte(s)
initial result of svn export --username anonsvn --password=anonsvn https://svn.smlnj-gforge.cs.uchicago.edu/svn/diderot/branches/vis15/src/tests/
vec2{} initPosns = load("data/pos.nrrd");
int numOfParticles = length(initPosns);
input int iterMax = 200;  // maximum number of steps to run
field#1(2)[] img = ctmr ⊛ image("data/ddro.nrrd");
real heat = 0.25;
real rr = 0.02;
real alpha = 0.5;
real threshold = 0.0;       // convergence threshold
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 iter = 1;             // which iteration we're on

// ============================ begin setting up gridding of domain
vec2 xDom = [-0.2,1.2];
vec2 yDom = [-0.2,1.2];
real xSamples = floor((xDom[1] - xDom[0])/RR); // Lamont verify logic floor vs ceil here
real ySamples = floor((yDom[1] - yDom[0])/RR);
vec4 qWinDim = [xDom[0],xDom[1],yDom[0],yDom[1]]; // i.e. [XMIN, XMAX, YMIN, YMAX]  (required for the query function)
vec2 qGridDim = [xSamples,ySamples];   // how many grid cells you want in each direction for the uniform grid (required for the query function)
vec2 qCellDim = [(xDom[1] - xDom[0])/xSamples,(yDom[1] - yDom[0])/ySamples];      // the length in each direction for a cell (required for the query function)
// ============================ end setting up gridding of domain

function vec2 engfrc(real r) {
  return [(1 - r)^4, -4*(1 - r)^3];
}

strand Particle (int ii, vec2 pos0) {
    vec2 pos = pos0;
    real hh = hhInit;
    vec2 posOld1 = pos;   // remember last TWO positions
    vec2 posOld2 = pos;
    output vec2 outPos = pos;
    real stepMax = rr/2;
    real ww = 0.0;
    real penergy = ii;
    vec2 pforce = [0,0];
    real ienergy = ii;
    vec2 iforce = [0,0];
    update {
        if (iter >= iterMax)  {
            stabilize;
        }
        if (!inside(pos, img)) {
            die;
        }
        penergy = 0;
        pforce = [0,0];
        ienergy = img(pos);
        iforce = -∇img(pos);
        foreach (Particle p_j in sphere(rr)) {
            vec2 r_ij = (pos - p_j.pos)/rr;
            vec2 d_ij = normalize(r_ij);
            if (|r_ij| < 1) {
              // BUG: this line works ...
              //vec2 ef = [(1 - |r_ij|)^4, -4*(1 - |r_ij|)^3];
              // .. but this line doesn't
              vec2 ef = engfrc(|r_ij|);
              penergy += ef[0];
              pforce += - ef[1] * d_ij;
            }
        }
        pforce /= rr;     // smaller particles make larger forces

        real energy = (1-alpha)*penergy + alpha*ienergy;
        vec2 force = (1-alpha)*pforce + alpha*iforce;
        // update position based on force
        posOld2 = posOld1;   // shuffle saved positions down
        posOld1 = pos;
        if (energy > 0.0) {  // we have neighbors
            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
            pos += step;
            real travel = |pos - posOld1| + |posOld1 - posOld2|;
            if (travel > 0) {
                real okay = |pos - posOld2|/travel;
                hh *= lerp(0.9, 1 + heat, 0, okay, 1);
            }
        }
        pos = [clamp(0, 1, pos[0]), clamp(0, 1, pos[1])];
        outPos = pos;
    }
}

// can add statements in here to do some global computation
global{
  iter+=1;
  real isum = sum{P.ienergy | P in Particle.active};
  real psum = sum{P.penergy | P in Particle.active};
  print("isum = ", isum, "; psum = ", psum, "\n");
  real aa = isum/(psum+isum);
  alpha = (alpha + aa)/2;
  print("----> alpha = ", alpha, "\n");
  heat *= 0.97;
}

initially {Particle(ii, initPosns{ii})
           | ii in 0 .. numOfParticles-1 };


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