File: img.i

package info (click to toggle)
yorick-yutils 1.3.0-2
  • links: PTS
  • area: main
  • in suites: lenny
  • size: 604 kB
  • ctags: 17
  • sloc: makefile: 105; python: 12
file content (902 lines) | stat: -rw-r--r-- 32,388 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
/*
 * img.i --
 *
 *	Routines for dealing with images (i.e. 2D arrays) in Yorick.
 *
 *-----------------------------------------------------------------------------
 *
 * Copyright (c) 2000-2004, Eric Thibaut (Centre de Recherches
 *	Astronomiques de Lyon).
 *
 * This program is free software; you can redistribute it and/or  modify it
 * under the terms of the GNU General Public License  as  published  by the
 * Free Software Foundation; either version 2 of the License,  or  (at your
 * option) any later version.
 *
 * This program is distributed in the hope  that  it  will  be  useful, but
 * WITHOUT  ANY   WARRANTY;   without   even   the   implied   warranty  of
 * MERCHANTABILITY or  FITNESS  FOR  A  PARTICULAR  PURPOSE.   See  the GNU
 * General Public License for more details (to receive a  copy  of  the GNU
 * General Public License, write to the Free Software Foundation, Inc., 675
 * Mass Ave, Cambridge, MA 02139, USA).
 *
 * Routines:
 *	img_dims - returns dimension [WIDTH,HEIGHT] of an image
 *	img_plot - plot an image
 *	img_cbar - add a color to an image plot
 *	img_interpolate - bi-linear interpolation of an image
 *	img_extract_parallelogram_as_rectangle - (as its name says)
 *	img_max - get coordinates of maximum in a 2-D array
 *	img_fft_centered_at_max - recenter an image at its maximum according
 *		to fft indexing
 *	img_pad - expand an image
 *	img_paste - copy an image into another one at a given location
 *	img_photometry - measure integrated intensity inside a circular region
 *	img_flt_max - filter an image
 *	img_flt_flac - flipped-local-auto-convolution of an image
 *	ing_get_type - get image file type
 *	img_read, img_write - read/write PNM/JPEG/PNG/TIFF/FITS/GIF image
 *	img_tmpnam - get name of temporary file
 *
 * History:
 *	$Id: img.i,v 1.1.1.1 2007/12/11 23:55:12 frigaut Exp $
 *	$Log: img.i,v $
 *	Revision 1.1.1.1  2007/12/11 23:55:12  frigaut
 *	Initial Import - yorick-yutils
 *	
 *	Revision 1.4  2004/10/14 09:54:53  eric
 *	 - "img_protect_path" and "img_expand_path" removed and replaced by
 *	   their counterparts, "expand_file_name" and "expand_file_name",
 *	   in "utils.i" (which was already required by "img.i").
 *
 *	Revision 1.3  2004/10/13 10:18:32  eric
 *	 - New functions (img_read, img_write, img_get_type,
 *	   img_protect_path, img_expand_path, img_tmpnam) to
 *	   implement reading/writing of various image file formats
 *	   (PNM/PBM/PGM/PPM, JPEG, PNG, TIFF, FITS, and GIF).
 *
 *	Revision 1.2  2004/09/17 11:21:18  eric
 *	 - removed jped_read and jpeg_write (backup in jpeg_img.i) which
 *	   are provided by jpeg.i plugin in Yorick from version >= 1.6
 *
 *	Revision 1.1  2004/07/17 13:37:40  eric
 *	Initial revision
 *
 *-----------------------------------------------------------------------------
 */

require, "utils.i";

func img_dims(img)
/* DOCUMENT img_dims(img)
     Returns dimensions of image IMG: [WIDTH,HEIGHT].  IMG must be a 2
     dimensional array.

   SEE ALSO dimsof. */
{
  if (! is_array(img) || (dims = dimsof(img))(1) != 2)
    error, "expecting 2D array";
  return dims(2:3);
}

/*---------------------------------------------------------------------------*/

func img_plot(img, first=, scale=, cmin=, cmax=, top=)
/* DOCUMENT img_plot, img;
     Plot image (2D array) IMG using pli.  Keyword FIRST can be used to set
     coordinates  of  center  of  first  (lower  left)  pixel,  default  is
     FIRST=1.0 (i.e. same coordinates as Yorick's indexing rules).  Keyword
     SCALE  can  be used  to  set the  step  size  between adjacent  pixels
     (default is SCALE=1.0).   Keywords FIRST and/or SCALE may  have one or
     two values depending whether or not both axis have same value.

   SEE ALSO pli, img_dims. */
{
  dims = img_dims(img);
  if (is_void(scale)) scale = 1.0;
  if (is_void(first)) first = 1.0;
  ll = first - scale/2.0; // lower-left corner coordinates
  ur = ll + scale*dims;   // upper-right corner coordinates
  pli, img, ll(1), ll(0), ur(1), ur(0), top=top, cmin=cmin, cmax=cmax;
}

