__init__.py

"""
The scene processing part of the toolbox defines methodes
to transform an image into a place code in the sense
of Basten and Mallot (2010).

The scene is either

* a 4d numpy array, used when the image is a equirectangular \
projection, i.e. a panoramic image.
* a 3d numpy array, used when the viewing direction can not \
be mapped/projected on a regular image (e.g. insect eye).

We thus define the following for a scene:

image based scene (IBS)
  A classical image. Each pixel is viewed in a direction
  (elevation,azimuth) in a regular manner.
  In that case the scene is a 4d numpy array
  [elevation-index,azimuth-index,channel-index,1].

Omatidium based scene (OBS)
  In an ommatidia based scene, the viewing direction
  do not need to be regularally spaced.
  In that case the scene is a 3d numpy array
  [ommatidia-index, channel-index,1].

By extension a place-code is either image based or ommatidium based.
The number of dimension of an ib-place-code is always 4, and of an
ob-place-code always 3.

image based place-code (IBPC)
    A place code derived from IBS. Each pixel is viewed in a direction
    (elevation,azimuth) in a regular manner.
    In that case the scene is a 4d numpy array
    [elevation-index,azimuth-index,channel-index,component-index].

Omatidium based place-code (OBPC)
    A place code derived from OBS, the viewing direction
    do not need to be regularally spaced.
    In that case the scene is a 3d numpy array
    [ommatidia-index, channel-index,component-index].


Abusing the terminology of a place-code, a scene can be a place-code.
Therefore ibs and obs have 4 and 3 dimension, respectively.

.. todo:

   * implement optic flow vector

"""
import numpy as np
import pandas as pd
from scipy.ndimage import maximum_filter, minimum_filter
import processing.constants as prc
import processing.tools as prt


def scene(database, posorient=None, rowid=None):
    """ Return a scene at a position orientation or given rowid
 in a given database.

    :param database: a DataBaseLoad class \
    :param posorient:  a pandas Series with index: \
['x','y','z','alpha_0,'alpha_1,'alpha_2'] (default None, i.e. not used)
    :param rowid: a row identification integer for directly reading \
in the database (default None, i.e. not used).
    :returns: a scene [elevation, azimuth, channel, 1] or \
[ommatidia,channel,1].
    :rtype: np.ndarray

    .. literalinclude:: example/processing/scene.py
       :lines: 13-14

    .. plot:: example/processing/scene.py

    """
    if posorient is not None:
        assert isinstance(posorient, pd.Series),\
            'posorient should be a pandas Series'
    scene = database.read_image(posorient=posorient,
                                rowid=rowid)
    scene = scene[..., np.newaxis]
    return scene


def skyline(scene):
    """Return the average along the elevation of a scene

    :param scene: the scenery at a given location (a 4d numpy array)
    :returns: the skyline [1,azimuth,channel,1]
    :rtype: np.ndarray

    .. literalinclude:: example/processing/skyline.py
       :lines: 12-14

    .. plot:: example/processing/skyline.py

    """
    assert prt.is_ibpc(scene),\
        'scene should be image based to compute a skyline'
    skyline = scene.mean(axis=prc.__ibpc_indeces__['elevation'])
    return skyline[np.newaxis, :]


def mztest():
  return 1


def michelson_contrast(scene, size=3):
    """Return the michelson constrast

    .. math::

       \\frac{I_\\text{max}-I_\\text{min}}{I_\\text{max}+I_\\text{min}}

    with :math:`I_\\text{max}` and :math:`I_\\text{min}` representing the \
highest and lowest luminance in an image region around each pixel.

    :param scene: an image based scene
    :param size: the size of the region to calculate the maximum \
and minimum of the local image intensity
    :returns: the michelson-contrast
    :rtype: np.ndarray

    .. literalinclude:: example/processing/michelson_contrast.py
       :lines: 12-14

    .. plot:: example/processing/michelson_contrast.py

    """
    assert prt.is_ibpc(scene), \
        'scene should be image based to compute the michelson constrast'
    contrast = np.zeros_like(scene)
    for channel in range(scene.shape[prc.__ibpc_indeces__['channel']]):
        i_max = maximum_filter(scene[..., channel, 0],
                               size=size, mode='wrap')
        i_min = minimum_filter(scene[..., channel, 0],
                               size=size, mode='wrap')
        contrast[..., channel, 0] = (i_max - i_min) / (i_max + i_min)
    return contrast


