Hello,
I am new to julia and am working on some neuroimaging analysis. I have a
function that calculates the Jacobian of the deformation field at each
point in the deformation field. A qr decomposition is then done on the
Jacobian. The principle eigenvector of the strain portion of the QR
decomposition is rotated, stored in an array and returned.
I was not going to parallelize this function as I think it is easier to do
a pmap of this process over all of the samples. Currently, this is taking a
few minutes to run on one large deformation field (225,425,225,3). All the
time is being taken up in the last three lines of the for loop. Any advice
is greatly appreciated.
Thanks
Ryan
function calcJac(deformation::Array{Float32,4})
temp=Array(Float32,3,3)
r=Array(Float32,3,3)
s=Array(Float32,3,3)
ev=Array(Float32,3,3)
i,j,k=size(deformation)
final=Array(Float32,i,j,k,3)
for x in 3:i-3
for y in 3:j-3
for z in 3:k-3
#compute frist row of jacobain
temp[1,1]=(-deformation[x+2,y,z,1]+8*deformation[x+1,y,z,1]-8*deformation[x-1,y,z,1]+deformation[x-2,y,z,1])/12.0
temp[1,2]=(-deformation[x+2,y,z,2]+8*deformation[x+1,y,z,2]-8*deformation[x-1,y,z,2]+deformation[x-2,y,z,2])/12.0
temp[1,3]=(-deformation[x+2,y,z,3]+8*deformation[x+1,y,z,3]-8*deformation[x-1,y,z,3]+deformation[x-2,y,z,3])/12.0
#compute 2nd row of jacobian
temp[1,1]=(-deformation[x,y+2,z,1]+8*deformation[x,y+1,z,1]-8*deformation[x,y-1,z,1]+deformation[x,y-2,z,1])/12.0
temp[1,2]=(-deformation[x,y+2,z,2]+8*deformation[x,y+1,z,2]-8*deformation[x,y-1,z,2]+deformation[x,y-2,z,2])/12.0
temp[1,3]=(-deformation[x,y+2,z,3]+8*deformation[x,y+1,z,3]-8*deformation[x,y-1,z,3]+deformation[x,y-2,z,3])/12.0
#compute third row of jacobian
temp[1,1]=(-deformation[x,y,z+2,1]+8*deformation[x,y,z+1,1]-8*deformation[x,y,z-1,1]+deformation[x,y,z-2,1])/12.0
temp[1,2]=(-deformation[x,y,z+2,2]+8*deformation[x,y,z+1,2]-8*deformation[x,y,z-1,2]+deformation[x,y,z-2,2])/12.0
temp[1,3]=(-deformation[x,y,z+2,3]+8*deformation[x,y,z+1,3]-8*deformation[x,y,z-1,3]+deformation[x,y,z-2,3])/12.0
r,s=qr(temp)
ev=eigvecs(s)
final[x,y,z,:]=r*ev[:,1]
end
end
end
return final
end
