Usage #1:  3dBallMatch dataset [radius]

Usage #2:  3dBallMatch [options]

where the pitifully few options are:

 -input dataset  =  read this dataset

 -ball  radius   =  set the radius of the 3D ball to match (mm)

 -spheroid a b   =  match with a spheroid of revolution, with principal
                    axis radius of 'a' and secondary axes radii 'b'
                    ++ this option is considerably slower

* This program tries to find a good match between a ball (filled sphere)
  of the given radius (in mm) and a dataset. The goal is to find a crude
  approximate center of the brain quickly.

* The output can be used to re-center a dataset so that its coordinate
  origin is inside the brain and/or as a starting point for more refined
  3D alignment. Sample scripts are given below.

* The reason for this program is that not all brain images are even
  crudely centered by using the center-of-mass ('3dAllineate -cmass')
  as a starting point -- if the volume covered by the image includes
  a lot of neck or even shoulders, then the center-of-mass may be
  far from the brain.

* If you don't give a radius, the default is 72 mm, which is about the
  radius of an adult human brain/cranium. A larger value would be needed
  for elephant brain images. A smaller value for marmosets.

* For advanced use, you could try a prolate spheroid, using something like
     3dBallMatch -input Fred.nii -spheroid 90 70
  for a human head image (that was not skull stripped). This option is
  several times slower than the 'ball' option, as multiple spheroids have
  to be correlated with the input dataset.

* This program does NOT work well with datasets containing large amounts
  of negative values or background junk -- such as I've seen with animal
  MRI scans and CT scans. Such datasets will likely require some repair
  first, such as cropping (cf. 3dZeropad), to make this program useful.

* Frankly, this program may not be that useful for any purpose :(

* The output is text to stdout containing 3 triples of numbers, all on
  one line:
    i j k xs ys zs xd yd zd
    i j k    = index triple of the central voxel
    xs ys zs = values to use in '3drefit -dxxorigin' (etc.)
               to make (i,j,k) be at coordinates (x,y,z)=(0,0,0)
    xd yd zd = DICOM-order (x,y,z) coordinates of (i,j,k) in the
               input dataset

* The intention is that this output line be captured and then the
  appropriate pieces be used for some higher purpose.

Below is a script to process all the entries in a directory.


# optional: start a virtual X11 server
  set xdisplay = `count -dig 1 3 999 R1`
  echo " -- trying to start Xvfb :${xdisplay}"
  Xvfb :${xdisplay} -screen 0 1024x768x24 >& /dev/null &
  sleep 1
  set display_old = $DISPLAY
  setenv DISPLAY :${xdisplay}

# loop over all subjects
  foreach sss ( sub-?????_T1w.nii.gz )

# extract subject ID code
    set sub = `echo $sss | sed -e 's/sub-//' -e 's/_T1w.nii.gz//'`

# skip if already finished
    if ( -f $sub.match   ) continue
    if ( -f $sub.sag.jpg ) continue
    if ( -f $sub.cor.jpg ) continue

# run the program, save output to a file
    3dBallMatch $sss > $sub.match

# capture the output for use below
    set ijk = ( `cat $sub.match` )
    echo $sub $ijk

# run afni to make some QC images
    afni -DAFNI_NOSPLASH=YES                            \
         -DAFNI_NOPLUGINS=YES                           \
         -com "OPEN_WINDOW A.sagittalimage"             \
         -com "OPEN_WINDOW A.coronalimage"              \
         -com "SET_IJK $ijk[1-3]"                       \
         -com "SAVE_JPEG A.sagittalimage $sub.sag.jpg"  \
         -com "SAVE_JPEG A.coronalimage $sub.cor.jpg"   \
         -com "QUITT"                                   \

# end of loop over subject

# kill the virtual X11 server (if it was started above)
  sleep 1
  killall Xvfb

# make a movie of the sagittal slices
  im_to_mov -resize -prefix Bsag -npure 4 -nfade 0 *.sag.jpg
# make a movie of the coronal slices
  im_to_mov -resize -prefix Bcor -npure 4 -nfade 0 *.cor.jpg
exit 0

This script is an extension of the one above, where it uses
3dAllineate to align the human brain image to the MNI template,
guided by the initial point computed by 3dBallMatch. The output
of 3dAllineate is the coordinate of the center of the original
volume, in the first 3 values stored in '*Aparam.1D' file.
 * Note that the 3dAllineate step presumes that the input
   dataset is a T1-weighted volume. A different set of options would
   have to be used for an EPI (T2*-weighted) or T2-weighted volume.
 * This script worked pretty well for putting the crosshairs at
   the 'origin' of the brain -- near the anterior commissure.
   Of course, you will need to evaluate its performance yourself.


# optional: start Xvfb to avoid the AFNI GUI starting visibly
  set xdisplay = `count -dig 1 3 999 R1`
  echo " -- trying to start Xvfb :${xdisplay}"
  Xvfb :${xdisplay} -screen 0 1024x768x24 >& /dev/null &
  sleep 1
  set display_old = $DISPLAY
  setenv DISPLAY :${xdisplay}

# loop over datasets in the current directory
  foreach sss ( anat_sub?????.nii.gz )

# extract the subject identfier code (the '?????')
    set sub = `echo $sss | sed -e 's/anat_sub//' -e 's/.nii.gz//'`

# if 3dAllineate was already run on this, skip to next dataset
    if ( -f $sub.Aparam.1D ) continue

# find the 'center' voxel location with 3dBallMatch
    if ( ! -f $sub.match ) then
      echo "Running 3dBallMatch $sss"
      3dBallMatch $sss | tee $sub.match

# extract results from 3dBallMatch output
# in this case, we want the final triplet of coordinates
    set ijk = ( `cat $sub.match` )
# set shift range to be 55 mm about 3dBallMatch coordinates
    set  xd = $ijk[7] ; set xbot = `ccalc "${xd}-55"` ; set xtop = `ccalc "${xd}+55"`
    set  yd = $ijk[8] ; set ybot = `ccalc "${yd}-55"` ; set ytop = `ccalc "${yd}+55"`
    set  zd = $ijk[9] ; set zbot = `ccalc "${zd}-55"` ; set ztop = `ccalc "${zd}+55"`

# Align the brain image volume with 3dAllineate:
#  match to 'skull on' part of MNI template = sub-brick [1]
#  only save the parameters, not the final aligned dataset
    3dAllineate                                          \
      -base ~/abin/MNI152_2009_template_SSW.nii.gz'[1]'  \
      -source $sss                                       \
      -parang 1 $xbot $xtop                              \
      -parang 2 $ybot $ytop                              \
      -parang 3 $zbot $ztop                              \
      -prefix NULL -lpa                                  \
      -1Dparam_save $sub.Aparam.1D                       \
      -conv 3.666 -fineblur 3 -num_rtb 0 -norefinal -verb

# 1dcat (instead of cat) to strip off the comments at the top of the file
# the first 3 values in 'param' are the (x,y,z) shifts
# Those values could be used in 3drefit to re-center the dataset
    set param = ( `1dcat $sub.Aparam.1D` )

# run AFNI to produce the snapshots with crosshairs at
# the 3dBallMatch center and the 3dAllineate center
# - B.*.jpg = 3dBallMatch result in crosshairs
# - A.*.jpg = 3dAllineate result in crosshairs
    afni -DAFNI_NOSPLASH=YES                             \
         -DAFNI_NOPLUGINS=YES                            \
         -com "OPEN_WINDOW A.sagittalimage"              \
         -com "SET_IJK $ijk[1-3]"                        \
         -com "SAVE_JPEG A.sagittalimage B.$sub.sag.jpg"  \
         -com "SET_DICOM_XYZ $param[1-3]"                \
         -com "SAVE_JPEG A.sagittalimage A.$sub.sag.jpg" \
         -com "QUITT"                                    \

# End of loop over datasets

# stop Xvfb (only needed if it was started above)
  sleep 1
  killall Xvfb

# make movies from the resulting images
  im_to_mov -resize -prefix Bsag -npure 4 -nfade 0 B.[1-9]*.sag.jpg
  im_to_mov -resize -prefix Asag -npure 4 -nfade 0 A.[1-9]*.sag.jpg
exit 0

HOW IT WORKS (approximately)
1] Create the automask of the input dataset (as in 3dAutomask).
   + This is a 0/1 binary marking of outside/inside voxels.
   + Then convert it to a -1/+1 mask instead.

