Maybe slightly off-topic, but it looks like you're using temp as a (1, 3)
Array, rather than a (3, 3).
On Wednesday, November 19, 2014 10:09:49 AM UTC+9, Ryan Young wrote:
>
> 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
>
