Fixed a bug in handling scans with oblique acquisition, added tests

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

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
 
 

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)
 

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