Show code cell source

Hide code cell source

import warnings

warnings.filterwarnings(
    "ignore",
    message="torch.meshgrid: in an upcoming release,",
    category=UserWarning,
)

warnings.filterwarnings(
    "ignore",
    message="torch.range is deprecated and will be removed",
    category=UserWarning,
)

Run this notebook yourself!

Download the executed notebook: generate_windows.ipynb!

Generate Windows

Warning

This notebook requires the optional dependency plenoptic, which can be installed with pip.

This notebook provides tutorials on the most common ways of initializing and interacting with PoolingWindows, which constructs foveated windows and uses them to take weighted averages across an image.

import matplotlib.pyplot as plt
import matplotlib as mpl
import numpy as np
import torch

import fenestration as fen
import plenoptic as po

%load_ext autoreload
%autoreload 2
%matplotlib inline

Initializing PoolingWindows

Let’s begin by creating a PoolingWindows object for image size (256,256) and visualize the window contours that are created. Beyond the size of the image, the only required argument is scaling, which is defined as the ratio between a window’s radial full-width at half-maximum and its central eccentricity and thus controls the window size (see Choosing Scaling Values for more details).

We can also change a number of other optional parameters:

• min_ecc and max_ecc define the extent of the windows within the image to support conversion between pixels and degrees of visual angle. We do not place windows in the foveal region within the min_ecc ring.

• num_scales which controls the number of window scales generated.

• cache_dir for specifying a directory to cache the windows. If windows are already cached there, will load them instead of re-creating them. If they’re not present, will save them after creation.

• window_type which can be defined as "gaussian" or "cosine".

pw = fen.PoolingWindows(
  scaling=0.5,
  img_res=(256,256),
  min_ecc=0.5,
  max_ecc=15,
  window_type="gaussian"
)

Using PoolingWindows

To return the pooled averages of your input (an image, for example), use the forward method. Let’s say we are using a 2d (grayscale) image. PoolingWindows operates on 4d tensors (batch by channel by height by width; the convention for images in pytorch) or dictionaries of 4d tensors as input, so we will first have to unsqueeze the image until it is 4d.

img = torch.rand((1, 1, 256, 256), dtype=torch.float32)

Then we can call pw.forward() on the image! The output tensor(s) will be 3d with shape batch by channel by windows, where the batch and channel are handled independently.

Note

Note that pw(img) and pw.forward(img) are the same!

pooled = pw(img)
pooled.shape

torch.Size([1, 1, 1380])

This process is also equivalent to running window (subsetting the input into windows) followed by pool (computing a weighted average), but is more efficient without the intermediate products.

windowed = pw.window(img)
print(windowed.shape)
pooled = pw.pool(windowed)
print(pooled.shape)

torch.Size([1, 1, 1380, 256, 256])
torch.Size([1, 1, 1380])

Choosing Scaling Values

However, the scaling value used in previous examples was arbitrary. Let’s take a deeper dive into use cases of scaling, which is defined as the ratio of a window’s radial full-width at half-maximum to eccentricity.

As shown in Freeman, J., & Simoncelli, E. P. (2011). Metamers of the ventral stream. Nature Neuroscience., scaling values were adjusted to match specific behavioral thresholds for discriminating metameric stimuli. Additionally, in Broderick, W. F., Rufo, G., Winawer, J. & Simoncelli, E. P. (2023). Foveated metamers of the early visual system. eLife., the goal was to find the largest scaling value for generating metamers at which human and model discrimination performance was matched, and they showed how this critical value is impacted by image statistics, types of discrimination, and metamer synthesis initialization.

The fenestration.calculate.scaling function can be used to calculate a scaling value based on some desired properties. Say you are displaying images for an experiment and want to build pooling windows ranging from 1-10 degrees of eccentricity, tiling the radial space with either 5 or 10 windows. fenestration.calculate.scaling computes the scaling value needed to generate the corresponding PoolingWindows objects.

In the following plot, note how decreasing the scaling (while holding other arguments constant) leads to smaller windows and thus more angular wedges and eccentricity rings!

scaling_5win = fen.calculate.scaling(n_windows=5, min_ecc=1, max_ecc=10, std_dev=1)
scaling_10win = fen.calculate.scaling(n_windows=10, min_ecc=1, max_ecc=10, std_dev=1)
pw_5win = fen.PoolingWindows(scaling_5win, (256,256), min_ecc=1, max_ecc=10)
pw_10win = fen.PoolingWindows(scaling_10win, (256,256), min_ecc=1, max_ecc=10)

ax = pw_5win.plot_windows(subset=False);
ax.set_title(f"Scaling = {scaling_5win:.4f}");
ax = pw_10win.plot_windows(subset=False);
ax.set_title(f"Scaling = {scaling_10win:.4f}");

Visualizing PoolingWindows

PoolingWindows has a variety of helper functions that you can use to visualize its windows. You can also generate the eccentricity rings and angular wedges separately, for visualization or other purposes.

