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

  Program:   ORFEO Toolbox
  Language:  C++
  Date:      $Date$
  Version:   $Revision$


  Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
  See OTBCopyright.txt for details.


     This software is distributed WITHOUT ANY WARRANTY; without even
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
     PURPOSE.  See the above copyright notices for more information.

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


#include "itkMacro.h"
#include <iostream>
#include <complex>  //only for the isnan() test line 148

#include "otbMath.h"
#include "otbVectorImage.h"
#include "otbBandMathXImageFilter.h"
#include "otbImageFileWriter.h"

#include "itkImageRegionIteratorWithIndex.h"

int otbBandMathXImageFilterNew( int itkNotUsed(argc), char* itkNotUsed(argv) [])
{
  typedef double                                            PixelType;
  typedef otb::VectorImage<PixelType, 2>                     ImageType;
  typedef otb::BandMathXImageFilter<ImageType>               FilterType;

   FilterType::Pointer         filter       = FilterType::New();

   return EXIT_SUCCESS;
}


int otbBandMathXImageFilter( int itkNotUsed(argc), char* itkNotUsed(argv) [])
{

  typedef otb::VectorImage<double, 2>              ImageType;
  typedef ImageType::PixelType                      PixelType;
  typedef otb::BandMathXImageFilter<ImageType>      FilterType;

  unsigned int i;
  const unsigned int N = 100, D1=3, D2=1, D3=1;

  ImageType::SizeType size;
  size.Fill(N);
  ImageType::IndexType index;
  index.Fill(0);
  ImageType::RegionType region;
  region.SetSize(size);
  region.SetIndex(index);

  ImageType::Pointer image1 = ImageType::New();
  ImageType::Pointer image2 = ImageType::New();
  ImageType::Pointer image3 = ImageType::New();

  image1->SetLargestPossibleRegion( region );
  image1->SetBufferedRegion( region );
  image1->SetRequestedRegion( region );
  image1->SetNumberOfComponentsPerPixel(D1);
  image1->Allocate();

  image2->SetLargestPossibleRegion( region );
  image2->SetBufferedRegion( region );
  image2->SetRequestedRegion( region );
  image2->SetNumberOfComponentsPerPixel(D2);
  image2->Allocate();

  image3->SetLargestPossibleRegion( region );
  image3->SetBufferedRegion( region );
  image3->SetRequestedRegion( region );
  image3->SetNumberOfComponentsPerPixel(D3);
  image3->Allocate();

  typedef itk::ImageRegionIteratorWithIndex<ImageType> IteratorType;
  IteratorType it1(image1, region);
  IteratorType it2(image2, region);
  IteratorType it3(image3, region);

  for (it1.GoToBegin(), it2.GoToBegin(), it3.GoToBegin(); !it1.IsAtEnd(); ++it1, ++it2, ++it3)
  {
    ImageType::IndexType i1 = it1.GetIndex();
    ImageType::IndexType i2 = it2.GetIndex();
    ImageType::IndexType i3 = it3.GetIndex();

    it1.Get()[0] = i1[0] + i1[1] -50; it1.Get()[1] = i1[0] * i1[1] -50; it1.Get()[2] = i1[0] / (i1[1]+1)+5;
    it2.Get()[0] = i2[0] * i2[1];
    it3.Get()[0] = i3[0] + i3[1] * i3[1];

  }


  FilterType::Pointer         filter       = FilterType::New();
  std::cout << "Number Of Threads  :  " << filter->GetNumberOfThreads() << std::endl;


  filter->SetNthInput(0, image1);
  filter->SetNthInput(1, image2);
  filter->SetNthInput(2, image3, "canal3");

  filter->SetExpression("vcos(2 * pi * im1) div (2 * pi * bands(im2,{1,1,1}) + {3.38,3.38,3.38}) mult vsin(pi * bands(canal3,{1,1,1}))"); //Sub-test 1
  filter->SetExpression("im1b1 / im2b1"); //Sub-test 2 (Edge Effect Handling)
  filter->Update();

  //if (filter->GetNumberOfOutputs() != 2)

  ImageType::Pointer output1 = filter->GetOutput(0);
  ImageType::Pointer output2 = filter->GetOutput(1);

  std::cout << "\n---  Standard Use\n";
  std::cout << "Parsed Expression :   " << filter->GetExpression(0) << std::endl;



  //Sub-test 1
  IteratorType itoutput1(output1, region);

  for (it1.GoToBegin(), it2.GoToBegin(), it3.GoToBegin(), itoutput1.GoToBegin(); !it1.IsAtEnd(); ++it1, ++it2, ++it3, ++itoutput1)
    {
    ImageType::IndexType i1 = it1.GetIndex();
    ImageType::IndexType i2 = it2.GetIndex();
    ImageType::IndexType i3 = it3.GetIndex();
    PixelType px1(D1),px2(D2),px3(D3);

    px1[0] = ( i1[0] + i1[1] -50 ); px1[1] = ( i1[0] * i1[1] -50 ); px1[2] = ( i1[0] / (i1[1]+1)+5 );
    px2[0] = ( i2[0] * i2[1] );
    px3[0] = ( i3[0] + i3[1] * i3[1] );

    double result1 = itoutput1.Get()[0], result2 = itoutput1.Get()[1], result3 = itoutput1.Get()[2];
    double error1,error2,error3;


    double expected1 = vcl_cos( 2 * otb::CONST_PI * px1[0] ) / ( 2 * otb::CONST_PI * px2[0] + 3.38 ) * vcl_sin( otb::CONST_PI * px3[0] );
    double expected2 = vcl_cos( 2 * otb::CONST_PI * px1[1] ) / ( 2 * otb::CONST_PI * px2[0] + 3.38 ) * vcl_sin( otb::CONST_PI * px3[0] );
    double expected3 = vcl_cos( 2 * otb::CONST_PI * px1[2] ) / ( 2 * otb::CONST_PI * px2[0] + 3.38 ) * vcl_sin( otb::CONST_PI * px3[0] );

    /*std::cout << "Pixel_1 =  " << it1.Get()[0] << "     Pixel_2 =  " << it2.Get()[0] << "     Pixel_3 =  " << it3.Get()[0]<commit_msg>BUG: proper use of iterators<commit_after>/*=========================================================================

  Program:   ORFEO Toolbox
  Language:  C++
  Date:      $Date$
  Version:   $Revision$


  Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
  See OTBCopyright.txt for details.


     This software is distributed WITHOUT ANY WARRANTY; without even
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
     PURPOSE.  See the above copyright notices for more information.

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


#include "itkMacro.h"
#include <iostream>
#include <complex>  //only for the isnan() test line 148

#include "otbMath.h"
#include "otbVectorImage.h"
#include "otbBandMathXImageFilter.h"
#include "otbImageFileWriter.h"

#include "itkImageRegionIteratorWithIndex.h"

int otbBandMathXImageFilterNew( int itkNotUsed(argc), char* itkNotUsed(argv) [])
{
  typedef double                                            PixelType;
  typedef otb::VectorImage<PixelType, 2>                     ImageType;
  typedef otb::BandMathXImageFilter<ImageType>               FilterType;

   FilterType::Pointer         filter       = FilterType::New();

   return EXIT_SUCCESS;
}


int otbBandMathXImageFilter( int itkNotUsed(argc), char* itkNotUsed(argv) [])
{

  typedef otb::VectorImage<double, 2>              ImageType;
  typedef ImageType::PixelType                      PixelType;
  typedef otb::BandMathXImageFilter<ImageType>      FilterType;

  unsigned int i;
  const unsigned int N = 100, D1=3, D2=1, D3=1;

  ImageType::SizeType size;
  size.Fill(N);
  ImageType::IndexType index;
  index.Fill(0);
  ImageType::RegionType region;
  region.SetSize(size);
  region.SetIndex(index);

  ImageType::Pointer image1 = ImageType::New();
  ImageType::Pointer image2 = ImageType::New();
  ImageType::Pointer image3 = ImageType::New();

  image1->SetLargestPossibleRegion( region );
  image1->SetBufferedRegion( region );
  image1->SetRequestedRegion( region );
  image1->SetNumberOfComponentsPerPixel(D1);
  image1->Allocate();

  image2->SetLargestPossibleRegion( region );
  image2->SetBufferedRegion( region );
  image2->SetRequestedRegion( region );
  image2->SetNumberOfComponentsPerPixel(D2);
  image2->Allocate();

  image3->SetLargestPossibleRegion( region );
  image3->SetBufferedRegion( region );
  image3->SetRequestedRegion( region );
  image3->SetNumberOfComponentsPerPixel(D3);
  image3->Allocate();

  typedef itk::ImageRegionIteratorWithIndex<ImageType> IteratorType;
  IteratorType it1(image1, region);
  IteratorType it2(image2, region);
  IteratorType it3(image3, region);

  ImageType::PixelType val1, val2, val3;
  val1.SetSize(D1);
  val2.SetSize(D2);
  val3.SetSize(D3);

  for (it1.GoToBegin(), it2.GoToBegin(), it3.GoToBegin(); !it1.IsAtEnd(); ++it1, ++it2, ++it3)
  {
    ImageType::IndexType i1 = it1.GetIndex();
    ImageType::IndexType i2 = it2.GetIndex();
    ImageType::IndexType i3 = it3.GetIndex();

    val1[0] = i1[0] + i1[1] -50;
    val1[1] = i1[0] * i1[1] -50;
    val1[2] = i1[0] / (i1[1]+1)+5;
    val2[0] = i2[0] * i2[1];
    val3[0] = i3[0] + i3[1] * i3[1];

    it1.Set(val1);
    it2.Set(val2);
    it3.Set(val3);

  }


  FilterType::Pointer         filter       = FilterType::New();
  std::cout << "Number Of Threads  :  " << filter->GetNumberOfThreads() << std::endl;


  filter->SetNthInput(0, image1);
  filter->SetNthInput(1, image2);
  filter->SetNthInput(2, image3, "canal3");

  filter->SetExpression("vcos(2 * pi * im1) div (2 * pi * bands(im2,{1,1,1}) + {3.38,3.38,3.38}) mult vsin(pi * bands(canal3,{1,1,1}))"); //Sub-test 1
  filter->SetExpression("im1b1 / im2b1"); //Sub-test 2 (Edge Effect Handling)
  filter->Update();

  //if (filter->GetNumberOfOutputs() != 2)

  ImageType::Pointer output1 = filter->GetOutput(0);
  ImageType::Pointer output2 = filter->GetOutput(1);

  std::cout << "\n---  Standard Use\n";
  std::cout << "Parsed Expression :   " << filter->GetExpression(0) << std::endl;

  //Sub-test 1
  IteratorType itoutput1(output1, region);

  for (it1.GoToBegin(), it2.GoToBegin(), it3.GoToBegin(), itoutput1.GoToBegin(); !it1.IsAtEnd(); ++it1, ++it2, ++it3, ++itoutput1)
    {
    ImageType::IndexType i1 = it1.GetIndex();
    ImageType::IndexType i2 = it2.GetIndex();
    ImageType::IndexType i3 = it3.GetIndex();
    PixelType px1(D1),px2(D2),px3(D3);

    px1[0] = ( i1[0] + i1[1] -50 ); px1[1] = ( i1[0] * i1[1] -50 ); px1[2] = ( i1[0] / (i1[1]+1)+5 );
    px2[0] = ( i2[0] * i2[1] );
    px3[0] = ( i3[0] + i3[1] * i3[1] );

    double result1 = itoutput1.Get()[0], result2 = itoutput1.Get()[1], result3 = itoutput1.Get()[2];
    double error1,error2,error3;


    double expected1 = vcl_cos( 2 * otb::CONST_PI * px1[0] ) / ( 2 * otb::CONST_PI * px2[0] + 3.38 ) * vcl_sin( otb::CONST_PI * px3[0] );
    double expected2 = vcl_cos( 2 * otb::CONST_PI * px1[1] ) / ( 2 * otb::CONST_PI * px2[0] + 3.38 ) * vcl_sin( otb::CONST_PI * px3[0] );
    double expected3 = vcl_cos( 2 * otb::CONST_PI * px1[2] ) / ( 2 * otb::CONST_PI * px2[0] + 3.38 ) * vcl_sin( otb::CONST_PI * px3[0] );

    /*std::cout << "Pixel_1 =  " << it1.Get()[0] << "     Pixel_2 =  " << it2.Get()[0] << "     Pixel_3 =  " << it3.Get()[0]
{'fraction_non_alphanumeric': 0.19123314065510597, 'fraction_numerical': 0.04980732177263969, 'mean_word_length': 3.84640522875817, 'pattern_counts': {'":': 0, '<': 69, '<?xml version=': 0, '>': 69, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '6742441', 'n_tokens_mistral': 3959, 'n_tokens_neox': 3798, 'n_words': 745}