Gaussian Filter in Keras (code snippet)

Very often we need to perform basic vision operations on a computational graph like building a Laplacian pyramid or filter a tensor with a specific precalculated filter. 

Below i present a code snippet for building a fixed non-trainable gaussian filter in keras.

import keras
import numpy as np
import scipy.stats as st

def gaussian_filter_block(input_layer,
                          kernel_size=3,
                          strides=(1, 1),
                          dilation_rate=(1, 1),
                          padding="same",
                          activation=None,
                          trainable=False,
                          use_bias=False):
    """
    Build a gaussian filter block
    :return:
    """

    def _gaussian_kernel(kernlen=[21, 21], nsig=[3, 3]):
        """
        Returns a 2D Gaussian kernel array
        """
        assert len(nsig) == 2
        assert len(kernlen) == 2
        kern1d = []
        for i in range(2):
            interval = (2 * nsig[i] + 1.) / (kernlen[i])
            x = np.linspace(-nsig[i] - interval / 2., nsig[i] + interval / 2.,
                            kernlen[i] + 1)
            kern1d.append(np.diff(st.norm.cdf(x)))

        kernel_raw = np.sqrt(np.outer(kern1d[0], kern1d[1]))
        # divide by sum so they all add up to 1
        kernel = kernel_raw / kernel_raw.sum()
        return kernel

    # Initialise to set kernel to required value
    def kernel_init(shape, dtype):
        kernel = np.zeros(shape)
        kernel[:, :, 0, 0] = _gaussian_kernel([shape[0], shape[1]])
        return kernel

    return keras.layers.DepthwiseConv2D(
        kernel_size=kernel_size,
        strides=strides,
        padding=padding,
        depth_multiplier=1,
        dilation_rate=dilation_rate,
        activation=activation,
        use_bias=use_bias,
        trainable=trainable,
        depthwise_initializer=kernel_init,
        kernel_initializer=kernel_init)(input_layer)

from my open source project https://github.com/NikolasMarkou/multiscale_variational_autoencoder

Normally Distributed Random Numbers

So you have a random generator that produces uniformly distributed random numbers from 0 to 1 and you want to create a normal distributed random number generator.

        /// <summary>
        /// Returns a normally distributed random number
        /// </summary>
        /// <param name="mean">The mean value of the distribution</param>
        /// <param name="stdDev">The standard deviation of the distribution</param>
        /// <returns>A normally distributed random number</returns>
        public Double NextDouble(double mean, double stdDev)
        {
            Random rand = new Random();
            double u1 = rand.NextDouble(), u2 = rand.NextDouble();
            double randStdNormal = Math.Sqrt(-2.0 * Math.Log(u1)) *
                Math.Sin(2.0 * Math.PI * u2);   // random normal(0,1)
            double randNormal =
                mean + stdDev * randStdNormal;  // random normal(mean,stdDev^2)
            return randNormal;
        }

The only problem is that if you want to create lots of random numbers all the logarithms and square roots just bog down the machine especially if you want real time performance. So what do you do ?

Simple, trade memory for computation time. I created a simple class that is used like the Random class of the .NET framework but it internally caches normally distributed numbers upon initialization.

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace Emgu_Tracker_Advanced_V2.Source
{
    public class RandomNormalCached : Random
    {
        #region Variables
        protected int noCached = 101;
        /// <summary>
        /// Array of precomputed number with 0 mean and 1 variance
        /// </summary>
        protected double[] precomputedNumbers;
        /// <summary>
        /// Internal state counter
        /// </summary>
        protected int counter;
        #endregion

        #region Constructor
        public RandomNormalCached()
            : base()
        {
            counter = 0;
            precomputedNumbers = new double[noCached];
            for (int i = 0; i < precomputedNumbers.Length; i++)
                precomputedNumbers[i] = ComputeNormalDistribution(0, 1);
        }

        /// <summary>
        /// Constructor with custom number of cached random numbers
        /// </summary>
        /// <param name="_noCached">Number of cached random numbers</param>
        /// <param name="_seed">Random seed initializer</param>
        public RandomNormalCached(int _noCached, int _seed)
            : base(_seed)
        {
            counter = 0;
            noCached = _noCached;
            precomputedNumbers = new double[noCached];
            // compute a number of normal distributed number with zero mean , 
            // and deviation equal to one
            for (int i = 0; i < precomputedNumbers.Length; i++)
                precomputedNumbers[i] = ComputeNormalDistribution(0, 1);
        }

        public RandomNormalCached(int _seed)
            : base(_seed)
        {
            counter = 0;
            precomputedNumbers = new double[noCached];
            // compute a number of normal distributed number with zero mean , 
            // and deviation equal to one
            for (int i = 0; i < precomputedNumbers.Length; i++)
                precomputedNumbers[i] = ComputeNormalDistribution(0, 1);
        }
        #endregion

        #region Functions
        public Double NextDouble(double mean, double stdDev)
        {
            double result = mean + stdDev * precomputedNumbers[counter];
            counter = (counter + 1) % noCached;
            return result;
        }

        /// <summary>
        /// Returns a nornal distributed random number
        /// </summary>
        /// <param name="mean">The mean/expected value</param>
        /// <param name="stdDev">The standard deviation</param>
        /// <returns>Normal distributed normal number</returns>
        protected Double ComputeNormalDistribution(double mean, double stdDev)
        {
            double u1 = base.NextDouble(), u2 = base.NextDouble();
            double randStdNormal = Math.Sqrt(-2.0 * Math.Log(u1)) *
                Math.Sin(2.0 * Math.PI * u2);   // random normal(0,1)
            double randNormal =
                mean + stdDev * randStdNormal;  // random normal(mean,stdDev^2)
            return randNormal;
        }
        #endregion
    }
}