2] Create a -1/+1 mask for the ball [-1=outside, +1=inside],
   inside a rectangular box.

3] Convolve these 2 masks (using FFTs for speed).
   + Basically, this is moving the ball around, then adding up
     the voxel counts where the masks match sign (both positive
     means ball and dataset are both 'inside'; both negative
     means ball and dataset are both 'outside'), and subtracting
     off the voxel counts where the mask differ in sign
     (one is 'inside' and one is 'outside' == not matched).
   + That is, the convolution value is the sum of matched voxels
     minus the sum of mismatched voxels, at every location of
     offset (i,j,k) of the corner of the ball mask.
   + The ball mask is in a cube of side 2*radius, which has volume
     8*radius^3. The volume of the ball is 4*pi/3*radius^3, so the
     inside of the ball is about 4*pi/(3*8) = 52% of the volume of the cube
     -- that is, inside and outside voxels are (roughly) matched, so they
     have (approximately) equal weight.
   + Most of the CPU time is in the 3D FFTs required.

4] Find the centroid of the locations where the convolution
   is positive (matches win over non-matches) and at least 5%
   of the maximum convolution. This centroid gives (i,j,k).

Why the centroid? I found that the peak convolution location
is not very stable, as a lot of locations have results barely less
than the peak value -- it was more stable to average them together.

WHY 'ball' NOT 'sphere'?
 * Because a 'sphere' is a 2D object, the surface of the 3D object 'ball'.
 * Because my training was in mathematics, where precise terminology has
   been developed and honed for centuries.
 * Because I'm yanking your chain. Any other questions? No? Good.

By RWCox, September 2020 (the year it all fell apart).
Delenda est. Never forget.