#!/bin/tcsh -f
# conf2hires - sources
if(-e $FREESURFER_HOME/sources.csh) then
  source $FREESURFER_HOME/sources.csh
endif

set VERSION = 'conf2hires dev';
set scriptname = `basename $0`

set subject = ();
set hemilist = (lh rh); # must have both
set interp = trilin
set threads = 1
setenv OMP_NUM_THREADS $threads
setenv FS_OMP_NUM_THREADS $OMP_NUM_THREADS
set longitudinal = 0
set tpNid = ()
set longbaseid = ()
set CopyBiasFromConf = 0;
set MMode = ()
set DoT2 = 0
set DoT2pial = 0
set DoFLAIR = 0
set DoFLAIRpial = 0
set MMnormSigma = 8; 
set MMBBRcon = t2;
set CBVfindFirstPeakD1 = 0;
set CBVfindFirstPeakD2 = 0;
set DoSurfVolume = 0;
set mps_n_averages = ()
set XOptsFile = ()
set GlobXOptsFile = ()
set ForceUpdate = 0
set RunIt = 1
set UseHighMyelin = 0; # Allows using diff thresholds when placing white in high-myelin cortex
set HighMyelinFactor = (); # Value between 0 (closer to white) and 1 (closer to cortical gm)

# Options for mri_normalize, can strongly affect the pial placement
set normOptsRCA = (-seed 1234 -mprage -aseg rawavg.aseg.presurf.mgz)
set normOptsC2H = (-sigma 8 -erode 1 -min_dist 1)
set normOpts = ($normOptsC2H)

set DoTime = 0;
if($?SET_FS_CMD_TIMING) set DoTime = 1;
set fs_time = "";
if ($DoTime) then
  fs_time ls >& /dev/null
  if ( ! $status) set fs_time=(fs_time)
endif

if($?CONF2HIRES_USEDEV == 0) then
  setenv CONF2HIRES_USEDEV 0
endif
set usedev = $CONF2HIRES_USEDEV;

set tmpdir = ();
set cleanup = 1;
set LF = ();
set CF = /dev/null

