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