// ParallelTSHelper.cpp : Helper class for Parallel Tilt Series and task-oriented routines
//
//
// Copyright (C) 2007-2026 by the Regents of the University of
// Colorado.  See Copyright.txt for full notice of copyright and limitations.
//
// Author: Leo Crowder
//
#include "stdafx.h"
#include "SerialEM.h"
#include "ParallelTSDlg.h"
#include "NavHelper.h"
#include "EMscope.h"
#include "NavigatorDlg.h"
#include "ParticleTasks.h"
#include "ComplexTasks.h"
#include "CameraController.h"
#include "EMmontageController.h"
#include "SerialEMDoc.h"
#include "MultiShotDlg.h"
#include "ShiftManager.h"
#include "ShiftCalibrator.h"
#include "AutoTuning.h"
#include "FocusManager.h"
#include "TSController.h"
#include "SerialEMView.h"
#include "ParallelTSHelper.h"
#include "MacroProcessor.h"
#include ".\Utilities\SEMUtilities.h"
#include "Shared\b3dutil.h"
#include "Shared\icont.h"

#if defined(_DEBUG) && defined(_CRTDBG_MAP_ALLOC)
#define new DEBUG_NEW
#endif

CParallelTSHelper::CParallelTSHelper()
{
  mWinApp = (CSerialEMApp *)AfxGetApp();
  mShiftManager = mWinApp->mShiftManager;
  mScope = NULL;
  mCamera = NULL;
  mNavHelper = NULL;
  mParallelTSDlg = NULL;

  mISTargetIter = -1;
  mActionAtTarget = -1;
  mDoingISToTargets = false;
  mLastActionFailed = false;
  mPretilt = 0.f;
  mXpitch = 0.f;
  mStartIndex = 0;
  mAreaMapFileName = "";
  mAreaMapID = 0;
  mTargetMapFileName = "";
  mParTSitem = NULL;
  mSavedTSparamIndex = -1;
  mInitialStateSaved = false;
  mTiltDuringFit = 0.;
  mMagIndex = -1;
  mAreaMapMagInd = -1;
  mOldAdjustingXform.xpx = 0.f;
  mAdjustingXform.xpx = 0.f;
  mSavedRefocusThreshold = 0.f;
}

CParallelTSHelper::~CParallelTSHelper()
{

}

void CParallelTSHelper::Initialize(void)
{
  mScope = mWinApp->mScope;
  mCamera = mWinApp->mCamera;
  mNavHelper = mWinApp->mNavHelper;

  mParTSopts = mNavHelper->GetParTSOptions();
  mParallelTSDlg = mNavHelper->mParallelTSDlg;

  if (mWinApp->LowDoseMode()) {
    LowDoseParams *ldp = mWinApp->GetLowDoseParams() + RECORD_CONSET;
    mMagIndex = ldp->magIndex;
  } else {
    mMagIndex = mParTSopts->acqMagIndNonLD;
  }
}

// Busy function for task when doing IS
int CParallelTSHelper::ISToTargetsBusy()
{
  return (mCamera->CameraBusy() || mScope->StageBusy() ||
    mNavHelper->GetRealigning() || mWinApp->mFocusManager->DoingFocus()) ? 1 : 0;
}

// Clears data from vectors and resets members for a new run
void CParallelTSHelper::ClearTargets(bool autofocus)
{
  mStartIndex = 0;
  mAlignedToFirstISTarget = false;
  mInitialStateSaved = false;
  mISTargetPointIDs.clear();
  mSavedTargetIDs.clear();
  if (autofocus) {
    mISTargetSSX.clear();
    mISTargetSSY.clear();
    mISTargetDefocus.clear();
  } else {
    mISTargetISX.clear();
    mISTargetISY.clear();
    mPrevMapShiftX.clear();
    mPrevMapShiftY.clear();
    mPreviewMapIDs.clear();
    mPrevMapSectNums.clear();
    mPreRefineISX.clear();
    mPreRefineISY.clear();
  }
}

//Clear saved data, but keep target IDs and computed ISs so refining can be retried
void CParallelTSHelper::ClearSavedTargets()
{
  mStartIndex = 0;
  mAlignedToFirstISTarget = false;
  mInitialStateSaved = false;
  mSavedTargetIDs.clear();
  mPrevMapShiftX.clear();
  mPrevMapShiftY.clear();
  mPreviewMapIDs.clear();
  mPrevMapSectNums.clear();
  mPreRefineISX.clear();
  mPreRefineISY.clear();
}

// Begin procedure of image shifting to targets with given map IDs, and either getting the 
// autofocusing or taking Previews at each target
int CParallelTSHelper::StartShiftToTargets(IntVec targetMapIDs, int actionType)
{
  int err;
  CString mess;

  err = AppendNewTargets(targetMapIDs, mess);
  if (err) {
    AfxMessageBox(mess, MB_EXCLAME);
    return err;
  }
  mDoingISToTargets = true;
  mISTargetIter = 0;
  mAtTarget = false;
  mLastActionFailed = false;
  mActionAtTarget = actionType;
  mParTSopts = mNavHelper->GetParTSOptions();

  if (mWinApp->LowDoseMode()) {
    LowDoseParams *ldp = mWinApp->GetLowDoseParams() + RECORD_CONSET;
    mMagIndex = ldp->magIndex;
  } else {
    mMagIndex = mParTSopts->acqMagIndNonLD;
  }
  mNavHelper->SetParTSRefiningISMag(mMagIndex);

  if (mActionAtTarget == PARALLELTS_ACTION_AUTOFOCUS) {
    mess = "FITTING PLANE";
  } else if (mActionAtTarget == PARALLELTS_ACTION_PREVIEW) {
    mSavedMouseStage = mShiftManager->GetMouseMoveStage();
    mShiftManager->SetMouseMoveStage(false);
    mWinApp->UpdateWindowSettings();
    mess = "REFINING TARGETS";
  } else if (mActionAtTarget == PARALLELTS_ACTION_ADJUST) {
    mSavedMouseStage = mShiftManager->GetMouseMoveStage();
    mShiftManager->SetMouseMoveStage(false);
    mWinApp->UpdateWindowSettings();
    mess = "ADJUSTING IS";
  } else {
    AfxMessageBox("Invalid flag for Parallel Tilt Series action", MB_EXCLAME);
    return -1;
  }

  mWinApp->UpdateBufferWindows();
  mWinApp->SetStatusText(COMPLEX_PANE, mess);
  mWinApp->AddIdleTask(TASK_IS_TO_PARALLELTS_TARGET, 0, 0);
  return 0;
}

//Pauses the routine, does not restore original state so that more targets can be added.
void CParallelTSHelper::PauseParallelTSShift()
{
  mISTargetIter = -1;
  mDoingISToTargets = false;

  if (mActionAtTarget == PARALLELTS_ACTION_PREVIEW || 
    mActionAtTarget == PARALLELTS_ACTION_ADJUST) {
    if (mWinApp->mZoomedOverview)
      mWinApp->mZoomedOverview->CloseFrame();
    mParallelTSDlg->FinishRefineTargets(true);
  }
  mWinApp->SetStatusText(COMPLEX_PANE, "");
}

//Stops the routine, and restores original state of the scope.
void CParallelTSHelper::StopParallelTSShift(bool error)
{
  int numPoints;
  CMapDrawItem *item, *mapItem;
  
  mISTargetIter = -1;

  if (mInitialStateSaved) {
    if (mCamera->Acquiring()) {
      mCamera->SetImageShiftToRestore(mBaseISX, mBaseISY);
      if (mParTSopts->adjustBeamTilt)
        mCamera->SetBeamTiltToRestore(mBaseBeamTiltX, mBaseBeamTiltY);
      if (mAdjustBeamTilt)
        mCamera->SetAstigToRestore(mBaseAstigX, mBaseAstigY);
    } else {
      mScope->SetImageShift(mBaseISX, mBaseISY);
      if (mParTSopts->adjustBeamTilt)
        mScope->SetBeamTilt(mBaseBeamTiltX, mBaseBeamTiltY);
      if (mAdjustBeamTilt)
        mScope->SetObjectiveStigmator(mBaseAstigX, mBaseAstigY);
    }
    if (mActionAtTarget == PARALLELTS_ACTION_AUTOFOCUS)
      mWinApp->mFocusManager->SetRefocusThreshold(mSavedRefocusThreshold);
  }

  mNavHelper->SetParTSRefiningISMag(0);

  mDoingISToTargets = false;
  mAlignedToFirstISTarget = false;
  mDoNextShift = false;
  mInitialStateSaved = false;

  if (mActionAtTarget == PARALLELTS_ACTION_AUTOFOCUS) {
    numPoints = (int)mISTargetPointIDs.size();
    if (numPoints > 0) {
      item = mWinApp->mNavigator->FindItemWithMapID(
        mISTargetPointIDs[numPoints - 1], false);
      if (item && item->mDrawnOnMapID) {
        mapItem = mWinApp->mNavigator->FindItemWithMapID(item->mDrawnOnMapID);
        if (mapItem)
          mWinApp->mNavigator->DoLoadMap(true, mapItem, -1);
      }
    }

    if (!GetDebugOutput('N')) {
      for (int ind = 0; ind < numPoints; ind++) {
        item = mWinApp->mNavigator->FindItemWithMapID(mISTargetPointIDs[ind], false);
        if (item)
          mWinApp->mNavigator->ExternalDeleteItem(item,
            mWinApp->mNavigator->GetFoundItem());
      }
    }

    // Remove the current parallel tilt series item if fit points ran on different map
    if (mapItem && mapItem->mMapID != mAreaMapID) {
      mParTSitem = NULL;
    }
    mParallelTSDlg->FinishFitPlane();
    ClearTargets(true);
  } else if (mActionAtTarget == PARALLELTS_ACTION_PREVIEW || 
    mActionAtTarget == PARALLELTS_ACTION_ADJUST) {
    if (mWinApp->mZoomedOverview)
      mWinApp->mZoomedOverview->CloseFrame();
    mParallelTSDlg->FinishRefineTargets(!error);
    mShiftManager->SetMouseMoveStage(mSavedMouseStage);
    mWinApp->UpdateWindowSettings();
    if (error) {
      DeleteTargetMapsFromNav();
      ClearSavedTargets();
    }
  }
  mWinApp->SetStatusText(COMPLEX_PANE, "");
}

