<commit_before>/*=========================================================================

  Library   : Image Registration Toolkit (IRTK)
  Module    : $Id$
  Copyright : Imperial College, Department of Computing
              Visual Information Processing (VIP), 2008 onwards
  Date      : $Date$
  Version   : $Revision$
  Changes   : $Author$

=========================================================================*/

#include <irtkRegistration.h>

// Used as temporary memory for transformed intensities
extern irtkGreyImage *_tmpImage;

// Used as temporary memory for metric
extern irtkSimilarityMetric *_tmpMetricA, *_tmpMetricB;

// Used as lookup table for transformed coordinates up to level n-1.  This
// lookup table needs to be calculated only once for each image resolution
// level.
extern float *_mffdLookupTable;

// Used as lookup table for transformed coordinates including level n. This
// lookup table needs to be calculated each time a control point has been
// modified.
extern float *_affdLookupTable;

// Used as lookup table for the contribution of each control point. This
// lookup table needs to be calculated only once.
extern float *_localLookupTable;

void irtkImageFreeFormRegistration2D::GuessParameter()
{
  int i;
  double xsize, ysize, zsize, spacing;

  if ((_target == NULL) || (_source == NULL)) {
    cerr << "irtkImageFreeFormRegistration2D::GuessParameter: Target and source image not found" << endl;
    exit(1);
  }

  // Default parameters for registration
  _NumberOfLevels     = 3;
  _NumberOfBins       = 64;

  // Default parameters for optimization
  _SimilarityMeasure  = NMI;
  _OptimizationMethod = GradientDescent;
  _Epsilon            = 0.0001;

  // Read target pixel size
  _target->GetPixelSize(&xsize, &ysize, &zsize);

  // Use xsize as spacing
  spacing = xsize;

  // Default target parameters
  _TargetBlurring[0]      = GuessResolution(xsize, ysize) / 2.0;
  _TargetResolution[0][0] = GuessResolution(xsize, ysize);
  _TargetResolution[0][1] = GuessResolution(xsize, ysize);
  _TargetResolution[0][2] = zsize;

  for (i = 1; i < _NumberOfLevels; i++) {
    _TargetBlurring[i]      = _TargetBlurring[i-1] * 2;
    _TargetResolution[i][0] = _TargetResolution[i-1][0] * 2;
    _TargetResolution[i][1] = _TargetResolution[i-1][1] * 2;
    _TargetResolution[i][2] = zsize;
  }

  // Read source pixel size
  _source->GetPixelSize(&xsize, &ysize, &zsize);

  // Default source parameters
  _SourceBlurring[0]      = GuessResolution(xsize, ysize, zsize) / 2.0;
  _SourceResolution[0][0] = GuessResolution(xsize, ysize);
  _SourceResolution[0][1] = GuessResolution(xsize, ysize);
  _SourceResolution[0][2] = zsize;

  for (i = 1; i < _NumberOfLevels; i++) {
    _SourceBlurring[i]      = _SourceBlurring[i-1] * 2;
    _SourceResolution[i][0] = _SourceResolution[i-1][0] * 2;
    _SourceResolution[i][1] = _SourceResolution[i-1][1] * 2;
    _SourceResolution[i][2] = zsize;
  }

  // Default parameters for non-rigid registration
  _Lambda1            = 0;
  _Lambda2            = 0;
  _Lambda3            = 0;
  _DX                 =_target->GetX() * spacing / 10.0;
  _DY                 =_target->GetX() * spacing / 10.0;
  _DZ                 = 1;
  _Subdivision        = true;

  // Remaining parameters
  for (i = 0; i < _NumberOfLevels; i++) {
    _NumberOfIterations[i] = 10;
    _NumberOfSteps[i]      = 4;
    _LengthOfSteps[i]      = _DX / 8.0 * pow(2.0, i);
  }

  // Try to guess padding by looking at voxel values in all eight corners of the volume:
  // If all values are the same we assume that they correspond to the padding value
  _TargetPadding = MIN_GREY;
  if ((_target->Get(_target->GetX()-1, 0, 0)                 == _target->Get(0, 0, 0)) &&
      (_target->Get(0, _target->GetY()-1, 0)                 == _target->Get(0, 0, 0)) &&
      (_target->Get(_target->GetX()-1, _target->GetY()-1, 0) == _target->Get(0, 0, 0))) {
    _TargetPadding = _target->Get(0, 0, 0);
  }
}

