JuliaAstro / FITSIO.jl

Flexible Image Transport System (FITS) file support for Julia
http://juliaastro.org/FITSIO.jl/
MIT License
55 stars 29 forks source link

Added a method to read an image into a pre-allocated array #122

Closed jishnub closed 4 years ago

jishnub commented 4 years ago

The read method currently used to allocate an array of the correct type and size and read from disk into it. This leads to temporary allocation that might hamper performance if data is read in a loop. In this PR a new method has been added to the non-allocating read! that lets us read from disk into a pre-allocated array.

jishnub commented 4 years ago

Additionally, https://github.com/JuliaAstro/FITSIO.jl/blob/450ddb6d839b11f5a6b0f5f2a612463fcf4a8e71/src/libcfitsio.jl#L799-L801

Does data need to be an array? I've tried redefining the method to match an AbstractArray, which lets me read data into views of an array as

julia> a=FITS("abc.fits","w"); write(a,ones(3,3)); close(a); a=FITS("abc.fits","r");

julia> b=zeros(3,3);

julia> read!(a[1],view(b,:,1),:,1)
3-element view(::Array{Float64,2}, :, 1) with eltype Float64:
 1.0
 1.0
 1.0

julia> b
3×3 Array{Float64,2}:
 1.0  0.0  0.0
 1.0  0.0  0.0
 1.0  0.0  0.0

Similarly in fits_read_pix. The tests seem to pass with this change.

giordano commented 4 years ago

Reading to a view is interesting, I happen to have this pending pull request for Julia to make Base.read! accept views: https://github.com/JuliaLang/julia/pull/33046 (I should push it forward, by the way). AbstractArray is definitely too loose, even though we could just rely on missing methods to make it work, or you can look at the type I've used in my PR (Union{Array{T}, FastContiguousSubArray{T,N,<:Array{T}} where N}).

giordano commented 4 years ago

Yesterday I tested this PR locally and observed that read!(f1[1],b) does allocate something, even though the memory used is considerably less than read(f1[1]) and read! in the end is indeed slightly faster than plain read. I couldn't profile memory allocations, but they may well come from CFITSIO library.

jishnub commented 4 years ago

Reading to a view is interesting, I happen to have this pending pull request for Julia to make Base.read! accept views: JuliaLang/julia#33046 (I should push it forward, by the way). AbstractArray is definitely too loose, even though we could just rely on missing methods to make it work, or you can look at the type I've used in my PR (Union{Array{T}, FastContiguousSubArray{T,N,<:Array{T}} where N}).

Union{Array{T}, FastContiguousSubArray{T,N,<:Array{T}} where N} does seem to be the right way to implement this, I didn't know that FastContiguousSubArray exists.

kbarbary commented 4 years ago

Looks good to me!

giordano commented 4 years ago

Thank you so much!