Add image_fft and image_rfft

[tkoyama010]

Overview

Add vtkImageFFT and vtkImageRFFT wrapper.

Requires:

• implement complex dtype for pyvista_ndarray #2773
• Return correct scalar type for pyvista_ndarray #2770
• Do not add arrays when plotting #2763
• Feat/update bound actor #2441 (displaying frequencies without scaling the dataset)

---

Convert rgb to grayscale (see https://en.wikipedia.org/wiki/Grayscale#Luma_coding_in_video_systems).

r = image["JPEGImage"][:, 0]                                                                                                                                                                                
g = image["JPEGImage"][:, 1]                                                                                                                                                                                
b = image["JPEGImage"][:, 2]                                                                                                                                                                                
image.clear_data()                                                                                                                                                                                          
image["GrayScale"] = 0.299 * r + 0.587 * g + 0.114 * b                                                                                                                                                      
pv.global_theme.cmap = "gray"                                                                                                                                                                               
image.copy().plot(cpos="xy") 

It is also possible to apply filters to images. The following is the Fast Fourier Transformed image data.

fft = image.image_fft()                                                                                                                                                                                     
fft.copy().plot(cpos="xy", log_scale=True)                                                                                                                                                                  

Once Fast Fourier Transformed, images can also Reverse Fast Fourier Transformed.

rfft = fft.image_rfft()                                                                                                                                                                                     
rfft.copy().plot(cpos="xy")

Details

• None

[codecov]

Codecov Report

Merging #1953 (6f0e62f) into main (5f3da83) will increase coverage by 0.03%.
The diff coverage is 100.00%.

@@            Coverage Diff             @@
##             main    #1953      +/-   ##
==========================================
+ Coverage   94.03%   94.06%   +0.03%     
==========================================
  Files          76       76              
  Lines       16399    16495      +96     
==========================================
+ Hits        15420    15516      +96     
  Misses        979      979              

[tkoyama010]

WFM

[banesullivan]

I'd like to see if we can find a better example where we use the FFT to perform noise reduction or some other, more visible effect.

[banesullivan]

examples.download_gpr_data_array() would be a good dataset for this

[banesullivan]

I'd also be nice to show off a 3D example as well

[adeak]

I'd like to see if we can find a better example where we use the FFT to perform noise reduction or some other, more visible effect.

Does Perlin noise have structure in its Fourier spectrum corresponding to its base frequency? If so that could also be a fun demo.

[tkoyama010]

I'd like to see if we can find a better example where we use the FFT to perform noise reduction or some other, more visible effect.

I think we can learn from the OpenCV example: by adding vtkImageButterworthHighPass and vtkImageButterworthLowPass, we can add an example like the following.

[tkoyama010]

examples.download_gpr_data_array() would be a good dataset for this

Thanks. I will add it.

[tkoyama010]

Does Perlin noise have structure in its Fourier spectrum corresponding to its base frequency? If so that could also be a fun demo.

Sounds good! Thanks. I will add it too.

[RichardScottOZ]

I'd like to see if we can find a better example where we use the FFT to perform noise reduction or some other, more visible effect.

I think we can learn from the OpenCV example: by adding vtkImageButterworthHighPass and vtkImageButterworthLowPass, we can add an example like the following.

Built in high pass and low pass sounds good!

[tkoyama010]

Built in high pass and low pass sounds good!

That's right. I am trying to create an image of a geometrically simple shape in PyVista and create its outline.

[adeak]

Built in high pass and low pass sounds good!

That's right. I am trying to create an image of a geometrically simple shape in PyVista and create its outline.

Not sure if this helps, but perhaps a depth image of a simple scene might work.

[adeak]

I'm going to abandon what I tried for the Perlin noise FFT example, but here's what it is:

I really liked the gif that sweeps the cutoff of the frequency threshold for low_pass():

(side note: it's weird that this gif is called a "png", but nobody ever looks at this so it's fine.)

What this plot currently does is show an increasing amplitude, since we are discarding less and less from our input "signal". What I wanted to do is to scale up each frame so that we can emphasize the change in frequencies rather than a change in amplitude. But this won't look as nice as I'd hoped, because for low frequency cutoffs there's a sharp relative peak at the origin, which biases the mesh for the start of the animation (and preventing this would probably be too messy for an example snippet):

What I did was to change the helper function a bit to normalize the scalars:

def warp_low_pass_noise(cfreq):
    """Process the sampled FFT and warp by scalars."""
    output = sampled_fft.low_pass(cfreq, cfreq, cfreq).rfft()
    output['scalars'] = output.active_scalars.real
    # scale to fixed amplitude
    output['scalars'] = output['scalars'] * abs(sampled['scalars']).ptp() / abs(output['scalars']).ptp()
    return output.warp_by_scalar()

Like I said, I think this overall looks worse.

Side note 2: there's yet another weird thing going on if I do the scaling in-place. Initially I was getting hilariously huge numbers, think 1e+300 scale. But look at this small example:

>>> from pyvista import examples
>>> puppy = examples.download_puppy()
>>> puppy['JPEGImage'] = puppy['JPEGImage'].astype(float)  # uint8 -> float
>>> puppy['JPEGImage'] *= 2/3
Segmentation fault

It even happens with scalar data:

>>> import numpy as np
>>> from pyvista import examples
>>> puppy = examples.download_puppy()
>>> puppy.clear_data()
>>> puppy['scalars'] = np.arange(puppy.n_points, dtype=float)
>>> puppy['scalars'] *= 2/3
>>> puppy['scalars']
Segmentation fault

I have no idea what's going on here, but we really should figure this out before release... Assuming others can reproduce this and it's not my env being haunted. (Should be fixed in a different PR of course.)

To be clear this only happens with an in-place update:

>>> import numpy as np
>>> from pyvista import examples
>>> puppy = examples.download_puppy()
>>> puppy.clear_data()
>>> puppy['scalars'] = np.arange(puppy.n_points, dtype=float)
>>> puppy['scalars'] = puppy['scalars'] * 2 / 3
>>> puppy['scalars']
pyvista_ndarray([0.00000000e+00, 6.66666667e-01, 1.33333333e+00, ...,
                 1.27999800e+06, 1.27999867e+06, 1.27999933e+06])

I'm starting to think this is related to earlier mentions of VTK reference shenanigans. Data being garbage collected out from under a still-existing pyvista_ndarray, something like that...

[akaszynski]

Like I said, I think this overall looks worse.

I’d like to have both actually. Can you add it back?

[akaszynski]

I have no idea what's going on here, but we really should figure this out before release... Assuming others can reproduce this and it's not my env being haunted. (Should be fixed in a different PR of course.)

I'm unable to reproduce bug on this branch with

Reproduced the error locally. 100% agree this must be fixed before the next release.

[adeak]

Like I said, I think this overall looks worse.

I’d like to have both actually. Can you add it back?

I did what I could. I didn't want to repeat most of the last code block, so instead I tried to put both renders into the same figure. I'm sort of happy with the rendered result, but I'm not very happy with the code I wrote to get there; I think it should be cleaned up but I can't tell how. Improvements of any kind would be appreciated.

Here's how it roughly looks (give or take a document plot theme):

[adeak]

I'll want to look at the built docs for the example I edited, and someone will have to review my recent changes, but pending these this should be good to go.

[akaszynski]

I'm sort of happy with the rendered result

Rendered result looks excellent!

[tkoyama010]

Thank for improving this PR lots.

[adeak]

Thanks for the proofreading @tkoyama010 🚀 There are a few cases which I'd prefer to leave in the original state though, please see comments below. I'll leave it up to you to decide if you want to accept any of these suggestions, these are just my two cents as a fellow not-native-speaker :)

[adeak]

Thank you all for the patience, this is great 🎉

I didn't get an email about the latest comments and commit from @akaszynski, I hope this doesn't become a more general thing with github. I very much rely on email notifications...

Update: got one of at least 3 emails now, more than 2 hours late. I guess they'll arrive eventually.

[adeak]

In case anyone else is surprised by the dismissal of @akaszynski's review: it was auto-dismissed by this setting being enabled: #2245 (reply in thread)

[akaszynski]

In case anyone else is surprised by the dismissal of @akaszynski's review: it was auto-dismissed by this setting being enabled: #2245 (reply in thread)

I think this is a great idea to add; I'm quite guilty of getting an early approval and then adding whatever I see fit. This will at least help me in-check.