Перейти к содержимому
Compvision.ru
mrgloom

Сложение преобразований через remap

Автор: , в OpenCV

Recommended Posts

mrgloom    242

mrgloom
Опубликовано

вообщем есть 2 функции одна двигает углы изображения, другая делает преобразование подушка\бочка, хотелось бы объединить их сразу в 1 через remap чтобы выполнялось быстрей.вообще даже скорей вопрос как фун-ию warp_patch запихнуть в рамки mapx,mapy.

похожий вопрос тут обсуждается http://tech.groups.yahoo.com/group/OpenCV/message/40225

IplImage* warp_patch(IplImage* src, std::vector<coor> vec, bool crop) 

{

	CvPoint2D32f srcQuad[4], dstQuad[4];

	CvMat*       warp_matrix = cvCreateMat(3,3,CV_32FC1);


	IplImage* dst= NULL;

	IplImage* src_ext= NULL;


	//src координаты

	srcQuad[0].x = 0;             //src Top left

	srcQuad[0].y = 0;

	srcQuad[1].x = src->width-1;  //src Top right

	srcQuad[1].y = 0;

	srcQuad[2].x = 0;             //src Bottom left

	srcQuad[2].y = src->height-1;

	srcQuad[3].x = src->width-1;  //src Bot right

	srcQuad[3].y = src->height-1;


	//отрицательные значения cx налево

	//						 cy вверх


	if (vec.size()!=4) //проверка вектора на содержание 4 углов

	return NULL;


	int dx= abs(MIN(vec[0].x,vec[2].x))+MAX(vec[1].x,vec[3].x);

	int dy= abs(MIN(vec[0].y,vec[1].y))+MAX(vec[2].y,vec[3].y);

	dst= cvCreateImage(cvSize(src->width+dx,src->height+dy),IPL_DEPTH_8U,3);

	src_ext= cvCreateImage(cvSize(src->width+dx,src->height+dy),IPL_DEPTH_8U,3);//из-за требования на одинаковый размер src и dst в cvWarpPerspective

	cvSetImageROI(src_ext,cvRect(0,0,src->width,src->height));

	cvResetImageROI(src);

	cvCopy(src,src_ext);


	//1 сторона 0 2 

	if (vec[0].x<vec[2].x)

	{

		if(vec[0].x<0)

		{

			dstQuad[0].x= 0;

			dstQuad[2].x= vec[2].x+abs(vec[0].x);

		}

		else

		{

			dstQuad[0].x= 0;

			dstQuad[2].x= vec[2].x-vec[0].x;

		}

	}

	else

	{

		if(vec[0].x<0)

		{

			dstQuad[0].x= vec[0].x-(vec[0].x);

			dstQuad[2].x= 0;

		}

		else

		{

			dstQuad[0].x= vec[0].x-vec[2].x;

			dstQuad[2].x= 0;

		}

	}

	//2 сторона 0 1 

	if (vec[0].y<vec[1].y)

	{

		if(vec[0].y<0)

		{

			dstQuad[0].y= 0;

			dstQuad[1].y= vec[1].y+abs(vec[0].y);

		}

		else

		{

			dstQuad[0].y= 0;

			dstQuad[1].y= vec[1].y-vec[0].y;

		}

	}

	else

	{

		if(vec[0].y<0)

		{

			dstQuad[0].y= vec[0].y-(vec[0].y);

			dstQuad[1].y= 0;

		}

		else

		{

			dstQuad[0].y= vec[0].y-vec[1].y;

			dstQuad[1].y= 0;

		}

	}

	//3 сторона 1 3 

	dstQuad[1].x= abs(MIN(vec[0].x,vec[2].x))+src->width-1+vec[1].x;

	dstQuad[3].x= abs(MIN(vec[0].x,vec[2].x))+src->width-1+vec[3].x;


	//4 сторона 2 3 

	dstQuad[2].y= abs(MIN(vec[0].y,vec[1].y))+src->height-1+vec[2].y;

	dstQuad[3].y= abs(MIN(vec[0].y,vec[1].y))+src->height-1+vec[3].y;


	cvGetPerspectiveTransform(srcQuad,dstQuad,warp_matrix);

	cvWarpPerspective( src_ext, dst, warp_matrix );


	if (crop) // обрезается максимальный вписаный прямоугольник

	{

		int cx1= MAX(dstQuad[0].x,dstQuad[2].x);

		int cx2= MIN(dstQuad[1].x,dstQuad[3].x);

		int cy1= MAX(dstQuad[0].y,dstQuad[1].y);

		int cy2= MIN(dstQuad[2].y,dstQuad[3].y);


		IplImage* dst_crop= cvCreateImage(cvSize(cx2-cx1,cy2-cy1),IPL_DEPTH_8U,3);

		cvSetImageROI(dst,cvRect(cx1,cy1,cx2-cx1,cy2-cy1));

		cvCopy(dst,dst_crop);

		return dst_crop;

	}

	else

	{

		return dst;

	}

}

