NIREPEvaluator.cxx

Go to the documentation of this file.

#include "NIREPEvaluator.h"
#include "NIREPInverseConsistencyError.h"
#include "NIREPTransitivityError.h"
#include "NIREPIntensityVariance.h"

#include "NIREPDefinitions.h"

/* +----------------+
   | Public Methods |
   +----------------+ */

/*
  function:  NIREPEvaluator
  purpose:   constructor
  overloads: NIREPEvaluator(NIREPDataManager *dm)
*/
NIREPEvaluator::NIREPEvaluator(void)
{
  this->dataManager = NULL;
  //taskList = NULL;
  this->Init();
}


/*
  function:  NIREPEvaluator
  purpose:   Class constructor that initializes the data manager
  overloads: NIREPEvaluator(void) 
*/
NIREPEvaluator::NIREPEvaluator(NIREPDataManager *dm)
{
  this->dataManager = dm;
  this->Init();
}


/* 
  function: ~NIREPEvaluator
  purpose:  class destructor, destroying references to pointers
*/
NIREPEvaluator::~NIREPEvaluator(void){
  //delete dataManager;
  //if(taskList)
  //{
  //  delete taskList;
  //}
  //taskList = NULL;
}


/* +-----------------+
   | Private Methods |
   +-----------------+ */
  
/*
  function: Init
  purpose:  Initializes fList and all other parameters we have control over.
  output:   initialized fList
*/
void NIREPEvaluator::Init()
{
}


