halftone.diderot

#version 1.0
/* ==========================================
## halftone.diderot: 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 and understood 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.  This creates a
result similar to the [Electrostatic
Halftoning](http://www.mia.uni-saarland.de/Research/Electrostatic_Halftoning/index.shtml)
method of Schmaltz et al. The very regular hexagonal grid patterns that this
program generates by energy minimization in uniform areas may be sub-optimal
for artistic or signal-processing purposes; this program lacks the jittering
that better haltoning methods add to avoid this.  The regular sampling
generated by this program could help subsequent meshing.

The standard test for this kind of program is a linear ramp, which is
available via:

	#!
	#I due to non-deterministic parallelism and FP sensitivity
	../fs2d/fs2d-scl -which x -width 2 -size0 401 -size1 401 |
	  unu crop -min 0 100 -max M M-100 |
	  unu affine -1 - 1 0 1 |
	  unu pad -min -2 -2 -max M+2 M+2 -o img.nrrd
	rm -f out.nrrd
	#_ junk img.nrrd

The domain of this `img.nrrd` is [-1,1] along X, and [-0.5,0.5] along Y, and
which varies from 0 at X=-1 to 1 at X=1. Two extra samples are added at
the boundaries to ensure that even at the boundary, the `inside()` test
passes. With this proxy image in place, we can compile:

	diderotc --snapshot --exec halftone.diderot
	#prog halftone.diderot

We use snapshots to monitor the progress of computation.  To make `NN` initial
positions with random number seen `RNG` that fit within the ramp image domain:

	NN=300
	RNG=5
	echo 0 0 | unu pad -min 0 0 -max M $((NN-1)) |
	  unu 1op rand -s $RNG | unu affine 0 - 1 -1 1 -o vec2.nrrd
	echo 1 0.5 | unu 2op x vec2.nrrd - -o vec2.nrrd
	#_ junk vec2.nrrd

Then to run with snapshots saved every iteration (`-s 1`), but limiting the program
to 800 iterations (`-l 800`), as well as cleaning results from previous run:

	rm -f pos-????.{png,nrrd} pos.nrrd
	#^ ./halftone -s 2 -l 800 -radmm 0.04 1 -eps 0.0001 -pcp 2
	#_ if true; then  # (begin possible block comment)
	#_ echo "##### FIRST (smaller, faster) TEST; should converge within -l limit"
	#_ ./halftone -s 0 -l 800 -radmm 0.04 1 -eps 0.0001 -pcp 2

Running with this large value (0.04) of minimum radius (considering the image domain)
is good for giving a visual impression of how the particle system populates the
domain. Next, some `unu` hacking makes images of the evolving system, which
are then turned into an animated `ramp.gif` (compare to `[ramp-ref.gif](ramp-ref.gif)).

	SZ=200
	OV=2
	export NRRD_STATE_VERBOSE_IO=0
	#^ for PIIN in pos-????.nrrd; do
	#^    IIN=${PIIN#*-}; II=${IIN%.*}
	#_ for PIIN in pos.nrrd; do
	#_    II=lores
	   echo "post-processing $PIIN to pos-$II.png ... "
	   unu jhisto -i $PIIN -min -1 -0.5 -max 1 0.5 -b $((OV*SZ*2)) $((SZ*OV)) |
	     unu resample -s /$OV /$OV -k bspln5 -t float |
	     unu quantize -b 8 -min 0 -max $(echo "0.15 / ($OV * $OV)" | bc -l) -o pos-$II.png
	done
	#^ convert -delay 6 pos-*.png ramp.gif
	#_ # with a tolerance of 256 we're saying "anything goes" but the reason
	#_ # is still to generate some way of looking at how the system ended up
	#> pos-lores.png 256
	#_ # lo-res and blurred versions histograms of x and y positions
	#_ unu slice -i pos.nrrd -a 0 -p 1 |
	#_   unu histo -min -0.5 -max 0.5 -b 50 -t float |
	#_   unu resample -s x1 -k gauss:4,3 -b mirror -o pos-lores-yhisto.nrrd
	#_ unu slice -i pos.nrrd -a 0 -p 0 |
	#_   unu histo -min -1 -max 1 -b 100 -t float |
	#_   unu resample -s x1 -k gauss:4,3 -b mirror -o pos-lores-xhisto.nrrd
	#> pos-lores-?histo.nrrd 1.8
	#_ fi  # (end possible block comment)

On the other hand, to quantitatively check that the particle density is as it should be,
we run with many more particles (which can take a few minutes to finish). The following
uses `-radmm 0.004 1` (versus `-radmm 0.04 1` above) for the particle computation,
so that inter-particle distance may be a tenth of what it was previously (resulting in
up to 100 times greater particle density).  A histogram of the X positions, plotted
above the rasterized image of the particle positions, provides visual confirmation that
horizontal (X position) particle density approaches the linear ramp expected from the
linear ramp of the underlying image.

	#_ if true; then  # (begin possible block comment)
	rm -f {hp,pos}-????.{png,nrrd} pos.nrrd
	#_ NN=30000
	#_ RNG=5
	#_ echo 0 0 | unu pad -min 0 0 -max M $((NN-1)) |
	#_   unu 1op rand -s $RNG | unu affine 0 - 1 -1 1 -o vec2.nrrd
	#_ echo 1 0.5 | unu 2op x vec2.nrrd - -o vec2.nrrd
	#_ echo "##### SECOND (bigger, slower) TEST (will not converge) within -l limit"
	#||:
	#_ ./halftone -s 50 -l 150 -radmm 0.006 1 -eps 0.00004 -pcp 1
	#^ ./halftone -s 10 -l 800 -radmm 0.004 1 -eps 0.00004 -pcp 2
	SZ=200
	OV=2
	export NRRD_STATE_VERBOSE_IO=0
	#^ for PIIN in pos-????.nrrd; do
	#_ for PIIN in pos-0{000,050,100,150}.nrrd; do
	#_   if [[ ! -e $PIIN ]]; then
	#_     echo "$0: expected output file $PIIN does not exist, skipping"
	#_     continue;
	#_   fi
	   IIN=${PIIN#*-}; II=${IIN%.*}
	   echo "post-processing $PIIN to pos-$II.png ... "
	   unu jhisto -i $PIIN -min -1 -0.5 -max 1 0.5 -b $((OV*SZ*2)) $((SZ*OV)) |
	     unu resample -s /$OV /$OV -k bspln3 -t float |
	#^     unu quantize -b 8 -min 0 -max $(echo "2 / ($OV * $OV)" | bc -l) -o pos-$II.png
	#_     unu quantize -b 8 -min 0 -max $(echo "1 / ($OV * $OV)" | bc -l) -o pos-$II.png
	   unu slice -i $PIIN -a 0 -p 0 |
	     unu histo -min -1 -max 1 -b $((SZ/3)) |
	     unu dhisto -h $((SZ/3)) -nolog |
	     unu resample -s $((SZ*2)) = -k box |
	     unu join -i - pos-$II.png -a 1 -o hp-$II.png # combine histogram and position image
	done
	#^ convert -delay 6 hp-*.png hp.gif
	#_
	#_ mv hp-0150.png histo-pos-hires.png
	#_ # like above; tolerance 256 says "anything goes; may be informative"
	#> histo-pos-hires.png 256
	#_ junk hp-0{000,050,100}.png pos-0{000,050,100,150}.{png,nrrd}
	#_ unu slice -i pos-0150.nrrd -a 0 -p 1 |
	#_   unu histo -min -0.5 -max 0.5 -b 100 -t float |
	#_   unu resample -s x1 -k gauss:4,3 -b mirror -o pos-hires-yhisto.nrrd
	#_ unu slice -i pos-0150.nrrd -a 0 -p 0 |
	#_   unu histo -min -1 -max 1 -b 200 -t float |
	#_   unu resample -s x1 -k gauss:4,3 -b mirror -o pos-hires-xhisto.nrrd
	#> pos-hires-?histo.nrrd 4.5
	#_ fi  # (end possible block comment)

The top part of the resulting `hp.gif` (compare to [`hp-ref.gif`](hp-ref.gif)) shows the
expected convergence to a linear variation in particle density, even though there are
too many particles to individually distinguish in the rasterized image.

Finally, to have some fun with a picture of [Diderot
himself](https://en.wikipedia.org/wiki/Denis_Diderot), we invert the
image intensity and darken a bit (tending to have more space between
particles), recompile with the new proxy image (the array axis sizes changed),
and then run with new initial positions (without snapshots this time):

	#R   not clear that this reveals anything new
	unu 2op - 1 ../data/ddro.nrrd | unu gamma -g 0.75 -o img.nrrd
	diderotc --exec halftone.diderot
	NN=1000
	RNG=5
	echo 0 0 | unu pad -min 0 0 -max M $((NN-1)) |
	  unu 1op rand -s $RNG | unu affine 0 - 1 -1 1 -o vec2.nrrd
	./halftone -l 800 -radmm 0.01 0.2 -eps 0.000009 -pcp 1
	#_ unu jhisto -i pos.nrrd -min -1 -1 -max 1 1 -b 100 100 -o ddro-jhisto.nrrd
	#> ddro-jhisto.nrrd 0
	#^ echo 1 -1 |
	#^   unu 2op x pos.nrrd - |
	#^   unu save -f text |
	#^   sed -e s/$/\ 0.005\ 0\ 360\ arc\ closepath\ fill/ |
	#^   cat head.eps - tail.eps > ddro.eps
	#^ epstopdf ddro.eps

The `ddro.eps` is generated with the help of `head.eps` and `tail.eps`,
which assume that the points live inside the domain [-1,1]x[-1,1].
The final command for generating `ddro.pdf` is the
[`epstopdf`](https://www.ctan.org/pkg/epstopdf?lang=en) that comes
with some LaTeX installations, but there are many other ways
to convert from eps to pdf. The `ddro.pdf` result (compare to [`drro-ref.pdf`](ddro-ref.pdf))
should look something like the `ddro.png` produced by

	#R
	unu quantize -b 8 -i ../data/ddro.nrrd -o ddro.png

Note that the minimum particle radius specified to `-radmm` was `0.01`, hence
at their tightest packing particles were 0.01 away from each other (they sat
in each other's potential wells, at radius 0.01). Drawing each point with a
circle of radius 0.005 (via `sed` above) creates, where the field is highest
(where `ddro.nrrd` is darkest), a dense packing of mutually tangent circles,
which is visible in the output.

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

input image(2)[] img ("field in which to disperse particles") = image("img.nrrd");
field#0(2)[] F = tent ⊛ clamp(img);

input vec2{} ipos ("initial positions for all particles") = load("vec2.nrrd");
input vec2 radmm ("particle minimum > 0 and maximum radius") = [0.01, 0.1];
input real eps ("system convergence threshold, computed as mean movement (normalized by particle radius)") = 0.005;
input int pcp ("periodicity of population control (if greater than zero)") = 2;
input real hhInit ("initial step size for potential energy gradient descent") = 1;
input real searchScl ("scaling of current radius to get neighbor search radius") = 2.5;

/* A potential function phi(r) found via Mathematica (see findingPhi.nb in
   this directory): this is the lowest-degree piecewise polynomial phi(r) with
   phi(0)=1 and a potential well at phi(0.666)=-0.005, with C^3 continuity
   across the well and with 0 for r >= 1. Having a slight potential well is
   essential for the strategy, used by this program, of using inter-particle
   energy to create the desired packing density, rather than the interaction
   of particles and some ambient image-based potential field. */
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;
real phiWellRad = 0.666;

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

