/*=========================================================================
 *
 *  Copyright NumFOCUS
 *
 *  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
 *
 *         https://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.
 *
 *=========================================================================*/

// Insight classes
#include "itkImageRegionIterator.h"
#include "itkMath.h"

#include "itkScalarImageToCooccurrenceMatrixFilter.h"

int
itkScalarImageToCooccurrenceMatrixFilterTest(int, char *[])
{

  // Data definitions
  constexpr unsigned int IMGWIDTH = 5;
  constexpr unsigned int IMGHEIGHT = 5;
  constexpr unsigned int NDIMENSION = 2;


  //------------------------------------------------------
  // Create a simple test images
  //------------------------------------------------------
  using InputImageType = itk::Image<unsigned char, NDIMENSION>;

  using InputImageIterator = itk::ImageRegionIterator<InputImageType>;


  auto image = InputImageType::New();

  InputImageType::SizeType inputImageSize = { { IMGWIDTH, IMGHEIGHT } };

  InputImageType::IndexType index;
  index.Fill(0);
  InputImageType::RegionType region;

  region.SetSize(inputImageSize);
  region.SetIndex(index);

  //--------------------------------------------------------------------------
  // Set up the image first. It looks like:
  //  1 2 1 2 1
  //  1 2 1 2 1
  //  1 2 1 2 1
  //  1 2 1 2 1
  //  1 2 1 2 1
  //--------------------------------------------------------------------------

  image->SetRegions(region);
  image->Allocate();

  // setup the iterator
  InputImageIterator imageIt(image, image->GetBufferedRegion());

  imageIt.GoToBegin();

  for (unsigned int i = 0; i < 5; ++i)
  {
    for (unsigned int j = 0; j < 5; j++, ++imageIt)
    {
      imageIt.Set(j % 2 + 1);
    }
  }


  //--------------------------------------------------------------------------
  // Generate the histogram. The un-normalized histogram should look like this:
  //
  //     0 1  2 ...
  //     ------
  //  0 |0 0  0
  //  1 |0 24 20
  //  2 |0 20 16
  //  3 |0 0  0
  //  .
  //  .
  // with zeroes elsewhere.
  //--------------------------------------------------------------------------

  bool passed = true;

  try
  {

    using FilterType = itk::Statistics::ScalarImageToCooccurrenceMatrixFilter<InputImageType>;

    auto filter = FilterType::New();

    // Invoke update before adding an input. An exception should be
    // thrown.
    try
    {
      filter->Update();
      passed = false;
      std::cerr << "Failed to throw expected exception due to nullptr input: " << std::endl;
      return EXIT_FAILURE;
    }
    catch (const itk::ExceptionObject & excp)
    {
      std::cout << "Expected exception caught: " << excp << std::endl;
    }

    filter->ResetPipeline();

    if (filter->GetInput() != nullptr)
    {
      std::cerr << "GetInput() should return nullptr since the input is not set yet " << std::endl;
      passed = false;
    }

    if (filter->GetMaskImage() != nullptr)
    {
      std::cerr << "GetMaskImage() should return nullptr since the mask image is not set yet " << std::endl;
      passed = false;
    }


    // Invoke update with a nullptr input. An exception should be
    // thrown.
    filter->SetInput(nullptr);
    try
    {
      filter->Update();
      passed = false;
      std::cerr << "Failed to throw expected exception due to nullptr input: " << std::endl;
      return EXIT_FAILURE;
    }
    catch (const itk::ExceptionObject & excp)
    {
      std::cout << "Expected exception caught: " << excp << std::endl;
    }

    filter->ResetPipeline();

    if (filter->GetInput() != nullptr)
    {
      passed = false;
    }


    filter->SetInput(image);

    InputImageType::OffsetType      offset1 = { { 0, 1 } };
    InputImageType::OffsetType      offset2 = { { 1, 0 } };
    FilterType::OffsetVectorPointer offsetV = FilterType::OffsetVector::New();
    offsetV->push_back(offset1);
    offsetV->push_back(offset2);

    filter->SetOffsets(offsetV);


    filter->Update();


    const FilterType::HistogramType * hist = filter->GetOutput();

    //--------------------------------------------------------------------------
    // Test the histogram.
    //--------------------------------------------------------------------------

    // First make sure the bins are sized properly:

    float max = hist->GetBinMax(0, 255);
    float min = hist->GetBinMin(0, 255);

    if (itk::Math::NotAlmostEquals(max, 256) || itk::Math::NotAlmostEquals(min, 255))
    {
      std::cerr << "Error" << std::endl;
      std::cerr << "The calculated bin sizes are incorrect" << std::endl;
      std::cerr << "Expected [255, 256), got [" << min << ", " << max << ')' << std::endl << std::endl;
      passed = false;
    }

    // Now make sure the contents of the bins are correct:
    using IndexType = FilterType::HistogramType::IndexType;
    IndexType one_one(hist->GetMeasurementVectorSize());
    IndexType one_two(hist->GetMeasurementVectorSize());
    IndexType two_one(hist->GetMeasurementVectorSize());
    IndexType two_two(hist->GetMeasurementVectorSize());

    one_one[0] = 1;
    one_one[1] = 1;

    one_two[0] = 1;
    one_two[1] = 2;

    two_one[0] = 2;
    two_one[1] = 1;

    two_two[0] = 2;
    two_two[1] = 2;

    float ooF, otF, toF, ttF, totalF;
    ooF = hist->GetFrequency(one_one);
    otF = hist->GetFrequency(one_two);
    toF = hist->GetFrequency(two_one);
    ttF = hist->GetFrequency(two_two);
    totalF = hist->GetTotalFrequency();

    if (itk::Math::NotAlmostEquals(ooF, 24.0f) || itk::Math::NotAlmostEquals(ttF, 16.0f) ||
        itk::Math::NotAlmostEquals(otF, 20.0f) || itk::Math::NotAlmostEquals(toF, 20.0f) ||
        itk::Math::NotAlmostEquals(ooF + ttF + otF + toF, totalF))
    {
      std::cerr << "Error:" << std::endl;
      std::cerr << "The histogram was calculated incorrectly" << std::endl;
      std::cerr << "Expected 24, 16, 20, 20, 80 got " << ooF << ", " << ttF << ", " << otF << ", " << toF << ", "
                << totalF << std::endl
                << std::endl;
      passed = false;
    }

    //--------------------------------------------------------------------------
    // Test the histogram with normalization
    //--------------------------------------------------------------------------


    auto filter0 = FilterType::New();

    filter0->SetInput(image);
    filter0->SetOffsets(offsetV);
    filter0->NormalizeOn();

    if (filter0->GetNormalize() != true)
    {
      std::cerr << "Normalize boolean is not set correctly";
      passed = false;
    }

    filter0->Update();
    const FilterType::HistogramType * hist0 = filter0->GetOutput();

    ooF = hist0->GetFrequency(one_one);
    otF = hist0->GetFrequency(one_two);
    toF = hist0->GetFrequency(two_one);
    ttF = hist0->GetFrequency(two_two);

    if ((ooF - 24 / 80.) > 0.001 || (ttF - 16 / 80.) > 0.001 || (otF - 20 / 80.) > 0.001 || (toF - 20 / 80.) > 0.001 ||
        (ooF + ttF + otF + toF - 1) > 0.001)
    {
      std::cerr << "Error:" << std::endl;
      std::cerr << "The histogram was calculated incorrectly" << std::endl;
      std::cerr << "Expected 0.3, 0.2, 0.25, 0.25 got " << ooF << ", " << ttF << ", " << otF << ", " << toF << std::endl
                << std::endl;
      passed = false;
    }


    //--------------------------------------------------------------------------
    // Generate some variant histograms and test them
    //--------------------------------------------------------------------------

    // First a histogram with 2 bins per axis
    auto filter2 = FilterType::New();

    filter2->SetInput(image);
    InputImageType::OffsetType offset3 = { { 0, 1 } };

    filter2->SetOffset(offset3);
    filter2->SetNumberOfBinsPerAxis(2);

    if (filter2->GetNumberOfBinsPerAxis() != 2)
    {
      std::cerr << "GetNumberOfBinsPerAxis() is not returning the expected value" << std::endl;
      passed = false;
    }

    filter2->Update();
    const FilterType::HistogramType * hist2 = filter2->GetOutput();

    IndexType zero_zero(hist2->GetMeasurementVectorSize());
    zero_zero[0] = 0;
    zero_zero[1] = 0;

    float zzF;
    zzF = hist2->GetFrequency(zero_zero);
    totalF = hist2->GetTotalFrequency();

    if (itk::Math::NotAlmostEquals(zzF, 40.0f) || itk::Math::NotAlmostEquals(zzF, totalF))
    {
      std::cerr << "Error:" << std::endl;
      std::cerr << "The degenerate histogram was calculated incorrectly" << std::endl;
      std::cerr << "Expected 40, 40 got " << zzF << ", " << totalF << std::endl << std::endl;
      passed = false;
    }


    // Next a histogram with a smaller range.
    auto filter3 = FilterType::New();

    filter3->SetInput(image);
    InputImageType::OffsetType offset4 = { { 1, 1 } };

    filter3->SetOffset(offset4);

    filter3->SetPixelValueMinMax(1, 2);

    if (filter3->GetMin() != 1)
    {
      std::cerr << "Error: " << std::endl;
      std::cerr << "GetMin() is not returning the value that is expected: 1" << std::endl;
      passed = false;
    }

    if (filter3->GetMax() != 2)
    {
      std::cerr << "Error: " << std::endl;
      std::cerr << "GetMin() is not returning the value that is expected: 2" << std::endl;
      passed = false;
    }


    filter3->SetNumberOfBinsPerAxis(2);

    filter3->Update();
    const FilterType::HistogramType * hist3 = filter3->GetOutput();

    IndexType zero_one(hist3->GetMeasurementVectorSize());
    IndexType one_zero(hist3->GetMeasurementVectorSize());

    zero_one[0] = 0;
    zero_one[1] = 1;

    one_zero[0] = 1;
    one_zero[1] = 0;

    float zoF, ozF;
    zzF = hist3->GetFrequency(zero_zero);
    zoF = hist3->GetFrequency(zero_one);
    ozF = hist3->GetFrequency(one_zero);
    ooF = hist3->GetFrequency(one_one);
    totalF = hist3->GetTotalFrequency();

    if (itk::Math::NotAlmostEquals(zzF, 0.0f) || itk::Math::NotAlmostEquals(zoF, 16.0f) ||
        itk::Math::NotAlmostEquals(ozF, 16.0f) || itk::Math::NotAlmostEquals(ooF, 0.0f) ||
        itk::Math::NotAlmostEquals(zzF + zoF + ozF + ooF, totalF))
    {
      std::cerr << "Error:" << std::endl;
      std::cerr << "The small size histogram was calculated incorrectly" << std::endl;
      std::cerr << "Expected 0, 16, 16, 0, 32 got " << zzF << ", " << zoF << ", " << ozF << ", " << ooF << ", "
                << totalF << std::endl
                << std::endl;
      passed = false;
    }

    // Next a histogram with a truncated range.
    auto filter4 = FilterType::New();

    filter4->SetInput(image);
    filter4->SetOffsets(offsetV);

    if (filter4->GetOffsets() != offsetV)
    {
      std::cerr << "Error: " << std::endl;
      std::cerr << "GetOffsets() is not returning the expected offset vector: " << std::endl;
      passed = false;
    }

    filter4->SetPixelValueMinMax(0, 1);
    filter4->SetNumberOfBinsPerAxis(2);

    filter4->Update();
    const FilterType::HistogramType * hist4 = filter4->GetOutput();

    zzF = hist4->GetFrequency(zero_zero);
    zoF = hist4->GetFrequency(zero_one);
    ozF = hist4->GetFrequency(one_zero);
    ooF = hist4->GetFrequency(one_one);
    totalF = hist4->GetTotalFrequency();

    if (itk::Math::NotAlmostEquals(zzF, 0.0f) || itk::Math::NotAlmostEquals(zoF, 0.0f) ||
        itk::Math::NotAlmostEquals(ozF, 0.0f) || itk::Math::NotAlmostEquals(ooF, 24.0f) ||
        itk::Math::NotAlmostEquals(zzF + zoF + ozF + ooF, totalF))
    {
      std::cerr << "Error:" << std::endl;
      std::cerr << "The truncated range histogram was calculated incorrectly" << std::endl;
      std::cerr << "Expected 0, 0, 0, 24, 24 got " << zzF << ", " << zoF << ", " << ozF << ", " << ooF << ", " << totalF
                << std::endl
                << std::endl;
      passed = false;
    }


    if (!passed)
    {
      std::cerr << "Test failed" << std::endl;
      return EXIT_FAILURE;
    }
    else
    {
      std::cerr << "Test succeeded" << std::endl;
      return EXIT_SUCCESS;
    }
  }
  catch (const itk::ExceptionObject & err)
  {
    std::cerr << "ExceptionObject caught !" << std::endl;
    std::cerr << err << std::endl;
    std::cerr << "Test failed" << std::endl;
    return EXIT_FAILURE;
  }
}
