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

SCM Repository

[diderot] View of /tests/vis15-bugs/halftone-bug2.diderot
ViewVC logotype

View of /tests/vis15-bugs/halftone-bug2.diderot

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4929 - (download) (annotate)
Thu Feb 16 20:17:03 2017 UTC (2 years, 5 months ago) by jhr
File size: 7865 byte(s)
  remove svnexecutable property from files that shouldn't be executable
/* ==========================================
Half-toning via interacting repulsive particles

This example is based on the [`sphere.diderot`](../sphere) example,
which in turn is based on the [`circle.diderot`](../circle) example;
both of these should be read first.  While in those cases a surface
(circle or sphere) is sampled uniformly, in this case the flat surface
of an image is sampled non-uniformly, according to the image
intensity.

... TODO: example of running program and looking at output

========================================== */

input image(2)[] img ("data in which to find features") = image("halftone-bug2-img.nrrd");
field#0(2)[] F = ctmr ⊛ clamp(img);

input vec2{} ipos ("initial positions for all particles") = load("halftone-bug2-vec2.nrrd");
input vec2 radmm ("particle minimum > 0 and maximum radius") = [0.01, 0.1];
input real eps ("system convergence threshold, computed as the coefficient-of-variation of normalized distances to nearest neighbors") = 0.05;
input int pcp ("periodicity of population control (if greater than zero)") = 2;
input real hhInit ("initial step size for potential energy gradient descent") = 1;

real newDist = ((0.666));  // how far away to birth, as multiple of radius
real stepMax = ((0.666));  // limit on travel per iter, as multiple of radius
int iter = 0;         // which iteration we're on

// potential function phi(r) found via Mathematica: this is the lowest-degree
// piecewise polynomial with phi(0)=1 and a potential well at phi(0.666)=-0.005.
// phi(r) is C^3 continuous across the well and with 0 at r=1.
function real phi(real r) =
1 + r*(-4.248582222661734 + r*(5.991287388535527 + (-2.864766048887071 + 0.06790595504541913*r)*r))
if r < 0.666 else
-0.005 + ((0.4482053856358993 + (-2.6838645846460842 + (6.026642031390917 - 4.811690244623482*(-0.666 + r))*(-0.666 + r))*(-0.666 + r))*(-0.666 + r))*(-0.666 + r)
if r < 1 else 0;
function real phi'(real r) =
-4.248582222661734 + r*(11.982574777071054 + (-8.594298146661213 + 0.2716238201816765*r)*r)
if r < 0.666 else
-15.42591894092919 + 0.4482053856358993*(-1.332 + 2*r) + r*(-2.6838645846460842*(-3.996 + 3*r) + 6.026642031390917*(5.322672 + r*(-7.992 + 4*r)) - 4.811690244623482*(-5.90816592 + r*(13.30668 + r*(-13.32 + 5*r))))
if r < 1 else 0;

// Energy and force from particle (with radius rad) at vec2 x
function real enr(vec2 x, real rad) = rad*rad*phi(|x|/rad);
function vec2 frc(vec2 x, real rad) = rad*rad*phi'(|x|/rad) * (1/rad) * x/|x|;

// From a given vec2, find a random-ish value uniformly sampling [0,1)
function real posrnd(vec2 v, real rad) {
   vec2 p = 10000*v/rad;
   return fmod(|fmod(p[0],1)| + |fmod(p[1],1)|, 1);
}

// Is this an iteration in which to do population control?
function bool pcIter() = (pcp > 0 && iter > 0 && 0 == iter % pcp);
function bool birth() = ( (iter/pcp) % 2 == 0);

function real radius(real ff) =
   min(radmm[1], radmm[0]*sqrt(1/clamp(0.0000001, 1, ff)));


/* The particle is initialized at position pos0, with initial stepsize hh0.
   The first set of particles gets hhInit for the initial stepsize, but new
   particles created by population control benefit from getting the stepsize
   that was adaptively learned by the parent. */