double irtkImageFreeFormRegistration2D::Evaluate()
{
  // Image coordinates
  int i, j;
  // World coordinates
  double x, y, z;
  // Pointer to reference data
  irtkGreyPixel *ptr2target;
  irtkGreyPixel *ptr2tmp;
  float *ptr;

  // Print debugging information
  this->Debug("irtkImageFreeFormRegistration2D::Evaluate");

  // Initialize histogram
  _metric->Reset();

  // Loop over all voxels in the target (reference) volume
  ptr2target = _target->GetPointerToVoxels();
  ptr2tmp    = _tmpImage->GetPointerToVoxels();
  ptr        = _mffdLookupTable;;
  for (j = 0; j < _target->GetY(); j++) {
    for (i = 0; i < _target->GetX(); i++) {
      // Check whether reference point is valid
      if (*ptr2target >= 0) {
        x = i;
        y = j;
        z = 0;
        _target->ImageToWorld(x, y, z);
        _affd->LocalDisplacement(x, y, z);
        x += ptr[0];
        y += ptr[1];
        _source->WorldToImage(x, y, z);
        // Check whether transformed point is inside volume
        if ((x > _source_x1) && (x < _source_x2) &&
            (y > _source_y1) && (y < _source_y2)) {
          // Add sample to metric
          *ptr2tmp =  round(_interpolator->EvaluateInside(x, y, 0));
          _metric->Add(*ptr2target, *ptr2tmp);
        } else {
          *ptr2tmp = -1;
        }
      }
      // Increment pointers to next voxel
      ptr2tmp++;
      ptr2target++;
      ptr += 3;
    }
  }

  // Evaluate similarity measure
  double similarity = _metric->Evaluate();

  // Add penalty for smoothness
  if (this->_Lambda1 > 0) {
    similarity += this->_Lambda1*this->SmoothnessPenalty();
  }
  // Add penalty for volume preservation
  if (this->_Lambda2 > 0) {
    similarity += this->_Lambda2*this->VolumePreservationPenalty();
  }
  // Add penalty for topology preservation
  if (this->_Lambda3 > 0) {
    similarity += this->_Lambda3*this->TopologyPreservationPenalty();
  }

  // Return similarity measure + penalty terms
  return similarity;
}