void barrel_pincusion_dist(IplImage* img, double kx, double ky)

{


	IplImage* mapx = cvCreateImage( cvGetSize(img), IPL_DEPTH_32F, 1 );

	IplImage* mapy = cvCreateImage( cvGetSize(img), IPL_DEPTH_32F, 1 );


	int w= img->width;

	int h= img->height;

	int Cx= w/2;

	int Cy= h/2; 


	for (int y = 0; y < h; y++)

	{

		int ty= y-Cy;

		for (int x = 0; x < w; x++)

		{

			int tx= x-Cx;

			int rt= tx*tx+ty*ty;

			double u= tx*(1+kx*rt)+Cx;

			cvSetReal2D(mapx, y, x, u);

		}

	}


	for (int y = 0;y < h; y++)

	{

		int ty= y-Cy;

		for (int x = 0; x < w; x++) 

		{


			int tx= x-Cx;

			int rt= tx*tx+ty*ty;

			double v= ty*(1+ky*rt)+Cy;

			cvSetReal2D(mapy, y, x, v);

		}

	}


	IplImage* temp = cvCloneImage(img);

    cvRemap( temp, img, mapx, mapy ); 

	cvReleaseImage(&temp);

	cvReleaseImage(&mapx);

	cvReleaseImage(&mapy);

}

и еще может есть что побыстрей cvSetReal2D?

Поделиться сообщением

Ссылка на сообщение
Поделиться на других сайтах

Nuzhny    243

Nuzhny
Опубликовано

Отвечу на последний вопрос темы - над первым надо ещё подумать. Быстрее так:

void barrel_pincusion_dist(IplImage* img, double kx, double ky)

{

	IplImage* mapx = cvCreateImage( cvGetSize(img), IPL_DEPTH_32F, 1 );

	IplImage* mapy = cvCreateImage( cvGetSize(img), IPL_DEPTH_32F, 1 );


	int w= img->width;

	int h= img->height;

	int Cx= w/2;

	int Cy= h/2; 


	float* pbuf = (float*)mapx->imageData;

	for (int y = 0; y < h; y++)

	{

		int ty= y-Cy;

		for (int x = 0; x < w; x++)

		{

			int tx= x-Cx;

			int rt= tx*tx+ty*ty;


			*pbuf = (float)(tx*(1+kx*rt)+Cx);

			++pbuf;

		}

	}


	pbuf = (float*)mapx->imageData;

	for (int y = 0;y < h; y++)

	{

		int ty= y-Cy;

		for (int x = 0; x < w; x++) 

		{


			int tx= x-Cx;

			int rt= tx*tx+ty*ty;


			*pbuf = (float)(ty*(1+ky*rt)+Cy);

			++pbuf;

		}

	}


	IplImage* temp = cvCloneImage(img);

	cvRemap( temp, img, mapx, mapy ); 

	cvReleaseImage(&temp);

	cvReleaseImage(&mapx);

	cvReleaseImage(&mapy);

}

Поделиться сообщением

Ссылка на сообщение
Поделиться на других сайтах

mrgloom    242

mrgloom
Опубликовано

только там вычисляется как double, а тут использован указатель float. (c указателями double почему то не работает)

мне кажется надо как то warp_patch переписать под mapx,mapy и сложить с mapx,mapy которые получаются от barrel_pincusion_dist.

(только не уверен, что просто сложить будет правильно)

Поделиться сообщением

