Source

yt.go / read_data.go

Full commit
package main

import (
	// stdlib
	"fmt"
    "image"
    //"image/draw"
    //"image/color"
    "image/png"
    "math"
    "os"
    "time"

	// local
	"github.com/sbinet/go-hdf5/pkg/hdf5"
)

const (

)

type PixelArray struct {
    Px []float64
    Py []float64
    Pdx []float64
    Pdy []float64
    Vals []float64
}

type ImageBuffer struct {
    Nx int
    Ny int
    Buffer []float64
    XMin float64
    XMax float64
    YMin float64
    YMax float64
}

func (img *ImageBuffer) Fill(pix PixelArray) () {
    var px, py, pdx, pdy, val float64
    var lypx, rypx, lxpx, rxpx float64
    var overlap1, overlap2 float64
    var lc, lr, rc, rr int
    width := img.XMax - img.XMin
    px_dx := (float64(width) / float64(img.Nx))
    height := img.YMax - img.YMin
    px_dy := (float64(height) / float64(img.Ny))
    ipx_dx := 1.0 / px_dx
    ipx_dy := 1.0 / px_dy
    for pi := range pix.Px {
        px = pix.Px[pi]
        pdx = pix.Pdy[pi]
        if (px + pdx < img.XMin || px - pdx > img.XMax) {
            continue
        }
        py = pix.Py[pi]
        pdy = pix.Pdy[pi]
        if (py + pdy < img.YMin || py - pdy > img.YMax) {
            continue
        }
        val = pix.Vals[pi]
        lc = int(math.Max((px-pdx-img.XMin)*ipx_dx, 0.0))
        lr = int(math.Max((py-pdy-img.YMin)*ipx_dy, 0.0))
        rc = int(math.Min((px+pdx-img.XMin)*ipx_dx, float64(img.Nx)))
        rr = int(math.Min((py+pdy-img.YMin)*ipx_dy, float64(img.Ny)))
        for xi := lr; xi < rr; xi++ {
            lypx = px_dy * float64(xi) + img.YMin
            rypx = px_dy * float64(xi+1) + img.YMin
            overlap2 = (math.Min(rypx, py+pdy) - math.Max(lypx, py-pdy))*ipx_dy
            for yi := lc; yi < rc; yi++ {
                lxpx = px_dx * float64(yi) + img.XMin
                rxpx = px_dx * float64(yi+1) + img.XMin
                overlap1 = (math.Min(rxpx, px+pdx) - math.Max(lxpx, px-pdx))*ipx_dx
                if (overlap1 < 0.0 || overlap2 < 0.0) {
                    continue;
                }
                img.Buffer[xi * img.Ny + yi] = val
            }
        }
    }
    return
}

func NewImageBuffer(Nx, Ny int, XMin, XMax, YMin, YMax float64) (img *ImageBuffer) {
    img = new(ImageBuffer)
    img.Nx = Nx
    img.Ny = Ny
    img.XMin = XMin
    img.XMax = XMax
    img.YMin = YMin
    img.YMax = YMax
    img.Buffer = make([]float64, Nx*Ny)
    return
}

func FillElements(f *hdf5.File, field string)(val[] float64) {
    name := fmt.Sprintf("/Projections/0/%s", field)
    ds, err := f.OpenDataSet(name)
    if err != nil {
        fmt.Printf("** error opening %s ** %s\n", name, err)
        return
    }
    dss := ds.Space()
    nelements := dss.SimpleExtentNPoints()
    fmt.Printf("** %d elements in %s **\n", nelements, name)
    val = make([]float64, nelements)
    ds.Read(val, hdf5.T_NATIVE_DOUBLE)
    return
}

func (img *ImageBuffer) MapToImage(mi, ma float64) (out *image.RGBA){
    out = image.NewRGBA(image.Rect(0, 0, img.Nx, img.Ny))
    var v float64
    var cv uint8
    for pi := 0; pi < img.Ny*img.Nx ; pi++ {
        v = img.Buffer[pi]
        cv = uint8(255*(v - mi) / (ma - mi))
        for ci := 0; ci < 3; ci ++ {
            out.Pix[pi*4+ci] = cv
        }
        out.Pix[pi*4+3] = 255
    }
    return
}

func main() {
    f, err := hdf5.OpenFile("redshift0058.yt", hdf5.F_ACC_RDONLY)
    if err != nil {
        fmt.Printf("** error opening ** %s\n", err)
        return
    }
    var pixels PixelArray
    pixels.Px = FillElements(f, "px")
    pixels.Py = FillElements(f, "py")
    pixels.Pdx = FillElements(f, "pdx")
    pixels.Pdy = FillElements(f, "pdy")
    pixels.Vals = FillElements(f, "Density_Density")
    img := NewImageBuffer(1024, 1024, 0.0, 1.0, 0.0, 1.0)
    t1 := time.Now()
    N := 100
    for i := 0; i < N; i++ {
        img.Fill(pixels)
    }
    dt := time.Since(t1)
    spf := dt.Seconds() / float64(N)
    fps := 1.0 / spf
    fmt.Printf("Took %0.3e seconds each call for %0.3e fps\n", spf, fps)
    var mi, ma = 1e30, -1e30
    for i := range img.Buffer {
        mi = math.Min(mi, img.Buffer[i])
        ma = math.Max(ma, img.Buffer[i])
        img.Buffer[i] = math.Log10(img.Buffer[i])
    }
    fmt.Printf("Min value in Buffer: %0.3e\n", mi)
    fmt.Printf("Max value in Buffer: %0.3e\n", ma)
    out := img.MapToImage(math.Log10(mi), math.Log10(ma))
    file, _ := os.Create("test.png")
    png.Encode(file, out)
    file.Close()
}