double irtkImageFreeFormRegistration2D::EvaluateDerivative(int index, double step)
{
  float *ptr;
  double bi, bj, bk, dx, dy, p[3];
  int i, j, i1, i2, j1, j2, k1, k2, dim;
  irtkPoint p1, p2;
  irtkGreyPixel *ptr2target, *ptr2tmp;

  // Print debugging information
  this->Debug("irtkImageFreeFormRegistration2D::EvaluateDerivative(int, double)");

  // Initialize metrics for forward and backward derivative steps
  _tmpMetricA->Reset(_metric);
  _tmpMetricB->Reset(_metric);

  // Calculate bounding box of control point in world coordinates
  _affd->BoundingBox(index, p1, p2);
  _target->WorldToImage(p1);
  _target->WorldToImage(p2);

  // Calculate bounding box of control point in image coordinates
  _affd->BoundingBox(_target, index, i1, j1, k1, i2, j2, k2, this->_SpeedupFactor);

  // Calculate incremental changes in lattice coordinates when looping
  // over target
  dx = (FFDLOOKUPTABLESIZE-1)/(p2._x-p1._x);
  dy = (FFDLOOKUPTABLESIZE-1)/(p2._y-p1._y);

  // Calculate whether this DOF corresponds to x, y or z-displacement<commit_msg>fixed a invalid pointer for caused by speadup factor in 2D FFD registration<commit_after>/*=========================================================================

  Library   : Image Registration Toolkit (IRTK)
  Module    : $Id$
  Copyright : Imperial College, Department of Computing
              Visual Information Processing (VIP), 2008 onwards
  Date      : $Date$
  Version   : $Revision$
  Changes   : $Author$

=========================================================================*/

#include <irtkRegistration.h>

// Used as temporary memory for transformed intensities
extern irtkGreyImage *_tmpImage;

// Used as temporary memory for metric
extern irtkSimilarityMetric *_tmpMetricA, *_tmpMetricB;

// Used as lookup table for transformed coordinates up to level n-1.  This
// lookup table needs to be calculated only once for each image resolution
// level.
extern float *_mffdLookupTable;

// Used as lookup table for transformed coordinates including level n. This
// lookup table needs to be calculated each time a control point has been
// modified.
extern float *_affdLookupTable;

// Used as lookup table for the contribution of each control point. This
// lookup table needs to be calculated only once.
extern float *_localLookupTable;

void irtkImageFreeFormRegistration2D::GuessParameter()
{
  int i;
  double xsize, ysize, zsize, spacing;

  if ((_target == NULL) || (_source == NULL)) {
    cerr << "irtkImageFreeFormRegistration2D::GuessParameter: Target and source image not found" << endl;
    exit(1);
  }

  // Default parameters for registration
  _NumberOfLevels     = 3;
  _NumberOfBins       = 64;

  // Default parameters for optimization
  _SimilarityMeasure  = NMI;
  _OptimizationMethod = GradientDescent;
  _Epsilon            = 0.0001;

  // Read target pixel size
  _target->GetPixelSize(&xsize, &ysize, &zsize);

  // Use xsize as spacing
  spacing = xsize;

  // Default target parameters
  _TargetBlurring[0]      = GuessResolution(xsize, ysize) / 2.0;
  _TargetResolution[0][0] = GuessResolution(xsize, ysize);
  _TargetResolution[0][1] = GuessResolution(xsize, ysize);
  _TargetResolution[0][2] = zsize;

  for (i = 1; i < _NumberOfLevels; i++) {
    _TargetBlurring[i]      = _TargetBlurring[i-1] * 2;
    _TargetResolution[i][0] = _TargetResolution[i-1][0] * 2;
    _TargetResolution[i][1] = _TargetResolution[i-1][1] * 2;
    _TargetResolution[i][2] = zsize;
  }

  // Read source pixel size
  _source->GetPixelSize(&xsize, &ysize, &zsize);

  // Default source parameters
  _SourceBlurring[0]      = GuessResolution(xsize, ysize, zsize) / 2.0;
  _SourceResolution[0][0] = GuessResolution(xsize, ysize);
  _SourceResolution[0][1] = GuessResolution(xsize, ysize);
  _SourceResolution[0][2] = zsize;

  for (i = 1; i < _NumberOfLevels; i++) {
    _SourceBlurring[i]      = _SourceBlurring[i-1] * 2;
    _SourceResolution[i][0] = _SourceResolution[i-1][0] * 2;
    _SourceResolution[i][1] = _SourceResolution[i-1][1] * 2;
    _SourceResolution[i][2] = zsize;
  }

  // Default parameters for non-rigid registration
  _Lambda1            = 0;
  _Lambda2            = 0;
  _Lambda3            = 0;
  _DX                 =_target->GetX() * spacing / 10.0;
  _DY                 =_target->GetX() * spacing / 10.0;
  _DZ                 = 1;
  _Subdivision        = true;

  // Remaining parameters
  for (i = 0; i < _NumberOfLevels; i++) {
    _NumberOfIterations[i] = 10;
    _NumberOfSteps[i]      = 4;
    _LengthOfSteps[i]      = _DX / 8.0 * pow(2.0, i);
  }

  // Try to guess padding by looking at voxel values in all eight corners of the volume:
  // If all values are the same we assume that they correspond to the padding value
  _TargetPadding = MIN_GREY;
  if ((_target->Get(_target->GetX()-1, 0, 0)                 == _target->Get(0, 0, 0)) &&
      (_target->Get(0, _target->GetY()-1, 0)                 == _target->Get(0, 0, 0)) &&
      (_target->Get(_target->GetX()-1, _target->GetY()-1, 0) == _target->Get(0, 0, 0))) {
    _TargetPadding = _target->Get(0, 0, 0);
  }
}

double irtkImageFreeFormRegistration2D::Evaluate()
{
  // Image coordinates
  int i, j;
  // World coordinates
  double x, y, z;
  // Pointer to reference data
  irtkGreyPixel *ptr2target;
  irtkGreyPixel *ptr2tmp;
  float *ptr;

  // Print debugging information
  this->Debug("irtkImageFreeFormRegistration2D::Evaluate");

  // Initialize histogram
  _metric->Reset();

  // Loop over all voxels in the target (reference) volume
  ptr2target = _target->GetPointerToVoxels();
  ptr2tmp    = _tmpImage->GetPointerToVoxels();
  ptr        = _mffdLookupTable;;
  for (j = 0; j < _target->GetY(); j++) {
    for (i = 0; i < _target->GetX(); i++) {
      // Check whether reference point is valid
      if (*ptr2target >= 0) {
        x = i;
        y = j;
        z = 0;
        _target->ImageToWorld(x, y, z);
        _affd->LocalDisplacement(x, y, z);
        x += ptr[0];
        y += ptr[1];
        _source->WorldToImage(x, y, z);
        // Check whether transformed point is inside volume
        if ((x > _source_x1) && (x < _source_x2) &&
            (y > _source_y1) && (y < _source_y2)) {
          // Add sample to metric
          *ptr2tmp =  round(_interpolator->EvaluateInside(x, y, 0));
          _metric->Add(*ptr2target, *ptr2tmp);
        } else {
          *ptr2tmp = -1;
        }
      }
      // Increment pointers to next voxel
      ptr2tmp++;
      ptr2target++;
      ptr += 3;
    }
  }

  // Evaluate similarity measure
  double similarity = _metric->Evaluate();

  // Add penalty for smoothness
  if (this->_Lambda1 > 0) {
    similarity += this->_Lambda1*this->SmoothnessPenalty();
  }
  // Add penalty for volume preservation
  if (this->_Lambda2 > 0) {
    similarity += this->_Lambda2*this->VolumePreservationPenalty();
  }
  // Add penalty for topology preservation
  if (this->_Lambda3 > 0) {
    similarity += this->_Lambda3*this->TopologyPreservationPenalty();
  }

  // Return similarity measure + penalty terms
  return similarity;
}

double irtkImageFreeFormRegistration2D::EvaluateDerivative(int index, double step)
{
  float *ptr;
  double bi, bj, bk, dx, dy, p[3];
  int i, j, i1, i2, j1, j2, k1, k2, dim;
  irtkPoint p1, p2;
  irtkGreyPixel *ptr2target, *ptr2tmp;

  // Print debugging information
  this->Debug("irtkImageFreeFormRegistration2D::EvaluateDerivative(int, double)");

  // Initialize metrics for forward and backward derivative steps
  _tmpMetricA->Reset(_metric);
  _tmpMetricB->Reset(_metric);

  // Calculate bounding box of control point in world coordinates
  _affd->BoundingBox(index, p1, p2);
  _target->WorldToImage(p1);
  _target->WorldToImage(p2);

  // Calculate bounding box of control point in image coordinates
  _affd->BoundingBox(_target, index, i1, j1, k1, i2, j2, k2, 1.0 / _SpeedupFactor);

  // Calculate incremental changes in lattice coordinates when looping
  // over target
  dx = (FFDLOOKUPTABLESIZE-1)/(p2._x-p1._x);
  dy = (FFDLOOKUPTABLESIZE-1)/(p2._y-p1._y);

  // Calculate whether this DOF corresponds to x, y or z-displacement
{'fraction_non_alphanumeric': 0.14324401368301026, 'fraction_numerical': 0.023588939566704676, 'mean_word_length': 3.1032163742690058, 'pattern_counts': {'":': 0, '<': 27, '<?xml version=': 0, '>': 104, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '13371', 'n_tokens_mistral': 4858, 'n_tokens_neox': 4531, 'n_words': 1244}