Ссылка на сообщение
Поделиться на других сайтах

mrgloom    242

mrgloom
Опубликовано

еще вопрос как распараллелить warp_patch (имеется ввиду либо подача изображения частями в разные треды, либо как то еще)

т.е. даже этот вот кусок


        cvWarpPerspective( src_ext, dst, warp_matrix );

Поделиться сообщением

Ссылка на сообщение
Поделиться на других сайтах

mrgloom    242

mrgloom
Опубликовано

вот нашел, получается все завязано на ремап и надо произвести постобработку на конечном изображении, которая возникает из-за деления на куски.

void warpPerspective( const Mat& src, Mat& dst, const Mat& M0, Size dsize,

                      int flags, int borderType, const Scalar& borderValue )

{

    dst.create( dsize, src.type() );

    CV_Assert( dst.data != src.data );


    const int BLOCK_SZ = 32;

    short XY[BLOCK_SZ*BLOCK_SZ*2], A[BLOCK_SZ*BLOCK_SZ];

    double M[9];

    Mat matM(3, 3, CV_64F, M);

    int interpolation = flags & INTER_MAX;

    if( interpolation == INTER_AREA )

        interpolation = INTER_LINEAR;


    CV_Assert( (M0.type() == CV_32F || M0.type() == CV_64F) && M0.rows == 3 && M0.cols == 3 );

    M0.convertTo(matM, matM.type());


    if( !(flags & WARP_INVERSE_MAP) )

         invert(matM, matM);


    int x, y, x1, y1, width = dst.cols, height = dst.rows;


    int bh0 = std::min(BLOCK_SZ/2, height);

    int bw0 = std::min(BLOCK_SZ*BLOCK_SZ/bh0, width);

    bh0 = std::min(BLOCK_SZ*BLOCK_SZ/bw0, height);


    for( y = 0; y < height; y += bh0 )

    {

        for( x = 0; x < width; x += bw0 )

        {

            int bw = std::min( bw0, width - x);

            int bh = std::min( bh0, height - y);


            Mat _XY(bh, bw, CV_16SC2, XY), matA;

            Mat dpart(dst, Rect(x, y, bw, bh));


            for( y1 = 0; y1 < bh; y1++ )

            {

                short* xy = XY + y1*bw*2;

                double X0 = M[0]*x + M[1]*(y + y1) + M[2];

                double Y0 = M[3]*x + M[4]*(y + y1) + M[5];

                double W0 = M[6]*x + M[7]*(y + y1) + M[8];


                if( interpolation == INTER_NEAREST )

                    for( x1 = 0; x1 < bw; x1++ )

                    {

                        double W = W0 + M[6]*x1;

                        W = W ? 1./W : 0;

                        int X = saturate_cast<int>((X0 + M[0]*x1)*W);

                        int Y = saturate_cast<int>((Y0 + M[3]*x1)*W);

                        xy[x1*2] = (short)X;

                        xy[x1*2+1] = (short)Y;

                    }

                else

                {

                    short* alpha = A + y1*bw;

                    for( x1 = 0; x1 < bw; x1++ )

                    {

                        double W = W0 + M[6]*x1;

                        W = W ? INTER_TAB_SIZE/W : 0;

                        int X = saturate_cast<int>((X0 + M[0]*x1)*W);

                        int Y = saturate_cast<int>((Y0 + M[3]*x1)*W);

                        xy[x1*2] = (short)(X >> INTER_BITS);

                        xy[x1*2+1] = (short)(Y >> INTER_BITS);

                        alpha[x1] = (short)((Y & (INTER_TAB_SIZE-1))*INTER_TAB_SIZE +

                                (X & (INTER_TAB_SIZE-1)));

                    }

                }

            }


            if( interpolation == INTER_NEAREST )

                remap( src, dpart, _XY, Mat(), interpolation, borderType, borderValue );

            else

            {

                Mat matA(bh, bw, CV_16U, A);

                remap( src, dpart, _XY, matA, interpolation, borderType, borderValue );

            }

        }

    }

}

Поделиться сообщением

Ссылка на сообщение
Поделиться на других сайтах

Nuzhny    243

Nuzhny
Опубликовано

