Разработка системы автоматического поиска объектов на изображении



Приложение А

Руководство программиста

По итогам данной работы было реализовано два алгоритма. Алгоритм поиска максимума кросскорреляционной функции и метод SURF. Первый алгоритм реализован на C++, для хранения данных о изображениях была выбрана библиотека OpenCV, для быстрого вычисления дискретного преобразования Фурье была использована математическая библиотека FFTW. Метод реализован в виде одной процедуры На вход функции подаются изображение эталон, сцена и изображение для результата В качестве результата выдается изображение функция яркости которого соответствует функции кросскорреляции.

void phase_correlation( IplImage *ref, IplImage *tpl, IplImage *poc )

{

int i, j, k;

double tmp;

int width = ref?>width;

int height = ref?>height;

int step = ref?>widthStep;

int fft_size = width * height;

uchar *ref_data = ( uchar* ) ref?>imageData;

uchar *tpl_data = ( uchar* ) tpl?>imageData;

double *poc_data = ( double* )poc?>imageData;

fftw_complex *img1 = ( fftw_complex* )fftw_malloc( sizeof( fftw_complex ) * width * height );

fftw_complex *img2 = ( fftw_complex* )fftw_malloc( sizeof( fftw_complex ) * width * height );

fftw_complex *res = ( fftw_complex* )fftw_malloc( sizeof( fftw_complex ) * width * height );

fftw_plan fft_img1 = fftw_plan_dft_1d( width * height, img1, img1, FFTW_FORWARD, FFTW_ESTIMATE );

fftw_plan fft_img2 = fftw_plan_dft_1d( width * height, img2, img2, FFTW_FORWARD, FFTW_ESTIMATE );

fftw_plan ifft_res = fftw_plan_dft_1d( width * height, res, res, FFTW_BACKWARD, FFTW_ESTIMATE );

for( i = 0, k = 0 ; i < height ; i++ ) {

for( j = 0 ; j < width ; j++, k++ ) {

img1[k][0] = ( double )ref_data[i * step + j];

img1[k][1] = 0.0;

img2[k][0] = ( double )tpl_data[i * step + j];

img2[k][1] = 0.0;

}

}//подготовка к вычислению преобразования Фурье конец

fftw_execute( fft_img1 );//вычисление преобразования Фурье

fftw_execute( fft_img2 );

IplImage *fti1=cvCreateImage(cvSize(ref?>width,ref?>height),ref?>depth,ref?>nChannels);

//вычисление свертки

for( i = 0; i < fft_size ; i++ ){

res[i][0] = ( img2[i][0] * img1[i][0] ) ? ( img2[i][1] * ( ?img1[i][1] ) );

res[i][1] = ( img2[i][0] * ( ?img1[i][1] ) ) + ( img2[i][1] * img1[i][0] );

tmp = sqrt( pow( res[i][0], 2.0 ) + pow( res[i][1], 2.0 ) );

res[i][0] /= tmp;

res[i][1] /= tmp;

}

fftw_execute(ifft_res);// вычисление обратного преобразования Фурье

for( i = 0 ; i < fft_size ; i++ ) {

poc_data[i] = res[i][0] / ( double )fft_size;

}

fftw_destroy_plan( fft_img1 );

fftw_destroy_plan( fft_img2 );

fftw_destroy_plan( ifft_res );

fftw_free( img1 );

fftw_free( img2 );

fftw_free( res );

}

Метод SURF был реализован на С# с использованием объектно-ориентированного подхода. Входе работы над реализацией были реализованы следующие классы: IntegralImage, IPoint, SurfDescriptor, FastHessian.

Класс IPoint

Данный класс предназначен для хранения информации об особой точке, как положение, дескриптор, ориентация, длина дескриптора и другие значения.

public class IPoint{

public IPoint(){

orientation = 0;

}

public float x, y;

//позиция на изображении

public float scale;

//масштаб особой точки

public float response;

public float orientation;

//ориентация градиента

public int laplacian;

public int descriptorLength;

//длиннадескриптора

public float [] descriptor = null;

public void SetDescriptorLength(int Size){

descriptorLength = Size;

descriptor = new float[Size];

}

}

Класс IntegralImage

Данный класс предназначен для представления изображений в виде удобном для применения преобразований.

