diff --git a/examples/resample.jl b/examples/resample.jl index b878323..4b96fa4 100644 --- a/examples/resample.jl +++ b/examples/resample.jl @@ -49,7 +49,13 @@ args = parse_commandline() if !isnothing(args["transform"]) - tfm=Minc2.load_transforms(args["transform"]) + + if endswith(args["transform"],".xfm") + tfm=Minc2.load_transforms(args["transform"]) + else # assume that it's a grid + tfm=Minc2.GridTransform(Minc2.read_volume(args["transform"],store=Float64)) + end + if args["invert"] tfm=Minc2.inv(tfm) end diff --git a/src/minc_hl.jl b/src/minc_hl.jl index 26d0518..4de4ab4 100644 --- a/src/minc_hl.jl +++ b/src/minc_hl.jl @@ -717,27 +717,27 @@ function calculate_jacobian!( # calculate scaling matrix from the voxel to world matrix # to compensate for the step size (makes no difference on 1x1x1 voxels) _,step,_ = decompose(out_v2w) - sc = det(diagm(Base.inv.(step))) + sc = prod(Base.inv.(step)) # First step: generate vector field of transformations vector_field = Array{T}(undef, 3, size(out_vol)...) - @inbounds @simd for c in CartesianIndices(size(vector_field)[2:end]) + @inbounds for c in CartesianIndices(size(vector_field)[2:end]) orig = Minc2.transform_point(out_v2w, c ) dst = Minc2.transform_point(tfm, orig; max_iter,ftol) - vector_field[:,c] .= dst #.- orig + vector_field[:,c] .= dst end # Second step: calculate jacobian determinant vector_field_itp = extrapolate( interpolate( vector_field, (NoInterp(),interp,interp,interp)), Flat()) - J=MMatrix{3,3,T}(undef) - @inbounds @simd for c in CartesianIndices(out_vol) - for i in 1:3 - J[:,i].=Interpolations.gradient( vector_field_itp, i, Tuple(c)...) - end + @inbounds for c in CartesianIndices(out_vol) + J = mapreduce(hcat, 1:3) do i + gradient(vector_field_itp, i, Tuple(c)...) + end :: SMatrix{3, 3, T, 9} + out_vol[c] = det( J ) * sc end