siox.cpp revision 411b8ce8881e841bd236bd6aab966cbbd166bd39
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani/*
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani Copyright 2005, 2006 by Gerald Friedland, Kristian Jantz and Lars Knipping
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani Conversion to C++ for Inkscape by Bob Jamison
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani Licensed under the Apache License, Version 2.0 (the "License");
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani you may not use this file except in compliance with the License.
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani You may obtain a copy of the License at
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani http://www.apache.org/licenses/LICENSE-2.0
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani Unless required by applicable law or agreed to in writing, software
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani distributed under the License is distributed on an "AS IS" BASIS,
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani See the License for the specific language governing permissions and
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani limitations under the License.
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani */
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani#include "siox.h"
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani#include <math.h>
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani#include <stdarg.h>
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani#include <map>
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani#include <algorithm>
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahaninamespace org
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani{
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahaninamespace siox
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani{
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani//########################################################################
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani//# C L A B
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani//########################################################################
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani/**
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani * Convert integer A, R, G, B values into an pixel value.
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani */
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahanistatic unsigned long getRGB(int a, int r, int g, int b)
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani{
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani if (a<0) a=0;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani else if (a>255) a=255;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani if (r<0) r=0;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani else if (r>255) r=255;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani if (g<0) g=0;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani else if (g>255) g=255;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani if (b<0) b=0;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani else if (b>255) b=255;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani return (a<<24)|(r<<16)|(g<<8)|b;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani}
15411c0cb1192799b37ec8f25d6f30e8d7292fc6David Herrmann
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani/**
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani * Convert float A, R, G, B values (0.0-1.0) into an pixel value.
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani */
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahanistatic unsigned long getRGB(float a, float r, float g, float b)
dc75168823540076b354135f6e2de7a9a978fbcaZbigniew Jędrzejewski-Szmek{
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani return getRGB((int)(a * 256.0),
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani (int)(r * 256.0),
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani (int)(g * 256.0),
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani (int)(b * 256.0));
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani}
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani//#########################################
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani//# Root approximations for large speedup.
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani//# By njh!
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani//#########################################
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahanistatic const int ROOT_TAB_SIZE = 16;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahanistatic float cbrt_table[ROOT_TAB_SIZE +1];
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahanidouble CieLab::cbrt(double x)
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani{
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani double y = cbrt_table[int(x*ROOT_TAB_SIZE )]; // assuming x \in [0, 1]
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani y = (2.0 * y + x/(y*y))/3.0;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani y = (2.0 * y + x/(y*y))/3.0; // polish twice
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani return y;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani}
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahanistatic float qn_table[ROOT_TAB_SIZE +1];
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahanidouble CieLab::qnrt(double x)
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani{
b3ff20978a5e9c5a9fdc693128f29988b4070c8aSusant Sahani double y = qn_table[int(x*ROOT_TAB_SIZE )]; // assuming x \in [0, 1]
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani double Y = y*y;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani y = (4.0*y + x/(Y*Y))/5.0;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani Y = y*y;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani y = (4.0*y + x/(Y*Y))/5.0; // polish twice
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani return y;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani}
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahanidouble CieLab::pow24(double x)
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani{
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani double onetwo = x*qnrt(x);
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani return onetwo*onetwo;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani}
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahanistatic bool _clab_inited_ = false;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahanivoid CieLab::init()
15411c0cb1192799b37ec8f25d6f30e8d7292fc6David Herrmann{
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani if (!_clab_inited_)
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani {
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani cbrt_table[0] = pow(float(1)/float(ROOT_TAB_SIZE*2), 0.3333);
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani qn_table[0] = pow(float(1)/float(ROOT_TAB_SIZE*2), 0.2);
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani for(int i = 1; i < ROOT_TAB_SIZE +1; i++)
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani {
dc75168823540076b354135f6e2de7a9a978fbcaZbigniew Jędrzejewski-Szmek cbrt_table[i] = pow(float(i)/float(ROOT_TAB_SIZE), 0.3333);
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani qn_table[i] = pow(float(i)/float(ROOT_TAB_SIZE), 0.2);
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani }
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani _clab_inited_ = true;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani }
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani}
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani/**
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani * Construct this CieLab from a packed-pixel ARGB value
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani */
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant SahaniCieLab::CieLab(unsigned long rgb)
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani{
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani init();
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani int ir = (rgb>>16) & 0xff;
dc75168823540076b354135f6e2de7a9a978fbcaZbigniew Jędrzejewski-Szmek int ig = (rgb>> 8) & 0xff;
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani int ib = (rgb ) & 0xff;
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani float fr = ((float)ir) / 255.0;
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani float fg = ((float)ig) / 255.0;
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani float fb = ((float)ib) / 255.0;
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani //printf("fr:%f fg:%f fb:%f\n", fr, fg, fb);
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani if (fr > 0.04045)
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani //fr = (float) pow((fr + 0.055) / 1.055, 2.4);
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani fr = (float) pow24((fr + 0.055) / 1.055);
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani else
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani fr = fr / 12.92;
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani if (fg > 0.04045)
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani //fg = (float) pow((fg + 0.055) / 1.055, 2.4);
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani fg = (float) pow24((fg + 0.055) / 1.055);
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani else
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani fg = fg / 12.92;
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani if (fb > 0.04045)
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani //fb = (float) pow((fb + 0.055) / 1.055, 2.4);
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani fb = (float) pow24((fb + 0.055) / 1.055);
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani else
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani fb = fb / 12.92;
dc75168823540076b354135f6e2de7a9a978fbcaZbigniew Jędrzejewski-Szmek
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani // Use white = D65
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani const float x = fr * 0.4124 + fg * 0.3576 + fb * 0.1805;
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani const float y = fr * 0.2126 + fg * 0.7152 + fb * 0.0722;
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani const float z = fr * 0.0193 + fg * 0.1192 + fb * 0.9505;
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani float vx = x / 0.95047;
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani float vy = y;
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani float vz = z / 1.08883;
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani //printf("vx:%f vy:%f vz:%f\n", vx, vy, vz);
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani if (vx > 0.008856)
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani //vx = (float) pow(vx, 0.3333);
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani vx = (float) cbrt(vx);
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani else
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani vx = (7.787 * vx) + (16.0 / 116.0);
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani if (vy > 0.008856)
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani //vy = (float) pow(vy, 0.3333);
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani vy = (float) cbrt(vy);
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani else
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani vy = (7.787 * vy) + (16.0 / 116.0);
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani if (vz > 0.008856)
dc75168823540076b354135f6e2de7a9a978fbcaZbigniew Jędrzejewski-Szmek //vz = (float) pow(vz, 0.3333);
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani vz = (float) cbrt(vz);
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani else
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani vz = (7.787 * vz) + (16.0 / 116.0);
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani C = 0;
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani L = 116.0 * vy - 16.0;
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani A = 500.0 * (vx - vy);
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani B = 200.0 * (vy - vz);
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani}
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani/**
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani * Return this CieLab's value converted to a packed-pixel ARGB value
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani */
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahaniunsigned long CieLab::toRGB()
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani{
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani float vy = (L + 16.0) / 116.0;
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani float vx = A / 500.0 + vy;
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani float vz = vy - B / 200.0;
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani float vx3 = vx * vx * vx;
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani float vy3 = vy * vy * vy;
dc75168823540076b354135f6e2de7a9a978fbcaZbigniew Jędrzejewski-Szmek float vz3 = vz * vz * vz;
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani if (vy3 > 0.008856)
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani vy = vy3;
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani else
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani vy = (vy - 16.0 / 116.0) / 7.787;
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani if (vx3 > 0.008856)
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani vx = vx3;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani else
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani vx = (vx - 16.0 / 116.0) / 7.787;
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani if (vz3 > 0.008856)
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani vz = vz3;
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani else
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani vz = (vz - 16.0 / 116.0) / 7.787;
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani vx *= 0.95047; //use white = D65
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani vz *= 1.08883;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani float vr =(float)(vx * 3.2406 + vy * -1.5372 + vz * -0.4986);
dc75168823540076b354135f6e2de7a9a978fbcaZbigniew Jędrzejewski-Szmek float vg =(float)(vx * -0.9689 + vy * 1.8758 + vz * 0.0415);
49699bac94d24b444274f91f85c82e6fad04d029Susant Sahani float vb =(float)(vx * 0.0557 + vy * -0.2040 + vz * 1.0570);
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani if (vr > 0.0031308)
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani vr = (float)(1.055 * pow(vr, (1.0 / 2.4)) - 0.055);
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani else
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani vr = 12.92 * vr;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani if (vg > 0.0031308)
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani vg = (float)(1.055 * pow(vg, (1.0 / 2.4)) - 0.055);
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani else
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani vg = 12.92 * vg;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani if (vb > 0.0031308)
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani vb = (float)(1.055 * pow(vb, (1.0 / 2.4)) - 0.055);
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani else
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani vb = 12.92 * vb;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani return getRGB(0.0, vr, vg, vb);
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani}
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani/**
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani * Squared Euclidian distance between this and another color
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani */
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahanifloat CieLab::diffSq(const CieLab &other)
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani{
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani float sum=0.0;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani sum += (L - other.L) * (L - other.L);
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani sum += (A - other.A) * (A - other.A);
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani sum += (B - other.B) * (B - other.B);
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani return sum;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani}
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani/**
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani * Computes squared euclidian distance in CieLab space for two colors
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani * given as RGB values.
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani */
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahanifloat CieLab::diffSq(unsigned int rgb1, unsigned int rgb2)
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani{
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani CieLab c1(rgb1);
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani CieLab c2(rgb2);
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani float euclid = c1.diffSq(c2);
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani return euclid;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani}
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani/**
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani * Computes squared euclidian distance in CieLab space for two colors
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani * given as RGB values.
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani */
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahanifloat CieLab::diff(unsigned int rgb0, unsigned int rgb1)
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani{
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani return (float) sqrt(diffSq(rgb0, rgb1));
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani}
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani//########################################################################
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani//# T U P E L
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani//########################################################################
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani/**
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani * Helper class for storing the minimum distances to a cluster centroid
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani * in background and foreground and the index to the centroids in each
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani * signature for a given color.
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani */
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahaniclass Tupel {
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahanipublic:
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani Tupel()
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani {
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani minBgDist = 0.0f;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani indexMinBg = 0;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani minFgDist = 0.0f;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani indexMinFg = 0;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani }
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani Tupel(float minBgDistArg, long indexMinBgArg,
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani float minFgDistArg, long indexMinFgArg)
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani {
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani minBgDist = minBgDistArg;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani indexMinBg = indexMinBgArg;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani minFgDist = minFgDistArg;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani indexMinFg = indexMinFgArg;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani }
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani Tupel(const Tupel &other)
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani {
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani minBgDist = other.minBgDist;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani indexMinBg = other.indexMinBg;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani minFgDist = other.minFgDist;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani indexMinFg = other.indexMinFg;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani }
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani Tupel &operator=(const Tupel &other)
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani {
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani minBgDist = other.minBgDist;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani indexMinBg = other.indexMinBg;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani minFgDist = other.minFgDist;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani indexMinFg = other.indexMinFg;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani return *this;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani }
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani virtual ~Tupel()
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani {}
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani float minBgDist;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani long indexMinBg;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani float minFgDist;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani long indexMinFg;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani };
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani//########################################################################
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani//# S I O X I M A G E
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani//########################################################################
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani/**
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani * Create an image with the given width and height
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani */
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant SahaniSioxImage::SioxImage(unsigned int widthArg, unsigned int heightArg)
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani{
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani init(width, height);
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani}
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani/**
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani * Copy constructor
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani */
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant SahaniSioxImage::SioxImage(const SioxImage &other)
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani{
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani pixdata = NULL;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani cmdata = NULL;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani assign(other);
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani}
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani/**
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani * Assignment
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani */
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant SahaniSioxImage &SioxImage::operator=(const SioxImage &other)
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani{
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani assign(other);
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani return *this;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani}
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani/**
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani * Clean up after use.
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani */
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant SahaniSioxImage::~SioxImage()
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani{
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani if (pixdata) delete[] pixdata;
9c8e3101ceef00342263d57dce3fa61956824985Lennart Poettering if (cmdata) delete[] cmdata;
9c8e3101ceef00342263d57dce3fa61956824985Lennart Poettering}
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani/**
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani * Returns true if the previous operation on this image
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani * was successful, else false.
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani */
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahanibool SioxImage::isValid()
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani{
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani return valid;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani}
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani/**
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani * Sets whether an operation was successful, and whether
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani * this image should be considered a valid one.
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani * was successful, else false.
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani */
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahanivoid SioxImage::setValid(bool val)
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani{
9c8e3101ceef00342263d57dce3fa61956824985Lennart Poettering valid = val;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani}
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani/**
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani * Set a pixel at the x,y coordinates to the given value.
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani * If the coordinates are out of range, do nothing.
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani */
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahanivoid SioxImage::setPixel(unsigned int x,
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani unsigned int y,
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani unsigned int pixval)
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani{
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani if (x > width || y > height)
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani return;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani unsigned long offset = width * y + x;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani pixdata[offset] = pixval;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani}
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani/**
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani * Set a pixel at the x,y coordinates to the given r, g, b values.
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani * If the coordinates are out of range, do nothing.
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani */
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahanivoid SioxImage::setPixel(unsigned int x, unsigned int y,
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani unsigned int a,
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani unsigned int r,
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani unsigned int g,
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani unsigned int b)
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani{
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani if (x > width || y > height)
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani return;
9c8e3101ceef00342263d57dce3fa61956824985Lennart Poettering unsigned long offset = width * y + x;
9c8e3101ceef00342263d57dce3fa61956824985Lennart Poettering unsigned int pixval = ((a << 24) & 0xff000000) |
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani ((r << 16) & 0x00ff0000) |
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani ((g << 8) & 0x0000ff00) |
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani ((b ) & 0x000000ff);
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani pixdata[offset] = pixval;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani}
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani/**
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani * Get a pixel at the x,y coordinates given. If
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani * the coordinates are out of range, return 0;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani */
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahaniunsigned int SioxImage::getPixel(unsigned int x, unsigned int y)
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani{
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani if (x > width || y > height)
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani return 0L;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani unsigned long offset = width * y + x;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani return pixdata[offset];
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani}
e7a2419a2ae2a8f56a3e2840f8d623d2a449277aDavid Herrmann
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani/**
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani * Return the image data buffer
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani */
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahaniunsigned int *SioxImage::getImageData()
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani{
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani return pixdata;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani}
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani/**
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani * Set a confidence value at the x,y coordinates to the given value.
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani * If the coordinates are out of range, do nothing.
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani */
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahanivoid SioxImage::setConfidence(unsigned int x,
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani unsigned int y,
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani float confval)
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani{
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani if (x > width || y > height)
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani return;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani unsigned long offset = width * y + x;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani cmdata[offset] = confval;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani}
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani/**
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani * Get a confidence valueat the x,y coordinates given. If
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani * the coordinates are out of range, return 0;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani */
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahanifloat SioxImage::getConfidence(unsigned int x, unsigned int y)
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani{
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani if (x > width || y > height)
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani return 0.0;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani unsigned long offset = width * y + x;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani return cmdata[offset];
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani}
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani/**
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani * Return the confidence data buffer
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani */
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahanifloat *SioxImage::getConfidenceData()
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani{
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani return cmdata;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani}
9c8e3101ceef00342263d57dce3fa61956824985Lennart Poettering
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani/**
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani * Return the width of this image
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani */
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahaniint SioxImage::getWidth()
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani{
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani return width;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani}
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani/**
7d4866548d028489d84c39af1bb9206842a77b2bDavid Herrmann * Return the height of this image
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani */
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahaniint SioxImage::getHeight()
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani{
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani return height;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani}
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani/**
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani * Initialize values. Used by constructors
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani */
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahanivoid SioxImage::init(unsigned int widthArg, unsigned int heightArg)
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani{
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani valid = true;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani width = widthArg;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani height = heightArg;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani imageSize = width * height;
9c8e3101ceef00342263d57dce3fa61956824985Lennart Poettering pixdata = new unsigned int[imageSize];
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani cmdata = new float[imageSize];
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani for (unsigned long i=0 ; i<imageSize ; i++)
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani {
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani pixdata[i] = 0;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani cmdata[i] = 0.0;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani }
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani}
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani/**
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani * Assign values to that of another
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani */
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahanivoid SioxImage::assign(const SioxImage &other)
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani{
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani if (pixdata) delete[] pixdata;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani if (cmdata) delete[] cmdata;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani valid = other.valid;
ad1ad5c8e36ea795034fcdac660b15d7c141d55bSusant Sahani width = other.width;
height = other.height;
imageSize = width * height;
pixdata = new unsigned int[imageSize];
cmdata = new float[imageSize];
for (unsigned long i=0 ; i<imageSize ; i++)
{
pixdata[i] = other.pixdata[i];
cmdata[i] = other.cmdata[i];
}
}
/**
* Write the image to a PPM file
*/
bool SioxImage::writePPM(const std::string fileName)
{
FILE *f = fopen(fileName.c_str(), "wb");
if (!f)
return false;
fprintf(f, "P6 %d %d 255\n", width, height);
for (unsigned int y=0 ; y<height; y++)
{
for (unsigned int x=0 ; x<width ; x++)
{
unsigned int rgb = getPixel(x, y);
//unsigned int alpha = (rgb>>24) & 0xff;
unsigned int r = ((rgb>>16) & 0xff);
unsigned int g = ((rgb>> 8) & 0xff);
unsigned int b = ((rgb ) & 0xff);
fputc((unsigned char) r, f);
fputc((unsigned char) g, f);
fputc((unsigned char) b, f);
}
}
fclose(f);
return true;
}
#ifdef HAVE_GLIB
/**
* Create an image from a GdkPixbuf
*/
SioxImage::SioxImage(GdkPixbuf *buf)
{
if (!buf)
return;
unsigned int width = gdk_pixbuf_get_width(buf);
unsigned int height = gdk_pixbuf_get_height(buf);
init(width, height); //DO THIS NOW!!
guchar *pixldata = gdk_pixbuf_get_pixels(buf);
int rowstride = gdk_pixbuf_get_rowstride(buf);
int n_channels = gdk_pixbuf_get_n_channels(buf);
//### Fill in the cells with RGB values
int row = 0;
for (unsigned int y=0 ; y<height ; y++)
{
guchar *p = pixldata + row;
for (unsigned int x=0 ; x<width ; x++)
{
int r = (int)p[0];
int g = (int)p[1];
int b = (int)p[2];
int alpha = (int)p[3];
setPixel(x, y, alpha, r, g, b);
p += n_channels;
}
row += rowstride;
}
}
/**
* Create a GdkPixbuf from this image
*/
GdkPixbuf *SioxImage::getGdkPixbuf()
{
bool has_alpha = true;
int n_channels = has_alpha ? 4 : 3;
guchar *pixdata = (guchar *)
malloc(sizeof(guchar) * width * height * n_channels);
if (!pixdata)
return NULL;
int rowstride = width * n_channels;
GdkPixbuf *buf = gdk_pixbuf_new_from_data(pixdata,
GDK_COLORSPACE_RGB,
has_alpha, 8, width, height,
rowstride, NULL, NULL);
//### Fill in the cells with RGB values
int row = 0;
for (unsigned int y=0 ; y < height ; y++)
{
guchar *p = pixdata + row;
for (unsigned x=0 ; x < width ; x++)
{
unsigned int rgb = getPixel(x, y);
p[0] = (rgb >> 16) & 0xff;//r
p[1] = (rgb >> 8) & 0xff;//g
p[2] = (rgb ) & 0xff;//b
if ( n_channels > 3 )
{
p[3] = (rgb >> 24) & 0xff;//a
}
p += n_channels;
}
row += rowstride;
}
return buf;
}
#endif /* GLIB */
//########################################################################
//# S I O X
//########################################################################
//##############
//## PUBLIC
//##############
/**
* Confidence corresponding to a certain foreground region (equals one).
*/
const float Siox::CERTAIN_FOREGROUND_CONFIDENCE=1.0f;
/**
* Confidence for a region likely being foreground.
*/
const float Siox::FOREGROUND_CONFIDENCE=0.8f;
/**
* Confidence for foreground or background type being equally likely.
*/
const float Siox::UNKNOWN_REGION_CONFIDENCE=0.5f;
/**
* Confidence for a region likely being background.
*/
const float Siox::BACKGROUND_CONFIDENCE=0.1f;
/**
* Confidence corresponding to a certain background reagion (equals zero).
*/
const float Siox::CERTAIN_BACKGROUND_CONFIDENCE=0.0f;
/**
* Construct a Siox engine
*/
Siox::Siox()
{
sioxObserver = NULL;
init();
}
/**
* Construct a Siox engine
*/
Siox::Siox(SioxObserver *observer)
{
init();
sioxObserver = observer;
}
/**
*
*/
Siox::~Siox()
{
cleanup();
}
/**
* Error logging
*/
void Siox::error(char *fmt, ...)
{
char msgbuf[256];
va_list args;
va_start(args, fmt);
vsnprintf(msgbuf, 255, fmt, args);
va_end(args) ;
#ifdef HAVE_GLIB
g_warning("Siox error: %s\n", msgbuf);
#else
fprintf(stderr, "Siox error: %s\n", msgbuf);
#endif
}
/**
* Trace logging
*/
void Siox::trace(char *fmt, ...)
{
char msgbuf[256];
va_list args;
va_start(args, fmt);
vsnprintf(msgbuf, 255, fmt, args);
va_end(args) ;
#ifdef HAVE_GLIB
g_message("Siox: %s\n", msgbuf);
#else
fprintf(stdout, "Siox: %s\n", msgbuf);
#endif
}
/**
* Progress reporting
*/
bool Siox::progressReport(float percentCompleted)
{
if (!sioxObserver)
return true;
bool ret = sioxObserver->progress(percentCompleted);
if (!ret)
{
trace("User selected abort");
keepGoing = false;
}
return ret;
}
/**
* Extract the foreground of the original image, according
* to the values in the confidence matrix.
*/
SioxImage Siox::extractForeground(const SioxImage &originalImage,
unsigned int backgroundFillColor)
{
trace("### Start");
init();
keepGoing = true;
SioxImage workImage = originalImage;
//# fetch some info from the image
width = workImage.getWidth();
height = workImage.getHeight();
pixelCount = width * height;
image = workImage.getImageData();
cm = workImage.getConfidenceData();
labelField = new int[pixelCount];
trace("### Creating signatures");
//#### create color signatures
std::vector<CieLab> knownBg;
std::vector<CieLab> knownFg;
CieLab *imageClab = new CieLab[pixelCount];
for (unsigned long i=0 ; i<pixelCount ; i++)
{
float conf = cm[i];
unsigned int pix = image[i];
CieLab lab(pix);
imageClab[i] = lab;
if (conf <= BACKGROUND_CONFIDENCE)
knownBg.push_back(lab);
else if (conf >= FOREGROUND_CONFIDENCE)
knownFg.push_back(lab);
}
/*
std::vector<CieLab> imageClab;
for (int y = 0 ; y < workImage.getHeight() ; y++)
for (int x = 0 ; x < workImage.getWidth() ; x++)
{
float cm = workImage.getConfidence(x, y);
unsigned int pix = workImage.getPixel(x, y);
CieLab lab(pix);
imageClab.push_back(lab);
if (cm <= BACKGROUND_CONFIDENCE)
knownBg.push_back(lab); //note: uses CieLab(rgb)
else if (cm >= FOREGROUND_CONFIDENCE)
knownFg.push_back(lab);
}
*/
if (!progressReport(10.0))
{
error("User aborted");
workImage.setValid(false);
delete[] imageClab;
delete[] labelField;
return workImage;
}
trace("knownBg:%d knownFg:%d", knownBg.size(), knownFg.size());
std::vector<CieLab> bgSignature ;
if (!colorSignature(knownBg, bgSignature, 3))
{
error("Could not create background signature");
workImage.setValid(false);
delete[] imageClab;
delete[] labelField;
return workImage;
}
if (!progressReport(30.0))
{
error("User aborted");
workImage.setValid(false);
delete[] imageClab;
delete[] labelField;
return workImage;
}
std::vector<CieLab> fgSignature ;
if (!colorSignature(knownFg, fgSignature, 3))
{
error("Could not create foreground signature");
workImage.setValid(false);
delete[] imageClab;
delete[] labelField;
return workImage;
}
//trace("### bgSignature:%d", bgSignature.size());
if (bgSignature.size() < 1)
{
// segmentation impossible
error("Signature size is < 1. Segmentation is impossible");
workImage.setValid(false);
delete[] imageClab;
delete[] labelField;
return workImage;
}
if (!progressReport(30.0))
{
error("User aborted");
workImage.setValid(false);
delete[] imageClab;
delete[] labelField;
return workImage;
}
// classify using color signatures,
// classification cached in hashmap for drb and speedup purposes
trace("### Analyzing image");
std::map<unsigned int, Tupel> hs;
unsigned int progressResolution = pixelCount / 10;
for (unsigned int i=0; i<pixelCount; i++)
{
if (i % progressResolution == 0)
{
float progress =
30.0 + 60.0 * (float)i / (float)pixelCount;
//trace("### progress:%f", progress);
if (!progressReport(progress))
{
error("User aborted");
delete[] imageClab;
delete[] labelField;
workImage.setValid(false);
return workImage;
}
}
if (cm[i] >= FOREGROUND_CONFIDENCE)
{
cm[i] = CERTAIN_FOREGROUND_CONFIDENCE;
}
else if (cm[i] <= BACKGROUND_CONFIDENCE)
{
cm[i] = CERTAIN_BACKGROUND_CONFIDENCE;
}
else // somewhere in between
{
bool isBackground = true;
std::map<unsigned int, Tupel>::iterator iter = hs.find(i);
if (iter != hs.end()) //found
{
Tupel tupel = iter->second;
isBackground = tupel.minBgDist <= tupel.minFgDist;
}
else
{
CieLab lab = imageClab[i];
float minBg = lab.diffSq(bgSignature[0]);
int minIndex = 0;
for (unsigned int j=1; j<bgSignature.size() ; j++)
{
float d = lab.diffSq(bgSignature[j]);
if (d<minBg)
{
minBg = d;
minIndex = j;
}
}
Tupel tupel(0.0f, 0, 0.0f, 0);
tupel.minBgDist = minBg;
tupel.indexMinBg = minIndex;
float minFg = 1.0e6f;
minIndex = -1;
for (unsigned int j = 0 ; j < fgSignature.size() ; j++)
{
float d = lab.diffSq(fgSignature[j]);
if (d < minFg)
{
minFg = d;
minIndex = j;
}
}
tupel.minFgDist = minFg;
tupel.indexMinFg = minIndex;
if (fgSignature.size() == 0)
{
isBackground = (minBg <= clusterSize);
// remove next line to force behaviour of old algorithm
//error("foreground signature does not exist");
//delete[] labelField;
//workImage.setValid(false);
//return workImage;
}
else
{
isBackground = minBg < minFg;
}
hs[image[i]] = tupel;
}
if (isBackground)
cm[i] = CERTAIN_BACKGROUND_CONFIDENCE;
else
cm[i] = CERTAIN_FOREGROUND_CONFIDENCE;
}
}
delete[] imageClab;
trace("### postProcessing");
//#### postprocessing
smooth(cm, width, height, 0.333f, 0.333f, 0.333f); // average
normalizeMatrix(cm, pixelCount);
erode(cm, width, height);
keepOnlyLargeComponents(UNKNOWN_REGION_CONFIDENCE, 1.0/*sizeFactorToKeep*/);
//for (int i=0; i < 2/*smoothness*/; i++)
// smooth(cm, width, height, 0.333f, 0.333f, 0.333f); // average
normalizeMatrix(cm, pixelCount);
for (unsigned int i=0; i < pixelCount; i++)
{
if (cm[i] >= UNKNOWN_REGION_CONFIDENCE)
cm[i] = CERTAIN_FOREGROUND_CONFIDENCE;
else
cm[i] = CERTAIN_BACKGROUND_CONFIDENCE;
}
keepOnlyLargeComponents(UNKNOWN_REGION_CONFIDENCE, 1.5/*sizeFactorToKeep*/);
fillColorRegions();
dilate(cm, width, height);
if (!progressReport(100.0))
{
error("User aborted");
delete[] labelField;
workImage.setValid(false);
return workImage;
}
//#### We are done. Now clear everything but the background
for (unsigned long i = 0; i<pixelCount ; i++)
{
float conf = cm[i];
if (conf < FOREGROUND_CONFIDENCE)
image[i] = backgroundFillColor;
}
delete[] labelField;
trace("### Done");
keepGoing = false;
return workImage;
}
//##############
//## PRIVATE
//##############
/**
* Initialize the Siox engine to its 'pristine' state.
* Performed at the beginning of extractForeground().
*/
void Siox::init()
{
limits[0] = 0.64f;
limits[1] = 1.28f;
limits[2] = 2.56f;
float negLimits[3];
negLimits[0] = -limits[0];
negLimits[1] = -limits[1];
negLimits[2] = -limits[2];
clusterSize = sqrEuclidianDist(limits, 3, negLimits);
}
/**
* Clean up any debris from processing.
*/
void Siox::cleanup()
{
}
/**
* Stage 1 of the color signature work. 'dims' will be either
* 2 for grays, or 3 for colors
*/
void Siox::colorSignatureStage1(CieLab *points,
unsigned int leftBase,
unsigned int rightBase,
unsigned int recursionDepth,
unsigned int *clusterCount,
const unsigned int dims)
{
unsigned int currentDim = recursionDepth % dims;
CieLab point = points[leftBase];
float min = point(currentDim);
float max = min;
for (unsigned int i = leftBase + 1; i < rightBase ; i++)
{
point = points[i];
float curval = point(currentDim);
if (curval < min) min = curval;
if (curval > max) max = curval;
}
//Do the Rubner-rule split (sounds like a dance)
if (max - min > limits[currentDim])
{
float pivotPoint = (min + max) / 2.0; //average
unsigned int left = leftBase;
unsigned int right = rightBase - 1;
//# partition points according to the dimension
while (true)
{
while ( true )
{
point = points[left];
if (point(currentDim) > pivotPoint)
break;
left++;
}
while ( true )
{
point = points[right];
if (point(currentDim) <= pivotPoint)
break;
right--;
}
if (left > right)
break;
point = points[left];
points[left] = points[right];
points[right] = point;
left++;
right--;
}
//# Recurse and create sub-trees
colorSignatureStage1(points, leftBase, left,
recursionDepth + 1, clusterCount, dims);
colorSignatureStage1(points, left, rightBase,
recursionDepth + 1, clusterCount, dims);
}
else
{
//create a leaf
CieLab newpoint;
newpoint.C = rightBase - leftBase;
for (; leftBase < rightBase ; leftBase++)
{
newpoint.add(points[leftBase]);
}
//printf("clusters:%d\n", *clusters);
if (newpoint.C != 0)
newpoint.mul(1.0 / (float)newpoint.C);
points[*clusterCount] = newpoint;
(*clusterCount)++;
}
}
/**
* Stage 2 of the color signature work
*/
void Siox::colorSignatureStage2(CieLab *points,
unsigned int leftBase,
unsigned int rightBase,
unsigned int recursionDepth,
unsigned int *clusterCount,
const float threshold,
const unsigned int dims)
{
unsigned int currentDim = recursionDepth % dims;
CieLab point = points[leftBase];
float min = point(currentDim);
float max = min;
for (unsigned int i = leftBase+ 1; i < rightBase; i++)
{
point = points[i];
float curval = point(currentDim);
if (curval < min) min = curval;
if (curval > max) max = curval;
}
//Do the Rubner-rule split (sounds like a dance)
if (max - min > limits[currentDim])
{
float pivotPoint = (min + max) / 2.0; //average
unsigned int left = leftBase;
unsigned int right = rightBase - 1;
//# partition points according to the dimension
while (true)
{
while ( true )
{
point = points[left];
if (point(currentDim) > pivotPoint)
break;
left++;
}
while ( true )
{
point = points[right];
if (point(currentDim) <= pivotPoint)
break;
right--;
}
if (left > right)
break;
point = points[left];
points[left] = points[right];
points[right] = point;
left++;
right--;
}
//# Recurse and create sub-trees
colorSignatureStage2(points, leftBase, left,
recursionDepth + 1, clusterCount, threshold, dims);
colorSignatureStage2(points, left, rightBase,
recursionDepth + 1, clusterCount, threshold, dims);
}
else
{
//### Create a leaf
unsigned int sum = 0;
for (unsigned int i = leftBase; i < rightBase; i++)
sum += points[i].C;
if ((float)sum >= threshold)
{
float scale = (float)(rightBase - leftBase);
CieLab newpoint;
for (; leftBase < rightBase; leftBase++)
newpoint.add(points[leftBase]);
if (scale != 0.0)
newpoint.mul(1.0 / scale);
points[*clusterCount] = newpoint;
(*clusterCount)++;
}
}
}
/**
* Main color signature method
*/
bool Siox::colorSignature(const std::vector<CieLab> &inputVec,
std::vector<CieLab> &result,
const unsigned int dims)
{
unsigned int length = inputVec.size();
if (length < 1) // no error. just don't do anything
return true;
CieLab *input = (CieLab *) malloc(length * sizeof(CieLab));
if (!input)
{
error("Could not allocate buffer for signature");
return false;
}
for (unsigned int i=0 ; i < length ; i++)
input[i] = inputVec[i];
unsigned int stage1length = 0;
colorSignatureStage1(input, 0, length, 0,
&stage1length, dims);
unsigned int stage2length = 0;
colorSignatureStage2(input, 0, stage1length, 0,
&stage2length, length * 0.001, dims);
result.clear();
for (unsigned int i=0 ; i < stage2length ; i++)
result.push_back(input[i]);
free(input);
return true;
}
/**
*
*/
void Siox::keepOnlyLargeComponents(float threshold,
double sizeFactorToKeep)
{
for (unsigned long idx = 0 ; idx<pixelCount ; idx++)
labelField[idx] = -1;
int curlabel = 0;
int maxregion= 0;
int maxblob = 0;
// slow but easy to understand:
std::vector<int> labelSizes;
for (unsigned long i=0 ; i<pixelCount ; i++)
{
int regionCount = 0;
if (labelField[i] == -1 && cm[i] >= threshold)
{
regionCount = depthFirstSearch(i, threshold, curlabel++);
labelSizes.push_back(regionCount);
}
if (regionCount>maxregion)
{
maxregion = regionCount;
maxblob = curlabel-1;
}
}
for (unsigned int i=0 ; i<pixelCount ; i++)
{
if (labelField[i] != -1)
{
// remove if the component is to small
if (labelSizes[labelField[i]] * sizeFactorToKeep < maxregion)
cm[i] = CERTAIN_BACKGROUND_CONFIDENCE;
// add maxblob always to foreground
if (labelField[i] == maxblob)
cm[i] = CERTAIN_FOREGROUND_CONFIDENCE;
}
}
}
int Siox::depthFirstSearch(int startPos,
float threshold, int curLabel)
{
// stores positions of labeled pixels, where the neighbours
// should still be checked for processing:
//trace("startPos:%d threshold:%f curLabel:%d",
// startPos, threshold, curLabel);
std::vector<int> pixelsToVisit;
int componentSize = 0;
if (labelField[startPos]==-1 && cm[startPos]>=threshold)
{
labelField[startPos] = curLabel;
componentSize++;
pixelsToVisit.push_back(startPos);
}
while (pixelsToVisit.size() > 0)
{
int pos = pixelsToVisit[pixelsToVisit.size() - 1];
pixelsToVisit.erase(pixelsToVisit.end() - 1);
unsigned int x = pos % width;
unsigned int y = pos / width;
// check all four neighbours
int left = pos - 1;
if (((int)x)-1>=0 && labelField[left]==-1 && cm[left]>=threshold)
{
labelField[left]=curLabel;
componentSize++;
pixelsToVisit.push_back(left);
}
int right = pos + 1;
if (x+1 < width && labelField[right]==-1 && cm[right]>=threshold)
{
labelField[right]=curLabel;
componentSize++;
pixelsToVisit.push_back(right);
}
int top = pos - width;
if (((int)y)-1>=0 && labelField[top]==-1 && cm[top]>=threshold)
{
labelField[top]=curLabel;
componentSize++;
pixelsToVisit.push_back(top);
}
int bottom = pos + width;
if (y+1 < height && labelField[bottom]==-1
&& cm[bottom]>=threshold)
{
labelField[bottom]=curLabel;
componentSize++;
pixelsToVisit.push_back(bottom);
}
}
return componentSize;
}
/**
*
*/
void Siox::fillColorRegions()
{
for (unsigned long idx = 0 ; idx<pixelCount ; idx++)
labelField[idx] = -1;
//int maxRegion=0; // unused now
std::vector<int> pixelsToVisit;
for (unsigned long i=0; i<pixelCount; i++)
{ // for all pixels
if (labelField[i]!=-1 || cm[i]<UNKNOWN_REGION_CONFIDENCE)
{
continue; // already visited or bg
}
unsigned int origColor = image[i];
unsigned long curLabel = i+1;
labelField[i] = curLabel;
cm[i] = CERTAIN_FOREGROUND_CONFIDENCE;
// int componentSize = 1;
pixelsToVisit.push_back(i);
// depth first search to fill region
while (pixelsToVisit.size() > 0)
{
int pos = pixelsToVisit[pixelsToVisit.size() - 1];
pixelsToVisit.erase(pixelsToVisit.end() - 1);
unsigned int x=pos % width;
unsigned int y=pos / width;
// check all four neighbours
int left = pos-1;
if (((int)x)-1 >= 0 && labelField[left] == -1
&& CieLab::diff(image[left], origColor)<1.0)
{
labelField[left]=curLabel;
cm[left]=CERTAIN_FOREGROUND_CONFIDENCE;
// ++componentSize;
pixelsToVisit.push_back(left);
}
int right = pos+1;
if (x+1 < width && labelField[right]==-1
&& CieLab::diff(image[right], origColor)<1.0)
{
labelField[right]=curLabel;
cm[right]=CERTAIN_FOREGROUND_CONFIDENCE;
// ++componentSize;
pixelsToVisit.push_back(right);
}
int top = pos - width;
if (((int)y)-1>=0 && labelField[top]==-1
&& CieLab::diff(image[top], origColor)<1.0)
{
labelField[top]=curLabel;
cm[top]=CERTAIN_FOREGROUND_CONFIDENCE;
// ++componentSize;
pixelsToVisit.push_back(top);
}
int bottom = pos + width;
if (y+1 < height && labelField[bottom]==-1
&& CieLab::diff(image[bottom], origColor)<1.0)
{
labelField[bottom]=curLabel;
cm[bottom]=CERTAIN_FOREGROUND_CONFIDENCE;
// ++componentSize;
pixelsToVisit.push_back(bottom);
}
}
//if (componentSize>maxRegion) {
// maxRegion=componentSize;
//}
}
}
/**
* Applies the morphological dilate operator.
*
* Can be used to close small holes in the given confidence matrix.
*/
void Siox::dilate(float *cm, int xres, int yres)
{
for (int y=0; y<yres; y++)
{
for (int x=0; x<xres-1; x++)
{
int idx=(y*xres)+x;
if (cm[idx+1]>cm[idx])
cm[idx]=cm[idx+1];
}
}
for (int y=0; y<yres; y++)
{
for (int x=xres-1; x>=1; x--)
{
int idx=(y*xres)+x;
if (cm[idx-1]>cm[idx])
cm[idx]=cm[idx-1];
}
}
for (int y=0; y<yres-1; y++)
{
for (int x=0; x<xres; x++)
{
int idx=(y*xres)+x;
if (cm[((y+1)*xres)+x] > cm[idx])
cm[idx]=cm[((y+1)*xres)+x];
}
}
for (int y=yres-1; y>=1; y--)
{
for (int x=0; x<xres; x++)
{
int idx=(y*xres)+x;
if (cm[((y-1)*xres)+x] > cm[idx])
cm[idx]=cm[((y-1)*xres)+x];
}
}
}
/**
* Applies the morphological erode operator.
*/
void Siox::erode(float *cm, int xres, int yres)
{
for (int y=0; y<yres; y++)
{
for (int x=0; x<xres-1; x++)
{
int idx=(y*xres)+x;
if (cm[idx+1] < cm[idx])
cm[idx]=cm[idx+1];
}
}
for (int y=0; y<yres; y++)
{
for (int x=xres-1; x>=1; x--)
{
int idx=(y*xres)+x;
if (cm[idx-1] < cm[idx])
cm[idx]=cm[idx-1];
}
}
for (int y=0; y<yres-1; y++)
{
for (int x=0; x<xres; x++)
{
int idx=(y*xres)+x;
if (cm[((y+1)*xres)+x] < cm[idx])
cm[idx]=cm[((y+1)*xres)+x];
}
}
for (int y=yres-1; y>=1; y--)
{
for (int x=0; x<xres; x++)
{
int idx=(y*xres)+x;
if (cm[((y-1)*xres)+x] < cm[idx])
cm[idx]=cm[((y-1)*xres)+x];
}
}
}
/**
* Normalizes the matrix to values to [0..1].
*/
void Siox::normalizeMatrix(float *cm, int cmSize)
{
float max = *std::max(cm, cm + cmSize);
if (max<=0.0 || max==1.0)
return;
float alpha=1.00f/max;
premultiplyMatrix(alpha, cm, cmSize);
}
/**
* Multiplies matrix with the given scalar.
*/
void Siox::premultiplyMatrix(float alpha, float *cm, int cmSize)
{
for (int i=0; i<cmSize; i++)
cm[i]=alpha*cm[i];
}
/**
* Blurs confidence matrix with a given symmetrically weighted kernel.
*
* In the standard case confidence matrix entries are between 0...1 and
* the weight factors sum up to 1.
*/
void Siox::smooth(float *cm, int xres, int yres,
float f1, float f2, float f3)
{
for (int y=0; y<yres; y++)
{
for (int x=0; x<xres-2; x++)
{
int idx=(y*xres)+x;
cm[idx]=f1*cm[idx]+f2*cm[idx+1]+f3*cm[idx+2];
}
}
for (int y=0; y<yres; y++)
{
for (int x=xres-1; x>=2; x--)
{
int idx=(y*xres)+x;
cm[idx]=f3*cm[idx-2]+f2*cm[idx-1]+f1*cm[idx];
}
}
for (int y=0; y<yres-2; y++)
{
for (int x=0; x<xres; x++)
{
int idx=(y*xres)+x;
cm[idx]=f1*cm[idx]+f2*cm[((y+1)*xres)+x]+f3*cm[((y+2)*xres)+x];
}
}
for (int y=yres-1; y>=2; y--)
{
for (int x=0; x<xres; x++)
{
int idx=(y*xres)+x;
cm[idx]=f3*cm[((y-2)*xres)+x]+f2*cm[((y-1)*xres)+x]+f1*cm[idx];
}
}
}
/**
* Squared Euclidian distance of p and q.
*/
float Siox::sqrEuclidianDist(float *p, int pSize, float *q)
{
float sum=0.0;
for (int i=0; i<pSize; i++)
{
float v = p[i] - q[i];
sum += v*v;
}
return sum;
}
} // namespace siox
} // namespace org
//########################################################################
//# E N D O F F I L E
//########################################################################