Fourier Transform Textural Ordination in Python

1. Introduction

The fototex package implements the FOTO algorithm described in Textural ordination based on Fourier spectral decomposition: a method to analyze and compare landscape patterns (Pierre Couteron, Nicolas Barbier and Denis Gautier, 2006), which allows texture characterization and comparison. Though this algorithm has already been implemented in MATLAB and R, fototex enriches this collection by both implementing new features and dramatically enhancing performance.

2. Basic usage

The main module of fototex is foto, in which various classes might be used and combined to structure your own FOTO algorithm, suited to your own needs. Here is the list of the main available classes at the moment:

class FotoBase
class Batch(FotoBase)
class Sector(FotoBase)
class Foto(FotoBase)
class FotoBatch(Batch)
class FotoSector(Foto, Sector)
class FotoSectorBatch(FotoBatch, Sector)
  • FotoBase: base class describing the core of the FOTO algorithm
    • Batch: base subclass implementing image batch supply
    • Sector: base subclass implementing radial spectra anisotropy calculation
    • Foto: subclass implementing the typical FOTO algorithm, i.e. the isotropic FOTO algorithm on one given image

By using inheritance and multiple inheritance, you can then combine the different features between them. Typically, at the moment, those subclasses are the following:

  • FotoBatch: inherits from Batch and allows the FOTO algorithm to be implemented from multiple image batches
  • FotoSector: inherits from both Foto and Sector classes, and implements anisotropic FOTO
  • FotoSectorBatch: inherits from both FotoBatch and Sectorclasses, and implements anisotropic FOTO from multiple image batches

2.1 Main features

2.1.1 Sliding window method

The FotoBase gives you 2 methods for sliding window over an image or image batches. The first is the block method (method="block") that computes r-spectra for each window in the image block by block. The second is the moving method (method="moving_window") that computes r-spectra for a window sliding from West to East and North to South. This method requires far more CPU resources. In that case, the total number of windows will be the size of the input image. For example, an image that would require 150k block 11x11 windows would require about 20M windows in moving mode.

2.1.2 Isotropic vs. anisotropic R-spectra

Isotropic r-spectra does not look at the direction but rather averages the spectra over the entire azimuth (i.e. from 0° to 360°). Conversely, anisotropic calculation takes care of the direction, that is the angular sector for which the spectra must be averaged. In the main classes of the fototex package, everything that is not Sector based is isotropic, and everything that is Sector based is anisotropic. In the anisotropic case, the maximum number of sectors should somehow be given by the square pixel-based geometry of a window: 8 as the number of pixels surrounding the window center. However, it is possible to increase this number by considering nearest neighboring for specific radius values. In other words, it is not possible to get 12 bin sectors for the first 8 surrouding pixels, but the 4 "in-between" sectors that are not digitally represented might be assessed from the nearest available values.

2.1.3 Single image against set of image batches

The basic FOTO algorithm takes one single image as input. In order to take into account different textures in the training process from different places and images, you can use the Batch subclass. Everything that is Batch implements FOTO algorithms from multiple image batches, and everything without implements the normal algorithm.

2.1.4 In memory against H5 hard storage

Whether your machine has a lot of memory of the image size is not too large, the FotoBase subclasses and processes might be run in memory (in_memory=True). Conversely, if for any reason you are out of memory, there is still the possibility to run the algorithm using H5 hard temporary storage. In this case, be warned that io reading and writing will increase the CPU time. You may adjust the size of the chunks to be read from and written to the H5 temporary file by setting the data_chunk_size argument.

2.1.5 Multiprocessing

Everything in the fototex package is parallelized (except io reading/writing). The number of processes is available through nb_processes in the different FotoBase class and subclass methods. By default, this number is set to the maximum available number of CPU on your machine (nb_processes=mp.cpu_count()). In case you do not want to blow your computer off, you may want to reduce this number.

2.4 The FotoBase class and subclasses

2.4.1 The Foto class

This class implements the main FOTO algorithm, i.e. computes isotropic r-spectra and applies dimensionality reduction through Principal Component Analysis.

