Commit 7e8f0d6f authored by Daniel Scheffler's avatar Daniel Scheffler
Browse files

a couple of bug fixes and stability improvements

components.CoReg.COREG:
- __init__():
    - modified docstring (fmt_out keyword)
    - unsupported strings for output raster format are now properly catched
- _calc_shifted_cross_power_spectrum():
    - modified docstring
    - added optional figure of fft images
- _validate_ssim_improvement():
    - more stable code for handling not equal shapes of SSIM input arrays (seems to not be needed anymore)
    - bugfix for flagging tie points where SSIM kept the same (fixes issue of extensive point flagging during SSIM validity check)

components.CoReg_local.COREG_LOCAL:
- __init__():
    - modified docstring (fmt_out keyword)
    - unsupported strings for output raster format are now properly catched
- view_CoRegPoints(): added keyword 'showFig' for optinally muting figure output

components.DeShifter.DESHIFTER:
- __init__(): modified docstring (fmt_out keyword)
- correct_shifts(): bugfix for not performing resampling although the updated geotransform does not match he output grid

components.Geom_Quality_Grid.Geom_Quality_Grid:
- get_CoRegPoints_table():
    - bugfix for exception in case 'max_points' is given and actual number of points is below
    - added additional console output for number of found matches
- to_GCPList(): bugfix for crash in case no tie point passed all validity checks

components.Geom_Quality_Grid.TiePoint_Refiner:
- _RANSAC_outlier_detection():
    - bugfix for negative residual threshold)
    - empty coordinate arrays are now properly catched

