ISCE_INSAR/components/isceobj/Alos2Proc/runFrameMosaic.py

#
# Author: Cunren Liang
# Copyright 2015-present, NASA-JPL/Caltech
#

import os
import glob
import logging
import datetime
import numpy as np

import isceobj
from isceobj.Alos2Proc.Alos2ProcPublic import create_xml

logger = logging.getLogger('isce.alos2insar.runFrameMosaic')

def runFrameMosaic(self):
    '''mosaic frames
    '''
    if hasattr(self, 'doInSAR'):
        if not self.doInSAR:
            return

    catalog = isceobj.Catalog.createCatalog(self._insar.procDoc.name)
    self.updateParamemetersFromUser()

    referenceTrack = self._insar.loadTrack(reference=True)
    secondaryTrack = self._insar.loadTrack(reference=False)

    mosaicDir = 'insar'
    os.makedirs(mosaicDir, exist_ok=True)
    os.chdir(mosaicDir)

    numberOfFrames = len(referenceTrack.frames)
    if numberOfFrames == 1:
        import shutil
        frameDir = os.path.join('f1_{}/mosaic'.format(self._insar.referenceFrames[0]))
        if not os.path.isfile(self._insar.interferogram):
            os.symlink(os.path.join('../', frameDir, self._insar.interferogram), self._insar.interferogram)
        #shutil.copy2() can overwrite
        shutil.copy2(os.path.join('../', frameDir, self._insar.interferogram+'.vrt'), self._insar.interferogram+'.vrt')
        shutil.copy2(os.path.join('../', frameDir, self._insar.interferogram+'.xml'), self._insar.interferogram+'.xml')
        if not os.path.isfile(self._insar.amplitude):
            os.symlink(os.path.join('../', frameDir, self._insar.amplitude), self._insar.amplitude)
        shutil.copy2(os.path.join('../', frameDir, self._insar.amplitude+'.vrt'), self._insar.amplitude+'.vrt')
        shutil.copy2(os.path.join('../', frameDir, self._insar.amplitude+'.xml'), self._insar.amplitude+'.xml')

        # os.rename(os.path.join('../', frameDir, self._insar.interferogram), self._insar.interferogram)
        # os.rename(os.path.join('../', frameDir, self._insar.interferogram+'.vrt'), self._insar.interferogram+'.vrt')
        # os.rename(os.path.join('../', frameDir, self._insar.interferogram+'.xml'), self._insar.interferogram+'.xml')
        # os.rename(os.path.join('../', frameDir, self._insar.amplitude), self._insar.amplitude)
        # os.rename(os.path.join('../', frameDir, self._insar.amplitude+'.vrt'), self._insar.amplitude+'.vrt')
        # os.rename(os.path.join('../', frameDir, self._insar.amplitude+'.xml'), self._insar.amplitude+'.xml')

        #update track parameters
        #########################################################
        #mosaic size
        referenceTrack.numberOfSamples = referenceTrack.frames[0].numberOfSamples
        referenceTrack.numberOfLines = referenceTrack.frames[0].numberOfLines
        #NOTE THAT WE ARE STILL USING SINGLE LOOK PARAMETERS HERE
        #range parameters
        referenceTrack.startingRange = referenceTrack.frames[0].startingRange
        referenceTrack.rangeSamplingRate = referenceTrack.frames[0].rangeSamplingRate
        referenceTrack.rangePixelSize = referenceTrack.frames[0].rangePixelSize
        #azimuth parameters
        referenceTrack.sensingStart = referenceTrack.frames[0].sensingStart
        referenceTrack.prf = referenceTrack.frames[0].prf
        referenceTrack.azimuthPixelSize = referenceTrack.frames[0].azimuthPixelSize
        referenceTrack.azimuthLineInterval = referenceTrack.frames[0].azimuthLineInterval

        #update track parameters, secondary
        #########################################################
        #mosaic size
        secondaryTrack.numberOfSamples = secondaryTrack.frames[0].numberOfSamples
        secondaryTrack.numberOfLines = secondaryTrack.frames[0].numberOfLines
        #NOTE THAT WE ARE STILL USING SINGLE LOOK PARAMETERS HERE
        #range parameters
        secondaryTrack.startingRange = secondaryTrack.frames[0].startingRange
        secondaryTrack.rangeSamplingRate = secondaryTrack.frames[0].rangeSamplingRate
        secondaryTrack.rangePixelSize = secondaryTrack.frames[0].rangePixelSize
        #azimuth parameters
        secondaryTrack.sensingStart = secondaryTrack.frames[0].sensingStart
        secondaryTrack.prf = secondaryTrack.frames[0].prf
        secondaryTrack.azimuthPixelSize = secondaryTrack.frames[0].azimuthPixelSize
        secondaryTrack.azimuthLineInterval = secondaryTrack.frames[0].azimuthLineInterval

    else:
        #choose offsets
        if self.frameOffsetMatching:
            rangeOffsets = self._insar.frameRangeOffsetMatchingReference
            azimuthOffsets = self._insar.frameAzimuthOffsetMatchingReference
        else:
            rangeOffsets = self._insar.frameRangeOffsetGeometricalReference
            azimuthOffsets = self._insar.frameAzimuthOffsetGeometricalReference

        #list of input files
        inputInterferograms = []
        inputAmplitudes = []
        for i, frameNumber in enumerate(self._insar.referenceFrames):
            frameDir = 'f{}_{}'.format(i+1, frameNumber)
            inputInterferograms.append(os.path.join('../', frameDir, 'mosaic', self._insar.interferogram))
            inputAmplitudes.append(os.path.join('../', frameDir, 'mosaic', self._insar.amplitude))

        #note that track parameters are updated after mosaicking
        #mosaic amplitudes
        frameMosaic(referenceTrack, inputAmplitudes, self._insar.amplitude, 
            rangeOffsets, azimuthOffsets, self._insar.numberRangeLooks1, self._insar.numberAzimuthLooks1, 
            updateTrack=False, phaseCompensation=False, resamplingMethod=0)
        #mosaic interferograms
        (phaseDiffEst, phaseDiffUsed, phaseDiffSource, numberOfValidSamples) = frameMosaic(referenceTrack, inputInterferograms, self._insar.interferogram, 
            rangeOffsets, azimuthOffsets, self._insar.numberRangeLooks1, self._insar.numberAzimuthLooks1, 
            updateTrack=True, phaseCompensation=True, resamplingMethod=1)

        create_xml(self._insar.amplitude, referenceTrack.numberOfSamples, referenceTrack.numberOfLines, 'amp')
        create_xml(self._insar.interferogram, referenceTrack.numberOfSamples, referenceTrack.numberOfLines, 'int')

        catalog.addItem('frame phase diff estimated', phaseDiffEst[1:], 'runFrameMosaic')
        catalog.addItem('frame phase diff used', phaseDiffUsed[1:], 'runFrameMosaic')
        catalog.addItem('frame phase diff used source', phaseDiffSource[1:], 'runFrameMosaic')
        catalog.addItem('frame phase diff samples used', numberOfValidSamples[1:], 'runFrameMosaic')

        #update secondary parameters here
        #do not match for secondary, always use geometrical
        rangeOffsets = self._insar.frameRangeOffsetGeometricalSecondary
        azimuthOffsets = self._insar.frameAzimuthOffsetGeometricalSecondary
        frameMosaicParameters(secondaryTrack, rangeOffsets, azimuthOffsets, self._insar.numberRangeLooks1, self._insar.numberAzimuthLooks1)

    os.chdir('../')
    #save parameter file
    self._insar.saveProduct(referenceTrack, self._insar.referenceTrackParameter)
    self._insar.saveProduct(secondaryTrack, self._insar.secondaryTrackParameter)

    catalog.printToLog(logger, "runFrameMosaic")
    self._insar.procDoc.addAllFromCatalog(catalog)


