diff --git a/Project.toml b/Project.toml index 68c47c3..ba18270 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,8 @@ name = "MRCFile" uuid = "c18c01d3-0ab9-49c3-bd2d-8f2e64a2b7a5" authors = ["Seth Axen "] -version = "0.1.2" +version = "0.1.3" + [deps] CodecBzip2 = "523fee87-0ab8-5b00-afb7-3ecf72e48cfd" diff --git a/README.md b/README.md index 7e3cc04..cfd6592 100644 --- a/README.md +++ b/README.md @@ -1,4 +1,5 @@ # MRCFile.jl + ![Lifecycle](https://img.shields.io/badge/lifecycle-experimental-orange.svg) [![Build Status](https://github.com/sethaxen/MRCFile.jl/workflows/CI/badge.svg)](https://github.com/sethaxen/MRCFile.jl/actions?query=workflow%3ACI+branch%main) [![Coverage](https://codecov.io/gh/sethaxen/MRCFile.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/sethaxen/MRCFile.jl) @@ -6,12 +7,16 @@ [![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://sethaxen.github.io/MRCFile.jl/dev) [![Code Style: Blue](https://img.shields.io/badge/code%20style-blue-4495d1.svg)](https://github.com/invenia/BlueStyle) +## Description + MRCFile.jl implements the [MRC2014 format](https://www.ccpem.ac.uk/mrc_format/mrc2014.php) for storing image and volume data such as those produced by electron microscopy. It offers the ability to read, edit, and write MRC files, as well as utility functions for extracting useful information from the headers. The key type is `MRCData`, which contains the contents of the MRC file, accessible with `header` and `extendedheader`. `MRCData` is an `AbstractArray` whose elements are those of the data portion of the file and can be accessed or modified accordingly. +**Notice:** Since `1.1.3`, the `MRCData` is represent as row-major order to meet the requirement of [MRC2014 spec](https://www.ccpem.ac.uk/mrc_format/mrc2014.php). + ## Installation ```julia @@ -19,6 +24,14 @@ using Pkg Pkg.add("MRCFile") ``` +or install from the Julia package REPL, which can be accessed by pressing `]` from the Julia REPL: + +```julia +add MRCFile +``` + +See the [documentation](https://sethaxen.github.io/MRCFile.jl/stable/) for information on how to use MRCFile. + ## Example This example downloads a map of [TRPV1](https://www.emdataresource.org/EMD-5778) and animates slices taken through the map. @@ -50,7 +63,7 @@ gif(anim, "emd-$(emdid)_slices.gif", fps = 30) ![EMD-5778 slices](https://github.com/sethaxen/MRCFile.jl/blob/main/docs/src/assets/emd-5778_slices.gif) -# Reading a map as a memory-mapped array +## Reading a map as a memory-mapped array MRC files can be huge. It is convenient then to load their data as [memory-mapped arrays](https://docs.julialang.org/en/v1/stdlib/Mmap/). diff --git a/docs/make.jl b/docs/make.jl index c545d39..ff0b7bc 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -10,7 +10,12 @@ makedocs(; canonical="https://sethaxen.github.io/MRCFile.jl", assets=String[], ), - pages=["Home" => "index.md", "API" => "api.md"], + pages=[ + "Home" => "index.md", + "Documentation" => "documentation.md", + "Examples" => "examples.md", + "API" => "api.md", + ], ) deploydocs(; repo="github.com/sethaxen/MRCFile.jl", devbranch="main") diff --git a/docs/src/documentation.md b/docs/src/documentation.md new file mode 100644 index 0000000..8082f67 --- /dev/null +++ b/docs/src/documentation.md @@ -0,0 +1,102 @@ +# MRCFile documentation + +MRCFile.jl implements the [MRC2014 format](https://www.ccpem.ac.uk/mrc_format/mrc2014.php) for storing image and volume data such as those produced by electron microscopy. +It offers the ability to read, edit, and write MRC files, as well as utility functions for extracting useful information from the headers. + +The key type is `MRCData`, which contains the contents of the MRC file, accessible with `header` and `extendedheader`. +`MRCData` is an `AbstractArray` whose elements are those of the data portion of the file and can be accessed or modified accordingly. +The `data`, `header` and `extendedheader` are consistent with [mrcfile](https://github.com/ccpem/mrcfile). +Help on individual functions can be found in the API section or by using `?function name` from within Julia. + +## Opening MRC files + +MRC files can be opened using the `read()` or `read_mmap()` functions. These return an instance of the `MRCData` struct. + +```julia +using MRCFile + +dmap = read("/path/to/mrc/file.mrc", MRCData) +h = header(dmap) +exh = extendedheader(dmap) + +dmin, dmax = extrema(h) +drange = dmax - dmin +``` + +It is efficient to load the data as [memory-mapped arrays](https://docs.julialang.org/en/v1/stdlib/Mmap/). + +``` +# Open the file in memory-mapped mode +dmap = read_mmap("/path/to/mrc/file.mrc", MRCData) +``` + +## Handling compressed files + +All the functions above can also handle `.gz` or `.bz2` compressed files. + +```julia +using MRCFile + +dmap = read("/path/to/mrc/file.map.gz", MRCData) +# or +dmap = read("/path/to/mrc/file.map.bz2", MRCData) +``` + +## Write MRC file + +MRC files can be saved using the `write()` function. + +```julia +write("/path/to/mrc/file.mrc", dmap) +``` + +## Accessing the header + +The variables in header can be accessed using the following variable names, suppose `h=MRCHeader`: + +| Variable name | Access | Type | +| :-------------- | :------------ | :---------------- | +| NX | h.nx | Int32 | +| NY | h.ny | Int32 | +| NZ | h.nz | Int32 | +| MODE | h.mode | Int32 | +| NXSTART | h.nxstart | Int32 | +| NYSTART | h.nystart | Int32 | +| NZSTART | h.nzstart | Int32 | +| MX | h.mx | Int32 | +| MY | h.my | Int32 | +| MZ | h.mz | Int32 | +| CELLA_X | h.cella_x | Float32 | +| CELLA_Y | h.cella_y | Float32 | +| CELLA_Z | h.cella_z | Float32 | +| CELLB_alpha | h.cellb_alpha | Float32 | +| CELLB_beta | h.cellb_beta | Float32 | +| CELLB_gamma | h.cellb_gamma | Float32 | +| MAPC | h.mapc | Int32 | +| MAPR | h.mapc | Int32 | +| MAPS | h.mapc | Int32 | +| DMIN | h.dmin | Float32 | +| DMAX | h.dmax | Float32 | +| DMEAN | h.dmean | Float32 | +| ISPG | h.ispg | Int32 | +| NSYMBP | h.nsymbt | Int32 | +| EXTRA1 | h.extra1 | NTuple{8,UInt8} | +| EXYTYP | h.exttyp | String | +| NVERSION | h.nversion | Int32 | +| EXTRA2 | h.extra1 | NTuple{8,UInt8} | +| ORIGIN_X | h.origin_x | Float32 | +| ORIGIN_Y | h.origin_y | Float32 | +| ORIGIN_Z | h.origin_z | Float32 | +| MAP | h.map | String | +| MACHST | h.machst | NTuple{4,UInt8} | +| RMS | h.rms | Float32 | +| NLABL | h.nlabl | Int32 | +| LABEL | h.label | NTuple{10,String} | + +## Keeping the header and data in sync + +Update the header stored in `MRCData` from the data and its extended header. + +```julia +updateheader!(dmap) +``` diff --git a/docs/src/examples.md b/docs/src/examples.md new file mode 100644 index 0000000..0bac2e7 --- /dev/null +++ b/docs/src/examples.md @@ -0,0 +1,32 @@ +# MRCFile Example + +## Example 1 + +This example downloads a map of [TRPV1](https://www.emdataresource.org/EMD-5778) and animates slices taken through the map. + +To set-up this example, install FTPClient and Plots with + +```julia +using Pkg +Pkg.add("FTPClient") +Pkg.add("Plots") +``` + +```julia +using MRCFile, FTPClient, Plots + +emdid = 5778 # TRPV1 +ftp = FTP(hostname = "ftp.rcsb.org/pub/emdb/structures/EMD-$(emdid)/map") +dmap = read(download(ftp, "emd_$(emdid).map.gz"), MRCData) +close(ftp) +dmin, dmax = extrema(header(dmap)) +drange = dmax - dmin + +anim = @animate for xsection in eachmapsection(dmap) + plot(RGB.((xsection .- dmin) ./ drange)) +end + +gif(anim, "emd-$(emdid)_slices.gif", fps = 30) +``` + +![EMD-5778 slices](https://github.com/sethaxen/MRCFile.jl/blob/main/docs/src/assets/emd-5778_slices.gif) \ No newline at end of file diff --git a/docs/src/index.md b/docs/src/index.md deleted file mode 100644 index 9778506..0000000 --- a/docs/src/index.md +++ /dev/null @@ -1,7 +0,0 @@ -# MRCFile - -MRCFile.jl implements the [MRC2014 format](https://www.ccpem.ac.uk/mrc_format/mrc2014.php) for storing image and volume data such as those produced by electron microscopy. -It offers the ability to read, edit, and write MRC files, as well as utility functions for extracting useful information from the headers. - -The key type is `MRCData`, which contains the contents of the MRC file, accessible with `header` and `extendedheader`. -`MRCData` is an `AbstractArray` whose elements are those of the data portion of the file and can be accessed or modified accordingly. diff --git a/docs/src/index.md b/docs/src/index.md new file mode 120000 index 0000000..fe84005 --- /dev/null +++ b/docs/src/index.md @@ -0,0 +1 @@ +../../README.md \ No newline at end of file diff --git a/src/data.jl b/src/data.jl index c9f72cf..202546d 100644 --- a/src/data.jl +++ b/src/data.jl @@ -25,24 +25,26 @@ Create an array of the specified size. struct MRCData{T<:Number,N,EH,D} <: AbstractArray{T,N} header::MRCHeader extendedheader::EH - data::D + original_data::D +end +function ReadMRCData(header, extendedheader, data::AbstractArray{T,N}) where {T,N} + return MRCData{T,N,typeof(extendedheader),typeof(data)}(header, extendedheader, data) end function MRCData(header, extendedheader, data::AbstractArray{T,N}) where {T,N} + data = permutedims(data, reverse(1:N)) return MRCData{T,N,typeof(extendedheader),typeof(data)}(header, extendedheader, data) end function MRCData(header=MRCHeader(), extendedheader=MRCExtendedHeader()) data_size = size(header) - data_length = prod(data_size) dtype = datatype(header) dims = ndims(header) - s = ntuple(i -> data_size[i], dims) - data = Array{dtype,dims}(undef, s) - return MRCData(header, extendedheader, data) + original_data = Array{dtype,dims}(undef, reverse(data_size)) + return ReadMRCData(header, extendedheader, original_data) end function MRCData(size::NTuple{3,<:Integer}) header = MRCHeader() for i in 1:3 - setproperty!(header, (:nx, :ny, :nz)[i], size[i]) + setproperty!(header, (:nz, :ny, :nx)[i], size[i]) end return MRCData(header) end @@ -62,6 +64,13 @@ Get extended header. """ extendedheader(d::MRCData) = d.extendedheader +""" + data(d::MRCData) -> AbstractArray + +Get data. +""" +data(d::MRCData) = parent(d) + """ updateheader!(data::MRCData; statistics = true) @@ -74,7 +83,7 @@ function updateheader!(d::MRCData; statistics=true) # update size s = size(d) for i in eachindex(s) - setproperty!(h, (:nx, :ny, :nz)[i], s[i]) + setproperty!(h, (:nz, :ny, :nx)[i], s[i]) end # update space group @@ -104,7 +113,9 @@ function updateheader!(d::MRCData; statistics=true) end # Array overloads -@inline Base.parent(d::MRCData) = d.data +@inline function Base.parent(d::MRCData) + return PermutedDimsArray(d.original_data, reverse(1:ndims(d.original_data))) +end @inline Base.getindex(d::MRCData, idx::Int...) = getindex(parent(d), idx...) @@ -134,21 +145,21 @@ end Return an iterator over map rows. """ -eachmaprow(d::MRCData) = eachslice(d; dims=header(d).mapr) +eachmaprow(d::MRCData) = eachslice(d; dims=ndims(d) - header(d).mapr + 1) """ eachmapcol(d::MRCData) Return an iterator over columns. """ -eachmapcol(d::MRCData) = eachslice(d; dims=header(d).mapc) +eachmapcol(d::MRCData) = eachslice(d; dims=ndims(d) - header(d).mapc + 1) """ eachmapsection(d::MRCData) Return an iterator over sections. """ -eachmapsection(d::MRCData) = eachslice(d; dims=header(d).maps) +eachmapsection(d::MRCData) = eachslice(d; dims=ndims(d) - header(d).maps + 1) """ eachstackunit(d::MRCData) diff --git a/src/header.jl b/src/header.jl index 270c437..bd89365 100644 --- a/src/header.jl +++ b/src/header.jl @@ -136,7 +136,7 @@ end return offsets end -Base.size(h::MRCHeader) = (h.nx, h.ny, h.nz) +Base.size(h::MRCHeader) = (h.nz, h.ny, h.nx) Base.length(h::MRCHeader) = prod(size(h)) diff --git a/src/io.jl b/src/io.jl index b84d765..bbb4b54 100644 --- a/src/io.jl +++ b/src/io.jl @@ -19,9 +19,9 @@ function Base.read(io::IO, ::Type{T}; compress=:auto) where {T<:MRCData} header = read(newio, MRCHeader) extendedheader = read(newio, MRCExtendedHeader; header=header) d = MRCData(header, extendedheader) - read!(newio, d.data) + read!(newio, d.original_data) close(newio) - map!(bswaptoh(header.machst), d.data, d.data) + map!(bswaptoh(header.machst), d.original_data, d.original_data) return d end function Base.read(fn::AbstractString, ::Type{T}; compress=:auto) where {T<:MRCData} @@ -64,8 +64,8 @@ function read_mmap(io::IO, ::Type{MRCData}) head = read(io, MRCHeader) exthead = read(io, MRCExtendedHeader; header=head) arraytype = Array{datatype(head),ndims(head)} - data = Mmap.mmap(io, arraytype, size(head)[1:ndims(head)]) - return MRCData(head, exthead, data) + data = Mmap.mmap(io, arraytype, reverse(size(head))) + return ReadMRCData(head, exthead, data) end function read_mmap(path::AbstractString, T::Type{MRCData}) return open(path, "r") do io @@ -105,7 +105,7 @@ function Base.write( sz = write(newio, h) sz += write(newio, extendedheader(d)) T = datatype(h) - data = parent(d) + original_data = d.original_data fswap = bswapfromh(h.machst) buffer_size = div(buffer_size, sizeof(T)) if buffer === nothing @@ -120,15 +120,15 @@ function Base.write( # If `buffer` was provided as a parameter then `buffer_size` is redundant and # we must make sure that it matches `buffer`. buffer_size = length(buffer) - vlen = length(data) + vlen = length(original_data) vrem = vlen % buffer_size @inbounds @views begin if vrem != 0 - buffer[1:vrem] .= fswap.(T.(data[1:vrem])) + buffer[1:vrem] .= fswap.(T.(original_data[1:vrem])) sz += write(newio, buffer[1:vrem]) end for i in (vrem + 1):buffer_size:vlen - buffer .= fswap.(T.(data[i:(i + buffer_size - 1)])) + buffer .= fswap.(T.(original_data[i:(i + buffer_size - 1)])) sz += write(newio, buffer) end end diff --git a/test/consistency.jl b/test/consistency.jl index 970b455..efd24b7 100644 --- a/test/consistency.jl +++ b/test/consistency.jl @@ -6,12 +6,8 @@ using PythonCall numpy = pyimport("numpy") mrcfile = pyimport("mrcfile") -function compare_map(map_jl::Array{Float32,3}, map_py::Array{Float32,3}) - # dimensions permuted due to row-major storage in Python and col-major in Julia - # see https://github.com/sethaxen/MRCFile.jl/issues/10 - @test_broken map_jl == map_py - permuted_py = permutedims(map_py, (3, 2, 1)) - @test map_jl == permuted_py +function compare_map(map_jl::AbstractArray{Float32,3}, map_py::AbstractArray{Float32,3}) + @test map_jl == map_py end function compare_header(header_jl::MRCHeader, header_py::Py) @@ -60,7 +56,7 @@ end function compare_mrcfile(map_path::String) emd_jl = read(map_path, MRCData) - map_jl = emd_jl.data + map_jl = data(emd_jl) header_jl = header(emd_jl) exh_jl = extendedheader(emd_jl) diff --git a/test/data.jl b/test/data.jl index 2125b06..f44f4a8 100644 --- a/test/data.jl +++ b/test/data.jl @@ -2,11 +2,11 @@ @testset "MRCData(header, extendedheader, data)" begin h = MRCHeader() exh = MRCExtendedHeader() - data = Array{Float32,3}(undef, (2, 2, 2)) - d = MRCData(h, exh, data) + map_data = Array{Float32,3}(undef, (2, 2, 2)) + d = MRCData(h, exh, map_data) @test d.header === h @test d.extendedheader === exh - @test d.data === data + @test data(d) == map_data @test eltype(d) === Float32 @test ndims(d) === 3 end @@ -24,15 +24,16 @@ @testset "MRCData(size::NTuple{3,<:Integer})" begin d = MRCData((10, 20, 30)) @test size(d) == (10, 20, 30) - @test size(d.data) == (10, 20, 30) + @test size(data(d)) == (10, 20, 30) + @test size(d.original_data) == (30, 20, 10) @test ndims(d) == 3 - @test ndims(d.data) == 3 + @test ndims(d.original_data) == 3 @testset "header(d::MRCData)" begin h = header(d) - @test h.nx == 10 + @test h.nx == 30 @test h.ny == 20 - @test h.nz == 30 + @test h.nz == 10 end end @@ -42,10 +43,10 @@ @test size(d1) == size(d2) end - data = fill(1.0f0, (2, 2, 2)) + map_data = fill(1.0f0, (2, 2, 2)) h = MRCHeader() exh = MRCExtendedHeader() - d = MRCData(h, exh, data) + d = MRCData(h, exh, map_data) @testset "getindex" begin @test d[1, 1, 1] == 1.0f0 @@ -109,18 +110,18 @@ end @testset "IndexStyle" begin - @test Base.IndexStyle(typeof(d)) == Base.IndexStyle(typeof(d.data)) + @test Base.IndexStyle(typeof(d)) == Base.IndexStyle(typeof(d.original_data)) end end @testset "updateheader!(d::MRCData; statistics=true)" begin - data = fill(1.0f0, (2, 2, 2)) + map_data = fill(1.0f0, (2, 2, 2)) h = MRCHeader() @test h.nx == 0 @test h.ny == 0 @test h.nz == 0 exh = MRCExtendedHeader() - d = MRCData(h, exh, data) + d = MRCData(h, exh, map_data) updateheader!(d) @test d.header.dmin == 1.0f0 @test d.header.dmax == 1.0f0 diff --git a/test/header.jl b/test/header.jl index 6c08d29..2351e60 100644 --- a/test/header.jl +++ b/test/header.jl @@ -178,7 +178,7 @@ end @testset "size(::MRCHeader)" begin h = MRCHeader(; nx=10, ny=20, nz=30) - @test size(h) == (10, 20, 30) + @test size(h) == (30, 20, 10) end @testset "length(::MRCHeader)" begin diff --git a/test/io.jl b/test/io.jl index f0fe4bf..ff80059 100644 --- a/test/io.jl +++ b/test/io.jl @@ -8,9 +8,9 @@ h = header(emd3001) eh = extendedheader(emd3001) p = parent(emd3001) - @test h.nx === Int32(size(p)[1]) === Int32(73) + @test h.nx === Int32(size(p)[3]) === Int32(73) @test h.ny === Int32(size(p)[2]) === Int32(43) - @test h.nz === Int32(size(p)[3]) === Int32(25) + @test h.nz === Int32(size(p)[1]) === Int32(25) @test MRCFile.datatype(h.mode) === eltype(p) === Float32 @test h.nxstart === Int32(0) @test h.nystart === Int32(-21)