remap задействована в самом внутреннем цикле. Значит, можно разделить изображение по высоте на число потоков и отдавать каждому из них свою "полоску" исходного изображения.

Можно предложить такое решение:

void block_warpPerspective(int y_from, int y_to)

{

	for( y = y_from; y < y_to; y += bh0 )

	{

		for( x = 0; x < width; x += bw0 )

		{

			int bw = std::min( bw0, width - x);

			int bh = std::min( bh0, height - y);


			Mat _XY(bh, bw, CV_16SC2, XY), matA;

			Mat dpart(dst, Rect(x, y, bw, bh));


			for( y1 = 0; y1 < bh; y1++ )

			{

				short* xy = XY + y1*bw*2;

				double X0 = M[0]*x + M[1]*(y + y1) + M[2];

				double Y0 = M[3]*x + M[4]*(y + y1) + M[5];

				double W0 = M[6]*x + M[7]*(y + y1) + M[8];


				if( interpolation == INTER_NEAREST )

					for( x1 = 0; x1 < bw; x1++ )

					{

						double W = W0 + M[6]*x1;

						W = W ? 1./W : 0;

						int X = saturate_cast<int>((X0 + M[0]*x1)*W);

						int Y = saturate_cast<int>((Y0 + M[3]*x1)*W);

						xy[x1*2] = (short)X;

						xy[x1*2+1] = (short)Y;

					}

				else

				{

					short* alpha = A + y1*bw;

					for( x1 = 0; x1 < bw; x1++ )

					{

						double W = W0 + M[6]*x1;

						W = W ? INTER_TAB_SIZE/W : 0;

						int X = saturate_cast<int>((X0 + M[0]*x1)*W);

						int Y = saturate_cast<int>((Y0 + M[3]*x1)*W);

						xy[x1*2] = (short)(X >> INTER_BITS);

						xy[x1*2+1] = (short)(Y >> INTER_BITS);

						alpha[x1] = (short)((Y & (INTER_TAB_SIZE-1))*INTER_TAB_SIZE +

							(X & (INTER_TAB_SIZE-1)));

					}

				}

			}


			if( interpolation == INTER_NEAREST )

				remap( src, dpart, _XY, Mat(), interpolation, borderType, borderValue );

			else

			{

				Mat matA(bh, bw, CV_16U, A);

				remap( src, dpart, _XY, matA, interpolation, borderType, borderValue );

			}

		}

	}

}



void parallel_warpPerspective( const Mat& src, Mat& dst, const Mat& M0, Size dsize,

	int flags, int borderType, const Scalar& borderValue )

{

	dst.create( dsize, src.type() );

	CV_Assert( dst.data != src.data );


	const int BLOCK_SZ = 32;

	short XY[BLOCK_SZ*BLOCK_SZ*2], A[BLOCK_SZ*BLOCK_SZ];

	double M[9];

	Mat matM(3, 3, CV_64F, M);

	int interpolation = flags & INTER_MAX;

	if( interpolation == INTER_AREA )

		interpolation = INTER_LINEAR;


	CV_Assert( (M0.type() == CV_32F || M0.type() == CV_64F) && M0.rows == 3 && M0.cols == 3 );

	M0.convertTo(matM, matM.type());


	if( !(flags & WARP_INVERSE_MAP) )

		invert(matM, matM);


	int x, y, x1, y1, width = dst.cols, height = dst.rows;


	int bh0 = std::min(BLOCK_SZ/2, height);

	int bw0 = std::min(BLOCK_SZ*BLOCK_SZ/bh0, width);

	bh0 = std::min(BLOCK_SZ*BLOCK_SZ/bw0, height);


	// Число потоков

	int num_threads = 4;

	for (i = 0; i < num_threads; ++i)

	{

		// Типа старт потока с функцией block_warpPerspective и двумя её аргументами y_from, y_to

		start_thread(block_warpPerspective, (i * height) / num_threads, ((i + 1) * height) / num_threads);

	}

}

Вместо warpPerspective вызывать parallel_warpPerspective. Функции start_thread и block_warpPerspective - условные, необходимо заменить на конкретную реализацию. Надеюсь, что моя мысль ясна, но могу пояснить или набросать компилирующийся параллельный код с бустовскими потоками (только мне нужен какой-нибудь рабочий последовательный).

Поделиться сообщением

Ссылка на сообщение
Поделиться на других сайтах

mrgloom    242

mrgloom
Опубликовано 

remap задействована в самом внутреннем цикле. Значит, можно разделить изображение по высоте на число потоков и отдавать каждому из них свою "полоску" исходного изображения.

про полоску то понятно, но дело в том как эти полоски потом сложатся в картинку,т.к. между полосками нету интерполяции никакой.

Поделиться сообщением

Ссылка на сообщение
Поделиться на других сайтах

Nuzhny    243

Nuzhny
Опубликовано

Надо пробовать

Поделиться сообщением

Ссылка на сообщение
Поделиться на других сайтах

mrgloom    242

mrgloom
Опубликовано

вот написал код для распараллеливания бочки, remap пока не трогал, только вычисления.

попробуйте на тредах может будет какой то прирост.

void parallel_barrel(float* pbufx, float* pbufy,int Y, int w, int H, double kx, double ky, int Cx, int Cy)

{

	//через указатели (может хватит точности float?) double через указатели не работает

	for (int y = Y; y < Y+H; y++)

	{

		int ty= y-Cy;

		for (int x = 0; x < w; x++)

		{

			int tx= x-Cx;

			int rt= tx*tx+ty*ty;


			*pbufx = (float)(tx*(1+kx*rt)+Cx);

			++pbufx;

		}

	}


	for (int y = Y;y < Y+H; y++)

	{

		int ty= y-Cy;

		for (int x = 0; x < w; x++) 

		{

			int tx= x-Cx;

			int rt= tx*tx+ty*ty;


			*pbufy = (float)(ty*(1+ky*rt)+Cy);

			++pbufy;

		}

	}

}

//распаралеливаем по mapx mapy 

	IplImage* img= cvLoadImage("1.PNG");

	IplImage* mapx= cvCreateImage( cvGetSize(img), IPL_DEPTH_32F, 1 );

        IplImage* mapy= cvCreateImage( cvGetSize(img), IPL_DEPTH_32F, 1 );

	int w= img->width;

        int h= img->height;

        int Cx= w/2;

        int Cy= h/2; 

	const double kx= 3.8E-7;

	const double ky= 3.8E-7;

	float* pbufx = (float*)mapx->imageData;

	float* pbufy = (float*)mapy->imageData;


	//parallel section

	//вместо цикла поставить разделение на треды

	int N=4; //кол-во тредов

	int H= floor((float)h/N);

	int y=0;

	for (int i=0;i<N;++i)

	{

		parallel_barrel(pbufx, pbufy,y, w, H, kx, ky, Cx, Cy);

		if (i!=N-1)

		{

			pbufx+=H*w;

			pbufy+=H*w;

		}

		y+=H;

	}

        //end parallel section


	IplImage* temp = cvCloneImage(img);

        cvRemap( temp, img, mapx, mapy ); 

        cvReleaseImage(&temp);

        cvReleaseImage(&mapx);

        cvReleaseImage(&mapy);

	cvSaveImage("out.png",img);

Поделиться сообщением

Ссылка на сообщение
Поделиться на других сайтах

Nuzhny    243

Nuzhny
Опубликовано

Вот код с бустом:

#include <cstdlib>

#include <iostream>

#include <deque>


#include <boost/bind.hpp>

#include <boost/smart_ptr.hpp>

#include <boost/thread.hpp>


#include <cv.h>

#include <highgui.h>


////////////////////////////////////////////////////////////////////////////


void parallel_barrel(float* pbufx, float* pbufy, int Y, int w, int H, float kx, float ky, int Cx, int Cy)

{

	//через указатели (может хватит точности float?) float через указатели не работает

	for (int y = Y; y < Y+H; y++)

	{

		int ty= y-Cy;

		for (int x = 0; x < w; x++)

		{

			int tx= x-Cx;

			int rt= tx*tx+ty*ty;


			*pbufx = (float)(tx*(1+kx*rt)+Cx);

			++pbufx;

		}

	}


	for (int y = Y;y < Y+H; y++)

	{

		int ty= y-Cy;

		for (int x = 0; x < w; x++) 

		{

			int tx= x-Cx;

			int rt= tx*tx+ty*ty;


			*pbufy = (float)(ty*(1+ky*rt)+Cy);

			++pbufy;

		}

	}

}

