!-----------------------------------------------------------------!
!								  !
!								  !
!    			Batch Reproject				  !
!							          !
!    This script is used to batch reproject a set of files        !
!    to a projection of an input image file.                      !
!								  !
!    A list of file is necessary.  It is to include the file      !
!    with the source projection, names of input and output files. !
!								  !
!    Please refer to the Readme file for instructions on          !
!    running this script.					  !
!								  !
!-----------------------------------------------------------------!




global GeoInfo sNewGeoInfo
global string sDirectory
global string sOutputType

local mstring FileList
local mstring IOFileList
local int i, fd

Rem !!!!!!! This must be set to indicate the type of file to export 
Rem         and also the directory for the necessary files

sOutputType = "PIX"
sDirectory = NormalizePath( "C:\Program Files\PCI Geomatics\Geomatica_V103\user\\" )

print "Start conversion process"

FileList = TEXT$Import(sDirectory + "filelist.txt")


rem Use nearest neighbor resmapling
    RESAMPLE="near"
rem No rubber-sheeting between GCPs
    ORDER=0            
rem No GCPs input
    DBGC=   
rem Default is segment 1 on FILI, this is the georeferencing for the input 
    INGEO=    
rem Default is segment 1 on FILO, this is the georeferencing for the output
    OUTGEO=      
    MEMSIZE=


for i = 1 to F$Len(FileList)

    IOFileList = tokenize(FileList(i))
    if F$len(IOFileList) = 2 then
        print "Reproject files : " + FileList(i)
        call ReprojectImage(IOFileList(1), IOFileList(2))
    elseif F$len(IOFileList) = 1 then
Rem Read proper projection from a temporary file
        print "Set projection : " + FileList(i)
        fd = DBOpen(IOFileList(1), "r")
        if fd < 1 then
            print "-----Error openning file : " + IOFileList(1)
        endif
        call DBReadGeoInfo(fd, sNewGeoInfo)
        call DBClose(fd)
    else
        print "Error on line " + F$string(i) + " must be two file names"
    endif

endfor

define function ReprojectImage(sFileIn, sFileOut)

local int fd, i, nPCT
local byte PCT[768]
local int nLines, nPixels, nChannels
local string sDataType
local GeoInfo sGeoInfoIn
local double dfCorners(4)
local string sSourceTempFIle, sDestinationTempFile

rem Simple check to see if input file exists
    fd = DBOpen(sFileIn, "r")
    if fd < 1 then
        print "-----Error openning file : " + sFileIn
        return
    endif
    nChannels = DBChannels(fd)
    nLines = DBLines(fd)
    nPixels = DBPixels(fd)
    sDataType = DBChanType(fd, 1)
    nPCT = DBNextSeg(fd, "PCT", -1)
    if nPCT > -1 then
        call DBReadPCT(fd, nPCT, PCT)
    endif
    call DBClose(fd)

    sSourceTempFIle = sDirectory + "sSourceTempPix.pix"
    sDestinationTempFIle = sDirectory + "sDestinationTempPix.pix"

rem Ensure temp file are deleted
    call DBDelete("PIX", sSourceTempFile)
    call DBDelete("PIX", sDestinationTempFile)

rem Copy over the source images to a temp file
    FILI=sFileIn
    FILO=sSourceTempFile
    
    DBIW=
    DBVS=
    DBLUT=
    DBIB=
    DBIC = 1
    DBOC = 1 
    DBPCT = 
    FTYPE="PIX"
    for i = 2 to nChannels
        DBIC(i) = i
        DBOC(i) = i
    endfor
    run FEXPORT

rem Update the projection information on the source temp file 
    print "    Converting file : " + sFileIn
    fd = DBOpen(sSourceTempFile, "r+")
    if fd < 1 then
        print "-----Error openning file : " + sSourceTempFile
        return
    endif

    call DBReadGeoInfo(fd, sGeoInfoIn)

    call DBClose(fd)
         
rem Create new output temp file with same dimensions as input
    print "    Create new file : " + sFileOut
    call DBDelete(FType, sFileOut)
    fd = DBCreate(sDestinationTempFile, "PIX", nPixels, nLines, nChannels, sDataType, "")

rem Create new PCT segment and write the PCT
    if nPCT > -1 then
        nPCT = DBSegCreate(fd, "", "", "PCT")
        call DBWritePCT(fd, nPCT, PCT)
    endif

rem Convert input corners in projection to lon/lat
    dfCorners(1) = sGeoInfoIn.ULX
    dfCorners(2) = sGeoInfoIn.ULY
    dfCorners(3) = sGeoInfoIn.LRX
    dfCorners(4) = sGeoInfoIn.LRY

    call Reproject(sGeoInfoIn, sNewGeoInfo, 2, dfCorners)
    sNewGeoInfo.XOff = dfCorners(1)
    sNewGeoInfo.YOff = dfCorners(2)
    sNewGeoInfo.XSize = (dfCorners(3) - dfCorners(1))/nPixels
    sNewGeoInfo.YSize = (dfCorners(4) - dfCorners(2))/nLines
    sNewGeoInfo.XRot= 0
    sNewGeoInfo.YRot = 0

    call DBWriteGeoInfo(fd, sNewGeoInfo)
    call DBClose(fd)

rem Perform the reprojection
    FILI=sSourceTempFile
    FILO=sDestinationTempFile
    RUN REGPRO

rem Copy reprojected file to destination
    FOPTIONS=""
    FILI=sDestinationTempFIle
    FILO=sFileOut
    FTYPE=sOutputType
    for i = 1 to nChannels
        DBIC(i) = i
    endfor
    if nPCT > -1 then
        DBPCT = nPCT
    endif
    run FEXPORT
    FOPTIONS=""

rem Ensure temp file are deleted
    call DBDelete("PIX", sSourceTempFile)
    call DBDelete("PIX", sDestinationTempFile)
enddefine