strand particle (int gen_, vec2 pos0, real hh0) {
   int gen = gen_;
   vec2 pos = pos0;
   output vec3 posn = [pos[0], pos[1], 3];
   real hh = hh0;
   vec2 step = [0,0];   // step along force
   real rad = 0;        // zero isn't a valid value
   real closest = 1;    // distance to closest neighbor (as mult of radius)
   int ncount = 0;      // how many neighbors did we have
   real mvmt = 1;
   real energy=0;
   update {
      if (!inside(pos, F)) { die; }                  
      rad = radius(F(pos));
      // compute energy and forces on us from neighbors
      real energyLast = 0;
      vec2 force = [0,0];
      ncount = 0;
      foreach (particle P in sphere(2*rad)) {
         vec2 x = P.pos - pos;
         real rr = radius(F(lerp(pos, P.pos, 0.5)));
         energyLast += enr(x, rr);
         force += frc(x, rr);
         ncount += 1 if |x| < rr else 0;
      }

      /* Else have interacting neighbors; find step, limit step size */
      step = hh*force;
      if (|step| > stepMax*rad) {
         hh *= stepMax*rad/|step|;
         step = hh*force;
      }
      mvmt = lerp(mvmt, |step|/(rad*((0.666))), 0.5);

      vec2 posLast = pos;
      pos += step;
      energy = 0;
      closest = 1;
      ncount = 0;
      vec2 mno = [0,0]; // mean neighbor offset
      if (|pos| < 0.025 && (iter==5 || iter==6)) {
         print("\n(iter ", iter, ") HELLO (gen ", gen, ") at pos ", pos, " (|pos|=",
		|pos|, "), looking for neighbors with base radius ", rad, "\n");
      }
      foreach (particle P in sphere(2*rad)) {
         real rr = radius(F(lerp(pos, P.pos, 0.5)));
         vec2 x = P.pos - pos;
//         if (|pos| < 0.025 && (iter==5 || iter==6)) {
//            print("    neighbor at ", P.pos, " ==> halfway ", lerp(pos, P.pos, 0.5),
//		" F=", F(lerp(pos, P.pos, 0.5)), " ==> rr=", rr, "; |x|=", |x|,
//		"<" if |x| < rr else ">=", rr, "\n");
//         }
         if (|x| < rr) {
            energy += enr(x, rr);
            closest = min(closest, |x|/rr);
            mno += x;
            ncount += 1;
            if (|pos| < 0.025 && (iter==5 || iter==6)) {
               print("\n        REAL neighbor P at ", P.pos, " (|P.pos|==", |P.pos|,
		"==" if |P.pos|==0 else "!=", "0; P.gen ", P.gen, "); now, energy=",
		energy, ", closest=", closest, "; ncount=", ncount, "\n");
            }
         }
      }
      mno /= ncount;

      // Armijo-Goldstein sufficient decrease condition
      if (energy - energyLast > 0.5*(pos - posLast)•(-force)) {
         hh *= 0.5;     // energy didn't decrease as expected: backtrack
         pos = posLast; // try again next iteration
      } else {
         hh *= 1.05; // opportunistically increase step size
         // try to have between 6 and 8 neighbors
         if (pcIter()) {
            if (birth()) {
               if (energy < 0 && posrnd(pos, rad) < (6.0 - ncount)/6.0) {
                  real a = atan2(-mno[1],-mno[0]);
                  a += lerp(-1, 1, gen % 2)*π/6.0;
                  vec2 npos = pos + newDist*rad*[cos(a),sin(a)];
                  print("\n(iter ", iter, ") NEW particle(", gen+1, ",", npos, ",", hh, ")\n");
                  new particle(gen+1, npos, hh);
               }
            } else if (energy > 0 && posrnd(pos, rad) < (ncount - 8.0)/ncount) {
               if (|pos| < 0.1 && (iter==5 || iter==6)) {
                  print("(iter ", iter, ") DIE at pos ", pos, " (|pos|=", |pos|, "): energy ",
			energy, "; ncount ", ncount, "; rad ", rad, "\n");
                  print("     ", posrnd(pos, rad), " < ", (ncount - 8.0)/ncount, "\n");
               }
               die;
            }
         }
      }

      posn = [pos[0], pos[1], ncount];
      if (|pos| == 0) {
         print("(iter ", iter, ") WHOA pos exactly at 0!\n");
      }

      // Record actual step taken, in case used in global update
      step = pos - posLast;
   }
}

global {
   /* Compute coefficient-of-variation of distance to closest neighbor */
   real totenr = sum { P.energy | P in particle.all};
   real meancl = mean { P.closest | P in particle.all};
   real varicl = mean { (P.closest - meancl)^2 | P in particle.all};
   real covcl = sqrt(varicl) / meancl;
   print("\n=-=-=\n=-=-= iter ", iter,
         "; tot#=", numActive(),
         "; totenr=", totenr,
         "; COV(cl)=", covcl,
         "; mean(cl)=", meancl,
         "; mean(hh)=", mean { P.hh | P in particle.all},
         "; mean(ncount)=", mean { real(P.ncount) | P in particle.all},
         "; mean(mvmt)=", mean { P.mvmt | P in particle.all },
         "\n");
   iter += 1;
}

initially { particle(2, ipos[ii], hhInit) | ii in 0 .. length(ipos)-1 };

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