一次蛋疼的圆形图片圆心查找

先来看一下图

就是这样一张图，要算出两边圆形的圆心，方便校准镜头使用，刚开始试了很多方法，canny边缘查找、霍夫圆感觉都没办法找到上面的圆，更没办法用找到圆上的三点计算圆心了，后来由校准镜头第二步联想到可以让用户在图上选择三个点来进行计算，那就要用到 Matlab 了，不过肯定比直接用算法去找边好实现多了，来看 Matlab 的写法，里面用到了图片降噪、滤波、边缘查找、选点、画图，也算全面了，哈哈

function getCenter(inputRaw,  configJSON)

configStruct = loadjson(configJSON);

% 读取图片
rawCircle=imread('center.jpg');

i2=im2double(rawCircle);
% 灰度图
ihd=rgb2gray(i2);
% 降噪
[thr,sorh,keepapp]=ddencmp('den','wv',ihd);
ixc=wdencmp('gbl',ihd,'sym4',2,thr,sorh,keepapp);
%滤波
k2=medfilt2(ixc,[7 7]);
% 边缘查找
eprew=edge(k2,'canny');

if isfield(configStruct,'left1') && configStruct.left1 == 1
    outStruct = ...
        getRGGB(eprew,  'Left1');
    patchStruct.left.point1 = outStruct{1};
    patchStruct.point1 = outStruct{2};  
end
if isfield(configStruct,'left2') && configStruct.left2 == 1
    outStruct = ...
        getRGGB(eprew,  'Left2');
    patchStruct.left.point2 = outStruct{1};
    patchStruct.point2 = outStruct{2}; 
end
if isfield(configStruct,'left3') && configStruct.left3 == 1
    outStruct = ...
        getRGGB(eprew,  'Left3');
    patchStruct.left.point3 = outStruct{1};
    patchStruct.point3= outStruct{2};
end
if isfield(configStruct,'right1') && configStruct.right1 == 1
    outStruct = ...
        getRGGB(eprew,  'Right1');
    patchStruct.right.point1 = outStruct{1};
    patchStruct.point4 = outStruct{2};
end
if isfield(configStruct,'right2') && configStruct.right2 == 1
    outStruct = ...
        getRGGB(eprew,  'Right2');
    patchStruct.right.point2 = outStruct{1};
    patchStruct.point5 = outStruct{2};
end
if isfield(configStruct,'right3') && configStruct.right3 == 1
    outStruct = ...
            getRGGB(eprew,  'Right3');
    patchStruct.right.point3 = outStruct{1};
    patchStruct.point6 = outStruct{2};
end

draw(rawCircle,patchStruct);
saveFile(patchStruct);

end


function outStruct = getRGGB(varargin)

rawIn = varargin{1};
titleStr = varargin{2};

handle = figure('position',[0,0,2688,1344]);
set(gca,'position',[0,0,1,1]);
imagesc(rawIn)
title(titleStr)
[X, Y] = ginput(1);
X=fix(X);Y=fix(Y);
close(handle);

outStruct{1} =  sprintf('[%d,"%dj"]',X,Y);
outStruct{2} =  [X,Y];

end


% 写入json文件
function saveFile(struct)

fID = fopen('./centerPro.json','w');
fprintf(fID, '{\n');

fprintf(fID, '        "left1": %s,\n', struct.left.point1);
fprintf(fID, '        "left2": %s,\n', struct.left.point2);
fprintf(fID, '        "left3": %s,\n', struct.left.point3);

fprintf(fID, '        "right1": %s,\n', struct.right.point1);
fprintf(fID, '        "right2": %s,\n', struct.right.point2);
fprintf(fID, '        "right3": %s\n', struct.right.point3);

fprintf(fID, '}');
fclose(fID);

end

% 保存图片
function draw(rawCircle, struct)

x1 = [struct.point1(1) struct.point2(1) struct.point3(1) struct.point1(1)];
y1 = [struct.point1(2) struct.point2(2) struct.point3(2) struct.point1(2)];

x2 = [struct.point4(1) struct.point5(1) struct.point6(1) struct.point4(1)];
y2 = [struct.point4(2) struct.point5(2) struct.point6(2) struct.point4(2)];

handle = figure('position',[0,0,2688,1344]);
set(gca,'position',[0,0,1,1]);
imagesc(rawCircle)
hold on;
plot(x1,y1,'-.w');
plot(x2,y2,'-.w');
plot(x1,y1,'*r');
plot(x2,y2,'*r');
hold off;
F = getframe(gcf);
imwrite(F.cdata,'center1.jpg');
close(handle);

end

这样就能找到圆上可以用的三个点了，完美~
至于计算，我用的是 Python 写的一个计算小脚本，来看脚本

#!/usr/bin/python
# -*- coding: UTF-8 -*-
import sys,re
import json

def center(point1, point2, point3):
	x, y, z = point1[0]+complex(point1[1]), point2[0]+complex(point2[1]), point3[0]+complex(point3[1])
	w=z-x
	w/=y-x
	c=(x-y)*(w-abs(w)**2)/2j/w.imag-x
	return (abs(c.real),abs(c.imag))

if __name__ == "__main__":
	centerPro = json.load(open('centerPro.json'))
	left1 = centerPro['left1']
	left2 = centerPro['left2']
	left3 = centerPro['left3']
	right1 = centerPro['right1']
	right2 = centerPro['right2']
	right3 = centerPro['right3']
	
	left = center(left1, left2, left3)
	right = center(right1, right2, right3)

	str = 'left_x=%.3f\nleft_y=%.3f\nright_x=%.3f\nright_y=%.3f'%(left[0],left[1],right[0],right[1])
	print str

	with open('center.txt','w') as f:
		f.write(str)

这个比起上面的简单多了，花了一周就搞了个这。。。