/*
  function: Evaluate
  purpose:  evaluates the provided expression
  input:    command std::string containing the abbreviation or expression to
                    be retrieved (ex: _A3 or diff(_A1, _A2))
  output:   gec::SpatialData* containing the computed expression
*/
gec::SpatialData* NIREPEvaluator::Evaluate(const std::string& command, DisplayDescription * displayDescription)
{
  EvaluatorList * taskList = displayDescription->GetEvaluatorList();
  std::string cmd = taskList->GetEvaluatorCommand(command) ; // unwrap to see if it's empty

  gec::SpatialData* res = this->dataManager->GetLoadedSpatialDataListEntry(cmd, displayDescription);

  // if it's empty, it's not in the data manager so we must execute the expression
  if(res == NULL){            
    std::string fcn = taskList->GetTaskFunction(command); 
    std::vector<std::string> params = taskList->GetTaskParameters(command); 

    // Convert fcn to lower case for comparison
    std::transform(fcn.begin(), fcn.end(), fcn.begin(), ::tolower);

    bool CMDprocessed = false;

    if (fcn.compare(EvalCmd["TransformGrid"]) == 0){
      if(params.size() == 1) {     // ensure only four parameters
        //res = this->IntensityVariance(params,"ImageVariance");
        res = this->TransformGrid(params[0], displayDescription);
      } else {
        wxString temp = wxT("ERROR - Wrong number of parameters for command: ");
        temp.Append(wxString::Format (wxT("%s"), cmd.c_str()));
        wxLogError(temp);
        res = NULL;
      }
      CMDprocessed = true; // command was processed
    }

    if (fcn.compare(EvalCmd["VectorField"]) == 0){
      if(params.size() == 1) {     // ensure only four parameters
        //res = this->IntensityVariance(params,"ImageVariance");
        res = this->VectorField(params[0], displayDescription);
      } else {
        wxString temp = wxT("ERROR - Wrong number of parameters for command: ");
        temp.Append(wxString::Format (wxT("%s"), cmd.c_str()));
        wxLogError(temp);
        res = NULL;
      }
      CMDprocessed = true; // command was processed
    }

    if (fcn.compare(EvalCmd["IntensityVarianceImage"]) == 0){
      if(params.size() == 2) {     // ensure only four parameters
        //res = this->IntensityVariance(params,"ImageVariance");
        res = this->IntensityVariance(params,fcn, displayDescription);
      } else {
        wxString temp = wxT("ERROR - Wrong number of parameters for command: ");
        temp.Append(wxString::Format (wxT("%s"), cmd.c_str()));
        wxLogError(temp);
        res = NULL;
      }
      CMDprocessed = true; // command was processed
    }

    if (fcn.compare(EvalCmd["IntensityVarianceTransform"]) == 0){
      if(params.size() == 4){      // ensure only four parameters
        //res = this->IntensityVariance(params,"RegisteredImageVariance");
        res = this->IntensityVariance(params,fcn, displayDescription);
      } else {
        wxString temp = wxT("ERROR - Wrong number of parameters for command: ");
        temp.Append(wxString::Format (wxT("%s"), cmd.c_str()));
        wxLogError(temp);
        res = NULL;
      }
      CMDprocessed = true; // command was processed
    }

    if (fcn.compare(EvalCmd["TransitivityErrorImage"]) == 0){
      if(params.size() == 4){      // ensure only three parameters -- transAB, transBC, transCA, method
        res = this->TransitivityError(params[0],params[1],params[2],params[3], displayDescription);
      } else {
        wxString temp = wxT("ERROR - Wrong number of parameters for command: ");
        temp.Append(wxString::Format (wxT("%s"), cmd.c_str()));
        wxLogError(temp);
        res = NULL;
      }
      CMDprocessed = true; // command was processed
    }
    
    if (fcn.compare(EvalCmd["InverseConsistencyErrorImage"]) == 0){
      if(params.size() == 3){      // ensure only one parameter -- fwdTrans, revTrans, method
        res = this->InverseConsistencyError(params[0],params[1],params[2], displayDescription);
      } else {
        wxString temp = wxT("ERROR - Wrong number of parameters for command: ");
        temp.Append(wxString::Format (wxT("%s"), cmd.c_str()));
        wxLogError(temp);
        res = NULL;
      }
      CMDprocessed = true; // command was processed
    }
    
    if (fcn.compare(EvalCmd["InvertTransformation"]) == 0){
      if(params.size() == 1){      // ensure only one parameter -- transformation
        res = this->InvertTransformation(params[0], displayDescription);
      } else {
        wxString temp = wxT("ERROR - Wrong number of parameters for command: ");
        temp.Append(wxString::Format (wxT("%s"), cmd.c_str()));
        wxLogError(temp);
        res = NULL;
      }
      CMDprocessed = true; // command was processed
    }
    
    if (fcn.compare(EvalCmd["TransformImage"]) == 0){
      if(params.size() == 2){      // ensure only two parameters -- spatial data, transformation
        res = this->Transform(params[0],params[1], displayDescription);
      } else {
        wxString temp = wxT("ERROR - Wrong number of parameters for command: ");
        temp.Append(wxString::Format (wxT("%s"), cmd.c_str()));
        wxLogError(temp);
        res = NULL;
      }
      CMDprocessed = true; // command was processed
    }
    
    if (fcn.compare(EvalCmd["Jacobian"]) == 0){
      if(params.size() == 1){      // ensure only one parameter -- transformation
        res = this->Jacobian(params[0], displayDescription);
      } else {
        wxString temp = wxT("ERROR - Wrong number of parameters for command: ");
        temp.Append(wxString::Format (wxT("%s"), cmd.c_str()));
        wxLogError(temp);
        res = NULL;
      }
      CMDprocessed = true; // command was processed
    }
    
    if (fcn.compare(EvalCmd["Disp2"]) == 0){
      if(params.size() == 1){      // ensure only one parameter -- transformation
        res = this->Disp(params[0],2, displayDescription);
      } else {
        wxString temp = wxT("ERROR - Wrong number of parameters for command: ");
        temp.Append(wxString::Format (wxT("%s"), cmd.c_str()));
        wxLogError(temp);
        res = NULL;
      }
      CMDprocessed = true; // command was processed
    }
    
    if (fcn.compare(EvalCmd["Disp1"]) == 0){
      if(params.size() == 1){      // ensure only one parameter -- transformation
        res = this->Disp(params[0],1, displayDescription);
      } else {
        wxString temp = wxT("ERROR - Wrong number of parameters for command: ");
        temp.Append(wxString::Format (wxT("%s"), cmd.c_str()));
        wxLogError(temp);
        res = NULL;
      }
      CMDprocessed = true; // command was processed
    }
    
    if (fcn.compare(EvalCmd["Disp0"]) == 0){
      if(params.size() == 1){      // ensure only one parameter -- transformation
        res = this->Disp(params[0],0, displayDescription);
      } else {
        wxString temp = wxT("ERROR - Wrong number of parameters for command: ");
        temp.Append(wxString::Format (wxT("%s"), cmd.c_str()));
        wxLogError(temp);
        res = NULL;
      }
      CMDprocessed = true; // command was processed
    }
    
    if (fcn.compare(EvalCmd["Diff"]) == 0){
      if (params.size() == 2){      // ensure only two parameters -- image1, image2
        res = this->Diff(params[0],params[1], displayDescription);
      } else {
        wxString temp = wxT( "ERROR - Wrong number of parameters for command: " );
        temp.Append(wxString::Format (wxT("%s"), cmd.c_str()));
        wxLogError(temp);
        res = NULL;
      }
      CMDprocessed = true; // command was processed
    }
    
    if (fcn.compare(EvalCmd["Transformation"]) == 0){
      if(params.size() == 3){      // ensure only three parameters -- source, target, algorithm
        res = this->ReadTransformation(params[0],params[1],params[2],displayDescription);
      } else {
        wxString temp = wxT("ERROR - Wrong number of parameters for command: ");
        temp.Append(wxString::Format (wxT("%s"), cmd.c_str()));
        wxLogError(temp);
        res = NULL;
      }
      CMDprocessed = true; // command was processed
    }
    
    if (fcn.compare(EvalCmd["SpatialData"]) == 0){
      if (params.size() == 2){      // ensure only two parameters -- datasetID, dataID
        res = this->ReadSpatialData(params[0],params[1],displayDescription);
      } else {
        wxString temp = wxT("ERROR - Wrong number of parameters for command: ");
        temp.Append(wxString::Format (wxT("%s"), cmd.c_str()));
        wxLogError(temp);
        res = NULL;
      }
      CMDprocessed = true; // command was processed
    }
    if (!CMDprocessed){
      wxString temp = wxT("ERROR - Command not known: ");
      temp.Append(wxString::Format (wxT("%s"), cmd.c_str()));
      wxLogError(temp);
      res = NULL;
    }
  }
  return res;
}