// Applies image shift to get to the next target. 
// If returned error code is < 0, the whole procedure should be aborted. If > 0, then just
// the given target should be skippped.
int CParallelTSHelper::ISToNextTarget(int targetID, CString &err)
{
  double ISX, ISY, delX, delY, delBTX = 0., delBTY = 0., delAstigX = 0., delAstigY = 0.;
  float delay, transISX, transISY;
  int BTdelay, fromMagInd, mapID;
  CMapDrawItem *item = mWinApp->mNavigator->FindItemWithMapID(targetID, false);
  int navInd;
  ScaleMat st2is;
  bool doBacklash = mScope->GetAdjustForISSkipBacklash() <= 0;
  ComaVsISCalib *comaVsIS = mWinApp->mAutoTuning->GetComaVsIScal();
  ParallelTSOptions *parTSopt = mNavHelper->GetParTSOptions();
  CString mess;

  if (!item) {
    err.Format("The given target item no longer exists");
    return -3;
  }
  
  navInd = mWinApp->mNavigator->GetFoundItem();
  mapID = item->mDrawnOnMapID;
  CMapDrawItem *mapItem = mWinApp->mNavigator->FindItemWithMapID(mapID);
  if (!mapItem) {
    err.Format("The map on which items were drawn no longer exists");
    return -2;
  }

  fromMagInd = mAreaMapMagInd;
  if (fromMagInd < 0) {
    fromMagInd = mapItem->mMapMagInd;
  }

  st2is = MatMul(mShiftManager->StageToCamera(mWinApp->GetCurrentCamera(),
    fromMagInd), mShiftManager->CameraToIS(fromMagInd));
  SEMTrace('N', "stage to IS: mag index %d xpx %.3f xpy %.3f ypx %.3f ypy %.3f",
    fromMagInd, st2is.xpx, st2is.xpy, st2is.ypx, st2is.ypy);
  if (!st2is.xpx) {
    err.Format("There is no calibration to get from stage to IS coordinates at the "
      "given mag");
    return -1;
  }
  ApplyScaleMatrix(st2is, item->mStageX - mCenterStageX,
    item->mStageY - mCenterStageY, ISX, ISY);
  mShiftManager->TransferGeneralIS(fromMagInd, ISX, ISY, mMagIndex, delX, delY);

  mOldAdjustingXform.xpx = 0.f;
  if ((mActionAtTarget == PARALLELTS_ACTION_PREVIEW || 
    mActionAtTarget == PARALLELTS_ACTION_ADJUST) && mParTSopts->applyAdjustingXform) {
    if (CanAdjustISVectors(mAreaMapMagInd, false, mess)) {
      AdjustXformData *adjustData =
        mNavHelper->GetNearestAdjustingXform(mAreaMapMagInd);
      SEMTrace('N', "adj xform xpx=%.4f, xpy=%.4f, ypx=%.4f, ypy=%.4f",
        adjustData->adjustingXform.xpx, adjustData->adjustingXform.xpy,
        adjustData->adjustingXform.ypx, adjustData->adjustingXform.ypy);

      //Save applied adjusting transform to factor into calculation of new one
      mOldAdjustingXform = adjustData->adjustingXform;

      ApplyScaleMatrix(adjustData->adjustingXform, delX, delY, transISX, transISY);
      SEMTrace('N', "IS transformed from (%.4f, %.4f) to (%.4f, %.4f)",
        delX, delY, transISX, transISY);
      delX = transISX;
      delY = transISY;
    } else {
      SEMTrace('N', "Adjusting transform not applied: %s", mess);
    }
  }

  ISX = mCenterISX + delX;
  ISY = mCenterISY + delY;

  if (!mShiftManager->ImageShiftIsOK(ISX, ISY, FALSE)) {
    err.Format("Point at navigator index %d is beyond image shift limit and will "
      "be skipped", navInd);
    mWinApp->mNavigator->ExternalDeleteItem(item, navInd);
    return 1;
  }

  mCurISTargetItem = item;

  mScope->GetImageShift(mLastISX, mLastISY);

  // Compute the delay for the amount that IS will change
  delay = mShiftManager->ComputeISDelay(ISX - mLastISX, ISY - mLastISY);
  if (parTSopt->extraDelayFactor > 0.)
    delay *= parTSopt->extraDelayFactor;

  if (mParTSopts->adjustBeamTilt) {
    BTdelay = mWinApp->mAutoTuning->GetBacklashDelay();
    mWinApp->mParticleTasks->GetBTandAstigAdjustment(delX, delY, delBTX, delBTY,
      delAstigX, delAstigY, FALSE, GetDebugOutput('I'));
  }

  mScope->SetImageShift(ISX, ISY);
  mShiftManager->SetISTimeOut(delay);

  if (mParTSopts->adjustBeamTilt) {
    mWinApp->mAutoTuning->BacklashedBeamTilt(mCenterBeamTiltX + delBTX,
      mCenterBeamTiltY + delBTY, doBacklash);
    if (mAdjustBeamTilt)
      mWinApp->mAutoTuning->BacklashedStigmator(mCenterAstigX + delAstigX,
        mCenterAstigY + delAstigY, doBacklash);
  }

  mLastISX = ISX;
  mLastISY = ISY;

  return 0;
}

// Next task in the IS to target routine
void CParallelTSHelper::ISToTargetNextTask(int param)
{
  int index, mapID;
  float residual, windowSize;
  float focusLim = -20.;
  float shiftX, shiftY;
  int sizeX, sizeY;
  float FOVchange, FOVchangeThresh = 0.1f;
  bool mapSaved, openNewFile = false, openOldFile = false;
  EMimageBuffer *imBufs = mWinApp->GetImBufs();
  CString mess, str;
  ScaleMat st2ss = MatInv(mShiftManager->SpecimenToStage(1., 1.));
  NavAlignParams *alignParams = mNavHelper->GetNavAlignParams();
  ComaVsISCalib *comaVsIS = mWinApp->mAutoTuning->GetComaVsIScal();
  bool canAdjustIS = mShiftManager->GetFocusISCals()->GetSize() > 0 &&
    mShiftManager->GetFocusMagCals()->GetSize() > 0;

  if (!mDoingISToTargets || mISTargetIter < 0)
    return;

  if (mActionAtTarget == PARALLELTS_ACTION_AUTOFOCUS) {
    mess = "FITTING PLANE";
  } else if (mActionAtTarget == PARALLELTS_ACTION_PREVIEW) {
    mess = "REFINING TARGETS";
  } else if (mActionAtTarget == PARALLELTS_ACTION_ADJUST) {
    mess = "ADJUSTING IS";
  }
  mWinApp->SetStatusText(COMPLEX_PANE, mess);
  mParallelTSDlg->SetInstructionLine();

  mLastActionFailed = false;
  mapID = mISTargetPointIDs[mStartIndex + mISTargetIter];
  mCurISTargetItem = mWinApp->mNavigator->FindItemWithMapID(mapID, false);

  if (mISTargetIter == 0 && !mAlignedToFirstISTarget) { 
    if (!mCurISTargetItem) {
      StopParallelTSShift();
      return;
    }

    //Realign to first IS Target
    if (mNavHelper->RealignToItem(mCurISTargetItem, 0, 
      alignParams->resetISthresh, alignParams->maxNumResetIS, 0, 0,
      mParTSopts->extractVirtPrevs == 1 ? PREVIEW_CONSET : -1)) {
      AfxMessageBox("Realign to first IS target failed");
      mLastActionFailed = true;
      StopParallelTSShift(true);
      ClearTargets(true);
      return;
    }
    mAlignedToFirstISTarget = true;
    mDoNextShift = false;
    mWinApp->AddIdleTask(TASK_IS_TO_PARALLELTS_TARGET, 0, 0);
    return;
  }

  if (mAtTarget && mCurISTargetItem) {

    if (mActionAtTarget == PARALLELTS_ACTION_PREVIEW) {

      if (mParTSopts->extractVirtPrevs == 0) {
        imBufs->mImage->getShifts(shiftX, shiftY);
        imBufs->mImage->getSize(sizeX, sizeY);
        FOVchange = (1 - (sizeX - fabs(shiftX)) * (sizeY - fabs(shiftY)) / 
          (float)(sizeX * sizeY));

        if (FOVchange > FOVchangeThresh) {
          mCamera->InitiateCapture(PREVIEW_CONSET);
          mWinApp->AddIdleTask(TASK_IS_TO_PARALLELTS_TARGET, 0, 0);
          return;
        }
        mess = "";
        index = SaveTargetMap(mess, mapSaved);
        if (index) {
          AfxMessageBox(mess, MB_EXCLAME);
          StopParallelTSShift();
          return;
        }

        if (!mapSaved) {
          if (!mess.IsEmpty()) {
            str.Format("Map was not saved: %s", mess);
            AfxMessageBox(str, MB_EXCLAME);
          }
          return;
        }
      }
    }

    if (!mInitialStateSaved && SaveInitialState(mess)) {
      AfxMessageBox(mess, MB_EXCLAME);
      StopParallelTSShift(true);
      ClearTargets(true);
      return;
    }
 
    bool skipSave = mActionAtTarget == PARALLELTS_ACTION_AUTOFOCUS && mISTargetIter == 0;

    if ((mSavedTargetIDs.size() == 0 && (mActionAtTarget == PARALLELTS_ACTION_PREVIEW
      || mActionAtTarget == PARALLELTS_ACTION_ADJUST)) || skipSave) {
      mCenterStageX = mCurISTargetItem->mStageX;
      mCenterStageY = mCurISTargetItem->mStageY;
      mScope->GetImageShift(mCenterISX, mCenterISY);
      mTiltDuringFit = mScope->GetTiltAngle();
    }
    
    if (!skipSave && SaveTarget(mess)) {
      AfxMessageBox(mess, MB_EXCLAME);
      StopParallelTSShift(true);
      ClearSavedTargets();
      return;
    }

  } else if (mAtTarget && !mCurISTargetItem) {

    //Center item was deleted but none have been saved, cancel refining
    if (mSavedTargetIDs.size() == 0) {
      StopParallelTSShift(true);
      ClearTargets(true);
      return;
    }
    if (!mInitialStateSaved && SaveInitialState(mess)) {
      AfxMessageBox(mess, MB_EXCLAME);
      StopParallelTSShift(true);
      ClearTargets(true);
      return;
    }
    mDoNextShift = true;
  }

  if (mDoNextShift) {
    mAtTarget = false;

    if (mStartIndex + mISTargetIter >= (int)mISTargetPointIDs.size() - 1) {
      if (mActionAtTarget == PARALLELTS_ACTION_AUTOFOCUS) {

        // Run the least squares fit to get pretilt and X pitch angle
        index = FitPlane(mPretilt, mXpitch, residual, mess);
        if (index) {
          str.Format("Plane fit on points failed: %s", mess);
          AfxMessageBox(str, MB_EXCLAME);
        }
        StopParallelTSShift();
        ClearTargets(true);
        if (!index && mParallelTSDlg->IsOpen()) {
          mParallelTSDlg->UpdatePlaneParams(mPretilt, mXpitch);
          UpdateSpecAngles(mPretilt, mXpitch);
        }
      } else if (mActionAtTarget == PARALLELTS_ACTION_PREVIEW || 
        mActionAtTarget == PARALLELTS_ACTION_ADJUST) {
         PauseParallelTSShift();
      }
      return;
    }

    // Increment the target iteration to get next target
    //If some targets were done on a previous separate run, skip to the first new target
    mISTargetIter++;
    mapID = mISTargetPointIDs[mStartIndex + mISTargetIter];

    // Do the image shift to next target
    index = ISToNextTarget(mapID, mess);

    // If returned error code is < 0, the whole procedure should be aborted. 
    // If > 0, then just the given target should be skippped.
    if (index) {
      AfxMessageBox(mess, MB_EXCLAME);
      if (index < 0) {
        StopParallelTSShift(true);
        ClearSavedTargets();
        return;
      } else if (index > 0) {
        mLastActionFailed = true;
      }
    }

    // If failed at this point, skip action and shift to next point if continuing on
    mDoNextShift = mLastActionFailed;
  }

  if (!mDoNextShift) {

    //Queue to save this target and shift to next target on next iteration
    mAtTarget = true;
    mDoNextShift = true;
    
    // Go to the desired mag or low dose area before recording any values
    if (mWinApp->LowDoseMode() && mScope->GetLowDoseArea() != RECORD_CONSET) {
      mScope->GotoLowDoseArea(RECORD_CONSET);
    } else if (!mWinApp->LowDoseMode() && mScope->GetMagIndex() != mMagIndex &&
      !mScope->SetMagIndex(mMagIndex)) {
      AfxMessageBox("Failed to switch to the desired magnification", MB_EXCLAME);
      StopParallelTSShift(true);
      ClearSavedTargets();
      return;
    }

    //Get Image shift before next iteration
    mScope->GetImageShift(mLastISX, mLastISY);

    if (mActionAtTarget == PARALLELTS_ACTION_AUTOFOCUS) {
      if (mISTargetIter)  
        mWinApp->mFocusManager->AutoFocusStart(1, -2);
    } else if (mActionAtTarget == PARALLELTS_ACTION_PREVIEW ||
      mActionAtTarget == PARALLELTS_ACTION_ADJUST) {
      if (mActionAtTarget == PARALLELTS_ACTION_ADJUST || 
        !(mParTSopts->flags & PTSFLAG_SKIP_REFINE)) {
        CameraParameters *camParam = mWinApp->GetActiveCamParam();
        windowSize = 1.5f * (float)camParam->sizeX *
          mShiftManager->GetPixelSize(mWinApp->GetCurrentCamera(), mMagIndex);
        CMapDrawItem *mapItem = mWinApp->mNavigator->FindItemWithMapID(mAreaMapID);
        if (!mapItem) {
          AfxMessageBox("The Parallel Tilt Series area map was deleted", MB_EXCLAME);
          StopParallelTSShift(true);
          ClearSavedTargets();
          return;
        }
        mWinApp->mNavigator->DoLoadMap(true, mapItem, -1);
        EMimageBuffer *mapBuf = imBufs + mWinApp->mBufferManager->GetBufToReadInto();
        mWinApp->mMainView->GetItemImageCoords(mapBuf, mCurISTargetItem, shiftX, shiftY);
        mWinApp->CopyBufferToZoomedView(mapBuf, windowSize);
        mWinApp->CenterZoomedViewAtPoint(shiftX, shiftY);
        mWinApp->RestoreViewFocus();
      }
      mCamera->InitiateCapture(PREVIEW_CONSET);

      // Stop iterations to allow user to refine IS, unless they want to skip that
      if (!(mActionAtTarget == PARALLELTS_ACTION_PREVIEW &&
        mParTSopts->flags & PTSFLAG_SKIP_REFINE)) {
        mess.Format("Adjust image shift to align to the target. Click \"%s\" or press S "
          "key to save and go on.", mParallelTSDlg->GetSaveBtnText());
        mParallelTSDlg->SetInstructionLine(mess);
        return;
      }
    }
  }

  mWinApp->AddIdleTask(TASK_IS_TO_PARALLELTS_TARGET, 0, 0);
}

