Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Get initial traits in place for NeuroCore #1

Merged
merged 15 commits into from
Dec 27, 2019
Merged
Show file tree
Hide file tree
Changes from 1 commit
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Prev Previous commit
Next Next commit
Get initial documentation in place before exhaustive testing
  • Loading branch information
Tokazama committed Dec 12, 2019
commit adc0cd330f8ebde2efd69898ba6f9a72d1b9235f
10 changes: 8 additions & 2 deletions src/NeuroCore.jl
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,8 @@ export freqdim,
slice_end!,
slice_duration,
slice_duration!,
slice_timing,
slice_timing!,
echo_time,
echo_time!,
inversion_time,
Expand Down Expand Up @@ -104,8 +106,12 @@ export freqdim,
total_readout_time,
total_readout_time!
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I feel that this is a large number of names to export, especially given

  • A lot of these seem to be trivial wrappers around an internal dictionary-like interface.
  • There's a large number of other attributes which aren't covered here, and achieving full coverage of the very large number of DICOM attributes seems impractical.

Did you consider a more dictionary-like interface for some of these attributes? For example, perhaps allowing syntax such as

metadata(x).manufacturer

which is possible to make type stable if done carefully

(Achieving x.manufacturer is possible and tantalizing but likely undesirable due to the need to "steal" getproperty on some non-concrete supertype of x.)

I feel there may be two rough categories of metadata:

  1. Metadata which is directly required for working with the pixel data. For example, pixel spacings and coordinate systems.
  2. Metadata which is somewhat irrelevant or only indirectly relevant, such as the institution name.

The first type is a more limited set of metadata and we may plausibly get a minimal required set of functions for accessing this like ImageCore. The latter type of metadata is fairly open ended and I'd suggest it's better handled with some kind of dictionary like lookup. Does the images ecosystem also have a similar distinction?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There's a large number of other attributes which aren't covered here, and achieving full coverage of the very large number of DICOM attributes seems impractical.

I agree. That's why I've focused on what's in BIDS. It's still a lot but theoretically this can be broken up into smaller packages if it's found to be useful (eg., BIDSImage, BIDSEeg, BIDSData, etc).

The latter type of metadata is fairly open ended and I'd suggest it's better handled with some kind of dictionary like lookup.

A large portion of this is in fact a wrapper around a dictionary lookup but makes it type stable. JuliaImages uses properties to retrieve a dictionary with keytype == String and valtype == Any. ImageMeta has a nice some nice syntax that allows img[::String] to search the underlying property structure. However, that would restrict everything to have ImageMeta as the top level wrapper where it might not be appropriate (e.g., non arrays such as graphs, fiber tracks, etc.).

One example where this could be immediately useful is if you look at MRIReco.jl where there are distinct structures for raw acquisitional data, sequence parameters, etc. But very little of that is needed in NIfTI files. There's some overlap but it's not enough to agree upon a single structure for say acquisition parameters. This example doesn't even account for GIFTI, CIFTI, BED , etc. where there the data isn't an image or is a mesh instead.

All that being said, I'm certainly open to a better solution. I just don't think there will ever be a single structure everyone can agree on which as you pointed out means using a dictionary. If there's another way of guaranteeing that I get Float64 everytime I call img["EchoTime"] then we could just require that syntax.

Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If there's another way of guaranteeing that I get Float64 everytime I call img["EchoTime"] then we could just require that syntax.

Sure. Very basic sketch:

struct Metadata
     d::Dict{Symbol,Any}
end

function Base.getproperty(m::Metadata, s::Symbol)
    if s === :EchoTime
        getfield(m, :d)[:EchoTime]::Float64 # Inferrable
    else
        getfield(m, :d)[s]  # Fallback (uninferrable)
    end
end
julia> m = Metadata(Dict(:EchoTime=>1.0, :Foo=>"hi"))
Metadata(Dict{Symbol,Any}(:EchoTime => 1.0,:Foo => "hi"))

julia> m.EchoTime
1.0

julia> m.Foo
"hi"

julia> @code_warntype (m->m.EchoTime)(m)
Variables
  #self#::Core.Compiler.Const(var"#15#16"(), false)
  m::Metadata

Body::Float64
1 ─ %1 = Base.getproperty(m, :EchoTime)::Float64
└──      return %1

