2024年1月14日发(作者：)

- 由$${i_x}^2 times A = intint {{y^2}dxdy} $$$${i_y}^2 times A = intint {{x^2}dxdy} $$可以求得${i_x}^2$、${i_y}^2$####- 由$${a_x} = - {{{i_y}^2} over {{X_F}}}$$$${a_y} = - {{{i_x}^2} over {{Y_F}}}$$得到${X_F}$、${Y_F}$的值####- 迭代；作图----------# 代码## In "CoreFunction.m"####function CoreFunction(ImageLocation,M,N)%% Image initialImage = imread(ImageLocation);imshow(Image);title('OriginImage');hold off;Temp = input('This is the original any key to ');TwoBitImage = im2bw(Image);TwoBitImage = imresize(TwoBitImage,[400 400]);% black in 0;white in 1[X Y] = find(~TwoBitImage);DotInImage = [X Y];% Dots in black areaimshow(TwoBitImage);title('ImageInBlackAndWhite');Temp = input('This is the two-bit any key to ');EdgeImage = edge(TwoBitImage,'canny');% Edge in black area% black in 0;white in 1;

imshow(EdgeImage);title('EdgeImage');[X,Y]=find(EdgeImage);EdgeData =[X Y];Temp = input('This is the Edge any key to ');%% Edge data sortEdgeData = SortData(EdgeData(2:end,:),[],X(1),Y(1));%% CheckX = EdgeData(:,1);Y = EdgeData(:,2);plot(Y,X);%% FilterEdgeData = EdgeData(1:M:end,:);X = EdgeData(:,1);Y = EdgeData(:,2);plot(Y,X);hold on%% Confirm Poid(XingXin)% Coordinate of Xc = Sy / A;Sy = sum(DotInImage(:,2));A = length(DotInImage);Xc = (Sy / A);% Coordinate of Yc = Sx / A;Sx = sum(DotInImage(:,1));Yc = (Sx / A);% Scatterscatter(Yc,Xc)%% Confirm RadiusOfInertia^2(GuanXinBanJingDePingFang)% using Parallel Axis Fomula% Ix = Ixc + a^2*A%Ixc = sum(DotInImage(:,2).^2) - Yc^2*A;Ixc = sum((DotInImage(:,2)-Yc).^2);ix_2 = (Ixc/A);% Iy = Iyc + b^2*A

%Iyc = sum(DotInImage(:,1).^2) - Xc^2*A;Iyc = sum((DotInImage(:,1)-Xc).^2);iy_2 = (Iyc/A);%% Confirm ax and ay(QiuGeGeDianDeQieXianZaiXinXingZhuZhouDeJieJu)InterceptDataInXcYc = Intercept(Xc,Yc,X(1:end-2),Y(1:end-2),X(3:end) ...,Y(3:end));ax = InterceptDataInXcYc(:,1);ay = InterceptDataInXcYc(:,2);%% Calculate Xf and YfXf = -iy_2./ax;Yf = -ix_2./ay;%% Inverse Coordinate TransformationXf = (Xf + Xc);Yf = (Yf + Yc);R = SortData([Xf Yf],[],Xf(1),Yf(1));% N%Xf = R(1:N:end,1);Yf = R(1:N:end,2);% FilterDic = (Xf<0|Xf>400);Xf(Dic)=[];Yf(Dic)=[];Dic = (Yf<0|Yf>400);Xf(Dic)=[];Yf(Dic)=[];% CatinateXf = [Xf;Xf(1)];Yf = [Yf;Yf(1)];scatter(Yf,Xf,'.','k');set(gca,'ydir','reverse')end## In "Intercept.m"####

function Output = Intercept(Xc,Yc,X0,Y0,X1,Y1)%% Initial%% TranformX0 = X0 - Xc;Y0 = Y0 - Yc;X1 = X1 - Xc;Y1 = Y1 - Yc;%% CalculateK = (Y1-Y0)./(X1-X0);X = X0 - Y0./K;Y = -K.*X0 + Y0;%%Output = [X Y];Output = SortData(Output(2:end,:),[],Output(1,1),Output(1,2));end## In "SortData.m"####function R = SortData(OldData,NewData,X0,Y0)%% InitialX = OldData(:,1);Y = OldData(:,2);%% DistanceDistance = sqrt((X0-X).^2+(Y0-Y).^2);OldData = [X Y Distance];%% Sort and appendOldData = sortrows(OldData,3);X0 = OldData(1,1);Y0 = OldData(1,2);OldData = OldData(2:end,[1 2]);NewData = [NewData;[X0 Y0]];%% Recursionif(length(OldData) > 0)R = [SortData(OldData,NewData,X0,Y0);[X0 Y0]];else