https://github.com/QBioLab/CountmRNA.jl
Raw File
Tip revision: 27fa41066bf09a3c1ca00fd333f72755e481b16f authored by H.F on 13 April 2021, 10:57:05 UTC
upload running script
Tip revision: 27fa410
julia2ims.jl
using HDF5
using Dates 
#using Images
#https://portal.hdfgroup.org/display/HDF5/Introduction+to+HDF5
#https://imaris.oxinst.com/support/imaris-file-format
# https://github.com/imaris/ImarisWriter/blob/master/doc/Imaris5FileFormat.pdf
# https://docs.julialang.org/en/v1/stdlib/Dates/#Dates.format-Tuple{TimeType,AbstractString} 
# https://cpb-us-w2.wpmucdn.com/sites.wustl.edu/dist/f/1861/files/2019/02/02-imaris-image-properties-and-edit-2jj1avj.pdf
# https://github.com/tlambert03/imarispy/blob/master/imarispy/imaris.py

"""
Version Commit
0.1		hf, first verion
0.2     hf, update API to follow upstream HDF5.jl v0.15.4
TODO: genrate thumb and fix timestamps
      add function to append image
"""

"Save 4D UInt16 Array(x,y,z,t) to .ims file"
function  save2ims(array, fname::String="myfile.ims";
              subsamp=((1, 1, 1)),
              chunks=((128, 128, 5)),
              compression=2,
              thumbsize=256,
              dx=0.108, dz=0.5, dt_min=10)

    df = DateFormat("Y-mm-dd HH:MM:SS.sss")
    RecordingDate = DateTime(2020,5,7,14,19,22,0)

    nx, ny, nz, nt = size(array)
    nc = 1
    nr = 1     #nr = length(subsamp)

    # GROUPS of Imaris file
    GROUPS = [
        "DataSetInfo",
        "Thumbnail",
        "DataSetTimes",
        "DataSetInfo/Imaris",
        "DataSetInfo/Image",
        "DataSetInfo/TimeInfo"
    ]
    
    # Attributes of Imaris file
o   ATTRS = [
        ("/", "ImarisDataSet", "ImarisDataSet"),
        ("/", "ImarisVersion", "5.5.0"),
        ("/", "DataSetInfoDirectoryName", "DataSetInfo"),
        ("/", "ThumbnailDirectoryName", "Thumbnail"),
        ("/", "DataSetDirectoryName", "DataSet"),
        ("DataSetInfo/Imaris", "Version", "9.5"),
        ("DataSetInfo/Imaris", "ThumbnailMode", "thumbnailMIP"),
        ("DataSetInfo/Imaris", "ThumbnailSize", thumbsize),
        ("DataSetInfo/Image", "X", nx),
        ("DataSetInfo/Image", "Y", ny),
        ("DataSetInfo/Image", "Z", nz),
        ("DataSetInfo/Image", "NumberOfChannels", nc),
        ("DataSetInfo/Image", "Noc", nc),
        ("DataSetInfo/Image", "Unit", "um"),
        ("DataSetInfo/Image", "Description", "Hi"),
        ("DataSetInfo/Image", "MicroscopeModality", "Inverted Microscope"),
        ("DataSetInfo/Image", "RecordingDate", Dates.format(RecordingDate, df)),
        ("DataSetInfo/Image", "Name", "qblab"),
        ("DataSetInfo/Image", "ExtMin0", "0"),
        ("DataSetInfo/Image", "ExtMin1", "0"),
        ("DataSetInfo/Image", "ExtMin2", "0"),
        ("DataSetInfo/Image", "ExtMax0", nx * dx),
        ("DataSetInfo/Image", "ExtMax1", ny * dx),
        ("DataSetInfo/Image", "ExtMax2", nz * dz),
        ("DataSetInfo/Image", "LensPower", "40x"),
        ("DataSetInfo/TimeInfo", "DatasetTimePoints", nt),
        ("DataSetInfo/TimeInfo", "FileTimePoints", nt)
    ];

    COLORS = ("1 1 1", "1 0 1", "1 0 0", "0 0 1")
    for c in 0:nc-1
        grp = "DataSetInfo/Channel $c"
        push!(GROUPS, grp)
        push!(ATTRS, (grp, "ColorOpacity", 1))
        push!(ATTRS, (grp, "ColorMode", "BaseColor"))
        push!(ATTRS, (grp, "Color", COLORS[3]))
        push!(ATTRS, (grp, "GammaCorrection", 1))
        push!(ATTRS, (grp, "ColorRange", "0 255"))
        push!(ATTRS, (grp, "Name", "Channel $c"))
    end

    for t in 0:nt-1
        strr = Dates.format(RecordingDate + Dates.Minute(dt_min*t), df)
        #strr = "2020-05-07 {:02d}:{:02d}:{:02d}.000".format(h, m, s)
        push!(ATTRS, ("DataSetInfo/TimeInfo", "TimePoint$t", strr))
        #TODO: Imaris don't use these timestamps as time axis
    end

    h5open(fname, "w") do file
        for grp in GROUPS
            create_group(file, grp)
        end

        for info in ATTRS
            grp, key = info[1:2]
            value = "$(info[3])"
            h5a_write_S1(file[grp], key, value)
        end

        file["Thumbnail/Data"] = UInt8.(maximum(array[:,:, 1:20,1], dims=3).>>8)
        #file["Thumbnail/Data"] = UInt8.(maximum(rawview(array[:,:, 1:20,1]), dims=3).>>8)

        # add data
        for t in 0:nt-1
            for c in 0:nc-1
                data = array[:, :, :, t+1]
                for r in 0:nr-1
                    edges, count = build_histogram(data, 256)
                    grp = create_group(file, "/DataSet/ResolutionLevel $r/TimePoint $t/Channel $c/")
                    grp["Histogram"] = UInt64.(count[1:end])
                    h5a_write_S1(grp, "HistogramMin", "1000")
                    h5a_write_S1(grp, "HistogramMax", "$(edges[end])")
                    grp["Data", chunk= chunks , compress=compression] = data
                    #grp["Data", "chunk", (256,256,4)] = data
                    h5a_write_S1(grp, "ImageSizeX", "$(size(data)[1])")
                    h5a_write_S1(grp, "ImageSizeY", "$(size(data)[2])")
                    h5a_write_S1(grp, "ImageSizeZ", "$(size(data)[3])")
                end
            end
        end
    end
	fname
