Tubular regular positions (#63)

* First go on a tubular regular grid

* Tests for tubular regular grid

* Change file constructor

* Fix spelling

* Fix tests

* Try to fix Widows CI

* Create new dir

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

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

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

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