/* +--------------------+
   | Evaluator Commands |
   +--------------------+ */
 
/*
  function: ReadSpatialData
  purpose : reads spatial data from database
  input:    datasetID std::string containing the ID of the dataset (ex: "1")
            dataID    std::string containing the ID of the data (ex: "MRI")
  output  : gec::SpatialData* pointer of data 
*/
gec::SpatialData* NIREPEvaluator::ReadSpatialData(const std::string& datasetID,
                                                  const std::string& dataID, 
                                                  DisplayDescription * displayDescription)
{
  gec::SpatialData* res = NULL;

  // before we can ask if it's there, we must form the database manager command
  // const std::string cmd = evaluatorCommands.fListTwo.find(0)->second + "(" + datasetID + "," + dataID + ")"; // GEC removed
  const std::string cmd = EvalCmd["SpatialData"]  + "(" + datasetID + "," + dataID + ")"; // GEC added

  // first ask the database if it's got it in memory
  // note: request should increment the reference counter of cmd
  // gec::SpatialData* res = NIREPDataManager::request(cmd);
  res = this->dataManager->GetLoadedSpatialDataListEntry(cmd, displayDescription);

  // if it does not, request the DM load the image and then send it on its way
  if(res == NULL)
    if(this->dataManager->LoadSpatialData(datasetID, dataID, displayDescription))
      res = this->dataManager->GetLoadedSpatialDataListEntry(cmd, displayDescription);

  return res;
}


/*
  function: ReadTransformation
  purpose:  reads transformation from database
  input:    source std::string containing the name of the source (ex: "1")
            target std::string containing the name of the target (ex: "2")
            transformation std::string containing the name of the transformation
                           (ex: "SICLE_param1")
  output:   gec::SpatialData* pointer of data 
*/
gec::SpatialData* NIREPEvaluator::ReadTransformation(const std::string& source,
                                                     const std::string& target,
                                                     const std::string& transformation, 
                                                     DisplayDescription * displayDescription)
{
  gec::SpatialData* res = NULL;

  // before we can ask if it's there, we must form the database manager command
  //const std::string cmd =
  //  evaluatorCommands.fListTwo.find(1)->second+"(" + source + "," + target + "," + transformation + ")"; // GEC removed
  std::string cmd = EvalCmd["Transformation"]+"(" + source + "," + target + "," + transformation + ")"; // GEC added

  // first ask the database if it's got it in memory
  // note: request should increment the reference counter of cmd
  // gec::SpatialData* res = NIREPDataManager::request(cmd);
  res = this->dataManager->GetLoadedSpatialDataListEntry(cmd, displayDescription);

  // if it does not, request the DM load the image and then send it on its way
  if(res == NULL)
    if(this->dataManager->LoadTransformation(source, target, transformation, displayDescription))
      res = this->dataManager->GetLoadedSpatialDataListEntry(cmd, displayDescription);

  return res;
}


