mrgloom 242 Жалоба Опубликовано October 16, 2012 вот еще вопрос про пики и правильный сдвиг код на матлабе im1= imread('001_001.tif'); im2= imread('001_002.tif'); sz1= size(im1); sz2= size(im2); F1= fftshift(fft2(im1)); F2= fftshift(fft2(im2)); F = F1.*conj(F2); F = F./abs(F); mag = (ifft2(F)); [y x]= find((mag == (max(max(mag))))); img= zeros(sz1(1)+sz2(1),sz1(2)+sz2(2),1); img=cast(img,'uint8'); img(1:sz1(1),1:sz1(2),=im1; img(y:y+sz2(1)-1,x:x+sz2(2)-1,=im2; imshow(img);[/code] всё работает, но если второе изображение сдвигается ниже и левее, а в общем случае мы имеем 4 варианта как одно изображение может наложится на второе. например меняем местами im1 и im2 получаем другое смещение, но как узнать как смещать не имея информации к какому углу приставлять? Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
mrgloom 242 Жалоба Опубликовано October 23, 2012 то что я пытался собрать на коленке уже сделано в последнем матлабе в виде imregister(непонятно какую метрику они используют, но вроде бы даже код из ITK) во всяком случае что то типа Intensity-Based Automatic Image Registration http://www.mathworks.com/help/images/intensity-based-automatic-image-registration.html но проблема в том, что это подходит если надо совместить 2 изображения целиком, а мне же надо совместить 2 изображения с перекрытием в 15% Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
mrgloom 242 Жалоба Опубликовано December 11, 2012 задумался над вопросом как получить корреляцию без fft вот вроде бы код, но вопрос в том как делать "смещение картики" , тупо смещать за пределы ,а остальное заполнять черным? void test(IplImage* src,IplImage* temp) { int w= src->width; int h= src->height; uchar* s_p= (uchar*) src->imageData; uchar* t_p= (uchar*) temp->imageData; __int64 suma=0; __int64 sumb=0; for(int y=0;y<h;++y) { for(int x=0;x<w;++x) { suma+= s_p[x]; sumb+= t_p[x]; } s_p+=w; t_p+=w; } int meana = (int)(suma / (w * h)), meanb = (int)(sumb / (w * h)); s_p= (uchar*) src->imageData; t_p= (uchar*) temp->imageData; __int64 corr = 0, denoma = 0, denomb = 0; for(int y=0;y<h;++y) { for(int x=0;x<w;++x) { int da = (int)s_p[x] - meana, db = (int)t_p[x] - meanb; corr += da * db; denoma += da * da; denomb += db * db; } s_p+=w; t_p+=w; } double peak= corr / sqrt( (double)denoma * denomb); } Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
mrgloom 242 Жалоба Опубликовано December 11, 2012 попробовал так, но что то видимо не так, находит не тот пик. и еще непонятно если искать коэффициент корреляции у абсолютно черного изображения и тестового, то какое значение должно получится? double max_peak=0; Point pt; double per_x=0.5; double per_y=0.5; for(int cy=-h*per_y;cy<h*per_y;++cy) { for(int cx=-w*per_x;cx<w*per_x;++cx) { uchar* s_p= (uchar*) src->imageData; uchar* t_p= (uchar*) temp->imageData; __int64 suma=0; __int64 sumb=0; for(int y=0;y<h;++y) { for(int x=0;x<w;++x) { suma+= s_p[x];//по идее это константно if((y+cy)<h&&(x+cx)<w&& (y+cy)>0&&(x+cx)>0) sumb+= t_p[x+cx]; } s_p+=w; t_p+=w; } int meana = (int)(suma / (w * h)), meanb = (int)(sumb / (w * h)); s_p= (uchar*) src->imageData; t_p= (uchar*) temp->imageData; __int64 corr = 0, denoma = 0, denomb = 0; for(int y=0;y<h;++y) { for(int x=0;x<w;++x) { int da = (int)s_p[x] - meana; int db= 0; if((y+cy)<h&&(x+cx)<w&& (y+cy)>0&&(x+cx)>0) db = (int)t_p[x+cx] - meanb; //тут не очень понят надо заменять нулем? else db= -meanb; corr += da * db; denoma += da * da; denomb += db * db; } s_p+=w; t_p+=w; } double peak= corr / sqrt( (double)denoma * denomb); if (peak>max_peak) { pt.x= cx; pt.y= cy; max_peak=peak; } } } Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
mrgloom 242 Жалоба Опубликовано December 13, 2012 вот так вроде работает, во всяком случае на картинке отсюда http://nashruddin.com/phase-correlation-function-in-opencv.html но на картинках где перекрытие например 15-20% выдает что то не то. void test(IplImage* src,IplImage* temp) { /*IplImage* poc= cvCloneImage(src); cvSmooth(src,src); cvSmooth(temp,temp);*/ //если размеры не одинаковые, то надо дополнять нулями? //делаем пирамиду vector<IplImage> pyr1; pyr_down(src,pyr1); vector<IplImage> pyr2; pyr_down(temp,pyr2); for(int i=0;i<pyr1.size();++i) { std::ostringstream oss; oss <<"pyr/pyr1_"<<i<<".png"; std::cout << oss.str(); cvSaveImage(oss.str().c_str(),&pyr1[i]); oss.str(""); oss <<"pyr/pyr2_"<<i<<".png"; std::cout << oss.str(); cvSaveImage(oss.str().c_str(),&pyr2[i]); } //пробуем по всем масштабам и потом сравним //for(int i=0;i<pyr1.size();++i) { /*src= &pyr1[i]; temp= &pyr2[i];*/ src= &pyr1[pyr1.size()-1]; temp= &pyr2[pyr1.size()-1]; int w= src->width; int h= src->height; //допустим одно изображение неподвижно, а второе сдвигается double max_peak=0; Point pt; double per_x=0.5; double per_y=0.5; for(int cy=-h*per_y;cy<h*per_y;++cy) { for(int cx=-w*per_x;cx<w*per_x;++cx) { uchar* s_p= (uchar*) src->imageData; uchar* t_p= (uchar*) temp->imageData; __int64 suma=0; __int64 sumb=0; for(int y=0;y<h;++y) { for(int x=0;x<w;++x) { suma+= s_p[x];//по идее это константно if((y+cy)<h&&(x+cx)<w&& (y+cy)>0&&(x+cx)>0) { uchar* t= t_p+cy*w; sumb+= t[x+cx]; } } s_p+=w; t_p+=w; } int meana = (int)(suma / (w * h)), meanb = (int)(sumb / (w * h)); s_p= (uchar*) src->imageData; t_p= (uchar*) temp->imageData; __int64 corr = 0, denoma = 0, denomb = 0; for(int y=0;y<h;++y) { for(int x=0;x<w;++x) { int da = (int)s_p[x] - meana; int db= 0; if((y+cy)<h&&(x+cx)<w&& (y+cy)>0&&(x+cx)>0) { uchar* t= t_p+cy*w; db= (int)t[x+cx] - meanb; } else db= -meanb; corr += da * db; denoma += da * da; denomb += db * db; } s_p+=w; t_p+=w; } double peak= corr / sqrt( (double)denoma * denomb); if (peak>max_peak) { pt.x= cx; //теперь вроде работает, только почему то проблема со знаком pt.y= cy; max_peak=peak; } } } cout<<pt.x<<" "<<pt.y<<"\n"; } } еще насчет пирамиды непонятно на каком размере мы должны остановится? void pyr_down(IplImage* src,vector<IplImage>& pyr) { int w= src->width; int h= src->height; pyr.push_back(*src); int dx=128; int dy=128; while(w>dx&&h>dy) { w=w/2; h=h/2; IplImage* temp= cvCreateImage(cvSize(w,h),src->depth,src->nChannels); cvResize(src,temp); pyr.push_back(*temp); } } Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
Smorodov 579 Жалоба Опубликовано December 13, 2012 еще насчет пирамиды непонятно на каком размере мы должны остановится? на 2х2 максимум, а так задается обычно 3-4 итерации, там всего не так и много этажей пирамиды получается. Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
mrgloom 242 Жалоба Опубликовано December 13, 2012 попробовал с dft из opencv всё равно что то не так. Point opencv_dft(Mat& image1,Mat& image2) { int width = getOptimalDFTSize(max(image1.cols,image2.cols)); int height = getOptimalDFTSize(max(image1.rows,image2.rows)); Mat fft1(Size(width,height),CV_32F,Scalar(0)); Mat fft2(Size(width,height),CV_32F,Scalar(0)); for(int j=0; j<image1.rows; j++) for(int i=0; i<image1.cols; i++) fft1.at<float>(j,i) = image1.at<unsigned char>(j,i); for(int j=0; j<image2.rows; j++) for(int i=0; i<image2.cols; i++) fft2.at<float>(j,i) = image2.at<unsigned char>(j,i); dft(fft1,fft1,0,image1.rows); dft(fft2,fft2,0,image2.rows); mulSpectrums(fft1,fft2,fft1,0,true); idft(fft1,fft1); double maxVal; Point maxLoc; minMaxLoc(fft1,NULL,&maxVal,NULL,&maxLoc); int resX = (maxLoc.x<width/2) ? (maxLoc.x) : (maxLoc.x-width); int resY = (maxLoc.y<height/2) ? (maxLoc.y) : (maxLoc.y-height); return Point(resX,resY); } Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
Smorodov 579 Жалоба Опубликовано December 13, 2012 Так тут нужно комплексные изображения подавать наверное (иначе результат в каком-то непонятном формате), и четверти переставлять. Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
mrgloom 242 Жалоба Опубликовано December 13, 2012 в SVN версии opencv нашел еще 1 версию, но она несколько сложнее, но все равно не дает нужный результат) #include "stdafx.h" #include "cv.h" #include "highgui.h" using namespace cv; using namespace std; void magSpectrums( InputArray _src, OutputArray _dst) { Mat src = _src.getMat(); int depth = src.depth(), cn = src.channels(), type = src.type(); int rows = src.rows, cols = src.cols; int j, k; CV_Assert( type == CV_32FC1 || type == CV_32FC2 || type == CV_64FC1 || type == CV_64FC2 ); if(src.depth() == CV_32F) _dst.create( src.rows, src.cols, CV_32FC1 ); else _dst.create( src.rows, src.cols, CV_64FC1 ); Mat dst = _dst.getMat(); dst.setTo(0);//Mat elements are not equal to zero by default! bool is_1d = (rows == 1 || (cols == 1 && src.isContinuous() && dst.isContinuous())); if( is_1d ) cols = cols + rows - 1, rows = 1; int ncols = cols*cn; int j0 = cn == 1; int j1 = ncols - (cols % 2 == 0 && cn == 1); if( depth == CV_32F ) { const float* dataSrc = (const float*)src.data; float* dataDst = (float*)dst.data; size_t stepSrc = src.step/sizeof(dataSrc[0]); size_t stepDst = dst.step/sizeof(dataDst[0]); if( !is_1d && cn == 1 ) { for( k = 0; k < (cols % 2 ? 1 : 2); k++ ) { if( k == 1 ) dataSrc += cols - 1, dataDst += cols - 1; dataDst[0] = dataSrc[0]*dataSrc[0]; if( rows % 2 == 0 ) dataDst[(rows-1)*stepDst] = dataSrc[(rows-1)*stepSrc]*dataSrc[(rows-1)*stepSrc]; for( j = 1; j <= rows - 2; j += 2 ) { dataDst[j*stepDst] = (float)((double)dataSrc[j*stepSrc]*dataSrc[j*stepSrc] + (double)dataSrc[(j+1)*stepSrc]*dataSrc[(j+1)*stepSrc]); } if( k == 1 ) dataSrc -= cols - 1, dataDst -= cols - 1; } } for( ; rows--; dataSrc += stepSrc, dataDst += stepDst ) { if( is_1d && cn == 1 ) { dataDst[0] = dataSrc[0]*dataSrc[0]; if( cols % 2 == 0 ) dataDst[j1] = dataSrc[j1]*dataSrc[j1]; } for( j = j0; j < j1; j += 2 ) { dataDst[j] = (float)((double)dataSrc[j]*dataSrc[j] + (double)dataSrc[j+1]*dataSrc[j+1]); } } } else { const double* dataSrc = (const double*)src.data; double* dataDst = (double*)dst.data; size_t stepSrc = src.step/sizeof(dataSrc[0]); size_t stepDst = dst.step/sizeof(dataDst[0]); if( !is_1d && cn == 1 ) { for( k = 0; k < (cols % 2 ? 1 : 2); k++ ) { if( k == 1 ) dataSrc += cols - 1, dataDst += cols - 1; dataDst[0] = dataSrc[0]*dataSrc[0]; if( rows % 2 == 0 ) dataDst[(rows-1)*stepDst] = dataSrc[(rows-1)*stepSrc]*dataSrc[(rows-1)*stepSrc]; for( j = 1; j <= rows - 2; j += 2 ) { dataDst[j*stepDst] = dataSrc[j*stepSrc]*dataSrc[j*stepSrc] + dataSrc[(j+1)*stepSrc]*dataSrc[(j+1)*stepSrc]; } if( k == 1 ) dataSrc -= cols - 1, dataDst -= cols - 1; } } for( ; rows--; dataSrc += stepSrc, dataDst += stepDst ) { if( is_1d && cn == 1 ) { dataDst[0] = dataSrc[0]*dataSrc[0]; if( cols % 2 == 0 ) dataDst[j1] = dataSrc[j1]*dataSrc[j1]; } for( j = j0; j < j1; j += 2 ) { dataDst[j] = dataSrc[j]*dataSrc[j] + dataSrc[j+1]*dataSrc[j+1]; } } } // do batch sqrt to use SSE optimizations... cv::sqrt(dst, dst); } void divSpectrums( InputArray _srcA, InputArray _srcB, OutputArray _dst, int flags, bool conjB) { Mat srcA = _srcA.getMat(), srcB = _srcB.getMat(); int depth = srcA.depth(), cn = srcA.channels(), type = srcA.type(); int rows = srcA.rows, cols = srcA.cols; int j, k; CV_Assert( type == srcB.type() && srcA.size() == srcB.size() ); CV_Assert( type == CV_32FC1 || type == CV_32FC2 || type == CV_64FC1 || type == CV_64FC2 ); _dst.create( srcA.rows, srcA.cols, type ); Mat dst = _dst.getMat(); bool is_1d = (flags & DFT_ROWS) || (rows == 1 || (cols == 1 && srcA.isContinuous() && srcB.isContinuous() && dst.isContinuous())); if( is_1d && !(flags & DFT_ROWS) ) cols = cols + rows - 1, rows = 1; int ncols = cols*cn; int j0 = cn == 1; int j1 = ncols - (cols % 2 == 0 && cn == 1); if( depth == CV_32F ) { const float* dataA = (const float*)srcA.data; const float* dataB = (const float*)srcB.data; float* dataC = (float*)dst.data; float eps = FLT_EPSILON; // prevent div0 problems size_t stepA = srcA.step/sizeof(dataA[0]); size_t stepB = srcB.step/sizeof(dataB[0]); size_t stepC = dst.step/sizeof(dataC[0]); if( !is_1d && cn == 1 ) { for( k = 0; k < (cols % 2 ? 1 : 2); k++ ) { if( k == 1 ) dataA += cols - 1, dataB += cols - 1, dataC += cols - 1; dataC[0] = dataA[0] / dataB[0]; if( rows % 2 == 0 ) dataC[(rows-1)*stepC] = dataA[(rows-1)*stepA] / dataB[(rows-1)*stepB]; if( !conjB ) for( j = 1; j <= rows - 2; j += 2 ) { double denom = (double)dataB[j*stepB]*dataB[j*stepB] + (double)dataB[(j+1)*stepB]*dataB[(j+1)*stepB] + (double)eps; double re = (double)dataA[j*stepA]*dataB[j*stepB] + (double)dataA[(j+1)*stepA]*dataB[(j+1)*stepB]; double im = (double)dataA[(j+1)*stepA]*dataB[j*stepB] - (double)dataA[j*stepA]*dataB[(j+1)*stepB]; dataC[j*stepC] = (float)(re / denom); dataC[(j+1)*stepC] = (float)(im / denom); } else for( j = 1; j <= rows - 2; j += 2 ) { double denom = (double)dataB[j*stepB]*dataB[j*stepB] + (double)dataB[(j+1)*stepB]*dataB[(j+1)*stepB] + (double)eps; double re = (double)dataA[j*stepA]*dataB[j*stepB] - (double)dataA[(j+1)*stepA]*dataB[(j+1)*stepB]; double im = (double)dataA[(j+1)*stepA]*dataB[j*stepB] + (double)dataA[j*stepA]*dataB[(j+1)*stepB]; dataC[j*stepC] = (float)(re / denom); dataC[(j+1)*stepC] = (float)(im / denom); } if( k == 1 ) dataA -= cols - 1, dataB -= cols - 1, dataC -= cols - 1; } } for( ; rows--; dataA += stepA, dataB += stepB, dataC += stepC ) { if( is_1d && cn == 1 ) { dataC[0] = dataA[0] / dataB[0]; if( cols % 2 == 0 ) dataC[j1] = dataA[j1] / dataB[j1]; } if( !conjB ) for( j = j0; j < j1; j += 2 ) { double denom = (double)(dataB[j]*dataB[j] + dataB[j+1]*dataB[j+1] + eps); double re = (double)(dataA[j]*dataB[j] + dataA[j+1]*dataB[j+1]); double im = (double)(dataA[j+1]*dataB[j] - dataA[j]*dataB[j+1]); dataC[j] = (float)(re / denom); dataC[j+1] = (float)(im / denom); } else for( j = j0; j < j1; j += 2 ) { double denom = (double)(dataB[j]*dataB[j] + dataB[j+1]*dataB[j+1] + eps); double re = (double)(dataA[j]*dataB[j] - dataA[j+1]*dataB[j+1]); double im = (double)(dataA[j+1]*dataB[j] + dataA[j]*dataB[j+1]); dataC[j] = (float)(re / denom); dataC[j+1] = (float)(im / denom); } } } else { const double* dataA = (const double*)srcA.data; const double* dataB = (const double*)srcB.data; double* dataC = (double*)dst.data; double eps = DBL_EPSILON; // prevent div0 problems size_t stepA = srcA.step/sizeof(dataA[0]); size_t stepB = srcB.step/sizeof(dataB[0]); size_t stepC = dst.step/sizeof(dataC[0]); if( !is_1d && cn == 1 ) { for( k = 0; k < (cols % 2 ? 1 : 2); k++ ) { if( k == 1 ) dataA += cols - 1, dataB += cols - 1, dataC += cols - 1; dataC[0] = dataA[0] / dataB[0]; if( rows % 2 == 0 ) dataC[(rows-1)*stepC] = dataA[(rows-1)*stepA] / dataB[(rows-1)*stepB]; if( !conjB ) for( j = 1; j <= rows - 2; j += 2 ) { double denom = dataB[j*stepB]*dataB[j*stepB] + dataB[(j+1)*stepB]*dataB[(j+1)*stepB] + eps; double re = dataA[j*stepA]*dataB[j*stepB] + dataA[(j+1)*stepA]*dataB[(j+1)*stepB]; double im = dataA[(j+1)*stepA]*dataB[j*stepB] - dataA[j*stepA]*dataB[(j+1)*stepB]; dataC[j*stepC] = re / denom; dataC[(j+1)*stepC] = im / denom; } else for( j = 1; j <= rows - 2; j += 2 ) { double denom = dataB[j*stepB]*dataB[j*stepB] + dataB[(j+1)*stepB]*dataB[(j+1)*stepB] + eps; double re = dataA[j*stepA]*dataB[j*stepB] - dataA[(j+1)*stepA]*dataB[(j+1)*stepB]; double im = dataA[(j+1)*stepA]*dataB[j*stepB] + dataA[j*stepA]*dataB[(j+1)*stepB]; dataC[j*stepC] = re / denom; dataC[(j+1)*stepC] = im / denom; } if( k == 1 ) dataA -= cols - 1, dataB -= cols - 1, dataC -= cols - 1; } } for( ; rows--; dataA += stepA, dataB += stepB, dataC += stepC ) { if( is_1d && cn == 1 ) { dataC[0] = dataA[0] / dataB[0]; if( cols % 2 == 0 ) dataC[j1] = dataA[j1] / dataB[j1]; } if( !conjB ) for( j = j0; j < j1; j += 2 ) { double denom = dataB[j]*dataB[j] + dataB[j+1]*dataB[j+1] + eps; double re = dataA[j]*dataB[j] + dataA[j+1]*dataB[j+1]; double im = dataA[j+1]*dataB[j] - dataA[j]*dataB[j+1]; dataC[j] = re / denom; dataC[j+1] = im / denom; } else for( j = j0; j < j1; j += 2 ) { double denom = dataB[j]*dataB[j] + dataB[j+1]*dataB[j+1] + eps; double re = dataA[j]*dataB[j] - dataA[j+1]*dataB[j+1]; double im = dataA[j+1]*dataB[j] + dataA[j]*dataB[j+1]; dataC[j] = re / denom; dataC[j+1] = im / denom; } } } } void fftShift(InputOutputArray _out) { Mat out = _out.getMat(); vector<Mat> planes; split(out, planes); int xMid = out.cols >> 1; int yMid = out.rows >> 1; for(size_t i = 0; i < planes.size(); i++) { // perform quadrant swaps... Mat tmp; Mat q0(planes[i], Rect(0, 0, xMid, yMid)); Mat q1(planes[i], Rect(xMid, 0, xMid, yMid)); Mat q2(planes[i], Rect(0, yMid, xMid, yMid)); Mat q3(planes[i], Rect(xMid, yMid, xMid, yMid)); q0.copyTo(tmp); q3.copyTo(q0); tmp.copyTo(q3); q1.copyTo(tmp); q2.copyTo(q1); tmp.copyTo(q2); } merge(planes, out); } //даёт субпиксельную точность? Point2d weightedCentroid(InputArray _src, cv::Point peakLocation, cv::Size weightBoxSize) { Mat src = _src.getMat(); int type = src.type(); CV_Assert( type == CV_32FC1 || type == CV_64FC1 ); int minr = peakLocation.y - (weightBoxSize.height >> 1); int maxr = peakLocation.y + (weightBoxSize.height >> 1); int minc = peakLocation.x - (weightBoxSize.width >> 1); int maxc = peakLocation.x + (weightBoxSize.width >> 1); Point2d centroid; double sumIntensity = 0.0; // clamp the values to min and max if needed. if(minr < 0) { minr = 0; } if(minc < 0) { minc = 0; } if(maxr > src.rows - 1) { maxr = src.rows - 1; } if(maxc > src.cols - 1) { maxc = src.cols - 1; } if(type == CV_32FC1) { const float* dataIn = (const float*)src.data; dataIn += minr*src.cols; for(int y = minr; y <= maxr; y++) { for(int x = minc; x <= maxc; x++) { centroid.x += (double)x*dataIn[x]; centroid.y += (double)y*dataIn[x]; sumIntensity += (double)dataIn[x]; } dataIn += src.cols; } } else { const double* dataIn = (const double*)src.data; dataIn += minr*src.cols; for(int y = minr; y <= maxr; y++) { for(int x = minc; x <= maxc; x++) { centroid.x += (double)x*dataIn[x]; centroid.y += (double)y*dataIn[x]; sumIntensity += dataIn[x]; } dataIn += src.cols; } } sumIntensity += DBL_EPSILON; // prevent div0 problems... centroid.x /= sumIntensity; centroid.y /= sumIntensity; return centroid; } cv::Point2d phaseCorrelate(InputArray _src1, InputArray _src2, InputArray _window) { Mat src1 = _src1.getMat(); Mat src2 = _src2.getMat(); Mat window = _window.getMat(); CV_Assert( src1.type() == src2.type()); CV_Assert( src1.type() == CV_32FC1 || src1.type() == CV_64FC1 ); CV_Assert( src1.size == src2.size); if(!window.empty()) { CV_Assert( src1.type() == window.type()); CV_Assert( src1.size == window.size); } int M = getOptimalDFTSize(src1.rows); int N = getOptimalDFTSize(src1.cols); Mat padded1, padded2, paddedWin; //дополняем нулями. if(M != src1.rows || N != src1.cols) { copyMakeBorder(src1, padded1, 0, M - src1.rows, 0, N - src1.cols, BORDER_CONSTANT, Scalar::all(0)); copyMakeBorder(src2, padded2, 0, M - src2.rows, 0, N - src2.cols, BORDER_CONSTANT, Scalar::all(0)); if(!window.empty()) { copyMakeBorder(window, paddedWin, 0, M - window.rows, 0, N - window.cols, BORDER_CONSTANT, Scalar::all(0)); } } else { padded1 = src1; padded2 = src2; paddedWin = window; } Mat FFT1, FFT2, P, Pm, C; //можно и без хаминг виндоу? хотя что то на результат это не влияет. // perform window multiplication if available if(!paddedWin.empty()) { // apply window to both images before proceeding... multiply(paddedWin, padded1, padded1); multiply(paddedWin, padded2, padded2); } // execute phase correlation equation // Reference: http://en.wikipedia.org/wiki/Phase_correlation dft(padded1, FFT1, DFT_REAL_OUTPUT); //непонятные настройки? dft(padded2, FFT2, DFT_REAL_OUTPUT); mulSpectrums(FFT1, FFT2, P, 0, true); magSpectrums(P, Pm); divSpectrums(P, Pm, C, 0, false); // FF* / |FF*| (phase correlation equation completed here...) idft(C, C); // gives us the nice peak shift location... fftShift(C); // shift the energy to the center of the frame. // locate the highest peak Point peakLoc; minMaxLoc(C, NULL, NULL, NULL, &peakLoc); // get the phase shift with sub-pixel accuracy, 5x5 window seems about right here... //может это наоборот портит? Point2d t; //t = weightedCentroid(C, peakLoc, Size(5, 5)); t= peakLoc; // adjust shift relative to image center... Point2d center((double)src1.cols / 2.0, (double)src1.rows / 2.0); return (center - t); } void createHanningWindow(OutputArray _dst, cv::Size winSize, int type) { CV_Assert( type == CV_32FC1 || type == CV_64FC1 ); _dst.create(winSize, type); Mat dst = _dst.getMat(); int rows = dst.rows; int cols = dst.cols; int step = dst.step/dst.elemSize1(); if(dst.depth() == CV_32F) { float* dstData = dst.ptr<float>(); for(int i = 0; i < rows; i++) { double wr = 0.5 * (1.0f - cos(2.0f * CV_PI * (double)i / (double)(rows - 1))); for(int j = 0; j < cols; j++) { double wc = 0.5 * (1.0f - cos(2.0f * CV_PI * (double)j / (double)(cols - 1))); dstData[i*step + j] = (float)(wr * wc); } } } else { double* dstData = dst.ptr<double>(); for(int i = 0; i < rows; i++) { double wr = 0.5 * (1.0 - cos(2.0 * CV_PI * (double)i / (double)(rows - 1))); for(int j = 0; j < cols; j++) { double wc = 0.5 * (1.0 - cos(2.0 * CV_PI * (double)j / (double)(cols - 1))); dstData[i*step + j] = wr * wc; } } } // perform batch sqrt for SSE performance gains cv::sqrt(dst, dst); } int _tmain(int argc, _TCHAR* argv[]) { //Mat r1 = Mat::ones(Size(128, 128), CV_64F); //Mat r2 = Mat::ones(Size(128, 128), CV_64F); //double expectedShiftX = -10.0; //double expectedShiftY = -20.0; //// draw 10x10 rectangles @ (100, 100) and (90, 80) should see ~(-10, -20) shift here... //cv::rectangle(r1, Point(100, 100), Point(110, 110), Scalar(0, 0, 0), CV_FILLED); //cv::rectangle(r2, Point(90, 80), Point(100, 90), Scalar(0, 0, 0), CV_FILLED); /*Mat image1= imread("pyr/im1.PNG",0); Mat image2= imread("pyr/im2.PNG",0);*/ /*Mat image1= imread("pyr/small_1.jpg",0); Mat image2= imread("pyr/small_2.jpg",0);*/ Mat image1= imread("pyr/001_001.tif",0); Mat image2= imread("pyr/001_002.tif",0); Mat r1; image1.convertTo(r1,CV_64F/*C1*/); Mat r2; image2.convertTo(r2,CV_64F/*C1*/); Mat hann; createHanningWindow(hann, r1.size(), CV_64F/*C1*/); Point2d phaseShift1 = phaseCorrelate(r1, r2, hann); Point2d phaseShift2 = phaseCorrelate(r2, r1, hann); return 0; } Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
Smorodov 579 Жалоба Опубликовано December 13, 2012 Посмотрите книжку: Lindblad T., Kinser J.M. Image Processing Using Pulse-Coupled Neural Networks там есть немного по этой теме. Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
mrgloom 242 Жалоба Опубликовано December 14, 2012 Посмотрите книжку: Lindblad T., Kinser J.M. Image Processing Using Pulse-Coupled Neural Networks там есть немного по этой теме. там что то совсем в другую степь. посмотрел исходники matchtemplate код я до конца не понял. там почему то используется block-wise cross correlation http://simon.rozman.si/computers/dsp/block-wise-cross-correlation In some applications, one of the input vectors (signals) is much longer than another, or even of infinite length. For example: searching a transmitted signal’s echo in the received input stream in radar or sonar, or detecting an RDS signal inside a continuous radio transmission. In such applications one must compute the cross correlation of searched pattern and input stream in a block-by-block manner. типа если у нас очень большой сигнал и маленький темплейт, но все равно непонятно зачем нам N маленьких, если мы можем сделать 1 раз для изображения и дополненного темплейта. Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
mrgloom 242 Жалоба Опубликовано December 14, 2012 Так тут нужно комплексные изображения подавать наверное (иначе результат в каком-то непонятном формате), и четверти переставлять. http://en.wikipedia.org/wiki/Discrete_Fourier_transform#The_real-input_DFT вроде как определено и только для real части. а перестановка частей только для отображения и применеия фильтров вроде. ну вообщем судя по тому что c fftw и со старым интерфейсом opecnv работает, с точностью до каких то смещений, которые я не понимаю откуда возникают. int resX = (maxLoc.x<width/2) ? (maxLoc.x) : (maxLoc.x-width); int resY = (maxLoc.y<height/2) ? (maxLoc.y) : (maxLoc.y-height); пробелема либо в том что всё таки надо еще и комплексную занулёную часть подавать, либо в том что mulSpectrums делает что то не то. dft(fft1,fft1,0,image1.rows); dft(fft2,fft2,0,image2.rows); mulSpectrums(fft1,fft2,fft1,0,true); idft(fft1,fft1); конволюция всё таки плотно связана с кореляцией почти одно и то же. (возможно в этом ошибка? а если нет то где разница) void convolveDFT(const Mat& A, const Mat& B, Mat& C) { // reallocate the output array if needed C.create(abs(A.rows - B.rows)+1, abs(A.cols - B.cols)+1, A.type()); Size dftSize; // compute the size of DFT transform dftSize.width = getOptimalDFTSize(A.cols + B.cols - 1); dftSize.height = getOptimalDFTSize(A.rows + B.rows - 1); // allocate temporary buffers and initialize them with 0's Mat tempA(dftSize, A.type(), Scalar::all(0)); Mat tempB(dftSize, B.type(), Scalar::all(0)); // copy A and B to the top-left corners of tempA and tempB, respectively Mat roiA(tempA, Rect(0,0,A.cols,A.rows)); A.copyTo(roiA); Mat roiB(tempB, Rect(0,0,B.cols,B.rows)); B.copyTo(roiB); // now transform the padded A & B in-place; // use "nonzeroRows" hint for faster processing dft(tempA, tempA, 0, A.rows); dft(tempB, tempB, 0, B.rows); // multiply the spectrums; // the function handles packed spectrum representations well mulSpectrums(tempA, tempB, tempA); // transform the product back from the frequency domain. // Even though all the result rows will be non-zero, // we need only the first C.rows of them, and thus we // pass nonzeroRows == C.rows dft(tempA, tempA, DFT_INVERSE + DFT_SCALE, C.rows); // now copy the result back to C. tempA(Rect(0, 0, C.cols, C.rows)).copyTo©; // all the temporary buffers will be deallocated automatically } Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
mrgloom 242 Жалоба Опубликовано December 17, 2012 использовал версию опенцв 2.4.3 она работает на некоторых изображениях, а например на изображении где 20% перекрытия не всегда выдает не тот пик, что то типа (w,h)-peak или -peak. не пойму с чем это всё связано, какие то краевые проблемы или связанные с цикличностью? когда изображения разного размера я добавляю нулями, но в зависимости от размера (см. в коде) получаются разные результаты. #include "stdafx.h" #include <opencv.hpp> using namespace cv; using namespace std; int _tmain(int argc, _TCHAR* argv[]) { Mat image1= imread("1.png",0); Mat image2= imread("2.png",0); Mat r1; image1.convertTo(r1,CV_64F); Mat r2; image2.convertTo(r2,CV_64F); Point2d phaseShift1; Point2d phaseShift2; //check size if(r1.cols!=r2.cols||r1.rows!=r2.rows) { //image1 size 883x646 //image2 size 882x599 //gives peak -423 -1 (423 1) int n_cols= 1024; int n_rows= 1024; //gives peak 224 -1 (-224 1) int n_cols= max(r1.cols,r2.cols); int n_rows= max(r1.rows,r2.rows); Mat r1_pad; copyMakeBorder(r1,r1_pad,0,n_rows-r1.rows,0,n_cols-r1.cols, BORDER_CONSTANT, Scalar::all(0)); Mat r2_pad; copyMakeBorder(r2,r2_pad,0,n_rows-r2.rows,0,n_cols-r2.cols, BORDER_CONSTANT, Scalar::all(0)); //imwrite("b1.png",r1_pad); //imwrite("b2.png",r2_pad); //Mat hann; //createHanningWindow(hann, r1_pad.size(), CV_64F); phaseShift1 = phaseCorrelate(r1_pad, r2_pad/*, hann*/); phaseShift2 = phaseCorrelate(r2_pad, r1_pad/*, hann*/); int f=0; } else { //Mat hann; //createHanningWindow(hann, r1.size(), CV_64F); phaseShift1 = phaseCorrelate(r1, r2/*, hann*/); phaseShift2 = phaseCorrelate(r2, r1/*, hann*/); } return 0; } Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
mrgloom 242 Жалоба Опубликовано December 19, 2012 Опять же вопрос, если работаем с темплейтом(он меньше), то всё понятно двигаем темплейт по изображению считаем совпадение. а что если изображения только частично перекрываются и надо найти перекрытие? вот я написал просто попиксельную разницу и корреляцию для теста. void plain_corr(Mat_<uchar> src, Mat_<uchar> tmp) { int w_s= src.cols; int h_s= src.rows; int w_t= tmp.cols; int h_t= tmp.rows; //просто разность double min_peak=INT_MAX; Point pt; for(int y=0;y<h_s-h_t+1;++y) { for(int x=0;x<w_s-w_t+1;++x) { __int64 sum=0; for(int cy=0;cy<h_t;++cy) { for(int cx=0;cx<w_t;++cx) { sum+= abs(src[y+cy][x+cx]-tmp[cy][cx]); } } double diff= (double)sum/(w_t*h_t); if (diff<min_peak) { pt.x= x; pt.y= y; min_peak=diff; } } } cout<<pt.x<<" "<<pt.y<<"\n"; double max_peak=0; Point pt; for(int y=0;y<h_s-h_t+1;++y) { for(int x=0;x<w_s-w_t+1;++x) { __int64 suma=0; __int64 sumb=0; for(int cy=0;cy<h_t;++cy) { for(int cx=0;cx<w_t;++cx) { suma+= src[y+cy][x+cx]; sumb+= tmp[cy][cx]; } } int meana = (int)(suma/(w_t*h_t)), meanb = (int)(sumb/(w_t*h_t)); __int64 corr = 0, denoma = 0, denomb = 0; for(int cy=0;cy<h_t;++cy) { for(int cx=0;cx<w_t;++cx) { int da= src[y+cy][x+cx] - meana; int db= tmp[cy][cx] - meanb; corr += da * db; denoma += da * da; denomb += db * db; } } double peak= corr / sqrt( (double)denoma * denomb); if (peak>max_peak) { pt.x= x; pt.y= y; max_peak=peak; } } } cout<<pt.x<<" "<<pt.y<<"\n"; } Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
Smorodov 579 Жалоба Опубликовано December 26, 2012 Кусочек из книжки "GPU Computing Gems Emerald Edition" Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
mrgloom 242 Жалоба Опубликовано December 27, 2012 не понял что за kthLaw, но похоже это что то типа окна. например применяют Hamming window, но на практике это делает и лучше и хуже. кстати понял как делать если у нас не маленький темплейт, а надо найти перекрытие 2-х изображений. надо пройтись вторым изображением по первому полностью и считать корреляцию в области пересечения(там получается нормируется на кол-во пикселей и всё норм), но тут может быть проблема, когда область пересечения довольно маленькая, ну предельный случай 1 пиксель, т.е. может быть очень большой коэффициент корреляции для маленькой области. кстати еще по поводу пирамиды: получается пирамиду мы можем использовать только в том случае, если не используем FFT. ну и если темплейт соответственно не маленький. еще исходя из вышесказанного про мин. размер в пикселях получается, что и мин. размер изображения в пирамиде тоже должен быть ограничен. Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
mrgloom 242 Жалоба Опубликовано April 30, 2013 http://www.lfd.uci.edu/~gohlke/code/imreg.py.html версия на питоне обещающая Implements an FFT-based technique for translation, rotation and scale-invariant image registration по References ---------- (1) An FFT-based technique for translation, rotation and scale-invariant image registration. BS Reddy, BN Chatterji. IEEE Transactions on Image Processing, 5, 1266-1271, 1996 (2) An IDL/ENVI implementation of the FFT-based algorithm for automatic image registration. H Xiea, N Hicksa, GR Kellera, H Huangb, V Kreinovich. Computers & Geosciences, 29, 1045-1055, 2003. (3) Image Registration Using Adaptive Polar Transform. R Matungka, YF Zheng, RL Ewing. IEEE Transactions on Image Processing, 18(10), 2009. только там чтение файлов .hdr, но можно вроде и с обычными. возможно способ кривой и можно и без Image import Image im0 = numpy.asarray(Image.open("001.png")) протестировал, но там какие то пробелемы опять связанные с перескакиванием через сторону def translation(im0, im1): """Return translation vector to register images.""" shape = im0.shape f0 = fft2(im0) f1 = fft2(im1) ir = abs(ifft2((f0 * f1.conjugate()) / (abs(f0) * abs(f1)))) t0, t1 = numpy.unravel_index(numpy.argmax(ir), shape) if t0 > shape[0] // 2: t0 -= shape[0] if t1 > shape[1] // 2: t1 -= shape[1] return [t0, t1] всмыстле с этим if t0 > shape[0] // 2: t0 -= shape[0] if t1 > shape[1] // 2: t1 -= shape[1] и я так и не понял что делает unravel_index точнее понятно она вытаскивает максимумы, но как работает/применяется непонятно. Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
mrgloom 242 Жалоба Опубликовано May 16, 2013 допустим хотим произвести регистрацию изображений(частично наложить\пересечь 2 изображения) как начальное приближение берём пик который выдал FFT, затем хотим произвести аффинную трансформацию изображения(имеем 6 параметров), чтобы изображение наложилось оптимально. допустим берем минимайзер http://www.chokkan.org/software/liblbfgs/ и сразу возникает вопрос как посчитать градиент? (если функция у нас задана как сумма модулей разностей пикселей в области пересечения или как пик FFT) Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
Smorodov 579 Жалоба Опубликовано May 16, 2013 Численным методом: df/dx=(f(x+epsilon)-f(x-epsilon))/(2*epsilon); если нужны частные производные, то epsilon добавляется к одной компоненте (к той, по которой дифференцируем). Последовательно прогнав по всем компонентам x, получим градиент. ЗЫ: Из оптимизаторов есть еще один, который мне удалось без проблем собрать в VS2010: http://ab-initio.mit.edu/wiki/index.php/NLopt_on_Windows Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
Nuzhny 243 Жалоба Опубликовано May 16, 2013 и сразу возникает вопрос как посчитать градиент? (если функция у нас задана как сумма модулей разностей пикселей в области пересечения или как пик FFT) Не пробовал findTransformECC - новую функцию, которая находит разные типы преообразований между картинками? Я погонял на паре пример, ничего так результаты. О самом методе написано в нескольких книгах и статьях, в основном при описании стабилизации видео. Впрочем, в исходниках есть отсылка к статье, которую и имплементировали. Там как раз итерациионно находится минимум градиента между картинками (это в двух словах). 1 Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
mrgloom 242 Жалоба Опубликовано May 16, 2013 хотя возможно для вычисления градиента(якобиана), я должен просто посчитать частные производные, как тут http://ru.wikipedia.org/wiki/%D0%A7%D0%B0%D1%81%D1%82%D0%BD%D0%B0%D1%8F_%D0%BF%D1%80%D0%BE%D0%B8%D0%B7%D0%B2%D0%BE%D0%B4%D0%BD%D0%B0%D1%8F только как выбрать дельта x? просто маленькое типа 0.001 ? http://www.machinelearning.ru/wiki/index.php?title=%D0%92%D1%8B%D1%87%D0%B8%D1%81%D0%BB%D0%B5%D0%BD%D0%B8%D0%B5_%D0%BC%D0%B0%D1%82%D1%80%D0%B8%D1%86_%D0%AF%D0%BA%D0%BE%D0%B1%D0%B8_%D0%B8_%D0%93%D0%B5%D1%81%D1%81%D0%B5 тут это называется метод конечных разностей ну вроде тоже что и вы написали, так в чем проблема подавать в функцию минимизации только F и epsilon(т.е. почему они хотят чтобы пользователь это сверху их функции считал)? findTransformECC это она официально появилась? в 2.4.3 у себя не нашел. Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
Nuzhny 243 Жалоба Опубликовано May 16, 2013 У меня есть последняя релизная версия OpenCV для работы. И одна из репозитория для опытов, там уже есть эта функция. И пример её использования тоже. Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
Smorodov 579 Жалоба Опубликовано May 16, 2013 только как выбрать дельта x? просто маленькое типа 0.001 ? Здесь две крайности: 1) Слишком маленькое значение может вообще не дать изменения из-за дискретности изображения. И производная получится равной нулю. 2) Слишком большое уменьшит точность аппроксимации. почему они хотят чтобы пользователь это сверху их функции считал Аналитическое выражение производной точнее и намного быстрее считается. В NLOpt, например, есть методы, не требующие от пользователя градиентов. Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
mrgloom 242 Жалоба Опубликовано May 16, 2013 кстати попробовал NLOpt собрать под питон D:\nlopt-2.3-dll>python setup.py in stall running install running build running build_py creating build creating build\lib.win32-2.7 copying nlopt.py -> build\lib.win32-2.7 running build_ext building '_nlopt' extension creating build\temp.win32-2.7 creating build\temp.win32-2.7\Release C:\Program Files\Microsoft Visual Studio 9.0\VC\BIN\cl.exe /c /nologo /Ox /MD /W 3 /GS- /DNDEBUG -I. -IC:\Python27\lib\site-packages\numpy\core\include -IC:\Pyth on27\include -IC:\Python27\PC /Tpnlopt-python.cpp /Fobuild\temp.win32-2.7\Releas e\nlopt-python.obj nlopt-python.cpp d:\nlopt-2.3-dll\nlopt.hpp(187) : w arning C4530: C++ exception handler used, but unwind semantics are not enabled. Specify /EHsc nlopt-python.cpp(8330) : warning C4101: '_e' : unreferenced local variable nlopt-python.cpp(8335) : warning C4101: '_e' : unreferenced local variable nlopt-python.cpp(8399) : warning C4101: '_e' : unreferenced local variable nlopt-python.cpp(8404) : warning C4101: '_e' : unreferenced local variable C:\Program Files\Microsoft Visual Studio 9.0\VC\BIN\link.exe /DLL /nologo /INCRE MENTAL:NO /LIBPATH:C:\Python27\libs /LIBPATH:C:\Python27\PCbuild libnlopt-0.lib /EXPORT:init_nlopt build\temp.win32-2.7\Release\nlopt-python.obj /OUT:build\lib. win32-2.7\_nlopt.pyd /IMPLIB:build\temp.win32-2.7\Release\_nlopt.lib /MANIFESTFI LE:build\temp.win32-2.7\Release\_nlopt.pyd.manifest LINK : fatal error LNK1181: cannot open input file 'libnlopt-0.lib' error: command '"C:\Program Files\Microsoft Visual Studio 9.0\VC\BIN\link.exe"' failed with exit status 1181 видимо не может найти библиотеку? которую предполагается получить из dll? а при использовании CMakeLists.txt написало CMake Error at CMakeLists.txt:60 (FILE): file Internal CMake error when trying to open file: D:/nlopt-2.3-dll/configure.ac for reading. NLOPT: Version number ..0 found in configure.ac Configuring incomplete, errors occurred! Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
Smorodov 579 Жалоба Опубликовано May 16, 2013 Для VS там есть пара отдельных файлов, Alternatively, you can use the following files, provided by Benoit Scherrer (benoitscherrer ατ gmail.com), to compile NLopt from source on Windows (with the Microsoft compiler) using CMake: CMakeLists.txt and config.cmake.h.in Да и бинарники под питон там есть: http://www.lfd.uci.edu/~gohlke/pythonlibs/#nlopt 1 Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах