cmd.davinci.avirisgeom.8

http://speclab.cr.usgs.gov


      PROGRAM AVRECGEN
c  ***************************************************************
c  *  Generate scanline and pixel lookup images using AVIRIS     *
c  *     engineering data for input into AVRECTFY (to rectify    *
c  *     AVIRIS single band and cube data                        *
c  *                                                             *
c  *  set band to band number to copy single band from BIL file  *
c  *    band ignored when input file is a single band,fcb(13)=1  *
c  *                                                             *
c  *  To change flightline array length, change lines:           *
c  *    INTEGER*4 slpos2(1014,6544) and pxpos(1014,6554)          *
c  *    maxsl = 6144 + (2*offset)                                *
c  *  This size is set for 12 (512sl) AVIRIS Segments for now    *
c  *                                                             *
c  *  K. Eric Livo     Sept. 1997                                *
c  ***************************************************************
      INTEGER*4 pixoff(256), sloff(256)
      INTEGER*2 pixbuf(8192), slbuf(8192)
      INTEGER*2 slpost(8192), pxpost(8192)
      INTEGER*4 slposfcb(256), pxposfcb(256)
      INTEGER*2 slpos(1014,6744), slpos2(1014,6744), pxpos(1014,6744)
      INTEGER*4 maxsl, maxpix, cursl, curpix, offset
      INTEGER*4 isl, ipix
      INTEGER*4 topsl, botsl
      INTEGER*4 leftpx, rghtpx, tmppx
      INTEGER*4 curpxout, curslout
      INTEGER*4 impixmin, impixmax, imslmin, imslmax
      INTEGER*4 ifound, xdist, ydist
      REAL*4 rsl, rpix
      DATA pixoff,sloff/3*0,2*1,251*0,3*0,2*1,251*0/

c     oversize top, left, & right sides 200 pixels
c     oversize bottom 400 scanlines (output image grows downward)
      offset = 200
      maxsl = 6144 + (3*offset)
      maxpix = 614 + (2*offset)

      impixmin = 0
      impixmax = 0
      imslmin = 0
      imslmax = 0

      CALL REMLOG('=========== AVRECGEN ===========')

c     sloff(1)=0 open read only (pc: to implement); 
c     sloff(1)=-1 open read/write

c  assumes scanline and pixel offset images enter a single band each
c  open scanline offset image
      sloff(1) = 0
      write(*,*) 'Enter SCANLINE OFFSET FILENAME'
      CALL DISKIO(0,11,slbuf,sloff)
      IF(sloff(1).LT.0) CALL EXIT(3)

c  open pixel offset image
      pixoff(1) = 0
      write(*,*) 'Enter PIXEL OFFSET FILENAME'
      CALL DISKIO(0,12,pixbuf,pixoff)
      IF(pixoff(1).LT.0) CALL EXIT(3)

C  *****
c  assuming scanline and pixel-offset images are the same dimensions
c  for now, just register single band inputs, not testing for BIL input
C  *****

c  read inputs

      cursl = 0

1     CALL DISKIO(9,11,slbuf,sloff)
      IF(sloff(1).EQ.-1) GOTO 2
      IF(sloff(1).LT.-1) CALL EXIT(3)
      CALL UNPACK(slbuf,sloff)

      CALL DISKIO(9,12,pixbuf,pixoff)
      IF(pixoff(1).EQ.-1) GOTO 2
      IF(pixoff(1).LT.-1) CALL EXIT(3)
      CALL UNPACK(pixbuf,pixoff)

      cursl = cursl + 1

c  process inputs
c  X & Y projection offsets are truncated to 'OFFSET'(200 pixel) maximum

      do 10,curpix=1,sloff(2)

         rsl = (slbuf(curpix) + 0.5) / 10.0
         isl = int(rsl)
         curslout = cursl + offset + isl

         if(curslout .lt. 1) then
            write(*,*) 'output scanline projected above top of image',
     &                 curslout
            curslout = 1
         elseif(curslout .gt. maxsl) then
            write(*,*) 'output scanline projected below bottom of',
     &                 ' image', curslout
            curslout = maxsl
         endif

         rpix = (pixbuf(curpix) + 0.5) / 10.0
         ipix = int(rpix)
         curpxout = curpix+offset+ipix

         if(curpxout .lt. 1 ) then
            write(*,*) 'output pixel projected to left of image',
     &                 curpxout
            curpxout = 1
         elseif(curpxout .gt. maxpix) then
            write(*,*) 'output pixel projected to right of image',
     &                  curpxout
            curpxout = maxpix
         endif

c  set 1st data point min and max
         if (impixmin.eq.0) impixmin=curpxout
         if (impixmax.eq.0) impixmax=curpxout
         if (imslmin.eq.0) imslmin=curslout
         if (imslmax.eq.0) imslmax=curslout

c  find min and max pix's and sl's
         if (curpxout.lt.impixmin) impixmin=curpxout
         if (curpxout.gt.impixmax) impixmax=curpxout
         if (curslout.lt.imslmin) imslmin=curslout
         if (curslout.gt.imslmax) imslmax=curslout

c  assign original-sl, and original-sample
         slpos(curpxout, curslout) = cursl
         slpos2(curpxout, curslout) = cursl
         pxpos(curpxout, curslout) = curpix

10    continue
      goto 1

2     continue

c  open output: scanline lookup (sl's of orig unrectified image)
      slposfcb(1)=16
      slposfcb(2)=(impixmax-impixmin+1)
      slposfcb(13)=1
      write(*,*) 'Enter scanline LOOKUP output image'
      CALL DISKIO(0,14,slpost,slposfcb)
      IF(slposfcb(1).LT.0) CALL EXIT(3)

c  open output: pixel lookup (pix's of orig unrectified image)
      pxposfcb(1)=16
      pxposfcb(2)=(impixmax-impixmin+1)
      pxposfcb(13)=1
      write(*,*) 'Enter pixel LOOKUP output image'
      CALL DISKIO(0,15,pxpost,pxposfcb)
      IF(pxposfcb(1).LT.0) CALL EXIT(3)

c  find DN 0's line top and bottom (1st non-zero value)
c     and left and right
c note: slpos2(pix,sl) = 0 only where holes are formed through warping
c       original image to rectified output image

      do 20, cursl=imslmin,imslmax
         do 30, curpix=impixmin,impixmax
            if(slpos2(curpix,cursl) .eq. 0) then

               topsl = cursl
               ifound = 0
               do while((ifound.eq.0).and.(topsl.gt.imslmin))
                  topsl = topsl - 1
                  ifound = slpos2(curpix,topsl)
               enddo
               if(ifound.eq.0) topsl = 0

               botsl = cursl
               ifound = 0
               do while((ifound.eq.0).and.(botsl.lt.imslmax))
                  botsl = botsl + 1
                  ifound = slpos2(curpix,botsl)
               enddo
               if(ifound.eq.0) botsl = 0

               leftpx = curpix
               ifound = 0
               do while((ifound.eq.0).and.(leftpx.gt.impixmin))
                  leftpx = leftpx - 1
                  ifound = slpos2(leftpx,cursl)
               enddo
               if(ifound.eq.0) leftpx = 0

               rghtpx = curpix
               ifound = 0
               do while((ifound.eq.0).and.(rghtpx.lt.impixmax))
                  rghtpx = rghtpx + 1
                  ifound = slpos2(rghtpx,cursl)
               enddo
               if(ifound.eq.0) rghtpx = 0

               ydist = 0
               if((topsl.gt.0).and.(botsl.gt.0)) then
                  ydist = topsl - cursl
                  if(IABS(ydist).gt.(botsl-cursl)) ydist=botsl-cursl
               endif

               xdist = 0
               if((leftpx.gt.0).and.(rghtpx.gt.0)) then
                 xdist = leftpx - curpix
                 if(IABS(xdist).gt.(rghtpx-curpix)) xdist=rghtpx-curpix
               endif

               if((xdist.ne.0).and.(ydist.ne.0)) then
                  if(IABS(xdist) .le. IABS(ydist)) then
                     slpos(curpix,cursl) = slpos2(curpix+xdist,cursl)
                     pxpos(curpix,cursl) = pxpos(curpix+xdist,cursl)
                  else
                     slpos(curpix,cursl) = slpos2(curpix,cursl+ydist)
                     pxpos(curpix,cursl) = pxpos(curpix,cursl+ydist)
                  endif

               else

                  if(xdist.ne.0) then 
                     slpos(curpix,cursl) = slpos2(curpix+xdist,cursl)
                     pxpos(curpix,cursl) = pxpos(curpix+xdist,cursl)
                  endif

                  if(ydist.ne.0) then
                     slpos(curpix,cursl) = slpos2(curpix,cursl+ydist)
                     pxpos(curpix,cursl) = pxpos(curpix,cursl+ydist)
                  endif

               endif
c  NOTE: if no match, slpos remains set to sl=0 & pxpos set to pix=0

            endif
30       continue
20    continue



c  write registered output

      do 200, cursl=imslmin,imslmax
         do 210, curpix=impixmin,impixmax
            tmppx = curpix - impixmin + 1
            slpost(tmppx) = slpos(curpix,cursl)
            pxpost(tmppx) = pxpos(curpix,cursl)
210      continue

         CALL DISKIO(10,14,slpost,slposfcb)
         IF(slposfcb(1).LT.0) CALL EXIT(4)

         CALL DISKIO(10,15,pxpost,pxposfcb)
         IF(pxposfcb(1).LT.0) CALL EXIT(4)

200   continue

c  close all files

      CALL DISKIO(6,11,slbuf,sloff)
      CALL DISKIO(6,12,pixbuf,pixoff)
      CALL DISKIO(6,14,slpost,slposfcb)
      CALL DISKIO(6,15,pxpost,pxposfcb)

      STOP
      END

U.S. Geological Survey, a bureau of the U.S. Department of the Interior
This page URL= http://speclab.cr.usgs.gov/aviris.geom.software/avrecgen.f.html
This page is maintained by: Eric Livo elivo@speclab.cr.usgs.gov
Last modified December 16, 1998.