InsightSoftwareConsortium / ITKIsotropicWavelets

External Module for ITK, implementing Isotropic Wavelets and Riesz Filter for multiscale phase analysis.
Apache License 2.0
13 stars 11 forks source link

How to convert a “Swig Object of type 'itkImageCF3_Pointer *'” into a “itk.itkImagePython.itkImageCF3” #147

Closed mgcyung closed 3 years ago

mgcyung commented 3 years ago

I try to port the examples/runRieszWaveletAnalysis.cpp from C++ to python. And part of the code is as foloows

#********************************************************************************
# FFT padding
#********************************************************************************
import itk
import numpy as np
import nibabel as nib
import os

#********************************************************************************
# Types
#********************************************************************************
# parametres
Dimension = 3
highSubBands = 4
levels = 4
ImageType = itk.Image[itk.F, Dimension]
ComplexImageType = itk.Image[itk.complex[itk.F], Dimension]

#********************************************************************************
# image data
#********************************************************************************
# image volume
image_volume = np.random.rand(900, 900, 900).astype(np.float32)
input_volume = image_volume.astype(np.float32)
input_volume_itk = itk.GetImageFromArray(np.ascontiguousarray(input_volume))

#********************************************************************************
# FFT padding
#********************************************************************************
BoundaryConditionType = itk.ConstantBoundaryCondition[ImageType]
bounds = BoundaryConditionType()
bounds.SetConstant(itk.NumericTraits[itk.F].ZeroValue())

FFTPadFilterType = itk.FFTPadImageFilter[ImageType]
fftPadFilter = FFTPadFilterType.New()
fftPadFilter.SetInput(input_volume_itk)
fftPadFilter.SetBoundaryCondition(bounds)
fftPadFilter.Update()

input_volume_itk_padded = fftPadFilter.GetOutput()
# input_volume_itk_padding.GetLargestPossibleRegion().GetSize()

#********************************************************************************
# IsotropicWavelet
#********************************************************************************
# Generate data in frequency domain for IsotropicWavelet by FFT
FFTForwardFilterType = itk.ForwardFFTImageFilter[ImageType, ComplexImageType]
fftForwardFilter = FFTForwardFilterType.New()
fftForwardFilter.SetInput(input_volume_itk_padded)
fftForwardFilter.Update()
# InverseFFTFilterType = itk.InverseFFTImageFilter[ComplexImageType, ImageType]
output_volume_fft = fftForwardFilter.GetOutput()

# Construct ShannonIsotropicWavelet
WaveletType = itk.ShannonIsotropicWavelet[itk.F, Dimension,
                                          itk.Point[itk.D, Dimension]]
WaveletFilterBankType = itk.WaveletFrequencyFilterBankGenerator[ComplexImageType,
                                        WaveletType].New()
ForwardWaveletType = itk.WaveletFrequencyForward[ComplexImageType,
                                                 ComplexImageType, WaveletFilterBankType]
forwardWavelet = ForwardWaveletType.New()

# ShannonIsotropicWavelet transform
forwardWavelet.SetHighPassSubBands( highSubBands )
forwardWavelet.SetLevels( levels )
forwardWavelet.SetInput(output_volume_fft)
forwardWavelet.Update()
analysisWavelets = forwardWavelet.GetOutputs()

# Apply Monogenic signal to wavelet results
MonogenicSignalFrequencyFilterType = itk.MonogenicSignalFrequencyImageFilter[ComplexImageType]
VectorMonoOutputType = itk.VectorImage[itk.complex[itk.F], Dimension]
VectorInverseFFTType = itk.VectorInverseFFTImageFilter[VectorMonoOutputType, itk.VectorImage[itk.F, Dimension]]
PhaseAnalysisFilter = itk.PhaseAnalysisSoftThresholdImageFilter[itk.VectorImage[itk.F, Dimension], itk.Image[itk.F, Dimension]]

numberOfOutputs = forwardWavelet.GetNumberOfOutputs()
for i in range(numberOfOutputs):
    print('Output #: ' + str(i) + ' / ' + str(numberOfOutputs - 1))
    monoFilter = MonogenicSignalFrequencyFilterType.New()
    vecInverseFFT = VectorInverseFFTType.New()
    phaseAnalyzer = PhaseAnalysisFilter.New()
    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 );

There is a type mismatching in

monoFilter.SetInput(analysisWavelets[i])

In python2.7 the error message is

*** TypeError: in method 'itkImageToImageFilterICF3VICF3_SetInput', argument 2 of type 'itkImageCF3 const *'

And in python3.6 the error message is

*** TypeError: in method 'itkImageToImageFilterICF3VICF3_SetInput', argument 2 of type 'itkImageCF3 const *'
Additional information:
Wrong number or type of arguments for overloaded function 'itkImageToImageFilterICF3VICF3_SetInput'.
  Possible C/C++ prototypes are:
    itkImageToImageFilterICF3VICF3::SetInput(itkImageCF3 const *)
    itkImageToImageFilterICF3VICF3::SetInput(unsigned int,itkImageCF3 const *)

It seems “itkImageToImageFilterICF3VICF3::SetInput” needs an argument of type “itkImageCF3 const ”, but the “analysisWavelets[i]” is a “Swig Object of type 'itkImageCF3_Pointer '”. Is there a convertion from a “Swig Object of type 'itkImageCF3_Pointer '” to a “itkImageCF3 const ”, or other ways?

phcerdan commented 3 years ago

Hi @mgcyung, thanks a lot for the complete issue report.

This is related to the lack of python wrapping for std::vector<ImageTypePointer>, see #67. The workaround is easy, instead of using analysisWavelets[index] use forwardWavelet.GetOutput(index), which gives you a wrapped type.

Closing for now, feel free to open more issues with doubts/questions.

PS: An image of size 900, 900, 900 might be slow to analyze, I would recommend something like 128, 128, 128 for testing.

mgcyung commented 3 years ago

Thanks for the solution and suggestion.