func img_cbar(img, cmin=, cmax=, top=,
              levs=, labs=, nticks=, vert=, adjust=, vport=,
              format=, color=, width=, height=, ticklen=, font=)
/* DOCUMENT img_cbar, img;
         or img_cbar, cmin=cmin, cmax=cmax;
         
     Draw a color bar below the  current coordinate system.  If LEVS is not
     specified uses plfc_levs (set by previous call to plfc).  If COLORS is
     specified, it should have one  more value than LEVS, otherwise equally
     spaced colors are chosen, or  plfc_colors if plfc_levs was used.  With
     the VERT=1  keyword the color bar  appears to the left  of the current
     coordinate  system (vert=0  is default).   By default,  color_bar will
     attempt  to  label some  of  the  color  interfaces.  With  the  LABS
     keyword,  you can  force the  labelling algorithm  as  follows: LABS=0
     supresses all  labels, LABS=n forces  a label at every  n-th interface,
     LABS=[i,n]  forces  a  label  at  every n-th  interface  starting  from
     interface i (0<=i<=numberof(LEVS)).

     You    can    specify   the    viewport    coordinates   by    keyword
     VPORT=[xmin,xmax,ymin,ymax]; by default the  colorbar is drawn next to
     the current viewport.  You can use the ADJUST keyword  to move the bar
     closer to (adjust<0) or further from (adjust>0) the viewport.

     You  can specify  the string  format  for labels  with keyword  FORMAT
     (default "%g"), the font  type with keyword FONT (default "helvetica")
     and the font height with keyword HEIGHT (default 14 points).

     Keyword COLOR  can be  used to  specify the color  of the  labels, the
     ticks and the frame of the colorbar.  Default is foreground color.

     Keyword WIDTH can be  used to set the width of the  lines used to draw
     the frame and the ticks of the colorbar.

     Keyword TICKLEN  can be used to set  the lenght (in NDC  units) of the
     ticks.  Default is 0.005 NDC.
     
   SEE ALSO: plfc. */
{
  nil = string(0);
  if (is_void(cmin)) {
    if (is_void(img)) error, "keyword CMIN must be given";
    cmin = min(img);
  }
  if (is_void(cmax)) {
    if (is_void(img)) error, "keyword CMAX must be given";
    cmax = max(img);
  }
  cmin = double(cmin); /* make sure CMIN is double */
  cmax = double(cmax); /* make sure CMAX is double */

  if (is_void(top)) {
    /* get indices in colormap */
    crange = long(bytscl([0,1],cmin=0,cmax=1));
    ncolors = crange(2) - crange(1) + 1;
  } else {
    ncolors = top + 1;
  }

  if (structof(labs) == string) {
    nticks = numberof(labs);
    if (is_void(levs)) {
      levs = array(double, nticks);
      if (sread(labs, levs) != nticks)
        error, "cannot convert tick labels into values from LABS, use keyword LEVS";
    } else if (numberof(levs) != nticks) {
      error, "LABS and LEVS must have the same number of elements";
    }
  } else if (is_void(labs)) {
    if (is_void(levs)) {
      if (is_void(nticks)) nticks = 11;
      /* avoid rounding errors in span() */
      levs = cmin + ((cmax - cmin)/(nticks - 1))*indgen(0:nticks-1);
    } else {
      nticks = numberof(levs);
    }
    if (nticks) labs= swrite(format=(is_void(format)?"%g":format), levs);
  } else {
    error, "LABS must be nil or an array of strings";
  }
  
  if (is_void(font)) font= "helvetica";
  if (is_void(vport)) vport = viewport();    
  if (is_void(adjust)) adjust = 0.0;
  if (is_void(ticklen)) ticklen = 0.005;
  oldsys = plsys(0); /* _after_ calling viewport() */
  if (vert) {
    x0 = vport(2) + adjust + 0.022;
    x1 = x0 + 0.020;
    y0 = vport(3);
    y1 = vport(4);
    if (nticks) {
      x = x1(-::nticks);
      y = y0 + (y1 - y0)/(cmax - cmin)*(levs - cmin);
      if (ticklen) pldj, x, y, x + ticklen, y, legend=nil, color=color, width=width;
      else ticklen = 0.005;
      x += 2*ticklen;
      justify="LH";
    }
  } else {
    x0 = vport(1);
    x1 = vport(2);
    y0 = vport(3) - adjust - 0.045;
    y1 = y0 - 0.020;
    if (nticks) {
      x = x0 + (x1 - x0)/(cmax - cmin)*(levs - cmin);
      y = y1(-::nticks);
      if (ticklen) pldj, x, y, x, y - ticklen, legend=nil, color=color, width=width;
      else ticklen = 0.005;
      y -= 2*ticklen;
      justify="CT";
    }
  }
  for (i=1 ; i<=nticks ; ++i) {
    plt, labs(i), x(i), y(i), justify=justify, height=height, font=font, color=color;
  }
  colors = char(indgen(0:ncolors-1));
  colors = (vert ? colors(-,) : colors(,-));
  /* FIXME: there is a bug in Yorick, I have to make the colorbar at least 2 pixel wide to
     avoid division by zero in pli... */
  pli, colors, x0, y0, x1, y1;
  
  if (width) {
    plg, [y0,y0,y1,y1,y0], [x0,x1,x1,x0,x0], closed=0,
      width=width, color=color, legend=nil, marks=0, type=1;
  }
  plsys, oldsys;
}

