home *** CD-ROM | disk | FTP | other *** search
/ vsiftp.vmssoftware.com / VSIPUBLIC@vsiftp.vmssoftware.com.tar / FREEWARE / FREEWARE40.ZIP / xpaint-247 / texture.c < prev    next >
C/C++ Source or Header  |  1997-01-03  |  16KB  |  598 lines
  1. /*
  2.  * This file is part of XPaint.
  3.  *
  4.  * Copyright 1996 by Torsten Martinsen
  5.  * Copyright 1994 by Jer Johnson
  6.  *
  7.  * This is free software; you can redistribute it and/or modify
  8.  * it under the terms of the GNU General Public License as published by
  9.  * the Free Software Foundation; either version 2 of the License, or
  10.  * (at your option) any later version.
  11.  */
  12.  
  13. #include <X11/X.h>
  14. #include <X11/Intrinsic.h>
  15. #include <X11/Xlib.h>
  16.  
  17. #include <stdio.h>
  18. #include <math.h>
  19. #include <stdlib.h>
  20.  
  21. #include "xpaint.h"
  22. #include "misc.h"
  23. #include "image.h"
  24.  
  25. #ifndef M_PI
  26. #define M_PI            3.14159265358979323846
  27. #endif
  28.  
  29. #if FEATURE_FRACTAL
  30.  
  31. static void draw_grid(int x, int y, unsigned char col, unsigned char *grid);
  32. static void subdiv(int x1, int y1, int x2, int y2, unsigned char *grid);
  33. static void adjust(int x1, int y1, int x, int y, int x2, int y2,
  34.            unsigned char *grid);
  35. static void fourn(float data[], int nn[], int ndim, int isign);
  36. static void spectralsynth(float **x, unsigned int n, double h);
  37.  
  38. static int number, numcolors, width, height;
  39.  
  40. /*
  41.  * draw_plasma() and friends are stolen from Xtacy, written by Jer Johnson
  42.  * (mpython@gnu.ai.mit.edu) who says that his work in turn is based on OOZE.C,
  43.  * written by Jeff Clough.
  44.  */
  45. Image *
  46. draw_plasma(int w, int h)
  47. {
  48.     int i, n;
  49.     unsigned char *grid;
  50.     Image *output;
  51.     unsigned char *op;
  52.  
  53.  
  54.     /* customizable parameters */
  55.     numcolors = 100;
  56.     number = 10;
  57.  
  58.     width = w;
  59.  
  60.     n = numcolors / 5;
  61.     numcolors = n * 5;
  62.     output = ImageNewCmap(w, h, numcolors);
  63.     grid = output->data;
  64.     memset(grid, 0, w * h);
  65.  
  66.     /*
  67.      * Make points at (0,0) (0, h-1) (w-1, 0) (w-1, h-1) of random colors
  68.      */
  69.     draw_grid(0, 0, RANDOMI2(0, numcolors-1), grid);
  70.     draw_grid(w - 1, 0, RANDOMI2(0, numcolors-1), grid);
  71.     draw_grid(0, h - 1, RANDOMI2(0, numcolors-1), grid);
  72.     draw_grid(w - 1, h - 1, RANDOMI2(0, numcolors-1), grid);
  73.  
  74.     subdiv(0, 0, w - 1, h - 1, grid);
  75.  
  76.     /* Build the colour map */
  77.     op = output->cmapData;
  78.     for (i = 0; i < n; ++i) {
  79.     *op++ = 255;
  80.     *op++ = 255 * i / n;
  81.     *op++ = 0;
  82.     }
  83.     for (i = 0; i < n; ++i) {
  84.     *op++ = 255 - 255 * i / n;
  85.     *op++ = 255;
  86.     *op++ = 0;
  87.     }
  88.     for (i = 0; i < n; ++i) {
  89.     *op++ = 0;
  90.     *op++ = 255;
  91.     *op++ = 255 * i / n;
  92.     }
  93.     for (i = 0; i < n; ++i) {
  94.     *op++ = 255 * i / n;
  95.     *op++ = 0;
  96.     *op++ = 255;
  97.     }
  98.     for (i = 0; i < n; ++i) {
  99.     *op++ = 255;
  100.     *op++ = 0;
  101.     *op++ = 255 - 255 * i / n;
  102.     }
  103. #if 0
  104.     op = output->cmapData;
  105.     for (i = 0; i < 5 * n; ++i) {
  106.     printf("%d: %d %d %d\n", i, op[0], op[1], op[2]);
  107.     op += 3;
  108.     }
  109. #endif
  110.  
  111.     return output;
  112. }
  113.  
  114. static void 
  115. draw_grid(int x, int y, unsigned char col, unsigned char *grid)
  116. {
  117.     grid[x + y * width] = col;
  118. }
  119.  
  120. static void 
  121. subdiv(int x1, int y1, int x2, int y2, unsigned char *grid)
  122. {
  123.     int x, y;
  124.  
  125.     if ((x2 - x1 < 2) && (y2 - y1 < 2))
  126.     return;
  127.  
  128.     x = (x1 + x2) / 2;
  129.     y = (y1 + y2) / 2;
  130.     adjust(x1, y1, x, y1, x2, y1, grid);
  131.     adjust(x2, y1, x2, y, x2, y2, grid);
  132.     adjust(x1, y2, x, y2, x2, y2, grid);
  133.     adjust(x1, y1, x1, y, x1, y2, grid);
  134.  
  135.     if (grid[x + y * width] == 0) {
  136.     unsigned char spooge;
  137.  
  138.     spooge = ((int) grid[x1 + y1 * width] +
  139.           (int) grid[x1 + y2 * width] +
  140.           (int) grid[x2 + y1 * width] +
  141.           (int) grid[x2 + y2 * width]) / 4;
  142.     draw_grid(x, y, spooge, grid);
  143.     }
  144.     subdiv(x1, y1, x, y, grid);
  145.     subdiv(x, y1, x2, y, grid);
  146.     subdiv(x1, y, x, y2, grid);
  147.     subdiv(x, y, x2, y2, grid);
  148. }
  149.  
  150. static void 
  151. adjust(int x1, int y1, int x, int y, int x2, int y2, unsigned char *grid)
  152. {
  153.     int c, d, spoo;
  154.     int clr1, clr2, clr3, clr4, horz, vert;
  155.  
  156.  
  157.     if (grid[x + y * width] != 0)
  158.     return;
  159.  
  160.     clr1 = grid[x1 + y1 * width];
  161.     clr2 = grid[x2 + y2 * width];
  162.     clr3 = grid[x1 + y2 * width];
  163.     clr4 = grid[x2 + y1 * width];
  164.     horz = x2 - x1;
  165.     vert = y2 - y1;
  166.     d = hypot(horz, vert);
  167.  
  168. #if 1                /* 0 for band effect */
  169.     if ((spoo = RANDOMI2(0, 100)) > number) {
  170.     if (spoo % 2 == 1)
  171.         c = (((clr1 + clr2) / 2) + ((int) RANDOMI2(0, d))) % numcolors;
  172.     else
  173.         c = (((clr1 + clr2) / 2) - ((int) RANDOMI2(0, d))) % numcolors;
  174.     } else
  175. #endif
  176.     c = (clr1 + clr2 + clr3 + clr4) / 4;
  177.     if (c < 0)
  178.     c = -c;
  179.     else if (c == 0)
  180.     c = 1;
  181.     draw_grid(x, y, c, grid);
  182. }
  183.  
  184. /*
  185.  * Fractal landscape generator from ppmforge.c,
  186.  * written by John Walker (kelvin@Autodesk.com).
  187.  */
  188.  
  189. #define Real(v, x, y)  v[1 + (((x) * meshsize) + (y)) * 2]
  190. #define Imag(v, x, y)  v[2 + (((x) * meshsize) + (y)) * 2]
  191.  
  192. #define ELEMENTS(array) (sizeof(array)/sizeof((array)[0]))
  193.  
  194. /* customize: */
  195. static double fracdim = 1.8;    /* Fractal dimension - ]0,3] */
  196. static double powscale = 1.0;    /* Power law scaling exponent - ]0,3] */
  197. static int meshsize = 256;    /* FFT mesh size - generally >= max(w,h),
  198.                    but small numbers can give special effects */
  199. /* colours */
  200.  
  201. static unsigned char pgnd[][3] =
  202. {
  203.     {206, 205, 0},
  204.     {208, 207, 0},
  205.     {211, 208, 0},
  206.     {214, 208, 0},
  207.     {217, 208, 0},
  208.     {220, 208, 0},
  209.     {222, 207, 0},
  210.     {225, 205, 0},
  211.     {227, 204, 0},
  212.     {229, 202, 0},
  213.     {231, 199, 0},
  214.     {232, 197, 0},
  215.     {233, 194, 0},
  216.     {234, 191, 0},
  217.     {234, 188, 0},
  218.     {233, 185, 0},
  219.     {232, 183, 0},
  220.     {231, 180, 0},
  221.     {229, 178, 0},
  222.     {227, 176, 0},
  223.     {225, 174, 0},
  224.     {223, 172, 0},
  225.     {221, 170, 0},
  226.     {219, 168, 0},
  227.     {216, 166, 0},
  228.     {214, 164, 0},
  229.     {212, 162, 0},
  230.     {210, 161, 0},
  231.     {207, 159, 0},
  232.     {205, 157, 0},
  233.     {203, 156, 0},
  234.     {200, 154, 0},
  235.     {198, 152, 0},
  236.     {195, 151, 0},
  237.     {193, 149, 0},
  238.     {190, 148, 0},
  239.     {188, 147, 0},
  240.     {185, 145, 0},
  241.     {183, 144, 0},
  242.     {180, 143, 0},
  243.     {177, 141, 0},
  244.     {175, 140, 0},
  245.     {172, 139, 0},
  246.     {169, 138, 0},
  247.     {167, 137, 0},
  248.     {164, 136, 0},
  249.     {161, 135, 0},
  250.     {158, 134, 0},
  251.     {156, 133, 0},
  252.     {153, 132, 0},
  253.     {150, 132, 0},
  254.     {147, 131, 0},
  255.     {145, 130, 0},
  256.     {142, 130, 0},
  257.     {139, 129, 0},
  258.     {136, 128, 0},
  259.     {133, 128, 0},
  260.     {130, 127, 0},
  261.     {127, 127, 0},
  262.     {125, 127, 0},
  263.     {122, 127, 0},
  264.     {119, 127, 0},
  265.     {116, 127, 0},
  266.     {113, 127, 0},
  267.     {110, 128, 0},
  268.     {107, 128, 0},
  269.     {104, 128, 0},
  270.     {102, 127, 0},
  271.     {99, 126, 0},
  272.     {97, 124, 0},
  273.     {95, 122, 0},
  274.     {93, 120, 0},
  275.     {92, 117, 0},
  276.     {92, 114, 0},
  277.     {92, 111, 0},
  278.     {93, 108, 0},
  279.     {94, 106, 0},
  280.     {96, 104, 0},
  281.     {98, 102, 0},
  282.     {100, 100, 0},
  283.     {103, 99, 0},
  284.     {106, 99, 0},
  285.     {109, 99, 0},
  286.     {111, 100, 0},
  287.     {114, 101, 0},
  288.     {117, 102, 0},
  289.     {120, 103, 0},
  290.     {123, 102, 0},
  291.     {125, 102, 0},
  292.     {128, 100, 0},
  293.     {130, 98, 0},
  294.     {132, 96, 0},
  295.     {133, 94, 0},
  296.     {134, 91, 0},
  297.     {134, 88, 0},
  298.     {134, 85, 0},
  299.     {133, 82, 0},
  300.     {131, 80, 0},
  301.     {129, 78, 0}
  302. };
  303.  
  304.  
  305. /*      FOURN  --  Multi-dimensional fast Fourier transform
  306.  
  307.    Called with arguments:
  308.  
  309.    data       A  one-dimensional  array  of  floats  (NOTE!!!   NOT
  310.    DOUBLES!!), indexed from one (NOTE!!!   NOT  ZERO!!),
  311.    containing  pairs of numbers representing the complex
  312.    valued samples.  The Fourier transformed results     are
  313.    returned in the same array.
  314.  
  315.    nn         An  array specifying the edge size in each dimension.
  316.    THIS ARRAY IS INDEXED FROM  ONE,     AND  ALL  THE  EDGE
  317.    SIZES MUST BE POWERS OF TWO!!!
  318.  
  319.    ndim       Number of dimensions of FFT to perform.  Set to 2 for
  320.    two dimensional FFT.
  321.  
  322.    isign      If 1, a Fourier transform is done; if -1 the  inverse
  323.    transformation is performed.
  324.  
  325.    This  function  is essentially as given in Press et al., "Numerical
  326.    Recipes In C", Section 12.11, pp.  467-470.
  327.  */
  328.  
  329. static void 
  330. fourn(float data[], int nn[], int ndim, int isign)
  331. {
  332.     register int i1, i2, i3;
  333.     int i2rev, i3rev, ip1, ip2, ip3, ifp1, ifp2;
  334.     int ibit, idim, k1, k2, n, nprev, nrem, ntot;
  335.     float tempi, tempr;
  336.     double theta, wi, wpi, wpr, wr, wtemp;
  337.  
  338. #define SWAP(a,b) tempr=(a); (a) = (b); (b) = tempr
  339.  
  340.     ntot = 1;
  341.     for (idim = 1; idim <= ndim; idim++)
  342.     ntot *= nn[idim];
  343.     nprev = 1;
  344.     for (idim = ndim; idim >= 1; idim--) {
  345.     n = nn[idim];
  346.     nrem = ntot / (n * nprev);
  347.     ip1 = nprev << 1;
  348.     ip2 = ip1 * n;
  349.     ip3 = ip2 * nrem;
  350.     i2rev = 1;
  351.     for (i2 = 1; i2 <= ip2; i2 += ip1) {
  352.         if (i2 < i2rev) {
  353.         for (i1 = i2; i1 <= i2 + ip1 - 2; i1 += 2) {
  354.             for (i3 = i1; i3 <= ip3; i3 += ip2) {
  355.             i3rev = i2rev + i3 - i2;
  356.             SWAP(data[i3], data[i3rev]);
  357.             SWAP(data[i3 + 1], data[i3rev + 1]);
  358.             }
  359.         }
  360.         }
  361.         ibit = ip2 >> 1;
  362.         while (ibit >= ip1 && i2rev > ibit) {
  363.         i2rev -= ibit;
  364.         ibit >>= 1;
  365.         }
  366.         i2rev += ibit;
  367.     }
  368.     ifp1 = ip1;
  369.     while (ifp1 < ip2) {
  370.         ifp2 = ifp1 << 1;
  371.         theta = isign * (M_PI * 2) / (ifp2 / ip1);
  372.         wtemp = sin(0.5 * theta);
  373.         wpr = -2.0 * wtemp * wtemp;
  374.         wpi = sin(theta);
  375.         wr = 1.0;
  376.         wi = 0.0;
  377.         for (i3 = 1; i3 <= ifp1; i3 += ip1) {
  378.         for (i1 = i3; i1 <= i3 + ip1 - 2; i1 += 2) {
  379.             for (i2 = i1; i2 <= ip3; i2 += ifp2) {
  380.             k1 = i2;
  381.             k2 = k1 + ifp1;
  382.             tempr = wr * data[k2] - wi * data[k2 + 1];
  383.             tempi = wr * data[k2 + 1] + wi * data[k2];
  384.             data[k2] = data[k1] - tempr;
  385.             data[k2 + 1] = data[k1 + 1] - tempi;
  386.             data[k1] += tempr;
  387.             data[k1 + 1] += tempi;
  388.             }
  389.         }
  390.         wr = (wtemp = wr) * wpr - wi * wpi + wr;
  391.         wi = wi * wpr + wtemp * wpi + wi;
  392.         }
  393.         ifp1 = ifp2;
  394.     }
  395.     nprev *= n;
  396.     }
  397. }
  398. #undef SWAP
  399.  
  400. /*
  401.  * SPECTRALSYNTH -- Spectrally synthesised fractal motion in two dimensions.
  402.  */
  403. static void 
  404. spectralsynth(float **x, unsigned int n, double h)
  405. {
  406.     unsigned bl;
  407.     int i, j, i0, j0, nsize[3];
  408.     double rad, phase, rcos, rsin;
  409.     float *a;
  410.  
  411.     bl = ((((unsigned long) n) * n) + 1) * 2 * sizeof(float);
  412.     a = (float *) XtCalloc(bl, 1);
  413.     *x = a;
  414.  
  415.     for (i = 0; i <= n / 2; i++) {
  416.     for (j = 0; j <= n / 2; j++) {
  417.         phase = 2 * M_PI * (RANDOMI() & 0x7FFF) / 32767.0;
  418.         if (i != 0 || j != 0) {
  419.         rad = pow((double) (i * i + j * j), -(h + 1) / 2) * gauss();
  420.         } else {
  421.         rad = 0;
  422.         }
  423.         rcos = rad * cos(phase);
  424.         rsin = rad * sin(phase);
  425.         Real(a, i, j) = rcos;
  426.         Imag(a, i, j) = rsin;
  427.         i0 = (i == 0) ? 0 : n - i;
  428.         j0 = (j == 0) ? 0 : n - j;
  429.         Real(a, i0, j0) = rcos;
  430.         Imag(a, i0, j0) = -rsin;
  431.     }
  432.     }
  433.     Imag(a, n / 2, 0) = 0;
  434.     Imag(a, 0, n / 2) = 0;
  435.     Imag(a, n / 2, n / 2) = 0;
  436.     for (i = 1; i <= n / 2 - 1; i++) {
  437.     for (j = 1; j <= n / 2 - 1; j++) {
  438.         phase = 2 * M_PI * (RANDOMI() & 0x7FFF) / 32767.0;
  439.         rad = pow((double) (i * i + j * j), -(h + 1) / 2) * gauss();
  440.         rcos = rad * cos(phase);
  441.         rsin = rad * sin(phase);
  442.         Real(a, i, n - j) = rcos;
  443.         Imag(a, i, n - j) = rsin;
  444.         Real(a, n - i, j) = rcos;
  445.         Imag(a, n - i, j) = -rsin;
  446.     }
  447.     }
  448.  
  449.     nsize[0] = 0;
  450.     nsize[1] = nsize[2] = n;    /* Dimension of frequency domain array */
  451.     fourn(a, nsize, 2, -1);    /* Take inverse 2D Fourier transform */
  452. }
  453.  
  454. Image *
  455. draw_landscape(int w, int h, int clouds)
  456. {
  457.     float *a;
  458.     double rmin = HUGE_VAL, rmax = -HUGE_VAL, rmean, rrange;
  459.     Image *output = ImageNew(w, h);
  460.     unsigned char *op = output->data;
  461.     int x, y, i, j;
  462.     unsigned char *cp, *ap;
  463.     double *u, *u1;
  464.     unsigned int *bxf, *bxc;
  465.  
  466.     width = w;
  467.     height = h;
  468.  
  469.     spectralsynth(&a, meshsize, 3.0 - fracdim);
  470.  
  471.     /* Apply power law scaling if non-unity scale is requested. */
  472.  
  473.     if (powscale != 1.0)
  474.     for (i = 0; i < meshsize; i++)
  475.         for (j = 0; j < meshsize; j++) {
  476.         double r = Real(a, i, j);
  477.  
  478.         if (r > 0)
  479.             Real(a, i, j) = pow(r, powscale);
  480.         }
  481.     /* Compute extrema for autoscaling. */
  482.     for (i = 0; i < meshsize; i++)
  483.     for (j = 0; j < meshsize; j++) {
  484.         double r = Real(a, i, j);
  485.  
  486.         rmin = MIN(rmin, r);
  487.         rmax = MAX(rmax, r);
  488.     }
  489.  
  490.     rmean = (rmin + rmax) / 2;
  491.     rrange = (rmax - rmin) / 2;
  492.     for (i = 0; i < meshsize; i++)
  493.     for (j = 0; j < meshsize; j++)
  494.         Real(a, i, j) = (Real(a, i, j) - rmean) / rrange;
  495.  
  496.     u = (double *) XtMalloc((unsigned int) (width * sizeof(double)));
  497.     u1 = (double *) XtMalloc((unsigned int) (width * sizeof(double)));
  498.     bxf = (unsigned int *) XtMalloc((unsigned int) width *
  499.                     sizeof(unsigned int));
  500.     bxc = (unsigned int *) XtMalloc((unsigned int) width *
  501.                     sizeof(unsigned int));
  502.  
  503.     /* Prescale the grid points into intensities. */
  504.  
  505.     cp = (unsigned char *) XtMalloc(meshsize * meshsize + 1);
  506.     ap = cp;
  507.     for (i = 0; i < meshsize; i++)
  508.     for (j = 0; j < meshsize; j++)
  509.         *ap++ = (255.0 * (Real(a, i, j) + 1.0)) / 2.0;
  510.     /* Hack: make an extra array entry */
  511.     *ap = ap[-1];
  512.  
  513.     /* Fill the screen from the computed intensity grid by mapping screen points
  514.        onto the grid, then calculating each pixel value by bilinear
  515.        interpolation from the surrounding grid points.  (N.b. the pictures would
  516.        undoubtedly look better when generated with small grids if a more
  517.        well-behaved interpolation were used.)
  518.  
  519.        Before  we get started, precompute the line-level interpolation
  520.        parameters and store them in an array so we don't  have  to  do
  521.        this every time around the inner loop. */
  522.  
  523. #define UPRJ(a,size) ((a)/((size)-1.0))
  524.  
  525.     for (x = 0; x < width; x++) {
  526.     double bx = (meshsize - 1) * UPRJ(x, width);
  527.  
  528.     bxf[x] = floor(bx);
  529.     bxc[x] = bxf[x] + 1;
  530.     u[x] = bx - bxf[x];
  531.     u1[x] = 1 - u[x];
  532.     }
  533.  
  534.     for (y = 0; y < height; y++) {
  535.     double t, t1, by;
  536.     int byf, byc;
  537.  
  538.     by = (meshsize - 1) * UPRJ(y, height);
  539.     byf = floor(by) * meshsize;
  540.     byc = byf + meshsize;
  541.     t = by - floor(by);
  542.     t1 = 1 - t;
  543.  
  544.     if (clouds) {
  545.         /* Render the FFT output as clouds. */
  546.  
  547.         for (x = 0; x < width; x++) {
  548.         double r = t1 * u1[x] * cp[byf + bxf[x]] +
  549.         t * u1[x] * cp[byc + bxf[x]] +
  550.         t * u[x] * cp[byc + bxc[x]] +
  551.         t1 * u[x] * cp[byf + bxc[x]];
  552.         unsigned char w = (r > 127.0) ? (255 * ((r - 127.0) / 128.0)) : 0;
  553.  
  554.         *op++ = w;    /* blue trough white */
  555.         *op++ = w;
  556.         *op++ = 255;
  557.         }
  558.     } else
  559.         for (x = 0; x < width; x++) {
  560.         double r = t1 * u1[x] * cp[byf + bxf[x]] +
  561.         t * u1[x] * cp[byc + bxf[x]] +
  562.         t * u[x] * cp[byc + bxc[x]] +
  563.         t1 * u[x] * cp[byf + bxc[x]];
  564.         int ir, ig, ib;
  565.  
  566.         if (r >= 128) {
  567.             int ix = ((r - 128) * (ELEMENTS(pgnd) - 1)) / 127;
  568.  
  569.             /* Land area.  Look up colour based on elevation from
  570.                precomputed colour map table. */
  571.  
  572.             ir = pgnd[ix][0];
  573.             ig = pgnd[ix][1];
  574.             ib = pgnd[ix][2];
  575.         } else {
  576.             /* Water.  Generate clouds above water based on elevation.  */
  577.             ir = ig = r > 64 ? (r - 64) * 4 : 0;
  578.             ib = 255;
  579.         }
  580.  
  581.         *op++ = ir;
  582.         *op++ = ig;
  583.         *op++ = ib;
  584.         }
  585.     }
  586.  
  587.     XtFree(cp);
  588.     XtFree((char *) u);
  589.     XtFree((char *) u1);
  590.     XtFree((char *) bxf);
  591.     XtFree((char *) bxc);
  592.     XtFree((char *) a);
  593.  
  594.     return output;
  595. }
  596.  
  597. #endif                /* FEATURE_FRACTAL */