/* Energy and force from particle (with radius rad) at vec2 x. Unlike the
   sphere and circle examples, this takes the radius as an argument. */
function real enr(vec2 x, real rad) = phi(|x|/rad);
function vec2 frc(vec2 x, real 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?  If not, pcIter()
   returns 0. Otherwise, it returns 1 if we should birth new particles, and -1
   if we should kill off particles. This alternation is not because of any
   language limitations; it is just a strategy for aiding the population
   control heuristics. */
function int pcIter() = ((iter/pcp)%2)*2 - 1
                        if (pcp > 0 && iter > 0 && 0 == iter % pcp)
                        else 0;

/* Map from image intensity ff, assummed to be in [0,1], to effective particle
   radius. The principle is that doubling intensity should double sampling
   density, so radius should scale by sqrt(1/2). If ff is in [0,1], sqrt(1/ff)
   is minimized by ff=1; this is the where particles should be closest, so we
   scale by user-specified minimum radius radmm[0], and then clamp by maximum
   desired radius radmm[1]. Finally, we account for how a particle of radius
   1.0 actually has minimum-energy neighbors at phiWellRad. Other image
   effects (like gamma), or a more sophisticated accounting of how dots of
   different densities will produce an effective grayscale, can be
   incorporated into this radius function. */
function real radius(real ff) =
   min(radmm[1], radmm[0]*sqrt(1/clamp(0.000000001, 1, ff)))/phiWellRad;

/* 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) {
   output vec2 pos = pos0;
   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 {
      /* Limit the system to only work inside the image domain. */
      if (!inside(pos, F)) {
         die;
      }
      // rad is nominal radius for current position
      rad = radius(F(pos));

      /* compute energy and forces on us from neighbors. One challenge is that
         even though our radius might be small, which would imply that we don't
         have to look far for interacting neighbors, some of those neighbors
         may have a larger radius, and so we will find them only with a larger
         search radius. The searchScl factor anticipates this. If it is too
         high we are needlessly querying non-interacting neighbors; too low
         and we could be missing some neighbors, which impedes convergence. */
      real energyLast = 0;
      vec2 force = [0,0];
      ncount = 0;
      foreach (particle P in sphere(searchScl*rad)) {
         vec2 x = P.pos - pos;
         if (|x| == 0) {
            die; // same rationale as in ../sphere
         }
         /* Particle interactions must be symmetric (force on A from B must be
            opposite force on B from A).  When particles can vary in radius,
            not enforcing this symmetry can create non-physical
            perpetual motion flows from small to large radii. Here we enforce
            symmetry by learning radius from the average position, which is
            the same strategy as used in the tensor glyph packing (c.f.
            Sect 3.2 of http://people.cs.uchicago.edu/~glk/pubs/#VIS-2006 */
         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;
      }
      if (0 == ncount) {
         if (pcIter() > 0) {
            /* no neighbors, so let's make one, but in a random direction */
            real a = 2*π*posrnd(pos, rad);
            vec2 npos = pos + newDist*rad*[cos(a),sin(a)];
            new particle(gen+1, npos, hh);
         }
         // set closest to something in case used in global update
         closest = newDist;
         continue;
      }
      /* Else we have interacting neighbors; find step, limit step size */
      step = hh*force;
      if (|step| > stepMax*rad) {
         hh *= stepMax*rad/|step|;
         step = hh*force;
      }

      // Take step, find new energy
      vec2 posLast = pos;
      pos += step;
      energy = 0;
      closest = 1;
      ncount = 0;
      /* we compute a mean neighbor offset mno because, with the negative
         potential well, the direction of force is (unlike in the ../sphere
         example) a reliable indicator of where neighbors are more dense */
      vec2 mno = [0,0];
      foreach (particle P in sphere(searchScl*rad)) {
         real rr = radius(F(lerp(pos, P.pos, 0.5)));
         vec2 x = P.pos - pos;
         if (|x| < rr) {
            energy += enr(x, rr);
            closest = min(closest, |x|/rr);
            mno += x;
            ncount += 1;
         }
      }
      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.1;     // opportunistically increase step size
         /* Try to have between 6 and 8 neighbors; see comments in ../sphere
            for description of logic of probabilistic population control */
         if (pcIter() > 0) {
            if (energy < 0 && ncount < 6
                && posrnd(pos, rad) < (6.0 - ncount)/6.0) {
               real a = atan2(-mno[1],-mno[0]);
               a += lerp(-1, 1, gen % 2)*π/6.0;
               vec2 noff = newDist*rad*[cos(a),sin(a)];
               vec2 npos = pos + noff;
               /* make sure that not only is new position inside field (or
                  else it will be immediately killed off) but that past the
                  new position is also inside the field, to prevent a stream
                  of shortly-lived points falling off the edge of the field */
               if (inside(npos, F) && inside(pos + 2*noff, F)) {
                  /* try to avoid landing on an existing point */
                  int nnc = 0;
                  // note that sphere() can take center as first arg
                  foreach (particle P in sphere(npos, rad/100)) { nnc += 1; }
                  if (0 == nnc) {
                     /* no existing particle really close to new position */
                     new particle(gen+1, npos, hh);
                  }
               }
            }
         } else if (pcIter() < 0 && energy > 0 && ncount > 8
                    && posrnd(pos, rad) < (ncount - 8.0)/ncount) {
            die;
         }
      }

      /* Record information that may be used in next global update,
         or by other strands in the next iteration */
      step = pos - posLast;
      mvmt = lerp(mvmt, |step|/rad, 0.5); // should decay to 0
      rad = radius(F(pos));
   }
}

global {
   real totenr = sum { P.energy | P in particle.all};
   real meanmv = mean { P.mvmt | P in particle.all };
   print("=-=-= iter ", iter,
         "; tot#=", numActive(),
         "; totenr=", totenr,
         "; mean(hh)=", mean { P.hh | P in particle.all},
         "; mean(ncount)=", mean { real(P.ncount) | P in particle.all},
         "; mean(mvmt)=", meanmv,
         "\n");
   /* The variation of inter-particle distances seems to decrease the
      reliability of COV(closest) as a convergence indicator, so unlike
      ../sphere this instead uses a particle motion test */
   if (meanmv < eps) {
      print("Stabilizing ", numActive(), " points with mean(movement) ", meanmv,
            " < ", eps, " (iter ", iter, ")\n");
      stabilize;
   }
   iter += 1;
}

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