Draw random black-and-white forest Your task is to write program which will draw 800x600 black-and-white image with something resembling a forest.

Like this (it is dithered photo):



Rules


You are disallowed to use any existing images - you should generate image purely algorithmically
Use only 2 colors - black and white (no grayscale)
Each time program runs image should be new - random every time
One tree is not a forest (let say 5 trees minumum)
Special libraries for drawing trees/forests are disallowed
Answer with most votes wins
C: 3863 1144 1023 999 942 927

The original solution saves 2 pnm files per run (one with g appended, before dithering). Because the dithering wasn't beautiful for the first few lines, there is a hack in place to render more lines than needed, and crop during output.

The golfed solution has a simpler dithering and saves only the dithered image. 
(no warnings with gcc -std=c11 -pedantic -Wall -Wextra)

Example images from 3 original program runs and one run of the golfed version (last image):






Golfed version

  #include &lt;math.h&gt;
  #include &lt;time.h&gt;
  #include &lt;stdlib.h&gt;
  #include &lt;stdio.h&gt;
  #define D float
  #define R rand()/RAND_MAX
  #define Z(a,b,c) if(10.*R&gt;a)T(h,s,j+b+c*R,g);
  #define U a[y][x]
  #define V e[y+1][x
  #define F(x) for(i=0;i&lt;x;++i)
  int i,x,y,W=800,H=600;unsigned char a[600][800],c;D e[601][802],r,b,k,l,m,n,f,w,
  d,q=.01;void T(D h,D s,D j,D g){r=b=0;do{j+=.04*R-.02;h+=sin(j)*q;s+=cos(j)*q;b
  +=q;g+=q;f=525/d;r=.25-g/44;m=w*f+s*f+W/2;n=2*f-h*f+H/2;f*=r;for(y=n-f-2;y&lt;n+f+2
  ;++y)if(y&gt;=0&amp;&amp;y&lt;H)for(x=m-f-2;x&lt;m+f+2;++x)if(x&gt;=0&amp;&amp;x&lt;W)if(k=m-x,l=n-y,f&gt;sqrt(k*k
  +l*l))if(U&gt;d*3)U=d*3;}while(b&lt;10*r+12*r*R&amp;&amp;r&gt;q);if(r&gt;q){Z(2,.26,.35)Z(2,-.26,-
  .35)Z(7,-.05,.1)}}int main(){FILE* o=fopen("i","wb");srand(time(0));F(W*H){y=i/W
  ;a[y][i%W]=(y&lt;313)?255:6e3/(2*y-H);}F(200)w=1e2*R-60,d=80.*R+5,T(0,0,1.58,0);F(W
  *H){x=i%W;y=i/W;k=U+e[y][x+1];U=-(k&gt;0);l=(k-U)*.1;e[y][x+2]+=l*4;V]+=l*2;V+1]+=l
  *3;V+2]+=l;}fprintf(o,"P5 800 600 255 ");fwrite(a,1,W*H,o);}


