TA的每日心情 | 奋斗 2024-7-1 22:21 |
---|
签到天数: 2014 天 [LV.Master]伴坛终老
- 自我介绍
- 数学中国站长
群组: 数学建模培训课堂1 群组: 数学中国美赛辅助报名 群组: Matlab讨论组 群组: 2013认证赛A题讨论群组 群组: 2013认证赛C题讨论群组 |
2#
发表于 2005-6-9 23:04
|只看该作者
|
|邮箱已经成功绑定
例子
< >下面是一个定义数组的例子,我估计你定义不了是因为在栈里面无法定义,而在堆里面是可以定义的!</P>
9 ?( g( E% v- C$ p: g7 c* W< >#define T 100000</P>& |% _' Q h3 V4 L/ H& r
< >#include "stdlib.h"<BR>#include "math.h"<BR>#include "stdio.h"<BR>#include <string.h><BR>struct position<BR>{<BR> int x;<BR> int y;<BR>};<BR>//////////////////////////////////////////////////////////////////////////////<BR>//摘自《C常用算法程序集(第二版)》,实矩阵求逆的全选主高斯-约当法<BR>//<BR>//int brinv(a,n)<BR>//<BR>//本函数返回一个整形标志值。若返回的标志值为0,则表示矩阵A奇异,还输出信息<BR>//“err * * not inv”;若返回的标志值不为0,则表示正常返回。<BR>//形参说明:<BR>//a——双精度实型二维数组,体积为n*n。存放原矩阵A;返回时存放其逆矩阵A逆。<BR>//n——整形变量。矩阵的阶数。<BR>//////////////////////////////////////////////////////////////////////////////<BR>int brinv(double a[],int n)<BR>{ <BR> int *is,*js,i,j,k,l,u,v;<BR> double d,p;<BR> is=(int *)malloc((unsigned)n*sizeof(int));<BR> js=(int *)malloc((unsigned)n*sizeof(int));<BR> for (k=0; k<=n-1; k++)<BR> {<BR> d=0.0;<BR> for (i=k; i<=n-1; i++)<BR> for (j=k; j<=n-1; j++)<BR> {<BR> l=i*n+j; <BR> p=fabs(a[l]);<BR> if (p>d)<BR> {<BR> d=p;<BR> is[k]=i;<BR> js[k]=j;<BR> }<BR> }<BR> if (fabs(d)<0.0000001)<BR> {<BR> free(is);<BR> free(js);<BR> printf("err**not inv\n");<BR> return(0);<BR> }<BR> if (is[k]!=k)<BR> {<BR> for (j=0; j<=n-1; j++)<BR> { <BR> u=k*n+j;<BR> v=is[k]*n+j;<BR> p=a;<BR> a=a[v]; <BR> a[v]=p;<BR> }<BR> }<BR> if (js[k]!=k)<BR> {<BR> for (i=0; i<=n-1; i++)<BR> {<BR> u=i*n+k; <BR> v=i*n+js[k];<BR> p=a; <BR> a=a[v];<BR> a[v]=p;<BR> }<BR> }<BR> l=k*n+k;<BR> a[l]=1.0/a[l];<BR> for (j=0; j<=n-1; j++)<BR> {<BR> if (j!=k)<BR> {<BR> u=k*n+j;<BR> a=a*a[l];<BR> }<BR> }<BR> for (i=0; i<=n-1; i++)<BR> {<BR> if (i!=k)<BR> {<BR> for (j=0; j<=n-1; j++)<BR> {<BR> if (j!=k)<BR> {<BR> u=i*n+j;<BR> a=a-a[i*n+k]*a[k*n+j];<BR> }<BR> }<BR> }<BR> }<BR> for (i=0; i<=n-1; i++)<BR> {<BR> if (i!=k)<BR> {<BR> u=i*n+k;<BR> a=-a*a[l];<BR> }<BR> }<BR> }<BR> for (k=n-1; k>=0; k--)<BR> {<BR> if (js[k]!=k)<BR> {<BR> for (j=0; j<=n-1; j++)<BR> {<BR> u=k*n+j; <BR> v=js[k]*n+j;<BR> p=a; <BR> a=a[v];<BR> a[v]=p;<BR> }<BR> }<BR> if (is[k]!=k)<BR> {<BR> for (i=0; i<=n-1; i++)<BR> {<BR> u=i*n+k;<BR> v=i*n+is[k];<BR> p=a; <BR> a=a[v];<BR> a[v]=p;<BR> }<BR> }<BR> }<BR> free(is); <BR> free(js);<BR> return(1);<BR>}</P>, V$ [; }% |( l7 a3 O$ k6 c6 d
< >//////////////////////////////////////////////////////////////////////////////<BR>//摘自《C常用算法程序集(第二版)》,实矩阵相乘<BR>//<BR>//void brmul(a,b,m,n,k,c)<BR>//<BR>//形参说明:<BR>//a——双精度实型二维数组,体积为m*n。存放矩阵A的元素。<BR>//b——双精度实型二维数组,体积为n*k。存放矩阵B的元素。<BR>//m——整形变量。矩阵A的行数,也是乘积矩阵C=AB的行数。<BR>//n——整形变量。矩阵A的列数,也是矩阵B的行数。<BR>//k——整形变量。矩阵B的列数,也是乘积矩阵C=AB的列数。<BR>//c——双精度实型二维数组,体积为m*k。返回乘积矩阵C=AB的元素。<BR>//////////////////////////////////////////////////////////////////////////////<BR>void brmul(double a[],double b[],int m,int n,int k,double c[])<BR>{<BR> int i,j,l,u;<BR> for (i=0; i<=m-1; i++)<BR> {<BR> for (j=0; j<=k-1; j++)<BR> {<BR> u=i*k+j;<BR> c=0.0;<BR> for (l=0; l<=n-1; l++) c=c+a[i*n+l]*b[l*k+j];<BR> }<BR> }<BR> return;<BR>}</P>- u# g/ I9 k( K3 r( S. ?
< >//////////////////////////////////////////////////////////////////////////////<BR>void main()<BR>{<BR> FILE *para_fp,*in_fp,*out_fp,*fp;<BR> int n;// 特征空间的大小,即共n个特征<BR> int lines,columns;//对应数据场的行列值<BR> int m;//分类数<BR> int *nl;//一维数组,存放每类的样本点数量<BR> struct position **style;//二维数组,存放每类的样本点坐标<BR> double *xp,**x;//二维数组x,存放训练样本数据的数组;<BR> double *xtp,**xt;//数组x的逆,存放训练样本数据的数组;<BR> double *up,**u;//二维数组,存放每一类的均值<BR> double ***sigma,*sigmap,**sigmapp;//三维数组,存放每类的散度矩阵<BR> double ***anti_sigma,*anti_sigmap,**anti_sigmapp;//三维数组,存放每类的散度矩阵的逆矩阵<BR> double **sigma_temp,*sigma_tempp,**anti_sigma_temp,*anti_sigma_tempp;//临时用的矩阵<BR> int return_brinv;//矩阵求逆返回值,标示矩阵是否奇异<BR> double *g;//数组,存放像元x 与每一类的广义距离<BR> double min_g;<BR> int min_i;//存放最大g的i<BR> int datatype;//对应每一数据类型所占的字节数<BR> int i,j,k;//循环变量<BR> long count;//循环变量,扫描整个数据场<BR> char d[10];//原数据场的数据类型<BR> //以下为临时变量<BR> struct position style_temp;<BR> char pels_temp;<BR> char temp;<BR> short temp_0;<BR> int temp_1;<BR> long temp_2;<BR> float temp_3;<BR> double temp_4;</P>
: ]! L# B! S* L- V< > double g_temp0[1],*g_temp1p,**g_temp1;<BR>//////////////////////////////////////////////////////////////////////////////////////////////<BR>// 变量说明<BR>//<BR>// nl每类样本点数量 style每类的样本点坐标<BR>// | |<BR>// nl[0] style[0][0]---style[0][(nl[0])]<BR>// nl[1] style[1][0]---style[0][(nl[1])]<BR>//<BR>// …… ……<BR>//<BR>// nl[m] style[m][0]---style[0][(nl[m])] <BR>//<BR>//sigma[m][n][n]:存放所有类的sigma矩阵;<BR>//anti_sigma[m][n][n]:存放所有类的sigma矩阵的逆矩阵;<BR>//<BR>////////////////////////////////////////////////////////////////////////////////////////////////</P>
2 p! a R! b$ g9 |< ><BR>////////////读取特征空间的大小//////////////////////////////////<BR> if((para_fp=fopen("guangyi.par","rt"))==NULL)<BR> { <BR> return;<BR> }<BR> fscanf(para_fp,"%d",&n);//读取特征空间的大小<BR> fscanf(para_fp,"%d",&m);//预分类数<BR> //band=(int *)malloc((unsigned)n*sizeof(int));//存取n个波段号<BR> //for(i=0;i<n;i++) band=0;<BR> //for(i=0;i<n;i++)<BR> //{<BR> // fscanf(para_fp,"%d",&(band));<BR> //}<BR> fclose(para_fp);</P>% V& y7 f! w/ _& x
< > //for(i=0;i<n;i++)<BR> //{<BR> //}<BR>////////////读取数据场的行列值//////////////////////////////////<BR> if((para_fp=fopen("sourceimg.par","rt"))==NULL)<BR> { <BR> return;<BR> }<BR> fscanf(para_fp,"%d,%d",&lines,&columns);<BR> fclose(para_fp);</P>5 b) r. p: E/ R
< > printf("lines=%d,columns=%d\n",lines,columns);<BR> ///////////////////////////////////<BR> if((fp=fopen("leixing.par","rt"))==NULL)<BR> { <BR> return;<BR> }<BR> i=0;<BR> while(!feof(fp))<BR> { <BR> d=fgetc(fp);<BR> // printf("%c\n",d);<BR> i=i+1;<BR> }<BR> d[i-1]='\0';<BR> //fscanf(fp,"%s",d);</P>
$ I7 x; D1 t' b. k< > if(strcmp(d,"char")==0) <BR> {<BR> datatype=1;<BR> }<BR> if(strcmp(d,"short")==0)<BR> {<BR> datatype=2;<BR> }<BR> if(strcmp(d,"int")==0)<BR> {<BR> datatype=4;<BR> }<BR> if(strcmp(d,"long")==0)<BR> {<BR> datatype=4;<BR> }<BR> if(strcmp(d,"float")==0)<BR> {<BR> datatype=4;<BR> }<BR> if(strcmp(d,"double")==0)<BR> {<BR> datatype=8;<BR> }<BR> printf("datatype=%d\n",datatype);<BR> printf("%s",d);<BR>////////////读取样本数和样本区//////////////////////////////////<BR> style=(struct position **)malloc((unsigned)m*sizeof(struct position *));<BR> nl=(int*)malloc((unsigned)m*sizeof(int));</P>' i; o! a, y: R4 |" p8 q1 g
< > if((para_fp=fopen("style.par","rt"))==NULL)<BR> { <BR> free(style);<BR> free(nl);<BR> return;<BR> }<BR> for(i=0;i<m;i++)<BR> {<BR> fscanf(para_fp,"%d",&(nl));<BR> style=(struct position *)malloc((unsigned)(nl)*sizeof(struct position));<BR> for(j=0;j<(nl);j++) <BR> {<BR> style_temp.x=0;<BR> style_temp.y=0;<BR> style[j]=style_temp;<BR> }<BR> for(j=0;j<(nl);j++)<BR> {<BR> // 读入每一类的样本点的坐标<BR> fscanf(para_fp,"%d,%d",&(style_temp.x),&(style_temp.y));<BR> style[j]=style_temp;<BR> }<BR> }<BR> fclose(para_fp);</P>7 W0 j4 H. h# r E n
< > for(i=0;i<m;i++)<BR> {<BR> printf("\nClass %3d:%3d\n",i,nl);<BR> for(j=0;j<(nl);j++)<BR> {<BR> style_temp=style[j];<BR> // printf("(%d,%d)",style_temp.x,style_temp.y);<BR> }<BR> }<BR>////////////计算m个类别训练组的u,sigma,sigma逆//////////////////////////////////<BR> if((in_fp=fopen("sourceimg.dat","rb"))==NULL)<BR> { <BR> for(k=0;k<m;k++)<BR> {<BR> free(style[k]);<BR> }<BR> free(style);<BR> free(nl);<BR> return;<BR> }</P>" f5 h# n9 d9 T# Q
< > up=(double *)malloc((unsigned)n*m*sizeof(double));<BR> u=(double **)malloc((unsigned)n*sizeof(double*));<BR> for(k=0;k<n;k++)<BR> {<BR> u[k]=&(up[k*m]);<BR> }<BR> for(k=0;k<n;k++)<BR> {<BR> for(j=0;j<m;j++)<BR> {<BR> u[k][j]=0.0;<BR> }<BR> }</P>. g$ t+ t" n9 T: A! |" ]
< > sigmap=(double *)malloc((unsigned)m*n*n*sizeof(double));<BR> sigmapp=(double **)malloc((unsigned)m*n*sizeof(double *));<BR> sigma=(double ***)malloc((unsigned)m*sizeof(double **));</P>
/ V) N& @; q; r< > for(i=0;i<m;i++) <BR> {<BR> for(j=0;j<n;j++)<BR> {<BR> sigmapp[i*n+j]=&(sigmap[i*n*n+j*n]);<BR> }<BR> sigma=&(sigmapp[i*n]);<BR> }</P>1 b# R( O7 H. \' w9 h0 o
< > anti_sigmap=(double *)malloc((unsigned)m*n*n*sizeof(double));<BR> anti_sigmapp=(double **)malloc((unsigned)m*n*sizeof(double *));<BR> anti_sigma=(double ***)malloc((unsigned)m*sizeof(double **));</P>
8 J1 _1 i, c, ]' A2 g- X" m: I< > for(i=0;i<m;i++) <BR> {<BR> for(j=0;j<n;j++)<BR> {<BR> anti_sigmapp[i*n+j]=&(anti_sigmap[i*n*n+j*n]);<BR> }<BR> anti_sigma=&(anti_sigmapp[i*n]);<BR> }</P>8 d2 F/ g3 s) w
< > //对每一类进行循环,计算计算m个类别训练组的u,sigma,sigma逆<BR> for(i=0;i<m;i++)<BR> {<BR> sigma_tempp=(double *)malloc((unsigned)n*n*sizeof(double));<BR> sigma_temp=(double **)malloc((unsigned)n*sizeof(double*));<BR> for(k=0;k<n;k++)<BR> {<BR> sigma_temp[k]=&(sigma_tempp[k*n]);<BR> }</P>
+ U' B3 z5 N) k' O2 P& m3 ]< > anti_sigma_tempp=(double *)malloc((unsigned)n*n*sizeof(double));<BR> anti_sigma_temp=(double **)malloc((unsigned)n*sizeof(double*));<BR> for(k=0;k<n;k++)<BR> {<BR> anti_sigma_temp[k]=&(anti_sigma_tempp[k*n]);<BR> }</P>
2 C0 }9 ^" M1 }< > xp=(double *)malloc((unsigned)n*nl*sizeof(double));<BR> x=(double **)malloc((unsigned)n*sizeof(double*)); //训练样本<BR> for(k=0;k<n;k++)<BR> {<BR> x[k]=&(xp[k*nl]);<BR> }<BR> xtp=(double *)malloc((unsigned)nl*n*sizeof(double));<BR> xt=(double **)malloc((unsigned)nl*sizeof(double*)); //训练样本的转置<BR> for(j=0;j<nl;j++)<BR> {<BR> xt[j]=&(xtp[j*n]);<BR> }</P>; L0 u& t2 ^3 z3 y# b+ B
< ><BR>///////////////////读取矩阵X///////////////////////////////////<BR> for(j=0;j<nl;j++)<BR> { <BR> style_temp=style[j];<BR> for(k=0;k<n;k++)<BR> {<BR> //style_temp=style[j];<BR> if(k==0) <BR> { <BR> fseek(in_fp,(style_temp.x*lines+style_temp.y)*datatype,SEEK_SET);<BR> }<BR> else<BR> {<BR> fseek(in_fp,(lines*columns-1)*datatype,SEEK_CUR);<BR> }<BR> /////////////////从数据场文件中读入训练样本的数据<BR> if(1==datatype)<BR> {<BR> fread(&temp,datatype,1,in_fp);<BR> x[k][j]=(double)temp;<BR> }<BR> if(2==datatype)<BR> {<BR> fread(&temp_0,datatype,1,in_fp);<BR> x[k][j]=(double)temp_0;<BR> }<BR> if(strcmp(d,"int")==0)<BR> {<BR> fread(&temp_1,datatype,1,in_fp);<BR> x[k][j]=(double)temp_1;<BR> }<BR> if(strcmp(d,"long")==0)<BR> {<BR> fread(&temp_2,datatype,1,in_fp); <BR> x[k][j]=(double)temp_2;<BR> } <BR> if(strcmp(d,"float")==0)<BR> {<BR> fread(&temp_3,datatype,1,in_fp); <BR> x[k][j]=(double)temp_3;<BR> }<BR> if(8==datatype)<BR> {<BR> fread(&temp_4,datatype,1,in_fp); <BR> x[k][j]=(double)temp_4;<BR> }<BR> /////////////////////<BR> // fread(&temp,datatype,1,in_fp);<BR> // x[k][j]=(double)temp;<BR> // if(x[k][j]<0) x[k][j]+=256;<BR> }<BR> }<BR> // printf("\nWhen Class= %3d,X is:\n",i);<BR> for(k=0;k<n;k++)<BR> {<BR> for(j=0;j<nl;j++)<BR> {<BR> // printf("%10.5f",x[k][j]);<BR> }<BR> // printf("\n");<BR> }<BR>///////////////////求X的转置矩阵Xt////////////////////////////////<BR> for(k=0;k<n;k++)<BR> {<BR> for(j=0;j<nl;j++)<BR> {<BR> //计算训练样本的转置矩阵<BR> xt[j][k]=x[k][j];<BR> }<BR> }</P>
6 m; Z2 A x5 a- W< > printf("\nWhen Class= %3d,Xt is:\n",i);<BR> for(j=0;j<nl;j++)<BR> {<BR> for(k=0;k<n;k++)<BR> {<BR> // printf("%10.5f",xt[j][k]);<BR> }<BR> // printf("\n");<BR> }</P>
2 t9 Z+ l# z' e/ y, e< >///////////////////求矩阵u/////////////////////////////////////////<BR> printf("\nWhen Class= %3d,U is:\n",i);<BR> for(k=0;k<n;k++)<BR> {<BR> for(j=0;j<nl;j++)<BR> { <BR> u[k]+=x[k][j];<BR> }<BR> u[k]=(u[k])/(nl);//求每一类的样本均值<BR> // printf("%13.7lf\n",u[k]);<BR> }</P>0 y4 d, A5 D7 L0 x5 K) `& W& @6 [
< > // printf("\n");</P>3 Q2 b% g2 @, J$ w! n
< >///////////////////计算矩阵Sigma=X*Xt//////////////////////////////<BR> brmul(xp,xtp,n,nl,n,sigma_tempp);<BR> for (k=0;k<n;k++)<BR> for(j=0;j<n;j++)<BR> sigma_temp[k][j]=sigma_temp[k][j]/(nl-1); //求每一类的散度矩阵<BR> printf("\nWhen Class= %3d,Sigma is:\n",i);<BR> for(k=0;k<n;k++)<BR> {<BR> for(j=0;j<n;j++)<BR> {<BR> // printf("%13.3lf",sigma_temp[k][j]);<BR> }<BR> // printf("\n");<BR> }</P>
& x6 W5 v: F% B0 c7 t< ><BR>///////////////////计算矩阵Sigma的逆矩阵anti_sigma////////////////<BR> for(k=0;k<n;k++)<BR> {<BR> for(j=0;j<n;j++)<BR> {<BR> anti_sigma_temp[k][j]=sigma_temp[k][j];<BR> }<BR> }</P>2 x1 N- Q0 Y& M' N
< > return_brinv=brinv(anti_sigma_tempp,n);//求散度矩阵的逆矩阵</P>/ L$ z& H/ g% s P& [9 ~
< > if(return_brinv==0)//矩阵Sigma奇异,退出程序,重新选择参数<BR> {<BR> free(x);<BR> free(xp);<BR> free(xt);<BR> free(xtp);<BR> free(sigma_tempp);<BR> free(sigma_temp);<BR> free(anti_sigma_tempp);<BR> free(anti_sigma_temp);</P>
4 `/ m4 v- r# I6 u0 L) D, ^< > <BR> <BR> // free(p);<BR> free(nl);<BR> for(k=0;k<m;k++)<BR> {<BR> free(style[k]);<BR> }<BR> free(style);<BR> free(u);<BR> free(up);<BR> free(sigmap);<BR> free(sigmapp);<BR> free(sigma);<BR> free(anti_sigmap);<BR> free(anti_sigmapp);<BR> free(anti_sigma);<BR> return;<BR> }</P>
+ h- F6 e/ i( ~! F! Q( V< > printf("\nWhen Class= %3d,Anti_Sigma is:\n",i);<BR> for(k=0;k<n;k++)<BR> {<BR> for(j=0;j<n;j++)<BR> {<BR> // printf("%13.3e ",anti_sigma_temp[k][j]);<BR> }<BR> // printf("\n");<BR> }</P>3 s5 Z7 C2 P: m# K! }( e
<P><BR>// printf("\nWhen Class= %3d,Det Sigma is:%13.7lf\n",i,bsdet(sigma_tempp,n));</P>/ P4 ^, ~! ~" i7 Q, J
<P>/*<BR> double c[n][n];<BR> brmul(sigma_tempp,anti_sigma_tempp,n,n,n,c);<BR> printf("\n\n");<BR> for(k=0;k<n;k++)<BR> {<BR> for(j=0;j<n;j++)<BR> {<BR> printf("%13.3e ",c[k][j]);<BR> }<BR> printf("\n");<BR> }<BR>*/</P>
' W: c( s, d" ~" p) H3 E<P><BR>//////////////////////////将每类的Sigma,Anti_Sigma存放入三维数组//////<BR> for(k=0;k<n;k++)<BR> {<BR> for(j=0;j<n;j++)<BR> {<BR> sigma[k][j]=sigma_temp[k][j];<BR> anti_sigma[k][j]=anti_sigma_temp[k][j];<BR> }<BR> }<BR> free(x);<BR> free(xp);<BR> free(xt);<BR> free(xtp);<BR> free(sigma_tempp);<BR> free(sigma_temp);<BR> free(anti_sigma_tempp);<BR> free(anti_sigma_temp);<BR> }</P>6 t# o6 c$ `1 z
<P>/////////////////////////////////////////////////////////////////////////////<BR>//////////逐点扫描整个数据场<BR> if((out_fp=fopen("m_printf.dat","wb"))==NULL)<BR> { <BR> free(sigmap);<BR> free(sigmapp);<BR> free(sigma);<BR> free(anti_sigmap);<BR> free(anti_sigmapp);<BR> free(anti_sigma);<BR> return;<BR> }</P>2 ^# C# {. G8 b" V* c' {& h
<P> xp=(double *)malloc((unsigned)n*1*sizeof(double));<BR> x=(double **)malloc((unsigned)n*sizeof(double*));<BR> for(k=0;k<n;k++)<BR> {<BR> x[k]=&(xp[k*1]); //这里的数组x是逐个要扫描的每一个数据点<BR> }<BR> xtp=(double *)malloc((unsigned)1*n*sizeof(double));<BR> xt=(double **)malloc((unsigned)1*sizeof(double*));//这里的数组xt是逐个要扫描的每一个数据点的转置向量<BR> for(j=0;j<1;j++)<BR> {<BR> xt[j]=&(xtp[j*n]);<BR> }<BR> g=(double*)malloc((unsigned)m*sizeof(double));<BR> sigma_tempp=(double *)malloc((unsigned)n*n*sizeof(double));<BR> sigma_temp=(double **)malloc((unsigned)n*sizeof(double*));<BR> for(k=0;k<n;k++)<BR> {<BR> sigma_temp[k]=&(sigma_tempp[k*n]);<BR> }</P>
1 ^' y0 ?; r( ~7 x% a3 _<P> anti_sigma_tempp=(double *)malloc((unsigned)n*n*sizeof(double));<BR> anti_sigma_temp=(double **)malloc((unsigned)n*sizeof(double*));<BR> for(k=0;k<n;k++)<BR> {<BR> anti_sigma_temp[k]=&(anti_sigma_tempp[k*n]);<BR> }<BR> g_temp1p=(double *)malloc((unsigned)1*n*sizeof(double));<BR> g_temp1=(double **)malloc((unsigned)1*sizeof(double*));<BR> for(j=0;j<1;j++)<BR> {<BR> g_temp1[j]=&(g_temp1p[j*n]);<BR> }<BR>//读取数据场中的数据</P>. a3 t! j/ \4 t1 |% x4 z6 Q
<P> for(count=0;count<lines*columns;count++)<BR> {<BR> //读取一个数据点值<BR> for(k=0;k<n;k++)<BR> {<BR> if(k==0)<BR> {<BR> fseek(in_fp,count*datatype,SEEK_SET);<BR> }<BR> else<BR> {<BR> fseek(in_fp,(lines*columns-1)*datatype,SEEK_CUR);<BR> }<BR> ////////////////////根据原数据场的数据类型,读取每个数据点<BR> if(1==datatype)<BR> {<BR> fread(&temp,datatype,1,in_fp);<BR> x[k][0]=(double)temp;<BR> }<BR> if(2==datatype)<BR> {<BR> fread(&temp_0,datatype,1,in_fp);<BR> x[k][0]=(double)temp_0;<BR> }<BR> if(strcmp(d,"int")==0)<BR> {<BR> fread(&temp_1,datatype,1,in_fp);<BR> x[k][0]=(double)temp_1;<BR> }<BR> if(strcmp(d,"long")==0)<BR> {<BR> fread(&temp_2,datatype,1,in_fp); <BR> x[k][0]=(double)temp_2;<BR> }<BR> if(strcmp(d,"float")==0)<BR> {<BR> fread(&temp_3,datatype,1,in_fp); <BR> x[k][0]=(double)temp_3;<BR> }<BR> if(8==datatype)<BR> {<BR> fread(&temp_4,datatype,1,in_fp); <BR> x[k][0]=(double)temp_4;<BR> }<BR> /////////////////////<BR> //fread(&temp,datatype,1,in_fp);<BR> //x[k][0]=(double)temp;<BR> //if(x[k][0]<0) x[k][0]+=256;<BR> }</P>% f3 S& y' K' y3 f6 j
<P><BR> for(i=0;i<m;i++)<BR> {<BR>///////////////计算矩阵(x-u)以及其转置///////////////////////////////////<BR> for(k=0;k<n;k++)<BR> {<BR> x[k][0]-=u[k];<BR> xt[0][k]=x[k][0];<BR> }<BR>///////////////取出本类的矩阵Sigma及其逆矩阵Anti_Sigma///////////////////<BR> for(k=0;k<n;k++)<BR> {<BR> for(j=0;j<n;j++)<BR> {<BR> sigma_temp[k][j]=sigma[k][j];<BR> anti_sigma_temp[k][j]=anti_sigma[k][j];<BR> }<BR> }<BR>///////////////计算g,并求出最小的g//////////////////////////////////<BR> brmul(xtp,anti_sigma_tempp,1,n,n,g_temp1p);<BR> brmul(g_temp1p,xp,1,n,1,g_temp0);<BR> //计算数据点与每一类的距离<BR> g=g_temp0[0];<BR> if(i==0) <BR> {<BR> min_g=g;<BR> min_i=0;<BR> }<BR> else<BR> {<BR> if(min_g>g)<BR> {<BR> min_g=g;<BR> min_i=i;<BR> }<BR> }<BR> }<BR>///////////////将该数据点所属类别写入分类文件/////////////////////////////</P>* U8 [9 l6 w0 @9 _+ }& a! B. A
<P> if(min_g>T) min_i=m;////////////////如果min_g>T,划为拒绝类///////<BR> pels_temp=(char)(min_i*255/m);<BR> fwrite(&pels_temp,sizeof(char),1,out_fp);//保存最终的分类结果</P>9 v! H; V K: Q4 A5 P' D
<P> }<BR> fclose(out_fp);<BR> fclose(in_fp);</P>$ b3 o9 i' A' C2 ^
<P><BR> free(x);<BR> free(xp);<BR> free(xt);<BR> free(xtp);<BR> free(g);<BR> free(g_temp1p);<BR> free(g_temp1);<BR> free(sigma_tempp);<BR> free(sigma_temp);<BR> free(anti_sigma_tempp);<BR> free(anti_sigma_temp);</P>0 f* r9 W" r8 W% r. f, S
. ]5 g5 A$ y: f' F m6 o
! K2 v* n/ ]. h. I3 O3 ~+ m* J<P><BR>// free(p);<BR> free(nl);<BR> for(i=0;i<m;i++)<BR> {<BR> free(style);<BR> }<BR> free(style);<BR> free(u);<BR> free(up);<BR> free(sigmap);<BR> free(sigmapp);<BR> free(sigma);<BR> free(anti_sigmap);<BR> free(anti_sigmapp);<BR> free(anti_sigma);<BR>}<BR></P> |
|