Constructor:

Foto(image, band=None, method="block", in_memory=True, data_chunk_size=50000)

Required parameters:

  • image: (str) path to a valid GeoTIFF raster/image

Optional parameters:

  • band: (int) band number if multiband raster
  • method: (str) valid window method among {"block", "moving_window"}
  • in_memory: (bool) either implement FOTO algorithm (r-spectra and PCA) in memory or use HDF5 file storage on the fly
  • data_chunk_size: (int) if HDF5 storage is used, the number of data per chunk for reading/writing to file

2.4.2 The FotoSector class

The FotoSector class implements the calculation of anisotropic r-spectra, that is averaged r-spectra for given angular sectors. It is possible to give the desired orientation to those r-spectra, by setting the direction of the starting sector. By default, it is set to 0°, i.e. North (start_sector=0).

Constructor:

FotoSector(image, nb_sectors=6, start_sector=0, band=None, method="block", in_memory=True, 
           data_chunk_size=50000)

Required parameters:

  • image: (str) path to a valid GeoTIFF raster/image

Optional parameters:

  • nb_sectors: (int) number of sectors for which r-spectra must be computed (default: 6, max: Sector.max_nb_sectors = fototex.MAX_NB_SECTORS)
  • start_sector: (float, int) starting sector direction (default: 0°, i.e. North)
  • band: (int) band number if multiband raster
  • method: (str) valid window method among {"block", "moving_window"}
  • in_memory: (bool) either implement FOTO algorithm (r-spectra and PCA) in memory or use HDF5 file storage on the fly
  • data_chunk_size: (int) if HDF5 storage is used, the number of data per chunk for reading/writing to file

2.4.3 The FotoBatch class

2.5 The run method

The run method typically consists in computing r-spectra for each window in either the image or set of image batches and then applying dimensionality reduction through a Principal Component Analysis (PCA).

Usage:

Foto.run(window_size, nb_sampled_frequencies=None, standardize=False, keep_dc_component=False,
         at_random=False, batch_size=None, max_iter=1000, nb_processes=mp.cpu_count())

Required parameters:

  • window_size: (int) size of the sliding window that will be passed throughout the image

Optional parameters:

  • nb_sampled_frequencies: (int) number of frequencies sampled within window. If None, is inferred
  • standardize: (bool) if True, standardize window power spectrum density by window variance
  • keep_dc_component: (bool) if True, keep DC component (frequency 0) when computing r-spectra
  • at_random: (bool) when HDF5 storage is active (in_memory == False), if True use random batches for incremental PCA, otherwise use the whole dataset
  • batch_size: (int) when using incremental PCA at random, define the size of the batch
  • max_iter: (int) when using incremental PCA at random, maximum number of iterations before stopping the process
  • nb_processes: (int) number of processes for multiprocessing (by default, all CPU are used)

2.6 Examples

Import everything we need

In [11]:
import os
from fototex.foto import Foto, FotoSector, FotoBatch, FotoSectorBatch
from fototex.foto_tools import degrees_to_cardinal
import numpy as np
from matplotlib import pyplot as plt
from skimage import io

Let's read one small sample image to see what we can do with it:

In [57]:
image_file = os.path.join("tutorial", "pan_image_test.tif")
sample_image = io.imread(image_file)
In [19]:
plt.rcParams['figure.figsize'] = [15, 10]
img = plt.imshow(sample_image, cmap='gray', vmin=50, vmax=106)
plt.axis('off')
Out[19]:
(-0.5, 2240.5, 1665.5, -0.5)

2.6.1 Run the isotropic FOTO algorithm on a given image

Build a simple Foto instance, using the block method and applying the process in memory:

In [6]:
foto_obj = Foto(image_file, method="block", in_memory=True)

Now we can run it with a window size of 11:

In [7]:
foto_obj.run(window_size=11, keep_dc_component=True)
Retrieve isotropic R-spectra: 100%|██████████| 31.0k/31.0k [00:05<00:00, 5.91kit/s]  
Run Principal Component Analysis: 100%|██████████| 1.53k/1.53k [00:00<00:00, 15.9kit/s]