def frameMosaic(track, inputFiles, outputfile, rangeOffsets, azimuthOffsets, numberOfRangeLooks, numberOfAzimuthLooks, updateTrack=False, phaseCompensation=False, phaseDiffFixed=None, snapThreshold=None, resamplingMethod=0):
    '''
    mosaic frames
    
    track:                 track
    inputFiles:            input file list
    output file:           output mosaic file
    rangeOffsets:          range offsets
    azimuthOffsets:        azimuth offsets
    numberOfRangeLooks:    number of range looks of the input files
    numberOfAzimuthLooks:  number of azimuth looks of the input files
    updateTrack:           whether update track parameters
    phaseCompensation:     whether do phase compensation for each frame
    phaseDiffFixed:        if provided, the estimated value will snap to one of these values, which is nearest to the estimated one.
    snapThreshold:         this is used with phaseDiffFixed
    resamplingMethod:      0: amp resampling. 1: int resampling. 2: slc resampling
    '''
    import numpy as np

    from contrib.alos2proc_f.alos2proc_f import rect_with_looks
    from contrib.alos2proc.alos2proc import resamp
    from isceobj.Alos2Proc.runSwathMosaic import readImage
    from isceobj.Alos2Proc.runSwathMosaic import findNonzero
    from isceobj.Alos2Proc.Alos2ProcPublic import create_xml
    from isceobj.Alos2Proc.Alos2ProcPublic import find_vrt_file
    from isceobj.Alos2Proc.Alos2ProcPublic import find_vrt_keyword
    from isceobj.Alos2Proc.Alos2ProcPublic import computePhaseDiff
    from isceobj.Alos2Proc.Alos2ProcPublic import snap

    numberOfFrames = len(track.frames)
    frames = track.frames

    rectWidth = []
    rectLength = []
    for i in range(numberOfFrames):
        infImg = isceobj.createImage()
        infImg.load(inputFiles[i]+'.xml')
        rectWidth.append(infImg.width)
        rectLength.append(infImg.length)

    #convert original offset to offset for images with looks
    #use list instead of np.array to make it consistent with the rest of the code
    rangeOffsets1 = [i/numberOfRangeLooks for i in rangeOffsets]
    azimuthOffsets1 = [i/numberOfAzimuthLooks for i in azimuthOffsets]

    #get offset relative to the first frame
    rangeOffsets2 = [0.0]
    azimuthOffsets2 = [0.0]
    for i in range(1, numberOfFrames):
        rangeOffsets2.append(0.0)
        azimuthOffsets2.append(0.0)
        for j in range(1, i+1):
            rangeOffsets2[i] += rangeOffsets1[j]
            azimuthOffsets2[i] += azimuthOffsets1[j]

    #resample each frame
    rinfs = []
    for i, inf in enumerate(inputFiles):
        rinfs.append("{}_{}{}".format(os.path.splitext(os.path.basename(inf))[0], i, os.path.splitext(os.path.basename(inf))[1]))
        #do not resample first frame
        if i == 0:
            rinfs[i] = inf
        else:
            #no need to resample
            if (abs(rangeOffsets2[i] - round(rangeOffsets2[i])) < 0.0001) and (abs(azimuthOffsets2[i] - round(azimuthOffsets2[i])) < 0.0001):
                if os.path.isfile(rinfs[i]):
                    os.remove(rinfs[i])
                os.symlink(inf, rinfs[i])
                #all of the following use of rangeOffsets2/azimuthOffsets2 is inside int(), we do the following in case it is like
                #4.99999999999...
                rangeOffsets2[i] = round(rangeOffsets2[i])
                azimuthOffsets2[i] = round(azimuthOffsets2[i])

                infImg = isceobj.createImage()
                infImg.load(inf+'.xml')
                if infImg.getImageType() == 'amp':
                    create_xml(rinfs[i], infImg.width, infImg.length, 'amp')
                else:
                    create_xml(rinfs[i], infImg.width, infImg.length, 'int')
            else:
                infImg = isceobj.createImage()
                infImg.load(inf+'.xml')
                rangeOffsets2Frac = rangeOffsets2[i] - int(rangeOffsets2[i])
                azimuthOffsets2Frac = azimuthOffsets2[i] - int(azimuthOffsets2[i])

                if resamplingMethod == 0:
                    rect_with_looks(inf,
                                    rinfs[i],
                                    infImg.width, infImg.length,
                                    infImg.width, infImg.length,
                                    1.0, 0.0,
                                    0.0, 1.0,
                                    rangeOffsets2Frac, azimuthOffsets2Frac,
                                    1,1,
                                    1,1,
                                    'COMPLEX',
                                    'Bilinear')
                    if infImg.getImageType() == 'amp':
                        create_xml(rinfs[i], infImg.width, infImg.length, 'amp')
                    else:
                        create_xml(rinfs[i], infImg.width, infImg.length, 'int')

                elif resamplingMethod == 1:
                    #decompose amplitude and phase
                    phaseFile = 'phase'
                    amplitudeFile = 'amplitude'
                    data = np.fromfile(inf, dtype=np.complex64).reshape(infImg.length, infImg.width)
                    phase = np.exp(np.complex64(1j) * np.angle(data))
                    phase[np.nonzero(data==0)] = 0
                    phase.astype(np.complex64).tofile(phaseFile)
                    amplitude = np.absolute(data)
                    amplitude.astype(np.float32).tofile(amplitudeFile)

                    #resampling
                    phaseRectFile = 'phaseRect'
                    amplitudeRectFile = 'amplitudeRect'
                    rect_with_looks(phaseFile,
                                    phaseRectFile,
                                    infImg.width, infImg.length,
                                    infImg.width, infImg.length,
                                    1.0, 0.0,
                                    0.0, 1.0,
                                    rangeOffsets2Frac, azimuthOffsets2Frac,
                                    1,1,
                                    1,1,
                                    'COMPLEX',
                                    'Sinc')
                    rect_with_looks(amplitudeFile,
                                    amplitudeRectFile,
                                    infImg.width, infImg.length,
                                    infImg.width, infImg.length,
                                    1.0, 0.0,
                                    0.0, 1.0,
                                    rangeOffsets2Frac, azimuthOffsets2Frac,
                                    1,1,
                                    1,1,
                                    'REAL',
                                    'Bilinear')

                    #recombine amplitude and phase
                    phase = np.fromfile(phaseRectFile, dtype=np.complex64).reshape(infImg.length, infImg.width)
                    amplitude = np.fromfile(amplitudeRectFile, dtype=np.float32).reshape(infImg.length, infImg.width)
                    (phase*amplitude).astype(np.complex64).tofile(rinfs[i])

                    #tidy up
                    os.remove(phaseFile)
                    os.remove(amplitudeFile)
                    os.remove(phaseRectFile)
                    os.remove(amplitudeRectFile)
                    if infImg.getImageType() == 'amp':
                        create_xml(rinfs[i], infImg.width, infImg.length, 'amp')
                    else:
                        create_xml(rinfs[i], infImg.width, infImg.length, 'int')
                else:
                    resamp(inf,
                           rinfs[i],
                           'fake',
                           'fake',
                           infImg.width, infImg.length,
                           frames[i].swaths[0].prf,
                           frames[i].swaths[0].dopplerVsPixel,
                           [rangeOffsets2Frac, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
                           [azimuthOffsets2Frac, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
                    create_xml(rinfs[i], infImg.width, infImg.length, 'slc')

    #determine output width and length
    #actually no need to calculate in azimuth direction
    xs = []
    xe = []
    ys = []
    ye = []
    for i in range(numberOfFrames):
        if i == 0:
            xs.append(0)
            xe.append(rectWidth[i] - 1)
            ys.append(0)
            ye.append(rectLength[i] - 1)
        else:
            xs.append(0 - int(rangeOffsets2[i]))
            xe.append(rectWidth[i] - 1 - int(rangeOffsets2[i]))
            ys.append(0 - int(azimuthOffsets2[i]))
            ye.append(rectLength[i] - 1 - int(azimuthOffsets2[i]))

    (xmin, xminIndex) = min((v,i) for i,v in enumerate(xs))
    (xmax, xmaxIndex) = max((v,i) for i,v in enumerate(xe))
    (ymin, yminIndex) = min((v,i) for i,v in enumerate(ys))
    (ymax, ymaxIndex) = max((v,i) for i,v in enumerate(ye))

    outWidth = xmax - xmin + 1
    outLength = ymax - ymin + 1


    #prepare for mosaicing using numpy
    xs = [x-xmin for x in xs]
    xe = [x-xmin for x in xe]
    ys = [y-ymin for y in ys]
    ye = [y-ymin for y in ye]


    #compute phase offset
    if phaseCompensation:

        phaseDiffEst    = [0.0 for i in range(numberOfFrames)]
        phaseDiffUsed   = [0.0 for i in range(numberOfFrames)]
        phaseDiffSource = ['estimated' for i in range(numberOfFrames)]
        numberOfValidSamples = [0 for i in range(numberOfFrames)]
        #phaseDiffEst    = [0.0]
        #phaseDiffUsed   = [0.0]
        #phaseDiffSource = ['estimated']

        phaseOffsetPolynomials = [np.array([0.0])]
        for i in range(1, numberOfFrames):
            upperframe = np.zeros((ye[i-1]-ys[i]+1, outWidth), dtype=np.complex128)
            lowerframe = np.zeros((ye[i-1]-ys[i]+1, outWidth), dtype=np.complex128)
            #upper frame
            if os.path.isfile(rinfs[i-1]):
                upperframe[:,xs[i-1]:xe[i-1]+1] = readImage(rinfs[i-1], rectWidth[i-1], rectLength[i-1], 0, rectWidth[i-1]-1, ys[i]-ys[i-1], ye[i-1]-ys[i-1])
            else:
                upperframe[:,xs[i-1]:xe[i-1]+1] = readImageFromVrt(rinfs[i-1], 0, rectWidth[i-1]-1, ys[i]-ys[i-1], ye[i-1]-ys[i-1])
            #lower frame
            if os.path.isfile(rinfs[i]):
                lowerframe[:,xs[i]:xe[i]+1] = readImage(rinfs[i], rectWidth[i], rectLength[i], 0, rectWidth[i]-1, 0, ye[i-1]-ys[i])
            else:
                lowerframe[:,xs[i]:xe[i]+1] = readImageFromVrt(rinfs[i], 0, rectWidth[i]-1, 0, ye[i-1]-ys[i])
            #get a polynomial
            diff = np.sum(upperframe * np.conj(lowerframe), axis=0)
            (firstLine, lastLine, firstSample, lastSample) = findNonzero(np.reshape(diff, (1, outWidth)))
            #here i use mean value(deg=0) in case difference is around -pi or pi.
            #!!!!!there have been updates, now deg must be 0
            deg = 0
            p = np.polyfit(np.arange(firstSample, lastSample+1), np.angle(diff[firstSample:lastSample+1]), deg)

            #need to use a more sophisticated method to compute the mean phase difference
            (phaseDiffEst[i], numberOfValidSamples[i]) = computePhaseDiff(upperframe, lowerframe, coherenceWindowSize=9, coherenceThreshold=0.80)

            #snap phase difference to fixed values
            if phaseDiffFixed is not None:
                (outputValue, snapped) = snap(phaseDiffEst[i], phaseDiffFixed, snapThreshold)
                if snapped == True:
                    phaseDiffUsed[i] = outputValue
                    phaseDiffSource[i] = 'estimated+snap'
                else:
                    phaseDiffUsed[i] = phaseDiffEst[i]
                    phaseDiffSource[i] = 'estimated'
            else:
                phaseDiffUsed[i] = phaseDiffEst[i]
                phaseDiffSource[i] = 'estimated'

            #use new phase constant value
            p[-1] = phaseDiffUsed[i]

            phaseOffsetPolynomials.append(p)


            #check fit result
            DEBUG = False
            if DEBUG:
                #create a dir and work in this dir
                diffDir = 'frame_mosaic'
                os.makedirs(diffDir, exist_ok=True)
                os.chdir(diffDir)

                #dump phase difference
                diffFilename = 'phase_difference_frame{}-frame{}.int'.format(i, i+1)
                (upperframe * np.conj(lowerframe)).astype(np.complex64).tofile(diffFilename)
                create_xml(diffFilename, outWidth, ye[i-1]-ys[i]+1, 'int')

                #plot phase difference vs range
                import matplotlib.pyplot as plt
                x = np.arange(firstSample, lastSample+1)
                y = np.angle(diff[firstSample:lastSample+1])
                plt.plot(x, y, label='original phase difference')
                plt.plot(x, np.polyval(p, x), label='fitted phase difference')
                plt.legend()

                plt.minorticks_on()
                plt.tick_params('both', length=10, which='major')
                plt.tick_params('both', length=5, which='minor')

                plt.xlabel('Range Sample Number [Samples]')
                plt.ylabel('Phase Difference [Rad]')
                plt.savefig('phase_difference_frame{}-frame{}.pdf'.format(i, i+1))

                os.chdir('../')


    #mosaic file
    outFp = open(outputfile,'wb')
    for i in range(numberOfFrames):
        print('adding frame: {}'.format(i+1))

        #phase offset in the polynomials
        if phaseCompensation:
            cJ = np.complex64(1j)
            phaseOffset = np.ones(outWidth, dtype=np.complex64)
            for j in range(i+1):
                phaseOffset *= np.exp(cJ*np.polyval(phaseOffsetPolynomials[j], np.arange(outWidth)))

        #get start line number (starts with zero)
        if i == 0:
            ys1 = 0
        else:
            ys1 = int((ye[i-1]+ys[i])/2.0) + 1 - ys[i]
        #get end line number (start with zero)
        if i == numberOfFrames-1:
            ye1 = rectLength[i] - 1
        else:
            ye1 = int((ye[i]+ys[i+1])/2.0) - ys[i]

        #get image format
        inputimage = find_vrt_file(rinfs[i]+'.vrt', 'SourceFilename', relative_path=True)
        byteorder = find_vrt_keyword(rinfs[i]+'.vrt', 'ByteOrder')
        if byteorder == 'LSB':
            swapByte = False
        else:
            swapByte = True
        imageoffset = int(find_vrt_keyword(rinfs[i]+'.vrt', 'ImageOffset'))
        lineoffset = int(find_vrt_keyword(rinfs[i]+'.vrt', 'LineOffset'))

        #read image
        with open(inputimage,'rb') as fp:
            for j in range(ys1, ye1+1):
                fp.seek(imageoffset+j*lineoffset, 0)
                data = np.zeros(outWidth, dtype=np.complex64)
                if swapByte:
                    tmp = np.fromfile(fp, dtype='>f', count=2*rectWidth[i])
                    cJ = np.complex64(1j)
                    data[xs[i]:xe[i]+1] = tmp[0::2] + cJ * tmp[1::2]
                else:
                    data[xs[i]:xe[i]+1] = np.fromfile(fp, dtype=np.complex64, count=rectWidth[i])
                if phaseCompensation:
                    data *= phaseOffset
                data.astype(np.complex64).tofile(outFp)
    outFp.close()


    #delete files. DO NOT DELETE THE FIRST ONE!!!
    for i in range(numberOfFrames):
        if i == 0:
            continue
        os.remove(rinfs[i])
        os.remove(rinfs[i]+'.vrt')
        os.remove(rinfs[i]+'.xml')


    #update frame parameters
    if updateTrack:
        #mosaic size
        track.numberOfSamples = outWidth
        track.numberOfLines = outLength
        #NOTE THAT WE ARE STILL USING SINGLE LOOK PARAMETERS HERE
        #range parameters
        track.startingRange = frames[0].startingRange + (int(rangeOffsets2[0]) - int(rangeOffsets2[xminIndex])) * numberOfRangeLooks * frames[0].rangePixelSize
        track.rangeSamplingRate = frames[0].rangeSamplingRate
        track.rangePixelSize = frames[0].rangePixelSize
        #azimuth parameters
        track.sensingStart = frames[0].sensingStart
        track.prf = frames[0].prf
        track.azimuthPixelSize = frames[0].azimuthPixelSize
        track.azimuthLineInterval = frames[0].azimuthLineInterval

    if phaseCompensation:
        # estimated phase diff, used phase diff, used phase diff source
        return (phaseDiffEst, phaseDiffUsed, phaseDiffSource, numberOfValidSamples)


def frameMosaicParameters(track, rangeOffsets, azimuthOffsets, numberOfRangeLooks, numberOfAzimuthLooks):
    '''
    mosaic frames (this simplified version of frameMosaic to only update parameters)
    
    track:                 track
    rangeOffsets:          range offsets
    azimuthOffsets:        azimuth offsets
    numberOfRangeLooks:    number of range looks of the input files
    numberOfAzimuthLooks:  number of azimuth looks of the input files
    '''

    numberOfFrames = len(track.frames)
    frames = track.frames

    rectWidth = []
    rectLength = []
    for i in range(numberOfFrames):
        rectWidth.append(frames[i].numberOfSamples)
        rectLength.append(frames[i].numberOfLines)

    #convert original offset to offset for images with looks
    #use list instead of np.array to make it consistent with the rest of the code
    rangeOffsets1 = [i/numberOfRangeLooks for i in rangeOffsets]
    azimuthOffsets1 = [i/numberOfAzimuthLooks for i in azimuthOffsets]

    #get offset relative to the first frame
    rangeOffsets2 = [0.0]
    azimuthOffsets2 = [0.0]
    for i in range(1, numberOfFrames):
        rangeOffsets2.append(0.0)
        azimuthOffsets2.append(0.0)
        for j in range(1, i+1):
            rangeOffsets2[i] += rangeOffsets1[j]
            azimuthOffsets2[i] += azimuthOffsets1[j]

    #determine output width and length
    #actually no need to calculate in azimuth direction
    xs = []
    xe = []
    ys = []
    ye = []
    for i in range(numberOfFrames):
        if i == 0:
            xs.append(0)
            xe.append(rectWidth[i] - 1)
            ys.append(0)
            ye.append(rectLength[i] - 1)
        else:
            xs.append(0 - int(rangeOffsets2[i]))
            xe.append(rectWidth[i] - 1 - int(rangeOffsets2[i]))
            ys.append(0 - int(azimuthOffsets2[i]))
            ye.append(rectLength[i] - 1 - int(azimuthOffsets2[i]))

    (xmin, xminIndex) = min((v,i) for i,v in enumerate(xs))
    (xmax, xmaxIndex) = max((v,i) for i,v in enumerate(xe))
    (ymin, yminIndex) = min((v,i) for i,v in enumerate(ys))
    (ymax, ymaxIndex) = max((v,i) for i,v in enumerate(ye))

    outWidth = xmax - xmin + 1
    outLength = ymax - ymin + 1

    #update frame parameters
    #mosaic size
    track.numberOfSamples = outWidth
    track.numberOfLines = outLength
    #NOTE THAT WE ARE STILL USING SINGLE LOOK PARAMETERS HERE
    #range parameters
    track.startingRange = frames[0].startingRange + (int(rangeOffsets2[0]) - int(rangeOffsets2[xminIndex])) * numberOfRangeLooks * frames[0].rangePixelSize
    track.rangeSamplingRate = frames[0].rangeSamplingRate
    track.rangePixelSize = frames[0].rangePixelSize
    #azimuth parameters
    track.sensingStart = frames[0].sensingStart
    track.prf = frames[0].prf
    track.azimuthPixelSize = frames[0].azimuthPixelSize
    track.azimuthLineInterval = frames[0].azimuthLineInterval


def readImageFromVrt(inputfile, startSample, endSample, startLine, endLine):
    '''
    read a chunk of image
    the indexes (startSample, endSample, startLine, endLine) are included and start with zero

    memmap is not used, because it is much slower

    tested against readImage in runSwathMosaic.py
    '''
    import os
    from isceobj.Alos2Proc.Alos2ProcPublic import find_vrt_keyword
    from isceobj.Alos2Proc.Alos2ProcPublic import find_vrt_file

    inputimage = find_vrt_file(inputfile+'.vrt', 'SourceFilename', relative_path=True)
    byteorder = find_vrt_keyword(inputfile+'.vrt', 'ByteOrder')
    if byteorder == 'LSB':
        swapByte = False
    else:
        swapByte = True
    imageoffset = int(find_vrt_keyword(inputfile+'.vrt', 'ImageOffset'))
    lineoffset = int(find_vrt_keyword(inputfile+'.vrt', 'LineOffset'))

    data = np.zeros((endLine-startLine+1, endSample-startSample+1), dtype=np.complex64)
    with open(inputimage,'rb') as fp:
        #fp.seek(imageoffset, 0)
        #for i in range(endLine-startLine+1):
        for i in range(startLine, endLine+1):
            fp.seek(imageoffset+i*lineoffset+startSample*8, 0)
            if swapByte:
                tmp = np.fromfile(fp, dtype='>f', count=2*(endSample-startSample+1))
                cJ = np.complex64(1j)
                data[i-startLine] = tmp[0::2] + cJ * tmp[1::2]
            else:
                data[i-startLine] = np.fromfile(fp, dtype=np.complex64, count=endSample-startSample+1)
    return data