NIREPIntensityVariance.cxx

Go to the documentation of this file.

#include "NIREPIntensityVariance.h"

//Gets called from the Evaluator to compute the Statistic
TextTable NIREPIntensityVariance::ComputeStatistic(const std::vector<std::string>& arg, DisplayDescription * displayDescription)
{
  TextTable resultTextTable;
  this->Expand(arg);
  this->Compute(arg,displayDescription);

  std::vector<float> values;

  int *dims = result->GetVTKData()->GetDimensions();
  const int size = dims[0] * dims[1] * dims[2];

  float * resultData = static_cast<float*>(result->GetVTKData()->GetScalarPointer());
  for(int i=0;i<size;i++)
  {
    if(resultData[i] != 0) {
      int debugStop = 0;
    }
    values.push_back( resultData[i] );
  }
  
  int debugStop = 0;
  resultTextTable.Value.push_back(values);

  return resultTextTable;
}

//Gets called from the Evaluator to compute the image
gec::SpatialData* NIREPIntensityVariance::ComputeIVImage(const std::vector<std::string>& arg, DisplayDescription * displayDescription)
{
  this->Expand(arg);
  this->Compute(arg,displayDescription);
  return result;
}

//mean: u = E{x} = 1/N * sum(x_i) from i=1 to N, where N is the number of images
//Standard deviation: = E{x^2} = 1/N * sum( (x_i)^2 ) from i=1 to N, where N is the number of images
//Variance: E{ (x-u)^2 } = E{x^2} - u^2
void NIREPIntensityVariance::ImageVariance(const std::vector<std::string>& arg, DisplayDescription * displayDescription) {
    //intialize the mean image to all zeros and to the corresponding dimensions, spacing and we
    //want a float image to do computations.
    if(this->evaluator->ReadSpatialData(datasetID[0],arg[3],displayDescription) == NULL) {
      return;
    }
    gec::Image *mean = new gec::Image(gec::Image::FLOAT,
                                      dynamic_cast<gec::Image *>(this->evaluator->ReadSpatialData(datasetID[0],arg[2],displayDescription))->GetVTKData()->GetDimensions(),
                                      dynamic_cast<gec::Image *>(this->evaluator->ReadSpatialData(datasetID[0],arg[2],displayDescription))->GetVTKData()->GetSpacing());
    if(mean == NULL) {
      return;
    }
    //intialize the result image to the same parameters as mean
    result = new gec::Image( *mean );
    if(result == NULL) {
      return;
    }

    //get the scalar pointers because that is how we have to do the computations
    float *meanData = static_cast<float*>(mean->GetVTKData()->GetScalarPointer());
    float *resultData = static_cast<float*>(result->GetVTKData()->GetScalarPointer());

    //get the dimensions and spacing
    int *dims = mean->GetVTKData()->GetDimensions();
    const int size = dims[0] * dims[1] * dims[2];

    //run through a for loop for each data that the user wants to use 
    //to figure out the variance
    for(int h=0;h<datasetID.size();h++)
    {
      gec::Image *temp3 = dynamic_cast<gec::Image *>(this->evaluator->ReadSpatialData(datasetID[h],arg[2],displayDescription));
      if(temp3 == NULL) {
        return;
      }
      unsigned char *temp3Data = static_cast<unsigned char*>(temp3->GetVTKData()->GetScalarPointer());

      //run through a for loop for each data point in the images to do the computations
      for(int i=0;i<size;i++)
      {
          meanData[i] = meanData[i] + static_cast<float>(temp3Data[i]);
          resultData[i] = resultData[i] + static_cast<float>(temp3Data[i]) * static_cast<float>(temp3Data[i]);
      }

      //release the data from the data manager
      //this way we are conserving memory
      const std::string cmd =  "SpatialData(" + datasetID[h] + "," + arg[2] + ")";
      this->evaluator->GetDataManager()->ReleaseLoadedSpatialDataListEntry(cmd);
    }

    //Do our final computations
    mean->Mul(1.0f/static_cast<float>(datasetID.size()));
    mean->Mul(mean);
    result->Mul(1.0f/static_cast<float>(datasetID.size()));
    result->Sub(mean);
}


