对于你给出的图片，我根据之前提到的思路编写了以下程序。感觉还蛮准的。这个程序只能用于两圆相交的情况。当小圆在大圆内部时，需要在其中加以判定。仅供参考。

clear;clc

c = imread('1.jpg');

cc = c;

tic

[m n] = size(c);

for i=1:m

index = find(~c(i,:));

if length(index) > 2

c(i,index(1)+1:index(end)-1)= true;

end

end

[x y] = find(~c);

maxdist = 0;

for i=1:length(x)

for j=i+1:length(x)

tmpdist = (x(i)-x(j))^2+(y(i)-y(j))^2;

if tmpdist > maxdist

maxdist = tmpdist;

pair = [i j];

end

end

end

x1 = x(pair(1));y1 = y(pair(1));

x2 = x(pair(2));y2 = y(pair(2));

dist = (y2-y1)*(x-x1)-(x2-x1)*(y-y1);

index = find(dist == max(dist));

x3 = x(index(1));y3=y(index(1));

index = find(dist == min(dist));

x4 = x(index(1));y4=y(index(1));

dist1 = (x3-x1)^2+(y3-y1)^2+(x4-x1)^2+(y4-y1)^2;

dist2 = (x3-x2)^2+(y3-y2)^2+(x4-x2)^2+(y4-y2)^2;

if dist2 < dist1

tmpx = x2; tmpy = y2;

x2 = x1; y2 = y1;

x1 = tmpx; y1 = tmpy;

end

centerx1 = ((y3-y1)*(x4^2+y4^2-x1^2-y1^2)-(y4-y1)*(x3^2+y3^2-x1^2-y1^2))/((y3-y1)*(x4-x1)-(y4-y1)*(x3-x1))/2;

centery1 = ((x3-x1)*(x4^2+y4^2-x1^2-y1^2)-(x4-x1)*(x3^2+y3^2-x1^2-y1^2))/((x3-x1)*(y4-y1)-(x4-x1)*(y3-y1))/2;

r1 = sqrt((x1-centerx1)^2+(y1-centery1)^2);

dist = (x-centerx1).^2+(y-centery1).^2;

index = find(dist < (r1+3)^2);

x(index)=[];

y(index)=[];

dist = (y2-y1)*(x-x1)-(x2-x1)*(y-y1);

index = find(dist == max(dist));

x5 = x(index(1));y5=y(index(1));

index = find(dist == min(dist));

x6 = x(index(1));y6=y(index(1));

centerx2 = ((y5-y2)*(x6^2+y6^2-x2^2-y2^2)-(y6-y2)*(x5^2+y5^2-x2^2-y2^2))/((y5-y2)*(x6-x2)-(y6-y2)*(x5-x2))/2;

centery2 = ((x5-x2)*(x6^2+y6^2-x2^2-y2^2)-(x6-x2)*(x5^2+y5^2-x2^2-y2^2))/((x5-x2)*(y6-y2)-(x6-x2)*(y5-y2))/2;

r2 = sqrt((x2-centerx2)^2+(y2-centery2)^2);

toc

imshow(cc);

hold on;

plot(centery1,centerx1,'r+');

plot(centery1+r1*cos(pi*(0:100)/50),centerx1+r1*sin(pi*(0:100)/50),'b-')

plot(centery2,centerx2,'r+');

plot(centery2+r2*cos(pi*(0:100)/50),centerx2+r2*sin(pi*(0:100)/50),'b-')