// Save a map for the parallel tilt series area
int CParallelTSHelper::SaveAreaMap(CString &err)
{
  int index;
  bool wrongFile, saveable, openNewFile = false, openOldFile = false;
  CFile *cfile;
  CMapDrawItem *item;
  EMimageBuffer *imBuf = mWinApp->mMainView->GetActiveImBuf();
  EMimageBuffer *imBufs = mWinApp->GetImBufs();
  CString bufLetter;

  // Current active image buffer is a loaded map, so map already exists
  item = mWinApp->mNavigator->FindItemWithMapID(imBuf->mMapID);
  if (!item) {
    if (mWinApp->Montaging() && imBufs[1].mCaptured == BUFFER_MONTAGE_OVERVIEW) {
      imBuf = &imBufs[1];
      saveable = mWinApp->mStoreMRC != NULL;
      bufLetter = "B";
    } else {
      imBuf = &imBufs[0];
      saveable = mWinApp->mStoreMRC && 
        mWinApp->mBufferManager->IsBufferSavable(imBuf, mWinApp->mStoreMRC);
      bufLetter = "A";
    }
    
    if (imBuf->mMapID > 0 && mWinApp->mNavigator->FindItemWithMapID(imBuf->mMapID)) {
      err.Format("The image in buffer %s has already been saved to a map", bufLetter);
      return -1;
    }

    //If no open file, or open file is unusable, do not use it
    index = mWinApp->mDocWnd->StoreIndexFromName(mTargetMapFileName);
    wrongFile = !mTargetMapFileName.IsEmpty() && (index >= 0 ||
      index == mWinApp->mDocWnd->GetCurrentStore());

    if (wrongFile || !saveable) {

      //If Area map is defined and exists, switch to or open it. If not, open new file
      if (!mAreaMapFileName.IsEmpty() && UtilFileExists(mAreaMapFileName) != 0) {
        index = mWinApp->mDocWnd->StoreIndexFromName(mAreaMapFileName);
        if (index >= 0)
          mWinApp->mDocWnd->SetCurrentStore(index);
        else {
          index = mWinApp->mDocWnd->OpenOldMrcCFile(&cfile, mAreaMapFileName, false);
          if (index == MRC_OPEN_NOERR || index == MRC_OPEN_ADOC || index == MRC_OPEN_HDF)
            index = mWinApp->mDocWnd->OpenOldFile(cfile, mAreaMapFileName, index, true);
          if (index != MRC_OPEN_NOERR) {
            err.Format("Error opening old area map file");
            return 1;
          }
        }
      } else {
        if (mWinApp->mDocWnd->DoOpenNewFile()) {
          err.Format("New file was not opened");
          return 2;
        }
      }
    }

    if (imBuf->GetSaveCopyFlag() >= 0 && !mWinApp->Montaging())
      mWinApp->mDocWnd->SaveRegularBuffer();

    if (mWinApp->mNavigator->NewMap()) {
      err.Format("Error making a new area map");
      return 3;
    }

    index = mWinApp->mDocWnd->StoreIndexFromName(mAreaMapFileName);
    if (index >= 0) {
      mWinApp->mDocWnd->SetCurrentStore(index);
      mWinApp->mDocWnd->DoCloseFile();
    }

    item = mWinApp->mNavigator->GetCurrentItem();
  }
  mAreaMapID = item->mMapID;
  mAreaMapMagInd = item->mMapMagInd;
  mAreaMapFileName = item->mMapFile;
  mMappingTilt = item->mMapTiltAngle;
  return 0;
}

// Save a map centered on a target.
int CParallelTSHelper::SaveTargetMap(CString &err, bool &saved)
{
  int areaStore, tgtStore, store, index, numMaps;
  bool areaFileOpen, tgtFileOpen, fileExists;
  float shiftX, shiftY;
  CFile *cfile;
  CMapDrawItem *item;
  EMimageBuffer *imBuf = mWinApp->GetImBufs();
  saved = false;

  // check that mag of current buffer image matches acquire mag
  if (imBuf->mMagInd != mMagIndex) {
    err.Format("Target Maps must match the magnification that will be used to "
      "acquire images");
    return 0;
  }

  fileExists = !mTargetMapFileName.IsEmpty() && UtilFileExists(mTargetMapFileName);
  tgtStore = mWinApp->mDocWnd->StoreIndexFromName(mTargetMapFileName);
  numMaps = (int)mPreviewMapIDs.size();

  //If no open file, or open file is unusable, do not use it
  areaStore = mWinApp->mDocWnd->StoreIndexFromName(mAreaMapFileName);
  store = mWinApp->mDocWnd->GetCurrentStore();
  areaFileOpen = (areaStore >= 0 && areaStore == store);
  tgtFileOpen = (tgtStore >= 0 && tgtStore == store);

  if (tgtStore >= 0 && tgtStore != store) {
    mWinApp->mDocWnd->SetCurrentStore(tgtStore);
  } else if (!tgtFileOpen && (!mWinApp->mStoreMRC || areaFileOpen || numMaps > 0 ||
    !mWinApp->mBufferManager->IsBufferSavable(imBuf, mWinApp->mStoreMRC))) {
    if (fileExists) {
      index = mWinApp->mDocWnd->OpenOldMrcCFile(&cfile, mTargetMapFileName, false);
      if (index == MRC_OPEN_NOERR || index == MRC_OPEN_ADOC || index == MRC_OPEN_HDF)
        index = mWinApp->mDocWnd->OpenOldFile(cfile, mTargetMapFileName, index, true);
      if (index != MRC_OPEN_NOERR) {
        err.Format("Error opening map file");
        return 2;
      }
    } else if (numMaps == 0) {
      if (mWinApp->mDocWnd->DoOpenNewFile())
        return 0;
    } else if (numMaps > 0) {
    //If some preview maps were saved but map file was deleted, throw error
      err.Format("The Preview map file was deleted");
      return 3;
    }
  }

  mWinApp->mDocWnd->SaveRegularBuffer();

  if (mWinApp->mNavigator->NewMap()) {
    err.Format("Error making a new Preview map");
    return 0;
  }

  item = mWinApp->mNavigator->GetCurrentItem();
  if (numMaps == 0)
    mTargetMapFileName = item->mMapFile;
  mPreviewMapIDs.push_back(item->mMapID);
  mPrevMapSectNums.push_back(item->mMapSection);
  imBuf->mImage->getShifts(shiftX, shiftY);
  mPrevMapShiftX.push_back(shiftX);
  mPrevMapShiftY.push_back(shiftY);

  saved = true;
  return 0;
}

int CParallelTSHelper::SaveMap(CString &err)
{
  int index;
  bool wrongFile, saveable, openNewFile = false, openOldFile = false;
  EMimageBuffer *imBuf;
  EMimageBuffer *imBufs = mWinApp->GetImBufs();
  CString bufLetter;

  if (mWinApp->Montaging() && imBufs[1].mCaptured == BUFFER_MONTAGE_OVERVIEW) {
    imBuf = &imBufs[1];
    saveable = mWinApp->mStoreMRC != NULL;
    bufLetter = "B";
  }
  else {
    imBuf = &imBufs[0];
    saveable = mWinApp->mStoreMRC &&
      mWinApp->mBufferManager->IsBufferSavable(imBuf, mWinApp->mStoreMRC);
    bufLetter = "A";
  }

  if (imBuf->mMapID > 0 && mWinApp->mNavigator->FindItemWithMapID(imBuf->mMapID)) {
    err.Format("The image in buffer %s has already been saved to a map", bufLetter);
    return -1;
  }

  //If no open file, or open file is unusable, do not use it
  index = mWinApp->mDocWnd->StoreIndexFromName(mTargetMapFileName);
  wrongFile = !mTargetMapFileName.IsEmpty() && (index >= 0 ||
    index == mWinApp->mDocWnd->GetCurrentStore());

  if (wrongFile || !saveable) {
    if (mWinApp->mDocWnd->DoOpenNewFile()) {
      err.Format("New file was not opened");
      return 2;
    }
  }

  if (imBuf->GetSaveCopyFlag() >= 0 && !mWinApp->Montaging())
    mWinApp->mDocWnd->SaveRegularBuffer();

  if (mWinApp->mNavigator->NewMap()) {
    err.Format("Error making a new map");
    return 3;
  }
  return 0;
}