- updated __version__
parent d27bf911
......@@ -9,7 +9,7 @@ from .components import utilities
from .components import geometry
__author__ = 'Daniel Scheffler'
__version__= '2016-12-13_03'
__version__= '2017-01-03_01'
__all__=['COREG',
'COREG_LOCAL',
......
......@@ -129,7 +129,8 @@ class COREG(object):
- if None (default), the method correct_shifts() does not write to disk
- if 'auto': /dir/of/im1/<im1>__shifted_to__<im0>.bsq
:param fmt_out(str): raster file format for output file. ignored if path_out is None. can be any GDAL
compatible raster file format (e.g. 'ENVI', 'GeoTIFF'; default: ENVI)
compatible raster file format (e.g. 'ENVI', 'GeoTIFF'; default: ENVI). Refer to
http://www.gdal.org/formats_list.html to get a full list of supported formats.
:param out_crea_options(list): GDAL creation options for the output image,
e.g. ["QUALITY=80", "REVERSIBLE=YES", "WRITE_METADATA=YES"]
:param r_b4match(int): band of reference image to be used for matching (starts with 1; default: 1)
......@@ -196,7 +197,7 @@ class COREG(object):
self.params = dict([x for x in locals().items() if x[0] != "self"])
# assertions
assert fmt_out, "'%s' is not a valid GDAL driver code." %fmt_out
assert gdal.GetDriverByName(fmt_out), "'%s' is not a supported GDAL driver." % fmt_out
if match_gsd and out_gsd: warnings.warn("'-out_gsd' is ignored because '-match_gsd' is set.\n")
if out_gsd: assert isinstance(out_gsd, list) and len(out_gsd) == 2, 'out_gsd must be a list with two values.'
if data_corners_ref and not isinstance(data_corners_ref[0], list): # group if not [[x,y],[x,y]..] but [x,y,x,y,]
......@@ -737,10 +738,10 @@ class COREG(object):
def _calc_shifted_cross_power_spectrum(self, im0=None, im1=None, precision=np.complex64):
"""Calculates shifted cross power spectrum for quantifying x/y-shifts.
:param im0: reference image
:param im1: subject image to shift
:param precision: to be quantified as a datatype
:return: 2D-numpy-array of the shifted cross power spectrum
:param im0: reference image
:param im1: subject image to shift
:param precision: to be quantified as a datatype
:return: 2D-numpy-array of the shifted cross power spectrum
"""
im0 = im0 if im0 is not None else self.matchWin[:] if self.matchWin.imID=='ref' else self.otherWin[:]
......@@ -785,6 +786,9 @@ class COREG(object):
fft_arr0 = np.fft.fft2(in_arr0)
fft_arr1 = np.fft.fft2(in_arr1)
#GeoArray(fft_arr0.astype(np.float32)).show(figsize=(15,15))
#GeoArray(fft_arr1.astype(np.float32)).show(figsize=(15,15))
if self.v: print('forward FFTW: %.2fs' %(time.time() -time0))
eps = np.abs(fft_arr1).max() * 1e-15
......@@ -987,8 +991,8 @@ class COREG(object):
# -> so the shift vectors have to be inverted if shifts are applied to reference image
coreg_info = self._get_inverted_coreg_info() if self.otherWin.imID=='ref' else self.coreg_info
matchFull = self.ref if self.matchWin.imID=='ref' else self.shift
otherFull = self.ref if self.otherWin.imID=='ref' else self.shift
matchFull = self.ref if self.matchWin.imID=='ref' else self.shift
otherFull = self.ref if self.otherWin.imID=='ref' else self.shift
ds_results = DESHIFTER(otherFull, coreg_info,
band2process = otherFull.band4match+1,
clipextent = list(np.array(self.matchBox.boundsMap)[[0,2,1,3]]),
......@@ -1030,14 +1034,22 @@ class COREG(object):
## check if shapes of two images are equal (due to bug (?), in some cases otherWin_deshift_geoArr does not have
## the exact same dimensions as self.matchWin -> maybe bounds are handled differently by gdal.Warp)
if not self.matchWin.shape == otherWin_deshift_geoArr.shape:
if not self.matchWin.shape == otherWin_deshift_geoArr.shape: # FIXME this seems to be already fixed
warnings.warn('SSIM input array shapes are not equal! This issue seemed to be already fixed.. ')
matchFull = self.ref if self.matchWin.imID=='ref' else self.shift
matchWinData, matchWinGt, matchWinPrj = matchFull.get_mapPos(
list(np.array(self.matchBox.boundsMap)[[0, 2, 1, 3]]), self.matchWin.prj, rspAlg='cubic',
band2get=matchFull.band4match)
self.ssim_deshifted = calc_ssim(otherWin_deshift_geoArr[:], matchWinData, dynamic_range=dr)
# at the image edges it is possible that the size of the matchBox must be reduced in order to make array shapes match
if not matchWinData.shape == otherWin_deshift_geoArr.shape: # FIXME this seems to be already fixed
self.matchBox.buffer_imXY(-np.ceil(abs(self.x_shift_px)), -np.ceil(abs(self.y_shift_px)))
otherWin_deshift_geoArr = self._get_deshifted_otherWin()
matchWinData, matchWinGt, matchWinPrj = matchFull.get_mapPos(
list(np.array(self.matchBox.boundsMap)[[0, 2, 1, 3]]), self.matchWin.prj, rspAlg='cubic',
band2get=matchFull.band4match)
self.ssim_deshifted = calc_ssim(otherWin_deshift_geoArr[:], matchWinData, dynamic_range=dr)
if v:
GeoArray(matchWinData).show()
......@@ -1048,7 +1060,7 @@ class COREG(object):
print('Image similarity within the matching window (SSIM before/after correction): %.4f => %.4f'
% (self.ssim_orig, self.ssim_deshifted))
self.ssim_improved = self.ssim_orig < self.ssim_deshifted
self.ssim_improved = self.ssim_orig <= self.ssim_deshifted
# write win data to disk
#outDir = '/home/gfz-fe/scheffler/temp/ssim_debugging/'
......
......@@ -53,7 +53,8 @@ class COREG_LOCAL(object):
- if None (default), no output is written to disk
- if 'auto': /dir/of/im1/<im1>__shifted_to__<im0>.bsq
:param fmt_out(str): raster file format for output file. ignored if path_out is None. Can be any GDAL
compatible raster file format (e.g. 'ENVI', 'GeoTIFF'; default: ENVI)
compatible raster file format (e.g. 'ENVI', 'GeoTIFF'; default: ENVI). Refer to
http://www.gdal.org/formats_list.html to get a full list of supported formats.
:param out_crea_options(list): GDAL creation options for the output image,
e.g. ["QUALITY=80", "REVERSIBLE=YES", "WRITE_METADATA=YES"]
:param projectDir(str): name of a project directory where to store all the output results. If given,
......@@ -127,7 +128,7 @@ class COREG_LOCAL(object):
"""
# assertions
assert fmt_out, "'%s' is not a valid GDAL driver code." %fmt_out
assert gdal.GetDriverByName(fmt_out), "'%s' is not a supported GDAL driver." % fmt_out
if match_gsd and out_gsd: warnings.warn("'-out_gsd' is ignored because '-match_gsd' is set.\n")
if out_gsd: assert isinstance(out_gsd, list) and len(out_gsd) == 2, 'out_gsd must be a list with two values.'
......@@ -286,7 +287,7 @@ class COREG_LOCAL(object):
def view_CoRegPoints(self, attribute2plot='ABS_SHIFT', cmap=None, exclude_fillVals=True, backgroundIm='tgt',
hide_filtered=True, figsize=None, savefigPath='', savefigDPI=96):
hide_filtered=True, figsize=None, savefigPath='', savefigDPI=96, showFig=True):
"""Shows a map of the calculated quality grid with the target image as background.
:param attribute2plot: <str> the attribute of the quality grid to be shown (default: 'ABS_SHIFT')
......@@ -299,6 +300,7 @@ class COREG_LOCAL(object):
:param figsize: <tuple> size of the figure to be viewed, e.g. (10,10)
:param savefigPath:
:param savefigDPI:
:param showFig: <bool> whether to show or to hide the figure
:return:
"""
......@@ -366,11 +368,14 @@ class COREG_LOCAL(object):
cax = divider.append_axes("right", size="2%", pad=0.1) # create axis on the right; size =2% of ax; padding = 0.1 inch
plt.colorbar(points, cax=cax)
plt.show(block=True)
if savefigPath:
fig.savefig(savefigPath, dpi=savefigDPI)
if showFig and not self.q:
plt.show(block=True)
else:
plt.close(fig)
def view_CoRegPoints_folium(self, attribute2plot='ABS_SHIFT', cmap=None, exclude_fillVals=True):
warnings.warn(UserWarning('This function is still under construction and may not work as expected!'))
......
......@@ -35,7 +35,7 @@ class DESHIFTER(object):
:Keyword Arguments:
- path_out(str): /output/directory/filename for coregistered results
- fmt_out (str): raster file format for output file. ignored if path_out is None. can be any GDAL
compatible raster file format (e.g. 'ENVI', 'GeoTIFF'; default: ENVI)
compatible raster file format (e.g. 'ENVI', 'GeoTIFF'; default: ENVI)
- out_crea_options(list): GDAL creation options for the output image,
e.g. ["QUALITY=20", "REVERSIBLE=YES", "WRITE_METADATA=YES"]
- band2process (int): The index of the band to be processed within the given array (starts with 1),
......@@ -63,7 +63,7 @@ class DESHIFTER(object):
- q(bool): quiet mode (default: False)
"""
# TODO validate kwargs for typos
self.init_args = dict([x for x in locals().items() if x[0] != "self" and not x[0].startswith('__')])
self.init_kwargs = self.init_args['kwargs']
......@@ -213,17 +213,18 @@ class DESHIFTER(object):
t_start = time.time()
equal_prj = prj_equal(self.ref_prj,self.shift_prj)
if equal_prj and is_coord_grid_equal(self.shift_gt, *self.out_grid) and not self.GCPList:
if equal_prj and not self.GCPList and is_coord_grid_equal(self.updated_gt, *self.out_grid):
"""NO RESAMPLING NEEDED"""
#print('not_warping')
self.is_shifted = True
self.is_resampled = False
xmin,ymin,xmax,ymax = self._get_out_extent()
if self.cliptoextent: # TODO validate results!
if self.cliptoextent: # TODO validate results -> output extent does not seem to be the requested one!
# get shifted array
shifted_geoArr = GeoArray(self.im2shift[:],tuple(self.updated_gt), self.shift_prj)
# clip with target extent
# clip with target extent # FIXME does get_mapPos perform a resampling?
self.arr_shifted, self.updated_gt, self.updated_projection = \
shifted_geoArr.get_mapPos((xmin,ymin,xmax,ymax), self.shift_prj, fillVal=self.nodata,
band2get=self.band2process)
......@@ -281,6 +282,7 @@ class DESHIFTER(object):
# # TO DO implement output writer
# FIXME avoid reading the whole band if clip_extent is passed
#print('warping')
in_arr = self.im2shift[:,:,self.band2process] if self.band2process is not None and self.im2shift.ndim==3\
else self.im2shift[:]
......
......@@ -227,8 +227,8 @@ class Geom_Quality_Grid(object):
else:
crs = None
GDF = GeoDataFrame(index=range(len(geomPoints)),crs=crs,
columns=['geometry','POINT_ID','X_IM','Y_IM','X_UTM','Y_UTM'])
GDF = GeoDataFrame(index=range(len(geomPoints)),crs=crs,
columns=['geometry','POINT_ID','X_IM','Y_IM','X_UTM','Y_UTM'])
GDF ['geometry'] = geomPoints
GDF ['POINT_ID'] = range(len(geomPoints))
GDF.loc[:,['X_IM' ,'Y_IM' ]] = self.XY_points
......@@ -238,7 +238,7 @@ class Geom_Quality_Grid(object):
GDF = self._exclude_bad_XYpos(GDF)
# choose a random subset of points if a maximum number has been given
if self.max_points:
if self.max_points and len(GDF) > self.max_points:
GDF = GDF.sample(self.max_points).copy()
# equalize pixel grids in order to save warping time
......@@ -314,8 +314,13 @@ class Geom_Quality_Grid(object):
GDF = GDF.merge(records, on='POINT_ID', how="inner")
GDF = GDF.fillna(int(self.outFillVal))
if not self.q:
print("Found %s matches." % len(GDF[GDF.LAST_ERR == int(self.outFillVal)]))
# filter tie points according to given filter level
if self.tieP_filter_level>0:
if not self.q:
print('Performing validity checks...')
TPR = TiePoint_Refiner(GDF[GDF.ABS_SHIFT != self.outFillVal], q=self.q)
GDF_filt, new_columns = TPR.run_filtering(level=self.tieP_filter_level)
GDF = GDF.merge(GDF_filt[ ['POINT_ID']+new_columns], on='POINT_ID', how="outer")
......@@ -347,6 +352,10 @@ class Geom_Quality_Grid(object):
GDF = GDF[GDF.OUTLIER == False].copy()
avail_TP = len(GDF)
if not avail_TP:
# no point passed all validity checks
return []
if avail_TP>7000:
GDF = GDF.sample(7000)
warnings.warn('By far not more than 7000 tie points can be used for warping within a limited '
......@@ -646,8 +655,12 @@ class TiePoint_Refiner(object):
if th_checked:
th_too_strict = count_inliers < ideal_count # True if too less inliers remaining
# calculate new theshold using old increment
th_new = th+th_substract if th_too_strict else th-th_substract
# calculate new theshold using old increment (but ensure th_new>0 by adjusting increment if needed)
th_new = 0
while th_new <= 0:
th_new = th+th_substract if th_too_strict else th-th_substract
if th_new <= 0:
th_substract /=2
# check if calculated new threshold has been used before
th_already_checked = th_new in th_checked.keys()
......@@ -659,14 +672,19 @@ class TiePoint_Refiner(object):
# RANSAC call
# model_robust, inliers = ransac((src, dst), PolynomialTransform, min_samples=3,
model_robust, inliers = \
ransac((src_coords, est_coords), AffineTransform,
min_samples = 6,
residual_threshold = th,
max_trials = 2000,
stop_sample_num = int((min_inlier_percentage-tolerance) /100*src_coords.shape[0]),
stop_residuals_sum = int((max_outlier_percentage-tolerance)/100*src_coords.shape[0])
)
if src_coords.size and est_coords.size:
model_robust, inliers = \
ransac((src_coords, est_coords), AffineTransform,
min_samples = 6,
residual_threshold = th,
max_trials = 2000,
stop_sample_num = int((min_inlier_percentage-tolerance) /100*src_coords.shape[0]),
stop_residuals_sum = int((max_outlier_percentage-tolerance)/100*src_coords.shape[0])
)
else:
inliers = np.array([])
break
count_inliers = np.count_nonzero(inliers)
th_checked[th] = count_inliers / src_coords.shape[0] * 100
......@@ -679,9 +697,11 @@ class TiePoint_Refiner(object):
count_iter+=1
outliers = inliers == False
outliers = inliers == False if inliers.size else np.array([])
if len(GDF) < len(self.GDF):
if GDF.empty:
gs = GeoSeries([False]*len(self.GDF))
elif len(GDF) < len(self.GDF):
GDF['outliers'] = outliers
fullGDF = GeoDataFrame(self.GDF['POINT_ID'])
fullGDF = fullGDF.merge(GDF[['POINT_ID', 'outliers']], on='POINT_ID', how="outer")
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment