<commit_before>/*=========================================================================
 *
 *  Copyright Insight Software 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.
 *
 *=========================================================================*/

#include "itkForwardFFTImageFilter.h"
#include "itkInverseFFTImageFilter.h"
#include "itkWaveletFrequencyForward.h"
#include "itkWaveletFrequencyInverse.h"
#include "itkWaveletFrequencyFilterBankGenerator.h"
#include "itkHeldIsotropicWavelet.h"
#include "itkVowIsotropicWavelet.h"
#include "itkSimoncelliIsotropicWavelet.h"
#include "itkShannonIsotropicWavelet.h"

#include "itkMonogenicSignalFrequencyImageFilter.h"
#include "itkVectorInverseFFTImageFilter.h"
#include "itkPhaseAnalysisSoftThresholdImageFilter.h"
#include "itkZeroDCImageFilter.h"

#include "itkImage.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"

#include "itkTestingMacros.h"

#include <string>

// Visualize for dev/debug purposes. Set in cmake file. Requires VTK
#ifdef ITK_VISUALIZE_TESTS
#include "itkViewImage.h"
#include "itkNumberToString.h"
#endif

// 1. Wavelet analysis (forward) on input image.
// 2. Create a Monogenic Signal (from Riesz function ) on each wavelet output..
// 3. Do a PhaseAnalysis on each Monogenic Signal.
// 4. Wavelet reconstruction (inverse) using as coefficients the output of the PhaseAnalysis.
// Without applying reconstruction factors: ApplyReconstructionFactorOff()
// 5. The result of the reconstruction will be an image that uses phase information at each level/band for improving local structure information, and can also work as an equalizator of brightness.
template< unsigned int VDimension, typename TWaveletFunction >
int
runRieszWaveletPhaseAnalysisTest( const std::string& inputImage,
  const std::string &, // outputImage
  const unsigned int& inputLevels,
  const unsigned int& inputBands,
  const bool applySoftThreshold)
{
  const unsigned int Dimension = VDimension;

  typedef double                             PixelType;
  typedef itk::Image< PixelType, Dimension > ImageType;
  typedef itk::ImageFileReader< ImageType >  ReaderType;

  typename ReaderType::Pointer reader = ReaderType::New();
  reader->SetFileName( inputImage );
  reader->Update();

  typedef itk::ZeroDCImageFilter< ImageType > ZeroDCFilterType;
  typename ZeroDCFilterType::Pointer zeroDCFilter = ZeroDCFilterType::New();
  zeroDCFilter->SetInput( reader->GetOutput() );
  zeroDCFilter->Update();

  // Perform FFT on input image.
  typedef itk::ForwardFFTImageFilter< typename ZeroDCFilterType::OutputImageType >
    FFTForwardFilterType;
  typename FFTForwardFilterType::Pointer fftForwardFilter = FFTForwardFilterType::New();
  fftForwardFilter->SetInput( zeroDCFilter->GetOutput() );
  fftForwardFilter->Update();
  typedef typename FFTForwardFilterType::OutputImageType ComplexImageType;

  typedef itk::InverseFFTImageFilter< ComplexImageType, ImageType > InverseFFTFilterType;

  // Forward Wavelet
  typedef TWaveletFunction WaveletFunctionType;
  typedef itk::WaveletFrequencyFilterBankGenerator< ComplexImageType, WaveletFunctionType >
    WaveletFilterBankType;
  typedef itk::WaveletFrequencyForward< ComplexImageType, ComplexImageType, WaveletFilterBankType >
    ForwardWaveletType;
  typename ForwardWaveletType::Pointer forwardWavelet = ForwardWaveletType::New();
  unsigned int highSubBands = inputBands;
  unsigned int levels = inputLevels;
  forwardWavelet->SetHighPassSubBands( highSubBands );
  forwardWavelet->SetLevels( levels );
  forwardWavelet->SetInput( fftForwardFilter->GetOutput() );
  forwardWavelet->Update();
  typename ForwardWaveletType::OutputsType analysisWavelets =
    forwardWavelet->GetOutputs();

  // Apply Monogenic signal to wavelet results
  typedef itk::MonogenicSignalFrequencyImageFilter< ComplexImageType >
    MonogenicSignalFrequencyFilterType;
  typedef typename MonogenicSignalFrequencyFilterType::OutputImageType
    VectorMonoOutputType;
  typedef itk::VectorInverseFFTImageFilter< VectorMonoOutputType >
    VectorInverseFFTType;
  typedef itk::PhaseAnalysisSoftThresholdImageFilter< typename VectorInverseFFTType::OutputImageType >
    PhaseAnalysisFilter;

  typename ForwardWaveletType::OutputsType modifiedWavelets;
  unsigned int numberOfOutputs = forwardWavelet->GetNumberOfOutputs();
  for ( unsigned int i = 0; i < forwardWavelet->GetNumberOfOutputs(); ++i )
    {
    std::cout << "Output #: " << i << " / " << numberOfOutputs - 1 << std::endl;
    typename MonogenicSignalFrequencyFilterType::Pointer monoFilter =
      MonogenicSignalFrequencyFilterType::New();
    typename VectorInverseFFTType::Pointer vecInverseFFT =
      VectorInverseFFTType::New();
    typename PhaseAnalysisFilter::Pointer phaseAnalyzer =
      PhaseAnalysisFilter::New();
    typename FFTForwardFilterType::Pointer fftForwardPhaseFilter =
      FFTForwardFilterType::New();

    // Generate a monogenic signal (vector valued)
    monoFilter->SetInput( analysisWavelets[i] );
    monoFilter->Update();

    vecInverseFFT->SetInput( monoFilter->GetOutput() );
    vecInverseFFT->Update();

    phaseAnalyzer->SetInput( vecInverseFFT->GetOutput() );
    phaseAnalyzer->SetApplySoftThreshold( applySoftThreshold );
    // if( i == numberOfOutputs - 1 ) // TODO Held does not modify approx image, but it does not generate better results.
    //   {
    //   phaseAnalyzer->SetApplySoftThreshold( false );
    //   }
    phaseAnalyzer->Update();

    fftForwardPhaseFilter->SetInput( phaseAnalyzer->GetOutputCosPhase() );
    fftForwardPhaseFilter->Update();

    modifiedWavelets.push_back( fftForwardPhaseFilter->GetOutput() );
    modifiedWavelets.back()->DisconnectPipeline();
    }

#ifdef ITK_VISUALIZE_TESTS
  // Visualize and compare modified wavelets coefficients (and approx image)
  bool visualizeCoefficients = false;
  if ( visualizeCoefficients )
    {
    for ( unsigned int i = 0; i < forwardWavelet->GetNumberOfOutputs(); ++i )
      {
      itk::NumberToString< unsigned int > n2s;
      typename InverseFFTFilterType::Pointer inverseFFT = InverseFFTFilterType::New();
      inverseFFT->SetInput(analysisWavelets[i]);
      inverseFFT->Update();
      itk::Testing::ViewImage( inverseFFT->GetOutput(), "WaveletCoef: output #" + n2s(i) );
      inverseFFT->SetInput(modifiedWavelets[i]);
      inverseFFT->Update();
      itk::Testing::ViewImage( inverseFFT->GetOutput(), "WaveletCoef. PhaseAnalyzed #" + n2s(i) );
      }
    }
#endif

  typedef itk::WaveletFrequencyInverse< ComplexImageType, ComplexImageType, WaveletFilterBankType >
    InverseWaveletType;
  typename InverseWaveletType::Pointer inverseWavelet = InverseWaveletType::New();
  inverseWavelet->SetHighPassSubBands( highSubBands );
  inverseWavelet->SetLevels( levels );
  inverseWavelet->SetInputs( modifiedWavelets );
  // The coefficients are now phases, do not apply reconstrucction factors.
  inverseWavelet->ApplyReconstructionFactorsOff();
  inverseWavelet->Update();

  typename InverseFFTFilterType::Pointer inverseFFT = InverseFFTFilterType::New();
  inverseFFT->SetInput( inverseWavelet->GetOutput() );
  inverseFFT->Update();

#ifdef ITK_VISUALIZE_TESTS
  itk::Testing::ViewImage( reader->GetOutput(), "Input Image" );
  itk::Testing::ViewImage( inverseFFT->GetOutput(), "Inverse Wavelet" );
#endif

  // typedef itk::ImageFileWriter< typename InverseFFTFilterType::OutputImageType > WriterType;
  // typename WriterType::Pointer writer = WriterType::New();
  // writer->SetFileName( outputImage );
  // writer->SetInput( inverseFFT->GetOutput() );
  //
  // TRY_EXPECT_NO_EXCEPTION( writer->Update() );
  //
  return EXIT_SUCCESS;
}