////////////////////////////////////////////////////////////////////////////



int _tmain(int argc, _TCHAR* argv[])

{

	argc; argv;


	//распаралеливаем по mapx mapy 

	IplImage* img= cvLoadImage("1.PNG");

	IplImage* mapx= cvCreateImage( cvGetSize(img), IPL_DEPTH_32F, 1 );

	IplImage* mapy= cvCreateImage( cvGetSize(img), IPL_DEPTH_32F, 1 );

	int w= img->width;

	int h= img->height;

	int Cx= w/2;

	int Cy= h/2; 

	const float kx= 3.8E-7;

	const float ky= 3.8E-7;

	float* pbufx = (float*)mapx->imageData;

	float* pbufy = (float*)mapy->imageData;


	//parallel section

	//вместо цикла поставить разделение на треды

	int N = 4; //кол-во тредов

	int H= h / N;

	int y = 0;

	 boost::thread_group threads;

	for (int i = 0; i < N; ++i)

	{

		threads.add_thread(new boost::thread(boost::bind(parallel_barrel, pbufx, pbufy, y, w, H, kx, ky, Cx, Cy)));

		//parallel_barrel(pbufx, pbufy, y, w, H, kx, ky, Cx, Cy);


		if (i != N - 1)

		{

			pbufx += H * w;

			pbufy += H * w;

		}

		y += H;

	}

	//end parallel section


	threads.join_all();


	IplImage* temp = cvCloneImage(img);

	cvRemap(temp, img, mapx, mapy);

	cvReleaseImage(&temp);

	cvReleaseImage(&mapx);

	cvReleaseImage(&mapy);

	cvSaveImage("out.png",img);


	return 0;

}

Как видишь, пользоваться бустовскими потоками очень просто, код практически не изменился.

У меня дома комп - старенький Core 2 Duo, поэтому ускорения не заметил никакого.

Поделиться сообщением

Ссылка на сообщение
Поделиться на других сайтах

TroyashkA    12

TroyashkA
Опубликовано

Nuzhny, boost::thread_group напоминает ThreadPool из .NET

mrgloom, напомните, что делает remap()? Смещает координаты изображения? А назначение функции warp_patсh()?

Поделиться сообщением

Ссылка на сообщение
Поделиться на других сайтах

mrgloom    242

mrgloom
Опубликовано 

cvRemap( src, dst, mapx, mapy );
remap преобразует изображения из старых координат в новые+делает интерполяцию. сначала хочу распараллелить warpPerspective, а потом может дойду и до remap. вообще если брать тупо циклы то можно их распараллелить с помощью OpenMP тупо parallel for. (вообщем попробую) с бустом пока не работал, да и с тредами тоже, вроде есть какие то более простые вин треды? (чтобы буст не подлкючать) 
warp_patсh()

двигает углы изображения, в std::vector<coor> vec хранится смещение для 4 углов по х,y.

для Core 2 Duo наверно надо было 2 треда.

Поделиться сообщением

Ссылка на сообщение
Поделиться на других сайтах

Nuzhny    243

Nuzhny
Опубликовано

Можно и Вин-треды сделать. Там тоже несложно.

Сделай parallel_barrel как unsigned __stdcall parallel_barrel(void *pArguments) (pArguments - это и есть указатель на структуру с входными параметрами). Входные данные для неё упакуй в одну структуру и передавай в функцию _beginthreadex. Все хендлы потоков складывай в массив, который подавай в функцию ожидания завершения всех потоков WaitForMultipleObjects. Чуть сложнее и более громоздко чем буст.

Поделиться сообщением

Ссылка на сообщение
Поделиться на других сайтах

mrgloom    242

mrgloom
Опубликовано

с openmp оказалось как то не просто, непонятно как передавать правильный указатель pbuf в тред.

сейчас выдает кашу.

IplImage* barrel_dist(IplImage* img, double kx, double ky)

