d_radon.f90

线性radon变换

!  d_radon.f90 
!
!  FUNCTIONS:
!	d_radon      - Entry point of console application.
!
!	Example of displaying 'Hello World' at execution time.
!

!****************************************************************************
!
!  PROGRAM: d_radon
!
!  PURPOSE:  Entry point for 'Hello World' sample console application.
!
!****************************************************************************
	
	program d_radon
	implicit none

	!定义各变量
	integer :: m,n,k,h
	integer :: pk,th,xm,yn
	real,allocatable :: p(:),t(:),g(:,:),q(:,:),g_radon(:,:)
	character(len=80) :: filename1="g_data.txt"
	character(len=80) :: filename2="h_data.txt"
	character(len=80) :: filename3="data.txt"
	real,parameter :: field1=10
	real,parameter :: field2=20
	real,parameter :: field3=30
	real,parameter :: delta_x=1
	real,parameter :: delta_y=1
	real,parameter :: delta_p=0.02
	real,parameter :: delta_t=1
	real,parameter :: x_min=-50
	real,parameter :: y_min=0
	real,parameter :: p_min=-1
	real,parameter :: t_min=0
	real :: alpha,beta,round
	real :: sum


	!输入各采样数并申请空间
	write(*,*)"输入像数值！"
	read(*,*) xm,yn,pk,th
	allocate(g_radon(pk,th))
	allocate(g(pk,th))
	allocate(q(pk,th))
	allocate(p(pk))
	allocate(t(th))

	
	
	!定义原函数g(m,n)
	open(field1,file=filename1)
	do m=0,xm-1
		do n=0,yn-1
				g(m+1,n+1)=0
		enddo
	enddo		

	do m=0,xm-1
		round=0.4*(x_min+m*delta_x)+20
			if (int(round+0.5)>int(round)) then
				n=int(round)+1
			else if (int(round+0.5)<=int(round)) then
				n=int(round)
			endif
				
			if (n>=0 .and. n<yn) then
				g(m+1,n+1)=1
			endif	
		write(field1,"(E15.7)") g(m+1,n+1)			
	end do
	close(field1)				 
		

	!原函数在直线上的值为1，不在线上的为0
	open(field2,file=filename2)
	do m=0,xm-1
		do n=0,yn-1
			if (g(m+1,n+1)/=0) then
				q(m+1,n+1)=1
			else
				q(m+1,n+1)=0
			endif	
				write(field2,"(E15.7)") q(m+1,n+1)			
		end do
	end do
	close(field2)


	!radon变换
	open(field3,file=filename3)
	do k=0,pk-1                                      !all values of p
		do h=0,th-1                                    !all values of tau
			p(k+1)=p_min+k*delta_p                      !the digital slope at k
			alpha=p(k+1)*delta_x/delta_y                !calculate digital slope
			t(h+1)=t_min+h*delta_t                      !the digital tau at h
			beta=(p(k+1)*x_min+t(h+1)-y_min)/delta_y    !calculate digital offset
			sum=0.0                         
			do m=0,xm-1
				round=alpha*m+beta
					if (int(round+0.5)>int(round)) then
						n=int(round)+1
					else if (int(round+0.5)<=int(round)) then
						n=int(round)
					endif
				
					if (n>=0 .and. n<yn) then
						sum=sum+g(m+1,n+1)
					endif	
			end do
				g_radon(k+1,h+1)=delta_x*sum
				write(field3,"(E15.7)") g_radon(k+1,h+1)
		end do
	end do	
	close(field3)
	     
	end program d_radon