/*---------------------------------------------------------------------------*/
/* BI-LINEAR INTERPOLATION */

func img_interpolate(z, x, y)
/* DOCUMENT  img_interpolate(img, x, y)
     Returns image IMG interpolated (by bi-linear interpolation) at
     positions (X,Y).  The coordinates X and Y must be conformable and the
     result has dimensopn list dimsof(X,Y).  Note that coordinates run like
     Yorick indices, for instance (1,1) is the location of the lower-left
     pixel in the image.
     
   SEE ALSO: interp. */
{
  if (! is_array(z) || (dims = dimsof(z))(1) != 2) {
    error, "expecting 2-D array";
  }
  if ((width = dims(2)) < 2 || (height = dims(3)) < 2) {
    error, "array to interpolate must have at least 2 elements in every dimension";
  }
  x -= (i = ceil(x));
  i = long(i);
  if (min(i) < 1) {
    k = where(i < 1);
    i(k) = 1;
    x(k) = 0.0;
  }
  if (max(i) >= width) {
    k = where(i >= width);
    i(k) = width - 1;
    x(k) = 1.0;
  }
  y -= (j = ceil(y));
  j = long(j);
  if (min(j) < 1) {
    k = where(j < 1);
    j(k) = 1;
    y(k) = 0.0;
  }
  if (max(j) >= height) {
    k = where(j >= height);
    j(k) = height - 1;
    y(k) = 1.0;
  }
  i += (j - 1)*width; /* univariate index */
  return ((1.0 - x)*((1.0 - y)*z(i) + y*z(i + width)) +
          x*((1.0 - y)*z(i + 1) + y*z(i + (1 + width))));
}

func img_extract_parallelogram_as_rectangle(img, x1, y1, x2, y2, x3, y3, w, h)
/* DOCUMENT img_extract_parallelogram_as_rectangle(img, x1, y1, x2, y2, x3, y3,
                                                   width, height)
      Returns a WIDTH-by-HEIGHT rectangle obtained by bi-linear
      interpolation of a parallelogram region from image IMG.  The
      parallelogram is specified by the coordinates of 3 of its corner:
      (X1,Y1) is the upper-left corner, (X2,Y2) is the lower-left corner
      and (X3,Y3) is the lower-right corner.  Note that coordinates run
      like Yorick indices: (1,1) is the location of the lower-left pixel in
      the image.
      
   SEE ALSO: img_interpolate, LUsolve. */
{
  a = LUsolve([[1,0,1,0,1,0],
               [1,0,1,0,w,0],
               [h,0,1,0,1,0],
               [0,1,0,1,0,1],
               [0,1,0,1,0,w],
               [0,h,0,1,0,1]],
              [x1, y1, x2, y2, x3, y3]);
  u = double(indgen(w));
  v = double(indgen(h)(-,));
  return img_interpolate(img,
                         a(1) + a(2)*u + a(3)*v,
                         a(4) + a(5)*u + a(6)*v);
}

/*---------------------------------------------------------------------------*/

func img_max(img)
/* DOCUMENT img_max(img)
     Returns coordinates of first image maximum.

   SEE ALSO: img_dims, img_flt_max. */
{
  width = img_dims(img)(1);
  i = img(*)(mxx) - 1;
  return [i%width + 1, i/width + 1];
}

func img_fft_centered_at_max(img)
/* DOCUMENT img_fft_centered_at_max(img)
     Returns image IMG rolled so that its maximum is at coordinates (1,1).

   SEE ALSO: img_dims, img_max, img_flt_max, roll. */
{
  type = structof(img);
  img = roll(img, 1 - img_max(img));
  return type == structof(img) ? img : type(img);
}

/*---------------------------------------------------------------------------*/

