P22 virion structure, EMD-1220 rendering notes

P22 virion structure, EMD_1220 rendering notes

"For a physicist, it is indeed a great joy to learn how we can use beautiful mathematics to understand the real world." In Stephen Weinberg's "Without God" (9/25/08, The New York Review of Books)

Data source
Processing notes
Test rendering notes
Test renderings of the P22 virion EMD-1220 data set

Programs and code
Space Software
FSL
MATLAB code
Other structural questions
Homogenous coordinate transform matrices

2008-09-08 I just took a look at the EMD_1220 volume data set and it is very nice. My goal is to fully segment and color render, with cutaways, this volume data.

2008-09-16 The EMD_1220 volume data set was reconstructed from cryo-electron microscopy images. It represents, roughly speaking, a three dimensional estimate of electron densities of an average P22 virion particle.

These notes outline post-processing that I performed in order to do renderings of the structure and sub-structures represented in the volume.

2008-09-18 The post-processing has morphed into a somewhat separate project: make an average over axial n-fold rotations, where n is point dependent and point-wise derived. I considered this approach to "cleaning" the data (averaging out noise that doesn't have a rotational symmetry) while rendering EMD-1222 . I've tested the idea (see Processing notes, averaged rotations) and now can write the MATLAB code that implements it with precision.

Data source

2008-08-03 mcman12 wrote a comment on a short rendered animation loop of mine posted on YouTube, Virus structure, which was rendered from the EMD-1222 data set .

"for a higher resolution (better) reconstruction check out EMD-1220"

This EMD-1220 data set (a cryo-EM assymetric reconstruction, 3-D grayscale) is available (accessed 2008-09-08) at http://www.ebi.ac.uk/msd-srv/emsearch/atlas/1220_downloads.html, in the European Bioinformatics Institute 's Protein Data Bank Europe (PDBe). Also see the EM Data Bank , which provides three-dimensional density maps for non-viruses.

From the data set web page description:

The structure of an infectious P22 virion shows the signal for headful DNA packaging
Gabriel C. Lander, Liang Tang, Sherwood R. Casjens, Eddie B. Gilcrease, Peter Prevelige, Anton Poliakov, Clinton S. Potter, Bridget Carragher, and John E. Johnson:
The structure of an infectious p22 virion shows the signal for headful DNA packaging J Mol Biol (2006) 312, pp. 1791-1795 [PubMed entry 16709746]
3-Jun-2006
singleParticle, 17 Å (resolution determined by FSC at 0.5 cut-off

Processing notes

2008-09-08

- Convert from CCP4 map (.map/.xml) to Space volume (.vol) file format.
2008-09-08 MATLAB program convert_map_volume.m does this (see attached). First read the .xml file (as text) to find the volume dimensions (the EMD_1220.map is 256 x 256 x 256). A hardcoded black and white point must be chosen for each data set conversion. A histogram of final values is displayed that can help with this.
I'm not sure I captured the full dynamic range, particularly at the lowest values. Add gamma conversion control? Or redo with 32-bit NIFTI format?
MATLAB command:
>> volData = convert_map_volume( 'emd_1220.map', 'emd_1220_8-bit', [ 256 256 256 ] );

2008-09-14 Rewrote convert_map_volume.m for NIFTI 32-bit output.
MATLAB command:
>> volData = convert_map_volume( 'emd_1220.map', 'emd_1220', [ 256 256 256 ] );

- Allign the virion's raw data longitudinal axis with the volume's z-axis accurately.
- 2008-09-16 The approach I chose was to find a best fit longitudinal axis took advantage of the high radial 5-fold symmetry of the capsid. The capsid was first roughly alligned with the volume's z-axis (within about an angle of 1 px./80 px., just an eyeball resampling at this orientation in Space), and roughly segmented (a binary mask that included almost all the capsid, but not the interior or tail structures). The masked volume (just the capsid) was then averaged with its four roughly equivalent 2*pi/5 rotations. This average volume is axially symmetric, and served as a target for coregistration of the raw data.
- Rough (low resolution) segmentation of capsule head and contents (DNA).
2008-09-08 Cubic resampled (x2) and linearly resampled in an approximately axial orientation. Binary mask for substructures (capsid, inner head, outer head, DNA shells 1, 2, 3, and 4) using 2 and 3-D fill on the full binay image. The capsid has the primary orientation information, at least for any 5-fold symmetry structure. All using Space Software.
- 2008-09-13 After masking the high-value 8-bit volume with the rough capsid binary segmentation (using Space), FSL was used to make 5 copies, at n*2*pi/5 rotations.
- 2*pi/5 rotation transform matrix
- MATLAB calculation for xt and yt coordinate transformation matrix:
>> xt = xc + ( sqrt( xc*xc + yc*yc ) )*sin( (2*pi)/5 - atan( xc/yc ) )
>> yt = yc - ( sqrt( xc*xc + yc*yc ) )*cos( (2*pi)/5 - atan( xc/yc ) )   ,
where xc = 165 and yc = 163 are the center coordinates with respect to the corner origin.
- The full transformation matrix, in homogenous coordinates is (see homogeneous transformation matrix for details):
0.309017 -0.951057   0 269.0934
0.951057   0.309017   0   -44.2941
0              0                1      0
0                0                0      1

- Apply rotation transform matrix four times to the masked capsid volume.
FLIRT | Utils | Apply FLIRT transforms (ApplyXFM)
or
flirt -applyxfm -init -out -ref
- These 5 volumes were averaged to create capsid_mean.nii:
- FSL merge of the five capsid rotations to 4-D volume and collapse across the "time" (n=5) dimension:
fslmerge -t capsid_avg.nii.gz emd_1220_rot_0 emd_1220_rot_72 emd_1220_rot_144 emd_1220_rot_216 emd_1220_rot_288
fslmaths capsid_avg -Tmean capsid_mean
Note that this took a long time, apparently due to memory limitations. I think this time could be reduced using fslmaths -add, but would require a conversion to a floating point format, and adding each one sequentially. Adding in the 8-bit format caused value wrap-around (8-bit overflow).
I think that adding sequentially will be required for the the alligned image rotation averages, to maintain the bit depth without memory limitations.
[To Do: Screenshot of capsid and capsid_mean section.]
- An assymetric capsid slice (left, high contrast image of an angled slice, linearly interpolation resampling) compared with the coregisitration target (right, a 5-fold mean image, an orthogonal slice):Comparison of the nearly alligned roughly segmented capsid (left) and the 5-fold rotational average (right). At top is an axial slice through the vertex at top of capsid (magnified x2 relative to bottom images), where the mis-allignment is most pronounced.
- 2008-09-15 Coregister over-sampled original volume data to axially symmetric target.
- Raw data was cropped (to limit memory footprint) and over-sampled (tri-cubic) by a factor of x2, in Space Software.
The raw 256x256x256x4 byte = 64 MB volume is too large to resample x2³. So need to crop (in Space after Volume | Volume Info | Center, exact), both this raw data and the capsid weighting volume, emd_1220_8-bit_high-125-230_x2_capsid_weight.nii.gz.

The raw data crop at native resolution:

x: [ 41 : 212 ] (even width s.t. the resampled dimensions will be odd)

y: [ 35 : 218 ]

z: [ 39 : 256 ]
- A binary coregistration weighting volume was created (emd_1220_8-bit_high-125-230_x2_capsid_weight_cropped.nii) by segmenting the capsid at the original orientation, in a procedure identical to that above (rough segmentation of capsule head).
This was using a high values (125-230) 8-bit volume, but at the original orientation. This tilted orientation made it more difficult to segment the tail/capsid boundary, but it doesn't need to be perfect -- even a majority portion of the capsid would coregister well.
Weighting origin tweaked to exactly overlap raw data:
Note the value wrap flaw on the high luminance "needle" center at bottom-left. This might be due to a flaw in the Space tri-cubic resampling algorithm, or in the MATLAB conversion code. [ 2008-09-15 I can minimize the effect of this, but it seems to be a problem with Space's storage or interpretation of NIFTI range information.]
- FSL coregistration:
Target: capsid_mean.nii
Object: emd_1220_cropped_x2.nii
Weighting of object: emd_1220_8-bit_high-125-230_x2_capsid_weight_cropped.nii
Coregistered output: emd_1220_cropped_x2_alligned.nii
Normalized correlation method. Straight correlation method doesn't work, because the raw values (32-bit) are a very different range than the mean capsid (8-bit).
Tri-linear resampling.
flirt -in /emd_1220_cropped_x2.nii.gz -ref capsid_mean.nii.gz -out /emd_1220_cropped_x2_alligned.nii.gz -omat /emd_1220_cropped_x2_alligned.mat -bins 256 -cost normcorr -searchrx 25 35 -searchry -5 5 -searchrz -5 5 -dof 9 -inweight /emd_1220_8-bit_high-125-230_x2_capsid_weight_cropped.nii.gz -interp trilinear
- Final alligned volume, emd_1220_cropped_x2_alligned.nii, is 32-bit, 333 x 339 x 435 vx. (after post-cropping), origin at (164, 169, 217), 188 MB:
EMD-1220, orthogonal slices of alligned volume, high contrast.
- Averaged 5-fold rotation of capsid/DNA and 6-fold rotation of tail structures.

2008-09-16 Write a specialized MATLAB program for radial symmetries analysis and n-fold rotational averaging. Started rotation_average_volume.m.
2008-09-17 Here's a proof of concept image using a test version of rotation_average_image.m. The concentric DNA rings near the axial capsid center are nearly circular. By rotational averaging the central radial regularities are highlighted:
The original image (left) and an average with its pi/2 radian rotation about the center (right). The ring definition near the center is enhanced. Note that the capsid (pink at left) does not have a 2-fold symmetry on this plane, so its rotation has only a small overlap (pink at right).

DNA rings averaged over 5-fold (left) and 90-fold (right) rotations. The 5-fold rotation matches the capsid symmetry. Nine DNA rings are clear, but the central features are likely artifacts (noise, or different rotational symmetry).

2008-09-21 Tentative procedure to find center accurately and average over 2*pi/5 and 2*pi/6 rotations. I've eperimented with this procedure on 2-D slices, but only minor modification need to be made to do 3-D serial slices:
- Remove background (set background to marker symbol, a low value). See remove_background_peak_image for 2-D demo. [To Do: Don't remove background values from central features.]
- For 3-D volume, use Space 3-D fill tool to accomplish this. Make an 8-bit mask and apply to 32-bit volume.
- 2008-10-01 Did this process on emd_1220_cropped_x2_alligned.nii. Note that it is not needed for the rotational averaging algorithms (below). But for finding the center in single slices, it worked well to remove the contribution from noise, only taking into account the standard deviation of histogram values that were not on the peripheral noise.
- Exterior mask Overlay | New Blank Overlay | 8-bit with Edit | 3-D Fill (range = 1). Highest contrast on original, at a black point that best segregated interior an exterior. Fill in all exterior to the capsid (except for discontinuous blobs) without much attention to the exact background level. At lower tail planes the background values were contiguous, so used a wide pen tool to manually segment (circle) ROIs.
- Interior mask: Close original, make a new blank overlay to the exterior mask. Fill interior. There is some bleed into random noise contiguous with the interior, but it is not important for rotational averaging as at least n-1 points must be non-zero for an n-fold rotation.
- Overlay | Mask original (16-bit) with interior mask (8-bit, on top), and Overlay | Merge. Save as NIFTI because there is some BP/WP/Min/Max (Colors | Transforms and Histograms) interpretation bug with Spaceafter re-opening.
- Sub-pixel approximation of center by applying rotations (5 and/or 6-fold) an minimizing standard deviation. Iterative hill climbing.
>>rotation_average_image( 'test_remove_background-zeroback.png', [169 169], 5, true, 'test_' );
[2008-09-22 To Do: Ignore zero values (background), and require a minimum number of low deviation points. Ignore background point in net deviation calculation. Set map points to zero if not enough forground points on rotation.]
[2008-09-22 To Do: Accept an n-fold rotation map and only calculate deviation of rotation points for each corresponding n-fold rotation.]
- ~~Use center estimate to find best n-fold rotation by finding minimum standard deviation of each rotated point with each n-fold rotation.~~
>>rotation_nfold_image( 'test_remove_background-zeroback.png', [169.2383 170.0625], [5 6 10 12 15 18 24 25 30], true, '' );
[2008-09-22 T~~o Do: Jitter the center point for each point rotation. Save the maximal point for the n-fold rotation that was most stable with jitter. Also save the jitter maps.~~]
- 2008-09-24 Experimented with this a lot. This point by point method (evaluating best n-fold fit at each point) is not the way to go -- to noisy and not easy to integrate neighborhood and global n-fold information. See more simple approach below.
- ~~Repeat center finding, using the new n-fold rotation map. Different centers for different n?~~
- Rough 2 part segmentation of alligned volume by symmetry (5-fold for capsid and DNA, 6 or 12-fold for tail). Include some overlap at intersections. Find best fit n-fold rotation and center for all of one segmented slice, checking for other possible symmetries (circular symmetry, nearby odd symmetries). Repeat on all slices. Trim and blend two (or more) symmetry averaged volumes.
High contrast original (left) and composite after 5-fold rotational averaging of capsid (white on right) and 6-fold averaging of tail (pink on right). Some (unknown) parts of the 12-fold symmetric "bumps" at the tail periphery are artifacts of the rotational averaging. Ameliorate by excluding oulliers on annular rings before averaging.
- Used rotational center finding and averaging on several planes of the 16-bit alligned data with exterior noise zeroed (see Exterior mask above). It is beautiful, both 5 and 6-fold components on the critical boundary planes. The calculated centers were dead on for all capsid only planes, and only a possible +-.25 pixel shift of tail only planes, with little or no axis skewing.
- Currently used this program, with only one plane of interest and one symmetry hardcoded, to center find and rotationally average:
>> rotation_average_volume( 'emd_1220_cropped_x2_alligned_background_masked.nii', [170 165] );
[To Do: Example average images, derived center estimates.]
- 2008-10-01 Rotational symmetry decomposition.
- Extract an annulus that has 5 and 6-fold components from a single plane (radial components_image.m). Do component analysis manually and code annulus_components.m
- Average over m, n, ... -fold increments, excluding outliers. Hold in array with dimensions radius, angle, [1, m, n, ...].
- Filter high frequencies ( > m*n cycles/annulus) and normalize.
- Least squares fit of these models to the original. Find residual.
- 2008-10-06 Testing radial_components_image.m
- 2008-10-15
>> fPwr = radial_components_volume( 'emd_1220_cropped_x2_alligned.nii', [170 165] );
calls:     fpwower = radial_components_image( imPlane, xyCenter, false, num2str( iP ) );
% fPwr = radial_components_volume( strInputFilename, xyCenter )
% ------------------------------------------------------------------------------
%
%   Fourier decomposition of annuli about the center of a selected serial
%   grayscale z-planes of a NIFTI volume.
% Hardcoded information:

bShow = false;   % [true, false] If true, displays each raw slice.
bClear = false;   % [true, false] Unload volume from memory while analyzing each plane.

iPlanes = 1:435    % [ n : m ] z-coordinates of planes to be analyzed.
rRes = 2;

% fPwr = radial_components_image( strInputFilename, xyCenter, bShow, strOutputFilePrefix )
% ------------------------------------------------------------------------------
%
%   Fourier decomposition of annuli about the center of a grayscale image.
%
%   This version is specialized for analysis of the 5 and 6-fold symmetry
%   of the P22 virion (particularly the EMD-1220 data set), but it could be
%   modified for other annular symmetries.
%%%%%%%%%%%%%%%%%%%%%%%%
% Hardcoded information:

    bWrite = true;
    bDoReconstruct = true;
    rRes = 2;

% To Do: Calculate rmax from the given center and image size, use this
% parameter as a second hard maximum.
    rmax = 163;

    % Annular sampling rate; the number of sample points on each annulus.
    % Lower integer factors have some computational speed advantage.

    % This is longer than needed, and takes up considerable memory.
    np = 2*2*3*3*5*7*11     % = 13860, excludes 8 as a factor

%
%%%%%%
%     fPwr           - total Fourier power of annuli at each radius.
%                       1 - 1 cycle/annulus
%                       2 - 2*n, excluding mod 5 and mod 6
%                       3 - 3*n, excluding mod 5 and mod 6
%                       4 - 4*n, excluding mod 5 and mod 6
%                       5 - 5*n cycles/annulus, excluding 30, 60, ...
%                       6 - 6*n cycles/annulus, excluding 30, 60, ...
%                       7 - 7*n, excluding mod 5 and mod 6
%                       8 - all nonzero and not in 5 or 6-fold
%                       7 - [not used, empty]
%                       10 - zeroth (DC, or constant) component
%                       11 - 11*n, excluding mod 5 and mod 6
%                       12 - 12
%                       13 - 13*n, excluding mod 5 and mod 6
%                       14 - 2
%                       15 - 3
%                       16 - 4
%                       17 - 6
%                       18 - 8
%                       19 - 30*n
%                       20 - 60*n
>> serial_mat_to_nii( 'test_', 'Test_out.nii' )
% serial_mat_to_nii( strInputFilestem, strOuputFilename )
% ------------------------------------------------------------------------------
%
%   Read a series of 2-D matrices from .mat format into a 3-D volume and
%   save as a NIFTI volume.
%
%   Currently it is used for assembling the planes output by
%   radial_components_image.m.
%   This version is specific to a set of .mat files that
%   contained an image in a structure field .im56fold.

% Hardcoded information:

bShow   = false;   % [true, false] If true, displays each raw slice.
iPlanes = 1:435 % [ n : m ] z-coordinates of planes to be analyzed.
Before (left) and after (right) rotational symmetry filtering. High contrast to emphasize DNA packing structures. The 5 and 6-fold filtering parameters are just a rough guess here. There are some flaws, for example the blobs at far right are artifacts because I neglected to include overtones of 5-fold symmetries that are also overtones of 6-fold symetries, even though there are negligable 6-fold components at that location.
- Spectral powers above a threshold, for several sets of frequencies using spectral_segmentation_EMD_1220.m:
MATLAB command:
>> load emd-1220_RPwr
>> spectral_segmentation_EMD_1220( 'Cross_section_half.png', fPwr )
where emd-1220_RPwr.mat was generated by radial_components_volume.m (see above), containing the fPwr 3-D array.
%%%%%%%%%%%%%%%%%%%
% Hardcoded information:

bShow   = false;   % [true, false] If true, displays each raw slice.
bWrite = true;    % [true, false] If true, write segmentation image.
strOutFilestem = 'test_seg_fundamental_6'

iChannels = [ 21:33 ]
%   iChannels = [ 32 ]
%ndev = 5.0;     % Segmentation threshold, n standard deviations from median cut.
ndev = 6.0;     % Segmentation threshold, n standard deviations from median cut.

Reperesentative cross section (left) compared with spectral segmentation (right).
Red: f = 5 or 10 cycles/annulus
Green and cyan: f = 6 cycles/annulus
Blue: f = 12 cycles/annulus
Magenta: f = 5 and 12 cycles/annulus
Yellow: f = 6 and 12 cycles/annulus
White: Other frequencies, primarily 7 and 13 cycles/annulus, and overlap of 5, 6 and 12 cycles/annulus. The bottom right corner is an artifact due to constant (zero) filling where there was no image data.
- 2008-10-20 For second itteration of spectral density segmentation:
- Use radial_components_volume.m to generate the f = 0 (DC, or constant) component volume.
- Output in 16-bit, format instead of 32-bit, for convenience.
- Add f = 14, 17, 18 and 19 components to spectral powers array.
- Pre-filter super-Nyquist frequencies?
- Use a map of significant spectral power for generating component volumes: imPwrSigMap( r, z, f_a:f_z), where ( r_n, z_n, f_a:f_z )_n <= 1
2008-10-24 Done and tested. Need to add logic for overlapping regions -- what ratio of overtones to use.
- Test generation of another frequency and harmonics volume. fa = n*7, fa = n*13, or fa = 1 (fundamental frequency only).
- 2008-10-21 Made imPwrSigMap. Find clusters with at least one value == 1 and remove all other values (set to zero).
- Fill center point of reconstructions properly
- 2008-10-24 It is apparent that a large fraction of the non-5-6-12-fold power spikes (1-7-13-fold, larger white blobs in the spectral segmentation map above) are due to mis-estimation of the center(s) of rotation. Remember that the original (single) estimate was based on the capsid segmentation, all 5-fold.
Re-estimate the centers for the 5-fold, 6-12-fold and central (about circularly symmetric) regions independently. Use hill climbing and spectral decomposition. For which center is there the least residual power?
- 2008-11-03 Currently testing rotation_center_volume/image.m for center estimation by n-fold spectral power peak , about the longitudinal axis.
rotation_center_image.m calls radial_components_image.mwith the hardcoded parameter:
nCriterion = 2
- Full final segmentation and mask.
- Separate thresholds for some or all structures, as separate 8-bit images.

- Wedge mask. Apply to elements.
2008-09-10 Tested a 90 degree (volume corner) wedge, it works well. In Space Software, simpy create an all white volume and adjust its x/y origin (Volume|Volume Info) until its corner is near the origin of the virion volume object. Apply Overlays|Mask with this Overlay and Overlays|Merge All.
- Render combinations.

Test rendering notes

2008-09-10 Test renderings of the P22 virion EMD-1220 data set:

Also see larger version of the same in Volume Renderings gallery.

The first segmentations (see above "Rough (low resolution) segmentation" and "binary coregistration weighting volume") were performed on an 8-bit volume that was extracted using convert_map_volume.m hardcoded with:

nBlackPoint = 125;
nWhitePoint = 230;
emd_1220_8-bit_high-125-230.vol.gz

This dynamic range was very good for that purpose. But I went back and looked at the very low values to see if I was missing dim structure. These parameters, after on-line contrast/brightness selection, worked better for rendering of some fine structure:

nBlackPoint = 25;
nWhitePoint = 200;
emd_1220_8-bit_low-25-200.vol.gz

Programs and code

Space Software: My free Windows GUI for display, navigation, rendering, editing, and measurement of 3-D data sets.
Space was used for some of the processing (resampling, value scaling, etc), most of the image slices, and all of the rendering.

FMRIB Software Library (FSL): Free linux based analysis tools intended for for FMRI, MRI and DTI brain imaging data, but the coregistration and volume math tools are suitable for use with other volume data. GUI and command line tools. Used for coregistration and some rotation transforms and averaging.

2008-09-12
To avoid writing code that allows finding the axis of rotation accurately (see above, iterated least squares self-coregistration fit about capsule's 5-fold symmetric axis, and head's six-fold rotational axis), I decided to try a well tested coregistration program. The FLIRT tool in the Linux FSL toolbox does volume coregistration with a variety of options, and also has functions that can resample, transform and average volumes.
As I'm using a Windows operating system, I needed to install VMware Player to run the FSL virtual machine. This was a huge hassle, taking a lot of time to get everything working about right.
Space Software still has some compatibility issues (and a bug or two) with respect to the NIFTI volume format. In particular, the first value of the dim[8] header item gets written as 1 (short integer), while it should be 3 (the number of dimensions of the data set. The work-around is to manually set this value in the Space Volume | Volume Information | View/Edit Header dialog box. Also, the header length appears to by off by 4 bytes, but this didn't crash FSLView (v3.0) or FLIRT
4x4 transforms can by applied to volumes using the FLIRT utility (flirt -applyxfm, -init and -out). This transform program can be accessed with a FLIRT GUI using the Util | Apply FLIRT Transform (ApplyXFM) button (lower right).
MATLAB code: MATLAB is an expensive high-level language and interactive environment that enables fast coding with matrix based data. I wrote these MATLAB programs for several specific purposes.

convert_map_volume.m (MATLAB)
Reads 32-bit floating point volume data from a CCP4 format file, and writes the 3-D array into an 8-bit unsigned Space Volume Type 1 or a 32-bit single file format NIFTI(.nii) file. Both formats are readable by Space Software, and FSL and other programs read the NIFTI format.
[To Do: Add argument to convert_map_volume.m for desired output volume format. Document and include the NIFTI MATLAB tools I used: make_nii.m and save_nii.m
"Part of this file is copied and modified under GNU license from MRI_TOOLBOX developed by CNSP in Flinders University, Australia NIFTI data format can be found on: http://nifti.nimh.nih.gov
- Jimmy Shen (pls@rotman-baycrest.on.ca) "
2008-09-14 Rewrote convert_map_volume.m to include NIFTI 32-bit output option.
[To Do: Document and attach this new code.]

rotation_average_image.m (MATLAB)
Average n-fold rotations of an image about a near-central point. The location of the center is refined using a hill climbing algorithm that maximizes the standard deviation of values in the average image.
[To Do: Example usage.]
rotation_average_volume.m (MATLAB)
rotation_nfold_image.m (MATLAB)
remove_background_peak_image.m (MATLAB)
See Description of remove_background_peak_image.m for example image, usage, and output.
[ vPeak fFWHM_range imZeroBackground ] = remove_background_peak_image( strInputFilename, fFWHM, bShow )
For an image with a relatively uniform noisy background (composing a luminance histogram peak), find the peak luminance histogram hump and the range of values above a given fraction of the full-width at half-maximum FWHM. For example, with fFWHM = 1.5 the width of the range of values that is set to zero is 1.5 times the derived FWHM.
2008-09-24 I'm not going to use this for this project, although this is why wrote it. I think it will be easier to mask the backgound on the whole volume using the Space Software 3-D Fill tool. This avoids dealing with interior pixels/voxels that this program removes from the interior.

The homogeneous (coordinate) transformation matrix for 3D points

FSL, and many other graphics programs use a 4x4 spatial transfomation matrix of the following form. For the 3D case, this 4x4 matrix operator performs the rotation given by $R(\alpha,\beta,\gamma)$ , followed by a translation given by :

$\displaystyle T = {\small \begin{pmatrix}\cos\alpha \cos\beta & \cos\alpha \sin... ...ta \sin\gamma & \cos\beta \cos\gamma & z_t 0 & 0 & 0 & 1 \end{pmatrix} } .$

The order of operations is critical. The matrix in represents the following sequence of transformations:
1. Roll by $\gamma$
2. Pitch by $\beta$
3. Yaw by $\alpha$
4. Translate by