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