/*
  function: Diff
  purpose:  computes the absolute difference image1-image2
  input:    params std::vector<std::string> of the parameters
  output:   gec::SpatialData* pointer of image data 
*/
gec::SpatialData* NIREPEvaluator::Diff(const std::string& img1, const std::string& img2, DisplayDescription * displayDescription)
{
  EvaluatorList * taskList = displayDescription->GetEvaluatorList();
  gec::SpatialData* res = NULL;

  const gec::Image* im1 = dynamic_cast<const gec::Image *>(this->Evaluate(img1, displayDescription));
  if(im1==NULL) {
    return NULL;
  }
  const gec::Image* im2 = dynamic_cast<const gec::Image *>(this->Evaluate(img2, displayDescription));
  if(im2==NULL) {
    return NULL;
  }
  res = new gec::Image(*im1);
  dynamic_cast<gec::Image *>(res)->AbsDiff(im2);

  //std::string cmd = evaluatorCommands.fListTwo.find(2)->second+"(" + taskList->GetEvaluatorCommand(img1) + "," + taskList->GetEvaluatorCommand(img2) + ")"; // GEC removed
  std::string cmd = EvalCmd["Diff"]+"(" + taskList->GetEvaluatorCommand(img1) + "," + taskList->GetEvaluatorCommand(img2) + ")"; // GEC added
  this->dataManager->SetLoadedSpatialDataListEntry(cmd, res);

  return res;
}


/*
  function: Disp
  purpose:  Extracts ith-displacement field from a transformation
  input:    transformID Evaluator task ID (e.g. "_A1")
            index i-th index of transformation
  output:   gec::SpatialData* pointer of data 
*/
gec::SpatialData* NIREPEvaluator::Disp(const std::string& transformID, const unsigned int& index, DisplayDescription * displayDescription)
{
  EvaluatorList * taskList = displayDescription->GetEvaluatorList();
  gec::SpatialData* res = NULL;

  // obtain the transformation from the database
  const gec::Transformation* trans = dynamic_cast<const gec::Transformation *>(this->Evaluate(transformID, displayDescription));
  if(trans == NULL) {
    return NULL;
  }
  res = trans->GetDispImage(index);

  // NOTE: xDisp, yDisp and zDisp MUST be contiguous in the evaluatorTaskList for the next line to work!
  //std::string cmd = evaluatorCommands.fListTwo.find(3+index)->second + "(" + taskList->GetEvaluatorCommand(transformID) + ")"; // GEC removed
  std::string cmd;
  if (index == 0){
    cmd = EvalCmd["Disp0"] + "(" + taskList->GetEvaluatorCommand(transformID) + ")";  // GEC added
  } else if (index == 1){
    cmd = EvalCmd["Disp1"] + "(" + taskList->GetEvaluatorCommand(transformID) + ")";
  } else if (index == 2){
    cmd = EvalCmd["Disp2"] + "(" + taskList->GetEvaluatorCommand(transformID) + ")";
  }
  
  this->dataManager->SetLoadedSpatialDataListEntry(cmd, res);

  return res;
}


/*
  function: Jacobian
  purpose:  Computes the Jacobian of a transformation
  input:    transformID Evaluator task ID (e.g. "_A1")
  output:   gec::SpatialData* pointer of data 
*/
gec::SpatialData* NIREPEvaluator::Jacobian(const std::string& transformID, DisplayDescription * displayDescription)
{
  EvaluatorList * taskList = displayDescription->GetEvaluatorList();
  gec::SpatialData* res = NULL;

  // obtain the transformation from the database
  const gec::Transformation* trans = dynamic_cast<const gec::Transformation *>(this->Evaluate(transformID, displayDescription));
  if(trans == NULL) {
    return NULL;
  }

  //float min, max;
  //int min_x, min_y, min_z, max_x, max_y, max_z;
  gec::Jacobian* jac = new gec::Jacobian(trans);
  //res = jac->JacobianDiscrete3D(min, max, min_x, min_y, min_z, max_x, max_y, max_z, trans, 1.0F, 1.0F, 1.0F);
  res = jac->GetImage();
  delete jac;

  //std::string cmd = evaluatorCommands.fListTwo.find(6)->second+"(" + taskList->GetEvaluatorCommand(transformID) + ")";
  std::string cmd = EvalCmd["Jacobian"] + "(" + taskList->GetEvaluatorCommand(transformID) + ")";
  this->dataManager->SetLoadedSpatialDataListEntry(cmd, res);

  return res;
}


