linereconst.c

超声波成像算法

/*                        LineReconst                                  */
/*      侾儔僀儞忋偺憸嵞惗傪FFT傪梡偄偰崅懍偵寁嶼                      */
#include <math.h>
#include "Reconst.h"
void LineReconst(AINF array,int n_time,int n_band,int i_delay,int i_shot,
	COMPLEX16 **wal_f,PVect p, int n_depth, double *depths, int ib,
	double *h, COMPLEX16 *rec, COMPLEX16 *s)
/*	堷悢丗
/*		AINF array	丗傾儗僀偵娭偡傞僨乕僞
/*		int n_time	丗帪娫僒儞僾儖悢
/*		int n_band	丗廃攇悢僒儞僾儖悢
/*		int i_delay	丗A/D曄姺偺抶傟帪娫乮僒儞僾儖悢乯
/*		COMPLEX16 wal_f[i][j]丗i斣栚偺憲怣攇宍偺戞j僼乕儕僄惉暘
/*		PVect p		丗埵抲儀僋僩儖
/*		int n_depth	丗墱峴偒偺暘妱悢乮嬤帡惛搙傛傝寛掕乯
/*		double depths[]丗墱峴暘妱椞堟偺嫬奅偲拞怱偺嫍棧
/*		int ib		丗
/*		double h[]	: 嶌嬈梡攝楍乮悺朄n_time*4乯
/*			僄僐乕攇宍偺僼乕儕僄惉暘偑寁嶼偝傟偰栠傞
/*		COMPLEX16 rec[j]丗i斣栚偺庴怣攇宍偺戞j僼乕儕僄惉暘偼
/*				rec[n_band*i+j]偱梌偊傜傟傞
/*		COMPLEX16 s[]丗暋慺悢偺嵞惗憸
/*		              s[ib+j]丗1儔僀儞偺j斣栚偺夋慺
/************************************************************/
{
	int		i, i_depth, i_freq,
			i_start, i_end, i_middle, i_bias,index, nsamp;
	double	*h1_real,*h1_imag,*h2_real,*h2_imag;
	double	Krx, Kry, Krz, drc, y0,z0, phase, buf1=0;
	double  L=0,Lmax=0,amp1=1,amp2=1,noise=0;
	int		tau_min, tau_max, tau_s_min, tau_s_max,n_bias;
	double tau, tau_bias, dtheta,cos_th,sin_th, low=10.0;
	double	R_min,R_max,R_start, R_end, R_middle;
	
	COMPLEX16 cx;
	dtheta    =	2*M_PI/(double)n_time;
	tau_bias=(double)(i_delay)/2.0;
	h1_real=h;
	h1_imag=h+  n_time;
	h2_real=h+2*n_time;
	h2_imag=h+3*n_time;
	nsamp=(int)(SPRATE*CYNUM*array.tnum);
	R_min=depths[0];/* 憸嵞惗椞堟偺嵟彫墱峴偒*/
	R_max=depths[2*n_depth];/* 憸嵞惗椞堟偺嵟戝墱峴偒*/
	tau_min=(int)(R_min*SPRATE*2.0);/* 憸嵞惗椞堟偱掕傑傞抶墑嵟彫抣*/
	tau_max=(int)(R_max*SPRATE*2.0);/* 憸嵞惗椞堟偱掕傑傞抶墑嵟戝抣*/
	tau_s_min=i_delay;             /*應掕僔僗僥儉偱掕傑傞抶墑嵟彫抣*/
	tau_s_max=i_delay+n_time-nsamp /*應掕僔僗僥儉偱掕傑傞抶墑嵟戝抣*/
	   -2*(int)(SPRATE*(array.radRc+array.radTr));									
	if(tau_min>tau_s_max || tau_max<tau_s_min) return;
	n_bias=4*n_time;
	//n_depth=8
		for(i_depth=0;i_depth<n_depth;i_depth++){
		R_start=depths[2*i_depth];
		R_middle=depths[2*i_depth+1];
		R_end=depths[2*i_depth+2];
		i_end=(int)(R_end*SPRATE*2.0)-tau_min;
		i_middle=(int)(R_middle*SPRATE*2.0)-tau_min;
		i_start=(int)(R_start*SPRATE*2.0)-tau_min;
		if((i_start<tau_s_max-tau_min)&&(i_end>tau_s_min-tau_min)){
			p.range=R_middle;
			p.x=R_middle*sin(p.theta);
			y0=R_middle*cos(p.theta);
			z0=R_middle*cos(p.theta);
			p.y=y0*sin(p.phai);
			p.z=z0*cos(p.phai);
			WaveFront0(array,n_time,n_band,i_delay,i_shot, wal_f, p,h);
			array.rxp=array.rx;
			array.ryp=array.ry;
            for(i=0;i<n_time;i++){
				h2_real[i]=0.0; 
				h2_imag[i]=0.0;
			}
			
			for ( i=0; i<array.rnum; i++ )
			//for ( i=0; i<0; i++ )
			{//庴怣巕偵娭偡傞儖乕僾
				i_bias=i*n_band;
				/* 儗僔乕僶偐傜僞乕僎僢僩傑偱偺曽岦儀僋僩儖傪寁嶼*/
				Krx = p.x - (*array.rxp);
				Kry = p.y - (*array.ryp);
				Krz = p.z;
				//儗僔乕僶偐傜僞乕僎僢僩傑偱偺嫍棧傪寁嶼
				drc = DIST( Krx, Kry, Krz);
				//僞乕僎僢僩乣儗僔乕僶傑偱
				//偺嫍棧傪僒儞僾儖悢偱寁嶼(僨傿儗僀傪寁嶼乯
				tau = dtheta*(drc*SPRATE - tau_bias);
				
				for(i_freq=0;i_freq<n_band;i_freq++){
				//  phase=dtheta*tau*(double)i_freq;
					phase=tau*(double)i_freq;
					phase=fmod(phase,2*M_PI);
					cx.real=rec[i_bias+i_freq].real;
					cx.imag=rec[i_bias+i_freq].imag;
					cos_th=cos(phase);
					sin_th=sin(phase);
				
					h2_real[i_freq] += (cx.real*cos_th-cx.imag*sin_th);
					h2_imag[i_freq] += (cx.real*sin_th+cx.imag*cos_th);
				}//i-freq loop
				array.rxp++;
				array.ryp++;
			}// i(receiver)loop
		
			 for(i_freq=0;i_freq<n_band;i_freq++){
				cx.real=h1_real[i_freq];
				cx.imag=h1_imag[i_freq];
				h1_real[i_freq]= h2_real[i_freq]*cx.real
								+h2_imag[i_freq]*cx.imag;
				h1_imag[i_freq]=-h2_real[i_freq]*cx.imag
								+h2_imag[i_freq]*cx.real;
			}
			for(i_freq=n_band;i_freq<n_time;i_freq++){
				h1_real[i_freq]= 0.0;
				h1_imag[i_freq]= 0.0;
			}
			
			FFT(n_time,h1_real,h1_imag,1.0);
			
			for(i=i_start;i<i_end;i++){
				index=(i-i_delay+n_bias)%n_time;
				s[ib+i].real+=h1_real[index];
				s[ib+i].imag+=h1_imag[index];
			}
		}//end if
	}//i_depth loop
	if(tau_min<tau_s_min){ 
		for(i=0;i<tau_s_min-tau_min;i++){
			s[ib+i].real=0.0;
			s[ib+i].imag=0.0;
		}
	}
	if(tau_max>tau_s_max){
		for(i=0;i<tau_max-tau_s_max;i++){
			s[ib+tau_max-tau_min-i].real=0.0;
			s[ib+tau_max-tau_min-i].imag=0.0;
		}
		
	}

	return;
}