newreconstlut.c

超声波成像算法

/***************************************************************************/
/*	                	                                                   */
/*	    Walsh娭悢偱曄挷偝傟偨憲怣攇宍傪梡偄傞嶣憸僔僗僥儉偺                */
/*      3師尦憸嵞惗僾儘僌儔儉                            	               */
/*	    NewReconst.c				                                       */
/*      嬤幉嬤帡偵傛傝崅懍壔                                               */                         
/*	    			2001/ 9/4                                              */
/***************************************************************************/

#include "Reconst.h" 
#include <conio.h>
#include <stdlib.h>
void ReadParameter(AINF *array,DAQINF *daq);
void Sector2Vol(double *vsector, double *vol,DAQINF daq, VIEWINF view);

int main()
{
	PVect	p; //埵抲儀僋僩儖峔憿懱
	int		n_depth=3,n_range, n_range_max, n_sector=256;
	//n_depth:墱峴偒偺椞堟暘妱悢丆n_range丗墱峴偒暘妱悢丆n_sector丗墱峴偒偺夋慺悢
	double	r_min, r_max, depths[100];//嵟傕庤慜偺嫍棧丆嵟傕墦偄嫍棧丆墱峴偒偺嬫愗傝埵抲

	double	dth, *dely,buf1=0;
	double  *delta_xr,*delta_yr,*delta_zr,*delta_xt,*delta_yt,*delta_zt,
		    *detimet,*detimer,
		    *dist,*tang,*rang;
	double  L=0,Lmax=0,amp1=1,amp2=1,run=0,noise=0;
	int		shot, i_shot, nsht, method, nsamp,ib;
	int		i_theta,i_phai, n_theta=64;
	double	theta, theta_max, phai, d_theta, dist_l;
	int n_band; //n_band丗廃攇悢懷堟偺僒儞僾儖悢
	int		i_delay=0,ix_bias,iy_bias;
	COMPLEX16	*rec_f,**wal_f,//嶲徠梡攇宍偺僼乕儕僄曄姺丂悺朄n_rec*n_band
				*s;
	double		*vsector, *h;      //h丗嶌嬈梡攝楍丂悺朄4*NTIME
	int *ad_wave; 
	int j,Shots=1;

   	AINF	array;
	DAQINF	daq;
	VIEWINF view;
	char	fname1[STRLEN],*mname;
	FILE *fp;
	double CPU_time, start_time;
	double *vol;
	
	double pix_sum=0.0,max_pixel=0.0,avg_qave=0.0;
	int i;
	int n_time =NTIME;

	n_band=NTIME/4;
	view.n_sector=n_sector;
	Temp   = 20.0;/* 壏搙*/
	method = MOD2;/* 孞傝曉偟偺曽朄偼Mod2 */
	mname  = NAMMOD;
	nsht   = 1;
	i_shot=0;
//	onsoku  = 331.4+0.61*temp;		//嬻拞
	Onsoku  = 1513;					//悈拞
	

	/*A/D曄姺偝傟偨僄僐乕僨乕僞傪撉傒崬傓丂攝楍偺愭摢傊偺億僀儞僞
	ad_wave偵戙擖丂A/D奐巒傑偱偺抶墑僒儞僾儖悢i_delay傕僼傽僀儖偐傜撉傑傟傞*/
	Input_DATA_d("Input Image Size[mm] : ",&view.vol_size);
	view.vol_size *= 0.001;
	array.dl=0.75;	//傾儗僀偺慺巕悺朄
	ad_wave = fileread(&array,&daq);
	nsht = daq.nsht;
	i_delay = daq.ad_delay;
	Unitlength  = Onsoku/Spfrq;/* 侾僒儞僾儖帪娫偺娫偵壒偑恑傓嫍棧   */
	Lambda = Onsoku/Trfrq;
	dist_l = Distm/Lambda;// 拞怱廃攇悢偱偺攇挿偱惓婯壔

	view.distm = daq.ad_delay;
	view.n_theta=n_theta;

	dth     =	2*M_PI/SPRATE;/* 侾僒儞僾儖偁偨傝偺埵憡夞揮崁  */

	/****************** Walsh娭悢傪奿擺偡傞攝楍椞堟傪妋曐****************/

	Wal=calloc(array.tnum*array.tnum,sizeof(int));
	
    if (Wal==NULL){
	  puts("memory over!![wal]");
	  exit(0); 
	}

	n_range_max=NTIME;

	nsamp   = (int)(array.tnum*CYNUM*SPRATE);/* 憲怣攇宍宲懕帪娫乮僒儞僾儖*/
	detimet	= alloc_d(array.rnum); /*  僩儔儞僗儈僢僞偺僞僀儈儞僌僄儔乕   */
	detimer	= alloc_d(array.rnum); /*  儗僔乕僶偺僞僀儈儞僌僄儔乕         */
	rec_f	= alloc_complex(array.rnum*n_band*2);/////////////////////////////////////////
	//h		= alloc_d(nsht*n_theta*n_theta*NTIME*4);
	h   	= alloc_d(NTIME*4); /*     嶌嬈梡攝楍丂丂丂丂丂   */
	/*      儗僔乕僶偺庢傝晅偗埵抲岆嵎           */
	delta_xr= alloc_d(array.rnum);
	delta_yr= alloc_d(array.rnum);
	delta_zr= alloc_d(array.rnum);
	/*    僩儔儞僗儈僢僞偺庢傝晅偗埵抲岆嵎       */
	delta_xt= alloc_d(array.tnum);
	delta_yt= alloc_d(array.tnum);
	delta_zt= alloc_d(array.tnum);
	dist    = alloc_d(array.tnum*array.rnum);
	tang    = alloc_d(array.tnum);
	rang    = alloc_d(array.rnum);
	dely    = alloc_d(array.rnum*array.tnum);

	wal_f	= alloc_complex2d(array.tnum, NTIME);
	vsector	= alloc_d((n_theta+1)*(n_theta+1)*(n_sector+1));

/*----------A/D偝傟偨僄僐乕偺僼乕儕僄曄姺--------*/

	//FFTAdWave(array.rnum,NTIME, n_band, h, ad_wave, rec_f);

	putchar('\n');
	printf( "Number of Transmitter: %d\n",array.tnum);
	printf( "Number of Receiver   : %d\n",array.rnum);
	printf( "Distance             : %.2lf [m]\n",   Distm                );
	printf( "Temperature          : %.2lf [C]\n",   Temp                 );
	printf( "Carrier Frequency    : %.1lf [MHz]\n", Trfrq/1.0e6          );
	printf( "Sampling Frequency   : %.1lf [MHz]\n", Spfrq/1.0e6          );
	printf( "Wave Length          : %.4lf [mm]\n",  Unitlength*SPRATE*1000.0 );
	printf( "Delay for Sampling   : %d [word]\n",   i_delay                );
	printf( "Transmitting Method  : '%s'\n",        mname                );

	/* Walsh娭悢偑奿擺偝傟偨攝楍傪嶌傞 */
	walsh( Wal, array.tnum );
	walsh_ft(0, array.tnum,NTIME, n_band, wal_f);

	Input_DATA_d("Input View Angle [deg.]    : ",&theta_max);
	Input_DATA_d("Input r_min  [mm]    : ",&r_min);
	Input_DATA_d("Input r_max  [mm]    : ",&r_max);
	//r_min,r_max偼distm偱梌偊偨揰傪娷傓奿巕揰忋偵偁傞傛偆偵婯奿壔
	view.theta_max = theta_max*M_PI/180.0;
	r_min = Distm/Lambda - floor((Distm-r_min*1.0e-3)/Lambda);
	r_max = Distm/Lambda + floor((r_max*1.0e-3-Distm)/Lambda);
	n_range=(int)((r_max-r_min)*SPRATE*2.0);
	if(n_range>n_range_max) {
		n_range=n_range_max;
		r_min = Distm/Lambda - n_range/4.0/SPRATE;
		r_max = Distm/Lambda + n_range/4.0/SPRATE;
		printf("r_min & r_max are changed!\n");
	}
	printf("r_min=%lf   r_max=%lf\n",r_min,r_max);
	view.n_range = n_range;
	view.range = r_max-r_min;
	view.range_min = r_min;
	s		= alloc_complex(n_theta*n_theta*n_range);
	/*墱峴偒傪暘妱丅岆嵎乮攇挿/10乯偵婎偯偄偰暘妱偟丄暘妱悢n_depth
	偲椞堟偺嫬奅偲婎弨揰偺墱峴偒抣偺攝楍depths傪寁嶼
	depths[2*i+1]丗i斣栚偺椞堟偺婎弨揰丄
	depths[2*i]丄depths[2*i+2]丗i斣栚偺椞堟偺壓尷丄忋尷*/
	RangeDevision(array, r_min, r_max, &n_depth, depths);

	/*    嬻娫奿巕偺僗僥僢僾妏      */
	// 帇栰妏偼亇theta_max
	d_theta=2.0*theta_max*M_PI/180.0/(double)n_theta;    //n_theta; //僗僥僢僾妏[rad]
/********************************************************************************************/
/* 帇栰傪n_theta*n_theta屄偺曽岦偵暘妱 */
/* 曽岦偺悢偩偗儖乕僾偡傞 */
	shot=0;

//LUT僼傽僀儖僆乕僾儞
    Input_DATA_s("Input WaveFront LUT name : ",fname1);
    fp = fopen(fname1,"rb");
	if ( fp == NULL ) {
		fprintf( stderr, "File '%s' cannot open.\n", fname1 );
		exit(EXIT_FAILURE);
	}
    
	//LUT僼傽僀儖慡偰撉傒崬傒
    /*start_time=clock();
	for( i=0; i<nsht*n_theta*n_theta; i++){
	  fread( h, sizeof(double), n_band, fp );
      fread( h+n_time, sizeof(double), n_band, fp );
	}
    CPU_time=(double)(((clock()-start_time)/CLOCKS_PER_SEC));
	printf ("Read LUT Time = %g\n",CPU_time);*/
    
	//for(i_shot=0;i_shot<nsht;i_shot++)
    //FFTAdWave(array.rnum,NTIME, n_band, i_shot,h, ad_wave, rec_f);
	for(i_shot=0;i_shot<nsht;i_shot++){	
	printf("i_shot = %d\n",i_shot);
	FFTAdWave(array.rnum,NTIME, n_band, i_shot,h, ad_wave, rec_f);  
	    for ( i_theta=0; i_theta<n_theta; i_theta++ ) {
		//for ( i_theta=0; i_theta<64; i_theta++ ) {	
        start_time=clock();
		printf("i_theta,=%d\r",i_theta);
			theta=(i_theta - n_theta/2)*d_theta;
			p.theta=theta;
			ix_bias=i_theta*n_sector*n_theta;
			for( i_phai=0; i_phai<n_theta; i_phai++){ /********************************/
			//for( i_phai=0; i_phai<1; i_phai++){	
			    phai=(i_phai - n_theta/2)*d_theta;
				p.phai=phai;
				ib=i_theta*n_theta*n_range+i_phai*n_range;
				//LUT僼傽僀儖傪LineReconst偺倛偵撉傒崬傓
				fread( h, sizeof(double), n_band, fp );
	            fread( h+n_time, sizeof(double), n_band, fp );
					
				LineReconstLUT(array,NTIME,n_band,i_delay,i_shot,wal_f,
				               p,n_depth,depths,ib, h,rec_f, s);
				
				//h偺億僀儞僞堏摦
				//*(h+(n_band+n_time)*((i_phai+1)+(n_theta)*i_phai));

				iy_bias=i_phai*n_sector;
				j=0;
			} // i_phai loop
			CPU_time=(double)(((clock()-start_time)/CLOCKS_PER_SEC));
	        printf ("CPU Time = %g\n",CPU_time);
		}// i_theta loop
		//CPU_time=(double)(((clock()-start_time)/CLOCKS_PER_SEC));
	    //printf ("CPU Time = %g\n",CPU_time);
	}//i_shot loop
	for ( i_theta=0; i_theta<n_theta; i_theta++ ) {
		for( i_phai=0; i_phai<n_theta; i_phai++){
			ib=i_theta*n_theta*n_range+i_phai*n_range;
			LineDecimation(i_theta,i_phai,n_theta, n_theta,
							n_sector,n_range, ib, s, vsector);
		}
	}

//	CPU_time=(double)(((clock()-start_time)/CLOCKS_PER_SEC));
//	printf ("CPU Time = %g\n",CPU_time);

	Input_DATA_s("Input 2D Filename : ",fname1);
	MaxProjection(n_theta, n_theta, n_sector,"xz", vsector,fname1);

	view.nx=64;
	view.ny=64;
	view.nz=256;
	vol = (double *)calloc(view.nx*view.ny*view.nz,sizeof(double));
	
	Sector2Vol(vsector, vol,daq, view); 
	max_pixel=MaxFind(view.nx*view.ny*view.nz,vol);
	
	for(i=0;i<view.nx*view.ny*view.nz;i++)
		pix_sum+=vol[i];
	/****************************************/
	pix_sum/=((double)view.nx*view.ny*view.nz);
	printf("夋慺偺嵟戝抣=%lf\n",max_pixel);
	printf("夋慺偺暯嬒=%lf\n",pix_sum);
	avg_qave=10*(log10(pix_sum)-log10(max_pixel));
	printf("暯嬒婾憸儗儀儖=%lf [dB]\n",avg_qave);
	/*****************************************/
	Input_DATA_s("Input 3D Filename : ",fname1);
	strcat(fname1,".3d");
	fp = fopen( fname1, "wb" );
	if ( fp == NULL ) {
		fprintf( stderr, "File '%s' cannot open.\n", fname1 );
		exit(EXIT_FAILURE);
	}
	fwrite( vol, sizeof(double), view.nx*view.ny*view.nz, fp );	
	fclose(fp);
	return(0);
}