<commit_before>/*=========================================================================
 *
 *  Copyright RTK Consortium
 *
 *  Licensed under the Apache License, Version 2.0 (the "License");
 *  you may not use this file except in compliance with the License.
 *  You may obtain a copy of the License at
 *
 *         http://www.apache.org/licenses/LICENSE-2.0.txt
 *
 *  Unless required by applicable law or agreed to in writing, software
 *  distributed under the License is distributed on an "AS IS" BASIS,
 *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 *  See the License for the specific language governing permissions and
 *  limitations under the License.
 *
 *=========================================================================*/

#ifndef __rtkBackwardDifferenceDivergenceImageFilter_hxx
#define __rtkBackwardDifferenceDivergenceImageFilter_hxx
#include "rtkBackwardDifferenceDivergenceImageFilter.h"

#include <itkConstShapedNeighborhoodIterator.h>
#include <itkNeighborhoodInnerProduct.h>
#include <itkImageRegionIterator.h>
#include <itkImageRegionConstIterator.h>
#include <itkNeighborhoodAlgorithm.h>
#include <itkConstantBoundaryCondition.h>
#include <itkOffset.h>
#include <itkProgressReporter.h>

namespace rtk
{

template <class TInputImage, class TOutputImage>
BackwardDifferenceDivergenceImageFilter<TInputImage, TOutputImage>
::BackwardDifferenceDivergenceImageFilter()
{
  m_UseImageSpacing = true;

  // default boundary condition
  m_BoundaryCondition = new itk::ConstantBoundaryCondition<TInputImage>();
  m_IsBoundaryConditionOverriden = false;

  // default behaviour is to process all dimensions
  for (int dim = 0; dim < TInputImage::ImageDimension; dim++)
    {
    m_DimensionsProcessed[dim] = true;
    }
}

template <class TInputImage, class TOutputImage>
BackwardDifferenceDivergenceImageFilter<TInputImage, TOutputImage>
::~BackwardDifferenceDivergenceImageFilter()
{
  delete m_BoundaryCondition;
}

// This should be handled by an itkMacro, but it doesn't seem to work with pointer types
template <class TInputImage, class TOutputImage>
void
BackwardDifferenceDivergenceImageFilter<TInputImage, TOutputImage>
::SetDimensionsProcessed(bool* DimensionsProcessed)
{
  bool Modified=false;
  for (int dim=0; dim<TInputImage::ImageDimension; dim++)
    {
    if (m_DimensionsProcessed[dim] != DimensionsProcessed[dim])
      {
      m_DimensionsProcessed[dim] = DimensionsProcessed[dim];
      Modified = true;
      }
    }
  if(Modified) this->Modified();
}

template <class TInputImage, class TOutputImage>
void
BackwardDifferenceDivergenceImageFilter<TInputImage, TOutputImage>
::OverrideBoundaryCondition(itk::ImageBoundaryCondition< TInputImage >* boundaryCondition)
{
  delete m_BoundaryCondition;
  m_BoundaryCondition = boundaryCondition;
  m_IsBoundaryConditionOverriden = true;
}

template <class TInputImage, class TOutputImage>
void
BackwardDifferenceDivergenceImageFilter<TInputImage, TOutputImage>
::GenerateInputRequestedRegion()
{
  // call the superclass' implementation of this method
  Superclass::GenerateInputRequestedRegion();

  // get pointers to the input and output
  typename TInputImage::Pointer inputPtr = const_cast< TInputImage * >( this->GetInput() );
  typename TOutputImage::Pointer outputPtr = this->GetOutput();

  if ( !inputPtr || !outputPtr )
    {
    return;
    }

  // get a copy of the input requested region (should equal the output
  // requested region)
  typename TInputImage::RegionType inputRequestedRegion;
  inputRequestedRegion = inputPtr->GetRequestedRegion();

  // pad the input requested region by the operator radius
  inputRequestedRegion.PadByRadius(1);

  // crop the input requested region at the input's largest possible region
  if ( inputRequestedRegion.Crop( inputPtr->GetLargestPossibleRegion() ) )
    {
    inputPtr->SetRequestedRegion(inputRequestedRegion);
    return;
    }
  else
    {
    // Couldn't crop the region (requested region is outside the largest
    // possible region).  Throw an exception.

    // store what we tried to request (prior to trying to crop)
    inputPtr->SetRequestedRegion(inputRequestedRegion);

    // build an exception
    itk::InvalidRequestedRegionError e(__FILE__, __LINE__);
    e.SetLocation(ITK_LOCATION);
    e.SetDescription("Requested region is (at least partially) outside the largest possible region.");
    e.SetDataObject(inputPtr);
    throw e;
    }
}

template <class TInputImage, class TOutputImage>
void
BackwardDifferenceDivergenceImageFilter< TInputImage, TOutputImage>
::BeforeThreadedGenerateData()
{
  m_SpacingCoeffs= this->GetInput()->GetSpacing();
  if (m_UseImageSpacing == false) m_SpacingCoeffs.Fill(1);
}

template <class TInputImage, class TOutputImage>
void
BackwardDifferenceDivergenceImageFilter< TInputImage, TOutputImage>
::ThreadedGenerateData(const typename TInputImage::RegionType& outputRegionForThread, itk::ThreadIdType itkNotUsed(threadId))
{
  // Generate a list of indices of the dimensions to process
  std::vector<int> dimsToProcess;
  for (int dim = 0; dim < TInputImage::ImageDimension; dim++)
    {
    if(m_DimensionsProcessed[dim]) dimsToProcess.push_back(dim);
    }

  typename TOutputImage::Pointer output = this->GetOutput();
  typename TInputImage::ConstPointer input = this->GetInput();

  itk::ConstantBoundaryCondition<InputImageType> cbc;

  itk::ImageRegionIterator<TOutputImage> oit(output, outputRegionForThread);
  oit.GoToBegin();

  itk::Size<InputImageDimension> radius;
  radius.Fill(1);

  itk::ConstNeighborhoodIterator<TInputImage> iit(radius, input, outputRegionForThread);
  iit.GoToBegin();
  iit.OverrideBoundaryCondition(&cbc);

  itk::SizeValueType c = (itk::SizeValueType) (iit.Size() / 2); // get offset of center pixel
  itk::SizeValueType strides[TOutputImage::ImageDimension]; // get offsets to access neighboring pixels
  for (int dim=0; dim<TOutputImage::ImageDimension; dim++)
    {
    strides[dim] = iit.GetStride(dim);
    }

  float div;
  while(!oit.IsAtEnd())
    {
    div = 0;
    // Compute the local differences around the central pixel
    for (unsigned int k = 0; k < dimsToProcess.size(); k++)
      {
      div += (iit.GetPixel(c)[k] - iit.GetPixel(c - strides[dimsToProcess[k]])[k]) / m_SpacingCoeffs[dimsToProcess[k]];<commit_msg>BUG: Requested BoundaryCondition was not respected<commit_after>/*=========================================================================
 *
 *  Copyright RTK Consortium
 *
 *  Licensed under the Apache License, Version 2.0 (the "License");
 *  you may not use this file except in compliance with the License.
 *  You may obtain a copy of the License at
 *
 *         http://www.apache.org/licenses/LICENSE-2.0.txt
 *
 *  Unless required by applicable law or agreed to in writing, software
 *  distributed under the License is distributed on an "AS IS" BASIS,
 *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 *  See the License for the specific language governing permissions and
 *  limitations under the License.
 *
 *=========================================================================*/

#ifndef __rtkBackwardDifferenceDivergenceImageFilter_hxx
#define __rtkBackwardDifferenceDivergenceImageFilter_hxx
#include "rtkBackwardDifferenceDivergenceImageFilter.h"

#include <itkConstShapedNeighborhoodIterator.h>
#include <itkNeighborhoodInnerProduct.h>
#include <itkImageRegionIterator.h>
#include <itkImageRegionConstIterator.h>
#include <itkNeighborhoodAlgorithm.h>
#include <itkConstantBoundaryCondition.h>
#include <itkOffset.h>
#include <itkProgressReporter.h>

namespace rtk
{

template <class TInputImage, class TOutputImage>
BackwardDifferenceDivergenceImageFilter<TInputImage, TOutputImage>
::BackwardDifferenceDivergenceImageFilter()
{
  m_UseImageSpacing = true;

  // default boundary condition
  m_BoundaryCondition = new itk::ConstantBoundaryCondition<TInputImage>();
  m_IsBoundaryConditionOverriden = false;

  // default behaviour is to process all dimensions
  for (int dim = 0; dim < TInputImage::ImageDimension; dim++)
    {
    m_DimensionsProcessed[dim] = true;
    }
}

template <class TInputImage, class TOutputImage>
BackwardDifferenceDivergenceImageFilter<TInputImage, TOutputImage>
::~BackwardDifferenceDivergenceImageFilter()
{
  delete m_BoundaryCondition;
}

// This should be handled by an itkMacro, but it doesn't seem to work with pointer types
template <class TInputImage, class TOutputImage>
void
BackwardDifferenceDivergenceImageFilter<TInputImage, TOutputImage>
::SetDimensionsProcessed(bool* DimensionsProcessed)
{
  bool Modified=false;
  for (int dim=0; dim<TInputImage::ImageDimension; dim++)
    {
    if (m_DimensionsProcessed[dim] != DimensionsProcessed[dim])
      {
      m_DimensionsProcessed[dim] = DimensionsProcessed[dim];
      Modified = true;
      }
    }
  if(Modified) this->Modified();
}

template <class TInputImage, class TOutputImage>
void
BackwardDifferenceDivergenceImageFilter<TInputImage, TOutputImage>
::OverrideBoundaryCondition(itk::ImageBoundaryCondition< TInputImage >* boundaryCondition)
{
  delete m_BoundaryCondition;
  m_BoundaryCondition = boundaryCondition;
  m_IsBoundaryConditionOverriden = true;
}

template <class TInputImage, class TOutputImage>
void
BackwardDifferenceDivergenceImageFilter<TInputImage, TOutputImage>
::GenerateInputRequestedRegion()
{
  // call the superclass' implementation of this method
  Superclass::GenerateInputRequestedRegion();

  // get pointers to the input and output
  typename TInputImage::Pointer inputPtr = const_cast< TInputImage * >( this->GetInput() );
  typename TOutputImage::Pointer outputPtr = this->GetOutput();

  if ( !inputPtr || !outputPtr )
    {
    return;
    }

  // get a copy of the input requested region (should equal the output
  // requested region)
  typename TInputImage::RegionType inputRequestedRegion;
  inputRequestedRegion = inputPtr->GetRequestedRegion();

  // pad the input requested region by the operator radius
  inputRequestedRegion.PadByRadius(1);

  // crop the input requested region at the input's largest possible region
  if ( inputRequestedRegion.Crop( inputPtr->GetLargestPossibleRegion() ) )
    {
    inputPtr->SetRequestedRegion(inputRequestedRegion);
    return;
    }
  else
    {
    // Couldn't crop the region (requested region is outside the largest
    // possible region).  Throw an exception.

    // store what we tried to request (prior to trying to crop)
    inputPtr->SetRequestedRegion(inputRequestedRegion);

    // build an exception
    itk::InvalidRequestedRegionError e(__FILE__, __LINE__);
    e.SetLocation(ITK_LOCATION);
    e.SetDescription("Requested region is (at least partially) outside the largest possible region.");
    e.SetDataObject(inputPtr);
    throw e;
    }
}

template <class TInputImage, class TOutputImage>
void
BackwardDifferenceDivergenceImageFilter< TInputImage, TOutputImage>
::BeforeThreadedGenerateData()
{
  m_SpacingCoeffs= this->GetInput()->GetSpacing();
  if (m_UseImageSpacing == false) m_SpacingCoeffs.Fill(1);
}

template <class TInputImage, class TOutputImage>
void
BackwardDifferenceDivergenceImageFilter< TInputImage, TOutputImage>
::ThreadedGenerateData(const typename TInputImage::RegionType& outputRegionForThread, itk::ThreadIdType itkNotUsed(threadId))
{
  // Generate a list of indices of the dimensions to process
  std::vector<int> dimsToProcess;
  for (int dim = 0; dim < TInputImage::ImageDimension; dim++)
    {
    if(m_DimensionsProcessed[dim]) dimsToProcess.push_back(dim);
    }

  typename TOutputImage::Pointer output = this->GetOutput();
  typename TInputImage::ConstPointer input = this->GetInput();

  itk::ImageRegionIterator<TOutputImage> oit(output, outputRegionForThread);
  oit.GoToBegin();

  itk::Size<InputImageDimension> radius;
  radius.Fill(1);

  itk::ConstNeighborhoodIterator<TInputImage> iit(radius, input, outputRegionForThread);
  iit.GoToBegin();
  iit.OverrideBoundaryCondition(m_BoundaryCondition);

  itk::SizeValueType c = (itk::SizeValueType) (iit.Size() / 2); // get offset of center pixel
  itk::SizeValueType strides[TOutputImage::ImageDimension]; // get offsets to access neighboring pixels
  for (int dim=0; dim<TOutputImage::ImageDimension; dim++)
    {
    strides[dim] = iit.GetStride(dim);
    }

  float div;
  while(!oit.IsAtEnd())
    {
    div = 0;
    // Compute the local differences around the central pixel
    for (unsigned int k = 0; k < dimsToProcess.size(); k++)
      {
      div += (iit.GetPixel(c)[k] - iit.GetPixel(c - strides[dimsToProcess[k]])[k]) / m_SpacingCoeffs[dimsToProcess[k]];
{'fraction_non_alphanumeric': 0.11657025120849512, 'fraction_numerical': 0.002218876297646406, 'mean_word_length': 5.338523355097941, 'pattern_counts': {'":': 0, '<': 72, '<?xml version=': 0, '>': 84, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '2362691', 'n_tokens_mistral': 3780, 'n_tokens_neox': 3693, 'n_words': 936}