(LOG)Laplacian of Guassian & (DOH)Determinant of Hessian 斑点检测

1. 原理

考虑一个一维信号，该信号在t=1000时发生了阶跃，那么如何检测这个阶跃的位置呢？

图1分别用1阶高斯导数和2阶高斯导数对该信号进行卷积

图1-1. 一阶高斯导数（左）、二阶高斯导数（右）与原始信号卷积

可以看到，一阶高斯导的卷积结果在t=1000产生了一个极值，二阶高斯导的卷积在该位置是个过零点

LOG就是利用二阶高斯导数（也叫拉普拉斯变换）与原始信号（图像）卷积，通过检测局部极值获得角点。

需要注意的是：不同方差的二阶高斯导数得到的响应是不同的，当方差与信号波动的半径匹配时获得最大响应（图1-2）

图1-2. 方差为1的二阶高斯导数与不同半径的信号（左）卷积结果（右）

但此时我们还不能愉快的利用这个性质，因为未经处理的二阶高斯导数随方差增大会发生衰减，所以可能会出现这种情况（图1-3）

图1-3. 二阶高斯导数随方差增大的衰减现象

处理的方法叫规范化，很简单，就是在导数前乘上方差的平方。

公式推导：

(1)LOG

二阶高斯函数，规范化二阶高斯函数导数，

整理得，将导数对方差求导，取0值对应的极值，

这个就是信号半径与方差的关系。

(2)DOH

基本思路与LOG一样，只不过使用了Hessian矩阵：，对一张图计算xx，yy，xy三个方向的卷积，然后计算其加权行列式：。理论上，与LOG相比，DOH对细长结构的斑点有较好的抑制作用。

LOG斑点检测算法步骤：

（1）预定义一组方差值（因为不知道待检信号的尺度），对每个方差生成一个二阶高斯模板

注：这里取（DOH为，生成3个模板），有个放大差异的作用；模板尺寸取 

（2）对每个方差，将对应的高斯模板与原始信号做卷积（DOH需要将三个模板分别与原始图做卷积，然后计算其加权行列式），得到一组不同尺度的图像集

（3）对每个空间位置，比较其在图像集里26（3*3*3-1）个位置（图1-4）的值，如果为极值，则认为在该点有一个斑点

通过此时的尺度值可进一步得到该斑点的半径。

 

图1-4. 在26个位置里比较

2. 代码

