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.analysis.colour;
031
032import org.openimaj.image.FImage;
033import org.openimaj.image.MBFImage;
034import org.openimaj.image.colour.ColourSpace;
035import org.openimaj.image.combiner.ImageCombiner;
036
037/**
038 * Implementation of the CIE 2000 colour difference equation, 
039 * and a {@link ImageCombiner} to calculate a colour disparity map between
040 * two images.
041 * 
042 * @author Jonathon Hare (jsh2@ecs.soton.ac.uk)
043 */
044public class CIEDE2000 implements ImageCombiner<MBFImage, MBFImage, FImage> {
045        /**
046         * Calculate the colour difference value between two colours in lab space.
047         * @param lab1 first colour
048         * @param lab2 second colour
049         * @return the CIE 2000 colour difference
050         */
051        public static double calculateDeltaE(double [] lab1, double[] lab2) {
052                return calculateDeltaE(lab1[0],lab1[1],lab1[2],lab2[0],lab2[1],lab2[2]);
053        }
054        
055        /**
056         * Calculate the colour difference value between two colours in lab space.
057         * @param lab1 first colour
058         * @param lab2 second colour
059         * @return the CIE 2000 colour difference
060         */
061        public static float calculateDeltaE(float [] lab1, float[] lab2) {
062                return (float) calculateDeltaE(lab1[0],lab1[1],lab1[2],lab2[0],lab2[1],lab2[2]);
063        }
064        
065        /**
066         * Calculate the colour difference value between two colours in lab space.
067         * @param lab1 first colour
068         * @param lab2 second colour
069         * @return the CIE 2000 colour difference
070         */
071        public static float calculateDeltaE(Float [] lab1, Float[] lab2) {
072                return (float) calculateDeltaE(lab1[0],lab1[1],lab1[2],lab2[0],lab2[1],lab2[2]);
073        }
074        
075        /**
076         * Calculate the colour difference value between two colours in lab space.
077         * @param L1 first colour's L component
078         * @param a1 first colour's a component
079         * @param b1 first colour's b component
080         * @param L2 second colour's L component
081         * @param a2 second colour's a component
082         * @param b2 second colour's b component
083         * @return the CIE 2000 colour difference
084         */
085        public static double calculateDeltaE(double L1, double a1, double b1, double L2, double a2, double b2) {
086                double Lmean = (L1 + L2) / 2.0; //ok
087                double C1 =  Math.sqrt(a1*a1 + b1*b1); //ok
088                double C2 =  Math.sqrt(a2*a2 + b2*b2); //ok
089                double Cmean = (C1 + C2) / 2.0; //ok
090                
091                double G =  ( 1 - Math.sqrt( Math.pow(Cmean, 7) / (Math.pow(Cmean, 7) + Math.pow(25, 7)) ) ) / 2; //ok
092                double a1prime = a1 * (1 + G); //ok
093                double a2prime = a2 * (1 + G); //ok
094                
095                double C1prime =  Math.sqrt(a1prime*a1prime + b1*b1); //ok
096                double C2prime =  Math.sqrt(a2prime*a2prime + b2*b2); //ok
097                double Cmeanprime = (C1prime + C2prime) / 2; //ok 
098                
099                double h1prime =  Math.atan2(b1, a1prime) + 2*Math.PI * (Math.atan2(b1, a1prime)<0 ? 1 : 0);
100                double h2prime =  Math.atan2(b2, a2prime) + 2*Math.PI * (Math.atan2(b2, a2prime)<0 ? 1 : 0);
101                double Hmeanprime =  ((Math.abs(h1prime - h2prime) > Math.PI) ? (h1prime + h2prime + 2*Math.PI) / 2 : (h1prime + h2prime) / 2); //ok
102                
103                double T =  1.0 - 0.17 * Math.cos(Hmeanprime - Math.PI/6.0) + 0.24 * Math.cos(2*Hmeanprime) + 0.32 * Math.cos(3*Hmeanprime + Math.PI/30) - 0.2 * Math.cos(4*Hmeanprime - 21*Math.PI/60); //ok
104                
105                double deltahprime =  ((Math.abs(h1prime - h2prime) <= Math.PI) ? h2prime - h1prime : (h2prime <= h1prime) ? h2prime - h1prime + 2*Math.PI : h2prime - h1prime - 2*Math.PI); //ok
106                
107                double deltaLprime = L2 - L1; //ok
108                double deltaCprime = C2prime - C1prime; //ok
109                double deltaHprime =  2.0 * Math.sqrt(C1prime*C2prime) * Math.sin(deltahprime / 2.0); //ok
110                double SL =  1.0 + ( (0.015*(Lmean - 50)*(Lmean - 50)) / (Math.sqrt( 20 + (Lmean - 50)*(Lmean - 50) )) ); //ok
111                double SC =  1.0 + 0.045 * Cmeanprime; //ok
112                double SH =  1.0 + 0.015 * Cmeanprime * T; //ok
113                
114                double deltaTheta =  (30 * Math.PI / 180) * Math.exp(-((180/Math.PI*Hmeanprime-275)/25)*((180/Math.PI*Hmeanprime-275)/25));
115                double RC =  (2 * Math.sqrt(Math.pow(Cmeanprime, 7) / (Math.pow(Cmeanprime, 7) + Math.pow(25, 7))));
116                double RT =  (-RC * Math.sin(2 * deltaTheta));
117                
118                double KL = 1;
119                double KC = 1;
120                double KH = 1;
121                
122                double deltaE = Math.sqrt(
123                                ((deltaLprime/(KL*SL)) * (deltaLprime/(KL*SL))) +
124                                ((deltaCprime/(KC*SC)) * (deltaCprime/(KC*SC))) +
125                                ((deltaHprime/(KH*SH)) * (deltaHprime/(KH*SH))) +
126                                (RT * (deltaCprime/(KC*SC)) * (deltaHprime/(KH*SH)))
127                                );
128                                
129                return deltaE;
130        }
131
132        /**
133         * Compute the disparity map between two images. If the images
134         * are not in Lab colour space, then copies in lab space will
135         * be created.
136         * @param im1 The first image.
137         * @param im2 The second image.
138         * @return the disparity map between the colours in the two images.
139         */
140        public static FImage makeDisparityMap(MBFImage im1, MBFImage im2) {
141                if (im1.colourSpace != ColourSpace.CIE_Lab) {
142                        im1 = ColourSpace.convert(im1, ColourSpace.CIE_Lab);
143                }
144                
145                if (im2.colourSpace != ColourSpace.CIE_Lab) {
146                        im2 = ColourSpace.convert(im2, ColourSpace.CIE_Lab);
147                }
148                
149                FImage disparity = new FImage(im1.getWidth(), im1.getHeight());
150                for (int y=0; y<disparity.height; y++) {
151                        for (int x=0; x<disparity.width; x++) {
152                                disparity.pixels[y][x] = calculateDeltaE(im1.getPixel(x, y), im2.getPixel(x, y));
153                        }
154                }
155                
156                return disparity;
157        }
158        
159        @Override
160        public FImage combine(MBFImage image1, MBFImage image2) {
161                return makeDisparityMap(image2, image2);
162        }
163}