// Saves the state of the microscope before beginning the routine
int CParallelTSHelper::SaveInitialState(CString &err)
{
  int area;
  bool canAdjustIS = mShiftManager->GetFocusISCals()->GetSize() > 0 &&
    mShiftManager->GetFocusMagCals()->GetSize() > 0;
  double stageZ;
  float focusLim = -20.;
  EMimageBuffer *imBufs = mWinApp->GetImBufs();
  ComaVsISCalib *comaVsIS = mWinApp->mAutoTuning->GetComaVsIScal();

  area = mScope->GetLowDoseArea();

  // First time, record stage position and check for defocus

  // Go to the desired mag or low dose area before recording any values
  if (mWinApp->LowDoseMode()) {
    mScope->GotoLowDoseArea(RECORD_CONSET);
  } else {
    if (mScope->GetMagIndex() != mMagIndex &&
      !mScope->SetMagIndex(mMagIndex)) {
      err.Format("Failed to switch to the desired magnification");
      return 1;
    }
  }

  mScope->GetStagePosition(mBaseStageX, mBaseStageY, stageZ);
  if ((mActionAtTarget == PARALLELTS_ACTION_PREVIEW
    || mActionAtTarget == PARALLELTS_ACTION_ADJUST) &&
    !canAdjustIS && ((imBufs->mLowDoseArea && (imBufs->mConSetUsed == VIEW_CONSET ||
    imBufs->mConSetUsed == SEARCH_CONSET) && imBufs->mViewDefocus < focusLim) ||
    (mWinApp->LowDoseMode() && IS_AREA_VIEW_OR_SEARCH(area) &&
      mScope->GetLDViewDefocus(area) < focusLim))) {
    err.Format("It appears that the scope is in the View Low Dose\n"
      "area with a View defocus offset bigger than %.0f microns.\n\n"
      "This procedure should be run closer to focus.\n"
      "Press Abort to end it, or continue if you know what you are doing.", focusLim);
    AfxMessageBox(err, MB_EXCLAME);
  }

  mScope->GetImageShift(mBaseISX, mBaseISY);

  mAdjustBeamTilt = mParTSopts->adjustBeamTilt && comaVsIS->astigMat.xpx != 0. &&
    mNavHelper->GetSkipAstigAdjustment() <= 0;

  if (mParTSopts->adjustBeamTilt) {
    if (!mScope->GetBeamTilt(mBaseBeamTiltX, mBaseBeamTiltY)) {
      err.Format("Failed to get beam tilt.");
      return 2;
    }
    mCenterBeamTiltX = mBaseBeamTiltX;
    mCenterBeamTiltY = mBaseBeamTiltY;
    
    if (mAdjustBeamTilt) {
      if (!mScope->GetObjectiveStigmator(mBaseAstigX, mBaseAstigY)) {
        err.Format("Failed to get objective stigmator.");
        return 3;
      }
      mCenterAstigX = mBaseAstigX;
      mCenterAstigY = mBaseAstigY;
    }
  }

  if (mActionAtTarget == PARALLELTS_ACTION_AUTOFOCUS) {
    mSavedRefocusThreshold = mWinApp->mFocusManager->GetRefocusThreshold();
    mWinApp->mFocusManager->SetRefocusThreshold(1.);
  }
  
  mInitialStateSaved = true;
  return 0;
}

//Save info about current target. For plane fit, specimen coords and defocus. For target
//refinement, image shifts and Preview map info.
int CParallelTSHelper::SaveTarget(CString &err)
{
  int index, index2, mapID;
  float defocus, SX, SY, factor;
  CString mess;
  ScaleMat mat;
  NavAlignParams *alignParams = mNavHelper->GetNavAlignParams();
  bool canAdjustIS = mShiftManager->GetFocusISCals()->GetSize() > 0 &&
    mShiftManager->GetFocusMagCals()->GetSize() > 0;
  double ISX, ISY, stageX, stageY, stageZ, delX, delY;
  float ISlimit = 2.f * mWinApp->mShiftCalibrator->GetCalISOstageLimit();
  float focusLim = -20.;
  EMimageBuffer *imBufs = mWinApp->GetImBufs();
  ComaVsISCalib *comaVsIS = mWinApp->mAutoTuning->GetComaVsIScal();

  mapID = mCurISTargetItem->mMapID;
  
  mScope->GetStagePosition(stageX, stageY, stageZ);
  if (fabs(stageX - mBaseStageX) > ISlimit ||
    fabs(stageY - mBaseStageY) > ISlimit) {
    err.Format("The stage appears to have moved by %.3f in X and %.3f in Y\n"
      "(more than ISoffsetCalStageLimit = %f).\n\n"
      "There must be no stage movement during this procedure.",
      stageX - mBaseStageX, stageY - mBaseStageY, ISlimit);
    return 1;
  }
  
  // If taking Previews, save the image shift at the target
  if (mActionAtTarget == PARALLELTS_ACTION_PREVIEW || 
    mActionAtTarget == PARALLELTS_ACTION_ADJUST) {

    // Get the IS and save it
    mScope->GetImageShift(ISX, ISY);
    delX = ISX - mLastISX;
    delY = ISY - mLastISY;
    
    if (mActionAtTarget == PARALLELTS_ACTION_ADJUST) {
      SEMTrace('N', "Image shift adjusted by %.4f, %.4f", delX, delY);
    }

    mPreRefineISX.push_back(mLastISX);
    mPreRefineISY.push_back(mLastISY);
    mLastISX = ISX;
    mLastISY = ISY;
    
    mISTargetISX.push_back(ISX);
    mISTargetISY.push_back(ISY);
    mSavedTargetIDs.push_back(mCurISTargetItem->mMapID);
    

    //If center point, update stage position
    if (mSavedTargetIDs.size() == 1 && 
      (B3DABS(delX) > ISlimit || B3DABS(delY) > ISlimit)) {

      //Convert IS change to stage shift
      mat = MatMul(mShiftManager->IStoCamera(mMagIndex), 
        MatInv(mShiftManager->StageToCamera(mWinApp->GetCurrentCamera(), mMagIndex)));
      mShiftManager->ApplyScaleMatrix(mat, delX, delY, SX, SY);
      SEMTrace('N', "Shifted stage coords of PTS center pt by %.3f, %.3f", SX, SY);

      mCenterStageX += SX;
      mCenterStageY += SY;
      mCurISTargetItem->mStageX = mCenterStageX;
      mCurISTargetItem->mStageY = mCenterStageY;
      mCurISTargetItem->mPtX[0] = mCenterStageX;
      mCurISTargetItem->mPtY[0] = mCenterStageY;
      mWinApp->mNavigator->FindItemWithMapID(mapID, false);
      mWinApp->mNavigator->UpdateListString(mWinApp->mNavigator->GetFoundItem());
      mWinApp->mNavigator->Redraw();
    }

  } else if (mActionAtTarget == PARALLELTS_ACTION_AUTOFOCUS) {

    // If doing autofocus at targets, get the defocus
    index = mWinApp->mFocusManager->GetLastFailed() ? -1 : 0;
    index2 = mWinApp->mFocusManager->GetLastAborted();
    if (index2)
      index = index2;
    if (index) {
      PrintfToLog("WARNING: Autofocus failed on point at index %d with error type %d",
        mapID, index);
      mLastActionFailed = true;

    } else {
      defocus = (float)mScope->GetDefocus();
      factor = mWinApp->mFocusManager->AdjustedDefocusZFactor(-1);
      mISTargetDefocus.push_back(-defocus / factor);

      if (mISTargetDefocus.size() == 1) {
        mBaseDefocus = defocus;
      } else {
        mScope->SetDefocus(mBaseDefocus);
      }

      mat = MatInv(mShiftManager->SpecimenToStage(1., 1.));
      mShiftManager->ApplyScaleMatrix(mat, mCurISTargetItem->mStageX,
        mCurISTargetItem->mStageY, SX, SY);
      mISTargetSSX.push_back(SX);
      mISTargetSSY.push_back(SY);
      mSavedTargetIDs.push_back(mCurISTargetItem->mMapID);
    }

    //Get Image shift since we will go straight into shift to next target
    mScope->GetImageShift(mLastISX, mLastISY);
  }
  return 0;
}

// Measure z-height at navigator items in the specified group, and perform a least 
// squares fit to a plane to determine the pretilt and X pitch angle.
int CParallelTSHelper::FitPlane(float &pretilt, float &xPitch, float &residual,
  CString &mess)
{
  int numPoints, numToDrop;
  FloatVec residuals;
  float alpha;
  float dev;
  float xCoef, yCoef, zIntercept;
  int i;

  if ((int)mISTargetDefocus.size() < MIN_NUM_POINTS_TO_FIT_PLANE) {
    mess.Format("Defocus was measured at %d out of %d points, "
      "but %d are required to fit to plane",
      (int)mISTargetDefocus.size(), (int)mISTargetPointIDs.size() - 1, 
      MIN_NUM_POINTS_TO_FIT_PLANE);
    return 3;
  }

  //Run least squares fit on vectors of x,y,z coords
  numPoints = (int)mISTargetDefocus.size();
  lsFit2(&mISTargetSSX[0], &mISTargetSSY[0], &mISTargetDefocus[0], numPoints, &xCoef,
    &yCoef, &zIntercept);

  for (i = 0; i < numPoints; i++) {
    dev = mISTargetDefocus[i] -
      (mISTargetSSX[i] * xCoef + mISTargetSSY[i] * yCoef + zIntercept);
    residuals.push_back(B3DABS(dev));
    if (GetDebugOutput('N')) {
      PrintfToLog("%d  %.4f  %.4f  %.4f  %.4f", i, mISTargetSSX[i], mISTargetSSY[i],
        mISTargetDefocus[i], dev);
    }
  }

  if (numPoints > MIN_NUM_POINTS_TO_FIT_PLANE) {
    FloatVec outliers, outlierResid;
    float elimMin = 0.01f; //TODO what should this be?
    FindOutliersInResultList(residuals, elimMin, 1, outliers);
    IntVec dropInd;
    for (i = 0; i < numPoints; i++) {
      if (outliers[i] > 0) {
        outlierResid.push_back(residuals[i]);
        dropInd.push_back(i);
      }
    }

    //Remove outliers and redo the fit
    numToDrop = B3DMIN((int)dropInd.size(), numPoints - MIN_NUM_POINTS_TO_FIT_PLANE);

    if (numToDrop > 0) {
      rsSortIndexedFloats(&outlierResid[0], &dropInd[0], numToDrop);
      dropInd.resize(numToDrop);

      //Resort into increasing index order to easily drop elements starting from the end
      std::sort(dropInd.begin(), dropInd.end());
      for (i = numToDrop - 1; i >= 0; i--) {
        mISTargetSSX.erase(mISTargetSSX.begin() + dropInd[i]);
        mISTargetSSY.erase(mISTargetSSY.begin() + dropInd[i]);
        mISTargetDefocus.erase(mISTargetDefocus.begin() + dropInd[i]);
      }

      if (AssessPtsToFitPlane(mISTargetSSX, mISTargetSSY,
        mISTargetDefocus, mess)) {
        return 2;
      }

      lsFit2(&mISTargetSSX[0], &mISTargetSSY[0], &mISTargetDefocus[0],
        (int)mISTargetDefocus.size(), &xCoef, &yCoef, &zIntercept);
    }

  }

  //Measured angle from fit, which is current stage tilt plus specimen tilt
  alpha = atan(yCoef); 

  //Pretilt is the tilt angle that will make the specimen level
  pretilt = (float)mTiltDuringFit - alpha / (float)RADIANS_PER_DEGREE; 
  xPitch = -atan(xCoef / (cos(alpha) - yCoef * sin(alpha))) / (float)RADIANS_PER_DEGREE;

  if (GetDebugOutput('N')) {
    PrintfToLog("xcoef =  %.4f, ycoef =  %.4f, zintercept =  %.4f", xCoef, yCoef, 
      zIntercept);
    PrintfToLog("alpha =  %.4f, pretilt = %.4f, xPitch = %.4f", alpha, pretilt, xPitch);
  }

  return 0;
}

// Find the long axis and aspect ratio of the convex hull of a set of points. 
// This is used to determine if a set of points sufficiently span 2d space.
void CParallelTSHelper::ConvexHullLongAxis(FloatVec ptsX, FloatVec ptsY, 
  float *aspectRatio, float *longAxis, float anglePrecision)
{
  int numPoints = B3DMIN((int)ptsX.size(), (int)ptsY.size());
  FloatVec xHull, yHull;
  int nHull;
  float xcen, ycen;
  Icont *cont;

  xHull.resize(numPoints);
  yHull.resize(numPoints);
  convexBound(&ptsX[0], &ptsY[0], numPoints, 0., 0., &xHull[0], &yHull[0], &nHull, &xcen,
    &ycen, numPoints);

  cont = imodContourNew();
  cont->pts = B3DMALLOC(Ipoint, nHull);
  if (cont->pts) {
    cont->psize = nHull;
    for (int ind = 0; ind < nHull; ind++) {
      cont->pts[ind].x = xHull[ind];
      cont->pts[ind].y = yHull[ind];
      cont->pts[ind].z = 0.;
    }
  }
  imodContourLongAxis(cont, anglePrecision, aspectRatio, longAxis);
  imodContourDelete(cont);
}