Visualizing Angles and Eccentricities

If you want to just generate the eccentricity rings and angular wedges separately, you can also call create_pooling_windows. Here we will use scaling=2 and and image size of (256,256). We will also take advantage of plenoptic’s plotting function po.plot.imshow.

angle_w, ecc_w = fen.create_pooling_windows(2, (256, 256))
# only show first 8 eccentricity rings
fig = po.plot.imshow(ecc_w[:8].unsqueeze(0))
fig = po.plot.imshow(angle_w.unsqueeze(0))

It is also possible to reconstruct the full windows from the separate angle and eccentricity tensors, though note this may take up a lot of memory. We can visualize each of the windows by indexing the first dimension of windows. If you plan to use this output, see Checking Windows for the important normalization step!

windows = torch.einsum('ahw,ehw->eahw', [angle_w, ecc_w]).flatten(0, 1)
win = windows[0,:,:]
while win.ndim < 4:
    win = win.unsqueeze(0)
po.plot.imshow(win, cmap="gray");

Visualizing Window Contours

We can also view the contours of the windows at their intersection points. For raised-cosine windows, this occurs at an amplitude of 0.5; for gaussian windows, this is at half a standard deviation away from the maximum. The default for plot_windows is to only plot four angle window slices to save time and memory, though you can plot all contours by passing subset=False.

pw.plot_windows();

Visualizing Window Values

Now let’s generate a figure with a noisy gradient across the image. We can then use plot_window_values to display the average values within each window.

img = torch.rand((1, 1, 256, 256), dtype=torch.float32) * torch.range(1/256,1,1/256)
po.plot.imshow(img, cmap="gray");

Since we are “pooling” the input within each window, we now see a smooth gradient across the windows returned after averaging out the noise.

pw = fen.PoolingWindows(0.8, (256,256))
pw.plot_window_values(img, subset=False);

Understanding PoolingWindows

In addition to the main functionality of creating and visualizing windows, we also have some tools that calculate window sizes and check for proper normalization.

Summarizing Window Sizes

We also have a few additional helper functions for understanding the windows, including plotting the window widths (left) and window areas (right). Both of these figures show the window sizes along the y axes as the eccentricity increases along the x axes, measured in degrees of visual angle (calculated based on min_ecc and max_ecc).

The window widths figure depicts two measurements: width of the windows along the radial (long) axis and angular (short) axis. Each individual window’s size is defined by three measurements: ‘top’, ‘half’, and ‘full’ widths. Top is the width of the flat-top region of each window where the window’s value is 1 (only present for cosine windows); full is the width of the entire window; half is the width at the half-max value. To get the approximate area, we multiply the radial width against the corresponding angular width, then divide by \(\frac{\pi}{4}\).

fig, ax = plt.subplots(1,2, figsize=(10, 4))
pw.plot_window_widths(ax=ax[0]);
pw.plot_window_areas(ax=ax[1]);

If you would like a summary of the size and values associated with the pooling windows, you can call summarize_window_sizes for either pixels or degrees. Here, we use pprint to help us cleanly display dictionaries of sizes.

from pprint import pprint
summary = pw.summarize_window_sizes(units="pixels")
pprint(summary)
summary = pw.summarize_window_sizes(units="degrees")
pprint(summary)

{'max_window_scale_0_area': np.float64(3449.4738989100847),
 'max_window_scale_0_center': np.float64(117.1538328521535),
 'max_window_scale_0_fwhm': np.float64(93.72306628172281),
 'min_window_scale_0_area': np.float64(8.874642083985288),
 'min_window_scale_0_center': np.float64(5.942312588205071),
 'min_window_scale_0_fwhm': np.float64(4.753850070564057)}
{'max_window_area': np.float64(47.37131514006159),
 'max_window_center': np.float64(13.728964787361738),
 'max_window_fwhm': np.float64(10.983171829889391),
 'min_window_area': np.float64(0.12187466240824522),
 'min_window_center': np.float64(0.6963647564302817),
 'min_window_fwhm': np.float64(0.5570918051442254)}

Checking Windows

These windows have been designed to have two important properties:

1. They sample the image in such a way that we can use interpolation to estimate any intermediate values without aliasing. See Sampling and Aliasing for more details.

2. They should be normalized so that each window has an L1-norm of 1 and thus contributes equally to the PoolingWindows output. This improves optimization performance during metamer generation and is explained more below.

We can use plot_window_checks to check whether the windows have been normalized properly so that they have an L1-norm of 1 to ensure that each eccentricity contributes equally. The first row shows the L1-norm of the windows, the second shows the sum. Each row will have one plot and, if everything worked correctly, they should each look like a sigmoid function that runs from 1 for small eccentricities to 0 for high eccentricities, measured in degrees of visual angle. Note that if you only use create_pooling_windows, it does not include this normalizaton step.

pw.plot_window_checks();