func img_pad(img, .., bg=, just=)
/* DOCUMENT img_pad(img)
       -or- img_pad(img, dimlist)
       -or- img_pad(img, width, height)
     Pad an image to another size, one of:
       - the smallest square 2D array that contains the image
       - an 2D array with dimension list DIMLIST
       - a WIDTH-by-HEIGHT array

     The padding value can be specified with keyword BG (for "background"),
     the default is 0.
     
     The type of the result depends on the types of IMG and BG.

     The justification is set by keyword JUST:
       JUST =  0 or nil -> lower-left (the default)
               1        -> center
              -1        -> at corners to preserve FFT indexing

   SEE ALSO: img_paste, img_dims. */
{
  /* Get new dimension list. */
  local new1, new2;
  dims = img_dims(img);
  old1 = dims(1);
  old2 = dims(2);
  nargs = 0;
  while (more_args()) {
    arg = next_arg();
    ++nargs;
    if ((s = structof(arg)) != long && s != int && s != short && s != char)
      error, "invalid non-integer dimension list";
    if ((arg_dims = dimsof(arg))(1) == 0) {
      /* got a scalar argument */
      if (nargs == 1)      new1 = long(arg);
      else if (nargs == 2) new2 = long(arg);
      else error, "too many new dimensions";
    } else if (nargs==1 && arg_dims(1)==1 && arg_dims(2)==3 && arg(1)==2) {
      /* got a valid dimlist */
      new1 = arg(2);
      new2 = arg(3);
    } else error, "bad dimension list";
  }
  if (nargs==0) new1 = new2 = max(dims);
  else if (nargs == 1 && is_void(new2)) new2 = new1;
  if (new1<old1 || new2<old2) error, "cannot reduce size";

  /* produce new image */
  if (is_void(bg)) {
    new = array(structof(img), new1, new2);
  } else if (is_array(bg) && dimsof(bg)(1) == 0) {
    new = array(structof(img(1) + bg)(bg), new1, new2);
  } else error, "bad value for keyword BG";
  if (is_void(just) || just==0) {
    /* image will not be centered */
    new(1:old1, 1:old2) = img;
  } else if (just==1) {
    /* image will be centered */
    i1 = (new1-old1)/2;
    i2 = (new2-old2)/2;
    new(i1+1:old1+i1, i2+1:old2+i2) = img;
  } else if (just==-1) {
    /* preserve FFT indexing */
    h1 = old1/2;
    h2 = old2/2;
    if (h1 && h2) new(new1-h1+1:new1, new2-h2+1:new2) = img(1:h1, 1:h2);
    if (h1) new(new1-h1+1:new1, 1:old2-h2) = img(1:h1, h2+1:old2);
    if (h2) new(1:old1-h1, new2-h2+1:new2) = img(h1+1:old1, 1:h2);
    new(1:old1-h1, 1:old2-h2) = img(h1+1:old1, h2+1:old2);
  } else error, "bad value for keyword JUST";
  
  return new;
}

/*---------------------------------------------------------------------------*/

func img_paste(dst, x, y, src)
/* DOCUMENT img_paste(dst, x, y, src)
     Paste image (2D array) SRC  at location (X,Y) in DST.  The coordinates
     (X,Y) are the  integer indices of the location in  DST where the lower
     left pixel  of SRC will be pasted.   SRC is possibly clipped  so as to
     fit  into DST  (e.g.  X  and/or Y  may be  less or  equal  zero).  The
     operation is  done "in-place"  and the result  is returned.   The data
     type of DST is unchanged.
     
  SEE ALSO: img_pad, img_dims. */
{
  dst_dims = img_dims(dst);
  src_dims = img_dims(src);
  if (! is_integer_scalar(x)) error, "X must be an integer scalar";
  if (! is_integer_scalar(y)) error, "Y must be an integer scalar";
  x1 = (x0 = long(x)) + src_dims(1) - 1;
  y1 = (y0 = long(y)) + src_dims(2) - 1;

  /* insure that SRC will not be completely outside DST */
  if (x0<=dst_dims(1) && x1>=1 && y0<=dst_dims(2) && y1>=1) {
      /* maybe clip SRC along 1st dimension */
      if ((clip= x0<1)) {
        s0 = 2 - x0;
        x0 = 1;
      } else s0 = 1;
      if (x1>dst_dims(1)) {
        x1 = dst_dims(1);
        clip = 1n;
      }
      if (clip) src = src(s0:s0+x1-x0,);
      
      /* maybe clip SRC along 2nd dimension */
      if ((clip = y0<1)) {
        s0 = 2 - y0;
        y0 = 1;
      } else s0 = 1;
      if (y1>dst_dims(2)) {
        y1 = dst_dims(2);
        clip = 1n;
      }
      if (clip) src = src(,s0:s0+y1-y0);

      /* paste SRC in-place */
      dst(x0:x1, y0:y1) = src;
  }
  return dst;
}



/*---------------------------------------------------------------------------*/
/* MEASUREMENT */