// Assess if the provided group of points are sufficient for a 2D least squares fit
int CParallelTSHelper::AssessPtsToFitPlane(FloatVec &ptsX, FloatVec &ptsY, FloatVec &ptsZ, 
  CString &mess)
{
  int numPoints, centerPtInd;
  FloatVec finalPtsX, finalPtsY, finalPtsZ;
  IntVec indexVec, sortedIndexVec;
  float dist, cenX = 0.f, cenY = 0.f;
  float ratio, longAxis, distToCen;

  float minSpanOfPlanePts = 1.f; //TODO get this from parameters or something?
  float ISlimit = mShiftManager->GetRegularShiftLimit();

  numPoints = (int)ptsZ.size();
  if (numPoints < MIN_NUM_POINTS_TO_FIT_PLANE) {
    mess.Format("Only %d point(s) added but %d are required to fit to plane",
      numPoints, MIN_NUM_POINTS_TO_FIT_PLANE);
    return -1;
  }

  cenX = (VECTOR_MIN(ptsX) + VECTOR_MAX(ptsX)) / 2.f;
  cenY = (VECTOR_MIN(ptsY) + VECTOR_MAX(ptsY)) / 2.f;

  distToCen = 1.e10;
  for (int i = 0; i < numPoints; i++) {
    indexVec.push_back(i);
    dist = (cenX - ptsX[i]) * (cenX - ptsX[i]) + (cenY - ptsY[i]) * (cenY - ptsY[i]);
    if (i == 0 || dist < distToCen) {
      centerPtInd = i;
      distToCen = dist;
    }
  }

  //Place center point first in final vectors
  sortedIndexVec.resize(1, indexVec[centerPtInd]);
  finalPtsX.resize(1, ptsX[centerPtInd]);
  finalPtsY.resize(1, ptsY[centerPtInd]);

  // Remove points beyond image shift limit
  for (int i = 0; i < numPoints; i++) {
    if (i != centerPtInd) {
      if (fabs(ptsX[i] - ptsX[centerPtInd]) <= ISlimit &&
        fabs(ptsY[i] - ptsY[centerPtInd]) <= ISlimit) {
        sortedIndexVec.push_back(indexVec[i]);
        finalPtsX.push_back(ptsX[i]);
        finalPtsY.push_back(ptsY[i]);
        finalPtsZ.push_back(ptsZ[i]);
      }
    }
  }

  //Check that there are still enough remaining points
  numPoints = (int)sortedIndexVec.size();
  if (numPoints < MIN_NUM_POINTS_TO_FIT_PLANE) {
    mess.Format("Only %d point(s) within IS limits but %d are required to fit to plane",
      numPoints, MIN_NUM_POINTS_TO_FIT_PLANE);
    return -2;
  }

  //Check that points sufficiently span
  ConvexHullLongAxis(finalPtsX, finalPtsY, &ratio, &longAxis);
  if (longAxis / ratio < minSpanOfPlanePts) {
    mess.Format("Given points are not suitable for a least squares fit. "
      "Points should span at least %.1f microns in perpendicular directions",
      minSpanOfPlanePts);
    return 1;
  }

  ptsX = finalPtsX;
  ptsY = finalPtsY;
  ptsZ = finalPtsZ;

  return 0;
}

// Given a set of points, determine if image shifts from the center point to the other 
// points are within IS limits. By default the first item is the starting point, unless 
// sortedIndexVec is given in which case the center-most point will be the starting point.
int CParallelTSHelper::AssessISTargetShiftLimit(IntVec indexVec, IntVec &farInd, 
  IntVec *sortedIndexVec)
{
  int numPoints, newInd;
  FloatVec ptsX, ptsY;
  CMapDrawItem *item;
  MapItemArray *itemArr = mWinApp->mNavigator->GetItemArray();
  float stageZ, cenX, cenY;
  float ISlimit = mShiftManager->GetRegularShiftLimit();
  farInd.clear();

  numPoints = (int)indexVec.size();

  for (int i = 0; i < numPoints; i++) {
    item = itemArr->GetAt(indexVec[i]);
    ptsX.push_back(item->mStageX);
    ptsY.push_back(item->mStageY);
    if (i == 0)
      stageZ = item->mStageZ;
  }

  if (sortedIndexVec) {

    cenX = (VECTOR_MIN(ptsX) + VECTOR_MAX(ptsX)) / 2.f;
    cenY = (VECTOR_MIN(ptsY) + VECTOR_MAX(ptsY)) / 2.f;

    mWinApp->mNavigator->AddItemFromStagePositions(&cenX, &cenY, 1, stageZ,
      mParallelTSDlg->GetFitPlaneGroupID());
    itemArr = mWinApp->mNavigator->GetItemArray();
    newInd = (int)itemArr->GetSize() - 1;
    item = itemArr->GetAt(newInd);
    item->mNote = "[Temp Center Pt]";
    mWinApp->mNavigator->UpdateListString(newInd);
    mWinApp->mNavigator->Redraw();

    //Place center point first in vectors
    sortedIndexVec->resize(1, newInd);
    indexVec.insert(indexVec.begin(), newInd);
    ptsX.insert(ptsX.begin(), cenX);
    ptsY.insert(ptsY.begin(), cenY);
  } else {
    cenX = ptsX[0];
    cenY = ptsY[0];
  }

  // Remove points beyond image shift limit
  for (int i = 1; i < (int)indexVec.size(); i++) {
    if (fabs(ptsX[i] - cenX) <= ISlimit && fabs(ptsY[i] - cenY) <= ISlimit) {
      if (sortedIndexVec)
        sortedIndexVec->push_back(indexVec[i]);
    } else {
      farInd.push_back(indexVec[i]);
    }
  }

  return (int)farInd.size();
}

// Assess if the provided points are sufficient for a 2D least squares fit
int CParallelTSHelper::AssessPtsToFitPlane(IntVec indexVec, IntVec &sortedIndexVec, 
  CString &mess)
{
  int numPoints, err;
  FloatVec ptsX, ptsY;
  IntVec farInds;
  CMapDrawItem *item;
  MapItemArray *itemArr = mWinApp->mNavigator->GetItemArray();
  float ratio, longAxis;

  float minSpanOfPlanePts = 1.f; //TODO get this from parameters or something?
  float ISlimit = mShiftManager->GetRegularShiftLimit();

  numPoints = (int)indexVec.size();
  if (numPoints < MIN_NUM_POINTS_TO_FIT_PLANE) {
    mess.Format("Only %d point(s) added but %d are required to fit to plane",
      numPoints, MIN_NUM_POINTS_TO_FIT_PLANE);
    return -1;
  }

  err = AssessISTargetShiftLimit(indexVec, farInds, &sortedIndexVec);

  //Check that there are still enough remaining points
  numPoints = (int)sortedIndexVec.size();
  if (numPoints < MIN_NUM_POINTS_TO_FIT_PLANE) {
    mess.Format("Only %d point(s) within IS limits but %d are required to fit to plane",
      numPoints, MIN_NUM_POINTS_TO_FIT_PLANE);
    return -2;
  }

  for (int i = 0; i < (int)sortedIndexVec.size(); i++) {
    item = itemArr->GetAt(sortedIndexVec[i]);
    ptsX.push_back(item->mStageX);
    ptsY.push_back(item->mStageY);
  }

  //Check that points sufficiently span
  //These should be in specimen coords I think
  ConvexHullLongAxis(ptsX, ptsY, &ratio, &longAxis);
  if (longAxis / ratio < minSpanOfPlanePts) {
    mess.Format("Given points are not suitable for a least squares fit. "
      "Points should span at least %.1f microns in perpendicular directions",
      minSpanOfPlanePts);
    return 1;
  }

  return 0;
}

// When a new round of targets are added, this appends new targets to existing targets
int CParallelTSHelper::AppendNewTargets(IntVec targetMapIDs, CString &mess)
{
  int ind, jnd, size;
  bool append;

  size = (int)mISTargetPointIDs.size();

  if (size == 0) {
    mStartIndex = 0;
    mISTargetPointIDs = targetMapIDs;
    mAlignedToFirstISTarget = false;
    return 0;
  } 
  mStartIndex = size > 0 ? size - 1 : 0;

  for (ind = 0; ind < (int)targetMapIDs.size(); ind++) {
    append = true;
    for (jnd = 0; jnd < size; jnd++) {
      if (mISTargetPointIDs[jnd] == targetMapIDs[ind]) {
        append = false;
        break;
      }
    }
    if (append) {
      mISTargetPointIDs.push_back(targetMapIDs[ind]);
    }
  }

  return 0;
}

// Searches the navigator for items with IDs that have been saved during routine
int CParallelTSHelper::GetSavedTargetsInNav(IntVec *navInd, IntVec *indices)
{
  int ind, jnd, numPoints, numSaved;
  CMapDrawItem *item;
  MapItemArray *itemArr = mWinApp->mNavigator->GetItemArray();
  if (indices)
    indices->clear();
  navInd->clear();
  numSaved = (int)mSavedTargetIDs.size();
  for (ind = 0; ind < numSaved; ind++) {
    if (mWinApp->mNavigator->FindItemWithMapID(mSavedTargetIDs[ind], false)) {
      navInd->push_back(mWinApp->mNavigator->GetFoundItem());
    }
  }
  numPoints = (int)navInd->size();
  
  //Sort saved targets by their order in the navigator
  std::sort(navInd->begin(), navInd->end());

  if (indices) {
    for (ind = 0; ind < numPoints; ind++) {
      item = itemArr->GetAt(navInd->at(ind));
      for (jnd = 0; jnd < numSaved; jnd++) {
        if (item->mMapID == mSavedTargetIDs[jnd]) {
          indices->push_back(jnd);
          break;
        }
      }
    }
  }

  return numPoints;
}