/*
  function: Transform
  purpose:  Transforms a spatial data with a transformation
  input:    imageID Evaluator task ID ("e.g. "_A1")
            transformID Evaluator task ID (e.g. "_A3")
  output:   gec::SpatialData* pointer of data 
*/
gec::SpatialData* NIREPEvaluator::Transform(const std::string& imageID, const std::string& transformID, DisplayDescription * displayDescription)
{
  EvaluatorList * taskList = displayDescription->GetEvaluatorList();
  gec::SpatialData* res = NULL;

  gec::SpatialData* data = this->Evaluate(imageID, displayDescription);
  if( data == NULL ) {
    return NULL;
  }
  gec::Transformation* trans = dynamic_cast<gec::Transformation *>(this->Evaluate(transformID, displayDescription));
  if( trans == NULL ) {
    return NULL;
  }

  switch(data->GetDataType()) {
    case gec::SpatialData::IMAGE:
      res = trans->Transform(data, gec::LINEAR);
      break;
    case gec::SpatialData::OBJECTMAP:
      res = trans->Transform(data, gec::NEARESTNEIGHBOR);
      break;
    case gec::SpatialData::TRANSFORMATION:
      res = trans->Transform(data, gec::LINEAR);
      break;
    default:
      std::cerr << "ERROR: UNSUPPORTED TRANSFORMATION TYPE" << std::endl;
      break;
  }

  /*std::string cmd = evaluatorCommands.fListTwo.find(7)->second+"(" + taskList->GetEvaluatorCommand(imageID) +
                                                   "," + taskList->GetEvaluatorCommand(transformID) + ")";*/
  std::string cmd = EvalCmd["TransformImage"]+"(" + taskList->GetEvaluatorCommand(imageID) +
                                                   "," + taskList->GetEvaluatorCommand(transformID) + ")";
  this->dataManager->SetLoadedSpatialDataListEntry(cmd, res);

  return res;
}

/*
  function: Transform
  purpose:  Transforms a spatial data with a transformation
  input:    imageID Evaluator task ID ("e.g. "_A1")
            transformID Evaluator task ID (e.g. "_A3")
  output:   gec::SpatialData* pointer of data 
*/
//orientation, width, color for grid, color for background, 
gec::SpatialData* NIREPEvaluator::TransformGrid(const std::string& transformID, DisplayDescription * displayDescription)
{
  
  EvaluatorList * taskList = displayDescription->GetEvaluatorList();
  gec::SpatialData* res = NULL;

  gec::Transformation* trans = dynamic_cast<gec::Transformation *>(this->Evaluate(transformID, displayDescription));
  if( trans == NULL ) {
    return NULL;
  }

  int gridSpace = 8;
  const  int *dimension = trans->GetDimension();
  gec::Image grid =  gec::Image(gec::Image::UCHAR, dimension);
  unsigned char * dataPointer = static_cast<unsigned char *> (grid.GetVTKData()->GetScalarPointer() );


  


  for(int k=0; k < dimension[2]; k++) {
    for(int j=0; j< dimension[1]; j++) {
      for(int i=0; i< dimension[0]; i++) {
        //need to add in for the other two orientations

        //if(i%gridSpace == 0 || j%gridSpace == 0) {
        //   //value at this location is 255;
        //  dataPointer[grid.GetIndex(i,j,k)] = 255;
        //}

        if(j%gridSpace == 0 || k%gridSpace == 0) {
           //value at this location is 255;
          dataPointer[grid.GetIndex(i,j,k)] = 255;
        }

        //if(i%gridSpace == 0 || k%gridSpace == 0) {
        //   //value at this location is 255;
        //  dataPointer[grid.GetIndex(i,j,k)] = 255;
        //}
      }
    }
  }


 res = trans->Transform(&grid, gec::LINEAR);

  /*std::string cmd = evaluatorCommands.fListTwo.find(7)->second+"(" + taskList->GetEvaluatorCommand(imageID) +
                                                   "," + taskList->GetEvaluatorCommand(transformID) + ")";*/
  std::string cmd = EvalCmd["TransformGrid"]+"(" + taskList->GetEvaluatorCommand(transformID) + ")";
  this->dataManager->SetLoadedSpatialDataListEntry(cmd, res);

  return res;
}

