Home My Page Projects Code Snippets Project Openings diderot

# SCM Repository

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

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

Tue Sep 27 20:54:47 2016 UTC (2 years, 9 months ago) by glk
File size: 3827 byte(s)
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 };