Change to row-major order and modify docs

Project.toml

@@ -1,7 +1,8 @@
 name = "MRCFile"
 uuid = "c18c01d3-0ab9-49c3-bd2d-8f2e64a2b7a5"
 authors = ["Seth Axen <[email protected]>"]
-version = "0.1.2"
+version = "0.1.3"
+
 
 [deps]
 CodecBzip2 = "523fee87-0ab8-5b00-afb7-3ecf72e48cfd"

README.md

@@ -1,24 +1,37 @@
 # 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)
 [![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://sethaxen.github.io/MRCFile.jl/stable)
 [![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
 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/).

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")

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)
+```

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)

docs/src/index.md

@@ -0,0 +1 @@
+../../README.md

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)

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))
 

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

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)