diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index ad5dc637..88d4a4a6 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -25,12 +25,12 @@ jobs: arch: - x64 steps: - - uses: actions/checkout@v2 + - uses: actions/checkout@v4 - uses: julia-actions/setup-julia@v1 with: version: ${{ matrix.version }} arch: ${{ matrix.arch }} - - uses: actions/cache@v1 + - uses: actions/cache@v3 env: cache-name: cache-artifacts with: @@ -43,14 +43,14 @@ jobs: - uses: julia-actions/julia-buildpkg@v1 - uses: julia-actions/julia-runtest@v1 - uses: julia-actions/julia-processcoverage@v1 - - uses: codecov/codecov-action@v1 + - uses: codecov/codecov-action@v3 with: file: lcov.info docs: name: Documentation runs-on: ubuntu-latest steps: - - uses: actions/checkout@v2 + - uses: actions/checkout@v4 - uses: julia-actions/setup-julia@v1 with: version: '1' diff --git a/src/FrequencyFilter.jl b/src/FrequencyFilter.jl index 0bd1b4f8..ba6a0a05 100644 --- a/src/FrequencyFilter.jl +++ b/src/FrequencyFilter.jl @@ -173,12 +173,12 @@ function filterFrequenciesByStepsize!(indices::Vector{CartesianIndex{2}}, stepsi filter!(x -> in(x[1], stepIndices), indices) end -function sortFrequencies!(indices::Vector{CartesianIndex{2}}, f::MPIFile; numPeriodGrouping = 1, stepsize = 1, sortBySNR = false, sortByMixFactors = false) +function sortFrequencies(indices::Vector{CartesianIndex{2}}, f::MPIFile; numPeriodGrouping = 1, stepsize = 1, sortBySNR = false, sortByMixFactors = false) if sortBySNR && !sortByMixFactors indices = sortFrequenciesBySNR(indices, f, numPeriodGrouping = numPeriodGrouping, stepsize = stepsize) elseif !sortBySNR && sortByMixFactors indices = sortFrequenciesByMixFactors(indices, f, numPeriodGrouping = numPeriodGrouping) - else + elseif sortBySNR && sortByMixFactors error("Can not apply multiple sorting algorithms to frequencies") end return indices diff --git a/src/Positions/Positions.jl b/src/Positions/Positions.jl index 57b3ab6f..f8a53664 100644 --- a/src/Positions/Positions.jl +++ b/src/Positions/Positions.jl @@ -52,6 +52,8 @@ function Positions(params::Dict) positions = UniformRandomPositions(params) elseif typ == "ArbitraryPositions" positions = ArbitraryPositions(params) + elseif typ == "TubularRegularGridPositions" + positions = TubularRegularGridPositions(params) else throw(ErrorException("No grid found to load from dict $params")) end @@ -542,6 +544,162 @@ function convert(::Type{UniformRandomPositions}, N::Integer,fov::Vector,center:: return UniformRandomPositions(N,rand(UInt32),fov,center) end =# + +# Tubular cartesian grid +export TubularRegularGridPositions +mutable struct TubularRegularGridPositions{T} <: GridPositions where {T<:Unitful.Length} + shape::Vector{Int} + fov::Vector{T} + center::Vector{T} + "Main axis of the tube; only effective in 3D grids" + mainAxis::Int64 + "Radius-defining axis of the tube" + radiusAxis::Int64 +end + +function TubularRegularGridPositions(params::Dict) + shape = params["positionsShape"] + + fov = params["positionsFov"] + if !(eltype(fov) <: Quantity) + fov = fov.*Unitful.m + end + + center = params["positionsCenter"] + if !(eltype(center) <: Quantity) + center = center.*Unitful.m + end + + mainAxis = params["positionsMainAxis"] + radiusAxis = params["positionsRadiusAxis"] + + return TubularRegularGridPositions(shape, fov, center, mainAxis, radiusAxis) +end + +function TubularRegularGridPositions(file::HDF5.File) + shape = read(file, "/positionsShape") + fov = read(file, "/positionsFov")*Unitful.m + center = read(file, "/positionsCenter")*Unitful.m + mainAxis = read(file, "positionsMainAxis") + radiusAxis = read(file, "positionsRadiusAxis") + + return TubularRegularGridPositions(shape, fov, center, mainAxis, radiusAxis) +end + +length(grid::TubularRegularGridPositions) = length(filteredPositions(grid)) + +export radius +radius(grid::TubularRegularGridPositions) = grid.fov[grid.radiusAxis] / 2 + +function write(file::HDF5.File, positions::TubularRegularGridPositions) + write(file, "/positionsType", "TubularRegularGridPositions") + write(file, "/positionsShape", positions.shape) + write(file, "/positionsFov", Float64.(ustrip.(uconvert.(Unitful.m, positions.fov))) ) + write(file, "/positionsCenter", Float64.(ustrip.(uconvert.(Unitful.m, positions.center))) ) + write(file, "/positionsMainAxis", positions.mainAxis) + write(file, "/positionsRadiusAxis", positions.radiusAxis) +end + +function toDict(positions::TubularRegularGridPositions) + params = Dict{String,Any}() + params["positionsType"] = "TubularRegularGridPositions" + params["positionsShape"] = positions.shape + params["positionsFov"] = Float64.(ustrip.(uconvert.(Unitful.m, positions.fov))) + params["positionsCenter"] = Float64.(ustrip.(uconvert.(Unitful.m, positions.center))) + params["positionsMainAxis"] = positions.mainAxis + params["positionsRadiusAxis"] = positions.radiusAxis + return params +end + +function filteredPositions(grid::TubularRegularGridPositions) + if length(grid.shape) == 1 #Very ugly but improves compile time + cartIndices = CartesianIndices(tuple(grid.shape[1])) + return [Tuple(idx) for idx in cartIndices] # Applying a radius in 1D is not possible + elseif length(grid.shape) == 2 + cartIndices = CartesianIndices(tuple(grid.shape[1], grid.shape[2])) + return [Tuple(idx) for idx in cartIndices if norm(collect(Tuple(idx)) .- grid.shape./2)] + elseif length(grid.shape) == 3 + cartIndices = CartesianIndices(tuple(grid.shape[1], grid.shape[2], grid.shape[3])) + else + cartIndices = CartesianIndices(tuple(grid.shape...)) + end + + return [idx for idx in cartIndices if norm(collect(Tuple(idx)) .- grid.shape./2 .- 0.5) <= grid.shape[grid.radiusAxis]/2] +end + +function getindex(grid::TubularRegularGridPositions, i::Integer) + filteredPositions_ = filteredPositions(grid) + if i > length(filteredPositions_) || i < 1 + throw(BoundsError(grid, i)) + end + + idx = Tuple(filteredPositions_[i]) + + return ((-shape(grid).+(2 .*idx.-1))./shape(grid)).*fieldOfView(grid)./2 + fieldOfViewCenter(grid) +end + +function getindex(grid::TubularRegularGridPositions, idx::Vector{T}) where T<:Number + filteredPositions_ = filteredPositions(grid) + cartIdx = CartesianIndex(idx...) + if !(cartIdx in filteredPositions_) + throw(BoundsError(grid, idx)) + else + return ((-shape(grid).+(2 .*idx.-1))./shape(grid)).*fieldOfView(grid)./2 + fieldOfViewCenter(grid) + end +end + +function posToIdxFloat(grid::TubularRegularGridPositions, pos) + idx = 0.5 .* (shape(grid) .* ((pos .- fieldOfViewCenter(grid)) ./ + ( 0.5 .* fieldOfView(grid) ) .+ 1) .+ 1) + idx = [isnan(val) ? one(eltype(idx)) : val for val in idx] + return idx +end + +function posToIdx(grid::TubularRegularGridPositions, pos) + return round.(Int64, posToIdxFloat(grid, pos)) +end + +function posToLinIdx(grid::TubularRegularGridPositions, pos) + filteredPositions_ = filteredPositions(grid) + idx = posToIdx(grid, pos) + matchingIdx = [i for (i, idx_) in enumerate(filteredPositions_) if idx_ == CartesianIndex(idx...)] + if length(matchingIdx) < 1 + error("The position $pos is not valid.") + else + return first(matchingIdx) + end +end + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + # General functions for handling grids @@ -622,7 +780,7 @@ function loadTDesign(t::Int64, N::Int64, radius::S=10Unitful.mm, center::Vector{ push!(Ns,parse(Int,N)) end sort!(Ns) - @info "No spherical $t-Design with $N points availible!\nThere are spherical $t-Designs with following N:" Ns + @info "No spherical $t-Design with $N points available!\nThere are spherical $t-Designs with following N:" Ns throw(DomainError(1)) else ts = Int[] @@ -633,7 +791,7 @@ function loadTDesign(t::Int64, N::Int64, radius::S=10Unitful.mm, center::Vector{ end end sort!(ts) - @info "No spherical $t-Design availible!\n Choose another t." + @info "No spherical $t-Design available!\n Choose another t." throw(DomainError(1)) end end diff --git a/test/MDFInMemory.jl b/test/MDFInMemory.jl index e94f8d80..23dac34f 100644 --- a/test/MDFInMemory.jl +++ b/test/MDFInMemory.jl @@ -29,19 +29,19 @@ end fnSM1DBruker = joinpath(datadir,"BrukerStore","20170807_142514_Service_1_1","89") - fnMeasV1 = joinpath(tmpdir,"mdf","measurement_V1.mdf") - fnMeasV2 = joinpath(tmpdir,"mdf","measurement_V2.mdf") - fnSMV1 = joinpath(tmpdir,"mdf","systemMatrix_V1.mdf") - fnSMV2 = joinpath(tmpdir,"mdf","systemMatrix_V2.mdf") - fnSMV3 = joinpath(tmpdir,"mdf","systemMatrix_V3.mdf") - fnSMV4 = joinpath(tmpdir,"mdf","systemMatrix_V4.mdf") - fnSMV5 = joinpath(tmpdir,"mdf","systemMatrix_V5.mdf") - fnSMV6 = joinpath(tmpdir,"mdf","systemMatrix_V6.mdf") - fnSM1DV1 = joinpath(tmpdir,"mdf","systemMatrix1D_V1.mdf") - fnSM1DV2 = joinpath(tmpdir,"mdf","systemMatrix1D_V2.mdf") - - fnMeasV2_converted = joinpath(tmpdir,"mdf","measurement_V2_converted.mdf") - fnSMV2_converted = joinpath(tmpdir,"mdf","systemMatrix_V2_converted.mdf") + fnMeasV1 = joinpath(tmpdir,"mdfim","measurement_V1.mdf") + fnMeasV2 = joinpath(tmpdir,"mdfim","measurement_V2.mdf") + fnSMV1 = joinpath(tmpdir,"mdfim","systemMatrix_V1.mdf") + fnSMV2 = joinpath(tmpdir,"mdfim","systemMatrix_V2.mdf") + fnSMV3 = joinpath(tmpdir,"mdfim","systemMatrix_V3.mdf") + fnSMV4 = joinpath(tmpdir,"mdfim","systemMatrix_V4.mdf") + fnSMV5 = joinpath(tmpdir,"mdfim","systemMatrix_V5.mdf") + fnSMV6 = joinpath(tmpdir,"mdfim","systemMatrix_V6.mdf") + fnSM1DV1 = joinpath(tmpdir,"mdfim","systemMatrix1D_V1.mdf") + fnSM1DV2 = joinpath(tmpdir,"mdfim","systemMatrix1D_V2.mdf") + + fnMeasV2_converted = joinpath(tmpdir,"mdfim","measurement_V2_converted.mdf") + fnSMV2_converted = joinpath(tmpdir,"mdfim","systemMatrix_V2_converted.mdf") @testset "Conversion" begin @testset "Fields" begin diff --git a/test/Positions.jl b/test/Positions.jl index 638a60cb..275a2bd4 100644 --- a/test/Positions.jl +++ b/test/Positions.jl @@ -299,6 +299,65 @@ pospath = joinpath(tmpdir,"positions","Positions.h5") @test p == caG[i] end + @testset "Tubular regular grid positions" begin + grid = TubularRegularGridPositions([81, 81, 1], [40, 40 ,0]u"mm", [0, 0, 0]u"mm", 3, 1) + + params = Dict{String, Any}() + params["positionsType"] = "TubularRegularGridPositions" + params["positionsShape"] = [81, 81, 1] + params["positionsFov"] = [40, 40 ,0]u"mm" + params["positionsCenter"] = [0, 0, 0]u"mm" + params["positionsMainAxis"] = 3 + params["positionsRadiusAxis"] = 1 + gridByParams = TubularRegularGridPositions(params) + + @test grid == gridByParams + + gridByParamsGeneral = Positions(params) + @test grid == gridByParamsGeneral + + paramsFromGrid = toDict(grid) + @test params["positionsType"] == paramsFromGrid["positionsType"] + @test all(params["positionsShape"] .== paramsFromGrid["positionsShape"]) + @test all(ustrip.(u"m", params["positionsFov"]) .≈ paramsFromGrid["positionsFov"]) + @test all(ustrip.(u"m", params["positionsCenter"]) .≈ paramsFromGrid["positionsCenter"]) + @test params["positionsMainAxis"] == paramsFromGrid["positionsMainAxis"] + @test params["positionsRadiusAxis"] == paramsFromGrid["positionsRadiusAxis"] + + @test length(grid) == 5169 + @test MPIFiles.radius(grid) == 20u"mm" + + filename = joinpath(tmpdir, "TubularRegularGridPositionsTest.h5") + + if isfile(filename) + rm(filename) + end + + h5open(filename, "w") do f + write(f, grid) + end + + gridByFile = nothing + h5open(filename, "r") do f + gridByFile = TubularRegularGridPositions(f) + end + + rm(filename) + + @test gridByFile isa TubularRegularGridPositions + @test all(grid.shape .≈ gridByFile.shape) + @test all(grid.fov .≈ gridByFile.fov) + @test all(grid.center .≈ gridByFile.center) + @test grid.mainAxis == gridByFile.mainAxis + @test grid.radiusAxis == gridByFile.radiusAxis + + @test all(grid[1] .≈ [-2.962962962962963, -19.753086419753085, 0.0]u"mm") + @test all(grid[2000] .≈ [-9.382716049382715, -3.45679012345679, 0.0]u"mm") + + @test posToLinIdx(grid, [-2.962962962962963, -19.753086419753085, 0.0]u"mm") == 1 + @test posToLinIdx(grid, [-9.382716049382715, -3.45679012345679, 0.0]u"mm") == 2000 + end + @testset "Squared positions regression test" begin # When supplying a dict already equipped with Unitful units, the positions were squared params = Dict{String, Any}() diff --git a/test/runtests.jl b/test/runtests.jl index d21e0e9c..994af6cb 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -17,6 +17,7 @@ const tmpdir = @get_scratch!("tmp") @info "If you want to check the output of the tests, please head to $tmpdir." mkpath(joinpath(tmpdir,"mdf")) +mkpath(joinpath(tmpdir,"mdfim")) mkpath(joinpath(tmpdir,"positions")) mkpath(joinpath(tmpdir,"transferFunction"))