/*
 * Circle of radius R and center (X0,Y0) rounded to nearest pixel:
 *
 *         floor(sqrt((X - X0)^2 + (Y - Y0)^2) + 0.5) = R
 *   <==>  R <= sqrt((X - X0)^2 + (Y - Y0)^2) + 0.5 < R + 1
 *   <==>  R - 0.5 <= sqrt((X - X0)^2 + (Y - Y0)^2) < R + 0.5
 *
 * Disk of radius R and center (X0,Y0) rounded to nearest pixel:
 *
 *         floor(sqrt((X - X0)^2 + (Y - Y0)^2) + 0.5) <= R
 *   <==>  sqrt((X - X0)^2 + (Y - Y0)^2) + 0.5 < R + 1
 *   <==>  sqrt((X - X0)^2 + (Y - Y0)^2) <  R + 0.5
 *   <==>       (X - X0)^2 + (Y - Y0)^2  < (R + 0.5)^2
 *
 * Corresponding rectangular region of interest:
 *
 *  Y = Y0  ==>  abs(X - X0) < R + 0.5
 *              -(R + 0.5) < X - X0 < R + 0.5
 *              X0 - (R + 0.5) < X < X0 + (R + 0.5)
 *
 */

func img_photometry(img, x, y, r, bg=)
/* DOCUMENT img_photometry(img, x, y, r)
     Returns  sum of pixel  values in  image IMG  whithin circular  area of
     radius R and centered at (X,Y).  Arrays X, Y and R are in pixel units,
     they may have any geometry  but must be conformable.  Coordinates have
     the  same  origin  as  array  indices:  the first  pixel  in  IMG  has
     coordinates (1,1).
     
   SEE ALSO dimsof, avg. */
{
  /* Get dimensions of image and make X, Y and R arrays conformable. */
  if (! is_array(img) || (dims = dimsof(img))(1) != 2)
    error, "expecting 2D array";
  nx = dims(2);
  ny = dims(3);
  if (is_void((dims = dimsof(x, y, r)))) error, "X, Y and R not conformable";
  zero = array(double, dims);
  if (structof((x += zero)) != double ||
      structof((y += zero)) != double ||
      structof((r += zero)) != double) error, "bad data type";
  zero = [];

  rp = r + 0.5;
  i1 = max(1, long(ceil(x - rp)));
  i2 = min(nx, long(floor(x + rp)));
  j1 = max(1, long(ceil(y - rp)));
  j2 = min(ny, long(floor(y + rp)));
  rp *= rp;
  result = array((structof(img) == complex ? complex : double), dims);
  n = numberof(x);
  for (k=1 ; k<=n ; ++k) {
    if ((i1k = i1(k)) <= (i2k = i2(k)) &&
        (j1k = j1(k)) <= (j2k = j2(k))) {
      irng = i1k:i2k;
      jrng = j1k:j2k;
      u = indgen(irng) - x(k);
      v = indgen(jrng) - y(k);
      if (is_array((i = where(u*u + (v*v)(-,) < rp(k))))) {
        result(k) = avg(img(irng, jrng)(i));
      }
    }
  }
  if (! is_void(bg)) result -= bg;
  PI = 3.14159265358979323848;
  return (2*PI*r*r)*result;
}

/*---------------------------------------------------------------------------*/
/* FILTERS */

func img_flt_max(img, width, uniq=)
/* DOCUMENT img_flt_max(img, width)
     Return indices of pixels in image IMG that have the maximum value in a
     local  WIDTH-by-WIDTH box  (WIDTH must  be odd).   If keyword  UNIQ is
     true, a uniq maximum is  selected in every WIDTH-by-WIDTH box, in this
     case the intensity of the maxima will be in ascending order.

   SEE ALSO: img_flt_flac. */
{
  if (width%2 != 1) error, "width must be odd";
  w = (width-1)/2;
  if (w<=0) {
    if (w<0) error, "width must be non-negative";
    return img;
  }
  dims = img_dims(img);
  dim1 = dims(1);
  dim2 = dims(2);
  
  /* Filter image to replace pixel value by local maximum value. */
  if (w+1 >= dim1) {
    tmp = img(max,)(-:1:dim1,);
  } else {
    tmp = img;
    for (i=1 ; i<=dim1 ; ++i) tmp(i,) = img(max:max(i-w,1):min(i+w,dim1),);
  }
  if (w+1 >= dim2) {
    lmx = tmp(,max)(,-:1:dim2);
  } else {
    lmx = tmp;
    for (i=1 ; i<=dim2 ; ++i) lmx(,i) = tmp(,max:max(i-w,1):min(i+w,dim2));
  }
  tmp = [];
  
  /* Get indices of local maxima. */
  i = where(img >= lmx);
  if (! uniq || ! is_array(i)) return i;
  lmx = [];

  /* Sort maxima (because we already know that close maxima must have the same
     value) and remove multiples. */
  z = img(i);
  s = sort(z);
  i = i(s);
  x = 1 + (i-1)%dim1;
  y = 1 + (i-1)/dim1;
  z = z(s);
  if (numberof(z) > 1) {
    j = 1;
    k = where((z(dif)!=0) + (abs(x(dif)) > w) + (abs(y(dif)) > w));
    if (is_array(k)) grow, j, 1+k;
    i = i(j);
  }
  return i;
}