public class IntegralImage{

const float cR = .2989f;

const float cG = .5870f;

const float cB = .1140f;

internal float[,] Matrix;

public int Width, Height;

public float this[int y, int x]{

get { return Matrix[y, x]; }

set { Matrix[y, x] = value; }

}

private IntegralImage(int width, int height){

this.Width = width;

this.Height = height;27

this.Matrix = new float[height, width];

}

public static IntegralImage FromImage(Bitmap image){

IntegralImage pic = new IntegralImage(image.Width, image.Height);

float rowsum = 0;

for (int x = 0; x < image.Width; x++){

Color c = image.GetPixel(x, 0);

rowsum += (cR * c.R + cG * c.G + cB * c.B) / 255f;

pic[0, x] = rowsum;

}

for (int y = 1; y < image.Height; y++){

rowsum = 0;

for (int x = 0; x < image.Width; x++){

Color c = image.GetPixel(x, y);

rowsum += (cR * c.R + cG * c.G + cB * c.B) / 255f;

pic[y, x] = rowsum + pic[y ? 1, x];

}

}

return pic;

}

public float BoxIntegral(int row, int col, int rows, int cols){

int r1 = Math.Min(row, Height) ? 1;

int c1 = Math.Min(col, Width) ? 1;

int r2 = Math.Min(row + rows, Height) ? 1;

int c2 = Math.Min(col + cols, Width) ? 1;

float A = 0, B = 0, C = 0, D = 0;

if (r1 >= 0 && c1 >= 0) A = Matrix[r1, c1];

if (r1 >= 0 && c2 >= 0) B = Matrix[r1, c2];

if (r2 >= 0 && c1 >= 0) C = Matrix[r2, c1];

if (r2 >= 0 && c2 >= 0) D = Matrix[r2, c2];

return Math.Max(0, A ? B ? C + D);

}

public float HaarX(int row, int column, int size){

return BoxIntegral(row ? size / 2, column, size, size / 2)

? 1 * BoxIntegral(row ? size / 2, column ? size / 2, size, size / 2);

}

public float HaarY(int row, int column, int size){

return BoxIntegral(row, column ? size / 2, size / 2, size)

? 1 * BoxIntegral(row ? size / 2, column ? size / 2, size / 2, size);

}

}

Класс SurfDescriptor

Данный класс предназнчен непосредственно для вычисления дескрипторов особых точек.