{

        IplImage* mapx = cvCreateImage( cvGetSize(img), IPL_DEPTH_32F, 1 );

        IplImage* mapy = cvCreateImage( cvGetSize(img), IPL_DEPTH_32F, 1 );


        int w= img->width;

        int h= img->height;

        int Cx= w/2;

        int Cy= h/2; 


		const int n_threads=2;

		//omp_set_dynamic(0);      // запретить библиотеке openmp менять число потоков во время исполнения

		omp_set_num_threads(n_threads); // установить число потоков


        float* pbuf = (float*)mapx->imageData;

		int y=0;

		int x=0;

		int tx=0;

		int ty=0;

		int rt=0;

		#pragma omp parallel private(x,y,tx,ty,rt)

		{

			#pragma omp for schedule (static,h/n_threads) 

			for (int y = 0; y < h; y++)

			{

				//int t= omp_in_parallel();//

				//int f= omp_get_num_threads();//

				int ty= y-Cy;

				for (int x = 0; x < w; x++)

				{

					//int t= omp_in_parallel();//

					//int f= omp_get_num_threads();//

					int tx= x-Cx;

					int rt= tx*tx+ty*ty;


					*pbuf = (float)(tx*(1+kx*rt)+Cx);

					++pbuf;

				}

			}

		}


        pbuf = (float*)mapy->imageData;

		y=0;

		x=0;

		tx=0;

		ty=0;

		rt=0;

		#pragma omp parallel private(x,y,tx,ty,rt)

		{

			#pragma omp for schedule (static,h/n_threads) 

			for (int y = 0;y < h; y++)

			{

				//int f= omp_get_num_threads();//

				int ty= y-Cy;

				for (int x = 0; x < w; x++) 

				{

					//int f= omp_get_num_threads();//

					int tx= x-Cx;

					int rt= tx*tx+ty*ty;


					*pbuf = (float)(ty*(1+ky*rt)+Cy);

					++pbuf;

				}

			}

		}


        IplImage* temp = cvCloneImage(img);

        cvRemap( temp, img, mapx, mapy ); 

        cvReleaseImage(&temp);

        cvReleaseImage(&mapx);

        cvReleaseImage(&mapy);

		return img;

}

Поделиться сообщением

Ссылка на сообщение
Поделиться на других сайтах

mrgloom    242

mrgloom
Опубликовано

опять возвращаясь к теме

во-первых вопрос, возможно ли как то определить по конечной картинке порядок наложения дисторсий?

если знать начальную картинку и набор возможных дисторсий.

во-вторых

если поворот описывается с помощью

x'= x*cosa-y*sina;

y'= x*siny+y*cosa;

а подушка\бочка

x'= x+Const*(x*x+y*y)*x;

y'= y+Const*(x*x+y*y)*y;

то как описывается преобразование трапеции?

я так понимаю это частный случай Affine homography.

post-701-0-34467300-1332141206_thumb.png

но как его из матричного представления переписать?

Поделиться сообщением

Ссылка на сообщение
Поделиться на других сайтах

mrgloom    242

mrgloom
Опубликовано

хотя если если исходить из этого

image185.gif

то

x'= m1*x+m3*y+t1;

y'= m2*x+m4*y+t2;

Поделиться сообщением

Ссылка на сообщение
Поделиться на других сайтах

mrgloom    242

mrgloom
Опубликовано

вообщем я вспомнил это уже обсуждалось, и кажется пора вспомнить линал.

хотя все равно непонятно как матрица 2х1 [x,y] умножается на матрицу гомографии 3х3.

в итоге общее для матрицы гомографии

u= (m11*x + m12*y + m13)/(m31*x + m32*y + m33);

v= (m21*x + m22*y + m23)/(m31*x + m32*y + m33);

как раз для аффинных m31=m32=0 m33=1

u= (m11*x + m12*y + m13); //так что сверху я перепутал коэффициенты

v= (m21*x + m22*y + m23);

и обратное (возможно не оптимальное)

x = (-V * m33 * m12 + V * m32 * m13 + m23 * m12 - m22 * m13 + m22 * U * m33 - m23 * U * m32) / (m11 * m22 - m11 * V * m32 - m12 * m21 + m12 * V * m31 - U * m31 * m22 + U * m32 * m21);