func img_flt_flac(img, width)
/* DOCUMENT img_flt_flac(img, width)
     Compute ``flipped-local-auto-convolution'' (sic!)   of image IMG.  The
     output  pixel  value  is  the   local  convolution  --  in  a  box  of
     WIDTH-by-WIDTH pixels  (WIDTH must  be odd) --  of the input  image by
     itself rotated  by 180 degrees.  This  is usefull to  locate spikes in
     IMG that may have different shapes but are all nearly symmetrical with
     respect to both axis (i.e.   each spike has its shape nearly unchanged
     by a  180 degrees  rotation).  This processing  is a kind  of adaptive
     filtering.   WIDTH should  be roughly  as large  as the  typical spike
     width.  The impact of input noise  in the result is reduced by using a
     larger WIDTH. But  WIDTH should remain smaller than  half the smallest
     separation between  spikes.  The  computation time is  proportional to
     WIDTH^2.

   SEE ALSO: img_flt_max. */
{
  if (width%2 != 1) error, "width must be odd";
  w = (width-1)/2;
  dims = img_dims(img);
  dim1 = dims(1);
  dim2 = dims(2);
  if (max(dim1,dim2)<width) return array(double, dim1, dim2);
  if (structof(img)==string || structof(img(1)+0.0)!=double)
    error, "bad data type for IMG";
  if (structof(img)!=double) img = double(img);
  if (w==0) return img*img;

  // zero offset
  s = img(w+1:-w, w+1:-w);
  s = s*s/2.0;

  // non-zero offsets (use symmetry for speed up by a factor of 2)
  for (w2=0 ; w2<=w ; ++w2) {
    for (w1=(w2>0?-w:1) ; w1<=w ; ++w1) {
      s1 = img(w+1+w1 :  w1-w, w+1+w2 :  w2-w);
      s2 = img(w+1-w1 : -w1-w, w+1-w2 : -w2-w);
      s += s1*s2;
    }
  }

  out = array(double, dim1, dim2);
  out(w+1:-w, w+1:-w) = 2.0*s;
  return out;
}

/*---------------------------------------------------------------------------*/
/* READING/WRITING IMAGES */

