Ouch - that code caused rather a large amount of memory to be used for your
256 cubes - try this which allows you to break the job into feasible chunks
A 64cube takes 0.6sec, but 256cube takes 50sec per transform on my laptop
though (which is not good as the inner function of a registration!!!! :(
good luck anyhow
Martin
===============
function transform_image3d, image, rotation = rot, $
scale=scale,translate=translate, centre_rot=centre_rot, chunks=chunks,$
t3dmat=t3dmat
; transform a 3d image volume
s = size(image)
sx=s(1) & sy=s(2) & sz=s(3)
st = sx*sy*sz
imageT = image
; get transform matrix
if undefined(t3dmat) then begin
if undefined(rot) then rot =[0,0,0]
if undefined(centre_rot) then $
centre_rot=[(sx-1)/2.0,(sy-1)/2.0,(sz-1)/2.0]
if undefined(translate) then translate =[0,0,0]
if undefined(scale) then scale =[1,1,1]
t3d, /reset,trans= -centre_rot
t3d, rot=rot, trans= centre_rot + translate, scale=scale
t3dmat = !p.t
endif
; do transformation
if undefined(chunks) then chunks =1
iter = st/chunks
for i0 = 0L, (st-1), iter do begin
; account for possible odd last chunk
bufsize = iter < (st-i0)
;generate image coordinates
i = i0+lindgen(bufsize) ; temp array = vector indices
coords = [ [(i mod sx)],[((i / sx) mod (sy))],$
[(i /(sx*sy))],[replicate(1b, bufsize)]]
coords = (coords#t3dmat)
imageT[i0:i0+bufsize-1] = interpolate(image, coords(*,0), coords(*,1),$
coords(*,2), missing=0)
endfor
return, imageT
end
pro test,s, rot=rot, chunks=chunks
if undefined(s) then s = 256
if undefined(chunks) then chunks = 32
if undefined(rot) then rot = [20,10,45]
vol = rebin(bytscl(dist(s,s)),s,s,s, /sample)
t0= systime(1)
vol = transform_image3d(vol, rot = rot, chunks=chunks)
t1= systime(1)
print, "done in ", t1-t0, " sec"
; Display volume:
XVOLUME, vol
end
====================================
----------------------------------------
Martin Downing,
Clinical Research Physicist,
Orthopaedic RSA Research Centre,
Woodend Hospital, Aberdeen, AB15 6LS.
"Martin Downing" <martin.downing@ntlworld.com> wrote in message
news:Aabp7.24911$Pm5.5585206@news2-win.server.ntlworld.com.. .
> Hi Bob,
> Is this code any help or have I missed the point?
>
> ===========================================
> function transform_image3d, im, rotation = rot,
> scale=scale,translate=translate, centre_rot=centre_rot
> ; translate an image volume using interplote
> s = size(im)
> ; for clarity:
> sx=s(1)
> sy=s(2)
> sz=s(3)
> if undefined(rot) then rot =[0,0,0]
> if undefined(centre_rot) then centre_rot
=[(sx-1)/2.0,(sy-1)/2.0,(sz-1)/2.0]
[cut]
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
