001/**
002 * Copyright (c) 2011, The University of Southampton and the individual contributors.
003 * All rights reserved.
004 *
005 * Redistribution and use in source and binary forms, with or without modification,
006 * are permitted provided that the following conditions are met:
007 *
008 *   *  Redistributions of source code must retain the above copyright notice,
009 *      this list of conditions and the following disclaimer.
010 *
011 *   *  Redistributions in binary form must reproduce the above copyright notice,
012 *      this list of conditions and the following disclaimer in the documentation
013 *      and/or other materials provided with the distribution.
014 *
015 *   *  Neither the name of the University of Southampton nor the names of its
016 *      contributors may be used to endorse or promote products derived from this
017 *      software without specific prior written permission.
018 *
019 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
020 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
021 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
022 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
023 * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
024 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
025 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
026 * ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
027 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
028 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
029 */
030package org.openimaj.image.processing.convolution.filterbank;
031
032import static java.lang.Math.PI;
033import static java.lang.Math.cos;
034import static java.lang.Math.exp;
035import static java.lang.Math.pow;
036import static java.lang.Math.sin;
037import static java.lang.Math.sqrt;
038
039import org.openimaj.image.FImage;
040import org.openimaj.image.processing.convolution.FConvolution;
041import org.openimaj.image.processing.convolution.Gaussian2D;
042import org.openimaj.image.processing.convolution.LaplacianOfGaussian2D;
043import org.openimaj.math.util.FloatArrayStatsUtils;
044
045/**
046 * Implementation of a the filter bank described in: T. Leung and J. Malik.
047 * Representing and recognizing the visual appearance of materials using
048 * three-dimensional textons. IJCV, 2001
049 * 
050 * Inspired by the matlab implementation from
051 * http://www.robots.ox.ac.uk/~vgg/research/texclass/filters.html
052 * 
053 * @author Jonathon Hare (jsh2@ecs.soton.ac.uk)
054 */
055public class LeungMalikFilterBank extends FilterBank {
056        /**
057         * Default constructor with a filter support of 49 pixels
058         */
059        public LeungMalikFilterBank() {
060                this(49);
061        }
062
063        /**
064         * Construct with given support (filter size).
065         * 
066         * @param size
067         *            the filter size
068         */
069        public LeungMalikFilterBank(int size) {
070                super(makeFilters(size));
071        }
072
073        protected static FConvolution[] makeFilters(int size) {
074                final int nScales = 3;
075                final int nOrientations = 6;
076
077                final int NROTINV = 12;
078                final int NBAR = nScales * nOrientations;
079                final int NEDGE = nScales * nOrientations;
080                final int NF = NBAR + NEDGE + NROTINV;
081
082                final FConvolution F[] = new FConvolution[NF];
083
084                int count = 0;
085                for (int i = 1; i <= nScales; i++) {
086                        final float scale = (float) pow(sqrt(2), i);
087
088                        for (int orient = 0; orient < nOrientations; orient++) {
089                                final float angle = (float) (PI * orient / nOrientations);
090
091                                F[count] = new FConvolution(makeFilter(scale, 0, 1, angle, size));
092                                F[count + NEDGE] = new FConvolution(makeFilter(scale, 0, 2, angle, size));
093                                count++;
094                        }
095                }
096
097                count = NBAR + NEDGE;
098                for (int i = 1; i <= 4; i++) {
099                        final float scale = (float) pow(sqrt(2), i);
100
101                        F[count] = new FConvolution(normalise(Gaussian2D.createKernelImage(size, scale)));
102                        F[count + 1] = new FConvolution(normalise(LaplacianOfGaussian2D.createKernelImage(size, scale)));
103                        F[count + 2] = new FConvolution(normalise(LaplacianOfGaussian2D.createKernelImage(size, 3 * scale)));
104                        count += 3;
105                }
106
107                return F;
108        }
109
110        protected static FImage makeFilter(float scale, int phasex, int phasey, float angle, int size) {
111                final int hs = (size - 1) / 2;
112
113                final FImage filter = new FImage(size, size);
114                for (int y = -hs, j = 0; y < hs; y++, j++) {
115                        for (int x = -hs, i = 0; x < hs; x++, i++) {
116                                final float cos = (float) cos(angle);
117                                final float sin = (float) sin(angle);
118
119                                final float rx = cos * x - sin * y;
120                                final float ry = sin * x + cos * y;
121
122                                final float gx = gaussian1D(3 * scale, 0, rx, phasex);
123                                final float gy = gaussian1D(scale, 0, ry, phasey);
124
125                                filter.pixels[j][i] = gx * gy;
126                        }
127                }
128                return normalise(filter);
129        }
130
131        protected static float gaussian1D(float sigma, float mean, float x, int order) {
132                x = x - mean;
133                final float num = x * x;
134
135                final float variance = sigma * sigma;
136                final float denom = 2 * variance;
137                final float g = (float) (exp(-num / denom) / pow(PI * denom, 0.5));
138
139                switch (order) {
140                case 0:
141                        return g;
142                case 1:
143                        return -g * (x / variance);
144                case 2:
145                        return g * ((num - variance) / (variance * variance));
146                default:
147                        throw new IllegalArgumentException("order must be 0, 1 or 2.");
148                }
149        }
150
151        protected static FImage normalise(FImage f) {
152                final float mean = FloatArrayStatsUtils.mean(f.pixels);
153                f.subtractInplace(mean);
154                final float sumabs = FloatArrayStatsUtils.sumAbs(f.pixels);
155                return f.divideInplace(sumabs);
156        }
157}