julia> @code_warntype (m->m.Foo)(m)
Variables
  #self#::Core.Compiler.Const(var"#17#18"(), false)
  m::Metadata

Body::Any
1 ─ %1 = Base.getproperty(m, :Foo)::Any
└──      return %1

Then provide a generic function metadata(x) which for any x should return something with the same property-like interface as Metadata (it could literally be a concrete Metadata, or some other package-specific type).

Copy link

@c42f c42f Dec 17, 2019

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Clearly, that type assertion for ::Float64 is a simplification and you may want to dispatch to extra logic there.

The main point is that getproperty can be used to look up properties in a flexible backing dictionary and certain blessed property names can be treated in a special way which makes them inferrable.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I like that approach a lot. I would need to create a unique array type that reimplements all of the dictionary traits in JuliaImages that using String as the keytype though. That or convince JuliaImages to change to a Symbol keytype.

Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Cool, it will be interesting to see how it goes.

Note also that Metadata in my sketch could be implemented with actual fields for a subset of "expected" metadata, and a dict as a fallback for the rest. This would avoid the dict lookup for the set of blessed attributes, and would therefore be more efficient provided x keeps the Metadata in this format when it's constructed so that no work is required to compute metadata(x). Of course, you may want to call it MRIMetadata or some such, given most of these fields are very MRI-specific. Other types of x would have their own type with the same API but different blessed attributes.

struct MRIMetadata
     echo_time::Union{Float64,Nothing}
     # ...
     d::Dict{Symbol,Any}
end

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This leads me to a somewhat ignorant question. When I call sizeof on your structure it doesn't seem to matter if I use Nothing or Float64. Would it make sense to do something like this instead:

struct MRIMetadata{ET<:Union{Float64,Nothing}}
    echo_time::ET
   d::Dict{Symbol,Any}
end

This would prevent allocating memory for things that only have a couple fields that aren't nothing. I know this may be trivial right now, but if we want to have this metadata be flexible enough to work with something like individual fibers tracking algorithms, then it could get pretty expensive.

Copy link

@c42f c42f Dec 19, 2019

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think that parameterizing MRIMetadata on the presence or absence of metadata fields would be unnecessarily hard on the compiler and confer no practical advantage (it will need to create a new type each time the metadata changes). If you're worried about storage, you could go with a pure-Dict based solution instead. That would be simpler anyway so it might be for the best. Dict lookup is quite fast and is unlikely to be done in any inner analysis loop.

(Besides this, I would have imagined the relative sizeof(MRIMetdata) to sizeof(image_data) could easily be 1:1000 in which case it's not worth worrying about.)

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think that parameterizing MRIMetadata on the presence or absence of metadata fields would be unnecessarily hard on the compiler and confer no practical advantage

I should probably clarify. I imagined that most users wouldn't ever construct these and that it would serve as a common way of organizing and accessing metadata for those actually creating packages.

struct Metadata{ET<:Union{Float64,Nothing}}
    echo_time::ET
   d::Dict{Symbol,Any}
end

const EmptMetadata =Metadata{Nothing}

This would allow the person implementing an EEG file reader to not worry about echo time but still use the same interface. As you said, it may be simpler to just go with a dictionary.

I would have imagined the relative sizeof(MRIMetdata) to sizeof(image_data) could easily be 1:1000 in which case it's not worth worrying about.

I agree with this in so far as it pertains to images, but this should be applicable outside of MRI as well. For example, we could have a small connectome composed of a 13x13 array across 10,000 subjects. This is definitely the extreme but I'd hate to paint myself into a corner early on.

Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If it's worrying you, I'd suggest go for the Dict. It will be less code anyway :-)

But I also don't think you're painting yourself into a corner too much by adopting one implementation or another. The public interface to this stuff will be getproperty, which is allowed to present the metadata as if it came from an actual field. But whether there's an actual field present is an implementation detail.

I suspect what's most important / potentially annoying here is to settle on appropriate names for the metadata fields because these are very much the public interface. For example echo_time vs EchoTime. Probably the only viable thing is to stick with the original field names as defined in the BIDS standard, even though they conflict with the Julia naming conventions.



include("neuroproperty.jl")
include("properties.jl")
include("task_events.jl")
include("spatial_properties.jl")
include("modality_agnostic.jl")
include("hardware.jl")
include("sequence.jl")
include("time.jl")

