src/sdf2img.c

   1 /*****************************************************************************
   2 **  This is part of the CTSim program
   3 **  Copyright (C) 1983-2000 Kevin Rosenberg
   4 **
   5 **  $Id: sdf2img.c,v 1.2 2000/05/08 20:02:32 kevin Exp $
   6 **  $Log: sdf2img.c,v $
   7 **  Revision 1.2  2000/05/08 20:02:32  kevin
   8 **  ANSI C changes
   9 **
  10 **  Revision 1.1.1.1  2000/04/28 13:02:44  kevin
  11 **  Initial CVS import for first public release
  12 **
  13 **
  14 **
  15 **  This program is free software; you can redistribute it and/or modify
  16 **  it under the terms of the GNU General Public License (version 2) as
  17 **  published by the Free Software Foundation.
  18 **
  19 **  This program is distributed in the hope that it will be useful,
  20 **  but WITHOUT ANY WARRANTY; without even the implied warranty of
  21 **  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  22 **  GNU General Public License for more details.
  23 **
  24 **  You should have received a copy of the GNU General Public License
  25 **  along with this program; if not, write to the Free Software
  26 **  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  27 ******************************************************************************/
  28 /* FILE
  29  *   sdf2img.c                          Convert an SDF file to a viewable format image
  30  */
  31
  32 #include "ct.h"
  33
  34 #if HAVE_PNG
  35 void sdf2d_to_png (IMAGE *im, char *outfile, int nxcell, int nycell, double densmin, double densmax);
  36 #endif
  37 #if HAVE_GIF
  38 void sdf2d_to_gif (IMAGE *im, char *outfile, int nxcell, int nycell, double densmin, double densmax);
  39 #endif
  40 void sdf2d_to_pgm (IMAGE *im, char *outfile, int nxcell, int nycell, double densmin, double densmax);
  41 void sdf2d_to_pgmasc (IMAGE *im, char *outfile, int nxcell, int nycell, double densmin, double densmax);
  42
  43 enum { O_SCALE, O_MIN, O_MAX, O_AUTO, O_CENTER, O_STATS, O_FORMAT, O_LABELS,
  44        O_HELP, O_VERBOSE, O_VERSION, O_DEBUG };
  45
  46 static struct option my_options[] =
  47 {
  48   {"scale", 1, 0, O_SCALE},
  49   {"min", 1, 0, O_MIN},
  50   {"max", 1, 0, O_MAX},
  51   {"auto", 1, 0, O_AUTO},
  52   {"center", 1, 0, O_CENTER},
  53   {"format", 1, 0, O_FORMAT},
  54   {"labels", 0, 0, O_LABELS},
  55   {"stats", 0, 0, O_STATS},
  56   {"verbose", 0, 0, O_VERBOSE},
  57   {"debug", 0, 0, O_DEBUG},
  58   {"help", 0, 0, O_HELP},
  59   {"version", 0, 0, O_VERSION},
  60   {0, 0, 0, 0}
  61 };
  62
  63 #define O_AUTO_FULL   1
  64 #define O_AUTO_STD0_1 2
  65 #define O_AUTO_STD0_5 3
  66 #define O_AUTO_STD1   4
  67 #define O_AUTO_STD2   5
  68 #define O_AUTO_STD3   6
  69
  70 #define O_AUTO_FULL_STR  "full"
  71 #define O_AUTO_STD0_1_STR  "std0.1"
  72 #define O_AUTO_STD0_5_STR  "std0.5"
  73 #define O_AUTO_STD1_STR  "std1"
  74 #define O_AUTO_STD2_STR  "std2"
  75 #define O_AUTO_STD3_STR  "std3"
  76
  77 #define O_CENTER_MEAN      1
  78 #define O_CENTER_MODE      2
  79 #define O_CENTER_MEAN_STR  "mean"
  80 #define O_CENTER_MODE_STR  "mode"
  81
  82 #define O_FORMAT_GIF       1
  83 #define O_FORMAT_PNG       2
  84 #define O_FORMAT_PGM       3
  85 #define O_FORMAT_PGMASC    4
  86 #define O_FORMAT_GIF_STR   "gif"
  87 #define O_FORMAT_PNG_STR   "png"
  88 #define O_FORMAT_PGM_STR   "pgm"
  89 #define O_FORMAT_PGMASC_STR "pgmasc"
  90
  91 void
  92 usage (const char *program)
  93 {
  94   fprintf(stdout, "usage: %s sdfname outfile [OPTIONS]\n", kbasename(program));
  95   fprintf(stdout, "Convert SDF2D file to an image file\n");
  96   fprintf(stdout, "\n");
  97   fprintf(stdout, "     sdfname    Name of input SDF file\n");
  98   fprintf(stdout, "     outfile    Name of output file\n");
  99   fprintf(stdout, "     --format   Output format\n");
 100   fprintf(stdout, "         pgm    PGM (portable graymap) format (default)\n");
 101   fprintf(stdout, "         pgmasc PGM (portable graymap) ASCII format\n");
 102 #ifdef HAVE_PNG
 103   fprintf(stdout, "         png    PNG format\n");
 104 #endif
 105 #ifdef HAVE_GIF
 106   fprintf(stdout, "         gif    GIF format\n");
 107 #endif
 108   fprintf(stdout, "     --center   Center of window\n");
 109   fprintf(stdout, "         mode   Mode is center of window (default)\n");
 110   fprintf(stdout, "         mean   Mean is center of window\n");
 111   fprintf(stdout, "     --auto     Set auto window\n");
 112   fprintf(stdout, "         full   Use full window (default)\n");
 113   fprintf(stdout, "         std0.1 Use 0.1 standard deviation about center\n");
 114   fprintf(stdout, "         std0.5 Use 0.5 standard deviation about center\n");
 115   fprintf(stdout, "         std1   Use one standard deviation about center\n");
 116   fprintf(stdout, "         std2   Use two standard deviations about center\n");
 117   fprintf(stdout, "         std3   Use three standard deviations about center\n");
 118   fprintf(stdout, "     --scale    Scaling factor for output size\n");
 119   fprintf(stdout, "     --min      Set minimum intensity\n");
 120   fprintf(stdout, "     --max      Set maximum intensity\n");
 121   fprintf(stdout, "     --stats    Print image statistics\n");
 122   fprintf(stdout, "     --labels   Print image labels\n");
 123   fprintf(stdout, "     --debug    Set debug mode\n");
 124   fprintf(stdout, "     --verbose  Set verbose mode\n");
 125   fprintf(stdout, "     --version  Print version\n");
 126   fprintf(stdout, "     --help     Print this help message\n");
 127   exit(1);
 128 }
 129
 130 int
 131 main (int argc, char *const argv[])
 132 {
 133   IMAGE *im;
 134   double densmin = HUGE_VAL, densmax = -HUGE_VAL;
 135   char *in_file, *out_file;
 136   int opt_verbose = 0;
 137   int opt_scale = 1;
 138   int opt_auto = O_AUTO_FULL;
 139   int opt_set_max = 0;
 140   int opt_set_min = 0;
 141   int opt_stats = 0;
 142   int opt_debug = 0;
 143   int opt_center = O_CENTER_MODE;
 144   int opt_format = O_FORMAT_PGM;
 145   int opt_labels = 0;
 146
 147   while (1)
 148     {
 149       int c = getopt_long (argc, argv, "", my_options, NULL);
 150       char *endptr = NULL;
 151       char *endstr;
 152
 153       if (c == -1)
 154         break;
 155
 156       switch (c)
 157         {
 158         case O_MIN:
 159           opt_set_min = 1;
 160           densmin = strtod(optarg, &endptr);
 161           endstr = optarg + strlen(optarg);
 162           if (endptr != endstr)
 163             {
 164               fprintf(stderr, "Error setting --min to %s\n", optarg);
 165               usage(argv[0]);
 166               exit(1);
 167             }
 168           break;
 169         case O_MAX:
 170           opt_set_max = 1;
 171           densmax = strtod(optarg, &endptr);
 172           endstr = optarg + strlen(optarg);
 173           if (endptr != endstr)
 174             {
 175               fprintf(stderr, "Error setting --max to %s\n", optarg);
 176               usage(argv[0]);
 177               exit(1);
 178             }
 179           break;
 180         case O_SCALE:
 181           opt_scale = strtol(optarg, &endptr, 10);
 182           endstr = optarg + strlen(optarg);
 183           if (endptr != endstr)
 184             {
 185               fprintf(stderr,"Error setting --scale to %s\n", optarg);
 186               usage(argv[0]);
 187               exit(1);
 188             }
 189           break;
 190         case O_AUTO:
 191           if (strcmp(optarg, O_AUTO_FULL_STR) == 0)
 192             opt_auto = O_AUTO_FULL;
 193           else if (strcmp(optarg, O_AUTO_STD1_STR) == 0)
 194             opt_auto = O_AUTO_STD1;
 195           else if (strcmp(optarg, O_AUTO_STD0_5_STR) == 0)
 196             opt_auto = O_AUTO_STD0_5;
 197           else if (strcmp(optarg, O_AUTO_STD0_1_STR) == 0)
 198             opt_auto = O_AUTO_STD0_1;
 199           else if (strcmp(optarg, O_AUTO_STD2_STR) == 0)
 200             opt_auto = O_AUTO_STD2;
 201           else if (strcmp(optarg, O_AUTO_STD3_STR) == 0)
 202             opt_auto = O_AUTO_STD3;
 203           else
 204             {
 205               fprintf(stderr, "Invalid auto mode %s\n", optarg);
 206               usage(argv[0]);
 207               exit(1);
 208             }
 209                 break;
 210         case O_CENTER:
 211           if (strcmp(optarg, O_CENTER_MEAN_STR) == 0)
 212             opt_center = O_CENTER_MEAN;
 213           else if (strcmp(optarg, O_CENTER_MODE_STR) == 0)
 214             opt_center = O_CENTER_MODE;
 215           else
 216             {
 217               fprintf(stderr, "Invalid center mode %s\n", optarg);
 218               usage(argv[0]);
 219               exit(1);
 220             }
 221                 break;
 222         case O_FORMAT:
 223           if (strcmp(optarg, O_FORMAT_PGM_STR) == 0)
 224             opt_format = O_FORMAT_PGM;
 225           else if (strcmp(optarg, O_FORMAT_PGMASC_STR) == 0)
 226             opt_format = O_FORMAT_PGMASC;
 227 #if HAVE_PNG
 228           else if (strcmp(optarg, O_FORMAT_PNG_STR) == 0)
 229             opt_format = O_FORMAT_PNG;
 230 #endif
 231 #if HAVE_GIF
 232           else if (strcmp(optarg, O_FORMAT_GIF_STR) == 0)
 233             opt_format = O_FORMAT_GIF;
 234 #endif
 235           else {
 236               fprintf(stderr, "Invalid format mode %s\n", optarg);
 237               usage(argv[0]);
 238               exit(1);
 239             }
 240                 break;
 241         case O_LABELS:
 242           opt_labels = 1;
 243           break;
 244         case O_VERBOSE:
 245           opt_verbose = 1;
 246           break;
 247         case O_DEBUG:
 248           opt_debug = 1;
 249           break;
 250         case O_STATS:
 251           opt_stats = 1;
 252           break;
 253         case O_VERSION:
 254 #ifdef VERSION
 255           fprintf(stdout, "Version %s\n", VERSION);
 256 #else
 257           fprintf(stderr, "Unknown version number");
 258 #endif
 259           exit(0);
 260         case O_HELP:
 261         case '?':
 262           usage(argv[0]);
 263           exit(0);
 264         default:
 265           usage(argv[0]);
 266           exit(1);
 267         }
 268     }
 269
 270   if (optind + 2 != argc)
 271     {
 272       usage(argv[0]);
 273       exit(1);
 274     }
 275
 276   in_file = argv[optind];
 277   out_file = argv[optind+1];
 278
 279   if ((im = image_load (in_file)) == NULL)
 280     {
 281       fprintf(stderr, "File %s does not exist\n", in_file);
 282       exit(1);
 283     }
 284
 285   if (opt_labels)
 286     {
 287       int nlabels = im->dfp_2d->dfp->num_label;
 288       int i;
 289
 290       for (i = 0; i < nlabels; i++)
 291         {
 292           SDF_BLK blk;
 293           if (sdf_read_label(&blk.lab, i, im->dfp_2d->dfp) != SDF_OK)
 294             break;
 295           if (blk.lab.label_type == LT_HISTORY) {
 296             fprintf(stdout, "History: %s\n", blk.lab.label_str);
 297             fprintf(stdout, "  calc time = %.2f secs\n", blk.lab.calc_time);
 298             fprintf(stdout, "  Timestamp = %s\n", td_str_tmdt(&blk.lab.timedate));
 299           } else if (blk.lab.label_type == LT_TITLE) {
 300             fprintf(stdout, "Title: %s\n", blk.lab.label_str);
 301             fprintf(stdout, "  Timestamp = %s\n", td_str_tmdt(&blk.lab.timedate));
 302           } else if (blk.lab.label_type == LT_NOTE) {
 303             fprintf(stdout, "Note: %s\n", blk.lab.label_str);
 304             fprintf(stdout, "  Timestamp = %s\n", td_str_tmdt(&blk.lab.timedate));
 305           }
 306         }
 307     }
 308
 309   if (opt_stats || (! (opt_set_max && opt_set_min)))
 310     {
 311       double minfound = HUGE_VAL, maxfound = -HUGE_VAL;
 312       double mode = 0, mean = 0, stddev = 0, window = 0;
 313       double spread;
 314       int hist[256];
 315       int ibin, nbin = 256;
 316       int max_bin, max_binindex;
 317       double maxbin;
 318       int ix, iy;
 319
 320       maxbin = nbin - 1;
 321       for (ix = 0; ix < im->nx; ix++)
 322         {
 323           for (iy = 0; iy < im->ny; iy++)
 324             {
 325               if (im->v[ix][iy] > maxfound)
 326                 maxfound = im->v[ix][iy];
 327               if (im->v[ix][iy] < minfound)
 328                 minfound = im->v[ix][iy];
 329               mean += im->v[ix][iy];
 330             }
 331         }
 332       spread = maxfound - minfound;
 333       if (spread == 0)
 334         mode = minfound;
 335       else
 336         {
 337           for (ibin = 0; ibin < nbin; ibin++)
 338             hist[ibin] = 0;
 339           for (ix = 0; ix < im->nx; ix++)
 340             {
 341               for (iy = 0; iy < im->ny; iy++)
 342                 {
 343                   int b = (int) ((((im->v[ix][iy] - minfound) / spread) * (double) maxbin) + 0.5);
 344                   hist[b]++;
 345                 }
 346             }
 347           max_binindex = 0;
 348           max_bin = -1;
 349           for (ibin = 0; ibin < nbin; ibin++)
 350             {
 351               if (hist[ibin] > max_bin)
 352                 {
 353                   max_bin = hist[ibin];
 354                   max_binindex = ibin;
 355                 }
 356             }
 357           mode = (((double) max_binindex) * spread / ((double) maxbin)) + minfound;
 358         }
 359
 360       if (opt_auto == O_AUTO_FULL)
 361         {
 362           if (! opt_set_max)
 363             densmax = maxfound;
 364           if (! opt_set_min)
 365             densmin = minfound;
 366         }
 367       if (opt_stats || opt_auto != O_AUTO_FULL)
 368         {
 369           mean /= im->nx * im->ny;
 370           for (ix = 0; ix < im->nx; ix++)
 371             {
 372               for (iy = 0; iy < im->ny; iy++)
 373                 {
 374                   double diff = (im->v[ix][iy] - mean);
 375                   stddev += diff * diff;
 376                 }
 377             }
 378           stddev = sqrt(stddev / (im->nx * im->ny));
 379           if (opt_auto == O_AUTO_FULL)
 380             ;
 381           else if (opt_auto == O_AUTO_STD1)
 382             window = stddev;
 383           else if (opt_auto == O_AUTO_STD0_1)
 384             window = stddev * 0.1;
 385           else if (opt_auto == O_AUTO_STD0_5)
 386             window = stddev * 0.5;
 387           else if (opt_auto == O_AUTO_STD2)
 388             window = stddev * 2;
 389           else if (opt_auto == O_AUTO_STD3)
 390             window = stddev * 3;
 391           else
 392             {
 393               fprintf(stderr, "Internal Error: Invalid auto mode %d\n", opt_auto);
 394               exit(1);
 395             }
 396         }
 397       if (opt_stats)
 398         {
 399           fprintf(stdout,"nx=%d\n", im->nx);
 400           fprintf(stdout,"ny=%d\n", im->ny);
 401           fprintf(stdout,"min=%f\n", minfound);
 402           fprintf(stdout,"max=%f\n", maxfound);
 403           fprintf(stdout,"mean=%f\n", mean);
 404           fprintf(stdout,"mode=%f\n", mode);
 405           fprintf(stdout,"stddev=%f\n", stddev);
 406         }
 407       if (opt_auto != O_AUTO_FULL)
 408         {
 409           double center;
 410
 411           if (opt_center == O_CENTER_MODE)
 412             center = mode;
 413           else if (opt_center == O_CENTER_MEAN)
 414             center = mean;
 415           else
 416             {
 417               fprintf(stderr, "Internal Error: Invalid center mode %d\n", opt_center);
 418               exit(1);
 419             }
 420           if (! opt_set_max)
 421             densmax = center + window;
 422           if (! opt_set_min)
 423             densmin = center - window;
 424         }
 425     }
 426
 427   if (opt_stats)
 428     {
 429       fprintf(stdout,"min_disp=%f\n", densmin);
 430       fprintf(stdout,"max_disp=%f\n", densmax);
 431     }
 432
 433   if (opt_format == O_FORMAT_PGM)
 434     sdf2d_to_pgm (im, out_file, opt_scale, opt_scale, densmin, densmax);
 435   else if (opt_format == O_FORMAT_PGMASC)
 436     sdf2d_to_pgmasc (im, out_file, opt_scale, opt_scale, densmin, densmax);
 437 #if HAVE_PNG
 438   else if (opt_format == O_FORMAT_PNG)
 439     sdf2d_to_png (im, out_file, opt_scale, opt_scale, densmin, densmax);
 440 #endif
 441 #if HAVE_GIF
 442   else if (opt_format == O_FORMAT_GIF)
 443     sdf2d_to_gif (im, out_file, opt_scale, opt_scale, densmin, densmax);
 444 #endif
 445   else
 446     {
 447       fprintf(stderr, "Internal Error: Invalid format mode %d\n", opt_format);
 448       exit(1);
 449     }
 450   return (0);
 451 }
 452
 453
 454 void
 455 sdf2d_to_pgm (IMAGE *im, char *outfile, int nxcell, int nycell, double densmin, double densmax)
 456 {
 457   FILE *fp;
 458   int irow;
 459   unsigned char *rowp;
 460
 461   rowp = (unsigned char *) sys_alloc(im->nx * nxcell,"sdf2d_to_img");
 462   if (rowp == NULL)
 463     return;
 464
 465   if ((fp = fopen (outfile, "wb")) == NULL)
 466      return;
 467
 468   fprintf(fp, "P5\n");
 469   fprintf(fp, "%d %d\n", im->nx, im->ny);
 470   fprintf(fp, "255\n");
 471
 472   for (irow = im->ny - 1; irow >= 0; irow--)
 473     {
 474       int icol, ir;
 475
 476       for (icol = 0; icol < im->nx; icol++)
 477         {
 478           double dens;
 479           int p;
 480           int pos = icol * nxcell;
 481           dens = (im->v[icol][irow] - densmin) / (densmax - densmin);
 482           if (dens < 0)
 483             dens = 0;
 484           else if (dens > 1)
 485             dens = 1;
 486           for (p = pos; p < pos + nxcell; p++)
 487             {
 488               rowp[p] = (unsigned int) (dens * 255.);
 489             }
 490         }
 491       for (ir = 0; ir < nycell; ir++)
 492         {
 493           int ic;
 494           for (ic = 0; ic < im->nx * nxcell; ic++)
 495             fprintf(fp, "%c ", rowp[ic]);
 496         }
 497     }
 498   sys_free(rowp, "sdf2img");
 499
 500   fclose(fp);
 501 }
 502
 503 void
 504 sdf2d_to_pgmasc (IMAGE *im, char *outfile, int nxcell, int nycell, double densmin, double densmax)
 505 {
 506   FILE *fp;
 507   int irow;
 508   unsigned char *rowp;
 509
 510   rowp = (unsigned char *) sys_alloc(im->nx * nxcell,"sdf2d_to_img");
 511   if (rowp == NULL)
 512     return;
 513
 514   if ((fp = fopen (outfile, "wb")) == NULL)
 515      return;
 516
 517   fprintf(fp, "P2\n");
 518   fprintf(fp, "%d %d\n", im->nx, im->ny);
 519   fprintf(fp, "255\n");
 520
 521   for (irow = im->ny - 1; irow >= 0; irow--)
 522     {
 523       int icol, ir;
 524
 525       for (icol = 0; icol < im->nx; icol++)
 526         {
 527           double dens;
 528           int p;
 529           int pos = icol * nxcell;
 530           dens = (im->v[icol][irow] - densmin) / (densmax - densmin);
 531           if (dens < 0)
 532             dens = 0;
 533           else if (dens > 1)
 534             dens = 1;
 535           for (p = pos; p < pos + nxcell; p++)
 536             {
 537               rowp[p] = (unsigned int) (dens * 255.);
 538             }
 539         }
 540       for (ir = 0; ir < nycell; ir++)
 541         {
 542           int ic;
 543           for (ic = 0; ic < im->nx * nxcell; ic++)
 544             fprintf(fp, "%d ", rowp[ic]);
 545           fprintf(fp, "\n");
 546         }
 547     }
 548   sys_free(rowp, "sdf2img");
 549
 550   fclose(fp);
 551 }
 552
 553
 554 #ifdef HAVE_PNG
 555 void
 556 sdf2d_to_png (IMAGE *im, char *outfile, int nxcell, int nycell, double densmin, double densmax)
 557 {
 558   FILE *fp;
 559   png_structp png_ptr;
 560   png_infop info_ptr;
 561   int irow;
 562   unsigned char *rowp;
 563
 564   rowp = (unsigned char *) sys_alloc(im->nx * nxcell,"sdf2d_to_img");
 565   if (rowp == NULL)
 566     return;
 567
 568   if ((fp = fopen (outfile, "wb")) == NULL)
 569      return;
 570
 571   png_ptr = png_create_write_struct (PNG_LIBPNG_VER_STRING, NULL, NULL, NULL);
 572   if (! png_ptr)
 573     return;
 574
 575   info_ptr = png_create_info_struct(png_ptr);
 576   if (! info_ptr)
 577     {
 578       png_destroy_write_struct(&png_ptr, (png_infopp) NULL);
 579       fclose(fp);
 580       return;
 581     }
 582
 583   if (setjmp(png_ptr->jmpbuf))
 584     {
 585       png_destroy_write_struct(&png_ptr, &info_ptr);
 586       fclose(fp);
 587       return;
 588     }
 589
 590   png_init_io(png_ptr, fp);
 591
 592   png_set_IHDR(png_ptr, info_ptr, im->nx * nxcell, im->ny * nycell, 8, PNG_COLOR_TYPE_GRAY, PNG_INTERLACE_NONE, PNG_COMPRESSION_TYPE_BASE, PNG_FILTER_TYPE_DEFAULT);
 593
 594   png_write_info(png_ptr, info_ptr);
 595   for (irow = im->ny - 1; irow >= 0; irow--)
 596     {
 597       png_bytep row_pointer;
 598       int icol, ir;
 599
 600       for (icol = 0; icol < im->nx; icol++)
 601         {
 602           double dens;
 603           int p;
 604           int pos = icol * nxcell;
 605           dens = (im->v[icol][irow] - densmin) / (densmax - densmin);
 606           if (dens < 0)
 607             dens = 0;
 608           else if (dens > 1)
 609             dens = 1;
 610           for (p = pos; p < pos + nxcell; p++)
 611             {
 612               rowp[p] = (unsigned int) (dens * 255.);
 613             }
 614         }
 615       row_pointer = rowp;
 616       for (ir = 0; ir < nycell; ir++)
 617         {
 618           png_write_rows (png_ptr, &row_pointer, 1);
 619         }
 620     }
 621   sys_free(rowp, "sdf2img");
 622
 623   png_write_end(png_ptr, info_ptr);
 624   png_destroy_write_struct(&png_ptr, &info_ptr);
 625
 626   fclose(fp);
 627 }
 628 #endif
 629
 630 #ifdef HAVE_GD
 631 #include "gd.h"
 632 #define N_GRAYSCALE 256
 633 #endif
 634
 635 void
 636 sdf2d_to_gif (IMAGE *im, char *outfile, int nxcell, int nycell, double densmin, double densmax)
 637 {
 638 #ifdef HAVE_GD
 639   gdImagePtr gif;
 640   FILE *out;
 641   int gs_indices[N_GRAYSCALE];
 642   int i;
 643   int irow;
 644   int lastrow;
 645   unsigned char *rowp;
 646
 647   rowp = (unsigned char *) sys_alloc(im->nx * nxcell,"sdf2d_to_img");
 648   if (rowp == NULL)
 649     return;
 650
 651   gif = gdImageCreate(im->nx * nxcell, im->ny * nycell);
 652   for (i = 0; i < N_GRAYSCALE; i++)
 653     gs_indices[i] = gdImageColorAllocate(gif, i, i, i);
 654
 655   lastrow = im->ny * nycell - 1;
 656   for (irow = 0; irow < im->ny; irow++)
 657     {
 658       int icol, ir;
 659       int rpos = irow * nycell;
 660       for (ir = rpos; ir < rpos + nycell; ir++)
 661         {
 662           for (icol = 0; icol < im->nx; icol++)
 663             {
 664               double dens;
 665               int ic;
 666               int cpos = icol * nxcell;
 667               dens = (im->v[icol][irow] - densmin) / (densmax - densmin);
 668               if (dens < 0)
 669                 dens = 0;
 670               else if (dens > 1)
 671                 dens = 1;
 672               for (ic = cpos; ic < cpos + nxcell; ic++)
 673                 {
 674                   rowp[ic] = (unsigned int) (dens * (double) (N_GRAYSCALE - 1));
 675                   gdImageSetPixel(gif, ic, lastrow - ir, gs_indices[rowp[ic]]);
 676                 }
 677             }
 678         }
 679     }
 680   sys_free(rowp, "sdf2img");
 681
 682   if ((out = fopen(outfile,"w")) == NULL)
 683     {
 684       fprintf(stderr,"Error opening output file %s for writing\n", outfile);
 685       exit(1);
 686     }
 687   gdImageGif(gif,out);
 688   fclose(out);
 689   gdImageDestroy(gif);
 690 #else
 691   fprintf(stderr, "This version does not support GIF");
 692 #endif
 693 }
 694