Flower disk sampling for the thin lens

on 8 Sep 2011 by Mukund (@muks)

A few days ago, I was on my morning walk when I saw a bunch of flowers with the distinctive spiral florets pattern. The day before yesterday, it occured that this might be a good way to make jittered samples off a thin lens for effects like depth of field. Currently pretty much every open source raytracer seems to use the concentric map as described by Shirley. The problem is essentially one of space-filling a disk with many smaller equal-sized disks. A formula for the flower pattern was described by H. Vogel in A better way to construct the sunflower head. A sample program implementing this follows. I had a bit of a struggle automatically adjusting the size of florets depending on their number, but that's now done. The program is structured such that lengths are normalized and it can easily be incorporated inside a raytracer. It can generate random jittered samples inside florets. Samples can further be weighed as necessary. The program can be adapted to vary the size of the florets for importance sampling instead of weighing.

360 florets

Flower disk sampling (florets)

360 florets with jittered samples (in black)

Flower disk sampling

1000 florets with jittered samples (in black)

This is just a demonstration of the space filling, and won't be useful I suppose.

Flower disk sampling

Just the 1000 jittered samples (in black)

Flower disk sampling

10000 florets

Flower disk sampling

10000 florets with jittered samples (in black)

Flower disk sampling

Just the 10000 jittered samples (in black)

Flower disk sampling

Guess

Flower disk sampling

16 florets with jittered samples (in black)

Flower disk sampling

Just the 16 jittered samples (in black)

Flower disk sampling

Here's the code. It uses Cairo.

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <unistd.h>
#include <cairo/cairo.h>

#define SIZE         800
#define RADIUS_SCALE (SIZE / 2.1)
#define ORIGIN       (SIZE / 2)

static int show_florets = 1;
static int show_samples = 0;
static int inside = 0;
static int verbose = 0;
static int size = SIZE;
static float radius_scale = RADIUS_SCALE;
static float origin = ORIGIN;

static inline float
frandom (void)
{
  return random() / (RAND_MAX * 1.0f);
}

static inline float
golden_angle (void)
{
  return M_PI * (3.0 - sqrtf (5.0));
}

static void
draw_floret (cairo_t *c,
             int      i,
             float    t,
             float    r,
             float    s)
{
  float x, y;

  x = r * cosf (t);
  y = r * sinf (t);

  if (verbose)
    printf ("Drawing floret %d at theta = %f, radius = %f, "
            "floret size (radius) = %f, x = %f y = %f\n",
            i, t, r, s, x, y);

  if (show_florets) {
    /* Draw the floret in a random color */
    cairo_new_sub_path (c);
    cairo_arc (c,
               origin + (x * radius_scale),
               origin + (y * radius_scale),
               s * radius_scale,
               0.0, 2 * M_PI);
    cairo_set_source_rgba (c, frandom (), frandom (), frandom (), 0.5);
    cairo_fill (c);
  }

  if (show_samples) {
    float rx, ry;
    float rs;
    float phi;
    float l;

    /* Draw a tiny black dot for the random sample inside the floret. */

    rs = s / 8;

    /* Make sure the random samples are inside the flower disk (put them
       in the implicit eqn of a disk) */
    rx = ry = -2.0;
    while (((rx * rx) + (ry * ry)) >= 1.0) {
      phi = frandom() * 2 * M_PI;
      l = (s - rs) * frandom ();

      rx = x + l * cosf (phi);
      ry = y + l * sinf (phi);
    }

    cairo_new_sub_path (c);
    cairo_arc (c,
               origin + (rx * radius_scale),
               origin + (ry * radius_scale),
               rs * radius_scale,
               0.0, 2 * M_PI);
    cairo_set_source_rgb (c, 0.0, 0.0, 0.0);
    cairo_fill (c);
  }
}

static void
draw_flower (cairo_t *c,
             int      n)
{
  float t, ga;
  float s;
  float rm;
  int i;

  t = 2 * M_PI * frandom();
  ga = golden_angle ();
  s = 1 / sqrtf (sqrtf (2.0) * n);
  rm = 1 / sqrtf (n);

  if (inside) {
    /* Scale radius_scale such that the entire floret, not just its
       center, is inside the flower disk. */
    radius_scale *= 1 - s;
  }

  for (i = 0; i < n; i++) {
    float r;

    r = rm * sqrtf (i);

    draw_floret (c, i, t, r, s);

    t = fmodf (t + ga, 2 * M_PI);
  }
}

static void
show_usage (const char *program)
{
  fprintf (stderr,
           "Usage: %s [-n num_florets] [-w width] [-f] [-i] [-s] [-v] "
           "filename.png\n\n",
           program);
  fprintf (stderr,
           " -n: Number of florets (and hence samples) to draw\n"
           " -w: Width and height of output image\n"
           " -f: Don't show florets\n"
           " -i: Don't draw florets outside the flower disk "
                 "(default keeps centers inside)\n"
           " -s: Show the jittered sample points (in black)\n"
           " -v: Show messages\n"
           "\n"
           "Florets or samples or both must be enabled for draw.\n"
           "num_florets should be > 0.\n"
           "width should be > 16.\n"
           );

  exit (1);
}

int
main (int argc, char *argv[])
{
  int opt;
  const char *filename = "flower.png";
  int num_florets = 1000;
  cairo_surface_t *s;
  cairo_t *c;

  while ((opt = getopt (argc, argv, "fin:svw:")) != -1) {
    switch (opt) {
    case 'f':
      show_florets = 0;
      break;
    case 'i':
      inside = 1;
      break;
    case 's':
      show_samples = 1;
      break;
    case 'n':
      num_florets = atoi (optarg);
      break;
    case 'v':
      verbose = 1;
      break;
    case 'w':
      size = atoi (optarg);
      radius_scale = size / 2.1;
      origin = size / 2.0;
      break;
    default:
      show_usage (argv[0]);
    }
  }

  if ((optind >= argc) ||
      (num_florets < 1) ||
      (size < 16) ||
      (!show_florets && !show_samples))
    show_usage (argv[0]);

  filename = argv[optind];

  srandom (time (NULL));

  s = cairo_image_surface_create (CAIRO_FORMAT_RGB24,
                                  size, size);
  c = cairo_create (s);

  /* Fill background with white */
  cairo_set_source_rgb (c, 1.0, 1.0, 1.0);
  cairo_paint (c);

  /* Draw the lens disk outline */
  cairo_new_sub_path (c);
  cairo_arc (c, origin, origin, radius_scale, 0.0, 2 * M_PI);
  cairo_set_source_rgb (c, 0.0, 0.0, 0.0);
  cairo_stroke (c);

  /* Draw the florets */
  draw_flower (c, num_florets);

  cairo_destroy (c);

  cairo_surface_write_to_png (s, filename);
  cairo_surface_destroy (s);

  return 0;
}

Download flower.tar.gz for your making pleasure.

Update: Bugs in sample generation, and in the usage message were fixed.