//mean: u = E{x} = 1/N * sum(x_i) from i=1 to N, where N is the number of images
//Standard deviation: = E{x^2} = 1/N * sum( (x_i)^2 ) from i=1 to N, where N is the number of images
//Variance: E{ (x-u)^2 } = E{x^2} - u^2
void NIREPIntensityVariance::RegisteredImageVariance(const std::vector<std::string>& arg, DisplayDescription * displayDescription) {
    //need to read in corresponding mri image and call transform  

    //intialize the mean image to all zeros and to the corresponding dimensions, spacing and we
    //want a float image to do computations.
    if(this->evaluator->ReadSpatialData(datasetID[0],arg[3],displayDescription) == NULL) {
      return;
    }
    gec::Image *mean = new gec::Image(gec::Image::FLOAT,
                                      dynamic_cast<gec::Image *>(this->evaluator->ReadSpatialData(datasetID[0],arg[3],displayDescription))->GetVTKData()->GetDimensions(),
                                      dynamic_cast<gec::Image *>(this->evaluator->ReadSpatialData(datasetID[0],arg[3],displayDescription))->GetVTKData()->GetSpacing());
    if(mean == NULL) {
      return;
    }
    //intialize the result image to the same parameters as mean
    result = new gec::Image( *mean );
    if(result == NULL) {
      return;
    }
      

    //get the dimensions and spacing
    int *dims = mean->GetVTKData()->GetDimensions();
    const int size = dims[0] * dims[1] * dims[2];

    //get the scalar pointers because that is how we have to do the computations
    float * resultData = static_cast<float*>(result->GetVTKData()->GetScalarPointer());
    float * meanData = static_cast<float*>(mean->GetVTKData()->GetScalarPointer());

    //run through a for loop for each data that the user wants to use 
    //to figure out the variance
    for(int h=1;h<datasetID.size();h++)
    {
      gec::Transformation *temp3Trans = dynamic_cast<gec::Transformation *>(this->evaluator->ReadTransformation(datasetID[h],datasetID[0],arg[4],displayDescription));
      if(temp3Trans == NULL) {
        return;
      }
      gec::Image *temp3Image = new gec::Image(*dynamic_cast<gec::Image *>(this->evaluator->ReadSpatialData(datasetID[h],arg[3],displayDescription)));
      if(temp3Image == NULL) {
        return;
      }


      gec::Image *temp3 = dynamic_cast<gec::Image *>(temp3Trans->Transform(temp3Image, gec::LINEAR));
      if(temp3 == NULL) {
        return;
      }

      unsigned char *temp3Data = static_cast<unsigned char*>(temp3->GetVTKData()->GetScalarPointer());

      //run through a for loop for each data point in the images to do the computations
      for(int i=0;i<size;i++)
      {
          meanData[i] = meanData[i] + static_cast<float>(temp3Data[i]);
          resultData[i] = resultData[i] + static_cast<float>(temp3Data[i]) * static_cast<float>(temp3Data[i]);
      }

      //release the data from the data manager
      //this way we are conserving memory
      const std::string cmd =  "Transformation("+ datasetID[h] + "," + datasetID[0] + "," + arg[4] + ")";;
      this->evaluator->GetDataManager()->ReleaseLoadedSpatialDataListEntry(cmd);

      const std::string cmdTwo =  "SpatialData(" + datasetID[h] + "," + arg[3] + ")";
      this->evaluator->GetDataManager()->ReleaseLoadedSpatialDataListEntry(cmdTwo);
    }

    //Do our final computations
    mean->Mul(1.0f/static_cast<float>(datasetID.size()-1));
    mean->Mul(mean);
    result->Mul(1.0f/static_cast<float>(datasetID.size()-1));
    result->Sub(mean);
}

void NIREPIntensityVariance::Compute(const std::vector<std::string>& arg, DisplayDescription * displayDescription)
{
  if(!arg[0].compare("ImageVariance")) {
    this->ImageVariance(arg,displayDescription); 
  }
  else if(!arg[0].compare("RegisteredImageVariance")) {
    this->RegisteredImageVariance(arg,displayDescription);
  }
}

//Expand out the part that is something like: [14;4:5;8:11] to a vector
//that contains [14,4,5,8,9,10,11]
void NIREPIntensityVariance::Expand(const std::vector<std::string>& arg)
{
  //First we pull out the part that we want to expand from the argument list
  std::string stringToExpand;
  if(!arg[0].compare("RegisteredImageVariance")) {
    std::stringstream outTwo;
    if(atoi(arg[1].c_str()) < 10)
      outTwo<<"00";
    else if(atoi(arg[1].c_str())<100)
      outTwo<<"0";
    outTwo << atoi(arg[1].c_str());
    datasetID.push_back(outTwo.str());
    stringToExpand = arg[2];
    stringToExpand.erase(0,1);
    stringToExpand.erase(stringToExpand.find_last_of("]"),1);
    
  }
  else if(!arg[0].compare("ImageVariance")) {
    stringToExpand = arg[1];
    stringToExpand.erase(0,1);
    stringToExpand.erase(stringToExpand.find_last_of("]"),1);
  }
  

  //This is where we begin the process of expanding

  //first we have to have a char * to pass to strtok
  //Since we are using strings, and c_str() returns
  //a const char*, we are going to copy the data from
  //c_str to a new variable called strtokStringToExpand
  //which will be passed to strtok to tokenize bassed
  //on ";"
  char * strtokStringToExpand = new char[stringToExpand.size()+1]();
  memcpy(strtokStringToExpand,stringToExpand.c_str(),stringToExpand.size()+1);
  char* number = strtok (strtokStringToExpand,";");

  //keep tokenizing until we hit the end of the string
  while (number != NULL)
  {
    //Figure out if there is a : in the part that we
    //tokenized
    std::string numberTwo = number;
    size_t isColon = numberTwo.find_first_of(":");
    std::string firstNumber;
    std::string secondNumber;
    //if there is a : then we figure out the first
    //number to expand from and the second number to expand to
    if(isColon != std::string::npos) {
      firstNumber = std::string(numberTwo,0,isColon);
      secondNumber = std::string(numberTwo,isColon+1, numberTwo.size()-isColon+1);
    }
    //if there is not a : then we set the first and second number
    //to the same value
    else {
      firstNumber = numberTwo;
      secondNumber = numberTwo;
    }
    //figure out the true number of the first number in terms
    //of integer
    int trueNumber = atoi(firstNumber.c_str());
    for(int i=trueNumber; i<= atoi(secondNumber.c_str());i++)
    {
      //since when we store the values as three bits
      //in the database, we have to make sure that
      //the string has 3 bits based on the numbers
      //that is within the display descrption
      std::stringstream out;
      if(i < 10)
        out<<"00";
      else if(i<100)
        out<<"0";
      out << i;
      datasetID.push_back(out.str());

    }
    //get the next tokenized part
    number = strtok (NULL, ";");
  }
}