428 lines
14 KiB
Fortran
428 lines
14 KiB
Fortran
subroutine rect_with_looks(infile,outfile,ndac,nddn,nrac,nrdn,a,b,c,d,e,f,lac,ldn,lac0,ldn0,filetype,intstyle)
|
|
c
|
|
c*****************************************************************************
|
|
c**
|
|
c** FILE NAME: rect.f
|
|
c**
|
|
c** DATE WRITTEN: 27-Nov-98
|
|
c**
|
|
c** PROGRAMMER: P.A.Rosen
|
|
c**
|
|
c** FUNCTIONAL DESCRIPTION: This program adjusts an image
|
|
c** by affine transformation and interpolation
|
|
c**
|
|
c** UPDATE LOG: Cunren Liang, 24-APR-2015
|
|
c** 1. support file types: real, double
|
|
c** for double, after reading in, data are saved in float variables
|
|
c** for caculation, and data written out get values from float variables.
|
|
c** 2. support the case that input and output file sizes are different
|
|
c** 3. support affine transformation coefficients estimated from multilook
|
|
c** images. In this case, the program assumes that, the locations of the
|
|
c** offset estimations start from (0, 0). The images handled by this program
|
|
c** also start from (0, 0).
|
|
C**
|
|
c** Cunren LIANG, 03-JUN-2015
|
|
c** 1. there is 1 pixel offset in the output, which is corrected.
|
|
c**
|
|
c** Cunren Liang, 20-APR-2017
|
|
C** 1. support original image of multiple looks
|
|
c**
|
|
c**
|
|
c**
|
|
c*****************************************************************************
|
|
|
|
|
|
implicit none
|
|
|
|
c integer CMAX, RMAX
|
|
c parameter (CMAX = 7000, RMAX = 7200)
|
|
c real*4 rvs(0:2*CMAX-1,0:RMAX-1)
|
|
c complex*8 carr(0:20000)
|
|
c real*4 arr(0:40000)
|
|
|
|
c input of resampling
|
|
REAL*4, DIMENSION(:,:), ALLOCATABLE :: rvs
|
|
|
|
c variables for reading data
|
|
c For byte format *****GP 01-05******
|
|
INTEGER*1, DIMENSION(:), ALLOCATABLE :: barr
|
|
real*4, DIMENSION(:), ALLOCATABLE :: rarr
|
|
real*8, DIMENSION(:), ALLOCATABLE :: darr
|
|
COMPLEX*8, DIMENSION(:), ALLOCATABLE :: carr
|
|
|
|
c output of resampling
|
|
REAL*4, DIMENSION(:), ALLOCATABLE :: arr
|
|
|
|
real*4 pt1(3),pt2(3),pt3(3),pt4(3)
|
|
real*8 colval, rowval, ocolval, orowval, ocolval_tmp, orowval_tmp
|
|
real*8 ifrac, jfrac
|
|
real*8 a,b,c,d,e,f
|
|
real*4 interp
|
|
|
|
integer oi, oj, i, j, k, ift, iis
|
|
integer iargc, ndac, nddn, nrac, nrdn, ierr
|
|
integer nac
|
|
integer psize
|
|
integer lac, ldn
|
|
integer lac0, ldn0
|
|
|
|
character*180 fname, infile, outfile, intstyle, filetype
|
|
|
|
integer rdflen
|
|
character*255 rdfval
|
|
character*255 rdfcullsp,rdfdata
|
|
character*255 a_rdtmp
|
|
|
|
save rvs
|
|
|
|
c if(iargc() .eq. 0) then
|
|
c write(*,*) 'usage: rect rect.rdf'
|
|
c stop
|
|
c end if
|
|
|
|
c call getarg(1,fname)
|
|
|
|
c call rdf_init('ERRFILE=SCREEN')
|
|
c write(6,'(a)') 'Reading command file data...'
|
|
c call rdf_read(fname)
|
|
|
|
c a_rdtmp = rdfval('Input Image File Name','-')
|
|
c read(unit=a_rdtmp,fmt='(a)') infile
|
|
c a_rdtmp = rdfval('Output Image File Name','-')
|
|
c read(unit=a_rdtmp,fmt='(a)') outfile
|
|
c a_rdtmp = rdfval('Input Dimensions','-')
|
|
c read(unit=a_rdtmp,fmt=*) ndac, nddn
|
|
c a_rdtmp = rdfval('Output Dimensions','-')
|
|
c read(unit=a_rdtmp,fmt=*) nrac, nrdn
|
|
c a_rdtmp = rdfval('Affine Matrix Row 1','-')
|
|
c read(unit=a_rdtmp,fmt=*) a, b
|
|
c a_rdtmp = rdfval('Affine Matrix Row 2','-')
|
|
c read(unit=a_rdtmp,fmt=*) c, d
|
|
c a_rdtmp = rdfval('Affine Offset Vector','-')
|
|
c read(unit=a_rdtmp,fmt=*) e, f
|
|
c a_rdtmp = rdfval('Looks of the Offsets','-')
|
|
c read(unit=a_rdtmp,fmt=*) lac, ldn
|
|
c a_rdtmp = rdfval('Looks of the Image File','-')
|
|
c read(unit=a_rdtmp,fmt=*) lac0, ldn0
|
|
c a_rdtmp = rdfval('File Type','-')
|
|
c read(unit=a_rdtmp,fmt='(a)') filetype
|
|
c a_rdtmp = rdfval('Interpolation Method','-')
|
|
c read(unit=a_rdtmp,fmt='(a)') intstyle
|
|
|
|
c if(ndac .gt. CMAX) stop 'Increase column array dimension in rect'
|
|
c if(nddn .gt. RMAX) stop 'Increase row array dimension in rect'
|
|
|
|
ift = 0
|
|
psize = 8
|
|
if(index(filetype,'RMG') .ne. 0)then
|
|
ift = 1
|
|
psize = 8
|
|
write (*,*) 'Assuming RMG file type '
|
|
c For byte format *****GP 01-05******
|
|
elseif(index(filetype,'BYTE') .ne. 0)then
|
|
ift = 2
|
|
psize = 1
|
|
write (*,*) 'Assuming byte file type '
|
|
elseif(index(filetype,'REAL') .ne. 0)then
|
|
ift = 3
|
|
psize = 4
|
|
write (*,*) 'Assuming real*4 file type '
|
|
elseif(index(filetype,'DOUBLE') .ne. 0)then
|
|
ift = 4
|
|
psize = 8
|
|
write (*,*) 'Assuming double (real*8) file type '
|
|
else
|
|
write (*,*) 'Assuming complex file type '
|
|
endif
|
|
|
|
iis = 0
|
|
if(index(intstyle,'Bilinear') .ne. 0)then
|
|
iis = 1
|
|
write (*,*) 'Assuming Bilinear Interpolation '
|
|
elseif(index(intstyle,'Sinc') .ne. 0)then
|
|
iis = 2
|
|
write (*,*) 'Assuming Sinc Interpolation '
|
|
else
|
|
write (*,*) 'Assuming Nearest Neighbor '
|
|
end if
|
|
|
|
c input of resampling
|
|
if(ift .le. 1) then
|
|
ALLOCATE( rvs(0:2*ndac-1,0:nddn-1) )
|
|
else
|
|
ALLOCATE( rvs(0:ndac-1,0:nddn-1) )
|
|
end if
|
|
write(*,*) 'Allocated a map of dimension ',ndac,nddn
|
|
|
|
c variable for reading data
|
|
|
|
if(ndac .gt. nrac) then
|
|
nac = ndac
|
|
else
|
|
nac = nrac
|
|
end if
|
|
ALLOCATE( carr(0:2*nac-1) )
|
|
write(*,*) 'Allocated an array of dimension ',nac
|
|
c there is no need to allocate an array for rmg type
|
|
c For byte format *****GP 01-05******
|
|
ALLOCATE( barr(0:nac-1) )
|
|
write(*,*) 'Allocated an array of dimension ',nac
|
|
ALLOCATE( rarr(0:nac-1) )
|
|
write(*,*) 'Allocated an array of dimension ',nac
|
|
ALLOCATE( darr(0:nac-1) )
|
|
write(*,*) 'Allocated an array of dimension ',nac
|
|
|
|
|
|
c output of resampling
|
|
if(ift .le. 1) then
|
|
ALLOCATE( arr(0:2*nrac-1) )
|
|
else
|
|
ALLOCATE( arr(0:nrac-1) )
|
|
end if
|
|
write(*,*) 'Allocated array of dimension ',nrac
|
|
|
|
write (*,*) 'opening files ...'
|
|
|
|
c open files
|
|
open(11,file=infile,form='unformatted',
|
|
. access='direct',recl=psize*ndac,status='old')
|
|
open(12,file=outfile,form='unformatted',
|
|
. access='direct',recl=psize*nrac,status='unknown')
|
|
|
|
c forcing NN interpolation for byte format
|
|
c if(ift .eq. 2) then
|
|
c iis = 0
|
|
c end if
|
|
|
|
|
|
c read in the data
|
|
|
|
write (*,*) 'reading data ...'
|
|
|
|
if(ift .eq. 0) then
|
|
do j = 0 , nddn-1
|
|
if(mod(j,4096) .eq. 0) write (*,*) j
|
|
read(11,rec=j+1,iostat=ierr) (carr(k),k=0,ndac-1)
|
|
if(ierr .ne. 0) goto 999
|
|
do k = 0 , ndac -1
|
|
rvs(k,j) = real(carr(k))
|
|
rvs(k+ndac,j) = aimag(carr(k))
|
|
end do
|
|
end do
|
|
elseif(ift .eq. 1) then
|
|
do j = 0 , nddn-1
|
|
if(mod(j,4096) .eq. 0) write (*,*) j
|
|
read(11,rec=j+1,iostat=ierr) (rvs(k,j),k=0,2*ndac-1)
|
|
if(ierr .ne. 0) goto 999
|
|
end do
|
|
elseif(ift .eq. 2) then
|
|
do j = 0 , nddn-1
|
|
if(mod(j,4096) .eq. 0) write (*,*) j
|
|
read(11,rec=j+1,iostat=ierr) (barr(k),k=0,ndac-1)
|
|
if(ierr .ne. 0) goto 999
|
|
do k = 0 , ndac -1
|
|
rvs(k,j) = barr(k)
|
|
end do
|
|
end do
|
|
elseif(ift .eq. 3) then
|
|
do j = 0 , nddn-1
|
|
if(mod(j,4096) .eq. 0) write (*,*) j
|
|
read(11,rec=j+1,iostat=ierr) (rarr(k),k=0,ndac-1)
|
|
if(ierr .ne. 0) goto 999
|
|
do k = 0 , ndac -1
|
|
rvs(k,j) = rarr(k)
|
|
end do
|
|
end do
|
|
else
|
|
do j = 0 , nddn-1
|
|
if(mod(j,4096) .eq. 0) write (*,*) j
|
|
read(11,rec=j+1,iostat=ierr) (darr(k),k=0,ndac-1)
|
|
if(ierr .ne. 0) goto 999
|
|
do k = 0 , ndac -1
|
|
rvs(k,j) = darr(k)
|
|
end do
|
|
end do
|
|
end if
|
|
|
|
999 write (*,*) 'finished read ',j,' now interpolating ...'
|
|
|
|
c do the interpolation
|
|
|
|
do j = 0 , nrdn-1
|
|
if(mod(j,4096) .eq. 0) write (*,*) j
|
|
c rowval = dble(j)
|
|
c rowval = (rowval - (ldn - 1.0) / 2.0) / ldn
|
|
call look_coord_conv(ldn0, dble(j), ldn, rowval)
|
|
|
|
if(iis .eq. 0) then ! nearest neighbor
|
|
|
|
do i = 0 , nrac-1
|
|
c colval = dble(i)
|
|
c colval = (colval - (lac - 1.0) / 2.0) / lac
|
|
call look_coord_conv(lac0, dble(i), lac, colval)
|
|
ocolval_tmp = a * colval + b * rowval + e
|
|
orowval_tmp = c * colval + d * rowval + f
|
|
c ocolval = ocolval * lac + (lac - 1.0) / 2.0
|
|
c orowval = orowval * ldn + (ldn - 1.0) / 2.0
|
|
call look_coord_conv(ldn, orowval_tmp, ldn0, orowval)
|
|
call look_coord_conv(lac, ocolval_tmp, lac0, ocolval)
|
|
oi = nint(ocolval)
|
|
oj = nint(orowval)
|
|
if(.not.(oi .lt. 0 .or. oi .ge. ndac .or. oj .lt. 0 .or
|
|
$ . oj .ge. nddn)) then
|
|
arr(i) = rvs(oi,oj)
|
|
if(ift .le. 1) then
|
|
arr(i+nrac) = rvs(oi+ndac,oj)
|
|
end if
|
|
else
|
|
arr(i) = 0.
|
|
if(ift .le. 1) then
|
|
arr(i+nrac) = 0.
|
|
end if
|
|
end if
|
|
end do
|
|
|
|
elseif(iis. eq. 1) then ! bilinear interpolation
|
|
|
|
do i = 0 , nrac-1
|
|
c colval = dble(i)
|
|
c colval = (colval - (lac - 1.0) / 2.0) / lac
|
|
call look_coord_conv(lac0, dble(i), lac, colval)
|
|
ocolval_tmp = a * colval + b * rowval + e
|
|
orowval_tmp = c * colval + d * rowval + f
|
|
c ocolval = ocolval * lac + (lac - 1.0) / 2.0
|
|
c orowval = orowval * ldn + (ldn - 1.0) / 2.0
|
|
call look_coord_conv(ldn, orowval_tmp, ldn0, orowval)
|
|
call look_coord_conv(lac, ocolval_tmp, lac0, ocolval)
|
|
oi = nint(ocolval)
|
|
oj = nint(orowval)
|
|
ifrac = (ocolval - oi)
|
|
jfrac = (orowval - oj)
|
|
if(ifrac .lt. 0.d0) then
|
|
oi = oi - 1
|
|
ifrac = (ocolval - oi)
|
|
end if
|
|
if(jfrac .lt. 0.d0) then
|
|
oj = oj - 1
|
|
jfrac = (orowval - oj)
|
|
end if
|
|
if(.not.(oi .lt. 0 .or. oi .ge. ndac-1 .or. oj .lt. 0 .or
|
|
$ . oj .ge. nddn-1)) then
|
|
pt1(1) = 0.
|
|
pt1(2) = 0.
|
|
pt1(3) = rvs(oi,oj)
|
|
pt2(1) = 1.
|
|
pt2(2) = 0.
|
|
pt2(3) = rvs(oi+1,oj)
|
|
pt3(1) = 0.
|
|
pt3(2) = 1.
|
|
pt3(3) = rvs(oi,oj+1)
|
|
pt4(1) = 1.
|
|
pt4(2) = 1.
|
|
pt4(3) = rvs(oi+1,oj+1)
|
|
call bilinear(pt1,pt2,pt3,pt4,sngl(ifrac),sngl(jfrac),arr(i))
|
|
if(ift .le. 1) then
|
|
pt1(1) = 0.
|
|
pt1(2) = 0.
|
|
pt1(3) = rvs(oi+ndac,oj)
|
|
pt2(1) = 1.
|
|
pt2(2) = 0.
|
|
pt2(3) = rvs(oi+1+ndac,oj)
|
|
pt3(1) = 0.
|
|
pt3(2) = 1.
|
|
pt3(3) = rvs(oi+ndac,oj+1)
|
|
pt4(1) = 1.
|
|
pt4(2) = 1.
|
|
pt4(3) = rvs(oi+1+ndac,oj+1)
|
|
call bilinear(pt1,pt2,pt3,pt4,sngl(ifrac),sngl(jfrac),arr(i+nrac))
|
|
end if
|
|
else
|
|
arr(i) = 0.
|
|
if(ift .le. 1) then
|
|
arr(i+nrac) = 0.
|
|
end if
|
|
end if
|
|
end do
|
|
|
|
|
|
elseif(iis. eq. 2) then ! sinc interpolation
|
|
|
|
do i = 0 , nrac-1
|
|
c colval = dble(i)
|
|
c colval = (colval - (lac - 1.0) / 2.0) / lac
|
|
call look_coord_conv(lac0, dble(i), lac, colval)
|
|
ocolval_tmp = a * colval + b * rowval + e
|
|
orowval_tmp = c * colval + d * rowval + f
|
|
c ocolval = ocolval * lac + (lac - 1.0) / 2.0
|
|
c orowval = orowval * ldn + (ldn - 1.0) / 2.0
|
|
call look_coord_conv(ldn, orowval_tmp, ldn0, orowval)
|
|
call look_coord_conv(lac, ocolval_tmp, lac0, ocolval)
|
|
oi = nint(ocolval)
|
|
oj = nint(orowval)
|
|
ifrac = (ocolval - oi)
|
|
jfrac = (orowval - oj)
|
|
if(ifrac .lt. 0.d0) then
|
|
oi = oi - 1
|
|
ifrac = (ocolval - oi)
|
|
end if
|
|
if(jfrac .lt. 0.d0) then
|
|
oj = oj - 1
|
|
jfrac = (orowval - oj)
|
|
end if
|
|
|
|
! if(.not.(oi .lt. 4 .or. oi .ge. ndac-3 .or. oj .lt. 4 .or
|
|
! $ . oj .ge. nddn-3)) then
|
|
! I changed the upper sentence, as I have debug the array problem in interp.f, Cunren Liang, 03-JUN-2015
|
|
if(.not.(oi .lt. 4 .or. oi .ge. ndac-4 .or. oj .lt. 4 .or
|
|
$ . oj .ge. nddn-4)) then
|
|
arr(i) = interp(oi, oj, ifrac, jfrac, rvs, ndac, 0)
|
|
if(ift .le. 1) then
|
|
arr(i+nrac) = interp(oi, oj, ifrac, jfrac, rvs, ndac, ndac)
|
|
end if
|
|
else
|
|
arr(i) = 0.
|
|
if(ift .le. 1) then
|
|
arr(i+nrac) = 0.
|
|
end if
|
|
end if
|
|
end do
|
|
|
|
end if
|
|
|
|
if(ift .eq. 0) then
|
|
do k = 0 , nrac -1
|
|
carr(k) = cmplx(arr(k),arr(k+nrac))
|
|
end do
|
|
write(12,rec=j+1) (carr(k),k=0,nrac-1)
|
|
elseif(ift .eq. 1) then
|
|
write(12,rec=j+1) (arr(k),k=0,2*nrac-1)
|
|
elseif(ift .eq. 2) then
|
|
do k = 0 , nrac -1
|
|
barr(k) = arr(k)
|
|
end do
|
|
write(12,rec=j+1) (barr(k),k=0,nrac-1)
|
|
elseif(ift .eq. 3) then
|
|
do k = 0 , nrac -1
|
|
rarr(k) = arr(k)
|
|
end do
|
|
write(12,rec=j+1) (rarr(k),k=0,nrac-1)
|
|
else
|
|
do k = 0 , nrac -1
|
|
darr(k) = arr(k)
|
|
end do
|
|
write(12,rec=j+1) (darr(k),k=0,nrac-1)
|
|
end if
|
|
|
|
end do
|
|
|
|
DEALLOCATE(rvs)
|
|
DEALLOCATE(carr)
|
|
DEALLOCATE(barr)
|
|
DEALLOCATE(rarr)
|
|
DEALLOCATE(darr)
|
|
DEALLOCATE(arr)
|
|
|
|
close(unit=11)
|
|
close(unit=12)
|
|
end
|