#define MAXPOINTS 50000typedef struct myPoint{int x;int y;int radiu;}_point;void getLOGkernel(int halfWin, float sita, float** kernel);void cornerLOG(unsigned char* srcImg, _point* corner, int* count, int height, int width, float minScale, float maxScale, int stepScale);void converlution(unsigned char* srcImg, int** dstImg, int height, int width, float* kernel, int kernelWin);bool isMax(int* compImg, int compValue, _point compPoint, int kerWin, int imgWidth, bool isSameK);void cornerDOH(unsigned char* srcImg, _point* corner, int* count, int height, int width, float minScale, float maxScale, int stepScale);void getDOHkernel(int halfWin, float sita, float** kernelxx, float** kernelyy, float** kernelxy);void detValue(int** scaleimg, int* imgxx, int* imgyy, int* imgxy, float scale, int height, int width);void main(){_point corner[MAXPOINTS];int numCorner = 0;float minScale = 1;float maxScale = 20;float stepScale = 2;unsigned char* srcImg = NULL;cv::Mat img = cv::imread("../file/sunflowers.jpg", 0);int height = img.rows;int width = img.cols;srcImg = new unsigned char[height*width];for(int i=0; i<height; i++){uchar* ptr = img.ptr<uchar>(i);for(int j=0; j<width; j++)srcImg[i*width+j] = ptr[j];}#if 1cornerLOG(srcImg, corner, &numCorner, height, width, minScale, maxScale, stepScale);#elsecornerDOH(srcImg, corner, &numCorner, height, width, minScale, maxScale, stepScale);#endifdelete[] srcImg;cv::Mat showImg;cvtColor(img, showImg, CV_GRAY2BGR);for(int i=0; i<numCorner; i++){cv::Point cor(corner[i].x, corner[i].y);cv::circle(showImg, cor, corner[i].radiu, cv::Scalar(0,0,255));}cv::namedWindow("show");cv::imshow("show", showImg);cv::waitKey(0);}void getLOGkernel(int halfWin, float sita, float** kernel){int winSize = 2*halfWin+1;float tmp1, tmp2, sumValue = 0;float powsita = sita*sita;for(int i=-halfWin; i<=halfWin; i++){for(int j=-halfWin; j<=halfWin; j++){tmp1 = -1*(i*i+j*j)/(2*powsita);tmp2 = exp(tmp1)*(i*i+j*j-2*powsita);//exp(tmp1)*(1+tmp1)/(-1*powsita*powsita);sumValue += tmp2;(*kernel)[(i+halfWin)*winSize+(j+halfWin)] = tmp2;}}for(int i=0; i<winSize*winSize; i++)(*kernel)[i] -= sumValue/(winSize*winSize);}void converlution(unsigned char* srcImg, int** dstImg, int height, int width, float* kernel, int kernelWin){for(int i=0; i<height; i++){for(int j=0; j<width; j++){float sumValue = 0;int count = 0;for(int m=i-kernelWin; m<=i+kernelWin; m++){for(int n=j-kernelWin; n<=j+kernelWin; n++){if(m>=0 && m<height && n>=0 && n<width)sumValue += int(srcImg[m*width+n])*kernel[count];count++;}}sumValue *= 100;(*dstImg)[i*width+j] = int(sumValue);}}}bool isMax(int* compImg, int compValue, _point compPoint, int kerWin, int imgWidth, bool isSameK){for(int i=compPoint.y-kerWin; i<=compPoint.y+kerWin; i++){for(int j=compPoint.x-kerWin; j<=compPoint.x+kerWin; j++){if(isSameK && i==compPoint.y && j==compPoint.x)continue;if(abs(compValue) <= abs(compImg[i*imgWidth+j]))return false;}}return true;}void cornerLOG(unsigned char* srcImg, _point* corner, int* count, int height, int width, float minScale, float maxScale, int stepScale){int numScale = int((maxScale - minScale)/stepScale);int** scaledImg = new int*[numScale]; for(int i=0; i<numScale; i++)scaledImg[i] = new int[height*width];for(int k=0; k<numScale; k++){float scale = minScale+stepScale*k;int kernelWin = 3*scale;float *kernel = new float[(2*kernelWin+1)*(2*kernelWin+1)];getLOGkernel(kernelWin, scale, &kernel);converlution(srcImg, &(scaledImg[k]), height, width, kernel, kernelWin);delete[] kernel;}*count = 0;for(int i=1; i<height-1; i++){for(int j=1; j<width-1; j++){for(int k=1; k<numScale-1; k++){if((*count)>=MAXPOINTS){for(int m=0; m<numScale; m++)delete[] scaledImg[m];delete[] scaledImg;return;}_point cp;cp.x = j;cp.y = i;float scale = minScale+k*stepScale;cp.radiu = int(1.414*scale+0.5);if(isMax(scaledImg[k-1], scaledImg[k][i*width+j], cp, 1, width, false) &&isMax(scaledImg[k], scaledImg[k][i*width+j], cp, 1, width, true) &&isMax(scaledImg[k+1], scaledImg[k][i*width+j], cp, 1, width, false)){corner[(*count)++] = cp;break;}}}}for(int i=0; i<numScale; i++)delete[] scaledImg[i];delete[] scaledImg;}void getDOHkernel(int halfWin, float sita, float** kernelxx, float** kernelyy, float** kernelxy){int winSize = 2*halfWin+1;float tmp1, tmpxx, tmpyy, tmpxy;float sumValuexx = 0, sumValueyy = 0, sumValuexy = 0;float powsita = sita*sita;for(int i=-halfWin; i<=halfWin; i++){for(int j=-halfWin; j<=halfWin; j++){tmp1 = -1*(i*i+j*j)/(2*powsita);tmpxx = exp(tmp1)*(j*j-powsita);tmpyy = exp(tmp1)*(i*i-powsita);tmpxy = exp(tmp1)*i*j;sumValuexx += tmpxx;sumValueyy += tmpyy;sumValuexy += tmpxy;(*kernelxx)[(i+halfWin)*winSize+(j+halfWin)] = tmpxx;(*kernelyy)[(i+halfWin)*winSize+(j+halfWin)] = tmpyy;(*kernelxy)[(i+halfWin)*winSize+(j+halfWin)] = tmpxy;}}for(int i=0; i<winSize*winSize; i++){(*kernelxx)[i] -= sumValuexx/(winSize*winSize);(*kernelyy)[i] -= sumValueyy/(winSize*winSize);(*kernelxy)[i] -= sumValuexy/(winSize*winSize);}}void detValue(int** imgscale, int* imgxx, int* imgyy, int* imgxy, float scale, int height, int width){for(int i=0; i<height; i++){for(int j=0; j<width; j++){float xx = imgxx[i*width+j]/100.0;float yy = imgyy[i*width+j]/100.0;float xy = imgxy[i*width+j]/100.0;int vv = xx*yy-xy*xy;(*imgscale)[i*width+j] = int(vv/(scale*scale));}}}void cornerDOH(unsigned char* srcImg, _point* corner, int* count, int height, int width, float minScale, float maxScale, int stepScale){int numScale = int((maxScale - minScale)/stepScale);int** scaledImg = new int*[numScale]; for(int i=0; i<numScale; i++)scaledImg[i] = new int[height*width];for(int k=0; k<numScale; k++){float scale = minScale+stepScale*k;int kernelWin = 3*scale;float *kernelxx = new float[(2*kernelWin+1)*(2*kernelWin+1)];float *kernelyy = new float[(2*kernelWin+1)*(2*kernelWin+1)];float *kernelxy = new float[(2*kernelWin+1)*(2*kernelWin+1)];getDOHkernel(kernelWin, scale, &kernelxx, &kernelyy, &kernelxy);int** tmpImg = new int*[3]; for(int i=0; i<3; i++)tmpImg[i] = new int[height*width];converlution(srcImg, &(tmpImg[0]), height, width, kernelxx, kernelWin);converlution(srcImg, &(tmpImg[1]), height, width, kernelyy, kernelWin);converlution(srcImg, &(tmpImg[2]), height, width, kernelxy, kernelWin);detValue(&(scaledImg[k]), tmpImg[0], tmpImg[1], tmpImg[2], scale, height, width);for(int i=0; i<3; i++)delete[] tmpImg[i];delete[] tmpImg;delete[] kernelxx;delete[] kernelyy;delete[] kernelxy;}*count = 0;for(int i=1; i<height-1; i++){for(int j=1; j<width-1; j++){for(int k=1; k<numScale-1; k++){if((*count)>=MAXPOINTS){for(int m=0; m<numScale; m++)delete[] scaledImg[m];delete[] scaledImg;return;}_point cp;cp.x = j;cp.y = i;float scale = minScale+k*stepScale;cp.radiu = int(1.414*scale+0.5);if(isMax(scaledImg[k-1], scaledImg[k][i*width+j], cp, 1, width, false) &&isMax(scaledImg[k], scaledImg[k][i*width+j], cp, 1, width, true) &&isMax(scaledImg[k+1], scaledImg[k][i*width+j], cp, 1, width,false)){corner[(*count)++] = cp;break;}}}}for(int i=0; i<numScale; i++)delete[] scaledImg[i];delete[] scaledImg;}

效果图（左LOG，右DOH）