From dba7ff03dc692024673c2afc8c18914b6df04c63 Mon Sep 17 00:00:00 2001 From: "Vladimir S. FONOV" Date: Wed, 3 Jan 2024 18:03:11 -0500 Subject: [PATCH] Fixed a bug in handling scans with oblique acquisition, added tests --- src/geo_transform.jl | 8 ++++---- src/minc2_io.jl | 12 ++++++------ src/minc_hl.jl | 4 ++-- test/runtests.jl | 45 ++++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 57 insertions(+), 12 deletions(-) diff --git a/src/geo_transform.jl b/src/geo_transform.jl index d4e43fa..ed700a7 100644 --- a/src/geo_transform.jl +++ b/src/geo_transform.jl @@ -324,7 +324,7 @@ Apply affine transform to a point p::SVector{3,T}; _whatever...)::SVector{3,T} where {T} - return (p' * tfm.rot)' + tfm.shift + return tfm.rot*p + tfm.shift end @@ -468,7 +468,7 @@ Apply affine transform to CartesianIndices tfm::AffineTransform{T}, p::CartesianIndex{3}; _whatever...)::SVector{3,T} where {T} - ( SVector{3,T}(p[1]-1.0, p[2]-1.0, p[3]-1.0)' * tfm.rot)' + tfm.shift + tfm.rot*SVector{3,T}(p[1]-1.0, p[2]-1.0, p[3]-1.0) + tfm.shift end @@ -487,8 +487,8 @@ function decompose(rot, shift) # remove scaling dir_cos = f.U * f.Vt - step = diag(rot * Base.inv(dir_cos)) - start = permutedims(permutedims(shift) * Base.inv(dir_cos)) + step = diag(Base.inv(dir_cos) * rot) + start = Base.inv(dir_cos)*shift return start, step, dir_cos end diff --git a/src/minc2_io.jl b/src/minc2_io.jl index 3cbb703..289e295 100644 --- a/src/minc2_io.jl +++ b/src/minc2_io.jl @@ -707,16 +707,16 @@ function voxel_to_world(hdr::MincHeader)::AffineTransform{Float64} for i=1:3 aa=findfirst(isequal(DIM(i)), hdr.axis) # HACK, assumes DIM_X=1,DIM_Y=2 etc if hdr.dir_cos_valid[aa] - rot[i,1:3] = hdr.dir_cos[aa,1:3] + rot[1:3,i] = hdr.dir_cos[aa,1:3] else rot[i,i] = 1 end scales[i,i] = hdr.step[aa] start[i] = hdr.start[aa] end - origin = permutedims(permutedims(start) * rot) + origin = rot*start - return AffineTransform(scales*rot, origin) + return AffineTransform(rot*scales, origin) end """ @@ -757,7 +757,7 @@ function create_header_from_v2w( hdr.start[i] = 0 hdr.step[i] = 1.0 hdr.dir_cos_valid[i] = false - hdr.dir_cos[i,:] = [0,0,0] + hdr.dir_cos[i,:] = [0.0,0.0,0.0] hdr.axis[i] = DIM_VEC elseif time_dim && (i-vector_dim)==4 # time dimension hdr.step[i] = time_step @@ -766,7 +766,7 @@ function create_header_from_v2w( else # spatial dimension hdr.start[i] = start[i-vector_dim] hdr.step[i] = step[i-vector_dim] - hdr.dir_cos[i,:] .= dir_cos[i-vector_dim,:] + hdr.dir_cos[i,:] .= dir_cos[:, i-vector_dim] hdr.dir_cos_valid[i] = true hdr.axis[i] = DIM(i-vector_dim) end @@ -785,7 +785,7 @@ function voxel_to_world(h::VolumeHandle,ijk::Vector{Float64})::Vector{Float64} @assert(length(ijk)==3) xyz::Vector{Float64}=zeros(Float64,3) @minc2_check minc2_simple.minc2_voxel_to_world(h.x[],ijk,xyz) - return ijk + return xyz end diff --git a/src/minc_hl.jl b/src/minc_hl.jl index 6254d30..a9bfefa 100644 --- a/src/minc_hl.jl +++ b/src/minc_hl.jl @@ -501,7 +501,7 @@ function resample_volume!(in_vol::AbstractArray{T,3}, fill=0.0, ftol=1.0/80, max_iter=10) where {C, T, I, XFM<:AnyTransform} - + in_vol_itp = extrapolate( interpolate( in_vol, interp), fill) @simd for c in CartesianIndices(out_vol) orig = transform_point(v2w, c ) @@ -578,7 +578,7 @@ function resample_volume!( end end end - + resample_volume!(in_vol.vol, out_vol.vol, out_vol.v2w, inv(in_vol.v2w), itfm; interp, ftol, max_iter, fill) diff --git a/test/runtests.jl b/test/runtests.jl index 59b410a..9e833fa 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -71,6 +71,7 @@ using StaticArrays vol,hdr,store_hdr = Minc2.read_minc_volume_std(i,Float64) @test hdr.dims == [30,40,10] @test size(vol) == (30,40,10) + # TODO: test hdr too if contains(i,"_byte_") || contains(i,"_ubyte_") @test mean(vol) ≈ 35.63256359 atol=0.5 @@ -260,5 +261,49 @@ end @test Minc2.transform_point(Minc2.inv(xfm),SA_F64[1, 2, 3]) ≈ SA_F64[2, 1, 3] end + # TODO: come up with test for nonlinear transforms +end + + + +@testset "Test Spatial coordinates correctness" begin + ### TODO: check that all trasformation works as expected + @testset "Low level interface" begin + h = Minc2.open_minc_file("input/t1_z+_float_cor.mnc") + @test Minc2.voxel_to_world(h, [0.0,0.0,0.0]) ≈ [-21.425182690867980995,59.125105562514228552,-8.3176836593094876093] atol = 1e-9 + @test Minc2.voxel_to_world(h, [10.0,10.0,10.0]) ≈ [-11.533050841537136222,69.446014720535146125,3.4986096455089885637] atol = 1e-9 + Minc2.close_minc_file(h) + end + + @testset "Oblique volume" begin + vol = Minc2.read_volume("input/t1_z+_float_cor.mnc", store=Float64) + v2w = Minc2.voxel_to_world(vol) + + @test Minc2.transform_point(v2w, SA_F64[0,0,0] ) ≈ SA_F64[-21.425182690867980995, 59.125105562514228552, -8.3176836593094876093] atol=1e-6 + @test Minc2.transform_point(v2w, SA_F64[10,10,10] ) ≈ SA_F64[-11.533050841537136222,69.446014720535146125,3.4986096455089885637] atol=1e-6 + + start, step, dir_cos = Minc2.decompose(v2w) + @test start ≈ SA_F64[-18, 60, -10] atol=1e-4 + @test step ≈ SA_F64[ 1, 1, 1.2] atol=1e-5 + + @test dir_cos[:,1] ≈ [0.998021217040696, 0.0523040113746069, -0.0348990075895229] atol=1e-6 + @test dir_cos[:,2] ≈ [-0.0517200153907152, 0.998509297133945, 0.0174420051903491] atol=1e-6 + @test dir_cos[:,3] ≈ [0.0357599860692531, -0.0156019939220494, 0.999238610734185] atol=1e-6 + + # calculate COM in voxel coordinates + arr = Minc2.array(vol) + v_com = @SVector [reduce((x,y)->x+arr[y]*(y[i]-1), collect(CartesianIndices(arr));init=0.0)/sum(arr) for i in 1:3] + @test v_com ≈ SA_F64[13.72049578, 12.72096456, 4.525147114] atol=1e-3 + + w_com = Minc2.transform_point(v2w, v_com) + @test w_com ≈ SA_F64[-8.195582241, 72.46002233, -3.148594157] atol=1e-3 + + end + + @testset "AffineTransform" begin + xfms = Minc2.load_transforms("input/linear.xfm") + + end + # TODO: come up with test for nonlinear transforms end \ No newline at end of file