func writeImage(h *fits.Unit, name string) {
	n := len(h.Naxis)
	maxis := make([]int, n)
	img := image.NewGray16(image.Rect(0, 0, h.Naxis[0], h.Naxis[1]))
	prod := 1
	for k := 2; k < n; k++ {
		prod *= h.Naxis[k]
	}
	min, max := h.Stats()
	for i := 0; i < prod; i++ {
		l := i
		s := name
		for k := 2; k < n; k++ {
			maxis[k] = l % h.Naxis[k]
			l = l / h.Naxis[k]
			s += fmt.Sprintf("-%d", maxis[k])
		}
		for x := 0; x < h.Naxis[0]; x++ {
			for y := 0; y < h.Naxis[1]; y++ {
				maxis[0] = x
				maxis[1] = y
				if !h.Blank(maxis...) {
					v := uint16((h.FloatAt(maxis...) - min) / (max - min) * 65535)  
					img.SetGray16(x, h.Naxis[1]-y, color.Gray16{v})
				} else {
					img.SetGray16(x, h.Naxis[1]-y, color.Gray16{0})  
				}
			}
		}
		g, _ := os.Create(s + ".png")
		defer g.Close()
		png.Encode(g, img)
	}
}
writeImage generates PNG file(s) for image type HDUs
If h contains a two-dimensional image, a single 16-bit normalized PNG file is generated
For higher dimensional units, NAXIS3xNAXIS4x... different images (each NAXIS1xNAXIS2 in size) are generated
For example, if h.Naxis=[512, 512, 2, 3] and name is "test_0", the following images are written:

test_0-0.0.png  contains pixels [0,0,0,0] to [511,511,0,0]
test_0-1.0.png  contains pixels [0,0,1,0] to [511,511,1,0]
test_0-0.1.png  contains pixels [0,0,0,1] to [511,511,0,1]
test_0-1.1.png  contains pixels [0,0,1,1] to [511,511,1,1]
test_0-0.2.png  contains pixels [0,0,0,2] to [511,511,0,2]
test_0-1.2.png  contains pixels [0,0,1,2] to [511,511,1,2]
package main
import (
	"bytes"	
	"fmt"
	"image"
	"image/color"
	"image/png"
	"io/ioutil"
	"log"
	"net/http"
	"net/url"
	"os"
	"path"
	"strings"    
    "fits"
)
func main() {
	var units []*fits.Unit
	var name string
	if len(os.Args) == 1 {
		fmt.Println("usage: extract filename|url")
		os.Exit(1)
	}
	if strings.HasPrefix(os.Args[1], "http://") {  
		url, err := url.Parse(os.Args[1])
		if err != nil {
			log.Fatal(err)
		}
		name = path.Base(url.Path)
		res, err := http.Get(os.Args[1])
		if err != nil {
			log.Fatal(err)
		}
		buf, _ := ioutil.ReadAll(res.Body)  
		res.Body.Close()
		units, err = fits.Open(bytes.NewReader(buf))
		if err != nil {
			log.Fatal(err)
		}
	} else {  
		name = path.Base(os.Args[1])
		reader, err := os.Open(os.Args[1])
		if err != nil {
			log.Fatal(err)
		}
		defer reader.Close()
		units, err = fits.Open(reader)
		if err != nil {
			log.Fatal(err)
		}
	}
	j := strings.LastIndex(name, ".")  
	if j != -1 {
		name = string(name[:j])
	}
	for i, h := range units {  
		fmt.Printf("******************** Header %d ********************\n", i)
		for key, value := range h.Keys {  
			fmt.Println(key, value)
		}
		out := fmt.Sprintf("%s_%d", name, i)
		if h.HasImage() {  
			writeImage(h, out)
		} else if h.HasTable() {  
			if len(h.Naxis) == 1 {  
				writeArray(h, out)
			} else {  
				writeTable(h, out)
			}
		} else {
			fmt.Println("Unsupported Header Data Unit")
		}
	}
}
func writeArray(h *fits.Unit, name string) {
	g, _ := os.Create(name + ".dat")
	defer g.Close()
	for i := 0; i < h.Naxis[0]; i++ {
		fmt.Fprintln(g, h.FloatAt(i))
	}
}