// Converts saved IS targets or the given item to a parallel TS item. 
int CParallelTSHelper::ConvertToParTSItem(CString &err, CMapDrawItem *item)
{
  int ind, jnd;
  CMapDrawItem *mapItem;
  MapItemArray *itemArr = mWinApp->mNavigator->GetItemArray();
  IntVec navInd, indices;
  FloatVec delISX, delISY;
  int numPoints, numX, numY, numDef;
  float ptX, ptY;
  EMimageBuffer *imBuf;
  double ISXcen, ISYcen;

  mapItem = mWinApp->mNavigator->FindItemWithMapID(mAreaMapID);
  if (!mapItem) {
    err.Format("The Parallel Tilt Series area map was deleted");
    return -1;
  }

  mParTSParam.prevSectNums.clear();
  mParTSParam.xCoordInArea.clear();
  mParTSParam.yCoordInArea.clear();
  mParTSParam.xShiftInImage.clear();
  mParTSParam.yShiftInImage.clear();

  //Handle a regular pattern set up from multishot params
  if (item) {
    MultiShotParams *msPars = mNavHelper->GetMultiShotParams();
    if (item->IsPoint()) {
      mNavHelper->GetNumHolesFromParam(numX, numY, numDef);
      numPoints = mNavHelper->GetNumHolesForItem(item, numDef);
      if (item->mNumXholes > 0)
        numX = item->mNumXholes;
      if (item->mNumYholes > 0)
        numY = item->mNumYholes;
      mWinApp->mParticleTasks->GetHolePositions(delISX, delISY, indices, mMagIndex,
        mWinApp->GetCurrentCamera(), numX, numY, mMappingTilt, item);
      mParTSitem = item;
    } else if (item->IsPolygon()) {
      mNavHelper->FillPolygonWithMultiShot(item, delISX, delISY, err);
      numPoints = (int) delISX.size();
      ptX = item->mStageX;
      ptY = item->mStageY;
      mWinApp->mNavigator->AddItemFromStagePositions(&ptX, &ptY, 1, item->mStageZ, 0);
      mParTSitem = itemArr->GetAt(itemArr->GetSize() - 1);
    } else {
      err.Format("The item must be a point or polygon to finalize the area");
      return 1;
    }
    
    mParTSitem->mNumIStargets = (short)numPoints;
    mParTSitem->mIStargetsXY = new float[numPoints * 2];
    mParTSitem->mMagOfIStargets = msPars->holeMagIndex[msPars->doHexArray ? 1 : 0];

    for (ind = 0; ind < numPoints; ind++) {
      mParTSitem->mIStargetsXY[2 * ind] = delISX[ind];
      mParTSitem->mIStargetsXY[2 * ind + 1] = delISY[ind];
    }
    
  } else {

    // If skipping refinements and no maps saved, need to find image shifts for nav points
    if ((mParTSopts->flags & PTSFLAG_SKIP_REFINE) && mParTSopts->extractVirtPrevs != 0) {
      if (GetISVectors(mParallelTSDlg->GetTargetGroupID(), err)) {
        return 2;
      }
    } 

    //Handle custom targets places in arbitrary positions
    numPoints = GetSavedTargetsInNav(&navInd, &indices);
    if (numPoints < 2) {
      err.Format("At least two targets are required to create a parallel tilt series "
        "item");
      return 3;
    }

    //Center point information
    mParTSitem = itemArr->GetAt(navInd[0]);
    mParTSitem->mNumIStargets = (short)numPoints;
    mParTSitem->mIStargetsXY = new float[numPoints * 2];
    mParTSitem->mMagOfIStargets = (mParTSopts->flags & PTSFLAG_SKIP_REFINE) ? 
      mAreaMapMagInd : mMagIndex;

    // Get IS and Preview map for first item
    if (mParTSopts->extractVirtPrevs == 0) {
      item = mWinApp->mNavigator->FindItemWithMapID(mPreviewMapIDs[indices[0]]);
      if (item) {
        mParTSParam.firstPrevMapID = mPreviewMapIDs[indices[0]];
      } else {
        err.Format("The first target map item no longer exists. ");
        return 3;
      }
    } else {
      mParTSParam.firstPrevMapID = 0;
    }
    ISXcen = mISTargetISX[indices[0]];
    ISYcen = mISTargetISY[indices[0]];

    // If using extracts, load the area map
    if (mParTSopts->extractVirtPrevs == 1) {
      mWinApp->mNavigator->DoLoadMap(true, mapItem, -1);
      imBuf = &mWinApp->GetImBufs()[mWinApp->mBufferManager->GetBufToReadInto()];
    }

    for (ind = 0; ind < numPoints; ind++) {
      jnd = indices[ind];
      mParTSitem->mIStargetsXY[2 * ind] = (float)(mISTargetISX[jnd] - ISXcen);
      mParTSitem->mIStargetsXY[2 * ind + 1] = (float)(mISTargetISY[jnd] - ISYcen);
      
      if (mParTSopts->extractVirtPrevs == 0) {
        mParTSParam.prevSectNums.push_back(mPrevMapSectNums[jnd]);
        mParTSParam.xShiftInImage.push_back(mPrevMapShiftX[jnd]);
        mParTSParam.yShiftInImage.push_back(mPrevMapShiftY[jnd]);
      } else if (mParTSopts->extractVirtPrevs == 1) {
        item = itemArr->GetAt(navInd[ind]);
        mWinApp->mMainView->GetItemImageCoords(imBuf, item, ptX, ptY);
        mParTSParam.xCoordInArea.push_back(ptX);
        mParTSParam.yCoordInArea.push_back(ptY);
      }
    }

    //Delete the extra nav items now that all parameters are stored in one item
    DeleteTargetsFromNav(true);
    ClearTargets(false);

    //Close preview map file if open
    ind = mWinApp->mDocWnd->StoreIndexFromName(mTargetMapFileName);
    if (ind >= 0) {
      mWinApp->mDocWnd->SetCurrentStore(ind);
      mWinApp->mDocWnd->DoCloseFile();
    }

    item = mParTSitem;
  }

  mParTSParam.navID = mWinApp->mNavigator->MakeUniqueID();
  mParTSParam.preTilt = roundf(mPretilt * 100) / 100.f;
  mParTSParam.xPitchAngle = roundf(mXpitch * 100) / 100.f;
  mParTSParam.mappingTilt = roundf(mMappingTilt * 100) / 100.f;
  mParTSParam.maxTiltFromStart = mNavHelper->mParallelTSDlg->m_fMaxTilt;

  ParallelTSParam *parTSParam = new ParallelTSParam;
  *parTSParam = mParTSParam;
  mParTSitem->mParallelTSIndex = 
    (int)mWinApp->mNavigator->GetParallelTSArray()->Add(parTSParam);

  mWinApp->mNavigator->FindItemWithMapID(mParTSitem->mMapID, false);
  ind = mWinApp->mNavigator->GetFoundItem();
  mParTSitem->mNote.Format("%d TS targets", mParTSitem->mNumIStargets);
  mWinApp->mNavigator->UpdateListString(ind);
  mWinApp->mNavigator->SetSelectedItem(ind);

  mParallelTSDlg->SetInstructionLine("Select the desired tilt series file.");

  //Check the TS box. If already checked for some reason, uncheck it first.
  if (mWinApp->mNavigator->m_bTiltSeries) {
    mWinApp->mNavigator->m_bTiltSeries = false;
    mWinApp->mNavigator->UpdateData(false);
    mWinApp->mNavigator->OnCheckTiltSeries();
  }
  mWinApp->mNavigator->m_bTiltSeries = true;
  mWinApp->mNavigator->UpdateData(false);
  mWinApp->mNavigator->OnCheckTiltSeries();

  mSavedTSparamIndex = mParTSitem->mTSparamIndex;
  return 0;
}

// Gets the item from which TS parameters are being inherited, and returnss the nav index 
int CParallelTSHelper::GetTSparamItem(CMapDrawItem *&item)
{
  item = NULL;
  MapItemArray *itemArr = mWinApp->mNavigator->GetItemArray();
  int ind, navInd = -1;

  if (mSavedTSparamIndex >= 0) {
    for (ind = 0; ind < itemArr->GetSize(); ind++) {
      item = itemArr->GetAt(ind);
      if (item->mTSparamIndex == mSavedTSparamIndex) {
        navInd = ind;
        break;
      }
    }
  }
  return navInd;
}

// Opens the tilt series parameters window for the current PTS item
void CParallelTSHelper::UpdateTSParams()
{
  CMapDrawItem *item;
  int ind;
  
  if (!mParTSitem)
    return;
  item = mWinApp->mNavigator->FindItemWithMapID(mParTSitem->mMapID, false);
  ind = mWinApp->mNavigator->GetFoundItem();
  if (ind < 0) {
    return;
  }
  mWinApp->mNavigator->SetSelectedItem(ind, false);
  mWinApp->mNavigator->OnButTsparams();
  mSavedTSparamIndex = item->mTSparamIndex;
}

// Checks preconditions and applies the adjusting transform to a given set of multishot
// vectors
bool CParallelTSHelper::CanAdjustISVectors(int fromMag, bool multiShot, CString &mess)
{
  int fromMagCam = 0, toMagCam = 0;
  int camera = mWinApp->GetCurrentCamera();
  if (multiShot) {
    MultiShotParams *params = mNavHelper->GetMultiShotParams();
    return mNavHelper->AdjustMultiShotVectors(params,
      B3DCHOICE(params->useCustomHoles, 1, params->doHexArray ? -1 : 0), true, mess) == 0;
  } else {
    AdjustXformData *adjustData = mNavHelper->GetNearestAdjustingXform(fromMag);
    if (adjustData) {
      fromMagCam = MagForCamera(camera, adjustData->xformFromMag);
      toMagCam = MagForCamera(camera, adjustData->xformToMag);
    }
    if (!adjustData || fromMagCam == 0 || toMagCam == 0) {
      mess = "No adjusting transform available";
      return false;
    }
    mess.Format("Adjusting transform available from %dx to %dx", fromMagCam, toMagCam);
    return adjustData->xformFromMag == fromMag;
  }
}

// Get the ID of the center point in the target group, which is always the first target
int CParallelTSHelper::GetCenterPtID()
{
  if (mSavedTargetIDs.size() > 0)
    return mSavedTargetIDs[0];
  else
    return -1;
}

// Just compute the image shift vectors without any refinement
int CParallelTSHelper::GetISVectors(int groupID, CString &err)
{
  int ind, numPoints, numAcq;
  CMapDrawItem *item;
  MapItemArray *itemArr = mWinApp->mNavigator->GetItemArray();
  float cenStageX, cenStageY;
  double delX, delY;
  double ISX, ISY;
  CString label, mess;
  IntVec indexVec;
  ScaleMat st2is;
  bool canAdjustIS = mShiftManager->GetFocusISCals()->GetSize() > 0 &&
    mShiftManager->GetFocusMagCals()->GetSize() > 0;

  numPoints = mWinApp->mNavigator->CountItemsInGroup(groupID, label, mess, numAcq, 
    &indexVec);
  if (numPoints < 2) {
    err.Format("At least two targets are required to create a parallel tilt series item");
    return -1;
  }

  CMapDrawItem *mapItem = mWinApp->mNavigator->FindItemWithMapID(mAreaMapID);
  if (!mapItem) {
    err.Format("The map on which items were drawn no longer exists");
    return -2;
  }

  st2is = MatMul(mShiftManager->StageToCamera(mWinApp->GetCurrentCamera(), mAreaMapMagInd), 
    mShiftManager->CameraToIS(mAreaMapMagInd));
  if (!st2is.xpx) {
    err.Format("There is no calibration to get from stage to IS coordinates at the "
    "given mag");
    return 1;
  }

  for (ind = 0; ind < numPoints; ind++) {
    item = itemArr->GetAt(indexVec[ind]);
    if (ind == 0) {
      cenStageX = item->mStageX;
      cenStageY = item->mStageY;
      ISX = ISY = 0.f;
    } else {

      // Get IS from center point to the next target
      mShiftManager->ApplyScaleMatrix(st2is, item->mStageX - cenStageX,
        item->mStageY - cenStageY, ISX, ISY);
      mShiftManager->TransferGeneralIS(mAreaMapMagInd, ISX, ISY, mMagIndex, delX, delY);

      if (mParTSopts->applyAdjustingXform) {
        if (CanAdjustISVectors(mAreaMapMagInd, false, mess)) {
          AdjustXformData *adjustData = 
            mNavHelper->GetNearestAdjustingXform(mAreaMapMagInd);
          SEMTrace('N', "adj xform xpx=%.4f, xpy=%.4f, ypx=%.4f, ypy=%.4f",
            adjustData->adjustingXform.xpx, adjustData->adjustingXform.xpy,
            adjustData->adjustingXform.ypx, adjustData->adjustingXform.ypy);
          ApplyScaleMatrix(adjustData->adjustingXform, delX, delY, ISX, ISY);
          SEMTrace('N', "IS transformed from (%.4f, %.4f) to (%.4f, %.4f)",
            delX, delY, ISX, ISY);
        } else {
          SEMTrace('N', "Adjusting transform not applied: %s", mess);
        }
      }

      if (!mShiftManager->ImageShiftIsOK(ISX, ISY, FALSE)) {
        SEMTrace('N', "Target at navigator index %d is beyond image shift limit and will "
          "be skipped", indexVec[ind]);
        continue;
      }
    }

    mISTargetISX.push_back(ISX);
    mISTargetISY.push_back(ISY);
    mSavedTargetIDs.push_back(item->mMapID);
  }

  return 0;
}