local IMG_PNM, IMG_JPEG, IMG_PNG, IMG_TIFF, IMG_FITS, IMG_GIF;
func img_get_type(filename, type=, reading=)
/* DOCUMENT img_get_type(filename)
     Returns image type for file FILENAME.  If keyword READING is true,
     then FILENAME must exists and the image type is obtained from the file
     signature (the 4 first bytes in the file).  Otherwise, if keyword TYPE
     is true it must be one of the string: "pnm", "jpeg", "png", "tiff",
     "fits", or "gif". Finally if none of this keywords is set, the image
     type is guessed form the file extension.
     
     The returned value is one of:
       1 = IMG_PNM    portable anymap (PBM/PGM/PPM) image;
       2 = IMG_JPEG   JPEG image;
       3 = IMG_PNG,   portable network graphic image;
       4 = IMG_TIFF   TIFF image;
       5 = IMG_FITS   FITS (flexible image transport system) file;
       6 = IMG_GIF    GIF image;
     
   SEE ALSO: img_read, img_write. */
{
  if (reading) {
    /* Guess file type from 4 byte header, according to magic numbers:
     *   \377 \330      = JPEG
     *   \211 "PNG"     = PNG
     *   "GIF8"         = GIF
     *   "MM" \000 \052 = TIFF image data, big-endian
     *   "II" \052 \000 = TIFF image data, little-endian
     *   "P1" space     = ascii PBM (portable bitmap)
     *   "P2" space     = ascii PGM (portable graymap)
     *   "P3" space     = ascii PPM (portable pixmap)
     *   "P4" space     = raw PBM
     *   "P5" space     = raw PGM
     *   "P6" space     = raw PPM
     */
    magic = array(char, 4);
    n = _read(open(filename, "rb"), 0, magic);
    c1 = (n >= 1 ? magic(1) : char(0));
    c2 = (n >= 2 ? magic(2) : char(0));
    c3 = (n >= 3 ? magic(3) : char(0));
    c4 = (n >= 4 ? magic(4) : char(0));
    if (c1 == '\377' &&  c2 == '\330')
      return IMG_JPEG;
    if (c1 == '\211' &&  c2 == 'P' && c3 == 'N' && c4 == 'G')
      return IMG_PNG;
    if ((c1 == 'M' &&  c2 == 'M' && c3 == '\000' && c4 == '\052') ||
        (c1 == 'I' &&  c2 == 'I' && c3 == '\052' && c4 == '\000'))
      return IMG_TIFF;
    if (c1 == 'S' &&  c2 == 'I' && c3 == 'M' && c4 == 'P')
      return IMG_FITS;
    if (c1 == 'P' && c2 >= '1' && c2 <= '6' &&
        (c3 == '\n' || c3 == '\r' || c3 == ' ' || c3 == '\t'))
      return IMG_PNM;
    if (c1 == 'G' &&  c2 == 'I' && c3 == 'F' && c4 == '8')
      return IMG_GIF;
    error, "cannot guess image file format";
  }
  if (is_void(type)) {
    /* Guess file type from its extension. */
    ext = strpart(filename, -3:0);
    if (ext == ".png" || ext == ".PNG") return IMG_PNG;    
    if (ext == ".jpg" || ext == ".JPG") return IMG_JPEG;
    if (ext == ".tif" || ext == ".TIF") return IMG_TIFF;
    if (ext == ".fit" || ext == ".FIT" ||
        ext == ".fts" || ext == ".FTS") return IMG_FITS;
    if (ext == ".gif" || ext == ".GIF") return IMG_GIF;
    if (ext == ".pnm" || ext == ".PNM" ||
        ext == ".pbm" || ext == ".PBM" ||
        ext == ".pgm" || ext == ".PGM" ||
        ext == ".ppm" || ext == ".PPM") return IMG_PNM;
    ext = strpart(filename, -4:0);
    if (ext == ".jpeg" || ext == ".JPEG") return IMG_JPEG;
    if (ext == ".tiff" || ext == ".TIFF") return IMG_TIFF;
    if (ext == ".fits" || ext == ".FITS") return IMG_FITS;
    error, "cannot guess image file format from file extension";
  }
  if (type == "pnm") return IMG_PNM;
  if (type == "jpeg") return IMG_JPEG;
  if (type == "png") return IMG_PNG;    
  if (type == "tiff") return IMG_TIFF;
  if (type == "fits") return IMG_FITS;
  if (type == "gif") return IMG_GIF;
  error, "bad image type name, expecting: pnm, jpeg, png, tiff, fits or gif";
}
IMG_PNM = 1;
IMG_JPEG = 2;
IMG_PNG = 3;
IMG_TIFF = 4;
IMG_FITS = 5;
IMG_GIF = 6;

func img_read(filename, tmp=)
/* DOCUMENT img_read(filename)
     Returns image  read from FILENAME.   Supported image formats  are: PNM
     (PBM/PBM/PPM), JPEG,  PNG, TIFF,  FITS and GIF.   For some  formats, a
     temporary PNM  image needs  to be created;  the name of  the temporary
     file can be specified with keyword TMP.
     
   SEE ALSO: system, pnm_read, fits_read, img_get_type,
             img_tmpnam, expand_file_name, protect_file_name. */
{
  filename = expand_file_name(filename);
  if (! open(filename,"r",1)) error, "input file does not exists";
  type = img_get_type(filename, reading=1);
  if (type == IMG_PNM) {
    convert = 0; /* no needs for conversion */
  } else if (type == IMG_JPEG) {
    convert = "jpegtopnm";
  } else if (type == IMG_PNG) {
    convert = "pngtopnm";
  } else if (type == IMG_TIFF) {
    convert = "tifftopnm";
  } else if (type == IMG_FITS) {
    /* TO DO: read compressed FITS file */
    if (! is_func(fits_read)) require, "fits.i";
    return fits_read(filename);
  } else if (type == IMG_GIF) {
    convert = "giftopnm";
  }

  /* read image as PNM file, possibly after conversion */
  if (! is_func(pnm_read)) require, "pnm.i";  
  if (convert) {
    if (is_void(tmp)) tmp = img_tmpnam(filename);
    if (catch(-1)) {
      remove, tmp;
      error, catch_message;
    }
    system, swrite(format="%s %s 2>/dev/null >%s", convert,
                   protect_file_name(filename), protect_file_name(tmp));
    img = pnm_read(tmp);
    remove, tmp;
  } else {
    img = pnm_read(filename);
  }
  return img;
}

func img_write(img, filename, type=, cmin=, cmax=, tmp=, noflip=,
               quality=, optimize=, grey=, gray=, progressive=, comment=,
               smooth=, eps=,
               compression=, interlace=, /* transparent=, gamma=, */
               bitpix=)