/*
  function: Transform
  purpose:  Transforms a spatial data with a transformation
  input:    imageID Evaluator task ID ("e.g. "_A1")
            transformID Evaluator task ID (e.g. "_A3")
  output:   gec::SpatialData* pointer of data 
*/
//orientation, width, color for grid, color for background, scale factor (float, adjusts where end point is) 
gec::SpatialData* NIREPEvaluator::VectorField(const std::string& transformID, DisplayDescription * displayDescription)
{
  
  EvaluatorList * taskList = displayDescription->GetEvaluatorList();
  //gec::SpatialData* res = NULL;

  gec::Transformation* trans = dynamic_cast<gec::Transformation *>(this->Evaluate(transformID, displayDescription));
  if( trans == NULL ) {
    return NULL;
  }

  int gridSpace = 8;
  const  int *dimension = trans->GetDimension();
  gec::Image *res =  new gec::Image(gec::Image::UCHAR, dimension);

  unsigned char * dataPointer = static_cast<unsigned char *> (res->GetVTKData()->GetScalarPointer() );
  float * transDataPointer = static_cast<float *>(trans->GetVTKData()->GetScalarPointer());


  //Debugging, remove after we solve orientation issues


  //for(int k=0; k < dimension[2]; k++) {
  //  for(int j=0; j< dimension[1]; j++) {
  //    for(int i=0; i< dimension[0]; i++) {
  //         //value at this location is 256;
  //        transDataPointer[trans->GetIndex(i,j,k,0)] = 0;
  //        transDataPointer[trans->GetIndex(i,j,k,1)] = (7*j)/dimension[1];
  //        transDataPointer[trans->GetIndex(i,j,k,2)] = (0);
  //    }
  //  }
  //}
  //-----------End of debugging

  for(int k=0; k < dimension[2]; k++) {
    for(int j=0; j< dimension[1]; j++) {
      for(int i=0; i< dimension[0]; i++) {
        //if(i == 0) {
        //   dataPointer[res->GetIndex(i,j,k)] = 255;
        //}
        if(false)
        {
        }
        if( (i%gridSpace == 0) && (j%gridSpace == 0) ) {
          
          // value at this location is 256;
          float fdx = transDataPointer[trans->GetIndex(i,j,k,0)] * 1;// ( (dimension[0]-1) /3);
          float fdy = transDataPointer[trans->GetIndex(i,j,k,1)] * 1;//( (dimension[1]-1)/3);
          int intfdx = (fdx  ) + i;
          int intfdy = (fdy  ) + j;

          //y= mx + b
          //m = fdy/fdx
          if(intfdx < 0) {
            intfdx = 0;
            intfdy = float(j) - float(i) * (fdy/fdx);
          }
          else if( intfdx > dimension[0] -1) {
            intfdx = dimension[0] -1;
            intfdy = (float((dimension[0] -1))*(fdy/fdx))+ (float(j) - float(i) * (fdy/fdx) );
          }


          if(intfdy < 0) {
            float temp = (fdy/fdx);
            float tempTwo = float(i) * (fdy/fdx);
            float tempThree = ( float(-j) + float(i) * (fdy/fdx) );
            intfdx = ( float(-j) + float(i) * (fdy/fdx) ) * (fdx/fdy);
            intfdy = 0;
          }
          else if( intfdy > dimension[1] -1) {
            intfdx = ( (dimension[1] -1) + (float(-j) + float(i) * (fdy/fdx)) ) * (fdx/fdy);
            intfdy = dimension[1] -1;
          }

          int x0 = i;
          int y0 = j;
          int dx = intfdx - i;
          int dy = intfdy - j;
          int stepx = -1;
          int stepy = -1;
          if (dy < 0) { dy = -dy;  stepy = -1; } else { stepy = 1; }
          if (dx < 0) { dx = -dx;  stepx = -1; } else { stepx = 1; }
          dy <<= 1;                                                  // dy is now 2*dy
          dx <<= 1;                                                  // dx is now 2*dx


          dataPointer[res->GetIndex(x0,y0,k)] = 255;
          if (dx > dy) {
              int fraction = dy - (dx >> 1);                         // same as 2*dy - dx
              while (x0 != (intfdx)) {
                  if (fraction >= 0) {
                      y0 += stepy;
                      fraction -= dx;                                // same as fraction -= 2*dx
                  }
                  x0 += stepx;
                  fraction += dy;                                    // same as fraction -= 2*dy
                  dataPointer[res->GetIndex(x0,y0,k)] = 255;
              }
          } else {
              int fraction = dx - (dy >> 1);
              while (y0 != (intfdy)) {
                  if (fraction >= 0) {
                      x0 += stepx;
                      fraction -= dy;
                  }
                  y0 += stepy;
                  fraction += dx;
                  dataPointer[res->GetIndex(x0,y0,k)] = 255;
              }
          }

   //       function line(x0, y0, x1, y1)
   //dx := abs(x1-x0)
   //dy := abs(y1-y0) 
   //if x0 < x1 then sx := 1 else sx := -1
   //if y0 < y1 then sy := 1 else sy := -1
   //err := dx-dy
 
   //loop
   //  setPixel(x0,y0)
   //  if x0 = x1 and y0 = y1 exit loop
   //  e2 := 2*err
   //  if e2 > -dy then 
   //    err := err - dy
   //    x0 := x0 + sx
   //  end if
   //  if e2 <  dx then 
   //    err := err + dx
   //    y0 := y0 + sy 
   //  end if
   //end loop
          
        }
      }
    }
  }


 //res = trans->Transform(&grid, gec::LINEAR);
  //res = &grid;

  /*std::string cmd = evaluatorCommands.fListTwo.find(7)->second+"(" + taskList->GetEvaluatorCommand(imageID) +
                                                   "," + taskList->GetEvaluatorCommand(transformID) + ")";*/
  std::string cmd = EvalCmd["VectorField"]+"(" + taskList->GetEvaluatorCommand(transformID) + ")";
  this->dataManager->SetLoadedSpatialDataListEntry(cmd, res);

  return res;
}