end
46 changes: 46 additions & 0 deletions src/abstractarrays.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,46 @@

using AbstractIndices, ImageMetadata, ImageCore

const MetaIndicesArray{T,N,A,I} = IndicesArray{T,N,ImageMeta{T,N,A},I}

const NamedIndicesArray{L,T,N,A,I} = NamedDimsArray{L,T,N,IndicesArray{T,N,A,I}}

const ImageArray{L,T,N,A,I} = NamedIndicesArray{L,T,N,ImageMeta{T,N,A},I}

ImageCore.HasProperties(::Type{T}) where {T<:NamedDimsArray} = HasProperties(parent_type(T))

ImageCore.HasProperties(::Type{T}) where {T<:IndicesArray} = HasProperties(parent_type(T))


# ImageMetadata
ImageMetadata.properties(img::ImageArray) = properties(parent(img))
ImageMetadata.properties(img::MetaIndicesArray) = properties(parent(img))

function ImageMetadata.spatialproperties(img::ImageArray)
return ImageMetadata.@get img "spatialproperties" ["spacedirections"]
end

Base.delete!(img::ImageArray, propname::AbstractString) = delete!(properties(img), propname)

Base.get(img::ImageArray, k::AbstractString, default) = get(properties(img), k, default)

Base.haskey(img::ImageArray, k::AbstractString) = haskey(properties(img), k)

is_color_axis(::NamedTuple{(:color,)}) = true

first_pair(x::NamedTuple{names}) where {names} = NamedTuple{(first(names),)}((first(x),))

"""
is_color_axis(ni::Index) -> Bool

Determines whether a given axis refers to a color dimension.
"""
is_color_axis(::T) where {T} = is_color_axis(T)
is_color_axis(::Type{<:Index{:color}}) = true
is_color_axis(::Type{<:Index{name}}) where {name} = false

colordim(x) = find_axis(is_color_axis, x)

function ImageCore.spacedirections(img::ImageArray)
return getter(img, "spacedirections", axes_type(img), spacedirections(namedaxes(img)))
end
174 changes: 174 additions & 0 deletions src/coordinates.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,174 @@
struct CoordinateSystem{S} end

"""
UnkownSpace
"""
const UnkownSpace = CoordinateSystem{:unkown}()
Tokazama marked this conversation as resolved.
Show resolved Hide resolved


"""
ScannerSpace

In scanner space.
"""
const ScannerSpace = CoordinateSystem{:scanner}()


"""
AnatomicalSpace

equivalent of 'aligned' space in NIfTI standard.
"""
const AnatomicalSpace = CoordinateSystem{:anatomical}()

"""
TailarachSpace

Tailarach space
"""
const TailarachSpace = CoordinateSystem{:tailarach}()

"""
MNI152Space

MNI152 space
"""
const MNI152Space = CoordinateSystem{:MNI152}()

"""
coordinate_system(x)

Return the coordinate space that `x` is in.
"""
coordinate_system(x::Any) = getter(x, "coordinatespace", CoordinateSystem, UnkownSpace)

"""
coordinate_system!(x, val)

Set the coordinate space for `x` to `val`.
"""
function coordinate_system!(x::Any, val::CoordinateSystem)
setter!(x, "coordinate_system", val, CoordinateSystem)
end


# TODO
"""
slice_direction

Possible values: `i`, `j`, `k`, `i-`, `j-`, `k-` (the axis of the NIfTI data
along which slices were acquired, and the direction in which `SliceTiming` is
defined with respect to). `i`, `j`, `k` identifiers correspond to the first,
second and third axis of the data in the NIfTI file. A `-` sign indicates that
the contents of `SliceTiming` are defined in reverse order - that is, the first
entry corresponds to the slice with the largest index, and the final entry
corresponds to slice index zero. When present, the axis defined by
`SliceEncodingDirection` needs to be consistent with the ‘slice_dim’ field in
the NIfTI header. When absent, the entries in `SliceTiming` must be in the
order of increasing slice index as defined by the NIfTI header. |
"""
function slice_direction end