Original version

  #include &lt;math.h&gt;
  #include &lt;stdio.h&gt;
  #include &lt;stdlib.h&gt;

  #define W 800
  #define H 600
  #define SPEED 0.01
  #define HEIGHT 11.0

  #define R(m) ((double)(m) * rand() / RAND_MAX)
  #define RAD(deg) ((deg) / 180.0 * M_PI)
  #define LIMIT(x, min, max) ((x) &lt; (min) ? (min) : (x) &gt; (max) ? (max) : (x))

  void shade(void);
  void growTree(double dist, double side, double h, double s, double alpha, double grown);
  void plot(double dist, double side, double h, double s, double alpha, double diam);
  void dither(void);
  void writeImg(int dither);

  unsigned char img[H+10][W];
  double err[H+10+2][W+4];
  long tim;

  int main(void)
  {
     int i;
     tim = time(0);
     srand(tim);
     shade();
     for(i = 0; i &lt; 200; ++i)
     {
        growTree(5 + R(75), -60 + R(120), 0.0, 0.0, RAD(90), 0.0);
     }
     writeImg(0);
     dither();
     writeImg(1);
  }

  void shade(void)
  {
     int y;
     for(y = -10; y &lt; H; ++y)
     {
        double dist = H * 3.5 / (2 * y - H);
        unsigned char color = dist / 80 * 255;
        if(y &lt;= H / 2 || dist &gt; 80) color = 255;
        memset(img[y+10], color, W);
     }
  }

  void growTree(double dist, double side, double h, double s, double alpha, double grown)
  {
     double diam, branchLength = 0.0;

     do
     {
        alpha += R(RAD(3)) - RAD(1.5);
        h += sin(alpha) * SPEED;
        s += cos(alpha) * SPEED;
        branchLength += SPEED;
        grown += SPEED;
        diam = (1.0 - grown / HEIGHT) * 0.5;
        plot(dist, side, h, s, alpha, diam);
     } while(branchLength &lt; 5 * diam + R(6 * diam) &amp;&amp; diam &gt; 0.02);

     if(diam &gt; 0.02)
     {
        int br = 0;

        if(R(10) &gt; 2) br++,growTree(dist, side, h, s, alpha + RAD(15) + R(RAD(20)), grown);
        if(R(10) &gt; 2) br++,growTree(dist, side, h, s, alpha - RAD(15) - R(RAD(20)), grown);
        if(R(10) &lt; 2 || br == 0) growTree(dist, side, h, s, alpha - RAD(2.5) + R(RAD(5)), grown);
     }
  }

  void plot(double dist, double side, double h, double s, double alpha, double diam)
  {
     int x, y;
     double scale = H / 4.0 * 3.5 / dist;
     double x0 = side * scale + s * scale + W / 2.0;
     double y0 = H / 2.0 + 2.0 * scale - h * scale;
     diam *= scale;
     h *= scale;
     s *= scale;
     for(y = y0 - diam / 2 - 2; y &lt; y0 + diam / 2 + 2; ++y)
     {
        if(y &lt; -10 || y &gt;= H) continue;
        for(x = x0 - diam / 2 - 2; x &lt; x0 + diam / 2 + 2; ++x)
        {
           double dx, dy, d;
           if(x &lt; 0 || x &gt;= W) continue;
           dx = x0 - x;
           dy = y0 - y;
           d = diam / 2 - sqrt(dx * dx + dy * dy) + 0.5;
           if(d &gt; 0)
           {
              unsigned char color = dist / 80 * 255;
              if(img[y+10][x] &gt; color) img[y+10][x] = color;
           }
        }
     }
  }

  void dither(void)
  {
     int x0, x, y;
     for(y = -10; y &lt; H; ++y)
     {
        for(x0 = 0; x0 &lt; W; ++x0)
        {
           double error, oldpixel;
           unsigned char newpixel;
           if(y%2) x = W - 1 - x0;
           else x = x0;
           oldpixel = img[y+10][x] + err[y+10][x+2];
           newpixel = oldpixel &gt; 127 ? 255 : 0;
           img[y+10][x] = newpixel;
           error = oldpixel - newpixel;
           err[y+10  ][x+1+2*(1-y%2)] += error * 7 / 48;
           err[y+10  ][x+4*(1-y%2)] += error * 5 / 48;
           err[y+10+1][x  ] += error * 3 / 48;
           err[y+10+1][x+1] += error * 5 / 48;
           err[y+10+1][x+2] += error * 7 / 48;
           err[y+10+1][x+3] += error * 5 / 48;
           err[y+10+1][x+4] += error * 3 / 48;
           err[y+10+2][x  ] += error * 1 / 48;
           err[y+10+2][x+1] += error * 3 / 48;
           err[y+10+2][x+2] += error * 5 / 48;
           err[y+10+2][x+3] += error * 3 / 48;
           err[y+10+2][x+4] += error * 1 / 48;
        }
     }
  }

  void writeImg(int dither)
  {
     FILE* fp;
     char buffer[32];
     sprintf(buffer, "%ld%s.pnm", tim, dither ? "" : "g");
     fp = fopen(buffer, "wb");
     fprintf(fp, "P5\n%d %d\n255\n", W, H);
     fwrite(&amp;img[10][0], 1, W * H, fp);
     fclose(fp);
  }