/* DOCUMENT img_write, img, filename;   
     Writes image IMG into file FILENAME as PNM (PBM/PBM/PPM), JPEG, PNG,
     TIFF or FITS image.  Except for a FITS file, if pixel type of IMG is
     not 'char', the pixel values are scaled to unsigned bytes (0-255) with
     bytscl function (which see).  The image type can be specified by the
     keyword TYPE otherwise it is automatically guessed from FILENAME
     extension (see img_get_type).

   KEYWORDS
     Keywords common to all format:
       TYPE - Output image type, one of: "jpeg", "pnm", "png", "tiff"
           or "fits".
       TMP - Name of the temporary file (format PBM, PGM, or PPM, see
           pnm_write) to creates.  Default is FILENAME~NUMBER where NUMBER is
           the smallest integer such that no file with the same name already
           exists (note that under race conditions the name of the default
           temporary file is not guaranteed to be unique).
       NOFLIP - See pnm_write -- not used for FITS format.
       CMIN, CMAX - Optional lower/upper cutoff (see bytscl) -- not used
           for FITS format.

     Keywords for JPEG images:
       EPS - If true, "jpeg2ps" is used to generate an encapsulated
           PostScript (level 2) image named FILENAME.eps (the JPEG image is
           not removed).
       QUALITY - JPEG quality (default 75).
       OPTIMIZE - Creates optimized JPEG image.
       GRAY/GREY - Creates grayscale JPEG image.
       PROGRESSIVE - Creates a progressive JPEG file.
       COMMENT - Text comment.
       SMOOTH=0-100 - Smooth the input image to eliminate dithering
           noise, 0 (the default) means no smoothing.

     Keywords for PNG images:
       COMPRESSION=1-9 - Level of compression (default is 6).
       INTERLACE - Creates an interlaced PNG file.
       
     Keywords for FITS images:
       BITPIX=n - Bits-per-pixel value.

     
   SEE ALSO: pnm_write, bytscl, system, img_get_type,
             img_tmpnam, expand_file_name, protect_file_name. */
{
  filename = expand_file_name(filename);
  type = img_get_type(filename, type=type);
  if (type == IMG_PNM) {
    convert = 0; /* no needs for conversion */
  } else if (type == IMG_JPEG) {
    if (! is_void(quality)) quatlity = 75;
    convert = swrite(format="pnmtojpeg --quality=%d", quality);
    if (optimize) convert += " --optimize";
    if (grey || gray) convert += " --grayscale";
    if (progressive) convert += " --progressive";
    if (! is_void(comment)) convert += swrite(format=" --comment=\"%s\"",
                                              comment);
    if (! is_void(smooth)) convert += swrite(format=" --smooth=%d", smooth);
  } else if (type == IMG_PNG) {
    convert = "pnmtopng";
    if (! is_void(compression))
      convert += swrite(format=" -compression=%d", compression);
    if (interlace)
      convert += " -interlace";
  } else if (type == IMG_TIFF) {
    convert = "pnmtotiff";
  } else if (type == IMG_FITS) {
    /* TO DO: read compressed FITS file */
    if (! is_func(fits_write)) require, "fits.i";
    return fits_write(filename, img, overwrite=1, bitpix=bitpix);
  } else if (type == IMG_GIF) {
    error, "GIF format not supported for writting";
  }

  /* byte scale image */
  if (structof(img) != char) img =  bytscl(img, top=255, cmin=cmin, cmax=cmax);

  /* write image as PNM file then convert to output type */
  if (! is_func(pnm_write)) require, "pnm.i";
  if (convert) {
    if (is_void(tmp)) tmp = img_tmpnam(filename);
    if (catch(-1)) {
      remove, tmp;
      error, catch_message;
    }
    pnm_write, img, tmp, noflip;
    filename_p = protect_file_name(filename);
    system, swrite(format="%s %s 2>/dev/null >%s",
                   convert, protect_file_name(tmp), filename_p);
    remove, tmp;
  } else {
    pnm_write, img, filename, noflip;
  }
  if (eps && type == IMG_JPEG) {
    /* make an encapsulated poscript file with JPEG image */
    if (is_void(filename_p)) filename_p = protect_file_name(filename);
    system, swrite(format="jpeg2ps %s >%s.eps", filename_p, filename_p);
  }
}

/*---------------------------------------------------------------------------*/
/* UTILITIES */

func img_tmpnam(name)
/* DOCUMENT img_tmpnam(name)
     Return a  string in the form: NAME~#  where # is an  integer chosen so
     that file  NAME~# does not exists.   Beware that there  is no absolute
     warranty that the returned name is not used elsewhere (for instance if
     two programs run at the same time and call the same function) but this
     is highly improbable.  In order to limit the  probabilty of such clash
     to occur, an empty file named NAME~# is created by the function.
     
   SEE ALSO: open. */
{
  fmt = "%s~%d";
  i = 0;
  do {
    tmp = swrite(format=fmt, name, ++i);
  } while (open(tmp, "r", 1));
  open, tmp, "w"; /* create empty file */
  return tmp;
}

/*---------------------------------------------------------------------------*/