int
itkRieszWaveletPhaseAnalysisTest( int argc, char *argv[] )
{
  if ( argc < 7 || argc > 8 )
    {
    std::cerr << "Usage: " << argv[0]
              << " inputImage outputImage inputLevels inputBands waveletFunction Apply|NoApply [dimension]"
              << std::endl;
    return EXIT_FAILURE;
    }<commit_msg>TEST: Write output for RieszWaveletPhaseAnalysisTest.<commit_after>/*=========================================================================
 *
 *  Copyright Insight Software 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.
 *
 *=========================================================================*/

#include "itkForwardFFTImageFilter.h"
#include "itkInverseFFTImageFilter.h"
#include "itkWaveletFrequencyForward.h"
#include "itkWaveletFrequencyInverse.h"
#include "itkWaveletFrequencyFilterBankGenerator.h"
#include "itkHeldIsotropicWavelet.h"
#include "itkVowIsotropicWavelet.h"
#include "itkSimoncelliIsotropicWavelet.h"
#include "itkShannonIsotropicWavelet.h"

#include "itkMonogenicSignalFrequencyImageFilter.h"
#include "itkVectorInverseFFTImageFilter.h"
#include "itkPhaseAnalysisSoftThresholdImageFilter.h"
#include "itkZeroDCImageFilter.h"

#include "itkImage.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkCastImageFilter.h"
#include "itkNumberToString.h"

#include "itkTestingMacros.h"

#include <string>

// Visualize for dev/debug purposes. Set in cmake file. Requires VTK
#ifdef ITK_VISUALIZE_TESTS
#include "itkViewImage.h"
#endif

std::string
AppendToFilenameRiesz(const std::string& filename, const std::string & appendix)
{
  std::size_t foundDot = filename.find_last_of('.');
  return filename.substr( 0, foundDot ) + appendix + filename.substr( foundDot );
}

// 1. Wavelet analysis (forward) on input image.
// 2. Create a Monogenic Signal (from Riesz function ) on each wavelet output..
// 3. Do a PhaseAnalysis on each Monogenic Signal.
// 4. Wavelet reconstruction (inverse) using as coefficients the output of the PhaseAnalysis.
// Without applying reconstruction factors: ApplyReconstructionFactorOff()
// 5. The result of the reconstruction will be an image that uses phase information at each level/band for improving local structure information, and can also work as an equalizator of brightness.
template< unsigned int VDimension, typename TWaveletFunction >
int
runRieszWaveletPhaseAnalysisTest( const std::string& inputImage,
  const std::string & outputImage,
  const unsigned int& inputLevels,
  const unsigned int& inputBands,
  const bool applySoftThreshold)
{
  const unsigned int Dimension = VDimension;

  typedef double                             PixelType;
  typedef itk::Image< PixelType, Dimension > ImageType;
  typedef itk::ImageFileReader< ImageType >  ReaderType;

  itk::NumberToString< unsigned int > n2s;
  typename ReaderType::Pointer reader = ReaderType::New();
  reader->SetFileName( inputImage );
  reader->Update();

  typedef itk::ZeroDCImageFilter< ImageType > ZeroDCFilterType;
  typename ZeroDCFilterType::Pointer zeroDCFilter = ZeroDCFilterType::New();
  zeroDCFilter->SetInput( reader->GetOutput() );
  zeroDCFilter->Update();

  // Perform FFT on input image.
  typedef itk::ForwardFFTImageFilter< typename ZeroDCFilterType::OutputImageType >
    FFTForwardFilterType;
  typename FFTForwardFilterType::Pointer fftForwardFilter = FFTForwardFilterType::New();
  fftForwardFilter->SetInput( zeroDCFilter->GetOutput() );
  fftForwardFilter->Update();
  typedef typename FFTForwardFilterType::OutputImageType ComplexImageType;

  typedef itk::InverseFFTImageFilter< ComplexImageType, ImageType > InverseFFTFilterType;

  // Forward Wavelet
  typedef TWaveletFunction WaveletFunctionType;
  typedef itk::WaveletFrequencyFilterBankGenerator< ComplexImageType, WaveletFunctionType >
    WaveletFilterBankType;
  typedef itk::WaveletFrequencyForward< ComplexImageType, ComplexImageType, WaveletFilterBankType >
    ForwardWaveletType;
  typename ForwardWaveletType::Pointer forwardWavelet = ForwardWaveletType::New();
  unsigned int highSubBands = inputBands;
  unsigned int levels = inputLevels;
  forwardWavelet->SetHighPassSubBands( highSubBands );
  forwardWavelet->SetLevels( levels );
  forwardWavelet->SetInput( fftForwardFilter->GetOutput() );
  forwardWavelet->Update();
  typename ForwardWaveletType::OutputsType analysisWavelets =
    forwardWavelet->GetOutputs();

  // Apply Monogenic signal to wavelet results
  typedef itk::MonogenicSignalFrequencyImageFilter< ComplexImageType >
    MonogenicSignalFrequencyFilterType;
  typedef typename MonogenicSignalFrequencyFilterType::OutputImageType
    VectorMonoOutputType;
  typedef itk::VectorInverseFFTImageFilter< VectorMonoOutputType >
    VectorInverseFFTType;
  typedef itk::PhaseAnalysisSoftThresholdImageFilter< typename VectorInverseFFTType::OutputImageType >
    PhaseAnalysisFilter;

  typename ForwardWaveletType::OutputsType modifiedWavelets;
  unsigned int numberOfOutputs = forwardWavelet->GetNumberOfOutputs();
  for ( unsigned int i = 0; i < forwardWavelet->GetNumberOfOutputs(); ++i )
    {
    std::cout << "Output #: " << i << " / " << numberOfOutputs - 1 << std::endl;
    typename MonogenicSignalFrequencyFilterType::Pointer monoFilter =
      MonogenicSignalFrequencyFilterType::New();
    typename VectorInverseFFTType::Pointer vecInverseFFT =
      VectorInverseFFTType::New();
    typename PhaseAnalysisFilter::Pointer phaseAnalyzer =
      PhaseAnalysisFilter::New();
    typename FFTForwardFilterType::Pointer fftForwardPhaseFilter =
      FFTForwardFilterType::New();

    // Generate a monogenic signal (vector valued)
    monoFilter->SetInput( analysisWavelets[i] );
    monoFilter->Update();

    vecInverseFFT->SetInput( monoFilter->GetOutput() );
    vecInverseFFT->Update();

    phaseAnalyzer->SetInput( vecInverseFFT->GetOutput() );
    phaseAnalyzer->SetApplySoftThreshold( applySoftThreshold );
    // if( i == numberOfOutputs - 1 ) // TODO Held does not modify approx image, but it does not generate better results.
    //   {
    //   phaseAnalyzer->SetApplySoftThreshold( false );
    //   }
    phaseAnalyzer->Update();

    fftForwardPhaseFilter->SetInput( phaseAnalyzer->GetOutputCosPhase() );
    fftForwardPhaseFilter->Update();

    modifiedWavelets.push_back( fftForwardPhaseFilter->GetOutput() );
    modifiedWavelets.back()->DisconnectPipeline();
    }

#ifdef ITK_VISUALIZE_TESTS
  // Visualize and compare modified wavelets coefficients (and approx image)
  bool visualizeCoefficients = false;
  if ( visualizeCoefficients )
    {
    for ( unsigned int i = 0; i < forwardWavelet->GetNumberOfOutputs(); ++i )
      {
      typename InverseFFTFilterType::Pointer inverseFFT = InverseFFTFilterType::New();
      inverseFFT->SetInput(analysisWavelets[i]);
      inverseFFT->Update();
      itk::Testing::ViewImage( inverseFFT->GetOutput(), "WaveletCoef: output #" + n2s(i) );
      inverseFFT->SetInput(modifiedWavelets[i]);
      inverseFFT->Update();
      itk::Testing::ViewImage( inverseFFT->GetOutput(), "WaveletCoef. PhaseAnalyzed #" + n2s(i) );
      }
    }
#endif

  typedef itk::WaveletFrequencyInverse< ComplexImageType, ComplexImageType, WaveletFilterBankType >
    InverseWaveletType;
  typename InverseWaveletType::Pointer inverseWavelet = InverseWaveletType::New();
  inverseWavelet->SetHighPassSubBands( highSubBands );
  inverseWavelet->SetLevels( levels );
  inverseWavelet->SetInputs( modifiedWavelets );
  // The coefficients are now phases, do not apply reconstrucction factors.
  inverseWavelet->ApplyReconstructionFactorsOff();
  inverseWavelet->Update();

  typename InverseFFTFilterType::Pointer inverseFFT = InverseFFTFilterType::New();
  inverseFFT->SetInput( inverseWavelet->GetOutput() );
  inverseFFT->Update();

#ifdef ITK_VISUALIZE_TESTS
  itk::Testing::ViewImage( reader->GetOutput(), "Input Image" );
  itk::Testing::ViewImage( inverseFFT->GetOutput(), "Inverse Wavelet" );
#endif

  // Cast To Float for save as tiff.
  typedef itk::Image< float, Dimension > ImageFloatType;
  typedef itk::CastImageFilter< ImageType, ImageFloatType> CastFloatType;
  typename CastFloatType::Pointer caster = CastFloatType::New();
  caster->SetInput(inverseFFT->GetOutput());
  caster->Update();

  // typedef itk::ImageFileWriter< typename InverseFFTFilterType::OutputImageType > WriterType;
  typedef itk::ImageFileWriter< ImageFloatType > WriterType;
  typename WriterType::Pointer writer = WriterType::New();
  std::string appendString = "_" + n2s(inputLevels) + "_" + n2s(inputBands);
  std::string outputFile = AppendToFilenameRiesz(outputImage, appendString);
  writer->SetFileName( outputFile );
  writer->SetInput( caster->GetOutput() );

  TRY_EXPECT_NO_EXCEPTION( writer->Update() );
  //
  return EXIT_SUCCESS;
}

int
itkRieszWaveletPhaseAnalysisTest( int argc, char *argv[] )
{
  if ( argc < 7 || argc > 8 )
    {
    std::cerr << "Usage: " << argv[0]
              << " inputImage outputImage inputLevels inputBands waveletFunction Apply|NoApply [dimension]"
              << std::endl;
    return EXIT_FAILURE;
    }
{'fraction_non_alphanumeric': 0.11587114184516782, 'fraction_numerical': 0.0023050542531062013, 'mean_word_length': 5.226111305565278, 'pattern_counts': {'":': 0, '<': 78, '<?xml version=': 0, '>': 145, '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': '120632', 'n_tokens_mistral': 5372, 'n_tokens_neox': 5035, 'n_words': 1204}