#!/usr/bin/scsh -s

WARNING: this script can destroy files in the current working directory.
Run it in a directory which you've created specifically for this script.

This script is released into the public domain.
First publication:
http://www.bautforum.com/bad-astronomy-stories/68345-fighting-bad-astronomy-bad-astronomy.html

This script requires scsh and netpbm and gifsicle.  If you're running
Debian 4.0, "aptitude install scsh netpbm gifsicle" will install what
you need.  (You might also want to install scsh-doc.)

Start with cydonia1.tiff in the current working directory.

  cydonia1.tiff sources:
  ftp://nssdcftp.gsfc.nasa.gov/miscellaneous/planetary/mgs/cydonia1.tiff
    (sha1sum: 611771f7cec388c7640f6110eba8aa0012379ab1  cydonia1.tiff)
  http://windsor.gsfc.nasa.gov/miscellaneous/planetary/mgs/cydonia1.tiff
    (sha1sum: 611771f7cec388c7640f6110eba8aa0012379ab1  cydonia1.tiff)
  http://mpfwww.jpl.nasa.gov/mgs/target/CYD1/cydonia1.tif
    (sha1sum: 9d9204ad231666ab0a0c3ecca65cd64e4c28a7da  cydonia1.tif)

  cydonia1.tiff information:
  http://nssdc.gsfc.nasa.gov/planetary/mgs_cydonia.html
  http://mpfwww.jpl.nasa.gov/mgs/target/CYD1/

This script generates an animation which demonstrates the effects of the
following algorithm:  Choose a convolution kernel height, and call it h.
Create a matrix with h rows and 1 column.  Set the middle element of the
matrix to 1, and set all the other elements to 1/(1-h).  Convolve each
column in the image using this matrix.  Normalize contrast.

In the original "catbox convolution," h was 31.  The image had been scaled
down, to a swath width of 512 pixels, before the convolution was applied.

In the limit, as h approaches infinity, the algorithm above becomes
equivalent to the following even simpler algorithm:  Subtract, from each
pixel, the average of all pixels in the same column.  Normalize contrast.

!#

;;; This function creates a convolution kernel for pnmconvol.  This is not the
;;; catbox kernel, because unfortunately pnmconvol isn't general enough to use
;;; it directly.  Instead this is the kernel you would get if you started with
;;; the catbox kernel, and multiplied the middle row by 0, and multiplied every
;;; other row by -1.  Later we'll undo the multiplication by -1 (using
;;; pnminvert) and undo the multiplication by 0 (using ppmmix).

;;; to do: a simple scsh function could perform the catbox convolution
;;; directly, without needing pnmconvol and pnminvert and ppmmix.  And it could
;;; be faster, because it could take advantage of the fact that nearly all of
;;; the elements in the kernel have the same value.

(define (make-kernel height)
  (format #t "P2 1 ~d ~d~%" height (* (- height 1) 2))
  (let ((f (lambda ()
             (do ((i (/ (- height 1) 2) (- i 1)))
                 ((<= i 0))
               (format #t "~d~%" height)))))
    (f)
    (format #t "~d~%" (- height 1))
    (f)))


;;; to do: instead of redirecting output to these silly throw-away files,
;;; use pipelines.  Scsh and netpbm both work beautifully with pipelines.

(define (prelim)
  (&& (epf (tifftopnm cydonia1.tiff)
           (> a.ppm))
      (epf (ppmtopgm a.ppm)
           (> b.pgm))
      (epf (pnmdepth 65535 b.pgm)
           (> c.pgm))
      (epf (pnmscale -reduce 2 c.pgm)
           (> d.pgm))
      (epf (pnmflip -lr d.pgm)
           (> e.pgm))
      (epf (pamcut -top 1769 -height 1440 e.pgm)
           (> f.pgm))))


(define (per-height height)
  (let ((subtrahend (format #f "~a.pgm" height)))
    (and (cond ((number? height) ; not infinity
                (&& (epf (begin (make-kernel height))
                         (> k.pgm))
                    (epf (pnmconvol k.pgm f.pgm)
                         (> ,subtrahend))))
               (else ; infinity
                (&& (epf (pnmscale -width 512 -height 1 e.pgm)
                         (> l.pgm))
                    (epf (pnmscale -width 512 -height 1440 l.pgm)
                         (> ,subtrahend)))))
         (&& (epf (pnminvert ,subtrahend)
                  (> m.pgm))
             (epf (ppmmix "0.5" f.pgm m.pgm)
                  (> n.ppm))
             (epf (ppmtopgm n.ppm)
                  (> o.pgm))
             (epf (pamcut -top 480 -height 480 o.pgm)
                  (> p.pgm))
             (epf (pnmnorm -bpercent "0.005" -wpercent "0.005" p.pgm)
                  (> q.pgm))
             (epf (pnmdepth 255 q.pgm)
                  (> r.pgm))
             (epf (ppmlabel -colour white -x 239 -y 475 -text
                            ,(if (number? height) height "inf.") r.pgm)
                  (> s.ppm))
             (epf (ppmtogif s.ppm)
                  (> ,(format #f "~a.gif" height)))))))


(define (list-heights)
  (do ((n 0 (- n 1/2))
       (l '("inf")
          (cons (+ 1 (* 2 (inexact->exact (round (* #i480 (expt 2 n)))))) l)))
      ((< n -5) l)))


(define (gifsicle-args height)
  (list "--name" height (format #f "~a.gif" height)))


(let ((heights '(31 43 61 85 121 171 241 341 481 679 961 "inf")))
  (and (prelim)
       (every per-height heights)
       (run (gifsicle -l -V
                      -d165 ,@(gifsicle-args (car (reverse heights)))
                      -d33 ,@(apply append (map gifsicle-args
                                                (cdr (reverse (cdr heights)))))
                      -d165 ,@(gifsicle-args (car heights)))
            (> catbconv.gif))))