/*
  function: InvertTransformation
  purpose:  Computes the inverse of a transformation
  input:    transformID Evaluator task ID (e.g. "_A1")
  output:   gec::SpatialData* pointer of data 
*/
gec::SpatialData* NIREPEvaluator::InvertTransformation(const std::string& transformID, DisplayDescription * displayDescription)
{
  EvaluatorList * taskList = displayDescription->GetEvaluatorList();
  gec::SpatialData* res = NULL;

  // obtain the transformation from the database
  gec::Transformation* trans = dynamic_cast<gec::Transformation *>(this->Evaluate(transformID, displayDescription));
  if(trans == NULL) {
    return NULL;
  }
  // prepare invert transformation
  res = trans->Invert();

  //std::string cmd = evaluatorCommands.fListTwo.find(8)->second+"(" + taskList->GetEvaluatorCommand(transformID) + ")";
  std::string cmd = EvalCmd["InvertTransformation"] +"(" + taskList->GetEvaluatorCommand(transformID) + ")";
  this->dataManager->SetLoadedSpatialDataListEntry(cmd, res);

  return res;
}


/*
  function: InverseConsistencyError
  purpose:  Computes the inverse consistency error of a pair transformations
  input:    fwdTransID Evaluator task ID (e.g. "_A1")
            revTransID Evaluator task ID (e.g. "_A2")
            method     Computation method (e.g. inv/comp)
  output:   gec::SpatialData* pointer of data 
*/
gec::SpatialData* NIREPEvaluator::InverseConsistencyError(const std::string& fwdTransID,
                                                          const std::string& revTransID,
                                                          const std::string& method,
                                                          DisplayDescription * displayDescription)
{
  EvaluatorList * taskList = displayDescription->GetEvaluatorList();
  gec::SpatialData* res = NULL;

  gec::Transformation* u12 = dynamic_cast<gec::Transformation *>(this->Evaluate(fwdTransID, displayDescription));
  if(u12 == NULL) {
    return NULL;
  }
  NIREPInverseConsistencyError* ice = new NIREPInverseConsistencyError(this);
  if(method == "inv") {
    // || h_12(x) - h_21^-1(x) ||^2 = || u_12(x) - u_21^-1(x) ||^2
    gec::Transformation* u21i =
      dynamic_cast<gec::Transformation *>(this->InvertTransformation(revTransID, displayDescription));
    if(u21i == NULL) {
      return NULL;
    }
    res = ice->ComputeICEImage(u12, u21i);
  } else {
    // || h_21(h_12(x)) - x ||^2 = || h_21(x+u_12(x)) - x ||^2
    //                           = || u_12(x) + u_21(x+u_12(x)) ||^2
    gec::Transformation* u121 =
      dynamic_cast<gec::Transformation *>(this->Transform(fwdTransID, revTransID, displayDescription));
    if(u121 == NULL) {
      return NULL;
    }
    res = ice->ComputeICEImage(u121);
  }
  delete ice;

  /*std::string cmd = evaluatorCommands.fListTwo.find(9)->second+"(" + taskList->GetEvaluatorCommand(fwdTransID) +
                                                   "," + taskList->GetEvaluatorCommand(revTransID) +
                                                   "," + method + ")";*/
  std::string cmd = EvalCmd["InverseConsistencyErrorImage"]+"(" + taskList->GetEvaluatorCommand(fwdTransID) +
                                                   "," + taskList->GetEvaluatorCommand(revTransID) +
                                                   "," + method + ")";
  this->dataManager->SetLoadedSpatialDataListEntry(cmd, res);

  return res;
}


