[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: The inverse of Extract_Slice
- Subject: Re: The inverse of Extract_Slice
- From: David Foster <foster(at)bial1.ucsd.edu>
- Date: Wed, 30 Sep 1998 13:14:16 -0700
- Newsgroups: comp.lang.idl-pvwave
- Organization: Univ. of Calif San Diego
- References: <3610E99A.CEB3FC6A@one.net>
- Xref: news.doit.wisc.edu comp.lang.idl-pvwave:12326
Larry Busse wrote:
>
> I'm starting a project which requires the "reconstruction" of a 3D data
> set from a series of 2D slices THROUGH the volume. The slices are NOT
> uniformly spaced (either by translation or rotation). I was thinking
> of hacking the IDL routine extract_slice to make a new function
> "insert_slice" to write the known data into the volume at the
> appropriate coordinates. Then do some sort of interpolation to fill in
> the holes.
>
> I have two questions:
> Has anyone already written an insert_slice routine?
> Any better ideas for approaching this problem??
Larry,
Your approach sounds very reasonable. I just want to caution you
that the EXTRACT_SLICE routine in essence assumes that you have
isotropic voxels; if you extract a slice from a volume having
voxels of size (1.0mm, 1.0mm, 4.0mm) you end up with a very
"stretched" image. Also, EXTRACT_SLICE.PRO does *NOT* handle
90-degree rotations correctly in general...you have to apply
rotations in a particular order when 90-degree rotations lead
to an "exchange of axes" situation.
I wrote a routine called RESLICE.PRO to get around these problems; it
borrows some of its code from EXTRACT_SLICE.PRO. Here it is, in case
it is useful, along with its documentation file. There's probably
some extra stuff in the routine that you can just throw out.
RESLICE() also corrects for the 90-degree rotation problem,
but for now it only supports one 90-degree rotation.
Also, the routine INTERP_3D() is a wrapper to a C function that
performs interpolation that treats the value 0 differently; for
your purposes you could just use INTERPOLATE() in its place.
Hope this is helpful.
Dave
;=============== RESLICE.PRO ========================================
; RESLICE.PRO 3-31-98 DSFoster
;
; Routine to reslice brain into image which is perpendicular to the Z
plane
; defined by the rotation angles supplied, and a specified distance "up"
; this Z plane. First computes the X,Y,Z coordinates for each pixel in
; matrix/section form by inverse rotation/translation, then reads pixels
; from volume supplied . Must supply the original midpoints and
; angles of rotation, the desired Z level (in mm and relative to angles
; given), the image to fill.
; Rotation angles passed are in degrees or, using keyword RAD, in
radians,
; and are relative to the matrix-section coordinate system.
;
; NOTE: We have learned that images that are "256x192" are in fact
; transformed into 256x256 format somewhere before they are stored
; onto tape. This means that the voxel is isotropic in the X and Y
; direction; therefore no correction needs to be made regarding
; "stretched" pixels within each image. However, the variables
; necessary to compute this correction, should it ever be necessary,
; are still read from a parameter file using READ_PARAMETER() (FOV_CM,
; XDIM and YDIM)
;
; If you want to get slices at angles relative to the "brain-centered"
; coordinate system, just start with rotation angles for the brain
relative
; to the image coordinate system, and then add your desired offsets.
;
; Use the INTERP keyword to specify which method of interpolation to
use:
; 0: Nearest neighbor sampling (the default)
; 1: Trilinear interpolation (IDL's INTERPOLATE)
; 2: Trilinear interp with zero-correction (our INTERP_3D from
; IMAGE_PROC_IDL.SO sharable module.
;
; This routine was adapted from RESLICE.FOR, a Fortran routine which
; resliced volumes on the PC.
;
; The algorithm for extracting the slice was borrowed from
; EXTRACT_SLICE.PRO, an IDL User's Library function, to speed things up.
;
; Modifications
;
; 3-22-94 DSF If OUT_VAL=-1 is flag to set out_val to zero. Flag
required
; since saying "...,OUT_VAL=0" disables keyword.
; 10-24-94 DSF Fix bug that led to strange behavior at the top and
bottom
; of the volume (z=0,z=zdim).
; 6-01-95 DSF Enable interpolation algorithm from sharable C module
; IMAGE_PROC_IDL.SO . Calls routine INTERP_3D() using
; call_external().
; DSF Remove CUBIC keyword, change SAMPLE keyword to INTERP.
; 9-18-96 DSF Add Z_BASE keyword to account for a coordinate
; system where Z begins at 1 (not 0). This was necessary
; to correct the inconsistency in MRORIENT.PRO that the
; X,Y coordinates were 0-based but Z, which is in section
; units, is 1-based. This led to all reslicing being off
; by one section; this was discovered when we noticed that
; the resliced 3D's were 2mm "behind" the FSE's.
; 12-18-96 DSF Use updated code from EXTRACT_SLICE.PRO (IDL User
Routine)
; for defining VOL_IND array, to speed things up.
; DSF Added AXIS_ROTATION to enable the exchange of axis
; necessary when creating one plane from another (eg.
coronal
; from sagittal). You cannot encorporate this exchange into
; the three offsets because once the axes are exchanged
; (say by a Y rotation of 90 degrees) then subsequent
rotations
; are with respect to the wrong axes. This exchange of axes
; must be done after the three offsets.
; 4-01-97 DSF Call wrapper routine INTERP_3D.PRO which calls
; C interpolation routines using CALL_EXTERNAL() (when
using
; zero-corrected trilinear interpolation INTERP=2).
; 12-15-97 DSF Allow only one element of AXIS_ROTATION to be nonzero.
; 3-23-98 DSF MAJOR BUG FIX: didn't treat Coronal -> Sagittal
transformation
; correctly. Use of AXIS_ROTATION was for Sagittal ->
Coronal
; transform, and produced correct results only because the
; X and Z rotation angles saved to the .pa were switched
; when the .pa was created by MrOrient (I misunderstood how
; Jim's routine in TRANSFORM_UTILS.PRO was working. The
angles
; computed for Sag scan are relative to the Sag
orientation,
; NOT the Cor orientation!). The AXIS_ROTATION keyword
failed
; for Cor -> Sag because the 90-degree Y rotation exchanges
; the X and Z angles unless the Y rotation is done FIRST.
;
; Now only a single angle of absolute value near 90 is
supported.
; When one occurs, this rotation is applied first so that
the
; remaining angles are applied with respect to the correct
axes.
;
; It is assumed that the 3D .pa's have been "corrected", by
; having their X and Z rotation angles exchanged, before
; using this routine.
;
;****************************************************************************
FUNCTION reslice, volume, x0,y0,z0, rotx,roty,rotz, zlocation,
param_file, $
RADIANS=radians, OUT_VAL=out_val, FOV=fov, $
INTERVAL=interval, INTERP=interp, Z_BASE=z_base
; Check for correct number of required arguments
if (keyword_set(FOV) and keyword_set(INTERVAL)) then begin
argexp = 8
endif else begin
argexp = 9
endelse
if (n_params() ne argexp) then begin
message, 'Incorrect number of arguments (' + strtrim(argexp,2) + $
' expected)', /continue
return, -1
endif
if ( keyword_set( Z_BASE ) ) then begin ; Is Z coord 0 => ? or
1 => ?
if ( z_base ne 1 ) then begin
message, 'Value for keyword Z_BASE must be 1 when defined',
/continue
return, -1
endif
endif else begin
z_base = 0
endelse
if (not keyword_set(INTERP)) then $
interp = 0 ; Default to sampling
; Get dimension and data-type information on volume, and create
; array that will hold resliced image (unless INTERP=2).
info = size(volume)
if (info(0) ne 3) then begin
message, "First argument must be three-dimensional array",/CONTINUE
return, -1
endif
xdim = info(1) ; Get dimensions of volume
ydim = info(2)
zdim = info(3)
type = info(4) ; Data type
; Check data types and create image to return, unless INTERP=2.
; Currently only supports BYTE or INT volumes. Define constants.
case (type) of
1: begin
if ( interp ne 2 ) then $
slice = bytarr(xdim,ydim,/nozero)
zero = 0B
one = 1B
end
2: begin
if ( interp ne 2 ) then $
slice = intarr(xdim,ydim,/nozero)
zero = 0
one = 1
end
else: begin
message, "Volume data type must be BYTE or INT", /CONTINUE
return, -1
end
endcase
if (keyword_set(OUT_VAL)) then begin ; Initialize image matrix
if (out_val(0) eq -1) then begin ; Flag for value of 0
sample_out_val = zero ; (one's will be removed
interp_out_val = one ; from image if
interpolation
if ( interp ne 2 ) then $ ; is used)
slice(*) = zero
endif else begin
case (type) of
1: begin
sample_out_val = byte( out_val(0) )
interp_out_val = byte( out_val(0) )
if ( interp ne 2 ) then $
slice(*) = byte( out_val(0) )
end
2: begin
sample_out_val = fix( out_val(0) )
interp_out_val = fix( out_val(0) )
if ( interp ne 2 ) then $
slice(*) = fix( out_val(0) )
end
endcase
endelse
endif else begin
if ( interp ne 2 ) then $
slice(*) = zero
endelse
; Convert angular measures to degrees if /RADIANS specified
if (keyword_set(RADIANS)) then begin
rot_x = (rotx * !PI)/180.0
rot_y = (roty * !PI)/180.0
rot_z = (rotz * !PI)/180.0
endif else begin
rot_x = rotx
rot_y = roty
rot_z = rotz
endelse
; Get required information from file containing image parameters.
; If keywords FOV and INTERVAL are set with appropriate values already
; then skip this.
if (keyword_set(FOV) and keyword_set(INTERVAL)) then begin
fov_cm = float(fov)
interval = float(interval)
endif else begin
interval = READ_PARAMETER(param_file, "SECTION_INTERVAL", /FLT)
fov_cm = READ_PARAMETER(param_file, "FIELD_OF_VIEW", /FLT)
info1 = size(interval)
info2 = size(fov_cm)
if (info1(1) eq 2 or info2(1) eq 2) then begin ; If INT then
error!
message, 'Unable to get volume parameters from file ' +
param_file, $
/continue
print, ' Missing parameters: SECTION_INTERVAL, FIELD_OF_VIEW'
return,-1
endif
endelse
; We want to allow for fully anisotropic pixels, so we must standardize
; our units of measure for the X,Y,Z coordinates. To do this, we assume
; the units for the X direction to be the standard, and then adjust
; the Y and Z accordingly. For the Y direction this means simply
; multiplying by the ratio of the X dimension of the original images to
; the Y dimension. For the Z direction, we calculate the number of
; X-pixels/section = (mm/section) / (mm/pixel) ( = INTERVAL/ZMM_PIX )
zmm_pix = (fov_cm/float(xdim)) * 10.0 ; mm/pixel (using X units)
zcorr = interval/zmm_pix ; Convert from Z-section =>
X-pixel units
ycorr = float(xdim)/float(ydim) ; Convert from Y-pixel to X-pixel
units
; Since Z location is in millimeters adjust to "image" units
(0,1,2,3...).
; This is the position of the desired image plane relative to the
rotation
; angles provided.
zlocsec = zlocation / interval
; Create one-dimensional matrix which will hold the X,Y,Z coordinates
; for the slice; these coordinates will be transformed using !P.T .
; (Borrowed from IDL User's Library routine EXTRACT_SLICE.PRO . Code
; commented out is from earlier version -- speed enhancement.)
im_size = long( xdim ) * ydim
; Statement replaces next 5 commented lines below, and is a speed
enhancement.
vol_ind = [ [reform( (findgen(xdim) # replicate(1.0, ydim)), im_size )],
$
[reform( (replicate(1.0, xdim) # findgen(ydim)), im_size )],
$
[replicate( 0.0, im_size )], $
[replicate( 1.0, im_size )] ]
; index = lindgen( im_size )
; vol_ind = fltarr((im_size), 4, /NOZERO)
; vol_ind(*,3) = 1.0
vol_ind(*,2) = zlocsec * zcorr ; Z => X pixel units
; vol_ind(*,1) = float(long(index / xdim))
; vol_ind(*,0) = float(index - (long(vol_ind(*,1)) * xdim))
; Compute transformation matrix using T3D function.
; First translate coordinates to "center" of VOL_IND 2D image.
; Now do rotation to get coordinates about matrix-section
; axes. If the rotations involve an exchange-of-axes (angle
; about 90 degrees) then perform this rotation *first* so
; subsequent rotations are about the correct axes.
; Then translate about the brain origin to give final matrix/section
; coordinates (account for whether the Z origin coordinate is 0-based
; or 1-based). Scale the Z dimension to account for section thickness.
save_pt = !P.T
T3D, /RESET ; Reset, translate to image
center
T3D, TRANSLATE = [ -(float(xdim-1L)/2.0), -(float(ydim-1L)/2.0), 0.0 ]
; Allow only one angle of rotation near (+/-) 90 degrees.
; If one exists, apply that exchange-of-axis rotation first so
; that the remaining angles are applied to the appropriate axes.
abs_rot = abs( [rot_x, rot_y, rot_z] )
obtuse = where(abs_rot ge 75.0)
if (n_elements(obtuse) gt 1) then begin
msg = 'Only one rotation > 90-degrees supported'
message, msg, /continue
return, -1
endif else begin
case (obtuse(0)) of
-1: begin ; All acute angles
T3D, ROTATE = [0.0, 0.0, rot_z]
T3D, ROTATE = [0.0, rot_y, 0.0]
T3D, ROTATE = [rot_x, 0.0, 0.0]
end
0: begin ; Obtuse X rotation
T3D, ROTATE = [rot_x, 0.0, 0.0]
T3D, ROTATE = [0.0, rot_y, 0.0]
T3D, ROTATE = [0.0, 0.0, rot_z]
end
1: begin ; Obtuse Y rotation
T3D, ROTATE = [0.0, rot_y, 0.0]
T3D, ROTATE = [0.0, 0.0, rot_z]
T3D, ROTATE = [rot_x, 0.0, 0.0]
end
2: begin ; Obtuse Z rotation
T3D, ROTATE = [0.0, 0.0, rot_z]
T3D, ROTATE = [0.0, rot_y, 0.0]
T3D, ROTATE = [rot_x, 0.0, 0.0]
end
endcase
endelse
T3D, SCALE = [1.0, 1.0, 1.0/zcorr] ; Adjust for slice
thickness
case ( z_base ) of ; Does Z start at 0 or 1?
0: T3D, TRANSLATE = float( [x0, y0, z0] )
1: T3D, TRANSLATE = float( [x0, y0, (z0 - 1.0)] )
endcase
vol_ind(*,*) = vol_ind(*,*) # !P.T(*,*) ; Transform coordinates
; Now fill the resliced image array
case (interp) of
0: begin
slice(*) = volume(0 > vol_ind(*,0) < (xdim-1), $
0 > vol_ind(*,1) < (ydim-1), $
0 > vol_ind(*,2) < (zdim-1))
if (n_elements(out_val) ne 0) then begin
out_v = where((((vol_ind(*,0) LT 0.0) OR $
(vol_ind(*,0) GE xdim)) OR $
((vol_ind(*,1) LT 0.0) OR $
(vol_ind(*,1) GE ydim))) OR $
((vol_ind(*,2) LT 0.0) OR $
(vol_ind(*,2) GE zdim)))
if (out_v(0) GE 0L) then $
slice(index(out_v)) = sample_out_val
endif
end
1: begin
if (Keyword_set(out_val)) then begin
slice(*) = $
interpolate(volume, vol_ind(*,0), vol_ind(*,1),
vol_ind(*,2), $
missing = interp_out_val)
if (out_val(0) eq -1) then begin ; Flag for 0 so remove
1's
zeroes = where(slice eq 1)
slice(zeroes) = zero
endif
endif else begin
slice(*) = interpolate(volume, vol_ind(*,0), vol_ind(*,1), $
vol_ind(*,2))
endelse
end
2: begin
; Call wrapper function that uses CALL_EXTERNAL() to call
; routines from C module IMAGE_PROC_IDL.SO .
status = 1
slice = interp_3d( volume, vol_ind(*,0), vol_ind(*,1), $
vol_ind(*,2), /zero_corrected, status=status )
if ( status ne 0 ) then begin
message, 'Internal error within IMAGE_PROC_IDL.SO',
/continue
return, -1
endif else if (slice(0) lt 0) then begin
message, 'Error: INTERP_3D.PRO, module IMAGE_PROC_IDL.SO', $
/CONTINUE
print, ' Invalid argument or module not found.'
return, -1
endif
end
endcase
!P.T = save_pt
return, slice
END
;================== END OF RESLICE.PRO ==============================
=================== RESLICE.DOC =====================================
NAME:
RESLICE
PURPOSE:
This function returns a 2-D planar slice extracted from
3-D volumetric data. The slicing plane may be oriented at
any angle, and may pass through any desired location in the
volume. Nearest neighbor sampling is the default sampling
technique. The user may specify that trilinear interpolation
be used. The field-of-view (FOV) and distance between images
for the imaging sequence are required, and may be specified
using keywords or contained in a configuration file
(specified as an optional parameter).
CALLING SEQUENCE:
Slice = RESLICE(Vol, X_center, Y_center, Z_center, $
X_rot, Y_rot, Z_rot, Z_location [, param_file])
INPUTS:
Vol: The three dimensional volume of data to slice.
Data type : BYTE or INT.
X_center: The X coordinate (index) of the "origin" of the volume.
Data type : Any scalar numeric value (usually Long).
Y_center: The Y coordinate (index) of the "origin" of the volume.
Data type : Any scalar numeric value (usually Long).
Z_center: The Z coordinate (index) of the "origin" of the volume.
Data type : Any scalar numeric value (usually Long).
X_rot: The orientation (X rotation) of the slicing plane.
Before transformations, the slicing plane is parallel
to the X-Y plane. The slicing plane transformations
are performed in the following order :
1. Rotate Z_rot degrees about the Z axis.
2. Rotate Y_rot degrees about the Y axis.
3. Rotate X_rot degrees about the X axis.
4. Perform any 90-degree rotations that result
in an exchange of axes.
5. Translate the center of the plane to
X_center, Y_center, Z_center.
6. Move the image along the "new" Z axis
Z_LOCATION units.
Data type : Float.
Y_rot: The orientation (Y rotation) of the slicing plane.
Data type : Float.
Z_rot: The orientation (Z rotation) of the slicing plane.
Data type : Float.
Z_location: The location along the "new" Z direction (relative to
the rotation angles provided) from the origin where the
center of the image is located. If zero, the image is
centered at the origin; otherwise, it is centered a
distance out from this origin in the direction defined
by X_rot, Y_rot and Z_rot. The units are millimeters.
Data type: float.
Param_file: A file containing information regarding the field-of-
view of the images of the volume and the distance
between images. Data type: String. See the example
below.
KEYWORD PARAMETERS:
FOV: Set this and the INTERVAL keyword to avoid needing
a parameter file. This specifies the field-of-view
of the original images in the volume, which is the
width of the volume in the X dimension. Units are cm.
Data type: Float.
INTERP: Specifies the type of sampling to use to compute the
resliced images. Possible values are:
0: Nearest neighbor sampling (the default)
1: Trilinear interpolation using IDL's INTERPOLATE()
2: Trilinear interpolation which does not use the
value zero to compute interpolate values. Avoids
aliases near borders with background voxels having
value zero. Requires INTERP_3D() from sharable
C library module IMAGE_PROC_IDL.SO .
INTERVAL: The distance between sections. Set this and FOV to
avoid using a parameter file. Units are mm.
Data type: Float.
OUT_VAL: If OUT_VAL is set, then the portions of the returned
slice that lie outside the original volume are set to
the value passed to OUT_VAL. To set outlying pixels
to 0 set OUT_VAL=-1 (since OUT_VAL=0 "unsets" the
keyword!). Data type : Any scalar numeric value
(usually the same type as Vol).
RADIANS: Set this keyword to a non-zero value to indicate that
X_rot, Y_rot, and Z_rot are in radians. The default
is degrees. Data type : Int.
OUTPUTS:
This function returns the planar slice as a two dimensional
array with the same data type as Vol. The dimensions of the
returned array are the same as the X- and Y-dimensions of the
volume.
EXAMPLE:
Display an oblique slice through volumetric data.
; Create a volume from an image sequence
ret = MAKEVOL('/spare/im/892413*.mr', vol)
; Extract and display a slice. Outlying pixels to zero.
slice = RESLICE(vol, 128.0, 128.0, 17.0, 3.0, -20.5, 0.0, $
30.0, param_file, OUT_VAL=-1)
TVSCL, REBIN(slice, 400, 400)
; Command when not using parameter file:
slice = RESLICE(vol, 128.0, 128.0, 17.0, 3.0, -20.5, 0.0, $
30.0, OUT_VAL=-1, FOV=24.0, INTERVAL=4.0)
; Create a coronal image from a sagittal volume:
slice = RESLICE(vol, 128.0, 128.0, 62.0, 0.0, 0.0, 0.0, $
0.0, OUT_VAL=-1, FOV=24.0, INTERVAL=1.2, $
AXIS_ROTATION=[0.0, 90.0, 0.0]
EXAMPLE PARAMETER FILE:
SECTION_INTERVAL: 4.0;
FIELD_OF_VIEW: 24.0;
The parameter file must be in a format which can be read using
READ_PARAMETER routine.
See Also: MAKEVOL, READ_PARAMETER
--
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
David S. Foster Univ. of California, San Diego
Programmer/Analyst Brain Image Analysis Laboratory
foster@bial1.ucsd.edu Department of Psychiatry
(619) 622-5892 8950 Via La Jolla Drive, Suite 2240
La Jolla, CA 92037
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~