// Use the saved image shifts and adjustments to generate an adjusting transform
int CParallelTSHelper::ComputeAdjustingTransform(CString &err)
{
  int ind, jnd;
  IntVec navInd, indices;
  int numPoints;
  double ISXcen, ISYcen, delX, delY;
  double maxErr = 0., meanErr = 0.;
  CString label;
  std::vector<double> fromISX, fromISY, toISX, toISY, predErr;

  //Handle custom targets places in arbitrary positions
  numPoints = GetSavedTargetsInNav(&navInd, &indices);
  if (numPoints < MIN_NUM_POINTS_FOR_PTSADJUST) {
    err.Format("There must be at least %d points to create an adjusting transform.", 
      MIN_NUM_POINTS_FOR_PTSADJUST);
    return 1;
  }

  ISXcen = mISTargetISX[indices[0]];
  ISYcen = mISTargetISY[indices[0]];

  for (ind = 1; ind < numPoints; ind++) {
    jnd = indices[ind];
    fromISX.push_back(mPreRefineISX[jnd] - ISXcen);
    fromISY.push_back(mPreRefineISY[jnd] - ISYcen);
    toISX.push_back(mISTargetISX[jnd] - ISXcen);
    toISY.push_back(mISTargetISY[jnd] - ISYcen);
    delX = mISTargetISX[jnd] - mPreRefineISX[jnd];
    delY = mISTargetISY[jnd] - mPreRefineISY[jnd];
  }

  mAdjustingXform = mShiftManager->GetTransformFromISAdjustments(fromISX, fromISY, toISX, 
    toISY, predErr);

  if (mAdjustingXform.xpx == 0) {
    err.Format("A transform could not be determined from the given IS adjustments");
    return 2;
  }
  if (mOldAdjustingXform.xpx != 0) {
    mAdjustingXform = MatMul(mAdjustingXform, mOldAdjustingXform);
  }
  
  jnd = 0;
  for (ind = 0; ind < (int)predErr.size(); ind++) {
    meanErr += predErr[ind];
    SEMTrace('N', "Adj transform fit error = %.1f nm at navigator item #%d (%s)", 
      predErr[ind] * 1.e3, navInd[ind + 1] + 1,
      mWinApp->mNavigator->GetItemArray()->GetAt(navInd[ind + 1])->mLabel);
    if (predErr[ind] > maxErr) {
      maxErr = predErr[ind];

      // Shift by one because first item is center point for which adjustment is not done
      jnd = navInd[ind + 1] + 1;
      label = mWinApp->mNavigator->GetItemArray()->GetAt(navInd[ind + 1])->mLabel;
    }
  }
  meanErr /= numPoints;
  PrintfToLog("Adjusting Transform:");
  PrintfToLog("  Mean fit error = %.1f nm", meanErr * 1.e3);
  PrintfToLog("  Max fit error = %.1f nm at navigator item #%d (%s)", 
     maxErr * 1.e3, jnd, label);
  PrintfToLog("  xpx=%.4f, xpy=%.4f, ypx=%.4f, ypy=%.4f",
    mAdjustingXform.xpx, mAdjustingXform.xpy,
    mAdjustingXform.ypx, mAdjustingXform.ypy);
  return 0;
}

void CParallelTSHelper::SaveAdjustingTransform()
{
  AdjustXformData xformData;
  xformData.adjustingXform = mAdjustingXform;
  xformData.xformFromMag = mAreaMapMagInd;
  xformData.xformToMag = mMagIndex;
  xformData.xformMinuteTime = mWinApp->MinuteTimeStamp();
  mNavHelper->AddAdjustingXform(xformData);
}

// Deletes all saved parallelTS targets from the navigator
void CParallelTSHelper::DeleteTargetsFromNav(bool keepFirst)
{
  int numPoints, ind, jnd;
  IntVec navInd, indices;
  CMapDrawItem *item;
  MapItemArray *itemArr = mWinApp->mNavigator->GetItemArray();

  numPoints = GetSavedTargetsInNav(&navInd, &indices);
  mWinApp->mNavigator->m_bCollapseGroups = false;
  for (ind = numPoints - 1; ind >= 0; ind--) {
    if (keepFirst && ind == 0)
      break;
    item = mWinApp->mNavigator->FindItemWithMapID(mSavedTargetIDs[ind], false);
    jnd = mWinApp->mNavigator->GetFoundItem();
    if (item) {
      mWinApp->mNavigator->ExternalDeleteItem(item, jnd);
    }
  }
}

void CParallelTSHelper::DeleteTargetMapsFromNav() {
  CMapDrawItem *item;
  int ind, jnd;
  for (ind = 0; ind < (int)mPreviewMapIDs.size(); ind++) {
    item = mWinApp->mNavigator->FindItemWithMapID(mPreviewMapIDs[ind]);
    if (item) {
      jnd = mWinApp->mNavigator->GetFoundItem();
      mWinApp->mNavigator->ExternalDeleteItem(item, jnd);
    }
  }
}


void CParallelTSHelper::UpdateSpecAngles(float pretilt, float xPitch)
{
  CArray<ParallelTSParam*, ParallelTSParam*> *parsArr = 
    mWinApp->mNavigator->GetParallelTSArray();
  ParallelTSParam *pars;

  mPretilt = pretilt;
  mXpitch  = xPitch;
  
  if (mParTSitem && mParTSitem->mParallelTSIndex >= 0 && 
    mParTSitem->mParallelTSIndex < parsArr->GetSize()) {
    pars = parsArr->GetAt(mParTSitem->mParallelTSIndex);
    pars->preTilt = roundf(mPretilt * 100) / 100.f;
    pars->xPitchAngle = roundf(mXpitch * 100) / 100.f;
  }
}

void CParallelTSHelper::UpdateMaxTilt(float angle)
{
  CArray<ParallelTSParam*, ParallelTSParam*> *parsArr =
    mWinApp->mNavigator->GetParallelTSArray();
  ParallelTSParam *pars;

  if (mParTSitem && mParTSitem->mParallelTSIndex >= 0 &&
    mParTSitem->mParallelTSIndex < parsArr->GetSize()) {
    pars = parsArr->GetAt(mParTSitem->mParallelTSIndex);
    pars->maxTiltFromStart = angle;
  }
}

int CParallelTSHelper::PruneDeletedTargets()
{
  CMapDrawItem *item;
  int numDeleted = 0;
  if (mActionAtTarget == PARALLELTS_ACTION_PREVIEW && 
    (mParTSopts->flags & PTSFLAG_SKIP_REFINE) && mParTSopts->extractVirtPrevs != 0)
    return 0;

  for (int ind = 0; ind < (int)mSavedTargetIDs.size(); ind++) {
    item = mWinApp->mNavigator->FindItemWithMapID(mSavedTargetIDs[ind], false);
    if (!item) {
      numDeleted++;
      VEC_REMOVE_AT(mSavedTargetIDs, ind);
      if (mActionAtTarget == PARALLELTS_ACTION_ADJUST) {
        VEC_REMOVE_AT(mISTargetISX, ind);
        VEC_REMOVE_AT(mISTargetISY, ind);
        VEC_REMOVE_AT(mPreRefineISX, ind);
        VEC_REMOVE_AT(mPreRefineISY, ind);
      } else if (mActionAtTarget == PARALLELTS_ACTION_AUTOFOCUS) {
        VEC_REMOVE_AT(mISTargetSSX, ind);
        VEC_REMOVE_AT(mISTargetSSY, ind);
        VEC_REMOVE_AT(mISTargetDefocus, ind);
      } else if (mActionAtTarget == PARALLELTS_ACTION_PREVIEW) {
        VEC_REMOVE_AT(mISTargetISX, ind);
        VEC_REMOVE_AT(mISTargetISY, ind);
        VEC_REMOVE_AT(mPreRefineISX, ind);
        VEC_REMOVE_AT(mPreRefineISY, ind);
        if (mParTSopts->extractVirtPrevs == 0 && ind < (int)mPreviewMapIDs.size() ) {
          item = mWinApp->mNavigator->FindItemWithMapID(mPreviewMapIDs[ind], false);
          if (item) {
            mWinApp->mNavigator->ExternalDeleteItem(item,
              mWinApp->mNavigator->GetFoundItem());
          }
          VEC_REMOVE_AT(mPreviewMapIDs, ind);
          VEC_REMOVE_AT(mPrevMapSectNums, ind);
          VEC_REMOVE_AT(mPrevMapShiftX, ind);
          VEC_REMOVE_AT(mPrevMapShiftY, ind);
        }
      }
    }
  }
  return numDeleted;
}

