Permalink
Switch branches/tags
Nothing to show
Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
109 lines (93 sloc) 2.96 KB
module UCSF
#Package to read .ucsf binary file format (NMR) multidimensional spectra.
#Reference:
#https://www.cgl.ucsf.edu/home/sparky/manual/files.html#UCSFFormat
export read_ucsf,get_axis_ppm
immutable UCSF_Header
dims::Int32
num_comp::Int32
format_ver::Int32
end
immutable Axis_Header
atom::String
data_length::Int32
tile_length::Int32
spectrometer_freq::Float32
spectral_width::Float32
center::Float32
end
function read_main_header(io::IO)
t = true
seek(io,0)
data = Vector{UInt8}(180)
readbytes!(io,data,180)
@assert join(Char.(data[1:8])) == "UCSF NMR" "Wrong file format! (not UCSF)"
UCSF_Header( Int(data[11]), Int(data[12]), Int(data[14]) )
end
function read_axis_header(io::IO,h::UCSF_Header)
seek(io,180)
header = Vector{Axis_Header}(h.dims)
data = Vector{UInt8}(128)
for i = 1:h.dims
readbytes!(io,data,128)
#Reversal needed because of the big-endian thing
header[i] = Axis_Header( rstrip(join(Char.(data[1:6])),'\0'),
reinterpret(Int32,reverse(data[9:12]))[1],
reinterpret(Int32,reverse(data[17:20]))[1],
reinterpret(Float32,reverse(data[21:24]))[1],
reinterpret(Float32,reverse(data[25:28]))[1],
reinterpret(Float32,reverse(data[29:32]))[1])
end
header
end
move(u::UnitRange,offset::Integer) = (u.start+offset):(u.stop+offset)
function read_data(io::IO,h::Vector{Axis_Header})
dims = length(h)
seek(io,180+128*dims)
sl = [ Int64(i.data_length) for i in h ]
data = Vector{Float32}(prod(sl))
b = Vector{UInt8}(4)
for (k,i) = enumerate(1:4:(4*prod(sl)))
readbytes!(io,b,4)
data[k] = reinterpret(Float32,reverse(b))[1]
end
data
end
function arrange_tiles(data::Vector,h::Vector{UCSF.Axis_Header})
dims = length(h)
sl = reverse([ Int64(i.data_length) for i in h ])
tl = reverse([ Int64(i.tile_length) for i in h ])
num_tiles = [ ceil(Int64,sl[i]/tl[i]) for i in 1:dims ]
spectrum = Array{Float32,dims}(sl...)
for i = 1:prod(num_tiles)
tile_ind = ind2sub((num_tiles...),i)
data_ind = move(1:prod(tl),(i-1)*prod(tl))
spec_ind = [ move(1:tl[j],(tile_ind[j]-1)*tl[j]) for j in 1:dims ]
tile = reshape(data[data_ind],tl...)
spectrum[spec_ind...] = tile
end
spectrum
end
"""
read_ucsf(file::String)
Returns header objects and spectrum `Array` from a ucsf file format.
"""
function read_ucsf(file::String)
f = open(file)
u = read_main_header(f)
a = read_axis_header(f,u)
data = read_data(f,a)
spectrum = arrange_tiles(data,a)
(u,a,spectrum)
end
"""
`get_axis_ppm(a::Axis_Header)`
Returns atom and its frequency range from a `Axis_Header` object.
"""
function get_axis_ppm(a::Axis_Header)
δ = linspace(a.center + a.spectral_width/a.spectrometer_freq/2,
a.center - a.spectral_width/a.spectrometer_freq/2,
a.data_length)
(a.atom,δ)
end
end # module