home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1994 January / usenetsourcesnewsgroupsinfomagicjanuary1994.iso / sources / misc / volume26 / pbmplus / patch10dec91 / part03 < prev    next >
Encoding:
Text File  |  1991-12-14  |  55.2 KB  |  1,693 lines

  1. Newsgroups: comp.sources.misc
  2. From: jef@well.sf.ca.us (Jef Poskanzer)
  3. Subject:  v26i108:  pbmplus - Extended Portable Bitmap Toolkit, Patch10dec91, Part03/05
  4. Message-ID: <1991Dec15.014640.20443@sparky.imd.sterling.com>
  5. X-Md4-Signature: 9a3f6109d610901d81f976509fa1525e
  6. Date: Sun, 15 Dec 1991 01:46:40 GMT
  7. Approved: kent@sparky.imd.sterling.com
  8.  
  9. Submitted-by: jef@well.sf.ca.us (Jef Poskanzer)
  10. Posting-number: Volume 26, Issue 108
  11. Archive-name: pbmplus/patch10dec91/part03
  12. Environment: UNIX
  13. Patch-To: pbmplus: Volume 23, Issue 36-59
  14.  
  15. #!/bin/sh
  16. # do not concatenate these parts, unpack them in order with /bin/sh
  17. # file pgm/pgmcrater.c continued
  18. #
  19. if test ! -r _shar_seq_.tmp; then
  20.     echo 'Please unpack part 1 first!'
  21.     exit 1
  22. fi
  23. (read Scheck
  24.  if test "$Scheck" != 3; then
  25.     echo Please unpack part "$Scheck" next!
  26.     exit 1
  27.  else
  28.     exit 0
  29.  fi
  30. ) < _shar_seq_.tmp || exit 1
  31. if test ! -f _shar_wnt_.tmp; then
  32.     echo 'x - still skipping pgm/pgmcrater.c'
  33. else
  34. echo 'x - continuing file pgm/pgmcrater.c'
  35. sed 's/^X//' << 'SHAR_EOF' >> 'pgm/pgmcrater.c' &&
  36. X    pages 31 and 32 of:
  37. X
  38. X    Peitgen, H.-O., and Saupe, D. eds., The Science Of Fractal
  39. X        Images, New York: Springer Verlag, 1988.
  40. X
  41. X    The  mathematical  technique  used    to calculate crater radii that
  42. X    obey the proper area law distribution from a uniformly distributed
  43. X    pseudorandom sequence was developed by Rudy Rucker.
  44. X
  45. X    Permission    to  use, copy, modify, and distribute this software and
  46. X    its documentation  for  any  purpose  and  without    fee  is  hereby
  47. X    granted,  without any conditions or restrictions.  This software is
  48. X    provided "as is" without express or implied warranty.
  49. X
  50. X                PLUGWARE!
  51. X
  52. X    If you like this kind of stuff, you may also enjoy "James  Gleick's
  53. X    Chaos--The  Software"  for  MS-DOS,  available for $59.95 from your
  54. X    local software store or directly from Autodesk, Inc., Attn: Science
  55. X    Series,  2320  Marinship Way, Sausalito, CA 94965, USA.  Telephone:
  56. X    (800) 688-2344 toll-free or, outside the  U.S. (415)  332-2344  Ext
  57. X    4886.   Fax: (415) 289-4718.  "Chaos--The Software" includes a more
  58. X    comprehensive   fractal    forgery      generator    which    creates
  59. X    three-dimensional  landscapes  as  well as clouds and planets, plus
  60. X    five more modules which explore other aspects of Chaos.   The  user
  61. X    guide  of  more  than  200    pages includes an introduction by James
  62. X    Gleick and detailed explanations by Rudy Rucker of the  mathematics
  63. X    and algorithms used by each program.
  64. X
  65. */
  66. X
  67. #include "pgm.h"
  68. #include <math.h>
  69. X
  70. /* Definitions for obtaining random numbers. */
  71. X
  72. #define Cast(low, high) ((low)+((high)-(low)) * ((random() & 0x7FFF) / arand))
  73. X
  74. /*  Data types    */
  75. X
  76. typedef int Boolean;
  77. #define FALSE 0
  78. #define TRUE 1
  79. X
  80. #define V   (void)
  81. X
  82. /*  Display parameters    */
  83. X
  84. #define SCRX    screenxsize          /* Screen width */
  85. #define SCRY    screenysize          /* Screen height */
  86. #define SCRGAMMA 1.0              /* Display gamma */
  87. X
  88. /*  Local variables  */
  89. X
  90. #define ImageGamma  0.5           /* Inherent gamma of mapped image */
  91. X
  92. static int screenxsize = 256;          /* Screen X size */
  93. static int screenysize = 256;          /* Screen Y size */
  94. static double dgamma = SCRGAMMA;      /* Display gamma */
  95. static double arand = 32767.0;          /* Random number parameters */
  96. static long ncraters = 50000L;          /* Number of craters to generate */
  97. static double CdepthPower = 1.5;      /* Crater depth power factor */
  98. static double DepthBias = 0.707107;   /* Depth bias */
  99. X
  100. /*  INITSEED  --  Generate initial random seed, if needed.  */
  101. X
  102. static void initseed()
  103. {
  104. X    int i;
  105. X
  106. X    i = time((long *) 0) * 0xF37C;
  107. X    srandom(i);
  108. X    for (i = 0; i < 7; i++) {
  109. X    V random();
  110. X    }
  111. }
  112. X
  113. /*  GENCRATERS    --  Generate cratered terrain.    */
  114. X
  115. static void gencraters()
  116. {
  117. X    int i, j, x, y;
  118. X    long l;
  119. X    unsigned short *aux;
  120. X    int slopemin = -52, slopemax = 52;
  121. #define RGBQuant    255
  122. X    unsigned char *slopemap;   /* Slope to pixel map */
  123. X    gray *pixels;           /* Pixel vector */
  124. X
  125. #define Auxadr(x, y)  ((unsigned short *) (aux + ((((y)) * SCRX) + (x))))
  126. X
  127. X    /* Acquire the elevation array and initialise it to mean
  128. X       surface elevation. */
  129. X
  130. X    aux = (unsigned short *) malloc(SCRX * SCRY * sizeof(short));
  131. X    if (aux == (unsigned short *) 0) {
  132. X        pm_error("out of memory allocating elevation array");
  133. X    }
  134. X
  135. X    /* Acquire the elevation buffer and initialise to mean
  136. X       initial elevation. */
  137. X
  138. X    for (i = 0; i < SCRY; i++) {
  139. X    unsigned short *zax = aux + (((long) SCRX) * i);
  140. X
  141. X    for (j = 0; j < SCRX; j++) {
  142. X        *zax++ = 32767;
  143. X    }
  144. X    }
  145. X
  146. X    /* Every time we go around this loop we plop another crater
  147. X       on the surface.    */
  148. X
  149. X    for (l = 0; l < ncraters; l++) {
  150. X    double g;
  151. X    int cx = Cast(0.0, ((double) SCRX - 1)),
  152. X        cy = Cast(0.0, ((double) SCRY - 1)),
  153. X        gx, gy, x, y;
  154. X    unsigned long amptot = 0, axelev;
  155. X    unsigned int npatch = 0;
  156. X
  157. X
  158. X    /* Phase 1.  Compute the mean elevation of the impact
  159. X             area.  We assume the impact area is a
  160. X             fraction of the total crater size. */
  161. X
  162. X    /* Thanks, Rudy, for this equation  that maps the uniformly
  163. X       distributed    numbers  from    Cast   into   an   area-law
  164. X       distribution as observed on cratered bodies. */
  165. X
  166. X    g = sqrt(1 / (M_PI * (1 - Cast(0, 0.9999))));
  167. X
  168. X    /* If the crater is tiny, handle it specially. */
  169. X
  170. X    if (g < 3) {
  171. X
  172. X       /* Set pixel to the average of its Moore neighbourhood. */
  173. X
  174. X        for (y = max(0, cy - 1); y <= min(SCRY - 1, cy + 1); y++) {
  175. X        int sx = max(0, cx - 1);
  176. X        unsigned short *a = Auxadr(sx, y);
  177. X
  178. X        for (x = sx; x <= min(SCRX - 1, cx + 1); x++) {
  179. X            amptot += *a++;
  180. X            npatch++;
  181. X        }
  182. X        }
  183. X        axelev = amptot / npatch;
  184. X
  185. X        /* Perturb the mean elevation by a small random factor. */
  186. X
  187. X        x = (g >= 1) ? ((random() >> 8) & 3) - 1 : 0;
  188. X        *Auxadr(cx, cy) = axelev + x;
  189. X
  190. X        /* Jam repaint sizes to correct patch. */
  191. X
  192. X        gx = 1;
  193. X        gy = 0;
  194. X
  195. X    } else {
  196. X
  197. X        /* Regular crater.    Generate an impact feature of the
  198. X           correct size and shape. */
  199. X
  200. X        /* Determine mean elevation around the impact area. */
  201. X
  202. X        gx = max(2, (g / 3));
  203. X        gy = max(2, g / 3);
  204. X
  205. X        for (y = max(0, cy - gy); y <= min(SCRY - 1, cy + gy); y++) {
  206. X        int sx = max(0, cx - gx);
  207. X        unsigned short *a = Auxadr(sx, y);
  208. X
  209. X        for (x = sx; x <= min(SCRX - 1, cx + gx); x++) {
  210. X            amptot += *a++;
  211. X            npatch++;
  212. X        }
  213. X        }
  214. X        axelev = amptot / npatch;
  215. X
  216. X        gy = max(2, g);
  217. X        g = gy;
  218. X        gx = max(2, g);
  219. X
  220. X        for (y = max(0, cy - gy); y <= min(SCRY - 1, cy + gy); y++) {
  221. X        int sx = max(0, cx - gx);
  222. X        unsigned short *ax = Auxadr(sx, y);
  223. X        double dy = (cy - y) / (double) gy,
  224. X               dysq = dy * dy;
  225. X
  226. X        for (x = sx; x <= min(SCRX - 1, cx + gx); x++) {
  227. X            double dx = ((cx - x) / (double) gx),
  228. X               cd = (dx * dx) + dysq,
  229. X               cd2 = cd * 2.25,
  230. X               tcz = DepthBias - sqrt(fabs(1 - cd2)),
  231. X               cz = max((cd2 > 1) ? 0.0 : -10, tcz),
  232. X               roll, iroll;
  233. X            unsigned short av;
  234. X
  235. X            cz *= pow(g, CdepthPower);
  236. X            if (dy == 0 && dx == 0 && ((int) cz) == 0) {
  237. X               cz = cz < 0 ? -1 : 1;
  238. X            }
  239. X
  240. #define         rollmin 0.9
  241. X            roll = (((1 / (1 - min(rollmin, cd))) /
  242. X                 (1 / (1 - rollmin))) - (1 - rollmin)) / rollmin;
  243. X            iroll = 1 - roll;
  244. X
  245. X            av = (axelev + cz) * iroll + (*ax + cz) * roll;
  246. X            av = max(1000, min(64000, av));
  247. X            *ax++ = av;
  248. X        }
  249. X        }
  250. X     }
  251. X    if ((l % 5000) == 4999) {
  252. X        pm_message( "%ld craters generated of %ld (%ld%% done)",
  253. X        l + 1, ncraters, ((l + 1) * 100) / ncraters);
  254. X    }
  255. X    }
  256. X
  257. X    i = max((slopemax - slopemin) + 1, 1);
  258. X    slopemap = (unsigned char *) malloc(i * sizeof(unsigned char));
  259. X    if (slopemap == (unsigned char *) 0) {
  260. X        pm_error("out of memory allocating slope map");
  261. X    }
  262. X    for (i = slopemin; i <= slopemax; i++) {
  263. X
  264. X        /* Confused?   OK,  we're using the  left-to-right slope to
  265. X       calculate a shade based on the sine of  the    angle  with
  266. X       respect  to the vertical (light incident from the left).
  267. X       Then, with one exponentiation, we account for  both    the
  268. X       inherent   gamma   of   the     image    (ad-hoc),  and    the
  269. X       user-specified display gamma, using the identity:
  270. X
  271. X         (x^y)^z = (x^(y*z))             */
  272. X
  273. X    slopemap[i - slopemin] = i > 0 ?
  274. X        (128 + 127.0 *
  275. X        pow(sin((M_PI / 2) * i / slopemax),
  276. X               dgamma * ImageGamma)) :
  277. X        (128 - 127.0 *
  278. X        pow(sin((M_PI / 2) * i / slopemin),
  279. X               dgamma * ImageGamma));
  280. X    }
  281. X
  282. X    /* Generate the screen image. */
  283. X
  284. X    pgm_writepgminit(stdout, SCRX, SCRY, RGBQuant, FALSE);
  285. X    pixels = pgm_allocrow(SCRX);
  286. X
  287. X    for (y = 0; y < SCRY; y++) {
  288. X    unsigned short *ax = Auxadr(0, y);
  289. X    gray *pix = pixels;
  290. X
  291. X    for (x = 0; x < SCRX - 1; x++) {
  292. X        int j = ax[1] - ax[0];
  293. X        j = min(max(slopemin, j), slopemax);
  294. X        *pix++ = slopemap[j - slopemin];
  295. X        ax++;
  296. X    }
  297. X    pgm_writepgmrow(stdout, pixels, SCRX, RGBQuant, FALSE);
  298. X    }
  299. X    pm_close(stdout);
  300. X    pgm_freerow(pixels);
  301. X
  302. #undef Auxadr
  303. #undef Scradr
  304. X    free((char *) slopemap);
  305. X    free((char *) aux);
  306. }
  307. X
  308. /*  MAIN  --  Main program.  */
  309. X
  310. int main(argc, argv)
  311. X  int argc;
  312. X  char *argv[];
  313. {
  314. X    int i;
  315. X    Boolean gammaspec = FALSE, numspec = FALSE,
  316. X        widspec = FALSE, hgtspec = FALSE;
  317. X    char *usage = "[-number <n>] [-width|-xsize <w>]\n\
  318. X                  [-height|-ysize <h>] [-gamma <f>]";
  319. X
  320. X    DepthBias = sqrt(0.5);          /* Get exact value for depth bias */
  321. X
  322. X    pgm_init(&argc, argv);
  323. X
  324. X    i = 1;
  325. X    while ((i < argc) && (argv[i][0] == '-') && (argv[i][1] != '\0')) {
  326. X        if (pm_keymatch(argv[i], "-gamma", 2)) {
  327. X        if (gammaspec) {
  328. X                pm_error("already specified gamma correction");
  329. X        }
  330. X        i++;
  331. X            if ((i == argc) || (sscanf(argv[i], "%lf", &dgamma)  != 1))
  332. X        pm_usage(usage);
  333. X        if (dgamma <= 0.0) {
  334. X                pm_error("gamma correction must be greater than 0");
  335. X        }
  336. X        gammaspec = TRUE;
  337. X        } else if (pm_keymatch(argv[i], "-number", 2)) {
  338. X        if (numspec) {
  339. X                pm_error("already specified number of craters");
  340. X        }
  341. X        i++;
  342. X            if ((i == argc) || (sscanf(argv[i], "%ld", &ncraters) != 1))
  343. X        pm_usage(usage);
  344. X        if (ncraters <= 0) {
  345. X                pm_error("number of craters must be greater than 0!");
  346. X        }
  347. X        numspec = TRUE;
  348. X        } else if (pm_keymatch(argv[i], "-xsize", 2) ||
  349. X                   pm_keymatch(argv[i], "-width", 2)) {
  350. X        if (widspec) {
  351. X                pm_error("already specified a width/xsize");
  352. X        }
  353. X        i++;
  354. X            if ((i == argc) || (sscanf(argv[i], "%d", &screenxsize) != 1))
  355. X        pm_usage(usage);
  356. X        if (screenxsize <= 0) {
  357. X                pm_error("screen width must be greater than 0");
  358. X        }
  359. X        widspec = TRUE;
  360. X        } else if (pm_keymatch(argv[i], "-ysize", 2) ||
  361. X                   pm_keymatch(argv[i], "-height", 2)) {
  362. X        if (hgtspec) {
  363. X                pm_error("already specified a height/ysize");
  364. X        }
  365. X        i++;
  366. X            if ((i == argc) || (sscanf(argv[i], "%d", &screenysize) != 1))
  367. X        pm_usage(usage);
  368. X        if (screenxsize <= 0) {
  369. X                pm_error("screen height must be greater than 0");
  370. X        }
  371. X        hgtspec = TRUE;
  372. X    } else {
  373. X        pm_usage(usage);
  374. X    }
  375. X    i++;
  376. X    }
  377. X
  378. X    initseed();
  379. X    gencraters();
  380. X
  381. X    exit(0);
  382. }
  383. SHAR_EOF
  384. echo 'File pgm/pgmcrater.c is complete' &&
  385. chmod 0664 pgm/pgmcrater.c ||
  386. echo 'restore of pgm/pgmcrater.c failed'
  387. Wc_c="`wc -c < 'pgm/pgmcrater.c'`"
  388. test 10318 -eq "$Wc_c" ||
  389.     echo 'pgm/pgmcrater.c: original size 10318, current size' "$Wc_c"
  390. rm -f _shar_wnt_.tmp
  391. fi
  392. # ============= pnm/Imakefile ==============
  393. if test ! -d 'pnm'; then
  394.     echo 'x - creating directory pnm'
  395.     mkdir 'pnm'
  396. fi
  397. if test -f 'pnm/Imakefile' -a X"$1" != X"-c"; then
  398.     echo 'x - skipping pnm/Imakefile (File already exists)'
  399.     rm -f _shar_wnt_.tmp
  400. else
  401. > _shar_wnt_.tmp
  402. echo 'x - extracting pnm/Imakefile (Text)'
  403. sed 's/^X//' << 'SHAR_EOF' > 'pnm/Imakefile' &&
  404. /* Imakefile for pnm tools
  405. X *
  406. X * Copyright (C) 1991 Rainer Klute
  407. X *
  408. X * Permission to use, copy, modify, distribute, and sell this software and
  409. X * its documentation for any purpose is hereby granted without fee, provided
  410. X * that the above copyright notice appear in all copies and that both that
  411. X * copyright notice and this permission notice appear in supporting
  412. X * documentation, and that the copyright holder's name not be used in
  413. X * advertising or publicity pertaining to distribution of the software
  414. X * without specific, written prior permission. The copyright holder makes
  415. X * no representations about the suitability of this software for any
  416. X * purpose. It is provided "as is" without express or implied warranty.
  417. X */
  418. X
  419. #define LibPnm libpnm.a
  420. #define DepLibPnm LibPnm
  421. #include <../Pbmplus.tmpl>
  422. X
  423. #if BuildLibTiff
  424. X   CURRENTLIBS = $(LIBTIFF) $(LIBPNM) $(LIBPPM) $(LIBPGM) $(LIBPBM)
  425. CURRENTDEPLIBS = $(DEPLIBTIFF) $(DEPLIBPNM) $(DEPLIBPPM) $(DEPLIBPGM) $(DEPLIBPBM)
  426. X     INCLUDES  = -I.. -I$(PBMDIR) -I$(PGMDIR) -I$(PPMDIR) -I$(TIFFDIR)
  427. X       DEFINES = -DLIBTIFF
  428. X         MERGE = pnmmerge
  429. X      TIFFMAN1 = tifftopnm.1 pnmtotiff.1
  430. X      TIFFSRCS = tifftopnm.c pnmtotiff.c
  431. X      TIFFOBJS = tifftopnm.o pnmtotiff.o
  432. X      TIFFBINS = tifftopnm pnmtotiff
  433. #else
  434. X   CURRENTLIBS = $(LIBPNM) $(LIBPPM) $(LIBPGM) $(LIBPBM)
  435. CURRENTDEPLIBS = $(DEPLIBPNM) $(DEPLIBPPM) $(DEPLIBPGM) $(DEPLIBPBM)
  436. X     INCLUDES  = -I.. -I$(PBMDIR) -I$(PGMDIR) -I$(PPMDIR)
  437. X         MERGE = pnmmerge
  438. X      TIFFMAN1 = 
  439. X      TIFFSRCS = 
  440. X      TIFFOBJS = 
  441. X      TIFFBINS = 
  442. #endif
  443. X
  444. X          MAN1 = pnmarith.1 pnmcat.1 pnmconvol.1 pnmcrop.1 pnmcut.1 \
  445. X         pnmdepth.1 pnmenlarge.1 pnmfile.1 pnmflip.1 pnminvert.1 \
  446. X                 pnmnoraw.1 pnmpaste.1 pnmscale.1 pnmtile.1 pnmtops.1 \
  447. X                 pnmtorast.1 pnmtoxwd.1 rasttopnm.1 xwdtopnm.1 \
  448. X         pnmgamma.1 pnmrotate.1 pnmshear.1 \
  449. X         anytopnm.1 pnmindex.1 pnmmargin.1 pnmsmooth.1 \
  450. X         $(TIFFMAN1)
  451. X          MAN3 = libpnm.3
  452. X          MAN5 = pnm.5
  453. X
  454. X          SRCS = pnmarith.c pnmcat.c pnmconvol.c pnmcrop.c pnmcut.c \
  455. X                 pnmdepth.c pnmenlarge.c pnmfile.c pnmflip.c pnminvert.c \
  456. X                 pnmnoraw.c pnmpaste.c pnmscale.c pnmtile.c pnmtops.c \
  457. X                 pnmtorast.c pnmtoxwd.c rasttopnm.c xwdtopnm.c \
  458. X         pnmgamma.c pnmrotate.c pnmshear.c \
  459. X         $(TIFFSRCS)
  460. X
  461. X          OBJS = pnmarith.o pnmcat.o pnmconvol.o pnmcrop.o pnmcut.o \
  462. X                 pnmdepth.o pnmenlarge.o pnmfile.o pnmflip.o pnminvert.o \
  463. X                 pnmnoraw.o pnmpaste.o pnmscale.o pnmtile.o pnmtops.o \
  464. X                 pnmtorast.o pnmtoxwd.o rasttopnm.o xwdtopnm.o \
  465. X         pnmgamma.o pnmrotate.o pnmshear.o \
  466. X         $(TIFFOBJS)
  467. X
  468. X          BINS = pnmarith pnmcat pnmconvol pnmcrop pnmcut \
  469. X         pnmdepth pnmenlarge pnmfile pnmflip pnminvert \
  470. X         pnmnoraw pnmpaste pnmscale pnmtile pnmtops \
  471. X         pnmtorast pnmtoxwd rasttopnm xwdtopnm \
  472. X         pnmgamma pnmrotate pnmshear \
  473. X         $(TIFFBINS)
  474. X
  475. includes:: anytopnm.script pnmindex.script pnmmargin.script pnmsmooth.script
  476. X
  477. anytopnm.script:
  478. X    $(LN) anytopnm anytopnm.script
  479. X
  480. pnmindex.script:
  481. X    $(LN) pnmindex pnmindex.script
  482. X
  483. pnmmargin.script:
  484. X    $(LN) pnmmargin pnmmargin.script
  485. X
  486. pnmsmooth.script:
  487. X    $(LN) pnmsmooth pnmsmooth.script
  488. X
  489. AllTarget($(LIBPNM) $(BINS))
  490. X
  491. DependTarget()
  492. X
  493. NormalPbmplusProgramTarget(pnmarith)
  494. NormalPbmplusProgramTarget(pnmcat)
  495. NormalPbmplusProgramTarget(pnmconvol)
  496. NormalPbmplusProgramTarget(pnmcrop)
  497. NormalPbmplusProgramTarget(pnmcut)
  498. NormalPbmplusProgramTarget(pnmdepth)
  499. NormalPbmplusProgramTarget(pnmenlarge)
  500. NormalPbmplusProgramTarget(pnmfile)
  501. NormalPbmplusProgramTarget(pnmflip)
  502. NormalPbmplusProgramTarget(pnminvert)
  503. NormalPbmplusProgramTarget(pnmnoraw)
  504. NormalPbmplusProgramTarget(pnmpaste)
  505. NormalPbmplusProgramTarget(pnmscale)
  506. NormalPbmplusProgramTarget(pnmtile)
  507. NormalPbmplusProgramTarget(pnmtops)
  508. NormalPbmplusProgramTarget(pnmtorast)
  509. NormalPbmplusProgramTarget(pnmtoxwd)
  510. NormalPbmplusProgramTarget(rasttopnm)
  511. NormalPbmplusProgramTarget(xwdtopnm)
  512. NormalPbmplusMathProgramTarget(pnmgamma)
  513. NormalPbmplusMathProgramTarget(pnmrotate)
  514. NormalPbmplusMathProgramTarget(pnmshear)
  515. #if BuildLibTiff
  516. NormalPbmplusProgramTarget(tifftopnm)
  517. NormalPbmplusProgramTarget(pnmtotiff)
  518. #endif
  519. X
  520. NormalLibraryObjectRule()
  521. NormalLibraryTarget(pnm,libpnm1.o libpnm2.o libpnm3.o libpnm4.o)
  522. X
  523. #if InstallMerged
  524. NormalProgramTarget($(MERGE),$(MERGE).o $(OBJS),$(CURRENTDEPLIBS),$(CURRENTLIBS),-lm)
  525. #if InstallBinaries
  526. InstallProgram($(MERGE),$(PBMPLUSDIR)$(PBMPLUSBINDIR))
  527. #endif
  528. #endif
  529. X
  530. #if InstallBinaries
  531. InstallPbmplusPrograms($(BINS),$(PBMPLUSDIR)$(PBMPLUSBINDIR),$(INSTPGMFLAGS))
  532. InstallScript(anytopnm,$(PBMPLUSDIR)$(PBMPLUSBINDIR))
  533. InstallScript(pnmindex,$(PBMPLUSDIR)$(PBMPLUSBINDIR))
  534. InstallScript(pnmmargin,$(PBMPLUSDIR)$(PBMPLUSBINDIR))
  535. InstallScript(pnmsmooth,$(PBMPLUSDIR)$(PBMPLUSBINDIR))
  536. #endif
  537. X
  538. #if InstallManuals
  539. InstallMultipleMan($(MAN1),$(PBMPLUSDIR)$(PBMPLUSMANDIR)/man1)
  540. InstallMultipleMan($(MAN3),$(PBMPLUSDIR)$(PBMPLUSMANDIR)/man3)
  541. InstallMultipleMan($(MAN5),$(PBMPLUSDIR)$(PBMPLUSMANDIR)/man5)
  542. #endif
  543. X
  544. #if InstallLibraries
  545. InstallLibrary(pnm,$(PBMPLUSDIR)$(PBMPLUSLIBDIR))
  546. #endif
  547. X
  548. #if InstallIncludes
  549. InstallMultipleFlags(pnm.h,$(PBMPLUSDIR)$(PBMPLUSINCDIR),$(INSTINCFLAGS))
  550. #endif
  551. SHAR_EOF
  552. chmod 0664 pnm/Imakefile ||
  553. echo 'restore of pnm/Imakefile failed'
  554. Wc_c="`wc -c < 'pnm/Imakefile'`"
  555. test 5158 -eq "$Wc_c" ||
  556.     echo 'pnm/Imakefile: original size 5158, current size' "$Wc_c"
  557. rm -f _shar_wnt_.tmp
  558. fi
  559. # ============= ppm/Imakefile ==============
  560. if test ! -d 'ppm'; then
  561.     echo 'x - creating directory ppm'
  562.     mkdir 'ppm'
  563. fi
  564. if test -f 'ppm/Imakefile' -a X"$1" != X"-c"; then
  565.     echo 'x - skipping ppm/Imakefile (File already exists)'
  566.     rm -f _shar_wnt_.tmp
  567. else
  568. > _shar_wnt_.tmp
  569. echo 'x - extracting ppm/Imakefile (Text)'
  570. sed 's/^X//' << 'SHAR_EOF' > 'ppm/Imakefile' &&
  571. /* Imakefile for ppm tools
  572. X *
  573. X * Copyright (C) 1991 Rainer Klute
  574. X *
  575. X * Permission to use, copy, modify, distribute, and sell this software and
  576. X * its documentation for any purpose is hereby granted without fee, provided
  577. X * that the above copyright notice appear in all copies and that both that
  578. X * copyright notice and this permission notice appear in supporting
  579. X * documentation, and that the copyright holder's name not be used in
  580. X * advertising or publicity pertaining to distribution of the software
  581. X * without specific, written prior permission. The copyright holder makes
  582. X * no representations about the suitability of this software for any
  583. X * purpose. It is provided "as is" without express or implied warranty.
  584. X */
  585. X
  586. #define LibPpm libppm.a
  587. #define DepLibPpm LibPpm
  588. #include <../Pbmplus.tmpl>
  589. X
  590. X   CURRENTLIBS = $(LIBPPM) $(LIBPGM) $(LIBPBM)
  591. CURRENTDEPLIBS = $(DEPLIBPPM) $(DEPLIBPGM) $(DEPLIBPBM)
  592. X     INCLUDES  = -I.. -I$(PBMDIR) -I$(PGMDIR)
  593. X       DEFINES = -DRGB_DB=\"DefaultRGBDatabase\"
  594. X         MERGE = ppmmerge
  595. X
  596. X          MAN1 = giftoppm.1 gouldtoppm.1 ilbmtoppm.1 imgtoppm.1 mtvtoppm.1 \
  597. X         pcxtoppm.1 pgmtoppm.1 pi1toppm.1 picttoppm.1 \
  598. X         pjtoppm.1 ppmdither.1 ppmhist.1 ppmmake.1 \
  599. X         ppmquant.1 ppmrelief.1 ppmtoacad.1 ppmtogif.1 ppmtoicr.1 \
  600. X         ppmtoilbm.1 ppmtopcx.1 ppmtopgm.1 ppmtopi1.1 ppmtopict.1 \
  601. X         ppmtopj.1 ppmtopuzz.1 ppmtorgb3.1 ppmtosixel.1 \
  602. X         ppmtotga.1 ppmtouil.1 ppmtoxpm.1 ppmtoyuv.1 qrttoppm.1 \
  603. X         rawtoppm.1 rgb3toppm.1 sldtoppm.1 spctoppm.1 sputoppm.1 \
  604. X         tgatoppm.1 ximtoppm.1 xpmtoppm.1 yuvtoppm.1 \
  605. X         ppmforge.1 ppmpat.1 \
  606. X         ppmquantall.1
  607. X          MAN3 = libppm.3
  608. X          MAN5 = ppm.5
  609. X
  610. X          SRCS = giftoppm.c gouldtoppm.c ilbmtoppm.c imgtoppm.c mtvtoppm.c \
  611. X             pcxtoppm.c pgmtoppm.c pi1toppm.c picttoppm.c \
  612. X         pjtoppm.c ppmdither.c ppmhist.c ppmmake.c \
  613. X         ppmquant.c ppmrelief.c ppmtoacad.c ppmtogif.c ppmtoicr.c \
  614. X         ppmtoilbm.c ppmtopcx.c ppmtopgm.c ppmtopi1.c ppmtopict.c \
  615. X         ppmtopj.c ppmtopuzz.c ppmtorgb3.c ppmtosixel.c \
  616. X         ppmtotga.c ppmtouil.c ppmtoxpm.c ppmtoyuv.c qrttoppm.c \
  617. X         rawtoppm.c rgb3toppm.c sldtoppm.c spctoppm.c sputoppm.c \
  618. X             tgatoppm.c ximtoppm.c xpmtoppm.c yuvtoppm.c \
  619. X         ppmforge.c ppmpat.c
  620. X
  621. X          OBJS = giftoppm.o gouldtoppm.o ilbmtoppm.o imgtoppm.o mtvtoppm.o \
  622. X             pcxtoppm.o pgmtoppm.o pi1toppm.o picttoppm.o \
  623. X         pjtoppm.o ppmdither.o ppmhist.o ppmmake.o \
  624. X         ppmquant.o ppmrelief.o ppmtoacad.o ppmtogif.o ppmtoicr.o \
  625. X         ppmtoilbm.o ppmtopcx.o ppmtopgm.o ppmtopi1.o ppmtopict.o \
  626. X         ppmtopj.o ppmtopuzz.o ppmtorgb3.o ppmtosixel.o \
  627. X         ppmtotga.o ppmtouil.o ppmtoxpm.o ppmtoyuv.o qrttoppm.o \
  628. X         rawtoppm.o rgb3toppm.o sldtoppm.o spctoppm.o sputoppm.o \
  629. X             tgatoppm.o ximtoppm.o xpmtoppm.o yuvtoppm.o \
  630. X         ppmforge.o ppmpat.o
  631. X
  632. X          BINS = giftoppm gouldtoppm ilbmtoppm imgtoppm mtvtoppm \
  633. X         pcxtoppm pgmtoppm pi1toppm picttoppm \
  634. X         pjtoppm ppmdither ppmhist ppmmake \
  635. X         ppmquant ppmrelief ppmtoacad ppmtogif ppmtoicr \
  636. X         ppmtoilbm ppmtopcx ppmtopgm ppmtopi1 ppmtopict \
  637. X         ppmtopj ppmtopuzz ppmtorgb3 ppmtosixel \
  638. X         ppmtotga ppmtouil ppmtoxpm ppmtoyuv qrttoppm \
  639. X         rawtoppm rgb3toppm sldtoppm spctoppm sputoppm \
  640. X         tgatoppm ximtoppm xpmtoppm yuvtoppm \
  641. X         ppmforge ppmpat
  642. X
  643. includes:: ppmquantall.script
  644. X
  645. ppmquantall.script:
  646. X    $(LN) ppmquantall ppmquantall.script
  647. X
  648. AllTarget($(LIBPPM) $(BINS))
  649. X
  650. DependTarget()
  651. X
  652. NormalPbmplusProgramTarget(giftoppm)
  653. NormalPbmplusProgramTarget(gouldtoppm)
  654. NormalPbmplusProgramTarget(ilbmtoppm)
  655. NormalPbmplusProgramTarget(imgtoppm)
  656. NormalPbmplusProgramTarget(mtvtoppm)
  657. NormalPbmplusProgramTarget(pcxtoppm)
  658. NormalPbmplusProgramTarget(pgmtoppm)
  659. NormalPbmplusProgramTarget(pi1toppm)
  660. NormalPbmplusProgramTarget(picttoppm)
  661. NormalPbmplusProgramTarget(pjtoppm)
  662. NormalPbmplusProgramTarget(ppmdither)
  663. NormalPbmplusProgramTarget(ppmhist)
  664. NormalPbmplusProgramTarget(ppmmake)
  665. NormalPbmplusProgramTarget(ppmquant)
  666. NormalPbmplusProgramTarget(ppmrelief)
  667. NormalPbmplusProgramTarget(ppmtoacad)
  668. NormalPbmplusProgramTarget(ppmtogif)
  669. NormalPbmplusProgramTarget(ppmtoicr)
  670. NormalPbmplusProgramTarget(ppmtoilbm)
  671. NormalPbmplusProgramTarget(ppmtopcx)
  672. NormalPbmplusProgramTarget(ppmtopgm)
  673. NormalPbmplusProgramTarget(ppmtopi1)
  674. NormalPbmplusProgramTarget(ppmtopict)
  675. NormalPbmplusProgramTarget(ppmtopj)
  676. NormalPbmplusProgramTarget(ppmtopuzz)
  677. NormalPbmplusProgramTarget(ppmtorgb3)
  678. NormalPbmplusProgramTarget(ppmtosixel)
  679. NormalPbmplusProgramTarget(ppmtotga)
  680. NormalPbmplusProgramTarget(ppmtouil)
  681. NormalPbmplusProgramTarget(ppmtoxpm)
  682. NormalPbmplusProgramTarget(ppmtoyuv)
  683. NormalPbmplusProgramTarget(qrttoppm)
  684. NormalPbmplusProgramTarget(rawtoppm)
  685. NormalPbmplusProgramTarget(rgb3toppm)
  686. NormalPbmplusProgramTarget(sldtoppm)
  687. NormalPbmplusProgramTarget(spctoppm)
  688. NormalPbmplusProgramTarget(sputoppm)
  689. NormalPbmplusProgramTarget(tgatoppm)
  690. NormalPbmplusProgramTarget(ximtoppm)
  691. NormalPbmplusProgramTarget(xpmtoppm)
  692. NormalPbmplusProgramTarget(yuvtoppm)
  693. NormalPbmplusMathProgramTarget(ppmforge)
  694. NormalPbmplusMathProgramTarget(ppmpat)
  695. X
  696. NormalLibraryObjectRule()
  697. NormalLibraryTarget(ppm,libppm1.o libppm2.o libppm3.o libppm4.o libppm5.o)
  698. X
  699. #if InstallMerged
  700. NormalProgramTarget($(MERGE),$(MERGE).o $(OBJS),$(CURRENTDEPLIBS),$(CURRENTLIBS),-lm)
  701. #if InstallBinaries
  702. InstallProgram($(MERGE),$(PBMPLUSDIR)$(PBMPLUSBINDIR))
  703. #endif
  704. #endif
  705. X
  706. #if InstallBinaries
  707. InstallPbmplusPrograms($(BINS),$(PBMPLUSDIR)$(PBMPLUSBINDIR),$(INSTPGMFLAGS))
  708. InstallScript(ppmquantall,$(PBMPLUSDIR)$(PBMPLUSBINDIR))
  709. #endif
  710. X
  711. #if InstallManuals
  712. InstallMultipleMan($(MAN1),$(PBMPLUSDIR)$(PBMPLUSMANDIR)/man1)
  713. InstallMultipleMan($(MAN3),$(PBMPLUSDIR)$(PBMPLUSMANDIR)/man3)
  714. InstallMultipleMan($(MAN5),$(PBMPLUSDIR)$(PBMPLUSMANDIR)/man5)
  715. #endif
  716. X
  717. #if InstallLibraries
  718. InstallLibrary(ppm,$(PBMPLUSDIR)$(PBMPLUSLIBDIR))
  719. #endif
  720. X
  721. #if InstallIncludes
  722. InstallMultipleFlags(ppm.h,$(PBMPLUSDIR)$(PBMPLUSINCDIR),$(INSTINCFLAGS))
  723. #endif
  724. SHAR_EOF
  725. chmod 0664 ppm/Imakefile ||
  726. echo 'restore of ppm/Imakefile failed'
  727. Wc_c="`wc -c < 'ppm/Imakefile'`"
  728. test 5811 -eq "$Wc_c" ||
  729.     echo 'ppm/Imakefile: original size 5811, current size' "$Wc_c"
  730. rm -f _shar_wnt_.tmp
  731. fi
  732. # ============= ppm/ppmforge.1 ==============
  733. if test -f 'ppm/ppmforge.1' -a X"$1" != X"-c"; then
  734.     echo 'x - skipping ppm/ppmforge.1 (File already exists)'
  735.     rm -f _shar_wnt_.tmp
  736. else
  737. > _shar_wnt_.tmp
  738. echo 'x - extracting ppm/ppmforge.1 (Text)'
  739. sed 's/^X//' << 'SHAR_EOF' > 'ppm/ppmforge.1' &&
  740. .TH ppmforge 1 "25 October 1991"
  741. .IX ppmforge
  742. .IX fractals
  743. .IX clouds
  744. .IX planets
  745. .IX stars
  746. .SH NAME
  747. ppmforge - fractal forgeries of clouds, planets, and starry skies
  748. .SH SYNOPSIS
  749. .na
  750. .nh
  751. .B ppmforge
  752. .RB [ -clouds ]
  753. 'in +9n
  754. .RB [ -night ]
  755. .RB [ -dimension
  756. .IR  dimen ]
  757. .RB [ -hour
  758. .IR hour ]
  759. .RB [ -inclination|-tilt
  760. .IR angle ]
  761. .RB [ -mesh
  762. .IR size ]
  763. .RB [ -power
  764. .IR factor ]
  765. .RB [ -glaciers
  766. .IR level ]
  767. .RB [ -ice
  768. .IR level ]
  769. .RB [ -saturation
  770. .IR sat ]
  771. .RB [ -seed
  772. .IR seed ]
  773. .RB [ -stars
  774. .IR fraction ]
  775. .RB [ -xsize|-width
  776. .IR width ]
  777. .RB [ -ysize|-height
  778. .IR height ]
  779. .in -9n
  780. .ad
  781. .hy
  782. .SH DESCRIPTION
  783. .B ppmforge
  784. generates three kinds of ``random fractal forgeries,'' the term coined
  785. by Richard F. Voss of the IBM Thomas J. Watson Research Center for
  786. seemingly realistic pictures of natural objects generated by simple
  787. algorithms embodying randomness and fractal self-similarity.  The
  788. techniques used by
  789. .B ppmforge
  790. are essentially those
  791. given by Voss[1], particularly the technique of spectral synthesis
  792. explained in more detail by Dietmar Saupe[2].
  793. .PP
  794. The program generates two varieties of pictures: planets and clouds,
  795. which are just different renderings of data generated in an identical
  796. manner, illustrating the unity of the fractal structure of these very
  797. different objects.  A third type of picture, a starry sky, is
  798. synthesised directly from pseudorandom numbers.
  799. .PP
  800. The generation of planets or clouds begins with the preparation of an
  801. array of random data in the frequency domain.  The size of this
  802. array, the ``mesh size,'' can be set with the
  803. .B -mesh
  804. option; the larger the mesh the more realistic the pictures but the
  805. calculation time and memory requirement increases as the square of the
  806. mesh size.  The fractal dimension, which you can specify with the
  807. .B -dimension
  808. option, determines the roughness of the terrain on the planet or the
  809. scale of detail in the clouds.  As the fractal dimension is increased,
  810. more high frequency components are added into the random mesh.
  811. .PP
  812. Once the mesh is generated, an inverse two dimensional Fourier
  813. transform is performed upon it.  This converts the original random
  814. frequency domain data into spatial amplitudes.  We scale the real
  815. components that result from the Fourier transform into numbers from 0
  816. to 1 associated with each point on the mesh.  You can further
  817. modify this number by applying a ``power law scale'' to it with the
  818. .B -power
  819. option.   Unity scale
  820. leaves the numbers unmodified; a power scale of 0.5 takes the square
  821. root of the numbers in the mesh, while a power scale of 3 replaces the
  822. numbers in the mesh with their cubes.  Power law scaling is best
  823. envisioned by thinking of the data as representing the elevation of
  824. terrain; powers less than 1 yield landscapes with vertical scarps that
  825. look like glacially-carved valleys; powers greater than one make
  826. fairy-castle spires (which require large mesh sizes and high
  827. resolution for best results).
  828. .PP
  829. After these calculations, we have a array of the specified size
  830. containing numbers that range from 0 to 1.  The pixmaps are generated as
  831. follows:
  832. .TP 10
  833. .B Clouds
  834. A colour map is created that ranges from pure blue to white by
  835. increasing admixture (desaturation) of blue with white.  Numbers less
  836. than 0.5 are coloured blue, numbers between 0.5 and 1.0 are coloured
  837. with corresponding levels of white, with 1.0 being pure white.
  838. .TP
  839. .B Planet
  840. The mesh is projected onto a sphere.  Values less than 0.5 are treated
  841. as water and values between 0.5 and 1.0 as land.  The water areas are
  842. coloured based upon the water depth, and land based on its elevation.
  843. The random depth data are used to create clouds over the oceans.  An
  844. atmosphere approximately like the Earth's is simulated; its light
  845. absorption is calculated to create a blue cast around the limb of the
  846. planet.  A function that rises from 0 to 1 based on latitude is
  847. modulated by the local elevation to generate polar ice caps--high
  848. altitude terrain carries glaciers farther from the pole.  Based on the
  849. position of the star with respect to the observer, the apparent colour
  850. of each pixel of the planet is calculated by ray-tracing from the star
  851. to the planet to the observer and applying a lighting model that sums
  852. ambient light and diffuse reflection (for most planets ambient light
  853. is zero, as their primary star is the only source of illumination).
  854. Additional random data are used to generate stars around the planet.
  855. .TP
  856. .B Night
  857. A sequence of pseudorandom numbers is used to generate stars with a
  858. user specified density.
  859. .PP
  860. Cloud pictures always contain 256 or fewer colours and may be
  861. displayed on most colour mapped devices without further processing.
  862. Planet pictures often contain tens of thousands of colours which
  863. must be compressed with
  864. .B ppmquant
  865. or
  866. .B ppmdither
  867. before encoding in a colour mapped format.  If the display resolution is
  868. high enough,
  869. .B ppmdither
  870. generally produces better looking planets.
  871. .B ppmquant
  872. tends to create discrete colour bands, particularly in the oceans,
  873. which are unrealistic and distracting.  The number of colours in starry
  874. sky pictures generated with the
  875. .B -night
  876. option depends on the value specified for
  877. .BR -saturation .
  878. Small values limit the colour temperature distribution of the stars
  879. and reduce the number of colours in the image.
  880. If the
  881. .B -saturation
  882. is set to 0, none of the stars will be coloured and the resulting
  883. image will never contain more than 256 colours.
  884. Night sky pictures with many different star colours often look
  885. best when colour compressed by
  886. .B pnmdepth
  887. rather than
  888. .B ppmquant
  889. or
  890. .BR ppmdither .
  891. Try
  892. .I newmaxval
  893. settings of 63, 31, or 15 with
  894. .B pnmdepth
  895. to reduce the number of colours in the picture to 256 or fewer.
  896. .SH OPTIONS
  897. .TP 10
  898. .B -clouds
  899. Generate clouds.  A pixmap of fractal clouds is generated.  Selecting clouds
  900. sets the default for fractal dimension to 2.15 and power scale factor
  901. to 0.75.
  902. .TP
  903. .BI -dimension " dimen"
  904. Sets the fractal dimension to the specified
  905. .IR dimen ,
  906. which may be any floating point value between 0 and 3.  Higher fractal
  907. dimensions create more ``chaotic'' images, which require higher
  908. resolution output and a larger FFT mesh size to look good.  If no
  909. dimension is specified, 2.4 is used when generating planets and 2.15
  910. for clouds.
  911. .TP
  912. .BI -glaciers " level"
  913. The floating point
  914. .I level
  915. setting controls the extent to which terrain elevation causes ice to
  916. appear at lower latitudes.  The default value of 0.75 makes the polar
  917. caps extend toward the equator across high terrain and forms glaciers
  918. in the highest mountains, as on Earth.  Higher values make ice sheets
  919. that cover more and more of the land surface, simulating planets in the
  920. midst of an ice age.  Lower values tend to be boring, resulting in
  921. unrealistic geometrically-precise ice cap boundaries.
  922. .TP
  923. .BI -hour " hour"
  924. When generating a planet,
  925. .I hour
  926. is used as the ``hour angle at the central meridian.''  If you specify
  927. .BR "-hour 12" ,
  928. for example, the planet will be fully illuminated, corresponding to
  929. high noon at the longitude at the centre of the screen.  You can
  930. specify any floating point value between 0 and 24 for
  931. .IR hour ,
  932. but values which place most of the planet in darkness (0 to 4 and 20
  933. to 24) result in crescents which, while pretty, don't give you many
  934. illuminated pixels for the amount of computing that's required.  If no
  935. .B -hour
  936. option is specified, a random hour angle is chosen, biased so that
  937. only 25% of the images generated will be crescents.
  938. .TP
  939. .BI -ice " level"
  940. Sets the extent of the polar ice caps to the given floating point
  941. .IR level .
  942. The default level of 0.4 produces ice caps similar to those of the Earth.
  943. Smaller values reduce the amount of ice, while larger
  944. .B -ice
  945. settings create more prominent ice caps.  Sufficiently large values,
  946. such as 100 or more, in conjunction with small settings for
  947. .B -glaciers
  948. (try 0.1) create ``ice balls'' like Europa.
  949. .TP
  950. .BI -inclination|-tilt " angle"
  951. The inclination angle of the planet with regard to its primary star is
  952. set to
  953. .IR angle ,
  954. which can be any floating point value from -90 to 90.  The inclination
  955. angle can be thought of as specifying, in degrees, the ``season'' the
  956. planet is presently experiencing or, more precisely, the latitude at
  957. which the star transits the zenith at local noon.  If 0, the planet
  958. is at equinox; the star is directly overhead at the equator.
  959. Positive values represent summer in the northern hemisphere, negative
  960. values summer in the southern hemisphere.  The Earth's inclination
  961. angle, for example, is about 23.5 at the June solstice, 0 at the
  962. equinoxes in March and September, and -23.5 at the December solstice.
  963. If no inclination angle is specified, a random value between -21.6 and
  964. 21.6 degrees is chosen.
  965. .TP
  966. .BI -mesh " size"
  967. A mesh of
  968. .IR size " by " size
  969. will be used for the fast Fourier transform (FFT).  Note that memory
  970. requirements and computation speed increase as the square of
  971. .IR size ;
  972. if you double the mesh size, the program will use four times the
  973. memory and run four times as long.  The default mesh is 256x256, which
  974. produces reasonably good looking pictures while using half a megabyte
  975. for the 256x256 array of single precision complex numbers
  976. required by the FFT.  On machines with limited memory capacity, you
  977. may have to reduce the mesh size to avoid running out of RAM.
  978. Increasing the mesh size produces better looking pictures; the
  979. difference becomes particularly noticeable when generating high
  980. resolution images with relatively high fractal dimensions (between 2.2
  981. and 3).
  982. .TP
  983. .B -night
  984. A starry sky is generated.  The stars are created by the same algorithm
  985. used for the stars that surround planet pictures, but the output
  986. consists exclusively of stars.
  987. .TP
  988. .BI -power " factor"
  989. Sets the ``power factor'' used to scale elevations synthesised from
  990. the FFT to
  991. .IR factor ,
  992. which can be any floating point number greater than zero.  If no
  993. factor is specified a default of 1.2 is used if a planet is being
  994. generated, or 0.75 if clouds are selected by the
  995. .B -clouds
  996. option.  The result of the FFT image synthesis is an array of elevation
  997. values between 0 and 1.  A non-unity power factor exponentiates each
  998. of these elevations to the specified power.  For example, a power
  999. factor of 2 squares each value, while a power factor of 0.5 replaces
  1000. each with its square root.  (Note that exponentiating values between 0
  1001. and 1 yields values that remain within that range.)  Power factors
  1002. less than 1 emphasise large-scale elevation changes at the expense of
  1003. small variations.  Power factors greater than 1 increase the roughness
  1004. of the terrain and, like high fractal dimensions, may require a larger
  1005. FFT mesh size and/or higher screen resolution to look good.
  1006. .TP
  1007. .BI -saturation " sat"
  1008. Controls the degree of colour saturation of the stars that surround planet
  1009. pictures and fill starry skies created with the
  1010. .B -night
  1011. option.  The default value of 125 creates stars which resemble the sky
  1012. as seen by the human eye from Earth's surface.  Stars are dim; only
  1013. the brightest activate the cones in the human retina, causing colour
  1014. to be perceived.  Higher values of
  1015. .I sat
  1016. approximate the appearance of stars from Earth orbit, where better
  1017. dark adaptation, absence of skyglow, and the concentration of light
  1018. from a given star onto a smaller area of the retina thanks to the lack
  1019. of atmospheric turbulence enhances the perception of colour.  Values
  1020. greater than 250 create ``science fiction'' skies that, while pretty,
  1021. don't occur in this universe.
  1022. .TP
  1023. Thanks to the inverse square law combined with Nature's love of
  1024. mediocrity, there are many, many dim stars for every bright one.
  1025. This population relationship is accurately reflected in the skies
  1026. created by
  1027. .BR ppmforge .
  1028. Dim, low mass stars live much longer than bright massive stars,
  1029. consequently there are many reddish stars for every blue giant.  This
  1030. relationship is preserved by
  1031. .BR ppmforge .
  1032. You can reverse the proportion, simulating the sky as seen in a starburst
  1033. galaxy, by specifying a negative
  1034. .I sat
  1035. value.
  1036. .TP
  1037. .BI -seed " num"
  1038. Sets the seed for the random number generator to the integer
  1039. .IR num .
  1040. The seed used to create each picture is displayed on standard output (unless
  1041. suppressed with the
  1042. .B -quiet
  1043. option).  Pictures generated with the same seed will be identical.  If no
  1044. .B -seed
  1045. is specified, a random seed derived from the date and time will be
  1046. chosen.  Specifying an explicit seed allows you to re-render a picture
  1047. you particularly like at a higher resolution or with different viewing
  1048. parameters.
  1049. .TP
  1050. .BI -stars " fraction"
  1051. Specifies the percentage of pixels, in tenths of a percent, which will
  1052. appear as stars, either surrounding a planet or filling the entire
  1053. frame if
  1054. .B -night
  1055. is specified.  The default
  1056. .I fraction
  1057. is 100.
  1058. .TP
  1059. .BI -xsize|-width " width"
  1060. Sets the width of the generated image to
  1061. .I width
  1062. pixels.  The default width is 256 pixels.  Images must be at least as
  1063. wide as they are high; if a width less than the height is specified,
  1064. it will be increased to equal the height.  If you must have a long
  1065. skinny pixmap, make a square one with
  1066. .BR ppmforge ,
  1067. then use
  1068. .B pnmcut
  1069. to extract a portion of the shape and size you require.
  1070. .TP
  1071. .BI -ysize|-height " height"
  1072. Sets the height of the generated image to
  1073. .I height
  1074. pixels.  The default height is 256 pixels.  If the height specified
  1075. exceeds the width, the width will be increased to equal the height.
  1076. .PP
  1077. All flags can be abbreviated to their shortest unique prefix.
  1078. .SH BUGS
  1079. .PP
  1080. The algorithms require the output pixmap to be at least as wide as it
  1081. is high, and the width to be an even number of pixels.  These
  1082. constraints are enforced by increasing the size of the requested
  1083. pixmap if necessary.
  1084. .PP
  1085. You may have to reduce the FFT mesh size on machines with 16 bit
  1086. integers and segmented pointer architectures.
  1087. .SH "SEE ALSO"
  1088. .PD
  1089. .BR pnmcut (1),
  1090. .BR pnmdepth (1),
  1091. .BR ppmdither (1),
  1092. .BR ppmquant (1),
  1093. .BR ppm (5)
  1094. .TP 5
  1095. [1] 
  1096. Voss, Richard F., ``Random Fractal Forgeries,'' in Earnshaw
  1097. et. al., Fundamental Algorithms for Computer Graphics, Berlin:
  1098. Springer-Verlag, 1985.
  1099. .TP
  1100. [2]
  1101. Peitgen, H.-O., and Saupe, D. eds., The Science Of Fractal Images,
  1102. New York: Springer Verlag, 1988.
  1103. .ne 10
  1104. .SH AUTHOR
  1105. .RS 5
  1106. .nf
  1107. John Walker
  1108. Autodesk SA
  1109. Avenue des Champs-Montants 14b
  1110. CH-2074 MARIN
  1111. Suisse/Schweiz/Svizzera/Svizra/Switzerland
  1112. .PD 0
  1113. .TP 9
  1114. Usenet:
  1115. kelvin@Autodesk.com
  1116. .TP
  1117. Fax:
  1118. 038/33 88 15
  1119. .TP
  1120. Voice:
  1121. 038/33 76 33
  1122. .fi
  1123. .RE
  1124. .PD
  1125. .PP
  1126. Permission to use, copy, modify, and distribute this software and its
  1127. documentation for any purpose and without fee is hereby granted,
  1128. without any conditions or restrictions.  This software is provided ``as
  1129. is'' without express or implied warranty.
  1130. .PP
  1131. .B PLUGWARE!
  1132. If you like this kind of stuff, you may also enjoy ``James Gleick's
  1133. Chaos--The Software'' for MS-DOS, available for $59.95 from your
  1134. local software store or directly from Autodesk, Inc., Attn: Science
  1135. Series, 2320 Marinship Way, Sausalito, CA 94965, USA.  Telephone:
  1136. (800) 688-2344 toll-free or, outside the U.S. (415) 332-2344 Ext
  1137. 4886.  Fax: (415) 289-4718.  ``Chaos--The Software'' includes a more
  1138. comprehensive fractal forgery generator which creates
  1139. three-dimensional landscapes as well as clouds and planets, plus five
  1140. more modules which explore other aspects of Chaos.  The user guide of
  1141. more than 200 pages includes an introduction by James Gleick and
  1142. detailed explanations by Rudy Rucker of the mathematics and algorithms
  1143. used by each program.
  1144. SHAR_EOF
  1145. chmod 0664 ppm/ppmforge.1 ||
  1146. echo 'restore of ppm/ppmforge.1 failed'
  1147. Wc_c="`wc -c < 'ppm/ppmforge.1'`"
  1148. test 15385 -eq "$Wc_c" ||
  1149.     echo 'ppm/ppmforge.1: original size 15385, current size' "$Wc_c"
  1150. rm -f _shar_wnt_.tmp
  1151. fi
  1152. # ============= ppm/ppmforge.c ==============
  1153. if test -f 'ppm/ppmforge.c' -a X"$1" != X"-c"; then
  1154.     echo 'x - skipping ppm/ppmforge.c (File already exists)'
  1155.     rm -f _shar_wnt_.tmp
  1156. else
  1157. > _shar_wnt_.tmp
  1158. echo 'x - extracting ppm/ppmforge.c (Text)'
  1159. sed 's/^X//' << 'SHAR_EOF' > 'ppm/ppmforge.c' &&
  1160. /*
  1161. X
  1162. X        Fractal forgery generator for the PPM toolkit
  1163. X
  1164. X    Originally  designed  and  implemented    in December of 1989 by
  1165. X    John Walker as a stand-alone program for the  Sun  and    MS-DOS
  1166. X    under  Turbo C.  Adapted in September of 1991 for use with Jef
  1167. X        Poskanzer's raster toolkit.
  1168. X
  1169. X    References cited in the comments are:
  1170. X
  1171. X        Foley, J. D., and Van Dam, A., Fundamentals of Interactive
  1172. X        Computer  Graphics,  Reading,  Massachusetts:  Addison
  1173. X        Wesley, 1984.
  1174. X
  1175. X        Peitgen, H.-O., and Saupe, D. eds., The Science Of Fractal
  1176. X        Images, New York: Springer Verlag, 1988.
  1177. X
  1178. X        Press, W. H., Flannery, B. P., Teukolsky, S. A., Vetterling,
  1179. X        W. T., Numerical Recipes In C, New Rochelle: Cambridge
  1180. X        University Press, 1988.
  1181. X
  1182. X    Author:
  1183. X        John Walker
  1184. X        Autodesk SA
  1185. X        Avenue des Champs-Montants 14b
  1186. X        CH-2074 MARIN
  1187. X        Switzerland
  1188. X        Usenet: kelvin@Autodesk.com
  1189. X        Fax:    038/33 88 15
  1190. X        Voice:  038/33 76 33
  1191. X
  1192. X    Permission    to  use, copy, modify, and distribute this software and
  1193. X    its documentation  for  any  purpose  and  without    fee  is  hereby
  1194. X    granted,  without any conditions or restrictions.  This software is
  1195. X    provided "as is" without express or implied warranty.
  1196. X
  1197. X                PLUGWARE!
  1198. X
  1199. X    If you like this kind of stuff, you may also enjoy "James  Gleick's
  1200. X    Chaos--The  Software"  for  MS-DOS,  available for $59.95 from your
  1201. X    local software store or directly from Autodesk, Inc., Attn: Science
  1202. X    Series,  2320  Marinship Way, Sausalito, CA 94965, USA.  Telephone:
  1203. X    (800) 688-2344 toll-free or, outside the  U.S. (415)  332-2344  Ext
  1204. X    4886.   Fax: (415) 289-4718.  "Chaos--The Software" includes a more
  1205. X    comprehensive   fractal    forgery      generator    which    creates
  1206. X    three-dimensional  landscapes  as  well as clouds and planets, plus
  1207. X    five more modules which explore other aspects of Chaos.   The  user
  1208. X    guide  of  more  than  200    pages includes an introduction by James
  1209. X    Gleick and detailed explanations by Rudy Rucker of the  mathematics
  1210. X    and algorithms used by each program.
  1211. X
  1212. */
  1213. X
  1214. #include "ppm.h"
  1215. #include <math.h>
  1216. X
  1217. #ifndef M_PI
  1218. #define M_PI    3.14159265358979323846
  1219. #endif
  1220. #ifndef M_E
  1221. #define M_E    2.7182818284590452354
  1222. #endif
  1223. X
  1224. /* Definitions used to address real and imaginary parts in a two-dimensional
  1225. X   array of complex numbers as stored by fourn(). */
  1226. X
  1227. #define Real(v, x, y)  v[1 + (((x) * meshsize) + (y)) * 2]
  1228. #define Imag(v, x, y)  v[2 + (((x) * meshsize) + (y)) * 2]
  1229. X
  1230. /* Co-ordinate indices within arrays. */
  1231. X
  1232. #define X     0
  1233. #define Y     1
  1234. #define Z     2
  1235. X
  1236. /* Definition for obtaining random numbers. */
  1237. X
  1238. #define nrand 4               /* Gauss() sample count */
  1239. #define Cast(low, high) ((low)+(((high)-(low)) * ((random() & 0x7FFF) / arand)))
  1240. X
  1241. /* Utility definition to get an  array's  element  count  (at  compile
  1242. X   time).   For  example:  
  1243. X
  1244. X       int  arr[] = {1,2,3,4,5};
  1245. X       ... 
  1246. X       printf("%d", ELEMENTS(arr));
  1247. X
  1248. X   would print a five.  ELEMENTS("abc") can also be used to  tell  how
  1249. X   many  bytes are in a string constant INCLUDING THE TRAILING NULL. */
  1250. X   
  1251. #define ELEMENTS(array) (sizeof(array)/sizeof((array)[0]))
  1252. X
  1253. /*  Data types    */
  1254. X
  1255. typedef int Boolean;
  1256. #define FALSE 0
  1257. #define TRUE 1
  1258. X
  1259. #define V     (void)
  1260. X
  1261. /*  Display parameters    */
  1262. X
  1263. #define SCRSAT     255              /* Screen maximum intensity */
  1264. X
  1265. /*  Local variables  */
  1266. X
  1267. static double arand, gaussadd, gaussfac; /* Gaussian random parameters */
  1268. static double fracdim;              /* Fractal dimension */
  1269. static double powscale;           /* Power law scaling exponent */
  1270. static int meshsize = 256;          /* FFT mesh size */
  1271. static unsigned int rseed;          /* Current random seed */
  1272. static Boolean seedspec = FALSE;      /* Did the user specify a seed ? */
  1273. static Boolean clouds = FALSE;          /* Just generate clouds */
  1274. static Boolean stars = FALSE;          /* Just generate stars */
  1275. static int screenxsize = 256;          /* Screen X size */
  1276. static int screenysize = 256;          /* Screen Y size */
  1277. static double inclangle, hourangle;   /* Star position relative to planet */
  1278. static Boolean inclspec = FALSE;      /* No inclination specified yet */
  1279. static Boolean hourspec = FALSE;      /* No hour specified yet */
  1280. static double icelevel;           /* Ice cap theshold */
  1281. static double glaciers;           /* Glacier level */
  1282. static int starfraction;          /* Star fraction */
  1283. static int starcolour;              /* Star colour saturation */
  1284. X
  1285. /*    FOURN  --  Multi-dimensional fast Fourier transform
  1286. X
  1287. X    Called with arguments:
  1288. X
  1289. X       data       A  one-dimensional  array  of  floats  (NOTE!!!    NOT
  1290. X              DOUBLES!!), indexed from one (NOTE!!!   NOT  ZERO!!),
  1291. X              containing  pairs of numbers representing the complex
  1292. X              valued samples.  The Fourier transformed results    are
  1293. X              returned in the same array.
  1294. X
  1295. X       nn          An  array specifying the edge size in each dimension.
  1296. X              THIS ARRAY IS INDEXED FROM  ONE,    AND  ALL  THE  EDGE
  1297. X              SIZES MUST BE POWERS OF TWO!!!
  1298. X
  1299. X       ndim       Number of dimensions of FFT to perform.  Set to 2 for
  1300. X              two dimensional FFT.
  1301. X
  1302. X       isign      If 1, a Fourier transform is done; if -1 the  inverse
  1303. X              transformation is performed.
  1304. X
  1305. X        This  function  is essentially as given in Press et al., "Numerical
  1306. X        Recipes In C", Section 12.11, pp.  467-470.
  1307. */
  1308. X
  1309. static void fourn(data, nn, ndim, isign)
  1310. X  float data[];
  1311. X  int nn[], ndim, isign;
  1312. {
  1313. X    register int i1, i2, i3;
  1314. X    int i2rev, i3rev, ip1, ip2, ip3, ifp1, ifp2;
  1315. X    int ibit, idim, k1, k2, n, nprev, nrem, ntot;
  1316. X    float tempi, tempr;
  1317. X    double theta, wi, wpi, wpr, wr, wtemp;
  1318. X
  1319. #define SWAP(a,b) tempr=(a); (a) = (b); (b) = tempr
  1320. X
  1321. X    ntot = 1;
  1322. X    for (idim = 1; idim <= ndim; idim++)
  1323. X    ntot *= nn[idim];
  1324. X    nprev = 1;
  1325. X    for (idim = ndim; idim >= 1; idim--) {
  1326. X    n = nn[idim];
  1327. X    nrem = ntot / (n * nprev);
  1328. X    ip1 = nprev << 1;
  1329. X    ip2 = ip1 * n;
  1330. X    ip3 = ip2 * nrem;
  1331. X    i2rev = 1;
  1332. X    for (i2 = 1; i2 <= ip2; i2 += ip1) {
  1333. X        if (i2 < i2rev) {
  1334. X        for (i1 = i2; i1 <= i2 + ip1 - 2; i1 += 2) {
  1335. X            for (i3 = i1; i3 <= ip3; i3 += ip2) {
  1336. X            i3rev = i2rev + i3 - i2;
  1337. X            SWAP(data[i3], data[i3rev]);
  1338. X            SWAP(data[i3 + 1], data[i3rev + 1]);
  1339. X            }
  1340. X        }
  1341. X        }
  1342. X        ibit = ip2 >> 1;
  1343. X        while (ibit >= ip1 && i2rev > ibit) {
  1344. X        i2rev -= ibit;
  1345. X        ibit >>= 1;
  1346. X        }
  1347. X        i2rev += ibit;
  1348. X    }
  1349. X    ifp1 = ip1;
  1350. X    while (ifp1 < ip2) {
  1351. X        ifp2 = ifp1 << 1;
  1352. X        theta = isign * (M_PI * 2) / (ifp2 / ip1);
  1353. X        wtemp = sin(0.5 * theta);
  1354. X        wpr = -2.0 * wtemp * wtemp;
  1355. X        wpi = sin(theta);
  1356. X        wr = 1.0;
  1357. X        wi = 0.0;
  1358. X        for (i3 = 1; i3 <= ifp1; i3 += ip1) {
  1359. X        for (i1 = i3; i1 <= i3 + ip1 - 2; i1 += 2) {
  1360. X            for (i2 = i1; i2 <= ip3; i2 += ifp2) {
  1361. X            k1 = i2;
  1362. X            k2 = k1 + ifp1;
  1363. X            tempr = wr * data[k2] - wi * data[k2 + 1];
  1364. X            tempi = wr * data[k2 + 1] + wi * data[k2];
  1365. X            data[k2] = data[k1] - tempr;
  1366. X            data[k2 + 1] = data[k1 + 1] - tempi;
  1367. X            data[k1] += tempr;
  1368. X            data[k1 + 1] += tempi;
  1369. X            }
  1370. X        }
  1371. X        wr = (wtemp = wr) * wpr - wi * wpi + wr;
  1372. X        wi = wi * wpr + wtemp * wpi + wi;
  1373. X        }
  1374. X        ifp1 = ifp2;
  1375. X    }
  1376. X    nprev *= n;
  1377. X    }
  1378. }
  1379. #undef SWAP
  1380. X
  1381. /*  INITGAUSS  --  Initialise random number generators.  As given in
  1382. X           Peitgen & Saupe, page 77. */
  1383. X
  1384. static void initgauss(seed)
  1385. X  unsigned int seed;
  1386. {
  1387. X    /* Range of random generator */
  1388. X    arand = pow(2.0, 15.0) - 1.0;
  1389. X    gaussadd = sqrt(3.0 * nrand);
  1390. X    gaussfac = 2 * gaussadd / (nrand * arand);
  1391. X    srandom(seed);
  1392. }
  1393. X
  1394. /*  GAUSS  --  Return a Gaussian random number.  As given in Peitgen
  1395. X           & Saupe, page 77. */
  1396. X
  1397. static double gauss()
  1398. {
  1399. X    int i;
  1400. X    double sum = 0.0;
  1401. X
  1402. X    for (i = 1; i <= nrand; i++) {
  1403. X    sum += (random() & 0x7FFF);
  1404. X    }
  1405. X    return gaussfac * sum - gaussadd;
  1406. }
  1407. X
  1408. /*  SPECTRALSYNTH  --  Spectrally  synthesised    fractal  motion in two
  1409. X               dimensions.  This algorithm is given under  the
  1410. X               name   SpectralSynthesisFM2D  on  page  108  of
  1411. X               Peitgen & Saupe. */
  1412. X
  1413. static void spectralsynth(x, n, h)
  1414. X  float **x;
  1415. X  unsigned int n;
  1416. X  double h;
  1417. {
  1418. X    unsigned bl;
  1419. X    int i, j, i0, j0, nsize[3];
  1420. X    double rad, phase, rcos, rsin;
  1421. X    float *a;
  1422. X
  1423. X    bl = ((((unsigned long) n) * n) + 1) * 2 * sizeof(float);
  1424. X    a = (float *) calloc(bl, 1);
  1425. X    if (a == (float *) 0) {
  1426. X        pm_error("Cannot allocate %d x %d result array (%ld bytes).",
  1427. X       n, n, bl);
  1428. X    }
  1429. X    *x = a;
  1430. X
  1431. X    for (i = 0; i <= n / 2; i++) {
  1432. X    for (j = 0; j <= n / 2; j++) {
  1433. X        phase = 2 * M_PI * ((random() & 0x7FFF) / arand);
  1434. X        if (i != 0 || j != 0) {
  1435. X        rad = pow((double) (i * i + j * j), -(h + 1) / 2) * gauss();
  1436. X        } else {
  1437. X        rad = 0;
  1438. X        }
  1439. X        rcos = rad * cos(phase);
  1440. X        rsin = rad * sin(phase);
  1441. X        Real(a, i, j) = rcos;
  1442. X        Imag(a, i, j) = rsin;
  1443. X        i0 = (i == 0) ? 0 : n - i;
  1444. X        j0 = (j == 0) ? 0 : n - j;
  1445. X        Real(a, i0, j0) = rcos;
  1446. X        Imag(a, i0, j0) = - rsin;
  1447. X    }
  1448. X    }
  1449. X    Imag(a, n / 2, 0) = 0;
  1450. X    Imag(a, 0, n / 2) = 0;
  1451. X    Imag(a, n / 2, n / 2) = 0;
  1452. X    for (i = 1; i <= n / 2 - 1; i++) {
  1453. X    for (j = 1; j <= n / 2 - 1; j++) {
  1454. X        phase = 2 * M_PI * ((random() & 0x7FFF) / arand);
  1455. X        rad = pow((double) (i * i + j * j), -(h + 1) / 2) * gauss();
  1456. X        rcos = rad * cos(phase);
  1457. X        rsin = rad * sin(phase);
  1458. X        Real(a, i, n - j) = rcos;
  1459. X        Imag(a, i, n - j) = rsin;
  1460. X        Real(a, n - i, j) = rcos;
  1461. X        Imag(a, n - i, j) = - rsin;
  1462. X    }
  1463. X    }
  1464. X
  1465. X    nsize[0] = 0;
  1466. X    nsize[1] = nsize[2] = n;          /* Dimension of frequency domain array */
  1467. X    fourn(a, nsize, 2, -1);          /* Take inverse 2D Fourier transform */
  1468. }
  1469. X
  1470. /*  INITSEED  --  Generate initial random seed, if needed.  */
  1471. X
  1472. static void initseed()
  1473. {
  1474. X    int i = time((long *) 0) ^ 0xF37C;
  1475. X    V srandom(i);
  1476. X    for (i = 0; i < 7; i++)
  1477. X    V random();
  1478. X    rseed = random();
  1479. }
  1480. X
  1481. /*  TEMPRGB  --  Calculate the relative R, G, and B components    for  a
  1482. X         black    body  emitting    light  at a given temperature.
  1483. X         The Planck radiation equation is solved directly  for
  1484. X         the R, G, and B wavelengths defined for the CIE  1931
  1485. X         Standard    Colorimetric    Observer.      The    colour
  1486. X         temperature is specified in degrees Kelvin. */
  1487. X
  1488. static void temprgb(temp, r, g, b)
  1489. X  double temp;
  1490. X  double *r, *g, *b;
  1491. {
  1492. X    double c1 = 3.7403e10,
  1493. X       c2 = 14384.0,
  1494. X       er, eg, eb, es;
  1495. X
  1496. /* Lambda is the wavelength in microns: 5500 angstroms is 0.55 microns. */
  1497. X
  1498. #define Planck(lambda)  ((c1 * pow((double) lambda, -5.0)) /  \
  1499. X             (pow(M_E, c2 / (lambda * temp)) - 1))
  1500. X
  1501. X    er = Planck(0.7000);
  1502. X    eg = Planck(0.5461);
  1503. X    eb = Planck(0.4358);
  1504. #undef Planck
  1505. X
  1506. X    es = 1.0 / max(er, max(eg, eb));
  1507. X
  1508. X    *r = er * es;
  1509. X    *g = eg * es;
  1510. X    *b = eb * es;
  1511. }
  1512. X
  1513. /*  ETOILE  --    Set a pixel in the starry sky.    */
  1514. X
  1515. static void etoile(pix)
  1516. X  pixel *pix;
  1517. {
  1518. X    if ((random() % 1000) < starfraction) {
  1519. #define StarQuality    0.5          /* Brightness distribution exponent */
  1520. #define StarIntensity    8          /* Brightness scale factor */
  1521. #define StarTintExp    0.5          /* Tint distribution exponent */
  1522. X    double v = StarIntensity * pow(1 / (1 - Cast(0, 0.9999)),
  1523. X                       (double) StarQuality),
  1524. X           temp,
  1525. X           r, g, b;
  1526. X
  1527. X    if (v > 255) {
  1528. X        v = 255;
  1529. X    }
  1530. X
  1531. X    /* We make a special case for star colour  of zero in order to
  1532. X       prevent  floating  point  roundoff  which  would  otherwise
  1533. X       result  in  more  than  256 star colours.  We can guarantee
  1534. X       that if you specify no star colour, you never get more than
  1535. X       256 shades in the image. */
  1536. X
  1537. X    if (starcolour == 0) {
  1538. X        int vi = v;
  1539. X
  1540. X        PPM_ASSIGN(*pix, vi, vi, vi);
  1541. X    } else {
  1542. X        temp = 5500 + starcolour *
  1543. X              pow(1 / (1 - Cast(0, 0.9999)), StarTintExp) *
  1544. X                  ((random() & 7) ? -1 : 1);
  1545. X        /* Constrain temperature to a reasonable value: >= 2600K 
  1546. X           (S Cephei/R Andromedae), <= 28,000 (Spica). */
  1547. X        temp = max(2600, min(28000, temp));
  1548. X        temprgb(temp, &r, &g, &b);
  1549. X        PPM_ASSIGN(*pix, (int) (r * v + 0.499),
  1550. X                 (int) (g * v + 0.499),
  1551. X                 (int) (b * v + 0.499));
  1552. X    }
  1553. X    } else {
  1554. X    PPM_ASSIGN(*pix, 0, 0, 0);
  1555. X    }
  1556. }
  1557. X
  1558. /*  GENPLANET  --  Generate planet from elevation array.  */
  1559. X
  1560. static void genplanet(a, n)
  1561. X  float *a;
  1562. X  unsigned int n;
  1563. {
  1564. X    int i, j;
  1565. X    unsigned char *cp, *ap;
  1566. X    double *u, *u1;
  1567. X    unsigned int *bxf, *bxc;
  1568. X
  1569. #define RGBQuant    255
  1570. X    pixel *pixels;              /* Pixel vector */
  1571. X    pixel *rpix;              /* Current pixel pointer */
  1572. X
  1573. #define Atthick 1.03              /* Atmosphere thickness as a percentage
  1574. X                                         of planet's diameter */
  1575. X    double athfac = sqrt(Atthick * Atthick - 1.0);
  1576. X    double sunvec[3];
  1577. X
  1578. X    Boolean flipped = FALSE;
  1579. X    double shang, siang;
  1580. X
  1581. X    ppm_writeppminit(stdout, screenxsize, screenysize,
  1582. X             (pixval) RGBQuant, FALSE);
  1583. X
  1584. X    if (!stars) {
  1585. X    u = (double *) malloc((unsigned int) (screenxsize * sizeof(double)));
  1586. X    u1 = (double *) malloc((unsigned int) (screenxsize * sizeof(double)));
  1587. X    bxf = (unsigned int *) malloc((unsigned int) screenxsize *
  1588. X          sizeof(unsigned int));
  1589. X    bxc = (unsigned int *) malloc((unsigned int) screenxsize *
  1590. X          sizeof(unsigned int));
  1591. X
  1592. X    if (u == (double *) 0 || u1 == (double *) 0 ||
  1593. X        bxf == (unsigned int *) 0 || bxc == (unsigned int *) 0) {
  1594. X            pm_error("Cannot allocate %d element interpolation tables.",
  1595. X             screenxsize);
  1596. X    }
  1597. X
  1598. X    /* Compute incident light direction vector. */
  1599. X
  1600. X    shang = hourspec ? hourangle : Cast(0, 2 * M_PI);
  1601. X    siang = inclspec ? inclangle : Cast(-M_PI * 0.12, M_PI * 0.12);
  1602. X
  1603. X    sunvec[X] = sin(shang) * cos(siang);
  1604. X    sunvec[Y] = sin(siang);
  1605. X    sunvec[Z] = cos(shang) * cos(siang);
  1606. X
  1607. X    /* Allow only 25% of random pictures to be crescents */
  1608. X
  1609. X    if (!hourspec && ((random() % 100) < 75)) {
  1610. X        flipped = sunvec[Z] < 0 ? TRUE : FALSE;
  1611. X        sunvec[Z] = abs(sunvec[Z]);
  1612. X    }
  1613. X        pm_message("%s: -seed %d -dimension %.2f -power %.2f -mesh %d",
  1614. X            clouds ? "clouds" : "planet",
  1615. X        rseed, fracdim, powscale, meshsize);
  1616. X    if (!clouds) {
  1617. X        pm_message(
  1618. X               "        -inclination %.0f -hour %d -ice %.2f -glaciers %.2f",
  1619. X           (siang * (180.0 / M_PI)),
  1620. X           (int) (((shang * (12.0 / M_PI)) + 12 +
  1621. X          (flipped ? 12 : 0)) + 0.5) % 24, icelevel, glaciers);
  1622. X            pm_message("        -stars %d -saturation %d.",
  1623. X        starfraction, starcolour);
  1624. X    }
  1625. X
  1626. X    /* Prescale the grid points into intensities. */
  1627. X
  1628. X    cp = (unsigned char *) malloc(n * n);
  1629. X    if (cp == (unsigned char *) 0)
  1630. X        return;
  1631. X    ap = cp;
  1632. X    for (i = 0; i < n; i++) {
  1633. X        for (j = 0; j < n; j++) {
  1634. X        *ap++ = (255.0 * (Real(a, i, j) + 1.0)) / 2.0;
  1635. X        }
  1636. X    }
  1637. X
  1638. X    /* Fill the screen from the computed  intensity  grid  by  mapping
  1639. X       screen  points onto the grid, then calculating each pixel value
  1640. X       by bilinear interpolation from  the    surrounding  grid  points.
  1641. X       (N.b. the pictures would undoubtedly look better when generated
  1642. X       with small grids if    a  more  well-behaved  interpolation  were
  1643. X       used.)
  1644. X
  1645. X       Before  we get started, precompute the line-level interpolation
  1646. X           parameters and store them in an array so we don't  have  to  do
  1647. X       this every time around the inner loop. */
  1648. X
  1649. #define UPRJ(a,size) ((a)/((size)-1.0))
  1650. X
  1651. X    for (j = 0; j < screenxsize; j++) {
  1652. X        double bx = (n - 1) * UPRJ(j, screenxsize);
  1653. X
  1654. X        bxf[j] = floor(bx);
  1655. X        bxc[j] = bxf[j] + 1;
  1656. X        u[j] = bx - bxf[j];
  1657. X        u1[j] = 1 - u[j];
  1658. X    }
  1659. X    } else {
  1660. X        pm_message("night: -seed %d -stars %d -saturation %d.",
  1661. X        rseed, starfraction, starcolour);
  1662. X    }
  1663. X
  1664. X    pixels = ppm_allocrow(screenxsize);
  1665. X    for (i = 0; i < screenysize; i++) {
  1666. X    double t, t1, by, dy, dysq, sqomdysq, icet, svx, svy, svz,
  1667. X           azimuth;
  1668. X    int byf, byc, lcos;
  1669. X
  1670. X    if (!stars) {              /* Skip all this setup if just stars */
  1671. X        by = (n - 1) * UPRJ(i, screenysize);
  1672. X        dy = 2 * (((screenysize / 2) - i) / ((double) screenysize));
  1673. X        dysq = dy * dy;
  1674. X        sqomdysq = sqrt(1.0 - dysq);
  1675. X        svx = sunvec[X];
  1676. X        svy = sunvec[Y] * dy;
  1677. X        svz = sunvec[Z] * sqomdysq;
  1678. X        byf = floor(by) * n;
  1679. SHAR_EOF
  1680. true || echo 'restore of ppm/ppmforge.c failed'
  1681. fi
  1682. echo 'End of  part 3'
  1683. echo 'File ppm/ppmforge.c is continued in part 4'
  1684. echo 4 > _shar_seq_.tmp
  1685. exit 0
  1686. exit 0 # Just in case...
  1687. -- 
  1688. Kent Landfield                   INTERNET: kent@sparky.IMD.Sterling.COM
  1689. Sterling Software, IMD           UUCP:     uunet!sparky!kent
  1690. Phone:    (402) 291-8300         FAX:      (402) 291-4362
  1691. Please send comp.sources.misc-related mail to kent@uunet.uu.net.
  1692.