https://github.com/QBioLab/CountmRNA.jl
Tip revision: 27fa41066bf09a3c1ca00fd333f72755e481b16f authored by H.F on 13 April 2021, 10:57:05 UTC
upload running script
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