public class SurfDescriptor{

float[,] gauss25 = new float[7, 7] {

{0.02350693969273f,0.01849121369071f,0.01239503121241f,0.00708015417522f,0.00344628101733f,0.0014294584748

4f,0.00050524879060f},

{0.02169964028389f,0.01706954162243f,0.01144205592615f,0.00653580605408f,0.00318131834134f,0.0013195564846

1f,0.00046640341759f},

{0.01706954162243f,0.01342737701584f,0.00900063997939f,0.00514124713667f,0.00250251364222f,0.0010379998950

4f,0.00036688592278f},

{0.01144205592615f,0.00900063997939f,0.00603330940534f,0.

3f,0.00024593098864f},

{0.00653580605408f,0.00514124713667f,0.00344628101733f,0.00196854695367f,0.00095819467066f,0.0003974427754

6f,0.00014047800980f},

{0.00318131834134f,0.00250251364222f,0.00167748505986f,0.00095819467066f,0.00046640341759f,0.0001934561675

7f,0.00006837798818f},

{0.00131955648461f,0.00103799989504f,0.00069579213743f,0.00039744277546f,0.00019345616757f,0.0000802423124

7f,0.00002836202103f}

};

IntegralImage img;

public static void DecribeInterestPoints(List<IPoint> ipts, bool upright, bool extended, IntegralImageimg){

SurfDescriptor des = new SurfDescriptor();

des.DescribeInterestPoints(ipts, upright, extended, img);

}

public void DescribeInterestPoints(List<IPoint> ipts, bool upright, bool extended, IntegralImage img){

if (ipts.Count == 0) return;

this.img = img;

foreach (IPoint ip in ipts){

if (extended) ip.descriptorLength = 128;

else ip.descriptorLength = 64;

if (!upright) GetOrientation(ip);

GetDescriptor(ip, upright, extended);

}

}

void GetOrientation(IPoint ip){

const byte Responses = 109;

float[] resX = new float[Responses];

float[] resY = new float[Responses];

float[] Ang = new float[Responses];

int idx = 0;

int[] id = { 6, 5, 4, 3, 2, 1, 0, 1, 2, 3, 4, 5, 6 };

int X = (int)Math.Round(ip.x, 0);

int Y = (int)Math.Round(ip.y, 0);

int S = (int)Math.Round(ip.scale, 0);

for (int i = ?6; i <= 6; ++i){

for (int j = ?6; j <= 6; ++j){

if (i * i + j * j < 36){

float gauss = gauss25[id[i + 6], id[j + 6]];

resX[idx] = gauss * img.HaarX(Y + j * S, X + i * S, 4 * S);

resY[idx] = gauss * img.HaarY(Y + j * S, X + i * S, 4 * S);

Ang[idx] = (float)GetAngle(resX[idx], resY[idx]);

++idx;

}

}

}

float sumX, sumY, max = 0, orientation = 0;

float ang1, ang2;

float pi = (float)Math.PI;

for (ang1 = 0; ang1 < 2 * pi; ang1 += 0.15f){

ang2 = (ang1 + pi / 3f > 2 * pi ? ang1 ? 5 * pi / 3f : ang1 + pi / 3f);

sumX = sumY = 0;

for (int k = 0; k < Responses; ++k){

if (ang1 < ang2 && ang1 < Ang[k] && Ang[k] < ang2){

sumX += resX[k];

sumY += resY[k];

}

else if (ang2 < ang1 &&((Ang[k] > 0 && Ang[k] < ang2) || (Ang[k] > ang1 && Ang[k] < pi))){

sumX += resX[k];

sumY += resY[k];

}

}

if (sumX * sumX + sumY * sumY > max){

max = sumX * sumX + sumY * sumY;

orientation = (float)GetAngle(sumX, sumY);

}

}

ip.orientation = (float)orientation;

}

void GetDescriptor(IPoint ip, bool bUpright, bool bExtended){

int sample_x, sample_y, count = 0;

int i = 0, ix = 0, j = 0, jx = 0, xs = 0, ys = 0;

float dx, dy, mdx, mdy, co, si;

float dx_yn, mdx_yn, dy_xn, mdy_xn;

float gauss_s1 = 0f, gauss_s2 = 0f;

float rx = 0f, ry = 0f, rrx = 0f, rry = 0f, len = 0f;

float cx = ?0.5f, cy = 0f;

int X = (int)Math.Round(ip.x, 0);

int Y = (int)Math.Round(ip.y, 0);

int S = (int)Math.Round(ip.scale, 0);

ip.SetDescriptorLength(64);

if (bUpright){

co = 1;

si = 0;

}

else{

co = (float)Math.Cos(ip.orientation);

si = (float)Math.Sin(ip.orientation);

}

i = ?8;

while (i < 12){

j = ?8;

i = i ? 4;

cx += 1f;

cy = ?0.5f;

while (j < 12){

cy += 1f;

j = j ? 4;

ix = i + 5;

jx = j + 5;

dx = dy = mdx = mdy = 0f;

dx_yn = mdx_yn = dy_xn = mdy_xn = 0f;

xs = (int)Math.Round(X + (?jx * S * si + ix * S * co), 0);

ys = (int)Math.Round(Y + (jx * S * co + ix * S * si), 0);

dx = dy = mdx = mdy = 0f;

dx_yn = mdx_yn = dy_xn = mdy_xn = 0f;

for (int k = i; k < i + 9; ++k){

for (int l = j; l < j + 9; ++l){

sample_x = (int)Math.Round(X + (?l * S * si + k * S * co), 0);

sample_y = (int)Math.Round(Y + (l * S * co + k * S * si), 0);

gauss_s1 = Gaussian(xs ? sample_x, ys ? sample_y, 2.5f * S);

rx = (float)img.HaarX(sample_y, sample_x, 2 * S);

ry = (float)img.HaarY(sample_y, sample_x, 2 * S);

rrx = gauss_s1 * (?rx * si + ry * co);

rry = gauss_s1 * (rx * co + ry * si);

if (bExtended){

if (rry >= 0){

dx += rrx;

mdx += Math.Abs(rrx);

}

else{

dx_yn += rrx;

mdx_yn += Math.Abs(rrx);

}

if (rrx >= 0){

dy += rry;

mdy += Math.Abs(rry);

}

else{

dy_xn += rry;

mdy_xn += Math.Abs(rry);

}

}

else{

dx += rrx;

dy += rry;

mdx += Math.Abs(rrx);

mdy += Math.Abs(rry);

}

}

}

gauss_s2 = Gaussian(cx ? 2f, cy ? 2f, 1.5f);

ip.descriptor[count++] = dx * gauss_s2;

ip.descriptor[count++] = dy * gauss_s2;

ip.descriptor[count++] = mdx * gauss_s2;

ip.descriptor[count++] = mdy * gauss_s2;

if (bExtended){

ip.descriptor[count++] = dx_yn * gauss_s2;

ip.descriptor[count++] = dy_xn * gauss_s2;

ip.descriptor[count++] = mdx_yn * gauss_s2;

ip.descriptor[count++] = mdy_xn * gauss_s2;

}

len += (dx * dx + dy * dy + mdx * mdx + mdy * mdy

+ dx_yn + dy_xn + mdx_yn + mdy_xn) * gauss_s2 * gauss_s2;

j += 9;

}

i += 9;

}

len = (float)Math.Sqrt((double)len);

if (len > 0){

for (int d = 0; d < ip.descriptorLength; ++d){

ip.descriptor[d] /= len;

double GetAngle(float X, float Y){

if (X >= 0 && Y >= 0)

return Math.Atan(Y / X);

else if (X < 0 && Y >= 0)

return Math.PI ? Math.Atan(?Y / X);

else if (X < 0 && Y < 0)

return Math.PI + Math.Atan(Y / X);

else if (X >= 0 && Y < 0)

return 2 * Math.PI ? Math.Atan(?Y / X);

return 0;

float Gaussian(int x, int y, float sig){

float pi = (float)Math.PI;

return (1f / (2f * pi * sig * sig)) * (float)Math.Exp(?(x * x + y * y) / (2.0f * sig * sig));

float Gaussian(float x, float y, float sig){

float pi = (float)Math.PI;

return 1f / (2f * pi * sig * sig) * (float)Math.Exp(?(x * x + y * y) / (2.0f * sig * sig));

}

Размещено на Allbest.ru