Notebook

In [1]:

%load_ext autoreload
%autoreload 2

SLIP¶

a Simple Library for Image Processing¶

Table of content

What is the SLIP package?
Installing
Reading an image
Importing the library
Filtering
Whitening
Handling a database
An application: statistics of natural images

What is the SLIP package?¶

The SLIP library (https://github.com/bicv/SLIP) defines a simple object-oriented class for gray-scale image processing. Use it to create a SLIP object with a dedicated image size (and optionnaly some other useful parameters) - which you can use to apply common image processing routines to your images.

In [2]:

from SLIP import Image, imread
im = Image(pe='https://raw.githubusercontent.com/bicv/SLIP/master/default_param.py')
help(im)

Help on Image in module SLIP.SLIP object:

class Image(builtins.object)
 |  Image(pe='https://raw.githubusercontent.com/bicv/SLIP/master/default_param.py')
 |  
 |  This library collects different Image Processing tools.
 |  
 |  Fork me on https://github.com/bicv/SLIP !
 |  
 |  This library is used in other projects, in particular  for use with the ``LogGabor`` and ``SparseEdges`` libraries
 |  For more information check respective pages @
 |      - http://pythonhosted.org/LogGabor and
 |      - http://pythonhosted.org/SparseEdges
 |  
 |  Collects image processing routines for one given image size:
 |   - Some classical related to pure Fourier number crunching:
 |      - creating masks
 |      - normalize,
 |      - fourier_grid : defines a useful grid for generating filters in FFT
 |      - show_FT : displays the envelope and impulse response of a filter
 |      - invert : go to the other of the fourier transform
 |  - Some usual application of Fourier filtering:
 |      - trans : translation filter in Fourier space
 |      - whitening procedures
 |   - Some related to handling experiments:
 |      - load_in_database : loads a random image in a folder and
 |      - patch : takes a random patch of the correct size
 |  
 |  Methods defined here:
 |  
 |  FTfilter(self, image, FT_filter, full=False)
 |      Using the ``FTfilter`` function, it is easy to filter an image with a filter defined in Fourier space.
 |  
 |  __init__(self, pe='https://raw.githubusercontent.com/bicv/SLIP/master/default_param.py')
 |      Initializes the Image class
 |      
 |      May take as input:
 |      
 |      - a dictionary containing parameters
 |      - a ``ndarray`` (dimensions ``N_X`` and ``N_Y`` are guessed from this array)
 |      - a string representing a file or URL pointing to an image file
 |      - a string pointing to  a file or URL containing a dictionary of parameters (or simply the name of the file)
 |      - a ``NeuroTools.parameters.ParameterSet`` object containing parameters
 |      
 |      Parameters are
 |      
 |      - N_X and N_Y which are respectively the number of pixels in the vertical and horizontal dimensions respectively (MANDATORY)
 |      - optional parameters which are used in the various functions such as N_image when handling a database or the whitening parameters.
 |  
 |  dewhitening(self, white, preprocess=True, center=True, use_max=True, f_0=None)
 |      Returns the dewhitened image
 |  
 |  enveloppe_color(self, alpha)
 |  
 |  extract_patches_2d(self, image, patch_size, N_patches)
 |      Reshape a 2D image into a collection of patches
 |      
 |      redundant with self.patch, but similar call as
 |      https://github.com/scikit-learn/scikit-learn/blob/14031f6/sklearn/feature_extraction/image.py#L300
 |  
 |  fourier(self, image, full=True)
 |      Using the ``fourierr`` function, it is easy to retieve its Fourier transformation.
 |  
 |  fourier_grid(self)
 |      use that function to define a reference frame for envelopes in Fourier space.
 |      In general, it is more efficient to define dimensions as powers of 2.
 |  
 |  frequency_angle(self)
 |  
 |  frequency_radius(self)
 |  
 |  full_url(self, name_database)
 |  
 |  get_imagelist(self, exp, name_database='natural')
 |      returns an imagelist from a pickled database.
 |      
 |      If the stored imagelist does not exist, creates it.
 |      The ``exp`` string allows to tag the list for a particular experiment.
 |  
 |  get_pe(self, pe)
 |      Guesses the parameters from the init variable
 |      
 |      We perform a duck-typing to guess parameters from different possible sources.
 |      outputs a ParameterSet
 |  
 |  get_size(self, im)
 |  
 |  hist_radial_frequency(self, FT, N_f=20)
 |      A simple function to compute a radial histogram in different spatial frequency bands.
 |  
 |  imread(self, URL, resize=True)
 |  
 |  imshow(self, image, fig=None, ax=None, cmap=<matplotlib.colors.LinearSegmentedColormap object at 0x1204a5bb0>, axis=False, norm=True, center=True, xlabel='Y axis', ylabel='X axis', figsize=None, mask=False, vmin=-1, vmax=1)
 |      Plotting routine to show an image
 |      
 |      Place the [0,0] index of the array in the upper left  corner of the axes.
 |      Data limits for the axes. The default assigns zero-based row, column
 |      indices to the x, y centers of the pixels.
 |      Note that the convention for coordinates follows that of matrices: the
 |      origin is at the top left of the image, and coordinates are first the
 |      rows (vertical axis, going down) then the columns (horizontal axis,
 |      going right).
 |  
 |  init(self, mask_exponent=3.0)
 |      Initializes different convenient matrices for image processing.
 |      
 |      To be called when keeping the same Image object but changing the size of the image.
 |  
 |  init_logging(self, filename='debug.log', name='SLIP')
 |  
 |  invert(self, FT_image, full=False)
 |      # Fourier number crunching
 |  
 |  list_database(self, name_database)
 |      Returns a list of the files in a folder
 |  
 |  load_in_database(self, name_database, i_image=None, filename=None, verbose=True)
 |      Loads a random image from the database ``name_database``.
 |      
 |      The strategy is to pick one image in the folder using the list provided by the ``list_database``function.
 |      
 |      TODO: it would be useful to be able to load from standard databases such as http://www.cps.utexas.edu/natural_scenes/db.shtml
 |  
 |  low_pass(self, f_0, steepness)
 |      Returns the low_pass filter used by (Olshausen, 98)
 |      
 |      parameters from Atick (p.240)
 |      f_0 = 22 c/deg in primates: the full image is approx 45 deg
 |      alpha makes the aspect change (1=diamond on the vert and hor, 2 = anisotropic)
 |      
 |      from Olshausen 98 (p.11):
 |      f_0  = 200 cycles / image (512 x 512 images)
 |      in absolute coordinates : f_0 = 200 / 512 / 2
 |      
 |      steepness is to produce a "fairly sharp cutoff"
 |  
 |  make_imagelist(self, name_database, verbose=False)
 |      Makes a list of images with no repetition.
 |      
 |      Takes as an input the name of a database (the name of a folder in the ``datapath``),
 |      returns a list of the filenames along with the crop area.
 |  
 |  mkdir(self)
 |      Initializes two folders for storing intermediate matrices and images.
 |      
 |      To be called before any operation to store or retrieve a result or figure.
 |  
 |  normalize(self, image, center=True, use_max=True)
 |  
 |  olshausen_whitening_filt(self, f_0)
 |      Returns the whitening filter used by (Olshausen, 98)
 |  
 |  patch(self, name_database, i_image=None, filename=None, croparea=None, threshold=0.2, verbose=True, preprocess=True, center=True, use_max=True, do_whitening=False)
 |      takes a subimage of size s (a tuple)
 |      
 |      does not accept if energy is relatively below a threshold (flat image)
 |  
 |  pipeline(self, image, preprocess=True, center=True, use_max=True, do_whitening=False)
 |      pre-processing pipeline
 |  
 |  power_spectrum(self, image)
 |  
 |  preprocess(self, image)
 |      Returns the pre-processed image
 |      
 |      From raw pixelized images, we want to keep information that is relevent to the content of
 |      the objects in the image. In particular, we want to avoid:
 |      
 |          - information that would not be uniformly distributed when rotating the image. In
 |          particular, we discard information outside the unit disk in Fourier space, in particular
 |          above the Nyquist frequency,
 |          - information that relates to information of the order the size of the image. This
 |          involves discarding information at low-level frequencies.
 |      
 |      See https://laurentperrinet.github.io/sciblog/posts/2015-05-21-a-simple-pre-processing-filter-for-image-processing.html
 |      for more information.
 |  
 |  retina(self, df=0.07, sigma=0.5)
 |      A parametric description of the envelope of retinal processsing.
 |      See https://laurentperrinet.github.io/sciblog/posts/2015-05-21-a-simple-pre-processing-filter-for-image-processing.html
 |      for more information.
 |      
 |      In digital images, some of the energy in Fourier space is concentrated outside the
 |      disk corresponding to the Nyquist frequency. Let's design a filter with:
 |      
 |          - a sharp cut-off for radial frequencies higher than the Nyquist frequency,
 |          - times a smooth but sharp transition (implemented with a decaying exponential),
 |          - times a high-pass filter designed by one minus a gaussian blur.
 |      
 |      This filter is rotation invariant.
 |      
 |      Note that this filter is defined by two parameters:
 |          - one for scaling the smoothness of the transition in the high-frequency range,
 |          - one for the characteristic length of the high-pass filter.
 |      
 |      The first is defined relative to the Nyquist frequency (in absolute values) while the second
 |      is relative to the size of the image in pixels and is given in number of pixels.
 |  
 |  savefig(self, fig, fname, figpath='', formats=None, display=True)
 |  
 |  set_size(self, im)
 |      Re-initializes the Image class with  the size given in ``im``
 |      
 |      May take as input:
 |      
 |      - a numpy array,
 |      - a string representing a file or URL pointing to an image file
 |      - a tuple
 |      
 |      Updated parameters are
 |      
 |      - N_X and N_Y which are respectively the number of pixels in the vertical and horizontal dimensions respectively (MANDATORY)
 |  
 |  show_FT(self, FT_image, fig=None, figsize=None, a1=None, a2=None, axis=False, title=True, FT_title='Spectrum', im_title='Image', norm=True, vmin=-1.0, vmax=1.0)
 |  
 |  show_image_FT(self, image, FT_image, fig=None, figsize=None, a1=None, a2=None, axis=False, title=True, FT_title='Spectrum', im_title='Image', norm=True, vmin=-1.0, vmax=1.0)
 |  
 |  show_spectrum(self, image, fig=None, figsize=None, a1=None, a2=None, axis=False, title=True, FT_title='Spectrum', im_title='Image', norm=True, vmin=-1.0, vmax=1.0)
 |  
 |  trans(self, u, v)
 |  
 |  translate(self, image, vec, preshift=True)
 |      Translate image by vec (in pixels)
 |      
 |      Note that the convention for coordinates follows that of matrices: the origin is at the top left of the image, and coordinates are first the rows (vertical axis, going down) then the columns (horizontal axis, going right).
 |  
 |  whitening(self, image, f_0=None)
 |      Returns the whitened image
 |  
 |  whitening_filt(self, recompute=False, f_0=None)
 |      Returns the envelope of the whitening filter.
 |      
 |      if we chose one based on structural assumptions (``struct=True``)
 |          then we return a 1/f spectrum based on the assumption that the structure of images
 |          is self-similar and thus that the Fourier spectrum scales a priori in 1/f.
 |      
 |      elif we chose to learn,
 |          returns theaverage correlation filter in FT space.
 |      
 |          Computes the average power spectrum = FT of cross-correlation, the mean decorrelation
 |          is given for instance by (Attick, 92).
 |      
 |      else
 |          we return the parametrization based on Olshausen, 1996
 |  
 |  ----------------------------------------------------------------------
 |  Data descriptors defined here:
 |  
 |  __dict__
 |      dictionary for instance variables (if defined)
 |  
 |  __weakref__
 |      list of weak references to the object (if defined)

Installing¶

Requirements :

numpy
scipy
NeuroTools
ipython
matplotlib

To install them, use

pip install -U -r requirements.txt

Install using pip:

pip install -U SLIP

!pip install -r requirements.txt !pip install -U SLIP!pip install -U .

In [3]:

import os
import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt
fig_width = 12
figsize=(fig_width, .618*fig_width)
fopts = {'fontsize':18}
opts= {'cmap': plt.cm.gray, 'vmin':-1., 'vmax':1., 'interpolation':'nearest', 'origin':'upper'}

Reading an image¶

In [4]:

from SLIP import imread
help(imread)

Help on function imread in module SLIP.SLIP:

imread(URL, grayscale=True, rgb2gray=[0.2989, 0.587, 0.114])
    Loads whatever image. Returns a grayscale (2D) image.
    
    Note that the convention for coordinates follows that of matrices: the origin is at the top left of the image, and coordinates are first the rows (vertical axis, going down) then the columns (horizontal axis, going right).
    
    These scalar values correspond to the grayscale luminance: "The intensity of a pixel is expressed within a given range between a minimum and a maximum, inclusive. This range is represented in an abstract way as a range from 0 (total absence, black) and 1 (total presence, white), with any fractional values in between." This range is here between 0 and 1.
    
    If ``grayscale`` is True, a grayscale image is obtained by summing over channels following the formula:
    
    Y  = 0.2989 * R + 0.5870 * G + 0.1140 * B
    
    http://docs.opencv.org/2.4/modules/imgproc/doc/miscellaneous_transformations.html#cvtcolor
    which corresponds to the definition of luma:
    http://www.poynton.com/notes/colour_and_gamma/ColorFAQ.html#RTFToC11
    
    This function tries to guess at best the range and format.
    If that fails, returns a string with the error message.
    
    TODO: the above formula is an approximation of the official conversion:
    
        Y = 0.2126 * R + 0.7152 * G + 0.0722 * B
    
    in the linear RGB space.
    (see https://en.wikipedia.org/wiki/Grayscale#Colorimetric_.28luminance-preserving.29_conversion_to_grayscale)

In [5]:

list_images = ['http://static.prsa.pl/images/fad62735-b1c7-4ff4-811e-928b42fa1c89.jpg', 
                            'http://www.cosy.sbg.ac.at/~pmeerw/Watermarking/lena_color.gif', 
                            'https://upload.wikimedia.org/wikipedia/commons/e/e3/The_Horsehead_Nebula_IC434.jpg',
                            'https://upload.wikimedia.org/wikipedia/commons/thumb/a/a3/Filterstef.JPG/330px-Filterstef.JPG']
fig, ax = plt.subplots(1, len(list_images), figsize=figsize)
for i, URL in enumerate(list_images):
    ax[i].imshow(imread(URL), cmap=plt.gray())    

Importing the library¶

In [6]:

N_X, N_Y = 128, 128
N_X, N_Y = 129, 129 # testing if it works with an even size
im = Image((N_X, N_Y))
help(im.__init__)

Help on method __init__ in module SLIP.SLIP:

__init__(pe='https://raw.githubusercontent.com/bicv/SLIP/master/default_param.py') method of SLIP.SLIP.Image instance
    Initializes the Image class
    
    May take as input:
    
    - a dictionary containing parameters
    - a ``ndarray`` (dimensions ``N_X`` and ``N_Y`` are guessed from this array)
    - a string representing a file or URL pointing to an image file
    - a string pointing to  a file or URL containing a dictionary of parameters (or simply the name of the file)
    - a ``NeuroTools.parameters.ParameterSet`` object containing parameters
    
    Parameters are
    
    - N_X and N_Y which are respectively the number of pixels in the vertical and horizontal dimensions respectively (MANDATORY)
    - optional parameters which are used in the various functions such as N_image when handling a database or the whitening parameters.

Filtering¶

In [7]:

help(im.FTfilter)

Help on method FTfilter in module SLIP.SLIP:

FTfilter(image, FT_filter, full=False) method of SLIP.SLIP.Image instance
    Using the ``FTfilter`` function, it is easy to filter an image with a filter defined in Fourier space.

In [8]:

from SLIP import Image
sf_0 = 0.15
B_sf = 0.05
theta_0 = np.pi/2
B_theta = 0.15
loggabor = True

def envelope_radial(im, sf_0=sf_0, B_sf=B_sf, loggabor=loggabor, norm=True):
    if sf_0 == 0.: return 1.
    if loggabor:
        env = 1./im.f*np.exp(-.5*(np.log(im.f/sf_0)**2)/(np.log((sf_0+B_sf)/sf_0)**2))
        if norm: env /= np.sqrt((env**2).sum())
        return env
    else:
        return np.exp(-.5*(im.f - sf_0)**2/B_sf**2)

im = Image('default_param.py')
env = envelope_radial(im)

fig, ax = plt.subplots(1, 4, figsize=figsize)
for i, (f, label) in enumerate(zip([im.f_x, im.f_y, im.f, env], ['f_x', 'f_y', 'f_r', 'envelope'])):
    ax[i].matshow(f)
    ax[i].set_title(label)
    ax[i].set_xticks([])
    ax[i].set_yticks([])
    

In [9]:

def envelope_orientation(im, theta_0=theta_0, B_theta=B_theta, norm=True):
    env =  np.exp(np.cos(2*(im.f_theta-theta_0))/B_theta**2)
    if norm: env /= np.sqrt((env**2).sum())
    return env


env = envelope_radial(im) * envelope_orientation(im)

fig, ax = plt.subplots(1, 4, figsize=figsize)
for i, (f, label) in enumerate(zip([im.f_x, im.f_y, im.f, env], ['f_x', 'f_y', 'f_r', 'env'])):
    ax[i].matshow(f)
    ax[i].set_title(label)
    ax[i].set_xticks([])
    ax[i].set_yticks([])
    

In [10]:

# one can then easily generate a texture
theta0 = np.pi/2
Btheta = 0.15
theta_0 = [0, np.pi/4, np.pi/2]
B_theta = [0.1, 0.5, 1.]

def texture(env):
    return np.fft.fft2(np.fft.ifftshift(env * np.exp(1j * 2 * np.pi * np.random.rand(env.shape[0], env.shape[1])))).real

def impulse(env, phi=2 * np.pi):
    I = np.fft.fftshift(np.fft.ifft2(np.fft.ifftshift(env * np.exp(1j * phi))).real)
    I /= env.sum()
    return I

fig, ax = plt.subplots(1, 3, figsize=(fig_width, fig_width*6/18))
for i, (theta0_, label) in enumerate(zip(theta_0, [r'$\theta = 0 = \pi$', r'$\theta = \pi/4$', r'$\theta = \pi/2$']) ) :
    env = envelope_radial(im) * envelope_orientation(im, theta_0=theta0_, B_theta=Btheta)
    I = texture(env)
    ax[i].matshow(I, cmap=plt.cm.gray)
    ax[i].set_title(label)
    ax[i].set_xticks([])
    ax[i].set_yticks([])
plt.tight_layout()

In [11]:

fig, ax = plt.subplots(1, 3, figsize=(fig_width, fig_width*6/18))
for i, (Btheta_, label) in enumerate(zip(B_theta, [r'$B_\theta = 0.1$', r'$B_\theta = 0.5$', r'$B_\theta = 1.0$']) ) :
    env = envelope_radial(im) * envelope_orientation(im, theta_0=theta0, B_theta=Btheta_)
    I = texture(env)
    ax[i].matshow(I, cmap=plt.cm.gray)
    ax[i].set_title(label)
    ax[i].set_xticks([])
    ax[i].set_yticks([])
plt.tight_layout()

In [12]:

env_in = envelope_radial(im) * envelope_orientation(im)
env_V1 = envelope_radial(im) * envelope_orientation(im, theta_0=np.random.rand()*np.pi)

fig, ax = plt.subplots(1, 2, figsize=figsize)
for i, (f, label) in enumerate(zip([texture(env_in), impulse(env_V1)], [u'Texture', u'Gabor'])):
    ax[i].matshow(f, cmap=plt.cm.gray)
    ax[i].set_title(label)
    ax[i].set_xticks([])
    ax[i].set_yticks([])    

In [13]:

def convolve(image_in, image_V1):
    env_in = np.fft.fft2(image_in)
    env_V1 = np.fft.fft2(image_V1)
    return np.fft.fftshift(np.fft.ifft2((env_in*env_V1)).real)

R = convolve(texture(env_in), impulse(env_V1))
fig, ax = plt.subplots(figsize=figsize)
ax.matshow(R, cmap=plt.cm.gray)
_= plt.title(u"Convolution of a texture with a gabor")

In [14]:

# images of convolutions with differents angles
N_theta = 360
theta0 = np.pi/2
theta_0 = np.linspace(0., np.pi, 10)

for i, theta0_ in enumerate(theta_0) :
    env_in = envelope_radial(im) * envelope_orientation(im)
    env_V1 = envelope_radial(im) * envelope_orientation(im, theta_0=theta0_)
    R = convolve(texture(env_in), impulse(env_V1))
    fig, ax = plt.subplots(1, 3, figsize=figsize)
    for i, (f, label) in enumerate(zip([texture(env_in), impulse(env_V1), R], [u'Texture', u'Gabor', u'convolution'])):
        ax[i].matshow(f, cmap=plt.cm.gray)
        ax[i].set_title(label)
        ax[i].set_xticks([])
        ax[i].set_yticks([])

Whitening¶

In [15]:

help(im.whitening_filt)

Help on method whitening_filt in module SLIP.SLIP:

whitening_filt(recompute=False, f_0=None) method of SLIP.SLIP.Image instance
    Returns the envelope of the whitening filter.
    
    if we chose one based on structural assumptions (``struct=True``)
        then we return a 1/f spectrum based on the assumption that the structure of images
        is self-similar and thus that the Fourier spectrum scales a priori in 1/f.
    
    elif we chose to learn,
        returns theaverage correlation filter in FT space.
    
        Computes the average power spectrum = FT of cross-correlation, the mean decorrelation
        is given for instance by (Attick, 92).
    
    else
        we return the parametrization based on Olshausen, 1996

Testing the whitening strategy.

Generates 2 figures:

one which compares Olshausen's formulation and the actual statistical measure
one which shows the actual effect of whitening

In [16]:

import os
import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt
fig_width = 16
figsize=(fig_width, .618*fig_width)
fopts = {'fontsize':18}

In [17]:

from SLIP import Image
im = Image('default_param.py')

im.pe.white_n_learning = 100
K = im.whitening_filt(recompute=True)
f_bins, theta_bins, K_rot = im.hist_radial_frequency(K)
im.pe.white_n_learning = 0
K_ols = im.whitening_filt()
f_bins, theta_bins, K_ols_rot = im.hist_radial_frequency(K_ols)

fig = plt.figure(figsize=(fig_width, fig_width/3))
bord = 0.03
a = plt.axes([bord, bord, 1/3.-bord,  1-bord])
f_bins_m = .5*(f_bins[:-1]+f_bins[1:])
a.semilogy(f_bins_m, 1/K_rot.mean(axis=1),c='b', label='stat')
a.semilogy(f_bins_m, 1/K_ols_rot.mean(axis=1), c='k', ls='--', label='ols')
a.set_xlabel('radial frequency')
a.set_ylabel('energy')
a.set_title('mean spectrum')
a.legend(loc="upper right")

a = plt.axes([1/3.+bord, bord, 1/3.-2*bord , 1. -2*bord])
a.matshow(np.log(K))
a.axis('off')

a = plt.axes([2/3.+bord, bord, 1/3.-bord, 1-bord])
a.semilogy(f_bins_m, K_rot.mean(axis=1),c='b', label='stat')
a.semilogy(f_bins_m, K_ols_rot.mean(axis=1), c='k', ls='--', label='ols')
a.set_xlabel('radial frequency')
a.set_ylabel('energy')
a.set_title('decorrelation filter')
_ = a.legend(loc="lower right")

👾 Learning the whitening filter

In [18]:

from SLIP import Image
im = Image('default_param.py')

im.pe.white_n_learning = 100
for size in [64, 128, 256]:
    im = Image('default_param.py')
    im.pe.white_n_learning = 100
    image = im.imread('https://github.com/bicv/SLIP/raw/master/database/lena' + np.str(size) + '.png')
    print(' Testing that whitening effectively works... ')
    xcorr, xcorr_white = np.zeros((im.pe.N_X, im.pe.N_Y)), np.zeros((im.pe.N_X, im.pe.N_Y))
    for i_learning in range(im.pe.white_n_learning):
        image_patch, filename, croparea = im.patch(im.pe.white_name_database, verbose=False)
        xcorr += im.power_spectrum(image_patch)/im.pe.white_n_learning
        image_patch = im.whitening(image_patch)
        xcorr_white += im.power_spectrum(image_patch)/im.pe.white_n_learning

    middle = np.ceil(K.shape[0]/2)
    print('Figure whitening')
    fig = plt.figure(figsize=(fig_width, .618*fig_width))
    bord = .03
    # top left
    a = plt.axes([0,  1/2., 1/3., 1/2.])
    a.matshow(image, cmap=plt.gray())
    a.set_ylabel('image')
    a.axis('off')
    # top right
    a = plt.axes([2/3., 1/2., 1/3., 1/2.])
    white = im.whitening(image)
    a.matshow(white, cmap=plt.gray())
    a.set_ylabel('whitened image')
    a.axis('off')
    # middle
    a = plt.axes([1/3.+bord, .25, 1/3.-2*bord , 1/2.])
    impulse = np.zeros(image.shape)#(25,25))
    impulse[24,24] = 1
    K = im.whitening(impulse)[:49,:49]
    middle = int(np.ceil(K.shape[0]/2))
    K_middle = .5 * (K[middle,:] + K[middle+1,:])
    a.plot([24,24],[np.min(K_middle), np.max(K_middle)], c='k')
    a.plot([0, 49],[0,0], c='k')
    a.plot(K_middle, c='b')
    a.set_title('decorrelation filter')
    a.axis('off')
    # bottom left
    a = plt.axes([0,  0, 1/3., 1/2.])
    a.matshow(xcorr, cmap=plt.jet())#
    a.axis('off')
    #bottom right
    a = plt.axes([2/3., 0, 1/3., 1/2.])
    a.matshow(xcorr_white, cmap=plt.jet())
    a.axis('off')
    plt.show()

<ipython-input-18-e15f1b40aba2>:8: DeprecationWarning: `np.str` is a deprecated alias for the builtin `str`. Use `str` by itself, which is identical in behavior, to silence this warning. If you specifically wanted the numpy scalar type, use `np.str_` here.
  image = im.imread('https://github.com/bicv/SLIP/raw/master/database/lena' + np.str(size) + '.png')

 Testing that whitening effectively works... 
Figure whitening

<ipython-input-18-e15f1b40aba2>:8: DeprecationWarning: `np.str` is a deprecated alias for the builtin `str`. Use `str` by itself, which is identical in behavior, to silence this warning. If you specifically wanted the numpy scalar type, use `np.str_` here.
  image = im.imread('https://github.com/bicv/SLIP/raw/master/database/lena' + np.str(size) + '.png')

 Testing that whitening effectively works... 
Figure whitening

<ipython-input-18-e15f1b40aba2>:8: DeprecationWarning: `np.str` is a deprecated alias for the builtin `str`. Use `str` by itself, which is identical in behavior, to silence this warning. If you specifically wanted the numpy scalar type, use `np.str_` here.
  image = im.imread('https://github.com/bicv/SLIP/raw/master/database/lena' + np.str(size) + '.png')

 Testing that whitening effectively works... 
Figure whitening

In [19]:

import holoviews as hv
%load_ext holoviews.ipython
%output size=150 dpi=120
%opts Image (cmap='gray')
image = imread('https://github.com/bicv/SLIP/raw/master/database/yelmo512.png')
im.set_size(image)
image = im.preprocess(image)
hv.Image(im.normalize(image), value_dimensions=[hv.Dimension('Image', range=(-1,1))]).hist()

WARNING:param.Image02242: Setting non-parameter attribute value_dimensions=[Dimension('Image')] using a mechanism intended only for parameters

Out[19]:

In [20]:

#! whitening to balance the energy of evey frequency band
im.pe.white_n_learning = 0
white = im.whitening(image)
hv.Image(im.normalize(white)).hist()

Out[20]:

In [21]:

import imageio
imageio.imsave('database/yelmo512_w.png', white)

In [22]:

#! the filtering operation preserves infomation (none is lost...)
hv.Image(white - im.FTfilter(white, 1.)).hist()

Out[22]:

In [23]:

""" This is just to remember that we use a simpler fitering technique.

In the Attick LGN, the gain is changed according to an estimation of the SNR.

"""
print('Figure Atick')
image = plt.imread('https://github.com/bicv/LogGabor/raw/master/database/lena256.png')[:,:,0]
im.set_size(image)


contrasts = 1. / 2**np.arange(6)
freqs = np.linspace(0.,.49, 40)
size  = image.shape
x, y = np.mgrid[0:size[0],0:size[1]]

response = np.zeros((len(freqs), len(contrasts)))
for i_contrast, contrast in enumerate(contrasts):
    for i_freq, freq in enumerate(freqs):
        image = contrast * np.sin( 2* np.pi * x * freq)
        white = im.whitening(image)
        response[i_freq,i_contrast] = np.std(white)
response.shape

fig = plt.figure(figsize=(fig_width, 0.618*fig_width))
a = plt.subplot(111)
a.loglog(freqs,response,c='b')
plt.xlabel('frequency (cycles / pixel)')
_ = plt.axis('tight')

Figure Atick

defining input image as Lena¶

In [24]:

image = plt.imread('https://github.com/bicv/LogGabor/raw/master/database/lena256.png').mean(axis=-1)
print('Mean=', image.mean(), ', std=', image.std())
im = Image('default_param.py')
im.pe.N_X, im.pe.N_Y = image.shape
im.init()
image = im.normalize(image, center=False)
print('Mean=', image.mean(), ', std=', image.std())
image = im.normalize(image, center=True)
print('Mean=', image.mean(), ', std=', image.std())

Mean= 0.5018028 , std= 0.1809845
Mean= 0.5399145 , std= 0.19473016
Mean= -0.04846428 , std= 0.44376016

In [25]:

fig = plt.figure(figsize=figsize)
a1 = fig.add_subplot(121)
a2 = fig.add_subplot(122)
a1.imshow(im.normalize(image, center=True, use_max=True), **opts)
a1.set_xlabel('Y'), a1.set_ylabel('X')
a1.plot([100.], [50.], 'r*')
a1.text(100., 50., 'X=50., Y=100.')
a1.axis([0, image.shape[0], image.shape[1], 0])
a2.imshow(im.normalize(im.whitening(image), center=True, use_max=True), **opts)
a2.plot([10.], [150.], 'r*')
a2.text(10., 150., 'X=150., Y=10.')
a2.set_xlabel('Y')
v = a2.axis([0, image.shape[0], image.shape[1], 0])

In [26]:

fig = plt.figure(figsize=figsize)
a1 = fig.add_subplot(121)
a2 = fig.add_subplot(122)
a1.imshow(im.normalize(image, center=True, use_max=True), **opts)
a1.set_xlabel('Y'), a1.set_ylabel('X')
a1.axis([0, image.shape[0], image.shape[1], 0])
a2.imshow(im.normalize(im.dewhitening(im.whitening(image)), center=True, use_max=True), **opts)
a2.set_xlabel('Y')
v = a2.axis([0, image.shape[0], image.shape[1], 0])

An application: statistics of natural images¶

In [27]:

import holoviews as hv
%load_ext holoviews.ipython
%output size=150 dpi=120

The holoviews.ipython extension is already loaded. To reload it, use:
  %reload_ext holoviews.ipython

In [28]:

%opts Image (cmap='hot')
key_dims = [hv.Dimension('$f_y$', range=(0, 1)), hv.Dimension('$f_x$', range=(0, 1))]
s_f = hv.Image(im.f, group='radial frequency coordinates', key_dimensions=key_dims).hist()
s_theta = hv.Image(im.f_theta, group='radial frequency coordinates', key_dimensions=key_dims).hist()
s_f + s_theta

WARNING:param.Image03588: Setting non-parameter attribute key_dimensions=[Dimension('$f_y$'), Dimension('$f_x$')] using a mechanism intended only for parameters
WARNING:param.Image03618: Setting non-parameter attribute key_dimensions=[Dimension('$f_y$'), Dimension('$f_x$')] using a mechanism intended only for parameters

Out[28]:

In [29]:

im = Image('default_param.py')
im.pe.datapath = 'database/'
name_database = 'serre07_targets'
imagelist = im.make_imagelist(name_database=name_database)

In [30]:

F = np.zeros_like(im.f_x)
for filename, croparea in imagelist:
    image, filename_, croparea_ = im.patch(name_database, filename=filename, croparea=croparea, center=False)
    F += np.fft.fftshift(np.absolute(np.fft.fftn(image))**2)
F /= F.max()

In [31]:

%opts Image (cmap='hot')
key_dims = [hv.Dimension('$f_y$', range=(0, 1)), hv.Dimension('$f_x$', range=(0, 1))]
hv.Image(np.log(F), group='average spectral energy', key_dimensions=key_dims).hist()

WARNING:param.Image04668: Setting non-parameter attribute key_dimensions=[Dimension('$f_y$'), Dimension('$f_x$')] using a mechanism intended only for parameters

Out[31]:

In [32]:

print('F is maximum at :', np.unravel_index(np.argmax(F), dims=F.shape), ' (we expected (', F.shape[0]/2, ', ', F.shape[0]/2, ') )')

F is maximum at : (127, 128)  (we expected ( 128.0 ,  128.0 ) )

<ipython-input-32-d2ca63b8146a>:1: DeprecationWarning: 'shape' argument should be used instead of 'dims'
  print('F is maximum at :', np.unravel_index(np.argmax(F), dims=F.shape), ' (we expected (', F.shape[0]/2, ', ', F.shape[0]/2, ') )')

We now try to fit a power-law to one slice:

In [33]:

# !pip install lmfit

In [34]:

from lmfit.models import PowerLawModel
mod = PowerLawModel()
valid = ~(np.arange(F.shape[0])==F.shape[0]//2)
pars = mod.guess(F[valid, F.shape[0]//2], x=np.abs(im.f_x[valid, F.shape[0]//2]))
out  = mod.fit(F[valid, F.shape[0]//2], pars, x=np.abs(im.f_x[valid, F.shape[0]//2]))
print(out.fit_report(min_correl=0.25))

valid = ~(np.arange(F.shape[1])==F.shape[1]//2)
pars = mod.guess(F[F.shape[1]//2, valid], x=np.abs(im.f_y[F.shape[1]//2, valid]))
out  = mod.fit(F[F.shape[1]//2, valid], pars, x=np.abs(im.f_y[F.shape[1]//2, valid]))
print(out.fit_report(min_correl=0.25))

[[Model]]
    Model(powerlaw)
[[Fit Statistics]]
    # fitting method   = leastsq
    # function evals   = 51
    # data points      = 255
    # variables        = 2
    chi-square         = 0.02018965
    reduced chi-square = 7.9801e-05
    Akaike info crit   = -2404.18136
    Bayesian info crit = -2397.09884
[[Variables]]
    amplitude:  9.6345e-05 +/- 9.0392e-06 (9.38%) (init = 9.982334e-07)
    exponent:  -1.67031658 +/- 0.01724693 (1.03%) (init = -2.852235)
[[Correlations]] (unreported correlations are < 0.250)
    C(amplitude, exponent) =  0.998
[[Model]]
    Model(powerlaw)
[[Fit Statistics]]
    # fitting method   = leastsq
    # function evals   = 58
    # data points      = 255
    # variables        = 2
    chi-square         = 0.00187915
    reduced chi-square = 7.4275e-06
    Akaike info crit   = -3009.64072
    Bayesian info crit = -3002.55819
[[Variables]]
    amplitude:  4.6415e-05 +/- 2.8079e-06 (6.05%) (init = 8.104911e-07)
    exponent:  -1.66602869 +/- 0.01112401 (0.67%) (init = -2.720384)
[[Correlations]] (unreported correlations are < 0.250)
    C(amplitude, exponent) =  0.998

This proves that the spectrum of natural images falls as $\frac{1}{f^2}$. But is that true for all directions? Let's slice the Fourier spectrum along different orientations:

In [35]:

N_f = 12 #F.shape[0]/2 # making an histogram with N_f bins
f_bins = np.linspace(0., 0.5, N_f+1)
f_bins = np.logspace(-.5, 0, N_f+1, base=10)*0.5

N_orientations = 12 # making an histogram with N_f bins
theta_bins = np.linspace(0, np.pi, N_orientations, endpoint=False)

F_rot = np.zeros((N_f, N_orientations))
for i_theta in range(N_orientations):
    for i_f in range(N_f):
        f_slice = (f_bins[i_f] < im.f) *  ( im.f < f_bins[i_f+1])
        theta_slice = np.exp(np.cos(im.f_theta - theta_bins[i_theta])/(1.5*2*np.pi/N_orientations)**2)
        F_rot[i_f, i_theta] = (f_slice * theta_slice * F).sum()
        F_rot[i_f, i_theta] /= (f_slice * theta_slice).sum() # normalize by the integration area (numeric)
if np.isnan(F_rot).any(): print('Beware of the NaNs!')
F_rot /= F_rot.max()
hv.Image(np.log10(F_rot), group='average spectral energy', key_dimensions=key_dims).hist()

WARNING:param.Image04874: Setting non-parameter attribute key_dimensions=[Dimension('$f_y$'), Dimension('$f_x$')] using a mechanism intended only for parameters

Out[35]:

In [36]:

fig, axs = plt.subplots(1, 1, figsize=figsize) # fib
#axs.set_color_cycle(np.array([[0., 1., 1.]]) * np.abs(f_bins)[:, np.newaxis])
fit_exp = np.zeros(N_orientations)
for i_theta in range(N_orientations): 
    axs.loglog((f_bins[:-1]+f_bins[1:])/2, F_rot[:, i_theta], 'o', alpha=.3)
    pars = mod.guess(F_rot[:, i_theta], x=(f_bins[:-1]+f_bins[1:])/2)
    out  = mod.fit(F_rot[:, i_theta], pars, x=(f_bins[:-1]+f_bins[1:])/2)
    fit_exp[i_theta] = -out.params.get('exponent').value
axs.axis('tight')
axs.set_xlabel('spatial frequency', **fopts)
_ = axs.set_ylabel('Power', **fopts)
inset = fig.add_axes([0.6, 0.58, .3, .3], facecolor='w')
inset.bar(theta_bins, fit_exp, width=theta_bins[1]-theta_bins[0])
inset.set_xlim((0, np.pi))
inset.set_xlabel('orientation', **fopts)
inset.set_ylabel('exponent of the best-fit \n of a power-law', fontsize=10)
fit = ((f_bins[:-1]+f_bins[1:])/2)**-fit_exp.mean()
axs.loglog((f_bins[:-1]+f_bins[1:])/2, fit/fit.max(), '-', alpha=.7)
_ = inset.annotate('Average exponent = ' + str(fit_exp.mean()), 
                   xy = (0., fit_exp.mean()), xytext = (-4.5, -3), 
                   arrowprops=dict(connectionstyle="angle,angleA=90,angleB=-180,rad=100", arrowstyle="->", lw=4, fc='b', ec='k', alpha=.2), **fopts)

Notice that the average exponent is sligthly different than the one obtained above when taking just one line or column of the spectrum, as slices in orientation are different (think of slicing a pie in a partition).

more book keeping¶

In [38]:

%load_ext watermark
%watermark -i -h -m -v -p numpy,matplotlib,scipy,imageio,SLIP  -r -g -b

The watermark extension is already loaded. To reload it, use:
  %reload_ext watermark
2020-09-08T16:55:28+02:00

CPython 3.8.5
IPython 7.16.1

numpy 1.20.0.dev0+7d04e22
matplotlib 3.2.2
scipy 1.5.2
imageio 2.9.0
SLIP 20191113

compiler   : Clang 11.0.3 (clang-1103.0.32.62)
system     : Darwin
release    : 19.6.0
machine    : x86_64
processor  : i386
CPU cores  : 36
interpreter: 64bit
host name  : fortytwo
Git hash   : 2451fce2503e4e422db56dee38d5680ada3ac7a2
Git repo   : https://github.com/bicv/SLIP
Git branch : master