/*
  function: TransitivityError
  purpose:  Computes the transitivity error of transformations
  input:    transABID Evaluator task ID (ex: "_A1")
            transBCID Evaluator task ID (ex: "_A2")
            transCAID Evaluator task ID (ex: "_A3")
            method    Computation method (e.g. inv/comp)
  output:   gec::SpatialData* pointer of data 
*/
gec::SpatialData* NIREPEvaluator::TransitivityError(const std::string& transABID,
                                                    const std::string& transBCID,
                                                    const std::string& transACID,
                                                    const std::string& method,
                                                    DisplayDescription * displayDescription)
{
  EvaluatorList * taskList = displayDescription->GetEvaluatorList();
  gec::SpatialData* res = NULL;

  gec::Transformation* abc =
    dynamic_cast<gec::Transformation *>(this->Transform(transBCID, transABID, displayDescription));
  if(abc == NULL) {
    return NULL;
  }
  gec::Transformation* ac =
    dynamic_cast<gec::Transformation *>(this->Evaluate(transACID, displayDescription));
  if(ac == NULL) {
    return NULL;
  }

  NIREPTransitivityError* te = new NIREPTransitivityError(this);
  if(method == "inv") {
    // || h_12(h_23(x)) - h_13(x) ||^2
    // = || u_23(x) + u_12(x+u_23(x)) - u_13(x) ||^2
    res = te->ComputeTEImage(abc, ac);
  } else {
    // || h_31(h_12(h_23(x))) - x ||^2
    // = || u_23(x) + u_12(x+u_23(x)) + u_31(x+u_23(x)+u_12(x+u_23(x))) ||^2
    // Perhaps there's a better way to do this, but this works for now
    std::string cmd = EvalCmd["TransformImage"] + "("
                    + EvalCmd["TransformImage"] + "("
                    + taskList->GetEvaluatorCommand(transBCID) + ","
                    + taskList->GetEvaluatorCommand(transABID) + "),"
                    + taskList->GetEvaluatorCommand(transACID) + ")";
    gec::Transformation* abca =
      dynamic_cast<gec::Transformation *>(ac->Transform(abc, gec::LINEAR));
    if(abca == NULL) {
      return NULL;
    }
    this->dataManager->SetLoadedSpatialDataListEntry(cmd, abca);
    res = te->ComputeTEImage(abca);
  }
  delete te;

  //std::string cmd = evaluatorCommands.fListTwo.find(10)->second+"(" + taskList->GetEvaluatorCommand(transABID) +
  //                                                  "," + taskList->GetEvaluatorCommand(transBCID) +
  //                                                  "," + taskList->GetEvaluatorCommand(transACID) +
  //                                                  "," + method + ")";
  std::string cmd = EvalCmd["TransitivityErrorImage"]+"(" + taskList->GetEvaluatorCommand(transABID) +
                                                    "," + taskList->GetEvaluatorCommand(transBCID) +
                                                    "," + taskList->GetEvaluatorCommand(transACID) +
                                                    "," + method + ")";
  this->dataManager->SetLoadedSpatialDataListEntry(cmd, res);

  return res;
}

/*
  function: IntensityVariance
  purpose:  Computes the intensity variance of a range of transformations
  input:    arg vector of strings (ex "1,[2:3;4:8],014, MRI")
  output:   gec::SpatialData* pointer of data 
*/
gec::SpatialData* NIREPEvaluator::IntensityVariance(std::vector<std::string>& arg,
                                                    std::string version,
                                                    DisplayDescription * displayDescription)
{
  gec::SpatialData* res = NULL;

  NIREPIntensityVariance* iv = new NIREPIntensityVariance(this);
  arg.insert(arg.begin(),version);
  res = iv->ComputeIVImage(arg,displayDescription);
  std::string cmd;

  //if(!arg[0].compare("RegisteredImageVariance")) {
  // cmd = evaluatorCommands.fListTwo.find(11)->second+"(" + arg[1] +
  //                                                 "," + arg[2] +
  //                                                 "," + arg[3] + 
  //                                                 "," + arg[4] + ")";
  //} else if(!arg[0].compare("ImageVariance")) {
  //  cmd = evaluatorCommands.fListTwo.find(12)->second+"(" + arg[1] +
  //                                                  "," + arg[2] + ")";
  //}
  if(!arg[0].compare(EvalCmd["IntensityVarianceTransform"])) {
   cmd = EvalCmd["IntensityVarianceTransform"]+"(" + arg[1] +
                                               "," + arg[2] +
                                               "," + arg[3] + 
                                               "," + arg[4] + ")";
  } else if(!arg[0].compare(EvalCmd["IntensityVarianceImage"])) {
    cmd = EvalCmd["IntensityVarianceImage"]+"(" + arg[1] +
                                            "," + arg[2] + ")";
  }
  this->dataManager->SetLoadedSpatialDataListEntry(cmd, res);

  return res;
}