sector2vol.c

超声波成像算法

#include "Reconst.h" 
#include <conio.h>

void Sector2Vol(double *vsector, double *vol,DAQINF daq, VIEWINF view){
	int i_theta,i_phai,i_sector,ix,iy,iz,i_vol,i;
	int nx,ny,nz,n_theta,n_phai,n_sector,i1,a1,a2,a3;
	double x,y,z,dx,dy,dz,drange,theta_max,dtheta;
	double delta1,delta2,delta3,distl;
	double v000,v001,v010,v011,v100,v101,v110,v111;
	double theta,phai,sector,max;
	distl = Distm/Lambda;
	nx = view.nx;
	ny = view.ny;
	nz = view.nz;
	n_theta = view.n_theta;
	n_phai = n_theta;
	n_sector = view.n_sector;
	//MaxFind(n_theta*n_phai*n_sector,vsector);
	drange = view.range/view.n_sector;
	theta_max = view.theta_max;
	dtheta = 2.0*theta_max/(double)n_theta;
	dx = view.vol_size/Lambda/(double)nx;
	dy = view.vol_size/Lambda/(double)ny;
	dz = view.vol_size/Lambda/(double)nz;  
	i_vol=0;
	a1 = n_sector*n_theta;
	a2= n_sector;
	a3=1;
	//vol = (double *)calloc(nx*ny*nz,sizeof(double));
	for(iz=0;iz<nz;iz++){
		z=(distl+(double)(iz-nz/2)*dz);	
		for(ix=0;ix<nx;ix++){
			x=(double)(ix-nx/2)*dx; 
			for(iy=0;iy<ny;iy++){
				y=(double)(iy-ny/2)*dy;
				theta = (atan2(x, z)+theta_max)/dtheta;
				phai = (atan2(y, z)+theta_max)/dtheta;
				sector = (sqrt(x*x+y*y+z*z)-view.range_min)/drange;
				i_theta = (int)theta;
				i_phai = (int)phai;
				i_sector = (int)sector;	
				delta1 = theta - (double)i_theta;
				delta2 = phai - (double)i_phai;
				delta3 = sector- (double)i_sector;
				if(i_sector>=0&&i_sector<n_sector&&i_theta>=0&&i_theta<n_theta
					&&i_phai>=0&&i_phai<n_phai){
					i1=a1*i_theta+a2*i_phai+i_sector;
					v000 = vsector[i1];
					v001 = vsector[i1+a3];
					v010 = vsector[i1+a2];
					v011 = vsector[i1+a2+a3];
					v100 = vsector[i1+a1];
					v101 = vsector[i1+a1+a3];
					v110 = vsector[i1+a1+a2];
					v111 = vsector[i1+a1+a2+a3];
					vol[i_vol] = (1.0-delta1)*(1.0-delta2)*(1.0-delta3)*v000 +
							(1.0-delta1)*(1.0-delta2)*delta3*v001 +
							(1.0-delta1)*delta2*(1.0-delta3)*v010 +
							(1.0-delta1)*delta2*delta3*v011 +
							delta1*(1.0-delta2)*(1.0-delta3)*v100 +
							delta1*(1.0-delta2)*delta3*v101 +
							delta1*delta2*(1.0-delta3)*v110 +
							delta1*delta2*delta3*v111;
					}//end if  
					i_vol++;		
			}//iy loop
		}//ix loop
	}//iz  loop
	/*
	i_vol=0;
	for(ix=0;ix<nx;ix++){
		for(iy=0;iy<ny;iy++){
			for(iz=0;iz<nz;iz++){
				vol[ix+nx*iy+nx*ny*iz]=vsector[i_vol++];
			}
		}
	}
	*/
	max=MaxFind(view.nx*view.ny*view.nz,vol);
	for(i=0;i<view.nx*view.ny*view.nz;i++){
		vol[i] = vol[i]/max;
	}
	max=MaxFind(view.nx*view.ny*view.nz,vol);
}