Now we may save the reduced r-spectra as a RGB map:

In [31]:
foto_obj.save_rgb()
Write RGB output image to 'tutorial/pan_image_test_method=block_wsize=11_rgb.tif': 100%|██████████| 152/152 [00:00<00:00, 6.91kit/s]

Now we can open this image, read it and look at what it looks like:

In [53]:
rgb_image = io.imread(foto_obj.rgb_file)
#rgb_image[rgb_image==-999] = np.nan

# 2%-98% range for values in each channel 
value_min = [-1.42, -1.04, -1.29]
value_max = [6.28, 1.56, 1.22]

for vmin, vmax, band in zip(value_min, value_max, range(rgb_image.shape[2])):
    rgb_image[:,:,band] = (rgb_image[:,:,band] - vmin)/ (vmax - vmin)
In [54]:
plt.imshow(rgb_image)
plt.axis('off')
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).
Out[54]:
(-0.5, 203.5, 151.5, -0.5)

2.6.2 Run the anisotropic FOTO algorithm

Build a FotoSector instance, considering 8 sectors starting from 0° (North) and the block method:

In [58]:
foto_obj = FotoSector(image_file, nb_sectors=8, start_sector=0, method="block", in_memory=True)

Now, we may run it for the same window size as before:

In [59]:
foto_obj.run(window_size=11, keep_dc_component=True)
Retrieve anisotropic R-spectra: 100%|██████████| 31.0k/31.0k [00:07<00:00, 4.21kit/s]  
Run Principal Component Analysis: 100%|██████████| 98.1k/98.1k [00:00<00:00, 471kit/s] 

And save the various RGB maps (one per sector):

In [60]:
foto_obj.save_rgb()
Write RGB output image to 'tutorial/pan_image_test_method=block_wsize=11_sector=225_SW_rgb.tif': 100%|██████████| 152/152 [00:00<00:00, 7.64kit/s]
Write RGB output image to 'tutorial/pan_image_test_method=block_wsize=11_sector=270_W_rgb.tif': 100%|██████████| 152/152 [00:00<00:00, 5.73kit/s]
Write RGB output image to 'tutorial/pan_image_test_method=block_wsize=11_sector=315_NW_rgb.tif': 100%|██████████| 152/152 [00:00<00:00, 8.33kit/s]
Write RGB output image to 'tutorial/pan_image_test_method=block_wsize=11_sector=0_N_rgb.tif': 100%|██████████| 152/152 [00:00<00:00, 8.79kit/s]
Write RGB output image to 'tutorial/pan_image_test_method=block_wsize=11_sector=45_NE_rgb.tif': 100%|██████████| 152/152 [00:00<00:00, 8.92kit/s]
Write RGB output image to 'tutorial/pan_image_test_method=block_wsize=11_sector=90_E_rgb.tif': 100%|██████████| 152/152 [00:00<00:00, 7.89kit/s]
Write RGB output image to 'tutorial/pan_image_test_method=block_wsize=11_sector=135_SE_rgb.tif': 100%|██████████| 152/152 [00:00<00:00, 8.85kit/s]
Write RGB output image to 'tutorial/pan_image_test_method=block_wsize=11_sector=180_S_rgb.tif': 100%|██████████| 152/152 [00:00<00:00, 11.8kit/s]

Now, we can plot the 8 RGB images corresponding to as much sectors:

In [67]:
plt.rcParams['figure.figsize'] = [20, 30]
fig, axes = plt.subplots(4,2)
sector = 0
for ax, file in zip(axes.ravel(), foto_obj.rgb_file):
    img = io.imread(file)
    for vmin, vmax, band in zip(value_min, value_max, range(img.shape[2])):
        img[:,:,band] = (img[:,:,band] - vmin)/ (vmax - vmin)
    ax.imshow(img)
    ax.axis('off')
    ax.set_title('Sector %d° (%s)' % (sector, degrees_to_cardinal(sector)))
    sector += 360/8
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).