Commit 33574a48 authored by Daniel Scheffler's avatar Daniel Scheffler
Browse files

added first prototype of SSIM improvement validation

components.CoReg.COREG:
- added _validate_ssim_improvement() - not yet working
- _get_image_windows_to_match(): bugfix for using wrong bands

components.CoReg_local.COREG_LOCAL:
- view_CoRegPoints(): now used the correct band of the background image

components.DeShifter.DESHIFTER:
- revised __init__()
- correct_shifts(): revised condition for diciding if warping is needed

components.Geom_Quality_Grid.Geom_Quality_Grid:
- _get_spatial_shifts(): bugfix for returning match win size XY vice versa
parent f5f610a5
......@@ -463,13 +463,13 @@ class COREG(object):
rS, rE, cS, cE = GEO.get_GeoArrayPosition_from_boxImYX(self.matchWin.boxImYX)
assert np.array_equal(np.abs(np.array([rS,rE,cS,cE])), np.array([rS,rE,cS,cE])), \
'Got negative values in gdalReadInputs for %s.' %self.matchWin.imParams.imName
self.matchWin.data = self.matchWin.imParams.GeoArray[rS:rE,cS:cE, self.matchWin.imParams.band4match-1]
self.matchWin.data = self.matchWin.imParams.GeoArray[rS:rE,cS:cE, self.matchWin.imParams.band4match]
# otherWin per subset-read einlesen
rS, rE, cS, cE = GEO.get_GeoArrayPosition_from_boxImYX(self.otherWin.boxImYX)
assert np.array_equal(np.abs(np.array([rS,rE,cS,cE])), np.array([rS,rE,cS,cE])), \
'Got negative values in gdalReadInputs for %s.' %self.otherWin.imParams.imName
self.otherWin.data = self.otherWin.imParams.GeoArray[rS:rE, cS:cE, self.otherWin.imParams.band4match - 1]
self.otherWin.data = self.otherWin.imParams.GeoArray[rS:rE, cS:cE, self.otherWin.imParams.band4match]
if self.v:
print('Original matching windows:')
......@@ -481,7 +481,6 @@ class COREG(object):
# resample otherWin.data to the resolution of matchWin AND make sure the pixel edges are identical
# (in order to make each image show the same window with the same coordinates)
# rsp_algor = 5 if is_avail_rsp_average else 2 # average if possible else cubic # OLD
# TODO replace cubic resampling by PSF resampling - average resampling leads to sinus like distortions in the fft image that make a precise coregistration impossible. Thats why there is currently no way around cubic resampling.
tgt_xmin,tgt_xmax,tgt_ymin,tgt_ymax = self.matchWin.boundsMap
self.otherWin.data = warp_ndarray(self.otherWin.data,
......@@ -719,6 +718,80 @@ class COREG(object):
return x_intshift+x_subshift, y_intshift+y_subshift
def _validate_ssim_improvement(self):
# get image dynamic range
dr = max(self.ref.win.data.max(), self.shift.win.data.max()) - \
min(self.ref.win.data.max(), self.shift.win.data.max())
# compute ssim BEFORE shift correction
from py_tools_ds.ptds.similarity.raster import calc_ssim
ssim_before = calc_ssim(self.matchWin.data, self.otherWin.data, dynamic_range=dr)
print('SSIM before', ssim_before)
#ws = int(self.matchWin.imDimsYX[1]), int(self.matchWin.imDimsYX[0])
# get shifted GeoArray in the reference image grid
shifted_geoArr = copy(self.shift.GeoArray)
geotransform = list(shifted_geoArr.gt)
geotransform[0] += self.x_shift_map
geotransform[3] += self.y_shift_map
shifted_geoArr.gt = geotransform
tgt_xmin, tgt_xmax, tgt_ymin, tgt_ymax = self.ref.win.boundsMap
arr2warp = shifted_geoArr[:,:,self.shift.band4match] if shifted_geoArr.ndim==3 else shifted_geoArr[:] # FIXME dont read complete band
from py_tools_ds.ptds.io.raster.GeoArray import get_array_at_mapPos
sub_arr, sub_gt, sub_prj = get_array_at_mapPos(arr2warp, shifted_geoArr.gt, shifted_geoArr.prj, self.ref.prj,
mapBounds=(tgt_xmin, tgt_ymin, tgt_xmax, tgt_ymax),
mapBounds_prj=self.ref.prj,
out_gsd=(5,5), # FIXME stimmt das?
rspAlg='cubic')
# # otherWin per subset-read einlesen
# rS, rE, cS, cE = GEO.get_GeoArrayPosition_from_boxImYX(self.otherWin.boxImYX)
# assert np.array_equal(np.abs(np.array([rS, rE, cS, cE])), np.array([rS, rE, cS, cE])), \
# 'Got negative values in gdalReadInputs for %s.' % self.otherWin.imParams.imName
# data2warp = shifted_geoArr[rS:rE, cS:cE, self.otherWin.imParams.band4match]
#
#
# if self.v:
# print('Original matching windows:')
# ref_data, shift_data = (self.matchWin.data, self.otherWin.data) if self.grid2use=='ref' else \
# (self.otherWin.data, self.matchWin.data)
# PLT.subplot_imshow([ref_data, shift_data],[self.ref.title,self.shift.title], grid=True)
#
# otherWin_subgt = GEO.get_subset_GeoTransform(self.otherWin.gt, self.otherWin.boxImYX)
#
# # resample otherWin.data to the resolution of matchWin AND make sure the pixel edges are identical
# # (in order to make each image show the same window with the same coordinates)
# # TODO replace cubic resampling by PSF resampling - average resampling leads to sinus like distortions in the fft image that make a precise coregistration impossible. Thats why there is currently no way around cubic resampling.
# #tgt_xmin,tgt_xmax,tgt_ymin,tgt_ymax = self.matchWin.boundsMap
# sub_arr = warp_ndarray(data2warp,
# otherWin_subgt,
# self.shift.prj,
# self.ref.prj,
# out_gsd = (self.imfft_gsd, self.imfft_gsd),
# out_bounds = ([tgt_xmin, tgt_ymin, tgt_xmax, tgt_ymax]),
# rspAlg = self.rspAlg_calc,
# in_nodata = self.shift.nodata,
# CPUs = None if self.mp else 1,
# progress = False) [0]
#out_arr, out_gt, out_prj = \
# warp_ndarray(arr, arr_gt, arr_prj, out_prj=out_prj, out_bounds=mapBounds, out_bounds_prj=mapBounds_prj,
# in_nodata=fillVal, out_nodata=fillVal, rspAlg=rspAlg, out_gsd=out_gsd)
print(sub_arr.shape)
#GeoArray(sub_arr).show(figsize=(15,15))
ssim_after = calc_ssim(sub_arr, self.otherWin.data)
print('SSIM after', ssim_after)
def calculate_spatial_shifts(self):
if self.success is False: return None,None
......
......@@ -180,12 +180,14 @@ class COREG_LOCAL(object):
# get a map showing target image
if backgroundIm not in ['tgt','ref']: raise ValueError('backgroundIm')
backgroundIm = self.im2shift if backgroundIm=='tgt' else self.imref
fig, ax, map2show = backgroundIm.show_map(figsize=figsize, nodataVal=self.nodata[1], return_map=True)
fig, ax, map2show = backgroundIm.show_map(figsize=figsize, nodataVal=self.nodata[1], return_map=True,
band=self.COREG_obj.shift.band4match)
# fig, ax, map2show = backgroundIm.show_map_utm(figsize=(20,20), nodataVal=self.nodata[1], return_map=True)
plt.title(attribute2plot)
# transform all points of quality grid to LonLat
GDF = self.CoRegPoints_table.loc[self.CoRegPoints_table.X_SHIFT_M != self.outFillVal, ['geometry', attribute2plot]].copy() \
GDF = self.CoRegPoints_table.loc\
[self.CoRegPoints_table.X_SHIFT_M != self.outFillVal, ['geometry', attribute2plot]].copy() \
if exclude_fillVals else self.CoRegPoints_table.loc[:, ['geometry', attribute2plot]]
# get LonLat coordinates for all points
......
......@@ -38,7 +38,7 @@ class DESHIFTER(object):
compatible raster file format (e.g. 'ENVI', 'GeoTIFF'; default: ENVI)
- band2process (int): The index of the band to be processed within the given array (starts with 1),
default = None (all bands are processed)
- nodata(int, float): no data value of an image to be de-shifted
- nodata(int, float): no data value of the image to be de-shifted
- out_gsd (float): output pixel size in units of the reference coordinate system (default = pixel size
of the input array), given values are overridden by match_gsd=True
- align_grids (bool): True: align the input coordinate grid to the reference (does not affect the
......@@ -62,24 +62,16 @@ class DESHIFTER(object):
"""
# unpack args
self.im2shift = im2shift if isinstance(im2shift, GeoArray) else GeoArray(im2shift)
self.shift_prj = self.im2shift.projection
self.shift_gt = list(self.im2shift.geotransform)
self.GCPList = coreg_results['GCPList'] if 'GCPList' in coreg_results else None
if not self.GCPList:
mapI = coreg_results['updated map info']
self.updated_map_info = mapI if mapI else geotransform2mapinfo(self.shift_gt, self.shift_prj)
self.original_map_info = coreg_results['original map info']
self.updated_gt = mapinfo2geotransform(self.updated_map_info) if mapI else self.shift_gt
self.ref_gt = coreg_results['reference geotransform']
self.ref_grid = coreg_results['reference grid']
self.ref_prj = coreg_results['reference projection']
self.updated_projection = self.ref_prj
# unpack kwargs
self.path_out = kwargs.get('path_out' , None)
self.fmt_out = kwargs.get('fmt_out' , 'ENVI')
self.band2process = kwargs.get('band2process', None) # starts with 1 # FIXME warum?
self.nodata = kwargs.get('nodata' , self.im2shift.nodata)
self.nodata = kwargs.get('nodata' , None) # FIXME can be replaced by self.im2shift.nodata
self.align_grids = kwargs.get('align_grids' , False)
self.rspAlg = kwargs.get('resamp_alg' , 'cubic')
self.cliptoextent = kwargs.get('cliptoextent', True)
......@@ -88,8 +80,20 @@ class DESHIFTER(object):
self.v = kwargs.get('v' , False)
self.q = kwargs.get('q' , False) if not self.v else False # overridden by v
self.progress = kwargs.get('progress' , True) if not self.q else False # overridden by q
self.out_grid = self._get_out_grid(kwargs) # needs self.ref_grid, self.im2shift
self.out_gsd = [abs(self.out_grid[0][1]-self.out_grid[0][0]), abs(self.out_grid[1][1]-self.out_grid[1][0])] # xgsd, ygsd
self.im2shift.nodata = self.nodata
self.im2shift.q = self.q
self.shift_prj = self.im2shift.projection
self.shift_gt = list(self.im2shift.geotransform)
if not self.GCPList:
mapI = coreg_results['updated map info']
self.updated_map_info = mapI if mapI else geotransform2mapinfo(self.shift_gt, self.shift_prj)
self.original_map_info = coreg_results['original map info']
self.updated_gt = mapinfo2geotransform(self.updated_map_info) if mapI else self.shift_gt
self.updated_projection = self.ref_prj
self.out_grid = self._get_out_grid(kwargs) # needs self.ref_grid, self.im2shift
self.out_gsd = [abs(self.out_grid[0][1]-self.out_grid[0][0]), abs(self.out_grid[1][1]-self.out_grid[1][0])] # xgsd, ygsd
# assertions
assert self.rspAlg in _dict_rspAlg_rsp_Int.keys(), \
......@@ -168,7 +172,7 @@ class DESHIFTER(object):
def _get_out_extent(self):
if self.cliptoextent and self.clipextent is None:
self.clipextent = self.im2shift.footprint_poly.bounds
self.clipextent = self.im2shift.footprint_poly.envelope.bounds # FIXME
else:
xmin, xmax, ymin, ymax = self.im2shift.box.boundsMap
self.clipextent = xmin, ymin, xmax, ymax
......@@ -192,7 +196,7 @@ 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.align_grids:
if equal_prj and is_coord_grid_equal(self.shift_gt, *self.out_grid) and not self.align_grids and not self.GCPList:
# FIXME buggy condition:
# reconstructable with correct_spatial_shifts from GMS
#DS = DESHIFTER(geoArr, self.coreg_info,
......
......@@ -135,7 +135,7 @@ class Geom_Quality_Grid(object):
CR = COREG(global_shared_imref, global_shared_im2shift, multiproc=False, **coreg_kwargs)
CR.calculate_spatial_shifts()
CR_res = [int(CR.matchWin.imDimsYX[0]), int(CR.matchWin.imDimsYX[1]),
CR_res = [int(CR.matchWin.imDimsYX[1]), int(CR.matchWin.imDimsYX[0]),
CR.x_shift_px, CR.y_shift_px, CR.x_shift_map, CR.y_shift_map,
CR.vec_length_map, CR.vec_angle_deg]
......
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