t2 = m11 * V;

t7 = U * m31;

y = -(m11 * m23 - t2 * m33 - m13 * m21 + m13 * V * m31 - t7 * m23 + U * m33 * m21) / (m11 * m22 - t2 * m32 - m12 * m21 + m12 * V * m31 - t7 * m22 + U * m32 * m21);

Поделиться сообщением

Ссылка на сообщение
Поделиться на других сайтах

Smorodov    579

Smorodov
Опубликовано

там не [x;y], а [x;y;1] -- гомогенные координаты.

вот

http://en.wikipedia.org/wiki/Homogeneous_coordinates

и вот:

http://www.packtpub.com/article/opencv-estimating-projective-relations-images

Поделиться сообщением

Ссылка на сообщение
Поделиться на других сайтах

mrgloom    242

mrgloom
Опубликовано

какой то математический трюк

In affine transformations an extra dimension is added and all points are given a value of 1 for this extra dimension. The advantage of doing this is that then all of the euclidean transformations become linear transformations and can be represented using matrix multiplication.

projective geometry это как по русски и относительно линала называется?

Поделиться сообщением

Ссылка на сообщение
Поделиться на других сайтах

mrgloom    242

mrgloom
Опубликовано

вроде нашел, там даже книжки есть какие то.

http://ru.wikipedia.org/wiki/%D0%9F%D1%80%D0%BE%D0%B5%D0%BA%D1%82%D0%B8%D0%B2%D0%BD%D0%B0%D1%8F_%D0%B3%D0%B5%D0%BE%D0%BC%D0%B5%D1%82%D1%80%D0%B8%D1%8F

Поделиться сообщением

Ссылка на сообщение
Поделиться на других сайтах

mrgloom    242

mrgloom
Опубликовано

еще некоторые подробности про гомогенные координаты

http://andrewharvey4.wordpress.com/2008/09/29/xyzw-in-opengldirect3d-homogeneous-coordinates/

там ввод такой системы координат объясняется тем, что сдвиг можно будет сделать через умножение матриц(так же по идее мы если запишем перспективное искажение через матричное уравнение в гомогенных координатах, то такое преобразование будет линейным, хотя если мы запишем его через формулу преобразования линейным оно являться не будет), но всё равно это как то кажется недостаточным объяснением, а рекомендованный там курс MATH5785 Geometry я не нашел онлайн.

Поделиться сообщением

Ссылка на сообщение
Поделиться на других сайтах

Smorodov    579

Smorodov
Опубликовано

По поводу преобразований в разных системах, вроде была глава посвященная этой теме в "Multiple view geometry in computer vision".

Там же объяснялось, почему ввели гомогенную СК.

Поделиться сообщением

Ссылка на сообщение
Поделиться на других сайтах

mrgloom    242

mrgloom
Опубликовано

прочитал введение в книге Multiple view geometry in computer vision и всё так же и осталось впечатление, что введение 3-ей координаты не более чем какой то трюк.

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

http://ru.wikipedia.org/wiki/%D0%9F%D0%B5%D1%80%D1%81%D0%BF%D0%B5%D0%BA%D1%82%D0%B8%D0%B2%D0%B0

Поделиться сообщением

Ссылка на сообщение
Поделиться на других сайтах

Smorodov    579

Smorodov
Опубликовано

Перспектива появляется из-за того, что мы видим угловые размеры предметов, а не линейные.

И получается, чем дальше предмет, тем меньше его угловой размер.

По поводу однородных координат, насколько помню, эта единица, дает возможность не терять информацию о векторе, при его трансформациях, приводящих к равенству длины вектора нулю. А так да, это удобный трюк, как и любое представление координат.

Поделиться сообщением

Ссылка на сообщение
Поделиться на других сайтах

Создайте учётную запись или войдите для комментирования

Вы должны быть пользователем, чтобы оставить комментарий

Создать учётную запись

Зарегистрируйтесь для создания учётной записи. Это просто!

Зарегистрировать учётную запись

Войти

Уже зарегистрированы? Войдите здесь.

Войти сейчас
Перейти к списку тем OpenCV

  • Сейчас на странице   0 пользователей

    Нет пользователей, просматривающих эту страницу

×