void R2Image::
blendOtherImageTranslated(R2Image * otherImage)
{
  R2Image *output = new R2Image(*otherImage);
	std::vector<Feature> features = *(this->Harris(3));
  int searchSpaceXDim = this->Width() / 10; // half the search space dimension
  int searchSpaceYDim = this->Height() / 10;
  int windowDimension = 12; // half the window dimension

  std::vector<Feature>::iterator it;

  for(it=features.begin(); it != features.end(); it++) {
    int i, j, m, n;
    double min_ssd = std::numeric_limits<double>::max();
    int min_ssd_x = 0, min_ssd_y = 0;

    Feature ft = *it;

    // Loop through search space

    for(
      i = std::max(ft.centerX - searchSpaceXDim, windowDimension);
      i <= std::min(ft.centerX + searchSpaceXDim, this->Width() - windowDimension);
      i++
    ) {
      for(
        j = std::max(ft.centerY - searchSpaceYDim, windowDimension);
        j <= std::min(ft.centerY + searchSpaceYDim, this->Height() - windowDimension);
        j++
      ) {
        // For each pixel (i, j) in the search space
        double ssd = 0;

        // Calculate the SSD with the feature assuming (i, j) is the center of the new feature
        for(m=-1*windowDimension; m<=windowDimension; m++) {
          for(n=-1*windowDimension; n<=windowDimension; n++) {
            double oldLuminance = this->Pixel(ft.centerX + m, ft.centerY + n).Luminance();
            double newLuminance = otherImage->Pixel(i + m, j + n).Luminance();
            double diff = oldLuminance - newLuminance;
            ssd += diff * diff;
          }
        }

        // If the computed SSD is lower than the current minimum, set the current minimum to (i, j)
        if(ssd < min_ssd) {
          min_ssd = ssd;
          min_ssd_x = i;
          min_ssd_y = j;
        }
      }
    }

    // Line and box drawing
    for(m=-4; m<=4; m++) {
      for(n=-4; n<=4; n++) {
        output->Pixel(min_ssd_x+m, min_ssd_y+n).Reset(1, 1, 0, 1);
      }
    }

    int x1 = ft.centerX;
    int x2 = min_ssd_x;
    int y1 = ft.centerY;
    int y2 = min_ssd_y;
    int dx = x2 - x1;
    int dy = y2 - y1;

    if(dx != 0) { // avoid div by zero errors
      for(int x=std::min(x1, x2); x<=std::max(x1, x2); x++) {
        int y = int(std::round(y1 + (double(dy * (x-x1)) / double(dx))));
        output->Pixel(x, y).Reset(1, 1, 0, 1);
      }
    }
  }

  this->pixels = output->pixels;
  output->pixels = nullptr;
  delete output;
}