From 6d3811ec11ca8d88cde621548b56956fbd000f7e Mon Sep 17 00:00:00 2001 From: giangi Date: Thu, 20 Jun 2019 09:55:55 -0700 Subject: [PATCH] Added downsample_snapu factory to TopsProc/Factories.py Added runner run_downsample_unwrapper.py to TopsProc added donwsample_unwraper module to contrib git commit -m --- components/isceobj/TopsProc/Factories.py | 4 +- components/isceobj/TopsProc/SConscript | 2 +- .../TopsProc/run_downsample_unwrapper.py | 33 ++++ contrib/SConscript | 2 + contrib/downsample_unwrapper/SConscript | 35 ++++ contrib/downsample_unwrapper/__init__.py | 1 + .../downsample_unwrapper.py | 173 ++++++++++++++++++ .../downsample_unwrapper/run_unwrap_snaphu.py | 91 +++++++++ 8 files changed, 339 insertions(+), 2 deletions(-) create mode 100644 components/isceobj/TopsProc/run_downsample_unwrapper.py create mode 100644 contrib/downsample_unwrapper/SConscript create mode 100644 contrib/downsample_unwrapper/__init__.py create mode 100644 contrib/downsample_unwrapper/downsample_unwrapper.py create mode 100644 contrib/downsample_unwrapper/run_unwrap_snaphu.py diff --git a/components/isceobj/TopsProc/Factories.py b/components/isceobj/TopsProc/Factories.py index c2d9f1c..744bc91 100644 --- a/components/isceobj/TopsProc/Factories.py +++ b/components/isceobj/TopsProc/Factories.py @@ -50,7 +50,9 @@ def createUnwrapper(other, do_unwrap = None, unwrapperName = None, elif unwrapperName.lower() == 'snaphu': from .runUnwrapSnaphu import runUnwrap elif unwrapperName.lower() == 'snaphu_mcf': - from .runUnwrapSnaphu import runUnwrapMcf as runUnwrap + from .runUnwrapSnaphu import runUnwrap + elif unwrapperName.lower() == 'downsample_snaphu': + from .run_downsample_unwrapper import runUnwrap elif unwrapperName.lower() == 'icu': from .runUnwrapIcu import runUnwrap elif unwrapperName.lower() == 'grass': diff --git a/components/isceobj/TopsProc/SConscript b/components/isceobj/TopsProc/SConscript index 83123ec..7bba615 100644 --- a/components/isceobj/TopsProc/SConscript +++ b/components/isceobj/TopsProc/SConscript @@ -40,6 +40,6 @@ project = 'TopsProc' install = os.path.join(envisceobj['PRJ_SCONS_INSTALL'],package,project) -listFiles = ['__init__.py', 'Factories.py', 'TopsProc.py', 'runPreprocessor.py', 'runComputeBaseline.py', 'runVerifyDEM.py', 'runTopo.py', 'runSubsetOverlaps.py', 'runCoarseOffsets.py', 'runCoarseResamp.py', 'runOverlapIfg.py', 'runPrepESD.py', 'runESD.py', 'runRangeCoreg.py', 'runFineOffsets.py', 'runFineResamp.py', 'runBurstIfg.py', 'runMergeBursts.py', 'runFilter.py', 'runUnwrapGrass.py', 'runUnwrapIcu.py', 'runUnwrapSnaphu.py', 'runGeocode.py', 'runMergeSLCs.py', 'runDenseOffsets.py', 'runOffsetFilter.py', 'runOffsetGeocode.py', 'runCropOffsetGeo.py', 'runUnwrap2Stage.py', 'VRTManager.py', 'runVerifyGeocodeDEM.py', 'runIon.py'] +listFiles = ['__init__.py', 'run_downsample_unwrapper.py','Factories.py', 'TopsProc.py', 'runPreprocessor.py', 'runComputeBaseline.py', 'runVerifyDEM.py', 'runTopo.py', 'runSubsetOverlaps.py', 'runCoarseOffsets.py', 'runCoarseResamp.py', 'runOverlapIfg.py', 'runPrepESD.py', 'runESD.py', 'runRangeCoreg.py', 'runFineOffsets.py', 'runFineResamp.py', 'runBurstIfg.py', 'runMergeBursts.py', 'runFilter.py', 'runUnwrapGrass.py', 'runUnwrapIcu.py', 'runUnwrapSnaphu.py', 'runGeocode.py', 'runMergeSLCs.py', 'runDenseOffsets.py', 'runOffsetFilter.py', 'runOffsetGeocode.py', 'runCropOffsetGeo.py', 'runUnwrap2Stage.py', 'VRTManager.py', 'runVerifyGeocodeDEM.py', 'runIon.py'] envisceobj.Install(install,listFiles) envisceobj.Alias('install',install) diff --git a/components/isceobj/TopsProc/run_downsample_unwrapper.py b/components/isceobj/TopsProc/run_downsample_unwrapper.py new file mode 100644 index 0000000..d22dd11 --- /dev/null +++ b/components/isceobj/TopsProc/run_downsample_unwrapper.py @@ -0,0 +1,33 @@ +import sys +import isceobj +import os +from isceobj.TopsProc.runUnwrapSnaphu import runUnwrapMcf +from contrib.downsample_unwrapper.downsample_unwrapper import DownsampleUnwrapper +def runUnwrap(self,costMode = None,initMethod = None, defomax = None, initOnly = None): + #generate inputs from insar obj + inps = { + "flat_name":self._insar.filtFilename, + "unw_name":self._insar.unwrappedIntFilename, + "cor_name":self._insar.coherenceFilename, + "range_looks":self.numberRangeLooks, + "azimuth_looks":self.numberAzimuthLooks, + "data_dir":self._insar.mergedDirname + } + + du = DownsampleUnwrapper(inps) + #modify the filenames so it uses the downsampled versions + self._insar.filtFilename = du._dflat_name + self._insar.unwrappedIntFilename = du._dunw_name + self._insar.coherenceFilename = du._dcor_name + self.numberRangeLooks = int(du._resamp*du._range_looks) + self.numberAzimuthLooks = int(du._resamp*du._azimuth_looks) + + du.downsample_images(du._ddir,du._flat_name,du._cor_name,du._resamp) + runUnwrapMcf(self) + du.upsample_unw(du._ddir,du._flat_name,du._dunw_name,du._dccomp_name,upsamp=du._resamp,filt_sizes=(3,4)) + #put back the original values + self._insar.filtFilename = du._flat_name + self._insar.unwrappedIntFilename = du._unw_name + self._insar.coherenceFilename = du._cor_name + self.numberRangeLooks = int(du._range_looks) + self.numberAzimuthLooks = int(du._azimuth_looks) \ No newline at end of file diff --git a/contrib/SConscript b/contrib/SConscript index f509751..fb91670 100644 --- a/contrib/SConscript +++ b/contrib/SConscript @@ -65,6 +65,8 @@ frameUtils = os.path.join('frameUtils','SConscript') SConscript(frameUtils) unwUtils = os.path.join('UnwrapComp','SConscript') SConscript(unwUtils) +downsample_unwrapper = os.path.join('downsample_unwrapper','SConscript') +SConscript(downsample_unwrapper) if 'MOTIFLIBPATH' in envcontrib.Dictionary(): mdx = os.path.join('mdx','SConscript') diff --git a/contrib/downsample_unwrapper/SConscript b/contrib/downsample_unwrapper/SConscript new file mode 100644 index 0000000..38e8b96 --- /dev/null +++ b/contrib/downsample_unwrapper/SConscript @@ -0,0 +1,35 @@ +#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# +# Giangi Sacco +# NASA Jet Propulsion Laboratory +# California Institute of Technology +# (c) 2019 All Rights Reserved +# +#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + +#!/usr/bin/env python +import os + +Import('envcontrib') +envDSU = envcontrib.Clone() +project = 'downsample_unwrapper' +package = envDSU['PACKAGE'] +envDSU['PROJECT'] = project +Export('envDSU') + +initFile = '__init__.py' +if not os.path.exists(initFile): + fout = open(initFile,"w") + fout.write("#!/usr/bin/env python") + fout.close() + +install = os.path.join(envDSU['PRJ_SCONS_INSTALL'],package,project) + +helpList,installHelp = envDSU['HELP_BUILDER'](envDSU,'__init__.py',install) +envDSU.Install(installHelp,helpList) +envDSU.Alias('install',installHelp) + +listFiles = ['downsample_unwrapper.py','run_unwrap_snaphu.py',initFile] +envDSU.Install(install,listFiles) +envDSU.Alias('install',install) diff --git a/contrib/downsample_unwrapper/__init__.py b/contrib/downsample_unwrapper/__init__.py new file mode 100644 index 0000000..fef66b5 --- /dev/null +++ b/contrib/downsample_unwrapper/__init__.py @@ -0,0 +1 @@ +#!/usr/bin/env python \ No newline at end of file diff --git a/contrib/downsample_unwrapper/downsample_unwrapper.py b/contrib/downsample_unwrapper/downsample_unwrapper.py new file mode 100644 index 0000000..1f27ac9 --- /dev/null +++ b/contrib/downsample_unwrapper/downsample_unwrapper.py @@ -0,0 +1,173 @@ +import sys +import isceobj +import os +import json +import numpy as np +from scipy.ndimage import zoom,gaussian_filter +from .run_unwrap_snaphu import runUnwrap +from iscesys.Component.Component import Component + +class DownsampleUnwrapper(Component): + def __init__(self,inps): + self._dtype_map = {'cfloat':np.complex64,'float':np.float32,'byte':np.uint8} + self._ddir = inps['data_dir'] + self._resamp = inps['resamp'] if 'resamp' in inps else 2 + self._cor_name = inps['cor_name'] if 'cor_name' in inps else 'phsig.cor' + self._unw_name = inps['unw_name'] if 'unw_name' in inps else 'filt_topophase.flat' + self._flat_name = inps['flat_name'] if 'flat_name' in inps else 'filt_topophase.unw' + self._ccomp_name = self._unw_name + '.conncomp' + self._dcor_name = '{0:s}_{1:d}x{1:d}.cor'.format(self._cor_name.replace('.cor',''),self._resamp) + self._dunw_name = '{0:s}_{1:d}x{1:d}.unw'.format(self._unw_name.replace('.unw',''),self._resamp) + self._dflat_name = '{0:s}_{1:d}x{1:d}.flat'.format(self._flat_name.replace('.flat',''),self._resamp) + self._dccomp_name = self._dunw_name + '.conncomp' + self._earth_radius = inps['earth_radius'] if 'earth_radius' in inps else 6371000 + self._altitude = inps['altitude'] if 'altitude' in inps else 800000 + self._wavelength = inps['wavelength'] if 'wavelength' in inps else 0 + self._azimuth_looks = inps['azimuth_looks'] + self._range_looks = inps['range_looks'] + self._remove_downsampled = True if 'remove_downsampled' in inps else False + + def get_isce_image(self,itype,fname,width,length): + if itype == 'flat': + im = isceobj.createIntImage() + else: + im = isceobj.createImage() + im.dataType = 'FLOAT' + if itype == 'unw': + im.bands = 2 + im.scheme = 'BIL' + elif itype == 'ccomp': + im.dataType = 'BYTE' + im.filename = fname + im.extraFilename = fname + '.vrt' + im.width = width + im.length = length + return im + + def save_image(self,ddir,fname,img,itype): + dname = os.path.join(ddir,fname) + im = self.get_isce_image(itype,dname,img.shape[-1],img.shape[0]) + img.astype(self._dtype_map[im.dataType.lower()]).tofile(dname) + im.dump(dname + '.xml') + + def remove_downsampled(self,ddir,flat,unw,phsig,ccomp): + try: + #subprocess keeps changing API. just use system + cmd = 'rm -rf {}* {}* {}* {}*'.format(os.path.join(ddir,flat),os.path.join(ddir,unw),os.path.join(ddir,phsig),os.path.join(ddir,ccomp)) + os.system(cmd) + except Exception as e: + print(e) + + def downsample_images(self,ddir,flat,phsig,resamp): + img,im = self.load_image(ddir,flat) + _,co = self.load_image(ddir,phsig) + ims = [] + cos = [] + width = img.width + length = img.length + width -= width%resamp + length -= length%resamp + for i in range(resamp): + for j in range(resamp): + ims.append(im[i:length:resamp,j:width:resamp]) + cos.append(co[i:length:resamp,j:width:resamp]) + ims = np.array(ims).transpose([1,2,0]) + cos = np.array(cos).transpose([1,2,0]) + nims = ims.mean(2) + ncos = cos.mean(2) + self.save_image(ddir, self._dcor_name,ncos,'cor') + self.save_image(ddir, self._dflat_name,nims,'flat') + + + def load_image(self,ddir,fname): + img = isceobj.createImage() + img.load(os.path.join(ddir,fname + '.xml')) + dtype = self._dtype_map[img.dataType.lower()] + width = img.getWidth() + length = img.getLength() + im = np.fromfile(os.path.join(ddir,fname),dtype) + if img.bands == 1: + im = im.reshape([length,width]) + else:#the other option is the unw which is 2 bands BIL + im = im.reshape([length,img.bands,width]) + return img,im + + def upsample_unw(self,ddir,flat,unw,ccomp,upsamp=2,filt_sizes=(3,4)): + _,dunw = self.load_image(ddir,unw) + _,flati = self.load_image(ddir,flat) + _,dccomp = self.load_image(ddir,ccomp) + uccomp = zoom(dccomp,upsamp) + uccomp = np.round(uccomp).astype(np.uint8) + ph_unw = dunw[:,1,:] + amp = np.abs(flati) + ph_flat = np.angle(flati) + uph_unw = zoom(ph_unw,upsamp) + uph_size = uph_unw.shape + ph_size = ph_flat.shape + if uph_size[0] != ph_size[0] or uph_size[0] != ph_size[1]: + #the lost one pixel during downsampling/upsampling + nunw = np.zeros(ph_flat.shape) + nunw[:uph_size[0],:uph_size[1]] = uph_unw + if ph_size[1] > uph_size[1]: + nunw[-1,:-1] = uph_unw[-1,:] + if ph_size[0] > uph_size[0]: + nunw[:-1,-1] = uph_unw[:,-1] + uph_unw = nunw + funw = self.filter_image(uph_unw,ph_flat,filt_sizes) + ifunw = np.round((funw - ph_flat)/(2*np.pi)).astype(np.int16) + funw = ph_flat + 2*np.pi*ifunw + unw_out = np.stack([amp,funw],0).transpose([1,0,2]) + self.save_image(ddir,self._unw_name,unw_out,'unw') + self.save_image(ddir,self._ccomp_name,uccomp,'ccomp') + if self._remove_downsampled: + self.remove_downsampled(ddir, self._dflat_name,self._dunw_name,self._dccomp_name,self._dccomp_name) + + def filter_image(self,unw,wrp,filt_sizes): + im0 = np.round((unw - wrp)/(2*np.pi)) + img = wrp + 2*np.pi*im0 + for filter in filt_sizes: + if not isinstance(filter,tuple): + filter = (filter,filter) + img = gaussian_filter(img,filter,0) + return wrp + 2*np.pi*np.round((img - wrp)/(2*np.pi)) + + def run_snaphu(self): + range_looks = int(self._range_looks*self._resamp) + azimuth_looks = int(self._azimuth_looks*self._resamp) + inps_json = {'flat_name':self._dflat_name,'unw_name':self._dunw_name, + 'cor_name':self._dcor_name,'wavelength':self._wavelength, + 'range_looks':range_looks,'azimuth_looks':azimuth_looks, + 'earth_radius':self._earth_radius,'altitude':self._altitude} + runUnwrap(inps_json) + +def main(inps): + """ + Run the unwrapper with downsampling + inputs: + inps = dictionary with the following key,value + { + "flat_name":"filt_topophase.flat", + "unw_name":"filt_topophase.unw", + "cor_name":"phsig.cor", + "range_looks":7, + "azimuth_looks":3, + "wavelength":0.05546576, + "resamp":2, + "data_dir":'./' + } + The range and azimuth looks are w.r.t the original image. + """ + du = DownsampleUnwrapper(inps) + #du.downsample_images(du._ddir,du._flat_name,du._cor_name,du._resamp) + #du.run_snaphu() + du.upsample_unw(du._ddir,du._flat_name,du._dunw_name,du._dccomp_name,upsamp=du._resamp,filt_sizes=(3,4)) + + +if __name__ == '__main__': + inp_json = sys.argv[1] + ddir = sys.argv[2] + inps = json.load(open(inp_json)) + inps['data_dir'] = ddir + main(inps) + + \ No newline at end of file diff --git a/contrib/downsample_unwrapper/run_unwrap_snaphu.py b/contrib/downsample_unwrapper/run_unwrap_snaphu.py new file mode 100644 index 0000000..36f983f --- /dev/null +++ b/contrib/downsample_unwrapper/run_unwrap_snaphu.py @@ -0,0 +1,91 @@ +# +# Author: Piyush Agram +# Copyright 2016 +# + +import sys +import isceobj +from contrib.Snaphu.Snaphu import Snaphu +from isceobj.Constants import SPEED_OF_LIGHT +from isceobj.Planet.Planet import Planet +import os +import json +def runUnwrap(inps_json): + costMode = 'SMOOTH' + initMethod = 'MCF' + defomax = 2.0 + initOnly = True + if isinstance(inps_json,str): + inps = json.load(open(inps_json)) + elif isinstance(inps_json,dict): + inps = inps_json + else: + print('Expecting a json filename or a dictionary') + raise ValueError + wrapName = inps['flat_name'] + unwrapName = inps['unw_name'] + img = isceobj.createImage() + img.load(wrapName + '.xml') + width = img.getWidth() + earthRadius = inps['earth_radius'] + altitude = inps['altitude'] + corrfile = inps['cor_name'] + rangeLooks = inps['range_looks'] + azimuthLooks = inps['azimuth_looks'] + wavelength = inps['wavelength'] + azfact = 0.8 + rngfact = 0.8 + corrLooks = rangeLooks * azimuthLooks/(azfact*rngfact) + maxComponents = 20 + + snp = Snaphu() + snp.setInitOnly(initOnly) + snp.setInput(wrapName) + snp.setOutput(unwrapName) + snp.setWidth(width) + snp.setCostMode(costMode) + snp.setEarthRadius(earthRadius) + snp.setWavelength(wavelength) + snp.setAltitude(altitude) + snp.setCorrfile(corrfile) + snp.setInitMethod(initMethod) + snp.setCorrLooks(corrLooks) + snp.setMaxComponents(maxComponents) + snp.setDefoMaxCycles(defomax) + snp.setRangeLooks(rangeLooks) + snp.setAzimuthLooks(azimuthLooks) + snp.setCorFileFormat('FLOAT_DATA') + snp.prepare() + snp.unwrap() + + ######Render XML + outImage = isceobj.Image.createUnwImage() + outImage.setFilename(unwrapName) + outImage.setWidth(width) + outImage.setAccessMode('read') + outImage.renderVRT() + outImage.createImage() + outImage.finalizeImage() + outImage.renderHdr() + + #####Check if connected components was created + if snp.dumpConnectedComponents: + connImage = isceobj.Image.createImage() + connImage.setFilename(unwrapName+'.conncomp') + #At least one can query for the name used + connImage.setWidth(width) + connImage.setAccessMode('read') + connImage.setDataType('BYTE') + connImage.renderVRT() + connImage.createImage() + connImage.finalizeImage() + connImage.renderHdr() + + return + + +def main(inps): + runUnwrap(inps) + +if __name__ == '__main__': + sys.exit(main(sys.argv[1])) \ No newline at end of file