примером обрасти может быть например круг, наложенная на нулевую матрицу с значениями матричных эллементов внутри матрицы равными единичкам.
выделяя контур окружности(в общем случае любой области) я получаю набор матричных координат(x,y), представляющих замкнутый контур, в ранном случае это окружность. Далее я апроксимирую его отрезками с некой погрешностью используя метод наименьших квадратов(он то и работает очень долго)
Я использую следующую функцию, возвращающе декартовы координаты ломанной(апроксимирующей), описывающей фигуру(последовательность точек с матричными координатами struct sElem *elem, колличество точек nodes, pogreshn- погрешность с которой строится апроксимирующая ломанная, *new_nodes- колличество точек в апроксимированной ломанной)
struct Point_dec * posledov_aprox_kont_s_pomowu_mnk(struct sElem *elem, int nodes,double pogreshn,int *new_nodes)
{
int i,xi,rez,j,k;
double b,c,d,sx,sy,aa,bb,cc,aa_old,bb_old,cc_old;
double sx2,sy2,sxy;
double *a_end,*b_end,*c_end;int*nom;
struct Point_dec *new_polig,tp,point;
a_end=(double*)malloc(sizeof(double)*nodes);
b_end=(double*)malloc(sizeof(double)*nodes);
c_end=(double*)malloc(sizeof(double)*nodes);
nom =(int*)malloc(sizeof(int)*nodes);
j=0+1;
b=0;c=0;d=0;sx=0;sy=0;
rez=1;
xi=0;
for(i=0;i<nodes;i++)
{
j=i+1;
if(j<nodes)
{
sx=(double)(elem[i].x);
sy=(double)(elem[i].y);
sx2=(double)(elem[i].x*elem[i].x);
sy2=(double)(elem[i].y*elem[i].y);
sxy=(double)(elem[i].x*elem[i].y);
aa=0;bb=0;cc=0;
do
{
aa_old=aa;
bb_old=bb;
cc_old=cc;
metod_naimenshih_kvadratov_posledov(elem, i,
&rez,
j,pogreshn, &sx, &sy,&sx2,&sy2,&sxy,
&aa,
&bb,
&cc);
j++;
}
// для каждой найденного апроксимирующего отрезка последовательность точек проверяем, чтобы расстояние от отрезка
//до точки было меньше погрешности
while (rez==1 && j < nodes );
if(rez==0)
{// с откатом
nom[xi]=j-2;
a_end[xi]=aa_old;
b_end[xi]=bb_old;
c_end[xi]=cc_old;
/////проверим на расстояние до прямой
//double del;
//del=aa_old* elem[i].x+ bb_old*+cc_old*
xi++;
i=j-3;
}
else if (rez>0 /*&& j == nodes-1*/)
{// без отката
nom[xi]=j-1;
a_end[xi]=aa;
b_end[xi]=bb;
c_end[xi]=cc;
i=nom[xi]-1;
xi++;
}
}
}
// для замыкающего отрезка
nom[xi]=0;
a_end[xi]=elem[nodes-1].y-elem[0].y;
b_end[xi]=-elem[nodes-1].x+elem[0].x;
c_end[xi]=elem[nodes-1].x*elem[0].y-elem[0].x*elem[nodes-1].y;
xi++;
//// построение результирующей ломанной
new_polig=(struct Point_dec *)malloc(sizeof(struct Point_dec )*2*xi);
k=0;
for(i=0;i<xi;i++)
{
if(i<xi-1)
findp_line(a_end[i],b_end[i],c_end[i],a_end[i+1],b_end[i+1],c_end[i+1],&tp,&rez);
/*находит точку пересечения tp прямых a_end[i]*x+b_end[i]*y+c_end[i]=0 и a_end[i+1]*x+b_end[i+1]*y+c_end[i+1]=0*/
else
findp_line(a_end[i],b_end[i],c_end[i],a_end[0],b_end[0],c_end[0],&tp,&rez);
if(rez==1 )
{
point.x=(double)(elem[nom[i]].x);
point.y=(double)(elem[nom[i]].y);
if(dlinna(tp,point) <=pogreshn)
{
new_polig[k]=tp;
k++;
}
else
{
perpendikular_point(a_end[i],b_end[i],c_end[i],point,&tp);
new_polig[k]=tp;
if(i<xi-1)
perpendikular_point(a_end[i+1],b_end[i+1],c_end[i+1],point,&tp);/*функция с помощью которой мы находим точку лежащую на пересечении 2 прямых -прямой a_end[i+1]*x+b_end[i+1]*x+c_end[i+1]=0
и прямой, проходящей через точку point и перпендикулярно a_end[i+1]*x+b_end[i+1]*x+c_end[i+1]=0 */
else
perpendikular_point(a_end[0],b_end[0],c_end[0],point,&tp);
new_polig[k+1]=tp;
k++;
k++;
}
}
}
free(a_end);
free(c_end);
free(b_end);
free(nom);
(*new_nodes)=k;
return new_polig;
}
void metod_naimenshih_kvadratov_posledov( struct sElem *elements,
int n_home,/*номер начальной точки с кот строим апроксимирующий отрезок*/
int *rez,
int j,double pogreshn,/*погрешность с которой проводим апроксимацию*/
double *sx,
double *sy,
double *sx2,
double *sy2,
double *sxy,
double *aa,
double *bb,
double *cc)
{
double dd;
int i,kol,nom;
double A,B, phi[2],p[2],d5,d6;
(*sx)=(*sx)+(double)(elements[j].x);
(*sy)=(*sy)+(double)(elements[j].y);
(*sx2)=(*sx2)+(double)(elements[j].x*elements[j].x);
(*sy2)=(*sy2)+(double)(elements[j].y*elements[j].y);
(*sxy)=(*sxy)+(double)(elements[j].x*elements[j].y);
kol=j-n_home+1;
A=2*((kol)*(*sxy)-(*sx)*(*sy));
B=-(*sx)*(*sx)+(kol)*(*sx2)+(*sy)*(*sy)-(kol)*(*sy2);
if(fabs(A)>eps2 )
{
phi[0]=atan((B-sqrt(A*A+B*B))/A);
phi[1]=atan((B+sqrt(A*A+B*B))/A);
}
if(fabs(A)<eps2 && fabs(B)<eps2)
{
phi[0]=pi;
phi[1]=pi;
}
if(fabs(A)<eps2 && fabs(B)>eps2)
{
phi[0]=pi/2.;
phi[1]=pi;
}
p[0]=(cos(phi[0])*(*sy)+sin(phi[0])*(*sx))/(kol);
p[1]=(cos(phi[1])*(*sy)+sin(phi[1])*(*sx))/(kol);
d5=kol*sqr(p[0])- 2.* p[0]* sin(phi[0]) *(*sx) + sqr(sin(phi[0]))* (*sx2) -
2.* p[0]* cos(phi[0])*(*sy) + 2.* cos(phi[0])*sin(phi[0])*(*sxy) +
sqr(cos(phi[0]))*(*sy2);
d6=kol*sqr(p[1])- 2.* p[1]* sin(phi[1]) *(*sx) + sqr(sin(phi[1]))* (*sx2) -
2.* p[1]* cos(phi[1])*(*sy) + 2.* cos(phi[1])*sin(phi[1])*(*sxy) +
sqr(cos(phi[1]))*(*sy2);
if( d5<d6 )
nom=0;
else
nom=1;
(*aa)=sin(phi[nom]);
(*bb)=cos(phi[nom]);
(*cc)=-p[nom];
////// теперь найдем макс индекс j при котором a*x+b*y+c<pogreshn
///// проверим отклонение прямой от контрольных точек
(*rez)=1;
// найдем расстояние до например нулевой точки
for(i=n_home;i<=j;i++)
{
dd=fabs((*aa)*((double)elements[i].x)+(*bb)*((double)elements[i].y)+(*cc));
if(dd>pogreshn)
{
(*rez)=0;
break;
// уменьшаем колличество рассматриваемых точек пополам
}
}
}
Здесь
struct Point_dec
{
double x;
double y;
};- для координат в декарте, тобишь с плавающей точкой, для метода наименьших квадратов
struct sElem
{
int x;
int y;
};- структура для матричных координат