end


"Append image to exits imaris5 file"
#=
function ims_append(fname, img)
    # find orignal array size
    h5open(fname, "r+") do file
        # add data
        for t in 0:nt-1
            for c in 0:nc-1
                data = array[:, :, :, t+1]
                for r in 0:nr-1
                    edges, count = build_histogram(data, 256)
                    grp = create_group(file, "/DataSet/ResolutionLevel $r/TimePoint $t/Channel $c/")
                    grp["Histogram"] = UInt64.(count[1:end])
                    h5a_write_S1(grp, "HistogramMin", "1000")
                    h5a_write_S1(grp, "HistogramMax", "$(edges[end])")
                    grp["Data", chunk= chunks , compress=compression] = data
                    #grp["Data", "chunk", (256,256,4)] = data
                    h5a_write_S1(grp, "ImageSizeX", "$(size(data)[1])")
                    h5a_write_S1(grp, "ImageSizeY", "$(size(data)[2])")
                    h5a_write_S1(grp, "ImageSizeZ", "$(size(data)[3])")
                end
            end
        end
    end


end
=#

"Encode string with ASCII S1"
function h5a_write_S1(parent, key, value)
    # https://support.hdfgroup.org/HDF5/doc1.6/UG/11_Datatypes.html
    # https://support.hdfgroup.org/HDF5/doc/RM/PredefDTypes.html
    dspace = HDF5.dataspace((length(value),))
    try
        #attr = create_attribute(parent, key, datatype(HDF5.H5T_CSET_ASCII), dspace)
        #NOTE: datatype(obj) -> return the Datatype of obj
        #      HDF5.Datatype(id) -> return Datatype responding to id
        attr = create_attribute(parent, key, HDF5.Datatype(HDF5.H5T_C_S1), dspace)
        #write(attr, "$value")
        write_attribute(attr, HDF5.Datatype(HDF5.H5T_C_S1), "$value")
    finally
        close(dspace)
    end
    nothing
end

function str2int(str)
    local tmp =""
    for i in str
        tmp=tmp*i
    end
    parse(Int, tmp)
end

"Load imaris file to multi dimension image"
function loadims(imsname::String)
    local attr_t_len = h5readattr(imsname, "DataSetInfo/TimeInfo")["DatasetTimePoints"];
    local attr_x_len = h5readattr(imsname, "DataSetInfo/Image")["X"];
    local attr_y_len = h5readattr(imsname, "DataSetInfo/Image")["Y"];
    local attr_z_len = h5readattr(imsname, "DataSetInfo/Image")["Z"];
    local t_len  = str2int(attr_t_len)
    local x_depth = str2int(attr_x_len)
    local y_depth = str2int(attr_y_len)
    local z_depth = str2int(attr_z_len)
    #local img = zeros(Gray{N0f16}, x_depth, y_depth, z_depth, t_len);
    #local img = zeros(UInt16, x_depth, y_depth, z_depth, t_len);
    local img = Array{UInt16}(undef, x_depth, y_depth, z_depth, t_len);
    h5open(imsname, "r") do imsfile
        for i in 1:t_len
            data = read(imsfile, 
                        "DataSet/ResolutionLevel 0/TimePoint $(i-1)/Channel 0/Data")
            #img[:, :, :, i] = reinterpret(N0f16, data)
            img[:, :, :, i] = data
        end
    end
    return img
end

function loadims2(imsname::String)
    local attr_t_len = h5readattr(imsname, "DataSetInfo/TimeInfo")["DatasetTimePoints"];
    local t_len  = str2int(attr_t_len)
    local z_depth = 20
    local img = zeros(Gray{N0f16}, 512, 512, t_len*z_depth);
    h5open(imsname) do imsfile
        for i in 0:t_len-1
			data = imsfile["DataSet/ResolutionLevel 0/TimePoint $i/Channel 0/Data"]
            img[:, :, i*20+1:20*(i+1)] = reinterpret(N0f16,data[1:512,1:512,1:20]);
        end
    end
    img
end
back to top