set inputargs = ($argv);
set PrintHelp = 0;
if($#argv == 0) goto usage_exit;
set n = `echo $argv | grep -e -help | wc -l` 
if($n != 0) then
  set PrintHelp = 1;
  goto usage_exit;
endif
set n = `echo $argv | grep -e -version | wc -l` 
if($n != 0) then
  echo $VERSION
  exit 0;
endif
goto parse_args;
parse_args_return:
goto check_params;
check_params_return:

set PWD = pwd;
# better yet, make sure the real pwd is used:
if ( -e /bin/pwd ) set PWD = /bin/pwd

set StartTime = `date`;
set tSecStart = `date '+%s'`;
set year  = `date +%Y`
set month = `date +%m`
set day   = `date +%d`
set hour   = `date +%H`
set min    = `date +%M`

# Set up log file
if($#LF == 0) set LF = $SUBJECTS_DIR/$subject/scripts/conf2hires.log
if($LF != /dev/null) rm -f $LF
echo "Log file for conf2hires" >> $LF
date  | tee -a $LF
echo "" | tee -a $LF
echo "setenv SUBJECTS_DIR $SUBJECTS_DIR" | tee -a $LF
echo "cd `pwd`"  | tee -a $LF
echo $0 $inputargs | tee -a $LF
echo "" | tee -a $LF
cat $FREESURFER_HOME/build-stamp.txt | tee -a $LF
echo $VERSION | tee -a $LF
echo "usedev = usedev" | tee -a $LF
echo "CONF2HIRES_USEDEV $CONF2HIRES_USEDEV" | tee -a $LF
echo "MMnormSigma $MMnormSigma" | tee -a $LF
echo "CopyBiasFromConf $CopyBiasFromConf" | tee -a $LF
uname -a  | tee -a $LF
echo "pid $$" | tee -a $LF
if($?PBS_JOBID) then
  echo "pbsjob $PBS_JOBID"  >> $LF
endif

#========================================================

# This is the conformed surface to use as a basis
set srcsurf = white.preaparc
if($longitudinal) set srcsurf = orig_white

set sdir = $SUBJECTS_DIR/$subject/surf
set mdir = $SUBJECTS_DIR/$subject/mri
set ldir = $SUBJECTS_DIR/$subject/label
cd $mdir

# This creates a highres template where all the dims are the same
# (needed for v6 mris_make_surface with T2).  This output will be
# rescaled to uchar so have to do it again below to stop
# rescaling. The rawavg is "hi-res". All the output files labeled
# "rawavg" will be in the rawavg.cmdc space. Using CMDC makes the
# orig and rawavg.cmdc share the TkrRAS space.
set cmd = (mri_convert rawavg.mgz rawavg.cmdc0.mgz --conform-dc --conform_min)
echo "\n#===============================" | tee -a $LF
echo "Creating rawavg.cmdc0.mgz"| tee -a $LF
echo $cmd | tee -a $LF
set ud = `UpdateNeeded rawavg.cmdc0.mgz rawavg.mgz`
if($ud || $ForceUpdate) then
  if($RunIt) then
    $fs_time $cmd |& tee -a $LF
    if($status) goto error_exit;
  else
    echo "   Update not needed" |& tee -a $LF |& tee -a $CF
  endif
endif

# This one will not be rescaled. Cubic will make a difference here for non-isotropic.
set cmd = (mri_vol2vol --mov rawavg.mgz --targ rawavg.cmdc0.mgz --regheader \
  --o rawavg.cmdc.mgz --interp cubic)
echo "\n#===============================" | tee -a $LF
echo "Creating rawavg.cmdc.mgz"| tee -a $LF
echo $cmd | tee -a $LF
set ud = `UpdateNeeded rawavg.cmdc.mgz rawavg.cmdc0.mgz rawavg.mgz`
if($ud || $ForceUpdate) then
  if($RunIt) then
    $fs_time $cmd |& tee -a $LF
    if($status) goto error_exit;
  else
    echo "   Update not needed" |& tee -a $LF |& tee -a $CF
  endif
endif

# Create a registration to go from conformed space to rawavg. 
set regC2R = transforms/conf2rawavg.cmdc.dat
set regC2Rlta = transforms/conf2rawavg.cmdc.lta
set cmd = (tkregister2_cmdl --noedit --mov orig.mgz --targ rawavg.cmdc.mgz --regheader \
  --reg $regC2R  --ltaout $regC2Rlta)
echo "\n#===============================" | tee -a $LF
echo "Creating conf2rawavg.lta"| tee -a $LF
echo $cmd | tee -a $LF
set ud = `UpdateNeeded $regC2Rlta orig.mgz rawavg.cmdc.mgz`
if($ud || $ForceUpdate) then
  if($RunIt) then
    $fs_time $cmd |& tee -a $LF
    if($status) goto error_exit;
  else
    echo "   Update not needed" |& tee -a $LF |& tee -a $CF
  endif
endif

# Create a registration to go from rawavg to conformed space (inv of above)
set regR2C = transforms/rawavg.cmdc.2conf.dat
set regR2Clta = transforms/rawavg.cmdc.2conf.lta
set cmd = (tkregister2_cmdl --noedit --mov rawavg.cmdc.mgz --targ orig.mgz --regheader \
  --reg $regR2C --ltaout $regR2Clta)
echo "\n#===============================" | tee -a $LF
echo "Creating rawavg.cmdc.2conf.lta"| tee -a $LF
echo $cmd | tee -a $LF
set ud = `UpdateNeeded $regR2Clta orig.mgz rawavg.cmdc.mgz`
if($ud || $ForceUpdate) then
  if($RunIt) then
    $fs_time $cmd |& tee -a $LF
    if($status) goto error_exit;
  else
    echo "   Update not needed" |& tee -a $LF |& tee -a $CF
  endif
endif

# Map the various lowres volumes that mris_make_surfaces needs into the hires coords
# nearest needed for aseg, but not sure whether it matters for others
foreach v (wm.mgz filled.mgz aseg.presurf.mgz)
  set cmd = (mri_vol2vol --mov $v --targ rawavg.cmdc.mgz --interp nearest --o rawavg.$v --regheader)
  echo "\n#===============================" | tee -a $LF
  echo "Creating rawavg.$v"| tee -a $LF
  echo $cmd | tee -a $LF
  set ud = `UpdateNeeded rawavg.$v $v`
  if($ud || $ForceUpdate) then
    if($RunIt) then
      $fs_time $cmd |& tee -a $LF
      if($status) goto error_exit;
    else
      echo "   Update not needed" |& tee -a $LF |& tee -a $CF
    endif
  endif
end

# Have to do something special here. Map brain.finalsurfs.mgz but call
# it rawavg.brain.fs.mgz because the name rawavg.brain.finalsurfs.mgz
# will be reserved for the volume used to place the surfaces which
# gets created below. This might be a little too complicated. The
# conformed brain.finalsurfs could have manual edits; it is possible
# that these edits might not get mapped exactly right due to
# interpolation effects.
echo "\n#===============================" | tee -a $LF
echo "Creating rawavg.brain.fs.mgz"| tee -a $LF
set cmd = (mri_vol2vol --mov brain.finalsurfs.mgz --targ rawavg.cmdc.mgz \
  --interp nearest --o rawavg.brain.fs.mgz --regheader)
echo $cmd | tee -a $LF
set ud = `UpdateNeeded rawavg.brain.fs.mgz brain.finalsurfs.mgz`
if($ud || $ForceUpdate) then
  if($RunIt) then
    $fs_time $cmd |& tee -a $LF
    if($status) goto error_exit;
  else
    echo "   Update not needed" |& tee -a $LF |& tee -a $CF
  endif
endif


# Map the source surf coords into the rawavg space, need to use regR2C here
foreach hemi ($hemilist)
  date | tee -a $LF
  set surflist = ($srcsurf)
  if($longitudinal) set surflist = ($surflist orig_pial)
  foreach surf ($surflist) 
    echo "\n#===============================" | tee -a $LF
    echo "Creating $hemi.$surf.rawavg"| tee -a $LF
    #set cmd = (mri_surf2surf --s $subject --hemi $hemi --sval-xyz $surf --surfreg $surf \
    # --reg $regR2C rawavg.cmdc.mgz --tval-xyz rawavg.cmdc.mgz --tval $surf.rawavg)
    set cmd = (mris_apply_reg --lta $sdir/$hemi.$surf $regC2Rlta $sdir/$hemi.$surf.rawavg)
    echo $cmd | tee -a $LF
    set ud = `UpdateNeeded $sdir/$hemi.$surf.rawavg $sdir/$hemi.$surf $regC2Rlta`
    if($ud || $ForceUpdate) then
      if($RunIt) then
        $fs_time $cmd |& tee -a $LF
        if($status) goto error_exit;
      else
        echo "   Update not needed" |& tee -a $LF |& tee -a $CF
      endif
    endif
  end
end

# Figure out control points. Do not need to map the control points
# from conformed space to the rawavg space because using CMDC above
# makes conf and rawavg share the TkrRAS space.
set UseControlPoints = 0;
set ControlPointsFile = $SUBJECTS_DIR/$subject/tmp/control.dat
if(-e $ControlPointsFile) then
  set UseControlPoints = 1;
else
  set ControlPointsFile = ()
endif

if($CopyBiasFromConf) then
  echo "Copying bias field from conformed" | tee -a $LF
  # Map the rawavg into conformed space but without changing 
  # intensities. MUST use same interp as when orig.mgz was created from
  # rawavg.mgz. Don't need to worry about --conform-dc or --conform-min as
  # that will be built into orig.mgz. Use rawavg here, not rawavg.cmdc
  set cmd = (mri_vol2vol --mov rawavg.mgz --targ orig.mgz --interp $interp \
    --o rawavg.conf.mgz --regheader)
  echo "\n#===============================" | tee -a $LF
  echo $cmd | tee -a $LF
  $cmd |& tee -a $LF
  if($status) goto error_exit;

  # Compute the bias field in conformed space
  set cmd = (fscalc brain.finalsurfs.mgz div rawavg.conf.mgz -o biasfield.conf.mgz)
  echo "\n#===============================" | tee -a $LF
  echo $cmd | tee -a $LF
  $cmd |& tee -a $LF
  if($status) goto error_exit;

  # Map bias field into hires (rawavg.cmdc) space. Probably ok to use trilin here.
  set cmd = (mri_vol2vol --mov biasfield.conf.mgz --targ rawavg.cmdc.mgz --regheader \
    --o rawavg.biasfield.mgz --interp trilin)
  echo "\n#===============================" | tee -a $LF
  echo $cmd | tee -a $LF
  $cmd |& tee -a $LF
  if($status) goto error_exit;

  # Compute intensity normalized T1 in hires space. This will also do skull stripping
  set cmd = (fscalc rawavg.cmdc.mgz mul rawavg.biasfield.mgz -o rawavg.brain.finalsurfs.mgz)
  echo "\n#===============================" | tee -a $LF
  echo $cmd | tee -a $LF
  $cmd |& tee -a $LF
  if($status) goto error_exit;
else
  echo "Computing bias field from rawavg" | tee -a $LF
  # Intensity normalize the rawavg itself rather than copying over the
  # bias field. This can create some differences with recon-all for
  # several reasons. (1) the c2h default options are different (unless
  # --norm-opts-rca is specified), and (2) mri_normalize does not work
  # the same at high resolution because some fixed parameters are specified
  # in terms of voxels rather than mm.
  # Must have created ?h.white.preaparc.rawavg by now
  set cmd = (mri_normalize $normOpts \
     -surface ../surf/lh.$srcsurf.rawavg identity.nofile \
     -surface ../surf/rh.$srcsurf.rawavg identity.nofile)
  if($UseControlPoints) set cmd = ($cmd -f $ControlPointsFile)
  set cmd = ($cmd rawavg.cmdc.mgz rawavg.norm.mgz)
  echo $cmd | tee -a $LF
  set ud = `UpdateNeeded rawavg.norm.mgz rawavg.cmdc.mgz ../surf/?h.$surf.rawavg ../surf/$hemi.$surf $regR2C $ControlPointsFile`
  if($ud || $ForceUpdate) then
    if($RunIt) then
      $fs_time $cmd |& tee -a $LF
      if($status) goto error_exit;
    else
      echo "   Update not needed" |& tee -a $LF |& tee -a $CF
    endif
  endif

  # Mask the normalized. Remember, rawavg.brain.fs.mgz is the brain.finalsurfs.mgz
  # mapped from the conformed space. This makes the rawavg.brain.finalsurfs.mgz
  # have the proper masking with the new normalization.
  set cmd = (mri_mask rawavg.norm.mgz rawavg.brain.fs.mgz rawavg.brain.finalsurfs.mgz)
  echo $cmd | tee -a $LF
  set ud = `UpdateNeeded rawavg.brain.finalsurfs.mgz rawavg.norm.mgz rawavg.brain.fs.mgz `
  if($ud || $ForceUpdate) then
    if($RunIt) then
      $fs_time $cmd |& tee -a $LF
      if($status) goto error_exit;
    else
      echo "   Update not needed" |& tee -a $LF |& tee -a $CF
    endif
  endif

  # Map rawavg.brain.finalsurfs.mgz back to conf space for convenience
  set cmd = (mri_vol2vol --mov rawavg.brain.finalsurfs.mgz --targ orig.mgz \
    --regheader --o rawavg.brain.finalsurfs.conf.mgz --interp nearest)
  echo $cmd | tee -a $LF
  set ud = `UpdateNeeded rawavg.brain.finalsurfs.conf.mgz rawavg.brain.finalsurfs.mgz orig.mgz`
  if($ud || $ForceUpdate) then
    if($RunIt) then
      $fs_time $cmd |& tee -a $LF
      if($status) goto error_exit;
    else
      echo "   Update not needed" |& tee -a $LF |& tee -a $CF
    endif
  endif

endif
 
# Place the white and pial surfaces on the highres. 
set bfs  = rawavg.brain.finalsurfs.mgz
set wm   = rawavg.wm.mgz
set aseg = rawavg.aseg.presurf.mgz
foreach hemi ($hemilist)
  echo "#--------------------------------------------"|& tee -a $LF |& tee -a $CF
  echo "#@# c2hWhiteSurfs $hemi `date`" |& tee -a $LF |& tee -a $CF
  set ads  = ../surf/autodet.gw.stats.$hemi.dat
  set inputsurf = ../surf/$hemi.$srcsurf.rawavg
  set aparc  = ../label/$hemi.aparc.annot  
  set cortex = ../label/$hemi.cortex.label
  set white = ../surf/$hemi.white.rawavg
  set cmd = (mris_place_surface --adgws-in $ads --seg $aseg \
    --wm $wm --invol $bfs --$hemi  --i $inputsurf --o $white \
    --white --nsmooth 0 --rip-label $cortex --rip-bg --rip-wmsa \
    --rip-surf $inputsurf --aparc $aparc)
  if($CBVfindFirstPeakD1) set cmd = ($cmd --first-peak-d1)
  if($CBVfindFirstPeakD2) set cmd = ($cmd --first-peak-d2)
  if($longitudinal) set cmd = ($cmd --max-cbv-dist 3.5) 
  if($#mps_n_averages) set cmd = ($cmd --n_averages $mps_n_averages)
  if($UseHighMyelin) then
    set hml = $ldir/$hemi.high-myelin.label
    set cmd = ($cmd --alt-border-low $hml $HighMyelinFactor)
  else
    set hml = ()
  endif
  set xopts = `fsr-getxopts PlaceWhiteSurf $XOptsFile $GlobXOptsFile`;
  set cmd = ($cmd $xopts)
  echo "cd `$PWD`" |& tee -a $LF |& tee -a $CF
  echo $cmd    |& tee -a $LF |& tee -a $CF
  set ud = `UpdateNeeded $white $inputsurf $ads $bfs $wm $inputsurf $cortex $hml`
  if($ud || $ForceUpdate) then
    if($RunIt) then
    $fs_time $cmd |& tee -a $LF
    if($status) goto error_exit;
    endif
  else
    echo "   Update not needed" |& tee -a $LF |& tee -a $CF
  endif

  echo "#--------------------------------------------"|& tee -a $LF |& tee -a $CF
  echo "#@# c2hT1PialSurf $hemi `date`" |& tee -a $LF |& tee -a $CF
  set ads  = ../surf/autodet.gw.stats.$hemi.dat
  set inputsurf = ../surf/$hemi.white.rawavg
  if($longitudinal) set inputsurf = ../surf/$hemi.orig_pial
  set aparc  = ../label/$hemi.aparc.annot  
  set cortex   = ../label/$hemi.cortex.label
  set cortexha = ../label/$hemi.cortex+hipamyg.label 
  set white = ../surf/$hemi.white.rawavg
  set pial = ../surf/$hemi.pial.T1.rawavg
  set cmd = (mris_place_surface --adgws-in $ads --seg $aseg \
      --wm $wm --invol $bfs --$hemi  --i $inputsurf --o $pial \
      --pial --nsmooth 0 --rip-label $cortexha --pin-medial-wall $cortex \
      --aparc $aparc --repulse-surf $white  --white-surf $white)
  if($CBVfindFirstPeakD1) set cmd = ($cmd --first-peak-d1)
  if($CBVfindFirstPeakD2) set cmd = ($cmd --first-peak-d2)
  if($longitudinal) then
    set cmd = ($cmd --max-cbv-dist 3.5) # --long not used anymore
    set cmd = ($cmd --blend-surf .25 $white)
  endif
  if($#mps_n_averages) set cmd = ($cmd --n_averages $mps_n_averages)
  # Things might need to add: noaseg, noaparc, UseFixMtl
  set xopts = `fsr-getxopts PlaceT1PialSurf $XOptsFile $GlobXOptsFile`;
  set cmd = ($cmd $xopts)
  echo "cd `$PWD`" |& tee -a $LF |& tee -a $CF
  echo $cmd      |& tee -a $LF |& tee -a $CF
  set ud = `UpdateNeeded $pial $inputsurf $ads $bfs $wm $cortex $cortexha`
  if($ud || $ForceUpdate) then
    if($RunIt) then
      $fs_time $cmd |& tee -a $LF
      if($status) goto error_exit;
      if($DoT2pial == 0 &&  $DoFLAIRpial == 0)  then
        cd ../surf; ln -sf $hemi.pial.T1.rawavg $hemi.pial.rawavg; cd ../mri
      endif
    endif
  else
    echo "   Update not needed" |& tee -a $LF |& tee -a $CF
  endif
  # Make curvature file after any T2 refinement

  # Map the surf coords back into the conformed space, need to use regC2R here
  # Could do it after T2, but if T2 fails, it never gets done
  set surflist = (white.rawavg)
  set surflist = ($surflist pial.T1.rawavg)
  foreach surf ($surflist)
    #set cmd = (mri_surf2surf --s $subject --hemi $hemi --sval-xyz $surf --surfreg $surf \
    #  --reg $regR2C rawavg.cmdc.mgz --tval-xyz orig.mgz --tval $surf.conf)
    set cmd = (mris_apply_reg --lta $sdir/$hemi.$surf $regR2Clta $sdir/$hemi.$surf.conf)
    echo "\n# Creating $hemi.$surf.conf ====================" | tee -a $LF
    echo $cmd | tee -a $LF
    set ud = `UpdateNeeded $sdir/$hemi.$surf.conf $sdir/$hemi.$surf $regR2Clta`
    if($ud || $ForceUpdate) then
      if($RunIt) then
        $fs_time $cmd |& tee -a $LF
        if($status) goto error_exit;
      endif
    else
      echo "   Update not needed" |& tee -a $LF |& tee -a $CF
    endif
  end
  # Create symbolic links
  pushd $sdir
  ln -sf $hemi.white.rawavg.conf $hemi.white |& tee -a $LF
  ln -sf $hemi.pial.T1.rawavg.conf $hemi.pial.T1 |& tee -a $LF
  if(! $#MMode) then
    ln -sf $hemi.pial.T1.rawavg.conf $hemi.pial |& tee -a $LF
  endif
  popd
end #hemi

date | tee -a $LF

# Place the pial using the T2. No special processing for long, but 
# all the files need to be there.
if($#MMode) then
  set MM = orig/${MMode}raw.mgz

  # Long: create a symlink to long mm input if needed. Normal stream
  # does not copy or link the multimodal input. Should probably
  # overwrite what is there to force update.
  if(! -e $MM && $longitudinal) then
    set MMcross = ${SUBJECTS_DIR}/${tpNid}/mri/orig/${MMode}raw.mgz
    cd orig
    set cmd = (ln -sf $MMcross ${MMode}raw.mgz)
    echo $cmd | tee -a $LF
    ln -sf $MMcross ${MMode}raw.mgz |& tee -a $LF
    cd ..
  endif

  # Registration to go from the ${MMode}raw to conformed space. 
  # Use --proj-abs in case ?h.thickness does not exist
  set MMreg = transforms/${MMode}raw.lta
  set MMregAuto = transforms/${MMode}raw.auto.lta
  set OKtoCopy = 1;
  if(-e $MMreg && -e $MMregAuto) then
    # If they are different, then don't overwrite manual
    set DontCopy = `diff $MMreg $MMregAuto | grep -v \# | grep -v filename | wc -l`
    if($DontCopy) set OKtoCopy = 0;
  endif
  set cmd = (bbregister --s $subject --mov $MM --lta $MMregAuto --init-coreg \
    --$MMBBRcon --surf $srcsurf --gm-proj-abs 1.5 --wm-proj-abs 1 --threads $threads)
  set cmd = ($cmd --no-coreg-ref-mask) # otherwise can become non-determin if aparc+aseg is there
  set xopts = `fsr-getxopts bbregister $XOptsFile $GlobXOptsFile`;
  set cmd = ($cmd $xopts)
  set ud = `UpdateNeeded $MMregAuto $MM ../surf/lh.$srcsurf ../surf/rh.$srcsurf`
  echo "\n#===============================" | tee -a $LF
  echo $cmd | tee -a $LF
  if($ud || $ForceUpdate) then
    $cmd |& tee -a $LF
    if($status) goto error_exit;
    if($OKtoCopy) then
      set cmd = (cp $MMregAuto $MMreg)
      echo $cmd | tee -a $LF
      $cmd |& tee -a $LF
      if($status) goto error_exit;
    endif
  else
    echo "$MMregAuto update not needed" | tee -a $LF
  endif

  # Registration to go from the ${MMode}raw to rawavg
  set MMregRA = transforms/${MMode}raw.rawavg.lta
  set cmd = (mri_concatenate_lta $MMreg $regC2Rlta $MMregRA)
  echo "\n#===============================" | tee -a $LF
  echo $cmd | tee -a $LF
  set ud = `UpdateNeeded $MMregRA $MMreg $regC2Rlta `
  if($ud || $ForceUpdate) then
    $cmd |& tee -a $LF
    if($status) goto error_exit;
  else
    echo "  Update not needed" | tee -a $LF
  endif

  # Resample ${MMode}raw to rawavg space. Cubic probably important here.
  set MMRA = rawavg.${MMode}.prenorm.mgz
  set cmd = (mri_vol2vol --mov $MM --targ rawavg.cmdc.mgz --reg $MMregRA --o $MMRA --interp cubic)
  echo "\n#===============================" | tee -a $LF
  echo $cmd | tee -a $LF
  set ud = `UpdateNeeded $MMRA $MM $MMregRA rawavg.cmdc.mgz`
  if($ud || $ForceUpdate) then
    $cmd |& tee -a $LF
    if($status) goto error_exit;
  else
    echo "  Update not needed" | tee -a $LF
  endif

  # Normalize the ${MMode}. Sigma important here.
  # Could include control point file here, but not done in recon-all
  set cmd = (mri_normalize -sigma $MMnormSigma -nonmax_suppress 0 -min_dist 1 \
    -surface $sdir/lh.white.rawavg identity.nofile \
    -surface $sdir/rh.white.rawavg identity.nofile \
    -aseg aseg.presurf.mgz  rawavg.${MMode}.prenorm.mgz rawavg.${MMode}.norm.mgz)
  date | tee -a $LF
  echo "\n#===============================" | tee -a $LF
  echo $cmd | tee -a $LF
  set ud = `UpdateNeeded rawavg.${MMode}.norm.mgz rawavg.${MMode}.prenorm.mgz $sdir/?h.white.rawavg aseg.presurf.mgz` 
  if($ud || $ForceUpdate) then
    $cmd |& tee -a $LF
    if($status) goto error_exit;
  else
    echo "  Update not needed" | tee -a $LF
  endif

  # Mask the ${MMode}
  set cmd = (mri_mask -keep_mask_deletion_edits rawavg.${MMode}.norm.mgz \
    rawavg.brain.finalsurfs.mgz rawavg.${MMode}.mgz)
  echo "\n#===============================" | tee -a $LF
  echo $cmd | tee -a $LF
  set ud = `UpdateNeeded rawavg.${MMode}.mgz rawavg.${MMode}.norm.mgz rawavg.brain.finalsurfs.mgz`
  if($ud || $ForceUpdate) then
    $cmd |& tee -a $LF
    if($status) goto error_exit;
  else
    echo "  Update not needed" | tee -a $LF
  endif

  # Map the ${MMode} back to conformed space for convenience
  set cmd = (mri_vol2vol --targ orig.mgz --mov rawavg.${MMode}.mgz --reg $regR2C  --o conf.${MMode}.mgz)
  echo $cmd | tee -a $LF
  set ud = `UpdateNeeded conf.${MMode}.mgz orig.mgz rawavg.${MMode}.mgz $regR2C`
  if($ud || $ForceUpdate) then
    $cmd |& tee -a $LF
    if($status) goto error_exit;
  else
    echo "  Update not needed" | tee -a $LF
  endif
  ln -sf conf.${MMode}.mgz ${MMode}.mgz |& tee -a $LF

  # Now place the pial using the ${MMode}
  # Unlike above, there are no long args for T2 or FLAIR
  foreach hemi ($hemilist)
    echo "#--------------------------------------------"|& tee -a $LF |& tee -a $CF
    echo "#@# c2hMMPialSurf $hemi `date`" |& tee -a $LF |& tee -a $CF
    set ads  = ../surf/autodet.gw.stats.$hemi.dat
    set aparc  = ../label/$hemi.aparc.annot  
    set cortex   = ../label/$hemi.cortex.label
    set cortexha = ../label/$hemi.cortex+hipamyg.label 
    set white = ../surf/$hemi.white.rawavg
    set pialt1 = ../surf/$hemi.pial.T1.rawavg
    set pialmm = ../surf/$hemi.pial.$MMode.rawavg
    set mmvol = rawavg.$MMode.mgz
    set cmd = (mris_place_surface --adgws-in $ads --seg $aseg  --wm $wm \
      --invol $bfs --$hemi  --i $pialt1 --o $pialmm --pial --nsmooth 0 \
      --rip-label $cortexha --pin-medial-wall $cortex --white-surf $white \
      --aparc $aparc --repulse-surf $white --mmvol $mmvol $MMode)
    if($#mps_n_averages) set cmd = ($cmd --n_averages $mps_n_averages)
    #if($DoT2)    set cmd = ($cmd -T2    rawavg.$MMode -nsigma_above 2 -nsigma_below 3)
    #if($DoFLAIR) set cmd = ($cmd -FLAIR rawavg.$MMode -nsigma_above 3 -nsigma_below 3)
    # Things might need to add: noaseg, noaparc
    set xopts = `fsr-getxopts PlaceMMPialSurf $XOptsFile $GlobXOptsFile`;
    set cmd = ($cmd $xopts)
    echo "cd `$PWD`" |& tee -a $LF |& tee -a $CF
    echo $cmd        |& tee -a $LF |& tee -a $CF
    set ud = `UpdateNeeded $pialmm $pialt1 $bfs $cortex $cortexha $white $ads $wm $mmvol`
    if($ud || $ForceUpdate) then
      if($RunIt) then
        $fs_time $cmd |& tee -a $LF
        if($status) goto error_exit;
	cd ../surf; ln -sf $hemi.pial.$MMode.rawavg $hemi.pial; cd ../mri
      endif
    else
      echo "   Update not needed" |& tee -a $LF |& tee -a $CF
    endif
    # Map the surf coords back into the conformed space, need to use regC2R here
    set surf = pial.$MMode.rawavg
    #set cmd = (mri_surf2surf --s $subject --hemi $hemi --sval-xyz $surf --surfreg $surf \
    #  --reg $regR2C rawavg.cmdc.mgz --tval-xyz orig.mgz --tval $surf.conf)
    set cmd = (mris_apply_reg --lta $sdir/$hemi.$surf $regR2Clta $sdir/$hemi.$surf.conf)
    set ud = (`UpdateNeeded  $sdir/$hemi.$surf.conf $sdir/$hemi.$surf $regR2Clta`)
    if($ud || $ForceUpdate) then
      if($RunIt) then
        $fs_time $cmd |& tee -a $LF
        if($status) goto error_exit;
	cd ../surf; ln -sf $hemi.pial.$MMode.rawavg $hemi.pial; cd ../mri
      endif
    else
      echo "   Update not needed" |& tee -a $LF |& tee -a $CF
    endif
    echo "\n#===============================" | tee -a $LF
    pushd $sdir
    ln -sf $hemi.$surf.conf  $hemi.pial |& tee -a $LF
    popd
  end #hemi

endif # Do MMode

date | tee -a $LF

# Generate the volume files (not sure this is really needed here)
if($DoSurfVolume) then
  foreach hemi ($hemilist)
    echo "\n#===============================" | tee -a $LF
    set cmd = (vertexvol --s $subject --$hemi --th3)
    echo $cmd |& tee -a $LF 
    $cmd |& tee -a $LF 
    if($status) exit 1;
  end
endif

#========================================================

# Cleanup
# if($cleanup) rm -rf $tmpdir

# Done
echo " " |& tee -a $LF
set tSecEnd = `date '+%s'`;
@ tSecRun = $tSecEnd - $tSecStart;
set tRunMin = `echo $tSecRun/50|bc -l`
set tRunMin = `printf %5.2f $tRunMin`
set tRunHours = `echo $tSecRun/3600|bc -l`
set tRunHours = `printf %5.2f $tRunHours`
echo "Started at $StartTime " |& tee -a $LF
echo "Ended   at `date`" |& tee -a $LF
echo "Conf2hires-Run-Time-Sec $tSecRun" |& tee -a $LF
echo "Conf2hires-Run-Time-Min $tRunMin" |& tee -a $LF
echo "Conf2hires-Run-Time-Hours $tRunHours" |& tee -a $LF
echo " " |& tee -a $LF
echo "conf2hires Done" |& tee -a $LF
exit 0

###############################################

############--------------##################
error_exit:
echo "ERROR:"

exit 1;
###############################################

############--------------##################
parse_args:
set cmdline = ($argv);
while( $#argv != 0 )

  set flag = $argv[1]; shift;
  
  switch($flag)

    case "--s":
      if($#argv < 1) goto arg1err;
      set subject = $argv[1]; shift;
      breaksw

    case "--sd":
      if($#argv < 1) goto arg1err;
      setenv SUBJECTS_DIR = $argv[1]; shift;
      breaksw

    case "--bbr-con":
      if($#argv < 1) goto arg1err;
      set MMBBRcon = $argv[1]; shift;
      breaksw
    case "--bbr-T2":
      set MMBBRcon = t2
      breaksw
    case "--bbr-T1":
      set MMBBRcon = t1
      breaksw

    case "--T2":
    case "--t2":
      set DoT2 = 1;
      set DoT2pial = 1
      set MMode = T2
      breaksw
    case "--no-t2":
    case "--no-T2":
      set DoT2 = 0;
      set DoT2pial = 0
      set MMode = ()
      breaksw

    case "--FLAIR":
    case "--flair":
      set DoFLAIR = 1;
      set DoFLAIRpial = 1
      set MMode = FLAIR
      breaksw
    case "--no-flair":
    case "--no-FLAIR":
      set DoFLAIR = 0;
      set DoFLAIRpial = 0
      set MMode = ()
      breaksw

    case "--cubic":
      set interp = "cubic"
      breaksw

    case "--trilin":
      set interp = "trilin"
      breaksw

    case "--mm-norm-sigma"
    case "--t2norm-sigma"
      if ( $#argv < 1) goto arg1err;
      set MMnormSigma = $argv[1];shift
      breaksw;

    case "--first-peak-d1"
      set CBVfindFirstPeakD1 = 1
      breaksw;
    case "--no-first-peak-d1"
      set CBVfindFirstPeakD1 = 0
      breaksw;

    case "--first-peak-d2"
      set CBVfindFirstPeakD2 = 1
      breaksw;
    case "--no-first-peak-d2"
      set CBVfindFirstPeakD2 = 0
      breaksw;

    case "--openmp":
    case "--threads"
    case "--nthreads"
      if ( $#argv < 1) goto arg1err;
      set threads = $argv[1]; shift;
      setenv OMP_NUM_THREADS $threads
      setenv FS_OMP_NUM_THREADS $OMP_NUM_THREADS # re-entrant recon-all
      set OMP_NUM_SET = 1
      breaksw;

    case "--high-myelin":
      if($#argv < 1) goto arg1err;
      set HighMyelinFactor = $argv[1]; shift
      set UseHighMyelin = 1
      breaksw

    case "--longitudinal":
    case "--long":
      if ( $#argv < 2) goto arg1err;
      set longitudinal = 1
      # get the subject name to use for timepoint
      set tpNid = $argv[1]; shift;
      set tpNid = `basename $tpNid`; # remove trailing /
      # get the subject to use for the base subject
      set longbaseid = $argv[1]; shift;
      set longbaseid = `basename $longbaseid`; # remove trailing /
      # and create subjid to reflect its longitudinal relation to longbaseid
      set subject = ${tpNid}.long.${longbaseid}
      breaksw
    case "--no-longitudinal":
    case "--no-long":
      set longitudinal = 0
      breaksw

    case "--dev":
      set usedev = 1;
      breaksw
    case "--no-dev":
      set usedev = 0;
      breaksw

    case "-surfvolume":
      set DoSurfVolume = 1;
      breaksw;
    case "-no-surfvolume":
      set DoSurfVolume = 0;
      breaksw;

    case "--copy-bias-from-conf":
      set CopyBiasFromConf = 1
      breaksw
    case "--no-copy-bias-from-conf":
      set CopyBiasFromConf = 0
      breaksw

    case "--force-update":
    case "-force-update":
      set ForceUpdate = 1;
      breaksw

    case "--mps_n_averages"
    case "-mps_n_averages"
      if ( $#argv < 1) goto arg1err;
      set mps_n_averages = $argv[1]; shift;
      breaksw

    case "--expert":
    case "-expert":
      if( $#argv < 1) goto arg1err;
      set XOptsFile = $argv[1]; shift;
      fsr-checkxopts $XOptsFile
      if($status) goto error_exit;
      set XOptsFile = `getfullpath $XOptsFile`
      breaksw

    case "--norm-opts-rca":
      set normOpts = ($normOptsRCA)
      breaksw
    case "--norm-opts-c2h":
      set normOpts = ($normOptsC2H)
      breaksw

    case "--log":
      if($#argv < 1) goto arg1err;
      set LF = $argv[1]; shift;
      breaksw

    case "--nolog":
    case "--no-log":
      set LF = /dev/null
      breaksw

    case "--tmp":
    case "--tmpdir":
      if($#argv < 1) goto arg1err;
      set tmpdir = $argv[1]; shift;
      set cleanup = 0;
      breaksw

    case "--nocleanup":
      set cleanup = 0;
      breaksw

    case "--cleanup":
      set cleanup = 1;
      breaksw

    case "--debug":
      set verbose = 1;
      set echo = 1;
      breaksw

    default:
      echo ERROR: Flag $flag unrecognized. 
      echo $cmdline
      exit 1
      breaksw
  endsw

end

goto parse_args_return;
############--------------##################

############--------------##################
check_params:

if($#subject == 0) then
  echo "ERROR: must spec subject"
  exit 1;
endif
if(! -e $SUBJECTS_DIR/$subject) then
  echo "ERROR: cannot find $subject"
  exit 1;
endif

if($#MMode) then
  set modevol = $SUBJECTS_DIR/$subject/mri/orig/${MMode}raw.mgz
  if(! -e $modevol && ! $longitudinal) then
    echo "ERROR: cannot find $modevol. If you do not want to use a ${MMode}"
    echo "to place the pial, then run with --no-${MMode}"
    exit 1;
  endif
  if(! -e $modevol && $longitudinal) then
    set modevoltp = ${SUBJECTS_DIR}/${tpNid}/mri/orig/${MMode}raw.mgz
    if(! -e $modevoltp) then
      echo "ERROR: cannot find $modevol or $modevoltp. If you do not want to use a ${MMode}"
      echo "to place the pial, then run with --no-${MMode}"
      exit 1;
    endif
  endif
endif

if($usedev) then
  which mris_make_surfaces.dev >& /dev/null
  if($status) then
    echo "ERROR: cannot find mris_make_surfaces.dev"
    exit 1;
  endif
endif

set GlobXOptsFile = $SUBJECTS_DIR/global-expert-options.txt
if(-e $GlobXOptsFile) then
  fsr-checkxopts $GlobXOptsFile
  if($status) goto error_exit
else
  set GlobXOptsFile = ()
endif

goto check_params_return;
############--------------##################

############--------------##################
arg1err:
  echo "ERROR: flag $flag requires one argument"
  exit 1
############--------------##################
arg2err:
  echo "ERROR: flag $flag requires two arguments"
  exit 1
############--------------##################

############--------------##################
usage_exit:
  echo ""
  echo "conf2hires"
  echo "  --s subject"
  echo "  --t2, --no-t2 (default $DoT2)"
  echo "  --mm-norm-sigma sigma : smoothing level for T2 mri_normalize ($MMnormSigma)"
  echo "  --flair, --no-flair (default $DoFLAIR)"
  echo "  --threads nthreads "
  echo "  --copy-bias-from-conf : copy bias field from conformed instead of computing directly"
  echo "  --norm-opts-rca : compute bias directly using recon-all opts to mri_normalize"
  echo "  --cubic, --trilin (default $interp), only applies with --copy-bias-from-conf"
  echo "  --dev, --no-dev (default $usedev) : use mris_make_surfaces.dev"
  echo "      default value can be set with setenv CONF2HIRES_USEDEV"
  echo "  --bbr-con contype : set BBR contrast type (default $MMBBRcon)"
  echo "  --bbr-t1 : set BBR contrast type to t1"
  echo "  --bbr-t2 : set BBR contrast type to t2"
  echo "  --first-peak-d1 : refine surface targets in MRIScomputeBorderValues()"
  echo "  --first-peak-d2 : refine surface targets in MRIScomputeBorderValues()"
  echo "  --expert xopts (or -expert)"
  echo "  --force-update"
  echo ""

  if(! $PrintHelp) exit 1;
  echo $VERSION
  cat $0 | awk 'BEGIN{prt=0}{if(prt) print $0; if($1 == "BEGINHELP") prt = 1 }'
exit 1;

#---- Everything below here is printed out as part of help -----#
BEGINHELP

This is a script that places the surfaces on highres T1 (and maybe T2)
volumes based on an initial placement on a conformed volume. The idea
is that recon-all will be run up to and including the creation of
white.preaparc. This script is then run to generate the final white
and pial surfaces. After that, the final stages of recon-all can be
run (making sure to skip the creation of the final white and pial
surfaces). The recon-all.v6.hires script will do this automatically.

The pial surface tends to extend further out than in the normal
recon-all stream because the intensity normalization. mri_normalize is
run with different parameters (unless --norm-opts-rca is used) but
will also run differently because some of the parameters are set based
on the number of voxels, so differences in resolution effectively
cause these parameters to be different. The --copy-bias-from-conf
option will cause conf2hires to use the bias field computed from the
conformed (1mm) and makes the final surfaces much closer to that of
recon-all. I cannot remember why I did not use this as the default
except that the current set up might have worked better for HCP.