"""
sform_code(x) -> CoordinateSystem

Code describing the orientation of the image.

Should only be one of following (although others are allowed):

* UnkownSpace
* AnatomicalSpace
* TalairachSpace
* MNI152Space
"""
sform_code(x::Any) = getter(x, "sform_code", CoordinateSystem, UnkownSpace)

"""
sform_code!(x, val)

Set the `sform` coordinate space of `x` to `val`.
"""
sform_code!(x::Any, val::CoordinateSystem) = setter!(x, "sform_code", val, CoordinateSystem)

"""
qform_code(x) -> CoordinateSystem

Code describing the orientation of the image in the scanner.

Should only be one of the following (although others are allowed):

* UnkownSpace
* ScannerSpace
"""
qform_code(x::Any) = getter(x, "sform_code", CoordinateSystem, UnkownSpace)

"""
qform_code!(x, val)

Set the `qfrom` coordinate space of `x` to `val`.
"""
qform_code!(x::Any, val::CoordinateSystem) = setter!(x, "qform_code", val, CoordinateSystem)

"""
qform(x) -> MMatrix{4,4,Float64,16}
"""
qform(x::Any) = getter(x, "qform", MMatrix{4,4,Float64,16}, i->default_affinemat(i))

qform!(x::Any, val::AbstractMatrix) = setter!(x, "qform", val, MMatrix{4,4,Float64,16})

"""
sform(x) -> MMatrix{4,4,Float64,16}

The 4th column of the matrix is the offset of the affine matrix.
This is primarily included for the purpose of compatibility with DICOM formats, where the
"Image Position" stores the coordinates of the center of the first voxel
(see the [DICOM standard](http://dicom.nema.org/medical/dicom/current/output/chtml/part03/sect_C.7.6.2.html#sect_C.7.6.2.1.1) for more details;
Note, these values should be in interpreted as 'mm').
"""
sform(x::Any) = getter(x, "sform", MMatrix{4,4,Float64,16}, i->default_affinemat(i))

sform!(x::Any, val::AbstractMatrix) = setter!(x, "sform", val, MMatrix{4,4,Float64,16})


function default_affinemat(x::Any)
if sform_code(x) === UnkownSpace
if qform_code(x) === UnkownSpace
return _default_affinemat(x)
else
return qform(x)
end
else
return sform(x)
end
end

function _default_affinemat(x::Any)
_default_affinemat(spacedirections(x), pixelspacing(x))
end


function _default_affinemat(sd::NTuple{2,NTuple{2,T}}, ps::NTuple{2,T}) where T
MMatrix{4,4,Float64,16}(sd[1][1], sd[2][1], 0, 0,
sd[1][2], sd[2][2], 0, 0,
0, 0, 0, 0,
ps[1], ps[2], 0, 0)
end

function _default_affinemat(sd::NTuple{3,NTuple{3,T}}, ps::NTuple{3,T}) where T
MMatrix{4,4,Float64,16}(sd[1][1], sd[2][1], sd[3][1], 0,
sd[1][2], sd[2][2], sd[3][2], 0,
sd[1][3], sd[2][3], sd[3][3], 0,
ps[1], ps[2], ps[3], 0)
end

"""
affine_matrix(x) -> MMatrix{4,4,Float64,16}

Returns affine matrix. For an instance that returns `spacedirections` this is
the corresponding tuple converted to an array.
"""
affine_matrix(x::Any) = getter(x, "affine_matrix", MMatrix{4,4,Float64,16}, i -> default_affinemat(i))
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Did you consider using CoordinateTransformations.AffineMap for affine maps?

Personally I consider the use of 4x4 matrices for affine transformations to be an unfortunate failure of abstraction: applying them to lenght-3 vectors involves lifting those vectors into homogenous coordinates and back again, and this lifting provides extra complication without any real value.

As far as I know, the only time homogenous coordinates are really worthwhile is

  • In computer vision applications where the set of 4x4 matrices also supports projection and there are some neat algebraic identities which can help with things like camera calibration.
  • In languages which have builtin matrices but don't have free abstractions. Like matlab.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That's what I was originally using but then there were so many weird transformations that NIfTI required that I just used a 4x4 matrix. That being said, I think you're right about this and it would be better to use CoordinateTransformations.AffineMap.



affine_matrix!(x::Any, val::AbstractMatrix) = setter!(x, "affine_matrix", val, MMatrix{4,4,Float64,16})
Loading