home *** CD-ROM | disk | FTP | other *** search
/ PC Pro 2002 April / pcpro0402.iso / essentials / graphics / Gimp / gimp-src-20001226.exe / src / gimp / plug-ins / flame / rect.c < prev    next >
Encoding:
C/C++ Source or Header  |  2000-08-25  |  11.1 KB  |  371 lines
  1. /*
  2.    flame - cosmic recursive fractal flames
  3.    Copyright (C) 1992  Scott Draves <spot@cs.cmu.edu>
  4.  
  5.    This program is free software; you can redistribute it and/or modify
  6.    it under the terms of the GNU General Public License as published by
  7.    the Free Software Foundation; either version 2 of the License, or
  8.    (at your option) any later version.
  9.  
  10.    This program is distributed in the hope that it will be useful,
  11.    but WITHOUT ANY WARRANTY; without even the implied warranty of
  12.    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13.    GNU General Public License for more details.
  14.  
  15.    You should have received a copy of the GNU General Public License
  16.    along with this program; if not, write to the Free Software
  17.    Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
  18. */
  19.  
  20. #include "rect.h"
  21.  
  22. #include <string.h>
  23.  
  24.  
  25. /* for batch
  26.  *   interpolate
  27.  *   compute colormap
  28.  *   for subbatch
  29.  *     compute samples
  30.  *     buckets += cmap[samples]
  31.  *   accum += time_filter[batch] * log(buckets)
  32.  * image = filter(accum)
  33.  */
  34.  
  35.  
  36. typedef short bucket[4];
  37.  
  38. /* if you use longs instead of shorts, you 
  39.    get higher quality, and spend more memory */
  40.  
  41. #if 1
  42. typedef short accum_t;
  43. #define MAXBUCKET (1<<14)
  44. #define SUB_BATCH_SIZE 10000
  45. #else
  46. typedef long accum_t;
  47. #define MAXBUCKET (1<<30)
  48. #define SUB_BATCH_SIZE 10000
  49. #endif
  50.  
  51. typedef accum_t abucket[4];
  52.  
  53.  
  54.  
  55. /* allow this many iterations for settling into attractor */
  56. #define FUSE 15
  57.  
  58. /* clamp spatial filter to zero at this std dev (2.5 ~= 0.0125) */
  59. #define FILTER_CUTOFF 2.5
  60.  
  61. /* should be MAXBUCKET / (OVERSAMPLE^2) */
  62. #define PREFILTER_WHITE (MAXBUCKET>>4)
  63.  
  64.  
  65. #define bump_no_overflow(dest, delta, type) { \
  66.    type tt_ = dest + delta;            \
  67.    if (tt_ > dest) dest = tt_;                 \
  68. }
  69.  
  70. /* sum of entries of vector to 1 */
  71. void normalize_vector(v, n)
  72.    double *v;
  73.    int n;
  74. {
  75.    double t = 0.0;
  76.    int i;
  77.    for (i = 0; i < n; i++)
  78.       t += v[i];
  79.    t = 1.0 / t;
  80.    for (i = 0; i < n; i++)
  81.       v[i] *= t;
  82. }
  83.  
  84.  
  85. void render_rectangle(spec, out, out_width, field, nchan, progress)
  86.    frame_spec *spec;
  87.    unsigned char *out;
  88.    int out_width;
  89.    int field;
  90.    int nchan;
  91.    int progress(double);
  92. {
  93.    int i, j, k, nsamples, nbuckets, batch_size, batch_num, sub_batch;
  94.    bucket  *buckets;
  95.    abucket *accumulate;
  96.    point *points;
  97.    double *filter, *temporal_filter, *temporal_deltas;
  98.    double bounds[4], size[2], ppux, ppuy;
  99.    int image_width, image_height;    /* size of the image to produce */
  100.    int width, height;               /* size of histogram */
  101.    int filter_width;
  102.    int oversample = spec->cps[0].spatial_oversample;
  103.    int nbatches = spec->cps[0].nbatches;
  104.    bucket cmap[CMAP_SIZE];
  105.    int gutter_width;
  106.    int sbc;
  107.  
  108.    image_width = spec->cps[0].width;
  109.    if (field) {
  110.       image_height = spec->cps[0].height / 2;
  111.       if (field == field_odd)
  112.      out += nchan * out_width;
  113.       out_width *= 2;
  114.    } else
  115.       image_height = spec->cps[0].height;
  116.  
  117.    if (1) {
  118.       filter_width = (2.0 * FILTER_CUTOFF * oversample *
  119.               spec->cps[0].spatial_filter_radius);
  120.       /* make sure it has same parity as oversample */
  121.       if ((filter_width ^ oversample) & 1)
  122.      filter_width++;
  123.  
  124.       filter = (double *) malloc(sizeof(double) * filter_width * filter_width);
  125.       /* fill in the coefs */
  126.       for (i = 0; i < filter_width; i++)
  127.      for (j = 0; j < filter_width; j++) {
  128.         double ii = ((2.0 * i + 1.0) / filter_width - 1.0) * FILTER_CUTOFF;
  129.         double jj = ((2.0 * j + 1.0) / filter_width - 1.0) * FILTER_CUTOFF;
  130.         if (field)
  131.            jj *= 2.0;
  132.         filter[i + j * filter_width] =
  133.            exp(-2.0 * (ii * ii + jj * jj));
  134.      }
  135.  
  136.       normalize_vector(filter, filter_width * filter_width);
  137.    }
  138.    temporal_filter = (double *) malloc(sizeof(double) * nbatches);
  139.    temporal_deltas = (double *) malloc(sizeof(double) * nbatches);
  140.    if (nbatches > 1) {
  141.       double t;
  142.       /* fill in the coefs */
  143.       for (i = 0; i < nbatches; i++) {
  144.      t = temporal_deltas[i] = (2.0 * ((double) i / (nbatches - 1)) - 1.0)
  145.         * spec->temporal_filter_radius;
  146.      temporal_filter[i] = exp(-2.0 * t * t);
  147.       }
  148.  
  149.       normalize_vector(temporal_filter, nbatches);
  150.    } else {
  151.       temporal_filter[0] = 1.0;
  152.       temporal_deltas[0] = 0.0;
  153.    }
  154.  
  155. #if 0
  156.    for (j = 0; j < nbatches; j++)
  157.       fprintf(stderr, "%4f %4f\n", temporal_deltas[j], temporal_filter[j]);
  158.    fprintf(stderr, "\n");
  159. #endif
  160.    
  161.    /* the number of additional rows of buckets we put at the edge so
  162.       that the filter doesn't go off the edge */
  163.  
  164.    gutter_width = (filter_width - oversample) / 2;
  165.    height = oversample * image_height + 2 * gutter_width;
  166.    width  = oversample * image_width  + 2 * gutter_width;
  167.  
  168.    nbuckets = width * height;
  169.    if (1) {
  170.      static char *last_block = NULL;
  171.      static int last_block_size = 0;
  172.      int memory_rqd = (sizeof(bucket) * nbuckets +
  173.                sizeof(abucket) * nbuckets +
  174.                sizeof(point) * SUB_BATCH_SIZE);
  175.      if (memory_rqd > last_block_size) {
  176.        if (last_block != NULL)
  177.       free(last_block);
  178.        last_block = (char *) malloc(memory_rqd);
  179.        if (NULL == last_block) {
  180.       fprintf(stderr, "render_rectangle: cannot malloc %d bytes.\n", memory_rqd);
  181.       exit(1);
  182.        }
  183.        /* else fprintf(stderr, "render_rectangle: mallocked %dMb.\n", Mb); */
  184.  
  185.        last_block_size = memory_rqd;
  186.      }
  187.      buckets = (bucket *) last_block;
  188.      accumulate = (abucket *) (last_block + sizeof(bucket) * nbuckets);
  189.      points = (point *)  (last_block + (sizeof(bucket) + sizeof(abucket)) * nbuckets);
  190.    }
  191.  
  192.    memset((char *) accumulate, 0, sizeof(abucket) * nbuckets);
  193.    for (batch_num = 0; batch_num < nbatches; batch_num++) {
  194.       double batch_time;
  195.       double sample_density;
  196.       control_point cp;
  197.       memset((char *) buckets, 0, sizeof(bucket) * nbuckets);
  198.       batch_time = spec->time + temporal_deltas[batch_num];
  199.  
  200.       /* interpolate and get a control point */
  201.       interpolate(spec->cps, spec->ncps, batch_time, &cp);
  202.  
  203.       /* compute the colormap entries.  the input colormap is 256 long with
  204.      entries from 0 to 1.0 */
  205.       for (j = 0; j < CMAP_SIZE; j++) {
  206.      for (k = 0; k < 3; k++) {
  207. #if 1
  208.         cmap[j][k] = (int) (cp.cmap[(j * 256) / CMAP_SIZE][k] *
  209.                 cp.white_level);
  210. #else
  211.         /* monochrome if you don't have any cmaps */
  212.         cmap[j][k] = cp.white_level;
  213. #endif        
  214.      }
  215.      cmap[j][3] = cp.white_level;
  216.       }
  217.  
  218.       /* compute camera */
  219.       if (1) {
  220.     double t0, t1, shift = 0.0, corner0, corner1;
  221.     double scale;
  222.  
  223.     scale = pow(2.0, cp.zoom);
  224.     sample_density = cp.sample_density * scale * scale;
  225.     
  226.     ppux = cp.pixels_per_unit * scale;
  227.     ppuy = field ? (ppux / 2.0) : ppux;
  228.     switch (field) {
  229.     case field_both: shift =  0.0; break;
  230.     case field_even: shift = -0.5; break;
  231.     case field_odd:  shift =  0.5; break;
  232.     }
  233.     shift = shift / ppux;
  234.     t0 = (double) gutter_width / (oversample * ppux);
  235.     t1 = (double) gutter_width / (oversample * ppuy);
  236.     corner0 = cp.center[0] - image_width / ppux / 2.0;
  237.     corner1 = cp.center[1] - image_height / ppuy / 2.0;
  238.     bounds[0] = corner0 - t0;
  239.     bounds[1] = corner1 - t1 + shift;
  240.     bounds[2] = corner0 + image_width  / ppux + t0;
  241.     bounds[3] = corner1 + image_height / ppuy + t1 + shift;
  242.     size[0] = 1.0 / (bounds[2] - bounds[0]);
  243.     size[1] = 1.0 / (bounds[3] - bounds[1]);
  244.       }
  245.  
  246.       nsamples = (int) (sample_density * nbuckets /
  247.             (oversample * oversample));
  248.    
  249.       batch_size = nsamples / cp.nbatches;
  250.  
  251.       sbc = 0;
  252.       for (sub_batch = 0;
  253.        sub_batch < batch_size;
  254.        sub_batch += SUB_BATCH_SIZE) {
  255.  
  256.     if (progress&&!(sbc++&7))
  257.       (*progress)(0.5*sub_batch/(double)batch_size);
  258.  
  259.      /* generate a sub_batch_size worth of samples */
  260.      points[0][0] = random_uniform11();
  261.      points[0][1] = random_uniform11();
  262.      points[0][2] = random_uniform01();
  263.      iterate(&cp, SUB_BATCH_SIZE, FUSE, points);
  264.      
  265.      
  266.      /* merge them into buckets, looking up colors */
  267.      for (j = 0; j < SUB_BATCH_SIZE; j++) {
  268.         int k, color_index;
  269.         double *p = points[j];
  270.         bucket *b;
  271.         if (p[0] < bounds[0] ||
  272.         p[1] < bounds[1] ||
  273.         p[0] > bounds[2] ||
  274.         p[1] > bounds[3])
  275.            continue;
  276.         color_index = (int) (p[2] * CMAP_SIZE);
  277.         if (color_index < 0) color_index = 0;
  278.         else if (color_index > (CMAP_SIZE-1))
  279.            color_index = CMAP_SIZE-1;
  280.         b = buckets +
  281.            (int) (width * (p[0] - bounds[0]) * size[0]) +
  282.           width * (int) (height * (p[1] - bounds[1]) * size[1]);        
  283.         for (k = 0; k < 4; k++)
  284.            bump_no_overflow(b[0][k], cmap[color_index][k], short);
  285.      }
  286.       }
  287.       
  288.       if (1) {
  289.      double k1 =(cp.contrast * cp.brightness *
  290.              PREFILTER_WHITE * 268.0 *
  291.              temporal_filter[batch_num]) / 256;
  292.      double area = image_width * image_height / (ppux * ppuy);
  293.      double k2 = (oversample * oversample * nbatches) /
  294.         (cp.contrast * area * cp.white_level * sample_density);
  295.  
  296.      /* log intensity in hsv space */
  297.      for (j = 0; j < height; j++)
  298.         for (i = 0; i < width; i++) {
  299.            abucket *a = accumulate + i + j * width;
  300.            bucket *b = buckets + i + j * width;
  301.            double c[4], ls;
  302.            c[0] = (double) b[0][0];
  303.            c[1] = (double) b[0][1];
  304.            c[2] = (double) b[0][2];
  305.            c[3] = (double) b[0][3];
  306.            if (0.0 == c[3])
  307.          continue;
  308.            
  309.            ls = (k1 * log(1.0 + c[3] * k2))/c[3];
  310.            c[0] *= ls;
  311.            c[1] *= ls;
  312.            c[2] *= ls;
  313.            c[3] *= ls;
  314.  
  315.            bump_no_overflow(a[0][0], c[0] + 0.5, accum_t);
  316.            bump_no_overflow(a[0][1], c[1] + 0.5, accum_t);
  317.            bump_no_overflow(a[0][2], c[2] + 0.5, accum_t);
  318.            bump_no_overflow(a[0][3], c[3] + 0.5, accum_t);
  319.         }
  320.       }
  321.    }
  322.    /*
  323.     * filter the accumulation buffer down into the image
  324.     */
  325.    if (1) {
  326.       int x, y;
  327.       double t[4];
  328.       double g = 1.0 / spec->cps[0].gamma;
  329.       y = 0;
  330.       for (j = 0; j < image_height; j++) {
  331.      if (progress && !(j&15))
  332.        (*progress)(0.5+0.5*j/(double)image_height);
  333.      x = 0;
  334.      for (i = 0; i < image_width; i++) {
  335.         int ii, jj, a;
  336.         unsigned char *p;
  337.         t[0] = t[1] = t[2] = t[3] = 0.0;
  338.         for (ii = 0; ii < filter_width; ii++)
  339.            for (jj = 0; jj < filter_width; jj++) {
  340.           double k = filter[ii + jj * filter_width];
  341.           abucket *a = accumulate + x + ii + (y + jj) * width;
  342.           
  343.           t[0] += k * a[0][0];
  344.           t[1] += k * a[0][1];
  345.           t[2] += k * a[0][2];
  346.           t[3] += k * a[0][3];
  347.            }
  348.         p = out + nchan * (i + j * out_width);
  349.         a = 256.0 * pow((double) t[0] / PREFILTER_WHITE, g) + 0.5;
  350.         if (a < 0) a = 0; else if (a > 255) a = 255;
  351.         p[0] = a;
  352.         a = 256.0 * pow((double) t[1] / PREFILTER_WHITE, g) + 0.5;
  353.         if (a < 0) a = 0; else if (a > 255) a = 255;
  354.         p[1] = a;
  355.         a = 256.0 * pow((double) t[2] / PREFILTER_WHITE, g) + 0.5;
  356.         if (a < 0) a = 0; else if (a > 255) a = 255;
  357.         p[2] = a;
  358.         if (nchan > 3) {
  359.           a = 256.0 * pow((double) t[3] / PREFILTER_WHITE, g) + 0.5;
  360.           if (a < 0) a = 0; else if (a > 255) a = 255;
  361.           p[3] = a;
  362.         }
  363.         x += oversample;
  364.      }
  365.      y += oversample;
  366.       }
  367.    }
  368.    
  369.    free(filter);
  370. }
  371.