www.pudn.com > Image_segment.rar > ImageSegmentation.cpp
//Copyright (c) 2004-2005, Baris Sumengen
//All rights reserved.
//
// CIMPL Matrix Performance Library
//
//Redistribution and use in source and binary
//forms, with or without modification, are
//permitted provided that the following
//conditions are met:
//
// * No commercial use is allowed.
// This software can only be used
// for non-commercial purposes. This
// distribution is mainly intended for
// academic research and teaching.
// * Redistributions of source code must
// retain the above copyright notice, this
// list of conditions and the following
// disclaimer.
// * Redistributions of binary form must
// mention the above copyright notice, this
// list of conditions and the following
// disclaimer in a clearly visible part
// in associated product manual,
// readme, and web site of the redistributed
// software.
// * Redistributions in binary form must
// reproduce the above copyright notice,
// this list of conditions and the
// following disclaimer in the
// documentation and/or other materials
// provided with the distribution.
// * The name of Baris Sumengen may not be
// used to endorse or promote products
// derived from this software without
// specific prior written permission.
//
//THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT
//HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
//EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT
//NOT LIMITED TO, THE IMPLIED WARRANTIES OF
//MERCHANTABILITY AND FITNESS FOR A PARTICULAR
//PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
//CONTRIBUTORS BE LIABLE FOR ANY
//DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
//EXEMPLARY, OR CONSEQUENTIAL DAMAGES
//(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT
//OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
//DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
//HOWEVER CAUSED AND ON ANY THEORY OF
//LIABILITY, WHETHER IN CONTRACT, STRICT
//LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
//OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
//OF THIS SOFTWARE, EVEN IF ADVISED OF THE
//POSSIBILITY OF SUCH DAMAGE.
#include "./ImageSegmentation.h"
#include
#include
namespace ImageProcessing
{
// Gaussian Filter
Matrix Gaussian2D(int side, float sigma_x, float angle, float ratio)
{
Matrix temp(2*side+1, 2*side+1);
float sigma_y = sigma_x*ratio;
double sin_angle = sin(angle*PI/180);
double cos_angle = cos(angle*PI/180);
int sample;
double d;
if (sigma_x < 2.0 || sigma_y < 2.0)
{
// use denser sample for better filter quality
sample = 5;
d = 1.0/(2.0*sample+1.0);
}
else
{
sample = 0;
d = 1.0;
}
double T = (double) (2*sample+1)*(2*sample+1);
for (int i=0;i<2*side+1;i++)
{
for (int j=0;j<2*side+1;j++)
{
double sum_G = 0.0;
double y = (double) (j-side)-d*sample-d;
for (int sx=0;sx<2*sample+1;sx++)
{
y += d;
double x = (double) (i-side)-d*sample-d;
for (int sy=0;sy<2*sample+1;sy++)
{
x += d;
double x_new = x*cos_angle + y*sin_angle;
double y_new = -x*sin_angle + y*cos_angle;
double g = 1.0/(2.0*PI*sigma_x*sigma_y)*exp(-(x_new*x_new/(sigma_x*sigma_x)+y_new*y_new/(sigma_y*sigma_y))/2.0);
sum_G += g;
}
}
temp.ElemNC(j,i) = (float)(sum_G/T);
}
}
return temp;
}
Matrix Gaussian2D(int side, double sigma_x, double angle, double ratio)
{
Matrix temp(2*side+1, 2*side+1);
double sigma_y = sigma_x*ratio;
double sin_angle = sin(angle*PI/180);
double cos_angle = cos(angle*PI/180);
int sample;
double d;
if (sigma_x < 2.0 || sigma_y < 2.0)
{
// use denser sample for better filter quality
sample = 5;
d = 1.0/(2.0*sample+1.0);
}
else
{
sample = 0;
d = 1.0;
}
double T = (double) (2*sample+1)*(2*sample+1);
for (int i=0;i<2*side+1;i++)
{
for (int j=0;j<2*side+1;j++)
{
double sum_G = 0.0;
double y = (double) (j-side)-d*sample-d;
for (int sx=0;sx<2*sample+1;sx++)
{
y += d;
double x = (double) (i-side)-d*sample-d;
for (int sy=0;sy<2*sample+1;sy++)
{
x += d;
double x_new = x*cos_angle + y*sin_angle;
double y_new = -x*sin_angle + y*cos_angle;
double g = 1.0/(2.0*PI*sigma_x*sigma_y)*exp(-(x_new*x_new/(sigma_x*sigma_x)+y_new*y_new/(sigma_y*sigma_y))/2.0);
sum_G += g;
}
}
temp.ElemNC(j,i) = sum_G/T;
}
}
return temp;
}
// First Directional Derivative of Gaussian Filter
Matrix FDGaussian2D(int side, float sigma_x, float angle, float ratio)
{
Matrix temp(2*side+1, 2*side+1);
float sigma_y = sigma_x*ratio;
double sin_angle = sin(angle*PI/180);
double cos_angle = cos(angle*PI/180);
int sample;
double d;
if (sigma_x < 2.0 || sigma_y < 2.0)
{
// use denser sample for better filter quality
sample = 5;
d = 1.0/(2.0*sample+1.0);
}
else
{
sample = 0;
d = 1.0;
}
double T = (double) (2*sample+1)*(2*sample+1);
for (int i=0;i<2*side+1;i++)
{
for (int j=0;j<2*side+1;j++)
{
double sum_G = 0.0;
double y = (double) (j-side)-d*sample-d;
for (int sx=0;sx<2*sample+1;sx++)
{
y += d;
double x = (double) (i-side)-d*sample-d;
for (int sy=0;sy<2*sample+1;sy++)
{
x += d;
double x_new = x*cos_angle + y*sin_angle;
double y_new = -x*sin_angle + y*cos_angle;
double g = -x_new/(2.0*PI*sigma_x*sigma_y*sigma_x*sigma_x)*exp(-(x_new*x_new/(sigma_x*sigma_x)+y_new*y_new/(sigma_y*sigma_y))/2.0);
sum_G += g;
}
}
temp.ElemNC(j,i) = (float)(sum_G/T);
}
}
return temp;
}
Matrix FDGaussian2D(int side, double sigma_x, double angle, double ratio)
{
Matrix temp(2*side+1, 2*side+1);
double sigma_y = sigma_x*ratio;
double sin_angle = sin(angle*PI/180);
double cos_angle = cos(angle*PI/180);
int sample;
double d;
if (sigma_x < 2.0 || sigma_y < 2.0)
{
// use denser sample for better filter quality
sample = 5;
d = 1.0/(2.0*sample+1.0);
}
else
{
sample = 0;
d = 1.0;
}
double T = (double) (2*sample+1)*(2*sample+1);
for (int i=0;i<2*side+1;i++)
{
for (int j=0;j<2*side+1;j++)
{
double sum_G = 0.0;
double y = (double) (j-side)-d*sample-d;
for (int sx=0;sx<2*sample+1;sx++)
{
y += d;
double x = (double) (i-side)-d*sample-d;
for (int sy=0;sy<2*sample+1;sy++)
{
x += d;
double x_new = x*cos_angle + y*sin_angle;
double y_new = -x*sin_angle + y*cos_angle;
double g = -x_new/(2.0*PI*sigma_x*sigma_y*sigma_x*sigma_x)*exp(-(x_new*x_new/(sigma_x*sigma_x)+y_new*y_new/(sigma_y*sigma_y))/2.0);
sum_G += g;
}
}
temp.ElemNC(j,i) = sum_G/T;
}
}
return temp;
}
// Second Directional Derivative of Gaussian Filter
Matrix SDGaussian2D(int side, float sigma_x, float angle, float ratio)
{
Matrix temp(2*side+1, 2*side+1);
float sigma_y = sigma_x*ratio;
double sin_angle = sin(angle*PI/180);
double cos_angle = cos(angle*PI/180);
int sample;
double d;
if (sigma_x < 2.0 || sigma_y < 2.0)
{
// use denser sample for better filter quality
sample = 5;
d = 1.0/(2.0*sample+1.0);
}
else
{
sample = 0;
d = 1.0;
}
double T = (double) (2*sample+1)*(2*sample+1);
for (int i=0;i<2*side+1;i++)
{
for (int j=0;j<2*side+1;j++)
{
double sum_G = 0.0;
double y = (double) (j-side)-d*sample-d;
for (int sx=0;sx<2*sample+1;sx++)
{
y += d;
double x = (double) (i-side)-d*sample-d;
for (int sy=0;sy<2*sample+1;sy++)
{
x += d;
double x_new = x*cos_angle + y*sin_angle;
double y_new = -x*sin_angle + y*cos_angle;
double g = (x_new*x_new-sigma_x*sigma_x)/(2.0*PI*sigma_x*sigma_y*sigma_x*sigma_x*sigma_x*sigma_x)*exp(-(x_new*x_new/(sigma_x*sigma_x)+y_new*y_new/(sigma_y*sigma_y))/2.0);
sum_G += g;
}
}
temp.ElemNC(j,i) = (float)(sum_G/T);
}
}
return temp;
}
Matrix SDGaussian2D(int side, double sigma_x, double angle, double ratio)
{
Matrix temp(2*side+1, 2*side+1);
double sigma_y = sigma_x*ratio;
double sin_angle = sin(angle*PI/180);
double cos_angle = cos(angle*PI/180);
int sample;
double d;
if (sigma_x < 2.0 || sigma_y < 2.0)
{
// use denser sample for better filter quality
sample = 5;
d = 1.0/(2.0*sample+1.0);
}
else
{
sample = 0;
d = 1.0;
}
double T = (double) (2*sample+1)*(2*sample+1);
for (int i=0;i<2*side+1;i++)
{
for (int j=0;j<2*side+1;j++)
{
double sum_G = 0.0;
double y = (double) (j-side)-d*sample-d;
for (int sx=0;sx<2*sample+1;sx++)
{
y += d;
double x = (double) (i-side)-d*sample-d;
for (int sy=0;sy<2*sample+1;sy++)
{
x += d;
double x_new = x*cos_angle + y*sin_angle;
double y_new = -x*sin_angle + y*cos_angle;
double g = (x_new*x_new-sigma_x*sigma_x)/(2.0*PI*sigma_x*sigma_y*sigma_x*sigma_x*sigma_x*sigma_x)*exp(-(x_new*x_new/(sigma_x*sigma_x)+y_new*y_new/(sigma_y*sigma_y))/2.0);
sum_G += g;
}
}
temp.ElemNC(j,i) = sum_G/T;
}
}
return temp;
}
// Laplacian of Gaussian Filter
Matrix LOG(int side, float sigma_x, float angle, float ratio)
{
Matrix temp(2*side+1, 2*side+1);
float sigma_y = sigma_x*ratio;
double sin_angle = sin(angle*PI/180);
double cos_angle = cos(angle*PI/180);
int sample;
double d;
if (sigma_x < 2.0 || sigma_y < 2.0)
{
// use denser sample for better filter quality
sample = 5;
d = 1.0/(2.0*sample+1.0);
}
else
{
sample = 0;
d = 1.0;
}
double T = (double) (2*sample+1)*(2*sample+1);
for (int i=0;i<2*side+1;i++)
{
for (int j=0;j<2*side+1;j++)
{
double sum_G = 0.0;
double y = (double) (j-side)-d*sample-d;
for (int sx=0;sx<2*sample+1;sx++)
{
y += d;
double x = (double) (i-side)-d*sample-d;
for (int sy=0;sy<2*sample+1;sy++)
{
x += d;
double x_new = x*cos_angle + y*sin_angle;
double y_new = -x*sin_angle + y*cos_angle;
double g = (x_new*x_new+y_new*y_new-sigma_x*sigma_x-sigma_y*sigma_y)
/(2.0*PI*sigma_x*sigma_y*(sigma_x*sigma_x*sigma_x*sigma_x+sigma_y*sigma_y*sigma_y*sigma_y))
*exp(-(x_new*x_new/(sigma_x*sigma_x)+y_new*y_new/(sigma_y*sigma_y))/2.0);
sum_G += g;
}
}
temp.ElemNC(j,i) = (float)(sum_G/T);
}
}
return temp;
}
Matrix LOG(int side, double sigma_x, double angle, double ratio)
{
Matrix temp(2*side+1, 2*side+1);
double sigma_y = sigma_x*ratio;
double sin_angle = sin(angle*PI/180);
double cos_angle = cos(angle*PI/180);
int sample;
double d;
if (sigma_x < 2.0 || sigma_y < 2.0)
{
// use denser sample for better filter quality
sample = 5;
d = 1.0/(2.0*sample+1.0);
}
else
{
sample = 0;
d = 1.0;
}
double T = (double) (2*sample+1)*(2*sample+1);
for (int i=0;i<2*side+1;i++)
{
for (int j=0;j<2*side+1;j++)
{
double sum_G = 0.0;
double y = (double) (j-side)-d*sample-d;
for (int sx=0;sx<2*sample+1;sx++)
{
y += d;
double x = (double) (i-side)-d*sample-d;
for (int sy=0;sy<2*sample+1;sy++)
{
x += d;
double x_new = x*cos_angle + y*sin_angle;
double y_new = -x*sin_angle + y*cos_angle;
double g = (x_new*x_new+y_new*y_new-sigma_x*sigma_x-sigma_y*sigma_y)
/(2.0*PI*sigma_x*sigma_y*(sigma_x*sigma_x*sigma_x*sigma_x+sigma_y*sigma_y*sigma_y*sigma_y))
*exp(-(x_new*x_new/(sigma_x*sigma_x)+y_new*y_new/(sigma_y*sigma_y))/2.0);
sum_G += g;
}
}
temp.ElemNC(j,i) = sum_G/T;
}
}
return temp;
}
// Difference of offset Gaussians filter
Matrix DOOG2D(int side, float sigma_x, float offset, float angle, float ratio)
{
Matrix temp(2*side+1, 2*side+1);
float sigma_y = sigma_x*ratio;
double sin_angle = sin(angle*PI/180);
double cos_angle = cos(angle*PI/180);
int sample;
double d;
if (sigma_x < 2.0 || sigma_y < 2.0)
{
// use denser sample for better filter quality
sample = 5;
d = 1.0/(2.0*sample+1.0);
}
else
{
sample = 0;
d = 1.0;
}
double T = (double) (2*sample+1)*(2*sample+1);
for (int i=0;i<2*side+1;i++)
{
for (int j=0;j<2*side+1;j++)
{
double sum_G = 0.0;
double y = (double) (j-side)-d*sample-d;
for (int sx=0;sx<2*sample+1;sx++)
{
y += d;
double x = (double) (i-side)-d*sample-d;
for (int sy=0;sy<2*sample+1;sy++)
{
x += d;
double x_new = x*cos_angle + y*sin_angle;
double y_new = -x*sin_angle + y*cos_angle;
double g = 1/(2.0*PI*sigma_x*sigma_y)
*exp(-(x_new*x_new/(sigma_x*sigma_x)+y_new*y_new/(sigma_y*sigma_y))/2.0)
- 1/(2.0*PI*sigma_x*sigma_y)
*exp(-((x_new-offset)*(x_new-offset)/(sigma_x*sigma_x)+y_new*y_new/(sigma_y*sigma_y))/2.0);
sum_G += g;
}
}
temp.ElemNC(j,i) = (float)(sum_G/T);
}
}
return temp;
}
Matrix DOOG2D(int side, double sigma_x, double offset, double angle, double ratio)
{
Matrix temp(2*side+1, 2*side+1);
double sigma_y = sigma_x*ratio;
double sin_angle = sin(angle*PI/180);
double cos_angle = cos(angle*PI/180);
int sample;
double d;
if (sigma_x < 2.0 || sigma_y < 2.0)
{
// use denser sample for better filter quality
sample = 5;
d = 1.0/(2.0*sample+1.0);
}
else
{
sample = 0;
d = 1.0;
}
double T = (double) (2*sample+1)*(2*sample+1);
for (int i=0;i<2*side+1;i++)
{
for (int j=0;j<2*side+1;j++)
{
double sum_G = 0.0;
double y = (double) (j-side)-d*sample-d;
for (int sx=0;sx<2*sample+1;sx++)
{
y += d;
double x = (double) (i-side)-d*sample-d;
for (int sy=0;sy<2*sample+1;sy++)
{
x += d;
double x_new = x*cos_angle + y*sin_angle;
double y_new = -x*sin_angle + y*cos_angle;
double g = 1/(2.0*PI*sigma_x*sigma_y)
*exp(-(x_new*x_new/(sigma_x*sigma_x)+y_new*y_new/(sigma_y*sigma_y))/2.0)
- 1/(2.0*PI*sigma_x*sigma_y)
*exp(-((x_new-offset)*(x_new-offset)/(sigma_x*sigma_x)+y_new*y_new/(sigma_y*sigma_y))/2.0);
sum_G += g;
}
}
temp.ElemNC(j,i) = sum_G/T;
}
}
return temp;
}
// Difference of offset Gaussians filter (centered at the middle of two Gaussians)
Matrix DOOG2DCentered(int side, float sigma_x, float offset, float angle, float ratio)
{
Matrix temp(2*side+1, 2*side+1);
float sigma_y = sigma_x*ratio;
double sin_angle = sin(angle*PI/180);
double cos_angle = cos(angle*PI/180);
double half_offset = offset/2.0;
int sample;
double d;
if (sigma_x < 2.0 || sigma_y < 2.0)
{
// use denser sample for better filter quality
sample = 5;
d = 1.0/(2.0*sample+1.0);
}
else
{
sample = 0;
d = 1.0;
}
double T = (double) (2*sample+1)*(2*sample+1);
for (int i=0;i<2*side+1;i++)
{
for (int j=0;j<2*side+1;j++)
{
double sum_G = 0.0;
double y = (double) (j-side)-d*sample-d;
for (int sx=0;sx<2*sample+1;sx++)
{
y += d;
double x = (double) (i-side)-d*sample-d;
for (int sy=0;sy<2*sample+1;sy++)
{
x += d;
double x_new = x*cos_angle + y*sin_angle;
double y_new = -x*sin_angle + y*cos_angle;
double g = 1/(2.0*PI*sigma_x*sigma_y)
*exp(-((x_new+half_offset)*(x_new+half_offset)/(sigma_x*sigma_x)+y_new*y_new/(sigma_y*sigma_y))/2.0)
- 1/(2.0*PI*sigma_x*sigma_y)
*exp(-((x_new-half_offset)*(x_new-half_offset)/(sigma_x*sigma_x)+y_new*y_new/(sigma_y*sigma_y))/2.0);
sum_G += g;
}
}
temp.ElemNC(j,i) = (float)(sum_G/T);
}
}
return temp;
}
Matrix DOOG2DCentered(int side, double sigma_x, double offset, double angle, double ratio)
{
Matrix temp(2*side+1, 2*side+1);
double sigma_y = sigma_x*ratio;
double sin_angle = sin(angle*PI/180);
double cos_angle = cos(angle*PI/180);
double half_offset = offset/2.0;
int sample;
double d;
if (sigma_x < 2.0 || sigma_y < 2.0)
{
// use denser sample for better filter quality
sample = 5;
d = 1.0/(2.0*sample+1.0);
}
else
{
sample = 0;
d = 1.0;
}
double T = (double) (2*sample+1)*(2*sample+1);
for (int i=0;i<2*side+1;i++)
{
for (int j=0;j<2*side+1;j++)
{
double sum_G = 0.0;
double y = (double) (j-side)-d*sample-d;
for (int sx=0;sx<2*sample+1;sx++)
{
y += d;
double x = (double) (i-side)-d*sample-d;
for (int sy=0;sy<2*sample+1;sy++)
{
x += d;
double x_new = x*cos_angle + y*sin_angle;
double y_new = -x*sin_angle + y*cos_angle;
double g = 1/(2.0*PI*sigma_x*sigma_y)
*exp(-((x_new+half_offset)*(x_new+half_offset)/(sigma_x*sigma_x)+y_new*y_new/(sigma_y*sigma_y))/2.0)
- 1/(2.0*PI*sigma_x*sigma_y)
*exp(-((x_new-half_offset)*(x_new-half_offset)/(sigma_x*sigma_x)+y_new*y_new/(sigma_y*sigma_y))/2.0);
sum_G += g;
}
}
temp.ElemNC(j,i) = sum_G/T;
}
}
return temp;
}
// Filter image with these filters
Matrix FilterGaussian2D(Matrix& image, float sigma_x, float angle, float ratio)
{
float sigma_y = ratio*sigma_x;
float sigma_max = sigma_x>=sigma_y?sigma_x:sigma_y;
int side = (int)(3*sigma_max+0.5);
Matrix filt = Gaussian2D(side, sigma_x, angle, ratio);
return Conv2(image, filt);
}
Matrix FilterGaussian2D(Matrix& image, double sigma_x, double angle, double ratio)
{
double sigma_y = ratio*sigma_x;
double sigma_max = sigma_x>=sigma_y?sigma_x:sigma_y;
int side = (int)(3*sigma_max+0.5);
Matrix filt = Gaussian2D(side, sigma_x, angle, ratio);
return Conv2(image, filt);
}
// First derivative of Gaussian
Matrix FilterFDGaussian2D(Matrix& image, float sigma_x, float angle, float ratio)
{
float sigma_y = ratio*sigma_x;
float sigma_max = sigma_x>=sigma_y?sigma_x:sigma_y;
int side = (int)(3*sigma_max+0.5);
Matrix filt = FDGaussian2D(side, sigma_x, angle, ratio);
return Conv2(image, filt);
}
Matrix FilterFDGaussian2D(Matrix& image, double sigma_x, double angle, double ratio)
{
double sigma_y = ratio*sigma_x;
double sigma_max = sigma_x>=sigma_y?sigma_x:sigma_y;
int side = (int)(3*sigma_max+0.5);
Matrix filt = FDGaussian2D(side, sigma_x, angle, ratio);
return Conv2(image, filt);
}
// Second derivative of Gaussian
Matrix FilterSDGaussian2D(Matrix& image, float sigma_x, float angle, float ratio)
{
float sigma_y = ratio*sigma_x;
float sigma_max = sigma_x>=sigma_y?sigma_x:sigma_y;
int side = (int)(3*sigma_max+0.5);
Matrix filt = SDGaussian2D(side, sigma_x, angle, ratio);
return Conv2(image, filt);
}
Matrix FilterSDGaussian2D(Matrix& image, double sigma_x, double angle, double ratio)
{
double sigma_y = ratio*sigma_x;
double sigma_max = sigma_x>=sigma_y?sigma_x:sigma_y;
int side = (int)(3*sigma_max+0.5);
Matrix filt = SDGaussian2D(side, sigma_x, angle, ratio);
return Conv2(image, filt);
}
// Laplacian of Gaussian
Matrix FilterLOG(Matrix& image, float sigma_x, float angle, float ratio)
{
float sigma_y = ratio*sigma_x;
float sigma_max = sigma_x>=sigma_y?sigma_x:sigma_y;
int side = (int)(3*sigma_max+0.5);
Matrix filt = LOG(side, sigma_x, angle, ratio);
return Conv2(image, filt);
}
Matrix FilterLOG(Matrix& image, double sigma_x, double angle, double ratio)
{
double sigma_y = ratio*sigma_x;
double sigma_max = sigma_x>=sigma_y?sigma_x:sigma_y;
int side = (int)(3*sigma_max+0.5);
Matrix filt = LOG(side, sigma_x, angle, ratio);
return Conv2(image, filt);
}
// Difference of offset Gaussians
Matrix FilterDOOG2D(Matrix& image, float sigma_x, float offset, float angle, float ratio)
{
float sigma_y = ratio*sigma_x;
int side = (int)((3*sigma_x+offset>=3*sigma_y?3*sigma_x+offset:3*sigma_y)+0.5);
Matrix filt = DOOG2D(side, sigma_x, offset, angle, ratio);
return Conv2(image, filt);
}
Matrix FilterDOOG2D(Matrix& image, double sigma_x, double offset, double angle, double ratio)
{
double sigma_y = ratio*sigma_x;
int side = (int)((3*sigma_x+offset>=3*sigma_y?3*sigma_x+offset:3*sigma_y)+0.5);
Matrix filt = DOOG2D(side, sigma_x, offset, angle, ratio);
return Conv2(image, filt);
}
Matrix FilterDOOG2DCentered(Matrix& image, float sigma_x, float offset, float angle, float ratio)
{
float sigma_y = ratio*sigma_x;
int side = (int)((3*sigma_x+offset/2>=3*sigma_y?3*sigma_x+offset/2:3*sigma_y)+0.5);
Matrix filt = DOOG2DCentered(side, sigma_x, offset, angle, ratio);
return Conv2(image, filt);
}
Matrix FilterDOOG2DCentered(Matrix& image, double sigma_x, double offset, double angle, double ratio)
{
double sigma_y = ratio*sigma_x;
int side = (int)((3*sigma_x+offset/2>=3*sigma_y?3*sigma_x+offset/2:3*sigma_y)+0.5);
Matrix filt = DOOG2DCentered(side, sigma_x, offset, angle, ratio);
return Conv2(image, filt);
}
// Various Edge Detection schemes
Matrix NonMaximaSuppress(Matrix& edgesMain, Matrix& theta)
{
Matrix vectorX = Cos(theta);
Matrix vectorY = Sin(theta);
return NonMaximaSuppress(edgesMain, vectorX, vectorY);
}
Matrix NonMaximaSuppress(Matrix& edgesMain, Matrix& theta)
{
Matrix vectorX = Cos(theta);
Matrix vectorY = Sin(theta);
return NonMaximaSuppress(edgesMain, vectorX, vectorY);
}
Matrix NonMaximaSuppress(Matrix& edgesMain, Matrix& vectorX, Matrix& vectorY)
{
Matrix edges = edgesMain.Clone();
Matrix mask = NonMaximaMask(edges - Minimum(edges(":")), vectorX, vectorY);
for(int i=0;i NonMaximaSuppress(Matrix& edgesMain, Matrix& vectorX, Matrix& vectorY)
{
Matrix edges = edgesMain.Clone();
Matrix mask = NonMaximaMask(edges - Minimum(edges(":")), vectorX, vectorY);
for(int i=0;i NonMaximaMask(Matrix& edges, Matrix& vectorX, Matrix& vectorY)
{
Matrix mask(edges.Rows(), edges.Columns(), 1);
Matrix theta(edges.Rows(), edges.Columns(), 0);
Matrix mask15(edges.Rows(), edges.Columns(), 0);
Matrix mask26(edges.Rows(), edges.Columns(), 0);
Matrix mask37(edges.Rows(), edges.Columns(), 0);
Matrix mask48(edges.Rows(), edges.Columns(), 0);
for(int i=0;i= 0 && theta(j,i) < PI/4 )
{
mask15(j,i) = 1;
}
else if( theta(j,i) >= PI/4 && theta(j,i) < PI/2 )
{
mask26(j,i) = 1;
}
else if( theta(j,i) >= PI/2 && theta(j,i) < PI*3/4 )
{
mask37(j,i) = 1;
}
else if( theta(j,i) >= PI*3/4 && theta(j,i) < PI )
{
mask48(j,i) = 1;
}
}
}
//Case 1
for(int i=0;i NonMaximaMask(Matrix& edges, Matrix& vectorX, Matrix& vectorY)
{
Matrix mask(edges.Rows(), edges.Columns(), 1);
Matrix theta(edges.Rows(), edges.Columns(), 0);
Matrix mask15(edges.Rows(), edges.Columns(), 0);
Matrix mask26(edges.Rows(), edges.Columns(), 0);
Matrix mask37(edges.Rows(), edges.Columns(), 0);
Matrix mask48(edges.Rows(), edges.Columns(), 0);
for(int i=0;i= 0 && theta(j,i) < PI/4 )
{
mask15(j,i) = 1;
}
else if( theta(j,i) >= PI/4 && theta(j,i) < PI/2 )
{
mask26(j,i) = 1;
}
else if( theta(j,i) >= PI/2 && theta(j,i) < PI*3/4 )
{
mask37(j,i) = 1;
}
else if( theta(j,i) >= PI*3/4 && theta(j,i) < PI )
{
mask48(j,i) = 1;
}
}
}
//Case 1
for(int i=0;i=0 && y>=0)
{
return atan2(y,x);
}
else if(x<0 && y == 0)
{
return atan2(y,x) - PI;
}
else if(x<0 && y>0)
{
return atan2(y,x);
}
else if(x<=0 && y<0)
{
return PI + atan2(y,x) ;
}
else if(x>0 && y<0)
{
return PI + atan2(y,x) ;
}
else
{
cout << "Something went wrong with Direction(x,y) function.. Returning -1..\n" << endl;
return -1;
}
}
float Direction(float y, float x)
{
if(x>=0 && y>=0)
{
return atan2(y,x);
}
else if(x<0 && y == 0)
{
return atan2(y,x) - PI;
}
else if(x<0 && y>0)
{
return atan2(y,x);
}
else if(x<=0 && y<0)
{
return PI + atan2(y,x) ;
}
else if(x>0 && y<0)
{
return PI + atan2(y,x) ;
}
else
{
cout << "Something went wrong with Direction(x,y) function.. Returning -1..\n" << endl;
return -1;
}
}
// Edgeflow
// both grayscale and multi-valued
// Outout is two dimensional (x and y components of the vector field)
MatrixList EdgeflowVectorField(Matrix& image, int angles, float sigma, float offset, float ratio, bool normalized)
{
Vector radAngles(2*angles);
Vector cosAngles(2*angles);
Vector sinAngles(2*angles);
for(int i=0; i > energies(2*angles);
float sigma_y = ratio*sigma;
int side = (int)((3*sigma+offset>=3*sigma_y?3*sigma+offset:3*sigma_y)+0.5);
for(int i=0; i filt = DOOG2D(side, sigma, offset, 180+angle, ratio);
energies(i) = AbsI(Conv2(image, filt));
//sprintf(s, "%02d.0.bmp", i);
//ImWrite(energies(i),s);
filt = DOOG2D(side, sigma, offset, angle, ratio);
energies(i+angles) = AbsI(Conv2(image, filt));
//sprintf(s, "%02d.1.bmp", i);
//ImWrite(energies(i+angles),s);
//Matrix filt = DOOG2DCentered(side, sigma, offset, angle, ratio);
//char s[200];
//Matrix filtOut = Conv2(image, filt, Full);
//
//int offsetX = (int)(offset*cos(radAngles[i])/2+0.5);
//offsetX = offsetX==0?1:offsetX;
//int offsetY = (int)(offset*sin(radAngles[i])/2+0.5);
//offsetY = offsetY==0?1:offsetY;
//// cout << (fmod(2.0+cos((double)radAngles[i]),2.0)-1) << " : " << (fmod(2.0+sin((double)radAngles[i]),2.0)-1) << endl;
//energies(i) = Abs(filtOut.Slice(side+offsetY, image.Rows()+side+offsetY-1, side+offsetX, image.Columns()+side+offsetX-1));
//sprintf(s, "%02d.0.bmp", i);
//ImWrite(energies(i), s);
//
//energies(i+angles) = Abs(filtOut.Slice(side-offsetY, image.Rows()+side-offsetY-1, side-offsetX, image.Columns()+side-offsetX-1));
//sprintf(s, "%02d.1.bmp", i);
//ImWrite(energies(i+angles), s);
}
//pf.Toc();
// Sum energies over half circles
//pf.Tic();
Vector > sumEnergies(2*angles);
for(int i=0; i<2*angles; i++)
{
sumEnergies(i) = Matrix(image.Rows(), image.Columns(), 0.0f);
int start = i-angles/2;
int end = start+angles;
for(int j=start; j > probabilities(2*angles);
for(int i=0; i(image.Rows(), image.Columns());
probabilities(i+angles) = Matrix(image.Rows(), image.Columns());
Matrix total = sumEnergies(i) + sumEnergies(i+angles);
for(int r=0;r 1e-9)
{
probabilities(i).Elem(r,c) = sumEnergies(i).Elem(r,c)/total(r,c);
probabilities(i+angles).Elem(r,c) = sumEnergies(i+angles).Elem(r,c)/total(r,c);
}
else
{
probabilities(i).Elem(r,c) = 0.5;
probabilities(i+angles).Elem(r,c) = 0.5;
}
}
}
}
//pf.Toc();
//pf.Tic();
Matrix xFlow(image.Rows(), image.Columns());
Matrix yFlow(image.Rows(), image.Columns());
Matrix maxEnergy(image.Rows(), image.Columns());
for(int r=0;r=vl && v>=vr)
{
float orientation, strength;
if(vl+vr != 2*v)
{
orientation = 0.5f*(vl - vr)/(vl + vr - 2*v);
strength = v + 0.5f*( (vr - vl)*orientation
+ (vr + vl - 2*v)*orientation*orientation );
orientation = (float)((k+orientation)*PI/angles);
}
else
{
strength = v;
orientation = (float)(k*PI/angles);
}
dirX += cos(orientation)*strength;
dirY += sin(orientation)*strength;
//maxEnergy(r,c) += strength;
//count++;
}
}
//maxEnergy(r,c) /= (float)count;
float dir = sqrt(dirX*dirX+dirY*dirY);
dirX /= dir;
dirY /= dir;
float value;
int ang;
if(probabilities(0).Elem(r,c)>probabilities(angles).Elem(r,c))
{
value = probabilities(0).Elem(r,c);
ang = 0;
}
else
{
value = probabilities(angles).Elem(r,c);
ang = angles;
}
float maximum = value;
int maxIndex = ang;
float minimum = value;
int minIndex = ang;
int maxOrientation = 0;
int minOrientation = 0;
for(int k=1;kprobabilities(k+angles).Elem(r,c))
{
value = probabilities(k).Elem(r,c);
ang = k;
}
else
{
value = probabilities(k+angles).Elem(r,c);
ang = k+angles;
}
if(value > maximum)
{
maximum = value;
maxIndex = ang;
maxOrientation = k;
}
if(value < minimum)
{
minimum = value;
minIndex = ang;
minOrientation = k;
}
}
if(normalized)
{
xFlow(r,c) = dirX*probabilities(maxIndex)(r,c);
yFlow(r,c) = dirY*probabilities(maxIndex)(r,c);
}
else
{
xFlow(r,c) = dirX*sumEnergies(maxIndex)(r,c);
yFlow(r,c) = dirY*sumEnergies(maxIndex)(r,c);
}
maxEnergy(r,c) = sumEnergies(maxIndex)(r,c);
}
}
//if(normalized)
//{
// xFlow *= ToFloat(maxEnergy > (float)angles);
// yFlow *= ToFloat(maxEnergy > (float)angles);
//}
//ImWrite(xFlow, "xflow.bmp");
//ImWrite(yFlow, "yflow.bmp");
//ImWrite(maxEnergy, "maxEnergy.bmp");
return MatrixList(xFlow, yFlow);
}
//MatrixList EdgeflowVectorField(Matrix image, double sigma, double offset, double ratio, int angles)
//{
//}
MatrixList EdgeflowVectorField(MatrixList& image, int angles, float sigma, float offset, float ratio, bool normalized)
{
Vector radAngles(2*angles);
Vector cosAngles(2*angles);
Vector sinAngles(2*angles);
for(int i=0; i > energies(2*angles);
for(int i=0;i(image.Rows(), image.Columns(), 0);
}
float sigma_y = ratio*sigma;
int side = (int)((3*sigma+offset>=3*sigma_y?3*sigma+offset:3*sigma_y)+0.5);
for(int i=0; i filt = DOOG2D(side, sigma, offset, 180+angle, ratio);
energies(i) += AbsI(Conv2(image[c], filt));
//sprintf(s, "%02d.0.bmp", i);
//ImWrite(energies(i),s);
filt = DOOG2D(side, sigma, offset, angle, ratio);
energies(i+angles) += AbsI(Conv2(image[c], filt));
//sprintf(s, "%02d.1.bmp", i);
//ImWrite(energies(i+angles),s);
}
//Matrix filt = DOOG2DCentered(side, sigma, offset, angle, ratio);
//char s[200];
//Matrix filtOut = Conv2(image, filt, Full);
//
//int offsetX = (int)(offset*cos(radAngles[i])/2+0.5);
//offsetX = offsetX==0?1:offsetX;
//int offsetY = (int)(offset*sin(radAngles[i])/2+0.5);
//offsetY = offsetY==0?1:offsetY;
//// cout << (fmod(2.0+cos((double)radAngles[i]),2.0)-1) << " : " << (fmod(2.0+sin((double)radAngles[i]),2.0)-1) << endl;
//energies(i) = Abs(filtOut.Slice(side+offsetY, image.Rows()+side+offsetY-1, side+offsetX, image.Columns()+side+offsetX-1));
//sprintf(s, "%02d.0.bmp", i);
//ImWrite(energies(i), s);
//
//energies(i+angles) = Abs(filtOut.Slice(side-offsetY, image.Rows()+side-offsetY-1, side-offsetX, image.Columns()+side-offsetX-1));
//sprintf(s, "%02d.1.bmp", i);
//ImWrite(energies(i+angles), s);
}
//pf.Toc();
// Sum energies over half circles
//pf.Tic();
Vector > sumEnergies(2*angles);
for(int i=0; i<2*angles; i++)
{
sumEnergies(i) = Matrix(image.Rows(), image.Columns(), 0.0f);
int start = i-angles/2;
int end = start+angles;
for(int j=start; j > probabilities(2*angles);
for(int i=0; i(image.Rows(), image.Columns());
probabilities(i+angles) = Matrix(image.Rows(), image.Columns());
Matrix total = sumEnergies(i) + sumEnergies(i+angles);
for(int r=0;r 1e-9)
{
probabilities(i).Elem(r,c) = sumEnergies(i).Elem(r,c)/total(r,c);
probabilities(i+angles).Elem(r,c) = sumEnergies(i+angles).Elem(r,c)/total(r,c);
}
else
{
probabilities(i).Elem(r,c) = 0.5;
probabilities(i+angles).Elem(r,c) = 0.5;
}
}
}
}
//pf.Toc();
//pf.Tic();
Matrix xFlow(image.Rows(), image.Columns());
Matrix yFlow(image.Rows(), image.Columns());
Matrix maxEnergy(image.Rows(), image.Columns());
for(int r=0;r=vl && v>=vr)
{
float orientation, strength;
if(vl+vr != 2*v)
{
orientation = 0.5f*(vl - vr)/(vl + vr - 2*v);
strength = v + 0.5f*( (vr - vl)*orientation
+ (vr + vl - 2*v)*orientation*orientation );
orientation = (float)((k+orientation)*PI/angles);
}
else
{
strength = v;
orientation = (float)(k*PI/angles);
}
dirX += cos(orientation)*strength;
dirY += sin(orientation)*strength;
//maxEnergy(r,c) += strength;
//count++;
}
}
//maxEnergy(r,c) /= (float)count;
float dir = sqrt(dirX*dirX+dirY*dirY);
dirX /= dir;
dirY /= dir;
float value;
int ang;
if(probabilities(0).Elem(r,c)>probabilities(angles).Elem(r,c))
{
value = probabilities(0).Elem(r,c);
ang = 0;
}
else
{
value = probabilities(angles).Elem(r,c);
ang = angles;
}
float maximum = value;
int maxIndex = ang;
float minimum = value;
int minIndex = ang;
int maxOrientation = 0;
int minOrientation = 0;
for(int k=1;kprobabilities(k+angles).Elem(r,c))
{
value = probabilities(k).Elem(r,c);
ang = k;
}
else
{
value = probabilities(k+angles).Elem(r,c);
ang = k+angles;
}
if(value > maximum)
{
maximum = value;
maxIndex = ang;
maxOrientation = k;
}
if(value < minimum)
{
minimum = value;
minIndex = ang;
minOrientation = k;
}
}
if(normalized)
{
xFlow(r,c) = dirX*probabilities(maxIndex)(r,c);
yFlow(r,c) = dirY*probabilities(maxIndex)(r,c);
}
else
{
xFlow(r,c) = dirX*sumEnergies(maxIndex)(r,c);
yFlow(r,c) = dirY*sumEnergies(maxIndex)(r,c);
}
maxEnergy(r,c) = sumEnergies(maxIndex)(r,c);
}
}
//if(normalized)
//{
// xFlow *= ToFloat(maxEnergy > (float)angles);
// yFlow *= ToFloat(maxEnergy > (float)angles);
//}
//ImWrite(xFlow, "xflow.bmp");
//ImWrite(yFlow, "yflow.bmp");
//ImWrite(maxEnergy, "maxEnergy.bmp");
return MatrixList(xFlow, yFlow);
}
void block_write(Matrix& matrix, int x_ind, int y_ind, int *keys, int energy);
double my_atan2(double y, double x);
Matrix CreateFlowImage(Matrix& xFlow, Matrix& yFlow)
{
if(xFlow.Rows() != yFlow.Rows() || xFlow.Columns() != yFlow.Columns())
{
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Matrix dimensions does not match to each other!");
}
int display_map[8][15] ={37,38,39,40,41,42,43,33,23,13,12,51,59,67,66,
10,20,30,40,50,60,70,61,52,43,34,69,68,67,66,
13,22,31,40,49,58,67,59,51,43,34,57,47,37,28,
16,24,32,40,48,56,64,55,46,37,28,65,66,67,68,
43,42,41,40,39,38,37,47,57,67,68,29,21,13,14,
70,60,50,40,30,20,10,19,28,37,46,11,12,13,14,
67,58,49,40,31,22,13,21,29,37,46,23,33,43,52,
64,56,48,40,32,24,16,25,34,43,52,15,14,13,12};
int quant[17] = {0, 45, 45, 90, 90, 135, 135, 180, 180, 225, 225, 270, 270, 315, 315, 0, 0};
Matrix angles(xFlow.Rows(), xFlow.Columns());
Matrix raw_angles(xFlow.Rows(), xFlow.Columns());
Matrix flow(9*xFlow.Rows(), 9*xFlow.Columns(), 255);
Matrix energies = SqrtI(xFlow*xFlow + yFlow*yFlow);
float maximum = Maximum(energies(":"));
float minimum = Minimum(energies(":"));
Matrix scaled_energies(xFlow.Rows(), xFlow.Columns(), 0);
if(maximum-minimum > 1e-9)
{
scaled_energies = 255 - ToInt((energies-minimum)*255.0f/(maximum-minimum));
//scaled_energies = 255 - ToInt(energies*255.0f/maximum);
}
int disp_keys[15];
for(int i=0;i& matrix, int x_ind, int y_ind, int *keys, int energy){
for(int t=0;t<15;t++){
int x_loc = 9*x_ind + (keys[t] % 9);
int y_loc = 9*y_ind + (keys[t] / 9);
matrix(y_loc,x_loc) = energy;
}
}
double my_atan2(double y, double x){
if(x>=0 && y>=0){
return atan2(y,x)*180.0/PI;
}
else if(x>=0 && y<=0){
return ( 2*PI + atan2(y,x) )*180.0/PI;
}
else if(x<=0 && y>=0){
return atan2(y,x)*180.0/PI;
}
else if(x<=0 && y<=0){
return (2*PI + atan2(y,x) )*180.0/PI;
}
else
{
return -1;
}
}
Matrix HystThreshold(Matrix& m, float thHigh, float thLow)
{
Matrix th = (m > thHigh);
Matrix th2 = (m > thLow);
int count;
do{
count = 0;
for(int i=0;i0 && th.ElemNC(i-1,j)==1)
{
th.ElemNC(i,j) = 1;
count++;
}
if(i0 && th.ElemNC(i,j-1)==1)
{
th.ElemNC(i,j) = 1;
count++;
}
if(j0 && j>0 && th.ElemNC(i-1,j-1)==1)
{
th.ElemNC(i,j) = 1;
count++;
}
if(i0 && th.ElemNC(i+1,j-1)==1)
{
th.ElemNC(i,j) = 1;
count++;
}
if(i> 0 && j 0);
return th;
}
Matrix& PMAnisoDiff(Matrix& image, float K, int iterations)
{
float lambda = 0.25;
float K2Inv = 1/(K*K);
Matrix gradX(image.Rows(), image.Columns());
Matrix gradY(image.Rows(), image.Columns());
gradX(0, image.Rows()-1, 0, 0) = 0;
gradY(0, 0, 0, image.Columns()-1) = 0;
for(int k=0; k fluxX = gradX*Exp(-K2Inv*Abs(gradX));
Matrix fluxY = gradY*Exp(-K2Inv*Abs(gradY));
// Update image
for(int i=0;i& PMAnisoDiff(Matrix& image, double K, int iterations)
{
double lambda = 0.25;
double K2Inv = 1/(K*K);
Matrix gradX(image.Rows(), image.Columns());
Matrix gradY(image.Rows(), image.Columns());
gradX(0, image.Rows()-1, 0, 0) = 0;
gradY(0, 0, 0, image.Columns()-1) = 0;
for(int k=0; k fluxX = gradX*Exp(-K2Inv*Abs(gradX));
Matrix fluxY = gradY*Exp(-K2Inv*Abs(gradY));
// Update image
for(int i=0;i& PMAnisoDiff(MatrixList& image, float K, int iterations)
//{
//}
//MatrixList& PMAnisoDiff(MatrixList& image, double K, int iterations)
//{
//}
class CompareGrads
{
public:
CompareGrads(){}
~CompareGrads(){}
static Matrix gradMatrix;
static int Compare( Point px, Point py )
{
if(gradMatrix(px.y,px.x) > gradMatrix(py.y,py.x))
{
return 1;
}
else if(gradMatrix(px.y,px.x) == gradMatrix(py.y,py.x))
{
return 0;
}
else
{
return -1;
}
}
};
Matrix CompareGrads::gradMatrix;
Matrix GetEgdes(Matrix &grads, Matrix &thick, Matrix &suppressed)
{
Matrix edges(grads.Rows(), grads.Columns(), 0);
Vector LabelCount(grads.Rows()*grads.Columns(), 0);
Vector LabelAvg(grads.Rows()*grads.Columns(), 0.0f);
Matrix labels(grads.Rows(), grads.Columns(), 0);
Matrix dummy(grads.Rows(), grads.Columns(), 0);
int label = 0;
// Point struct needs to be defined.
queue outsQ;
queue tmpQ;
for (int x=0; x 0)
{
Point p = tmpQ.front();
tmpQ.pop();
for (int i=-1; i<=1; i++)
{
for (int j=-1; j<=1; j++)
{
if(i == 0 || j == 0)
{
if(p.x+i>=0 && p.x+i=0 && p.y+j tmp(grads.Rows(), grads.Columns(), 0);
//Hashtable Marked = new Hashtable();
bool greedy = false;
int counter = 0;
while(outsQ.size() > 0)
{
counter++;
//Marked.Clear();
int len = outsQ.size();
//cout << "Length: " << len << endl;
// copy queue to a list
vector psA;
for(int i=0; i A;
for (int i=-1; i<=1; i++)
{
for (int j=-1; j<=1; j++)
{
if(i == 0 || j == 0)
{
if(p.x+i>=0 && p.x+i=0 && p.y+j1?1:ac;
ac = ac<-1?-1:ac;
return acos(ac)*180/PI;
}
//The following code for RGB to Lab conversion is adapted from code
//written by Yossi Rubner and Mark Ruzon. Following comments belong to them.
//Baris Sumengen 09/09/2005
//
// Convert between RGB and CIE-Lab color spaces
// Uses ITU-R recommendation BT.709 with D65 as reference white.
// Yossi Rubner 2/24/98
//
// Mark Ruzon 3/18/99 -- Added thresholds to truncate parts of the space
// where changing a value has no visible effect on the color.
void RGB2Lab(double R, double G, double B, double &L, double &a, double &b)
{
const double BLACK = 20;
const double YELLOW = 70;
double X, Y, Z, fX, fY, fZ;
X = 0.412453*R + 0.357580*G + 0.180423*B;
Y = 0.212671*R + 0.715160*G + 0.072169*B;
Z = 0.019334*R + 0.119193*G + 0.950227*B;
X /= (255 * 0.950456);
Y /= 255;
Z /= (255 * 1.088754);
if(Y > 0.008856)
{
fY = pow(Y, 1.0/3.0);
L = 116.0*fY - 16.0;
}
else
{
fY = 7.787*Y + 16.0/116.0;
L = 903.3*Y;
}
if(X > 0.008856)
{
fX = pow(X, 1.0/3.0);
}
else
{
fX = 7.787*X + 16.0/116.0;
}
if(Z > 0.008856)
{
fZ = pow(Z, 1.0/3.0);
}
else
{
fZ = 7.787*Z + 16.0/116.0;
}
a = 500.0*(fX - fY);
b = 200.0*(fY - fZ);
if(L < BLACK)
{
a *= exp((L - BLACK) / (BLACK / 4));
b *= exp((L - BLACK) / (BLACK / 4));
L = BLACK;
}
if(b > YELLOW)
{
b = YELLOW;
}
}
void Lab2RGB(double L, double a, double b, double &R, double &G, double &B)
{
const double BLACK = 20;
const double YELLOW = 70;
double X, Y, Z, fX, fY, fZ;
double RR, GG, BB;
fY = pow((L + 16.0) / 116.0, 3.0);
if(fY < 0.008856)
{
fY = L / 903.3;
}
Y = fY;
if(fY > 0.008856)
{
fY = pow(fY, 1.0/3.0);
}
else
{
fY = 7.787 * fY + 16.0/116.0;
}
fX = a / 500.0 + fY;
if(fX > 0.206893)
{
X = pow(fX, 3.0);
}
else
{
X = (fX - 16.0/116.0) / 7.787;
}
fZ = fY - b /200.0;
if(fZ > 0.206893)
{
Z = pow(fZ, 3.0);
}
else
{
Z = (fZ - 16.0/116.0) / 7.787;
}
X *= (0.950456 * 255);
Y *= 255;
Z *= (1.088754 * 255);
RR = 3.240479*X - 1.537150*Y - 0.498535*Z;
GG = -0.969256*X + 1.875992*Y + 0.041556*Z;
BB = 0.055648*X - 0.204043*Y + 1.057311*Z;
R = (double)(RR < 0 ? 0 : RR > 255 ? 255 : RR);
G = (double)(GG < 0 ? 0 : GG > 255 ? 255 : GG);
B = (double)(BB < 0 ? 0 : BB > 255 ? 255 : BB);
}
MatrixList RGB2Lab(MatrixList input)
{
if(input.Planes() != 3)
{
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Length of MatrixList should be 3 for an RGB image!");
}
MatrixList tmp(3, input.Rows(),input.Columns());
for(int i=0;i Lab2RGB(MatrixList input)
{
if(input.Planes() != 3)
{
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("Length of MatrixList should be 3 for an Lab image!");
}
MatrixList tmp(3, input.Rows(),input.Columns());
for(int i=0;i SegmentEF(Matrix &im, bool normalized, float initScale, float scaleJump,
float endScale, float angleLimit, float ratioLimit, float smoothWeighting, float stopError, int accuracy)
{
float diag = sqrt(float(im.Rows()*im.Rows() + im.Columns()*im.Columns()));
float atomic = diag/1000.0;
cout << "Unit scale: " << atomic << " pixels" << endl << endl;
float scale = initScale*atomic;
cout << "Flow field at scale: " << scale << endl;
MatrixList flow = EdgeflowVectorField(im, 8, scale, scale, 1.0f, normalized);
endScale = endScale*atomic;
scaleJump = scaleJump*atomic;
scale += scaleJump;
while( scale <= endScale)
{
Matrix efMag = Sqrt(flow(0)*flow(0)+flow(1)*flow(1));
float magLimit = Maximum(efMag(":"))/ratioLimit;
cout << "Flow field at scale: " << scale << " pixels" << endl;
MatrixList tempflow = EdgeflowVectorField(im, 4, scale, scale, 1.0f, normalized);
Matrix tempMag = Sqrt(tempflow(0)*tempflow(0)+tempflow(1)*tempflow(1));
for (int i=0; i BB = Divergence(flow(0), flow(1));
float *tmp = poisson((-BB).Data(), BB.Columns(), BB.Rows());
Matrix CC(BB.Rows(), BB.Columns());
for (int j=0; j b(im.Rows()+6, im.Columns()+6,0);
b(3,im.Rows()+2,3,im.Columns()+2) = smoothWeighting*CC;
Matrix u(im.Rows()+6, im.Columns()+6,0);
u(3,im.Rows()+2,3,im.Columns()+2) = flow(0);
Matrix v(im.Rows()+6, im.Columns()+6,0);
v(3,im.Rows()+2,3,im.Columns()+2) = flow(1);
Matrix phi(im.Rows()+6, im.Columns()+6);
phi(3,im.Rows()+2,3,im.Columns()+2) = im;
Matrix delta(im.Rows()+6, im.Columns()+6);
Matrix delta2(im.Rows()+6, im.Columns()+6);
Matrix old = im.Clone();
PerfTimer pt200;
pt200.Tic();
float oldChange = 1e10;
int k=0;
while(1)
{
ExtendConst2D(phi,3);
Evolve2DKappa(phi, 1, 1, 1, 1, b, delta);
switch( accuracy )
{
case 1:
Evolve2DVectorENO1(phi, 1, 1, u, v, delta2);
break;
case 2:
Evolve2DVectorENO2(phi, 1, 1, u, v, delta2);
break;
case 3:
Evolve2DVectorENO3(phi, 1, 1, u, v, delta2);
break;
case 5:
Evolve2DVectorWENO(phi, 1, 1, u, v, delta2);
break;
default:
cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
Utility::RunTimeError("accuracy parameter can only be 1, 2, 3, or 5!");
}
float dt = Getdt2DVectorKappa(0.95f, 1, 1, 1, 1, u, v, b);
AddPhi2D(phi,delta,dt);
SubtractPhi2D(phi,delta2,dt);
if(k%200 == 199 ){
cout << "EF: " << k+1 << endl;
pt200.Toc();
char str[100];
sprintf(str, "phi_%03d.bmp", k);
for(int i=0; i 255)
{
phi(i) = 255;
}
}
Matrix diffused = phi.Slice(3,im.Rows()+2,3,im.Columns()+2);
// ImWrite(diffused, str);
float change = Sum(Vector(Abs(old-diffused)))/im.Numel();
cout << "Error: " << change << endl << endl;
if(change < stopError || k > 10000)
{
break;
}
if(change > oldChange)
{
Utility::Warning("!!!!! Instaability detected during diffusion... Exiting diffusion stage!");
break;
}
oldChange = change;
old = diffused.Clone();
pt200.Tic();
}
k++;
}
phi = phi.Slice(3,im.Rows()+2,3,im.Columns()+2);
for(int i=0; i 255)
{
phi(i) = 255;
}
}
ImWrite(phi, "diffused.png");
Matrix gx = FilterFDGaussian2D(phi, atomic, 0);
Matrix gy = FilterFDGaussian2D(phi, atomic, 90);
Matrix gradIm = Sqrt(SQR(gx) + SQR(gy));
Matrix gradSupp = NonMaximaSuppress(gradIm, gx, gy);
Matrix edges = GetEgdes(gradIm, gradIm > 1, gradSupp > 1);
return edges;
}
Matrix SegmentEF(MatrixList &im, bool normalized, float initScale, float scaleJump,
float endScale, float angleLimit, float ratioLimit, float smoothWeighting, float stopError, int accuracy)
{
float diag = sqrt(float(im.Rows()*im.Rows() + im.Columns()*im.Columns()));
float atomic = diag/1000;
cout << "Unit scale: " << atomic << " pixels" << endl << endl;
float scale = initScale*atomic;
cout << "Flow field at scale: " << scale << " pixels";
PerfTimer pt;
pt.Tic();
MatrixList flow = EdgeflowVectorField(im, 8, scale, scale, 1.0f, normalized);
double el = pt.Toc(false);
cout << " (" << el << " seconds)" << endl;
endScale = endScale*atomic;
scaleJump = scaleJump*atomic;
scale += scaleJump;
while( scale <= endScale)
{
Matrix efMag = Sqrt(flow(0)*flow(0)+flow(1)*flow(1));
float magLimit = Maximum(efMag(":"))/ratioLimit;
cout << "Flow field at scale: " << scale << " pixels";
pt.Tic();
MatrixList tempflow = EdgeflowVectorField(im, 4, scale, scale, 1.0f, normalized);
double el = pt.Toc(false);
cout << " (" << el << " seconds)" << endl;
Matrix tempMag = Sqrt(tempflow(0)*tempflow(0)+tempflow(1)*tempflow(1));
for (int i=0; i BB = Divergence(flow(0), flow(1));
float *tmp = poisson((-BB).Data(), BB.Columns(), BB.Rows());
Matrix CC(BB.Rows(), BB.Columns());
for (int j=0; j b(im.Rows()+6, im.Columns()+6,0);
b(3,im.Rows()+2,3,im.Columns()+2) = smoothWeighting*CC;
Matrix u(im.Rows()+6, im.Columns()+6,0);
u(3,im.Rows()+2,3,im.Columns()+2) = flow(0);
Matrix v(im.Rows()+6, im.Columns()+6,0);
v(3,im.Rows()+2,3,im.Columns()+2) = flow(1);
Vector< Matrix > phis(3);
phis[0] = Matrix(im.Rows()+6, im.Columns()+6);
phis[0](3,im.Rows()+2,3,im.Columns()+2) = im[0];
phis[1] = Matrix