/*
* Attempts to find outliers among the shifts between expected and actual position for
* the various targets; derived from findXfWithoutOutliers from tiltxcorr.cpp in IMOD,
* it tries dropping successively more points to find a consistent point where dropped
* points are statistically out of range compared to the rest
*/
void CParallelTSHelper::IdentifyParTSOutliers(float *xload, float *yload, float *xfound,
  float *yfound, int *dropping, int numFit, int maxDrop, float elimMin, float critProb,
  float absProbCrit, float &slope, float *intcp)
{
  FloatVec devs, xftmp, yftmp, xltmp, yltmp;
  IntVec devInds, topInds, tmpInds;
  float devAvg, devSd, maxErr;
  float dev, probPerPoint, absPerPoint, z, prob, gprob;
  float sigmaFromMean, sigmaFromSD, sigma;
  int indMax, from;
  int totalDrop = 0, numDrop = 0, ind, jdrop, numKeep;
  int lastDrop = 0;
  CString mess, str;

  devs.resize(numFit);
  if (intcp) {
    FitZvsYWithDropping(xload, yload, dropping, numFit, slope, *intcp, &devs[0],
      devAvg, devSd, maxErr, indMax);
    SEMTrace('N', "Initial fit Y vs Z avg slope %f  intcp %.3f dev %.3f  max %.3f at %d",
      slope, *intcp, devAvg, maxErr, indMax);
  } else {
    FitYSlopeWithDropping(xload, yload, xfound, yfound, dropping, numFit, slope, &devs[0],
      devAvg, devSd, maxErr, indMax);
    SEMTrace('N', "Initial fit exp to actual avg dev %.3f  max %.3f at %d", devAvg, maxErr,
      indMax);
  }
  if (maxErr < elimMin || numFit <= (intcp ? 3 : 2))
    return;

  // Sort the residuals from the full fit
  for (ind = 0; ind < numFit; ind++)
    devInds.push_back(ind);
  rsSortIndexedFloats(&devs[0], &devInds[0], numFit);
  probPerPoint = (float)pow(1. - critProb, 1. / numFit);
  absPerPoint = (float)pow(1. - absProbCrit, 1. / numFit);

  // Reload the data in that order
  for (ind = 0; ind < numFit; ind++) {
    from = devInds[ind];
    xltmp.push_back(xload[from]);
    yltmp.push_back(yload[from]);
    xftmp.push_back(xfound[from]);
    if (!intcp)
      yftmp.push_back(yfound[from]);
  }
  /*
  FitYSlopeWithDropping(&xltmp[0], &yltmp[0], &xftmp[0], &yftmp[0], dropping, numFit,
  slope, &devs[0], devAvg, devSd, maxErr, indMax);
  mess = "Refit with all ordered data";
  for (ind = 0; ind < numFit; ind++) {
  str.Format(" %.4f", devs[ind]);
  mess += str;
  }
  SEMTrace('1', "%s", (LPCTSTR)mess);*/

  for (jdrop = 1; jdrop <= maxDrop + 1; jdrop++) {
    for (ind = 0; ind < numFit; ind++)
      dropping[ind] = (ind < numFit - jdrop) ? 0 : 1;
    if (intcp)
      FitZvsYWithDropping(&xltmp[0], &yltmp[0], dropping, numFit,
        slope, *intcp, &devs[0], devAvg, devSd, maxErr, indMax);
    else
      FitYSlopeWithDropping(&xltmp[0], &yltmp[0], &xftmp[0], &yftmp[0], dropping, numFit,
        slope, &devs[0], devAvg, devSd, maxErr, indMax);
    //mess.Format("drop %d avg dev %.4f  LO devs: ", jdrop, devAvg);
    sigmaFromMean = devAvg / sqrt(8.f / 3.14159f);
    sigmaFromSD = devSd / sqrtf(3.f - 8.f / 3.14159f);
    sigma = B3DMAX(sigmaFromMean, sigmaFromSD);
    numKeep = 0;

    // load the devs of the dropped ones and sort them
    topInds.resize(jdrop);
    for (ind = numFit - jdrop; ind < numFit; ind++) {
      devs[ind - (numFit - jdrop)] = devs[ind];
      topInds[ind - (numFit - jdrop)] = ind - (numFit - jdrop);
    }
    if (jdrop > 1)
      rsSortIndexedFloats(&devs[0], &topInds[0], jdrop);

    // Evaluate the deviation of each dropped one
    for (ind = numFit - jdrop; ind < numFit; ind++) {
      dev = devs[topInds[ind - (numFit - jdrop)] + (numFit - jdrop)];
      //str.Format("%.4f ", dev);
      //mess += str;
      if (sigma < 0.1 * dev || sigma < 1.e-5)
        z = 10.;
      else
        z = dev / sigma;
      gprob = (float)(1. - 0.5 * (1. - errFunc(z / 1.414214)));

      prob = 2.f * (gprob - 0.5f) - sqrt(2.f / 3.14159f) * z * expf(-z * z / 2.f);
      if (prob < probPerPoint) {
        //PRINT4(dev, z, prob, "Inc nk");
        numKeep = numKeep + 1;
      }
      if (prob >= absPerPoint) {
        //PRINT4(dev, z, prob, "Set ND");
        numDrop = B3DMIN(maxDrop, B3DMAX(numDrop, numFit - ind));
      }
    }
    //SEMTrace('1', "%s", mess);

    /*
    * If all points are outliers, this is a candidate for a set to drop
    * When only the first point is kept, and all the rest of the points
    * were outliers on the previous round, then this is a safe place to
    * draw the line between good data and outliers.  In this case, set
    * numDrop; and at end take the biggest ndrop that fits these criteria
    */
    if (numKeep == 0)
      lastDrop = jdrop;
    if (numKeep == 1 && lastDrop == jdrop - 1 && lastDrop > 0)
      numDrop = lastDrop;

    // Maintain the indexes in the current order
    tmpInds.resize(jdrop);
    for (ind = 0; ind < jdrop; ind++)
      tmpInds[ind] = devInds[topInds[ind] + (numFit - jdrop)];
    for (ind = numFit - jdrop; ind < numFit; ind++)
      devInds[ind] = tmpInds[ind - (numFit - jdrop)];

    // Reload the dropped ones in the right order
    for (ind = numFit - jdrop; ind < numFit; ind++) {
      from = devInds[ind];
      xltmp[ind] = xload[from];
      yltmp[ind] = yload[from];
      xftmp[ind] = xfound[from];
      if (!intcp)
        yftmp[ind] = yfound[from];
    }
  }

  // Redo fit with final numDrop
  for (ind = 0; ind < numFit; ind++)
    dropping[ind] = (ind < numFit - numDrop) ? 0 : 1;
  if (intcp)
    FitZvsYWithDropping(&xltmp[0], &yltmp[0], dropping, numFit,
      slope, *intcp, &devs[0], devAvg, devSd, maxErr, indMax);
  else
    FitYSlopeWithDropping(&xltmp[0], &yltmp[0], &xftmp[0], &yftmp[0], dropping, numFit,
      slope, &devs[0], devAvg, devSd, maxErr, indMax);

  // Mark the selected number as dropped
  mess = "";
  for (ind = 0; ind < numFit; ind++) {
    dropping[devInds[ind]] = (ind < numFit - numDrop) ? 0 : 1;
    if (ind >= numFit - numDrop) {
      str.Format(" %d", devInds[ind] + 1);
      mess += str;
    }
  }
  if (!mess.IsEmpty())
    SEMTrace('N', "Final drop %d: pos %s", numDrop, (LPCTSTR)mess);
}

/*
* Finds the best least-squared fit for scaling between found and loaded Y values,
* computes error between found X/Y and loaded positions with this scaling for all
* points, and returns the mean and sd error for non-dropped points
*/
void CParallelTSHelper::FitYSlopeWithDropping(float *xload, float *yload, float *xfound,
  float *yfound, int *dropping, int numPos, float &slope, float *errors, float &mean,
  float &sd, float &maxErr, int &indMax)
{
  int ind, numFit = 0;
  float xsqsum = 0., xysum = 0., delx, dely;
  FloatVec errArr;
  for (ind = 0; ind < numPos; ind++) {
    if (!dropping[ind]) {
      xsqsum += yload[ind] * yload[ind];
      xysum += yfound[ind] * yload[ind];
      numFit++;
    }
  }
  slope = xysum / xsqsum;
  maxErr = -1.;
  for (ind = 0; ind < numPos; ind++) {
    delx = xfound[ind] - xload[ind];
    dely = yfound[ind] - slope * yload[ind];
    errors[ind] = sqrtf(delx * delx + dely * dely);
    if (!dropping[ind]) {
      errArr.push_back(errors[ind]);
      if (errors[ind] > maxErr) {
        maxErr = errors[ind];
        indMax = ind;
      }
    }
  }
  avgSD(&errArr[0], numFit, &mean, &sd, &delx);
}

/*
* Does simple linear fit of Z versus Y, dropping points marged in the dropping array
*/
void CParallelTSHelper::FitZvsYWithDropping(float *yload, float *zvec, int *dropping,
  int numPos, float &slope, float &intcp, float *errors, float &mean, float &sd,
  float &maxErr, int &indMax)
{
  FloatVec ytemp, ztemp, errArr;
  int ind;
  float ro;
  for (ind = 0; ind < numPos; ind++) {
    if (!dropping[ind]) {
      ytemp.push_back(yload[ind]);
      ztemp.push_back(zvec[ind]);
    }
  }
  lsFit(&ytemp[0], &ztemp[0], (int)ytemp.size(), &slope, &intcp, &ro);
  for (ind = 0; ind < numPos; ind++) {
    errors[ind] = B3DABS(yload[ind] * slope + intcp - zvec[ind]);
    if (!dropping[ind]) {
      errArr.push_back(errors[ind]);
      if (errors[ind] > maxErr) {
        maxErr = errors[ind];
        indMax = ind;
      }
    }
  }
  avgSD(&errArr[0], (int)ytemp.size(), &mean, &sd, &ro);
}

/*
* Fit to data on Y positions only at one target point to derive yZero and zZero values
* and return summed and mean errors at teh specified pretilt
*/
void CParallelTSHelper::FitParallelTSYonly(FloatVec &angles, FloatVec &yOnly, float preTilt,
  float &yZero, float &zZero, double &errSum, float &meanErr)
{
  FloatVec x1, x2, yy;
  int ind, numFit = (int)angles.size();
  for (ind = 0; ind < numFit; ind++) {
    x1.push_back(cosf(DTORFL * (angles[ind] - preTilt)));
    x2.push_back(-sinf(DTORFL * (angles[ind] - preTilt)));
    yy.push_back(yOnly[ind]);
  }
  lsFit2(&x1[0], &x2[0], &yy[0], numFit, &yZero, &zZero, NULL);
  errSum = 0.;
  for (ind = 0; ind < numFit; ind++)
    errSum += powf(yZero * x1[ind] + zZero * x2[ind] - yy[ind], 2.f);
  meanErr = (float)sqrt(errSum / numFit);
}

/*
* Fit to data on Y and Z at one target point to derive yZero and zZero values and
* return summed and mean errors at teh specified pretilt  (UNUSED)
*/
void CParallelTSHelper::FitParallelTSYwithZ(FloatVec &angles, FloatVec &yWithZ,
  FloatVec &zVec, float preTilt, float &yZero, float &zZero, double &errSum,
  float &meanErr)
{
  FloatVec x1, x2, yy;
  int ind, numFit = (int)angles.size();
  float cosAng, sinAng;
  for (ind = 0; ind < numFit; ind++) {
    cosAng = cosf(DTORFL * (angles[ind] - preTilt));
    sinAng = sinf(DTORFL * (angles[ind] - preTilt));
    x1.push_back(cosAng);
    x2.push_back(-sinAng);
    yy.push_back(yWithZ[ind]);
    x1.push_back(sinAng);
    x2.push_back(cosAng);
    yy.push_back(zVec[ind]);
  }
  lsFit2(&x1[0], &x2[0], &yy[0], numFit * 2, &yZero, &zZero, NULL);
  errSum = 0.;
  for (ind = 0; ind < numFit * 2; ind++)
    errSum += powf(yZero * x1[ind] + zZero * x2[ind] - yy[ind], 2.f);
  meanErr = (float)sqrt(errSum / (2 * numFit));
}

/*
* Identifies outliers in an arbitrary set of values either with a leave-one-out approach
* for < 5 items or using MAD median outlier routine.  Polarity should be 1 or -1 for
* just the outliers above or below the median, or 0 for outliers in either direction.
* If nonzero, elimMin is a threshold value for any elimination; it should be
* a lower limit if polarity is -1, an upper limit if polarity is 1, or values should
* be positive and negative (e.g., differences from a mean) to use elimMin with polarity 0
*/
void CParallelTSHelper::FindOutliersInResultList(FloatVec &errors, float elimMin,
  int polarity, FloatVec &outliers)
{
  int minForRobust = 5;
  float kcrit = 2.5f, semCrit = 2.f;
  float avg, sd, sem, maxSem = -1., numSem;
  int idrop, ind, itmp, maxDrop, numErr = (int)errors.size();
  float temp[8];

  CLEAR_RESIZE(outliers, float, numErr);
  if (elimMin != 0. && ((polarity > 0 && VECTOR_MAX(errors) < elimMin) ||
    (polarity < 0 && VECTOR_MIN(errors) > elimMin) ||
    (!polarity && VECTOR_MAX(errors) < elimMin) && VECTOR_MIN(errors) > -elimMin))
    return;
  if (numErr < minForRobust) {

    // Try dropping each point, getting stats with the rest, and computing # of SEMs from
    // the mean
    for (idrop = 0; idrop < numErr; idrop++) {
      itmp = 0;
      for (ind = 0; ind < numErr; ind++)
        if (ind != idrop)
          temp[itmp++] = errors[ind];
      avgSD(temp, numErr - 1, &avg, &sd, &sem);
      numSem = (errors[idrop] - avg) / sem;
      if ((polarity > 0 && numSem > maxSem) || (polarity < 0 && numSem < maxSem) ||
        (!polarity && fabs(numSem) > fabs(maxSem))) {
        maxSem = numSem;
        maxDrop = idrop;
      }
    }

    // If the worst one exceeds the crit, drop it
    if (fabs(maxSem) > semCrit && (!polarity ||
      BOOL_EQUIV(polarity > 0, maxSem > semCrit)))
      outliers[maxDrop] = 1.;

  } else {

    // Robust is easy if it works
    rsMadMedianOutliers(&errors[0], numErr, kcrit, &outliers[0]);
    if (polarity) {
      for (ind = 0;ind < numErr; ind++)
        if (!BOOL_EQUIV(polarity > 0, outliers[ind] > 0))
          outliers[ind] = 0.;
    }
  }
}