def contrast_weighted_nearness(scene, contrast_size=3, distance_channel=3):
    """Return the michelson contrast wheighted nearness

    :param scene: an image based scene
    :param contrast_size: the size of the region to calculate the maximum \
and minimum of the local image intensity in the michelson-contrast.
    :param distance_channel: the index of the distance-channel.

    .. literalinclude:: example/processing/contrast_weighted_nearness.py
       :lines: 12-14

    .. plot:: example/processing/contrast_weighted_nearness.py

    """
    assert prt.is_ibpc(scene), \
        'scene should be image based to compute the contrast weighted nearness'
    contrast = michelson_contrast(scene, size=contrast_size)
    distance = scene[..., distance_channel, 0]
    distance = distance[..., np.newaxis, np.newaxis]
    distance = np.tile(distance, (1, 1, scene.shape[-2], 1))
    return contrast / distance


def pcv(place_code, viewing_directions):
    """Place code vectors

    :param place_code: the place code at a given location (e.g. an ibs scene)
    :param viewing_directions: viewing direction of each pixel
    :returns: the place code vectors in cartesian coordinates
    :rtype: (np.ndarray)

    .. literalinclude:: example/processing/pcv.py
       :lines: 12-14

    .. plot:: example/processing/pcv.py

    """
    if prt.is_ibpc(place_code):
        component_dim = prc.__ibpc_indeces__['component']
        channel_dim = prc.__ibpc_indeces__['channel']
    elif prt.is_obpc(place_code):
        component_dim = prc.__obpc_indeces__['component']
        channel_dim = prc.__obpc_indeces__['channel']
    else:
        raise TypeError('place code should be either an ibpc or obpc')
    assert isinstance(viewing_directions, np.ndarray), \
        'viewing_directions should be a numpy array'
    assert place_code.shape[component_dim] == 1, \
        'the last dimension ({}) of the place-code should be 1'.format(
        place_code.shape[component_dim])
    elevation = viewing_directions[..., prc.__spherical_indeces__['elevation']]
    azimuth = viewing_directions[..., prc.__spherical_indeces__['azimuth']]
    unscaled_lv = prt.spherical_to_cartesian(elevation, azimuth, radius=1)
    scaled_lv = np.zeros_like(place_code)
    # (3,) -> (1,1,3) or (1,1,1,3) see numpy.tile
    scaled_lv = np.tile(scaled_lv, (unscaled_lv.shape[-1],))
    for channel_index in range(0, scaled_lv.shape[channel_dim]):
        radius = np.tile(place_code[..., channel_index, 0]
                         [..., np.newaxis], (scaled_lv.shape[-1],))
        scaled_lv[..., channel_index, :] = unscaled_lv * radius
    return scaled_lv


def apcv(place_code, viewing_directions):
    """Calculate the average scene vector

    :param place_code: the place code at a given location (e.g. an ibs scene)
    :param viewing_directions: viewing direction of each pixel
    :returns: the average place-code vector
    :rtype: (np.ndarray)

    .. literalinclude:: example/processing/apcv.py
       :lines: 12-14

    .. plot:: example/processing/apcv.py

    """
    scaled_lv = pcv(place_code, viewing_directions)
    if prt.is_ibpc(place_code):
        return (scaled_lv.sum(axis=0).sum(axis=0))[np.newaxis, np.newaxis, ...]
    elif prt.is_obpc(place_code):
        return (scaled_lv.sum(axis=0))[np.newaxis, ...]
    else:
        raise TypeError('place code is neither an ibpc nor obpc')


def optic_flow(place_code, viewing_directions, velocity):
    """NOT IMPLEMENTED"""
    raise NameError('Not Implemented')