Tie_Point_Grid.py 50.6 KB
Newer Older
1
2
# -*- coding: utf-8 -*-

3
4
# AROSICS - Automated and Robust Open-Source Image Co-Registration Software
#
5
# Copyright (C) 2017-2020  Daniel Scheffler (GFZ Potsdam, daniel.scheffler@gfz-potsdam.de)
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
#
# This software was developed within the context of the GeoMultiSens project funded
# by the German Federal Ministry of Education and Research
# (project grant code: 01 IS 14 010 A-C).
#
# This program is free software: you can redistribute it and/or modify it under
# the terms of the GNU Lesser General Public License as published by the Free
# Software Foundation, either version 3 of the License, or (at your option) any
# later version.
#
# This program is distributed in the hope that it will be useful, but WITHOUT
# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
# FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
# details.
#
# You should have received a copy of the GNU Lesser General Public License along
# with this program.  If not, see <http://www.gnu.org/licenses/>.

24
25
26
27
import collections
import multiprocessing
import os
import warnings
28
import time
29
30

# custom
31
from osgeo import gdal
32
import numpy as np
33
34
from geopandas import GeoDataFrame
from pandas import DataFrame, Series
35
from shapely.geometry import Point
36
37

# internal modules
38
from .CoReg import COREG
39
from py_tools_ds.geo.projection import isProjectedOrGeographic, isLocal, get_UTMzone
40
from py_tools_ds.io.pathgen import get_generic_outpath
41
from py_tools_ds.processing.progress_mon import ProgressBar
42
from py_tools_ds.geo.vector.conversion import points_to_raster
43
from py_tools_ds.io.vector.writer import write_shp
44
from geoarray import GeoArray
45

46
from .CoReg import GeoArray_CoReg  # noqa F401  # flake8 issue
47

48
__author__ = 'Daniel Scheffler'
49

50
global_shared_imref = None
51
52
53
global_shared_im2shift = None


54
def mp_initializer(imref, imtgt):
Daniel Scheffler's avatar
Daniel Scheffler committed
55
    """Declare global variables needed for self._get_spatial_shifts().
56
57
58
59
60
61
62
63
64

    :param imref:   reference image
    :param imtgt:   target image
    """
    global global_shared_imref, global_shared_im2shift
    global_shared_imref = imref
    global_shared_im2shift = imtgt


65
class Tie_Point_Grid(object):
Daniel Scheffler's avatar
Daniel Scheffler committed
66
67
68
69
70
71
72
73
74
75
    """
    The 'Tie_Point_Grid' class applies the algorithm to detect spatial shifts to the overlap area of the input images.

    Spatial shifts are calculated for each point in grid of which the parameters can be adjusted using keyword
    arguments. Shift correction performs a polynomial transformation using te calculated shifts of each point in the
    grid as GCPs. Thus 'Tie_Point_Grid' can be used to correct for locally varying geometric distortions of the target
    image.

    See help(Tie_Point_Grid) for documentation!
    """
76

77
    def __init__(self, COREG_obj, grid_res, max_points=None, outFillVal=-9999, resamp_alg_calc='cubic',
78
79
                 tieP_filter_level=3, outlDetect_settings=None, dir_out=None, CPUs=None, progress=True, v=False,
                 q=False):
Daniel Scheffler's avatar
Daniel Scheffler committed
80
        """Get an instance of the 'Tie_Point_Grid' class.
81

82
        :param COREG_obj(object):       an instance of COREG class
83
        :param grid_res:                grid resolution in pixels of the target image (x-direction)
84
        :param max_points(int):         maximum number of points used to find coregistration tie points
85
86
87
                                        NOTE: Points are selected randomly from the given point grid (specified by
                                        'grid_res'). If the point does not provide enough points, all available points
                                        are chosen.
Daniel Scheffler's avatar
Daniel Scheffler committed
88
        :param outFillVal(int):         if given the generated tie points grid is filled with this value in case
89
                                        no match could be found during co-registration (default: -9999)
90
91
        :param resamp_alg_calc(str)     the resampling algorithm to be used for all warping processes during calculation
                                        of spatial shifts
92
93
                                        (valid algorithms: nearest, bilinear, cubic, cubic_spline, lanczos, average,
                                                           mode, max, min, med, q1, q3)
94
                                        default: cubic (highly recommended)
95
        :param tieP_filter_level(int):  filter tie points used for shift correction in different levels (default: 3).
96
                                        NOTE: lower levels are also included if a higher level is chosen
97
                                            - Level 0: no tie point filtering
98
99
100
                                            - Level 1: Reliablity filtering - filter all tie points out that have a low
                                                reliability according to internal tests
                                            - Level 2: SSIM filtering - filters all tie points out where shift
101
102
                                                correction does not increase image similarity within matching window
                                                (measured by mean structural similarity index)
103
                                            - Level 3: RANSAC outlier detection
104
105
106
107
        :param outlDetect_settings      a dictionary with the settings to be passed to
                                        arosics.TiePointGrid.Tie_Point_Refiner. Available keys: min_reliability,
                                        rs_max_outlier, rs_tolerance, rs_max_iter, rs_exclude_previous_outliers,
                                        rs_timeout, q. See documentation there.
108
109
        :param dir_out(str):            output directory to be used for all outputs if nothing else is given
                                        to the individual methods
Daniel Scheffler's avatar
Daniel Scheffler committed
110
        :param CPUs(int):               number of CPUs to use during calculation of tie points grid
111
                                        (default: None, which means 'all CPUs available')
112
        :param progress(bool):          show progress bars (default: True)
113
114
        :param v(bool):                 verbose mode (default: False)
        :param q(bool):                 quiet mode (default: False)
115
        """
116
117
        if not isinstance(COREG_obj, COREG):
            raise ValueError("'COREG_obj' must be an instance of COREG class.")
118

119
        self.COREG_obj = COREG_obj  # type: COREG
120
121
122
123
        self.grid_res = grid_res
        self.max_points = max_points
        self.outFillVal = outFillVal
        self.rspAlg_calc = resamp_alg_calc
124
        self.tieP_filter_level = tieP_filter_level
125
        self.outlDetect_settings = outlDetect_settings or dict()
126
127
128
129
130
        self.dir_out = dir_out
        self.CPUs = CPUs
        self.v = v
        self.q = q if not v else False  # overridden by v
        self.progress = progress if not q else False  # overridden by q
131

132
133
134
        if 'q' not in self.outlDetect_settings:
            self.outlDetect_settings['q'] = self.q

135
136
        self.ref = self.COREG_obj.ref  # type: GeoArray_CoReg
        self.shift = self.COREG_obj.shift  # type: GeoArray_CoReg
137

138
        self.XY_points, self.XY_mapPoints = self._get_imXY__mapXY_points(self.grid_res)
139
140
141
142
        self._CoRegPoints_table = None  # set by self.CoRegPoints_table
        self._GCPList = None  # set by self.to_GCPList()
        self.kriged = None  # set by Raster_using_Kriging()

143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
    @property
    def mean_x_shift_px(self):
        return self.CoRegPoints_table['X_SHIFT_PX'][self.CoRegPoints_table['X_SHIFT_PX'] != self.outFillVal].mean()

    @property
    def mean_y_shift_px(self):
        return self.CoRegPoints_table['Y_SHIFT_PX'][self.CoRegPoints_table['Y_SHIFT_PX'] != self.outFillVal].mean()

    @property
    def mean_x_shift_map(self):
        return self.CoRegPoints_table['X_SHIFT_M'][self.CoRegPoints_table['X_SHIFT_M'] != self.outFillVal].mean()

    @property
    def mean_y_shift_map(self):
        return self.CoRegPoints_table['Y_SHIFT_M'][self.CoRegPoints_table['Y_SHIFT_M'] != self.outFillVal].mean()
158

159
160
    @property
    def CoRegPoints_table(self):
Daniel Scheffler's avatar
Daniel Scheffler committed
161
162
163
164
        """Return a GeoDataFrame containing all the results from coregistration for all points in the tie point grid.

        Columns of the GeoDataFrame: 'geometry','POINT_ID','X_IM','Y_IM','X_UTM','Y_UTM','X_WIN_SIZE', 'Y_WIN_SIZE',
                                     'X_SHIFT_PX','Y_SHIFT_PX', 'X_SHIFT_M', 'Y_SHIFT_M', 'ABS_SHIFT' and 'ANGLE'
165
        """
166
167
168
169
170
171
172
173
174
175
176
177
        if self._CoRegPoints_table is not None:
            return self._CoRegPoints_table
        else:
            self._CoRegPoints_table = self.get_CoRegPoints_table()
            return self._CoRegPoints_table

    @CoRegPoints_table.setter
    def CoRegPoints_table(self, CoRegPoints_table):
        self._CoRegPoints_table = CoRegPoints_table

    @property
    def GCPList(self):
Daniel Scheffler's avatar
Daniel Scheffler committed
178
        """Return a list of GDAL compatible GCP objects."""
179
180
181
182
        if self._GCPList:
            return self._GCPList
        else:
            self._GCPList = self.to_GCPList()
183
            return self._GCPList
184
185
186
187
188
189

    @GCPList.setter
    def GCPList(self, GCPList):
        self._GCPList = GCPList

    def _get_imXY__mapXY_points(self, grid_res):
Daniel Scheffler's avatar
Daniel Scheffler committed
190
191
192
        """Return a numpy array containing possible positions for coregistration tie points.

        NOTE: The returned positions are dependent from the given grid resolution.
193
194
195
196

        :param grid_res:
        :return:
        """
197
        if not self.q:
Daniel Scheffler's avatar
Daniel Scheffler committed
198
            print('Initializing tie points grid...')
199

200
201
        Xarr, Yarr = np.meshgrid(np.arange(0, self.shift.shape[1], grid_res),
                                 np.arange(0, self.shift.shape[0], grid_res))
202

203
204
        mapXarr = np.full_like(Xarr, self.shift.gt[0], dtype=np.float64) + Xarr * self.shift.gt[1]
        mapYarr = np.full_like(Yarr, self.shift.gt[3], dtype=np.float64) - Yarr * abs(self.shift.gt[5])
205

206
207
208
        XY_points = np.empty((Xarr.size, 2), Xarr.dtype)
        XY_points[:, 0] = Xarr.flat
        XY_points[:, 1] = Yarr.flat
209

210
211
212
        XY_mapPoints = np.empty((mapXarr.size, 2), mapXarr.dtype)
        XY_mapPoints[:, 0] = mapXarr.flat
        XY_mapPoints[:, 1] = mapYarr.flat
213

Daniel Scheffler's avatar
Daniel Scheffler committed
214
215
        assert XY_points.shape == XY_mapPoints.shape

216
        return XY_points, XY_mapPoints
217

218
    def _exclude_bad_XYpos(self, GDF):
Daniel Scheffler's avatar
Daniel Scheffler committed
219
        """Exclude all points outside of the image overlap area and where the bad data mask is True (if given).
220
221
222
223

        :param GDF:     <geopandas.GeoDataFrame> must include the columns 'X_UTM' and 'Y_UTM'
        :return:
        """
224
225
        from skimage.measure import points_in_poly  # import here to avoid static TLS ImportError

226
227
228
229
        # exclude all points outside of overlap area
        inliers = points_in_poly(self.XY_mapPoints,
                                 np.swapaxes(np.array(self.COREG_obj.overlap_poly.exterior.coords.xy), 0, 1))
        GDF = GDF[inliers].copy()
230
        # GDF = GDF[GDF['geometry'].within(self.COREG_obj.overlap_poly.simplify(tolerance=15))] # works but much slower
231

232
        assert not GDF.empty, 'No coregistration point could be placed within the overlap area. Check your input data!'
233

234
        # exclude all points where bad data mask is True (e.g. points on clouds etc.)
235
236
237
238
239
240
        orig_len_GDF = len(GDF)  # length of GDF after dropping all points outside the overlap polygon
        mapXY = np.array(GDF.loc[:, ['X_UTM', 'Y_UTM']])
        GDF['REF_BADDATA'] = self.COREG_obj.ref.mask_baddata.read_pointData(mapXY) \
            if self.COREG_obj.ref.mask_baddata is not None else False
        GDF['TGT_BADDATA'] = self.COREG_obj.shift.mask_baddata.read_pointData(mapXY) \
            if self.COREG_obj.shift.mask_baddata is not None else False
Daniel Scheffler's avatar
Daniel Scheffler committed
241
        GDF = GDF[(~GDF['REF_BADDATA']) & (~GDF['TGT_BADDATA'])]
242
        if self.COREG_obj.ref.mask_baddata is not None or self.COREG_obj.shift.mask_baddata is not None:
Daniel Scheffler's avatar
Daniel Scheffler committed
243
            if not self.q:
244
245
246
247
248
249
                if not GDF.empty:
                    print('With respect to the provided bad data mask(s) %s points of initially %s have been excluded.'
                          % (orig_len_GDF - len(GDF), orig_len_GDF))
                else:
                    warnings.warn('With respect to the provided bad data mask(s) no coregistration point could be '
                                  'placed within an image area usable for coregistration.')
250
251
252

        return GDF

253
254
    @staticmethod
    def _get_spatial_shifts(coreg_kwargs):
Daniel Scheffler's avatar
Daniel Scheffler committed
255
        # unpack
256
        pointID = coreg_kwargs['pointID']
257
258
        fftw_works = coreg_kwargs['fftw_works']
        del coreg_kwargs['pointID'], coreg_kwargs['fftw_works']
259

Daniel Scheffler's avatar
Daniel Scheffler committed
260
        # assertions
261
        assert global_shared_imref is not None
262
        assert global_shared_im2shift is not None
Daniel Scheffler's avatar
Daniel Scheffler committed
263
264

        # run CoReg
265
        CR = COREG(global_shared_imref, global_shared_im2shift, CPUs=1, **coreg_kwargs)
266
        CR.fftw_works = fftw_works
267
        CR.calculate_spatial_shifts()
Daniel Scheffler's avatar
Daniel Scheffler committed
268
269

        # fetch results
270
        last_err = CR.tracked_errors[-1] if CR.tracked_errors else None
271
        win_sz_y, win_sz_x = CR.matchBox.imDimsYX if CR.matchBox else (None, None)
272
273
274
        CR_res = [win_sz_x, win_sz_y, CR.x_shift_px, CR.y_shift_px, CR.x_shift_map, CR.y_shift_map,
                  CR.vec_length_map, CR.vec_angle_deg, CR.ssim_orig, CR.ssim_deshifted, CR.ssim_improved,
                  CR.shift_reliability, last_err]
275

276
        return [pointID] + CR_res
277

278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
    def _get_coreg_kwargs(self, pID, wp):
        return dict(
            pointID=pID,
            fftw_works=self.COREG_obj.fftw_works,
            wp=wp,
            ws=self.COREG_obj.win_size_XY,
            resamp_alg_calc=self.rspAlg_calc,
            footprint_poly_ref=self.COREG_obj.ref.poly,
            footprint_poly_tgt=self.COREG_obj.shift.poly,
            r_b4match=self.ref.band4match + 1,  # band4match is internally saved as index, starting from 0
            s_b4match=self.shift.band4match + 1,  # band4match is internally saved as index, starting from 0
            max_iter=self.COREG_obj.max_iter,
            max_shift=self.COREG_obj.max_shift,
            nodata=(self.COREG_obj.ref.nodata, self.COREG_obj.shift.nodata),
            force_quadratic_win=self.COREG_obj.force_quadratic_win,
            binary_ws=self.COREG_obj.bin_ws,
            v=False,  # otherwise this would lead to massive console output
            q=True,  # otherwise this would lead to massive console output
            ignore_errors=True
        )

299
    def get_CoRegPoints_table(self):
300
301
        assert self.XY_points is not None and self.XY_mapPoints is not None

302
303
304
305
        # create a dataframe containing 'geometry','POINT_ID','X_IM','Y_IM','X_UTM','Y_UTM'
        # (convert imCoords to mapCoords
        XYarr2PointGeom = np.vectorize(lambda X, Y: Point(X, Y), otypes=[Point])
        geomPoints = np.array(XYarr2PointGeom(self.XY_mapPoints[:, 0], self.XY_mapPoints[:, 1]))
306

307
308
309
        if isLocal(self.COREG_obj.shift.prj):
            crs = None
        elif isProjectedOrGeographic(self.COREG_obj.shift.prj) == 'geographic':
310
            crs = dict(ellps='WGS84', datum='WGS84', proj='longlat')
311
        elif isProjectedOrGeographic(self.COREG_obj.shift.prj) == 'projected':
312
            UTMzone = abs(get_UTMzone(prj=self.COREG_obj.shift.prj))
313
314
315
316
            south = get_UTMzone(prj=self.COREG_obj.shift.prj) < 0
            crs = dict(ellps='WGS84', datum='WGS84', proj='utm', zone=UTMzone, south=south, units='m', no_defs=True)
            if not south:
                del crs['south']
317
318
319
        else:
            crs = None

320
321
322
        GDF = GeoDataFrame(index=range(len(geomPoints)), crs=crs,
                           columns=['geometry', 'POINT_ID', 'X_IM', 'Y_IM', 'X_UTM', 'Y_UTM'])
        GDF['geometry'] = geomPoints
Daniel Scheffler's avatar
Daniel Scheffler committed
323
        GDF['POINT_ID'] = range(len(geomPoints))
324
325
        GDF.loc[:, ['X_IM', 'Y_IM']] = self.XY_points
        GDF.loc[:, ['X_UTM', 'Y_UTM']] = self.XY_mapPoints
326

327
328
        # exclude offsite points and points on bad data mask
        GDF = self._exclude_bad_XYpos(GDF)
329
330
331
        if GDF.empty:
            self.CoRegPoints_table = GDF
            return self.CoRegPoints_table
332

333
        # choose a random subset of points if a maximum number has been given
334
        if self.max_points and len(GDF) > self.max_points:
335
            GDF = GDF.sample(self.max_points).copy()
336

337
        # equalize pixel grids in order to save warping time
338
339
340
341
        if len(GDF) > 100:
            # NOTE: actually grid res should be also changed here because self.shift.xgsd changes and grid res is
            # connected to that
            self.COREG_obj.equalize_pixGrids()
342

343
        # validate reference and target image inputs
344
        assert self.ref.footprint_poly  # this also checks for mask_nodata and nodata value
345
        assert self.shift.footprint_poly
346
347
348

        # ensure the input arrays for CoReg are in memory -> otherwise the code will get stuck in multiprocessing if
        # neighboured matching windows overlap during reading from disk!!
349
350
        self.ref.cache_array_subset(
            [self.COREG_obj.ref.band4match])  # only sets geoArr._arr_cache; does not change number of bands
Daniel Scheffler's avatar
Daniel Scheffler committed
351
352
        self.shift.cache_array_subset([self.COREG_obj.shift.band4match])

353
        # get all variations of kwargs for coregistration
354
        list_coreg_kwargs = (self._get_coreg_kwargs(i, self.XY_mapPoints[i]) for i in GDF.index)  # generator
355
356

        # run co-registration for whole grid
357
        if self.CPUs is None or self.CPUs > 1:
Daniel Scheffler's avatar
CoReg:    
Daniel Scheffler committed
358
            if not self.q:
359
                cpus = self.CPUs if self.CPUs is not None else multiprocessing.cpu_count()
360
                print("Calculating tie point grid (%s points) using %s CPU cores..." % (len(GDF), cpus))
361

362
            with multiprocessing.Pool(self.CPUs, initializer=mp_initializer, initargs=(self.ref, self.shift)) as pool:
363
364
365
366
                if self.q or not self.progress:
                    results = pool.map(self._get_spatial_shifts, list_coreg_kwargs)
                else:
                    results = pool.map_async(self._get_spatial_shifts, list_coreg_kwargs, chunksize=1)
367
                    bar = ProgressBar(prefix='\tprogress:')
368
369
                    while True:
                        time.sleep(.1)
370
371
                        # this does not really represent the remaining tasks but the remaining chunks
                        # -> thus chunksize=1
Daniel Scheffler's avatar
Fix.    
Daniel Scheffler committed
372
                        # noinspection PyProtectedMember
373
                        numberDone = len(GDF) - results._number_left
374
                        if self.progress:
375
                            bar.print_progress(percent=numberDone / len(GDF) * 100)
376
                        if results.ready():
377
378
379
                            # <= this is the line where multiprocessing can freeze if an exception appears within
                            # COREG ans is not raised
                            results = results.get()
380
                            break
Daniel Scheffler's avatar
Daniel Scheffler committed
381

382
        else:
383
384
385
386
387
            # declare global variables needed for self._get_spatial_shifts()
            global global_shared_imref, global_shared_im2shift
            global_shared_imref = self.ref
            global_shared_im2shift = self.shift

Daniel Scheffler's avatar
CoReg:    
Daniel Scheffler committed
388
            if not self.q:
389
                print("Calculating tie point grid (%s points) on 1 CPU core..." % len(GDF))
390
391
392
            results = np.empty((len(geomPoints), 14), np.object)
            bar = ProgressBar(prefix='\tprogress:')
            for i, coreg_kwargs in enumerate(list_coreg_kwargs):
393
                if self.progress:
394
395
                    bar.print_progress((i + 1) / len(GDF) * 100)
                results[i, :] = self._get_spatial_shifts(coreg_kwargs)
396

397
        # merge results with GDF
398
399
400
401
402
403
        # NOTE: We use a pandas.DataFrame here because the geometry column is missing.
        #       GDF.astype(...) fails with geopandas>0.6.0 if the geometry columns is missing.
        records = DataFrame(results,
                            columns=['POINT_ID', 'X_WIN_SIZE', 'Y_WIN_SIZE', 'X_SHIFT_PX', 'Y_SHIFT_PX', 'X_SHIFT_M',
                                     'Y_SHIFT_M', 'ABS_SHIFT', 'ANGLE', 'SSIM_BEFORE', 'SSIM_AFTER',
                                     'SSIM_IMPROVED', 'RELIABILITY', 'LAST_ERR'])
404

Daniel Scheffler's avatar
Daniel Scheffler committed
405
406
        # merge DataFrames (dtype must be equal to records.dtypes; We need np.object due to None values)
        GDF = GDF.astype(np.object).merge(records.astype(np.object), on='POINT_ID', how="inner")
407
408
        GDF = GDF.replace([np.nan, None], int(self.outFillVal))  # fillna fails with geopandas==0.6.0
        GDF.crs = crs  # gets lost when using GDF.astype(np.object), so we have to reassign that
409

410
411
412
        if not self.q:
            print("Found %s matches." % len(GDF[GDF.LAST_ERR == int(self.outFillVal)]))

413
        # filter tie points according to given filter level
414
        if self.tieP_filter_level > 0:
415
416
            if not self.q:
                print('Performing validity checks...')
417
            TPR = Tie_Point_Refiner(GDF[GDF.ABS_SHIFT != self.outFillVal], **self.outlDetect_settings)
418
            GDF_filt, new_columns = TPR.run_filtering(level=self.tieP_filter_level)
419
            GDF = GDF.merge(GDF_filt[['POINT_ID'] + new_columns], on='POINT_ID', how="outer")
420

421
        GDF = GDF.replace([np.nan, None], int(self.outFillVal))  # fillna fails with geopandas==0.6.0
422

423
        self.CoRegPoints_table = GDF
424

425
426
427
428
        if not self.q:
            if GDF.empty:
                warnings.warn('No valid GCPs could by identified.')
            else:
429
430
                if self.tieP_filter_level > 0:
                    print("%d valid tie points remain after filtering." % len(GDF[GDF.OUTLIER.__eq__(False)]))
431

432
433
        return self.CoRegPoints_table

434
435
    def calc_rmse(self, include_outliers=False):
        # type: (bool) -> float
Daniel Scheffler's avatar
Daniel Scheffler committed
436
        """Calculate root mean square error of absolute shifts from the tie point grid.
437

Daniel Scheffler's avatar
Daniel Scheffler committed
438
        :param include_outliers:    whether to include tie points that have been marked as false-positives (if present)
439
440
        """
        tbl = self.CoRegPoints_table
441
        tbl = tbl if include_outliers else tbl[tbl['OUTLIER'] == 0].copy() if 'OUTLIER' in tbl.columns else tbl
442
443
444
445
446
447

        shifts = np.array(tbl['ABS_SHIFT'])
        shifts_sq = [i * i for i in shifts if i != self.outFillVal]

        return np.sqrt(sum(shifts_sq) / len(shifts_sq))

448
449
    def calc_overall_ssim(self, include_outliers=False, after_correction=True):
        # type: (bool, bool) -> float
Daniel Scheffler's avatar
Daniel Scheffler committed
450
        """Calculate the median value of all SSIM values contained in tie point grid.
451
452

        :param include_outliers:    whether to include tie points that have been marked as false-positives
453
        :param after_correction:    whether to compute median SSIM before correction or after
454
455
        """
        tbl = self.CoRegPoints_table
456
        tbl = tbl if include_outliers else tbl[tbl['OUTLIER'] == 0].copy()
457

458
459
        ssim_col = np.array(tbl['SSIM_AFTER' if after_correction else 'SSIM_BEFORE'])
        ssim_col = [i * i for i in ssim_col if i != self.outFillVal]
460

461
        return float(np.median(ssim_col))
462
463
464

    def plot_shift_distribution(self, include_outliers=True, unit='m', interactive=False, figsize=None, xlim=None,
                                ylim=None, fontsize=12, title='shift distribution'):
Daniel Scheffler's avatar
Daniel Scheffler committed
465
        # type: (bool, str, bool, tuple, list, list, int, str) -> tuple
Daniel Scheffler's avatar
Daniel Scheffler committed
466
        """Create a 2D scatterplot containing the distribution of calculated X/Y-shifts.
467
468
469
470
471
472
473
474
475
476

        :param include_outliers:    whether to include tie points that have been marked as false-positives
        :param unit:                'm' for meters or 'px' for pixels (default: 'm')
        :param interactive:         interactive mode uses plotly for visualization
        :param figsize:             (xdim, ydim)
        :param xlim:                [xmin, xmax]
        :param ylim:                [ymin, ymax]
        :param fontsize:            size of all used fonts
        :param title:               the title to be plotted above the figure
        """
477
478
        from matplotlib import pyplot as plt

479
480
        if unit not in ['m', 'px']:
            raise ValueError("Parameter 'unit' must have the value 'm' (meters) or 'px' (pixels)! Got %s." % unit)
481

482
483
484
        if self.CoRegPoints_table.empty:
            raise RuntimeError('Shift distribution cannot be plotted because no tie points were found at all.')

485
        tbl = self.CoRegPoints_table
Daniel Scheffler's avatar
Daniel Scheffler committed
486
        tbl = tbl[tbl['ABS_SHIFT'] != self.outFillVal]
487
        tbl_il = tbl[tbl['OUTLIER'] == 0].copy() if 'OUTLIER' in tbl.columns else tbl
Daniel Scheffler's avatar
Daniel Scheffler committed
488
        tbl_ol = tbl[tbl['OUTLIER']].copy() if 'OUTLIER' in tbl.columns else None
489
490
        x_attr = 'X_SHIFT_M' if unit == 'm' else 'X_SHIFT_PX'
        y_attr = 'Y_SHIFT_M' if unit == 'm' else 'Y_SHIFT_PX'
491
492
        rmse = self.calc_rmse(include_outliers=False)  # always exclude outliers when calculating RMSE
        figsize = figsize if figsize else (10, 10)
493
494
495
496

        if interactive:
            from plotly.offline import iplot, init_notebook_mode
            import plotly.graph_objs as go
Daniel Scheffler's avatar
Daniel Scheffler committed
497
            # FIXME outliers are not plotted
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518

            init_notebook_mode(connected=True)

            # Create a trace
            trace = go.Scatter(
                x=tbl_il[x_attr],
                y=tbl_il[y_attr],
                mode='markers'
            )

            data = [trace]

            # Plot and embed in ipython notebook!
            iplot(data, filename='basic-scatter')

            return None, None

        else:
            fig = plt.figure(figsize=figsize)
            ax = fig.add_subplot(111)

Daniel Scheffler's avatar
Daniel Scheffler committed
519
            if include_outliers and 'OUTLIER' in tbl.columns:
Daniel Scheffler's avatar
Daniel Scheffler committed
520
                ax.scatter(tbl_ol[x_attr], tbl_ol[y_attr], marker='+', c='r', label='false-positives')
521
522
523
524
525
526
527
528
529
530
531
532
            ax.scatter(tbl_il[x_attr], tbl_il[y_attr], marker='+', c='g', label='valid tie points')

            # set axis limits
            if not xlim:
                xmax = np.abs(tbl_il[x_attr]).max()
                xlim = [-xmax, xmax]
            if not ylim:
                ymax = np.abs(tbl_il[y_attr]).max()
                ylim = [-ymax, ymax]
            ax.set_xlim(xlim)
            ax.set_ylim(ylim)

Daniel Scheffler's avatar
Daniel Scheffler committed
533
            # add text box containing RMSE of plotted shifts
534
            xlim, ylim = ax.get_xlim(), ax.get_ylim()
535
536
            plt.text(xlim[1] - (xlim[1] / 20), -ylim[1] + (ylim[1] / 20),
                     'RMSE:  %s m / %s px' % (np.round(rmse, 2), np.round(rmse / self.shift.xgsd, 2)),
537
                     ha='right', va='bottom', fontsize=fontsize, bbox=dict(facecolor='w', pad=None, alpha=0.8))
538

Daniel Scheffler's avatar
Daniel Scheffler committed
539
            # add grid and increase linewidth of middle line
540
541
542
            plt.grid()
            xgl = ax.get_xgridlines()
            middle_xgl = xgl[int(np.median(np.array(range(len(xgl)))))]
Daniel Scheffler's avatar
Daniel Scheffler committed
543
            middle_xgl.set_linewidth(2)
544
545
546
            middle_xgl.set_linestyle('-')
            ygl = ax.get_ygridlines()
            middle_ygl = ygl[int(np.median(np.array(range(len(ygl)))))]
Daniel Scheffler's avatar
Daniel Scheffler committed
547
            middle_ygl.set_linewidth(2)
548
549
            middle_ygl.set_linestyle('-')

Daniel Scheffler's avatar
Daniel Scheffler committed
550
551
            # set title and adjust tick labels
            ax.set_title(title, fontsize=fontsize)
552
553
            [tick.label.set_fontsize(fontsize) for tick in ax.xaxis.get_major_ticks()]
            [tick.label.set_fontsize(fontsize) for tick in ax.yaxis.get_major_ticks()]
Daniel Scheffler's avatar
Daniel Scheffler committed
554
555
            plt.xlabel('x-shift [%s]' % 'meters' if unit == 'm' else 'pixels', fontsize=fontsize)
            plt.ylabel('y-shift [%s]' % 'meters' if unit == 'm' else 'pixels', fontsize=fontsize)
556

557
558
            # add legend with labels in the right order
            handles, labels = ax.get_legend_handles_labels()
Daniel Scheffler's avatar
Daniel Scheffler committed
559
560
            leg = plt.legend(reversed(handles), reversed(labels), fontsize=fontsize, loc='upper right', scatterpoints=3)
            leg.get_frame().set_edgecolor('black')
561

562
563
564
565
            plt.show()

            return fig, ax

566
    def dump_CoRegPoints_table(self, path_out=None):
567
        path_out = path_out if path_out else \
568
569
570
571
572
            get_generic_outpath(dir_out=self.dir_out,
                                fName_out="CoRegPoints_table_grid%s_ws(%s_%s)__T_%s__R_%s.pkl"
                                          % (self.grid_res, self.COREG_obj.win_size_XY[0],
                                             self.COREG_obj.win_size_XY[1], self.shift.basename,
                                             self.ref.basename))
573
574
575
        if not self.q:
            print('Writing %s ...' % path_out)
        self.CoRegPoints_table.to_pickle(path_out)
576

577
    def to_GCPList(self):
Daniel Scheffler's avatar
Daniel Scheffler committed
578
        # get copy of tie points grid without no data
Daniel Scheffler's avatar
Daniel Scheffler committed
579
580
581
582
583
        try:
            GDF = self.CoRegPoints_table.loc[self.CoRegPoints_table.ABS_SHIFT != self.outFillVal, :].copy()
        except AttributeError:
            # self.CoRegPoints_table has no attribute 'ABS_SHIFT' because all points have been excluded
            return []
584

585
        if getattr(GDF, 'empty'):  # GDF.empty returns AttributeError
586
587
            return []
        else:
588
            # exclude all points flagged as outliers
589
            if 'OUTLIER' in GDF.columns:
590
                GDF = GDF[GDF.OUTLIER.__eq__(False)].copy()
591
592
            avail_TP = len(GDF)

593
594
595
596
            if not avail_TP:
                # no point passed all validity checks
                return []

597
            if avail_TP > 7000:
598
599
600
                GDF = GDF.sample(7000)
                warnings.warn('By far not more than 7000 tie points can be used for warping within a limited '
                              'computation time (due to a GDAL bottleneck). Thus these 7000 points are randomly chosen '
601
                              'out of the %s available tie points.' % avail_TP)
602

603
604
605
            # calculate GCPs
            GDF['X_UTM_new'] = GDF.X_UTM + GDF.X_SHIFT_M
            GDF['Y_UTM_new'] = GDF.Y_UTM + GDF.Y_SHIFT_M
606
607
608
609
610
611
            GDF['GCP'] = GDF.apply(lambda GDF_row: gdal.GCP(GDF_row.X_UTM_new,
                                                            GDF_row.Y_UTM_new,
                                                            0,
                                                            GDF_row.X_IM,
                                                            GDF_row.Y_IM),
                                   axis=1)
612
613
614
            self.GCPList = GDF.GCP.tolist()

            return self.GCPList
615

616
    def test_if_singleprocessing_equals_multiprocessing_result(self):
617
618
        # RANSAC filtering always produces different results because it includes random sampling
        self.tieP_filter_level = 1
619

Daniel Scheffler's avatar
Daniel Scheffler committed
620
        self.CPUs = None
621
        dataframe = self.get_CoRegPoints_table()
622
        mp_out = np.empty_like(dataframe.values)
623
        mp_out[:] = dataframe.values
Daniel Scheffler's avatar
Daniel Scheffler committed
624
        self.CPUs = 1
625
        dataframe = self.get_CoRegPoints_table()
626
        sp_out = np.empty_like(dataframe.values)
627
628
        sp_out[:] = dataframe.values

629
        return np.array_equal(sp_out, mp_out)
630

631
632
    def _get_line_by_PID(self, PID):
        return self.CoRegPoints_table.loc[PID, :]
633

634
    def _get_lines_by_PIDs(self, PIDs):
635
636
637
638
        assert isinstance(PIDs, list)
        lines = np.zeros((len(PIDs), self.CoRegPoints_table.shape[1]))
        for i, PID in enumerate(PIDs):
            lines[i, :] = self.CoRegPoints_table[self.CoRegPoints_table['POINT_ID'] == PID]
639
640
        return lines

641
    def to_PointShapefile(self, path_out=None, skip_nodata=True, skip_nodata_col='ABS_SHIFT'):
642
        # type: (str, bool, str) -> None
Daniel Scheffler's avatar
Daniel Scheffler committed
643
644
645
        """Write the calculated tie points grid to a point shapefile (e.g., for visualization by a GIS software).

        NOTE: The shapefile uses Tie_Point_Grid.CoRegPoints_table as attribute table.
646
647
648

        :param path_out:        <str> the output path. If not given, it is automatically defined.
        :param skip_nodata:     <bool> whether to skip all points where no valid match could be found
649
        :param skip_nodata_col: <str> determines which column of Tie_Point_Grid.CoRegPoints_table is used to
650
651
                                identify points where no valid match could be found
        """
652
        GDF = self.CoRegPoints_table
653
654
655

        if skip_nodata:
            GDF2pass = GDF[GDF[skip_nodata_col] != self.outFillVal].copy()
656
657
        else:
            GDF2pass = GDF
658
            GDF2pass.LAST_ERR = GDF2pass.apply(lambda GDF_row: repr(GDF_row.LAST_ERR), axis=1)
Daniel Scheffler's avatar
Daniel Scheffler committed
659
660

        # replace boolean values (cannot be written)
661
662
        GDF2pass = GDF2pass.replace(False, 0).copy()  # replace booleans where column dtype is not np.bool but np.object
        GDF2pass = GDF2pass.replace(True, 1).copy()
Daniel Scheffler's avatar
Daniel Scheffler committed
663
664
665
        for col in GDF2pass.columns:
            if GDF2pass[col].dtype == np.bool:
                GDF2pass[col] = GDF2pass[col].astype(int)
666
667
668
669
670
671

        path_out = path_out if path_out else \
            get_generic_outpath(dir_out=os.path.join(self.dir_out, 'CoRegPoints'),
                                fName_out="CoRegPoints_grid%s_ws(%s_%s)__T_%s__R_%s.shp"
                                          % (self.grid_res, self.COREG_obj.win_size_XY[0],
                                             self.COREG_obj.win_size_XY[1], self.shift.basename, self.ref.basename))
Daniel Scheffler's avatar
CoReg:    
Daniel Scheffler committed
672
        if not self.q:
673
            print('Writing %s ...' % path_out)
674
675
        GDF2pass.to_file(path_out)

676
    def _to_PointShapefile(self, skip_nodata=True, skip_nodata_col='ABS_SHIFT'):  # pragma: no cover
677
678
679
680
        warnings.warn(DeprecationWarning(
            "'_tiepoints_grid_to_PointShapefile' is deprecated."  # TODO delete if other method validated
            " 'tiepoints_grid_to_PointShapefile' is much faster."))
        GDF = self.CoRegPoints_table
681
682
683
        GDF2pass = \
            GDF if not skip_nodata else \
            GDF[GDF[skip_nodata_col] != self.outFillVal]
684
        shapely_points = GDF2pass['geometry'].values.tolist()
685
686
687
        attr_dicts = [collections.OrderedDict(zip(GDF2pass.columns,
                                                  GDF2pass.loc[i].values))
                      for i in GDF2pass.index]
688

689
690
        fName_out = "CoRegPoints_grid%s_ws%s.shp" \
                    % (self.grid_res, self.COREG_obj.win_size_XY)
691
        path_out = os.path.join(self.dir_out, fName_out)
692
693
694
695
        write_shp(path_out,
                  shapely_points,
                  prj=self.COREG_obj.shift.prj,
                  attrDict=attr_dicts)
696

697
    def to_vectorfield(self, path_out=None, fmt=None, mode='md'):
Daniel Scheffler's avatar
Daniel Scheffler committed
698
        # type: (str, str, str) -> GeoArray
Daniel Scheffler's avatar
Daniel Scheffler committed
699
700
701
        """Save the calculated X-/Y-shifts to a 2-band raster file that can be used to visualize a vectorfield.

        NOTE: For example ArcGIS is able to visualize such 2-band raster files as a vectorfield.
702
703
704

        :param path_out:    <str> the output path. If not given, it is automatically defined.
        :param fmt:         <str> output raster format string
705
706
        :param mode:        <str> The mode how the output is written ('uv' or 'md'; default: 'md')
                                    'uv': outputs X-/Y shifts
707
708
                                    'md': outputs magnitude and direction
        """
709
710
        assert mode in ['uv', 'md'], "'mode' must be either 'uv' (outputs X-/Y shifts) or 'md' " \
                                     "(outputs magnitude and direction)'. Got %s." % mode
711
712
        attr_b1 = 'X_SHIFT_M' if mode == 'uv' else 'ABS_SHIFT'
        attr_b2 = 'Y_SHIFT_M' if mode == 'uv' else 'ANGLE'
713

714
715
716
        xshift_arr, gt, prj = points_to_raster(points=self.CoRegPoints_table['geometry'],
                                               values=self.CoRegPoints_table[attr_b1],
                                               tgt_res=self.shift.xgsd * self.grid_res,
717
                                               prj=self.CoRegPoints_table.crs.to_wkt(),
718
                                               fillVal=self.outFillVal)
719

720
721
722
        yshift_arr, gt, prj = points_to_raster(points=self.CoRegPoints_table['geometry'],
                                               values=self.CoRegPoints_table[attr_b2],
                                               tgt_res=self.shift.xgsd * self.grid_res,
723
                                               prj=self.CoRegPoints_table.crs.to_wkt(),
724
                                               fillVal=self.outFillVal)
725
726
727

        out_GA = GeoArray(np.dstack([xshift_arr, yshift_arr]), gt, prj, nodata=self.outFillVal)

728
729
730
731
732
        path_out = path_out if path_out else \
            get_generic_outpath(dir_out=os.path.join(self.dir_out, 'CoRegPoints'),
                                fName_out="CoRegVectorfield%s_ws(%s_%s)__T_%s__R_%s.tif"
                                          % (self.grid_res, self.COREG_obj.win_size_XY[0],
                                             self.COREG_obj.win_size_XY[1], self.shift.basename, self.ref.basename))
733
734
735
736
737

        out_GA.save(path_out, fmt=fmt if fmt else 'Gtiff')

        return out_GA

738
    def to_Raster_using_Kriging(self, attrName, skip_nodata=1, skip_nodata_col='ABS_SHIFT', outGridRes=None,
739
                                fName_out=None, tilepos=None, tilesize=500, mp=None):
740

741
        mp = False if self.CPUs == 1 else True
742
743
        self._Kriging_sp(attrName, skip_nodata=skip_nodata, skip_nodata_col=skip_nodata_col,
                         outGridRes=outGridRes, fName_out=fName_out, tilepos=tilepos)
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765

        # if mp:
        #     tilepositions = UTL.get_image_tileborders([tilesize,tilesize],self.tgt_shape)
        #     args_kwargs_dicts=[]
        #     for tp in tilepositions:
        #         kwargs_dict = {'skip_nodata':skip_nodata,'skip_nodata_col':skip_nodata_col,'outGridRes':outGridRes,
        #                        'fName_out':fName_out,'tilepos':tp}
        #         args_kwargs_dicts.append({'args':[attrName],'kwargs':kwargs_dict})
        #     # self.kriged=[]
        #     # for i in args_kwargs_dicts:
        #     #     res = self.Kriging_mp(i)
        #     #     self.kriged.append(res)
        #     #     print(res)
        #
        #     with multiprocessing.Pool() as pool:
        #        self.kriged = pool.map(self.Kriging_mp,args_kwargs_dicts)
        # else:
        #     self.Kriging_sp(attrName,skip_nodata=skip_nodata,skip_nodata_col=skip_nodata_col,
        #                     outGridRes=outGridRes,fName_out=fName_out,tilepos=tilepos)
        res = self.kriged if mp else None
        return res

766
767
    def _Kriging_sp(self, attrName, skip_nodata=1, skip_nodata_col='ABS_SHIFT', outGridRes=None,
                    fName_out=None, tilepos=None):
768
769
        GDF = self.CoRegPoints_table
        GDF2pass = GDF if not skip_nodata else GDF[GDF[skip_nodata_col] != self.outFillVal]
770

771
        X_coords, Y_coords, ABS_SHIFT = GDF2pass['X_UTM'], GDF2pass['Y_UTM'], GDF2pass[attrName]
772

773
        xmin, ymin, xmax, ymax = GDF2pass.total_bounds
774

775
776
        grid_res = outGridRes if outGridRes else int(min(xmax - xmin, ymax - ymin) / 250)
        grid_x, grid_y = np.arange(xmin, xmax + grid_res, grid_res), np.arange(ymax, ymin - grid_res, -grid_res)
777
778
779

        # Reference: P.K. Kitanidis, Introduction to Geostatistcs: Applications in Hydrogeology,
        #            (Cambridge University Press, 1997) 272 p.
780
        from pykrige.ok import OrdinaryKriging
781
782
        OK = OrdinaryKriging(X_coords, Y_coords, ABS_SHIFT, variogram_model='spherical', verbose=False)
        zvalues, sigmasq = OK.execute('grid', grid_x, grid_y, backend='C', n_closest_points=12)
783

784
        if self.CPUs is None or self.CPUs > 1:
785
            fName_out = fName_out if fName_out else \
786
                "Kriging__%s__grid%s_ws%s_%s.tif" % (attrName, self.grid_res, self.COREG_obj.win_size_XY, tilepos)
787
788
        else:
            fName_out = fName_out if fName_out else \
789
790
                "Kriging__%s__grid%s_ws%s.tif" % (attrName, self.grid_res, self.COREG_obj.win_size_XY)
        path_out = get_generic_outpath(dir_out=self.dir_out, fName_out=fName_out)
791
        # add a half pixel grid points are centered on the output pixels
792
793
794
795
        xmin = xmin - grid_res / 2
        # ymin = ymin - grid_res / 2
        # xmax = xmax + grid_res / 2
        ymax = ymax + grid_res / 2
796
797
798
799

        GeoArray(zvalues,
                 geotransform=(xmin, grid_res, 0, ymax, 0, -grid_res),
                 projection=self.COREG_obj.shift.prj).save(path_out)
800
801
802

        return zvalues

803
    def _Kriging_mp(self, args_kwargs_dict):
804
805
        args = args_kwargs_dict.get('args', [])
        kwargs = args_kwargs_dict.get('kwargs', [])
806

807
        return self._Kriging_sp(*args, **kwargs)
808
809


810
class Tie_Point_Refiner(object):
Daniel Scheffler's avatar
Daniel Scheffler committed
811
812
    """A class for performing outlier detection."""

Daniel Scheffler's avatar
Daniel Scheffler committed
813
    def __init__(self, GDF, min_reliability=60, rs_max_outlier=10, rs_tolerance=2.5, rs_max_iter=15,
814
                 rs_exclude_previous_outliers=True, rs_timeout=20, q=False):
Daniel Scheffler's avatar
Daniel Scheffler committed
815
        """Get an instance of Tie_Point_Refiner.
Daniel Scheffler's avatar
Daniel Scheffler committed
816

817
818
        :param GDF:                             GeoDataFrame like TiePointGrid.CoRegPoints_table containing all tie
                                                points to be filtered and the corresponding metadata
Daniel Scheffler's avatar
Daniel Scheffler committed
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
        :param min_reliability:                 <float, int> minimum threshold for previously computed tie X/Y shift
                                                reliability (default: 60%)
        :param rs_max_outlier:                  <float, int> RANSAC: maximum percentage of outliers to be detected
                                                (default: 10%)
        :param rs_tolerance:                    <float, int> RANSAC: percentage tolerance for max_outlier_percentage
                                                (default: 2.5%)
        :param rs_max_iter:                     <int> RANSAC: maximum iterations for finding the best RANSAC threshold
                                                (default: 15)
        :param rs_exclude_previous_outliers:    <bool> RANSAC: whether to exclude points that have been flagged as
                                                outlier by earlier filtering (default:True)
        :param rs_timeout:                      <float, int> RANSAC: timeout for iteration loop in seconds (default: 20)

        :param q:
        """
        self.GDF = GDF.copy()
        self.min_reliability = min_reliability
        self.rs_max_outlier_percentage = rs_max_outlier
        self.rs_tolerance = rs_tolerance
        self.rs_max_iter = rs_max_iter
        self.rs_exclude_previous_outliers = rs_exclude_previous_outliers
        self.rs_timeout = rs_timeout
        self.q = q
        self.new_cols = []
842
843
        self.ransac_model_robust = None

844
    def run_filtering(self, level=3):
845
846
        """Filter tie points used for shift correction.

847
        :param level:   tie point filter level (default: 3).
848
849
850
851
852
853
854
855
                        NOTE: lower levels are also included if a higher level is chosen
                            - Level 0: no tie point filtering
                            - Level 1: Reliablity filtering - filter all tie points out that have a low
                                reliability according to internal tests
                            - Level 2: SSIM filtering - filters all tie points out where shift
                                correction does not increase image similarity within matching window
                                (measured by mean structural similarity index)
                            - Level 3: RANSAC outlier detection
Daniel Scheffler's avatar
Daniel Scheffler committed
856
857
858

        :return:
        """
859
860
        # TODO catch empty GDF

861
        # RELIABILITY filtering
862
        if level > 0:
863
            marked_recs = self._reliability_thresholding()  # type: Series
864
865
            self.GDF['L1_OUTLIER'] = marked_recs
            self.new_cols.append('L1_OUTLIER')
Daniel Scheffler's avatar
Daniel Scheffler committed
866

867
868
869
870
871
872
873
874
875
876
            n_flagged = len(marked_recs[marked_recs])
            perc40 = np.percentile(self.GDF.RELIABILITY, 40)

            if n_flagged / len(self.GDF) > .7:
                warnings.warn(r"More than 70%% of the found tie points have a reliability lower than %s%% and are "
                              r"therefore marked as false-positives. Consider relaxing the minimum reliability "
                              r"(parameter 'min_reliability') to avoid that. For example min_reliability=%d would only "
                              r"flag 40%% of the tie points in case of your input data."
                              % (self.min_reliability, perc40))

877
            if not self.q:
878
                print('%s tie points flagged by level 1 filtering (reliability).'
879
                      % n_flagged)
Daniel Scheffler's avatar
Daniel Scheffler committed
880

881
        # SSIM filtering
882
        if level > 1:
883
884
            marked_recs = self._SSIM_filtering()
            self.GDF['L2_OUTLIER'] = marked_recs  # type: Series
885
            self.new_cols.append('L2_OUTLIER')
Daniel Scheffler's avatar
Daniel Scheffler committed
886

887
            if not self.q:
Daniel Scheffler's avatar
Daniel Scheffler committed
888
                print('%s tie points flagged by level 2 filtering (SSIM).' % (len(marked_recs[marked_recs])))
Daniel Scheffler's avatar
Daniel Scheffler committed
889

890
        # RANSAC filtering
891
        if level > 2:
Daniel Scheffler's avatar
Daniel Scheffler committed
892
            # exclude previous outliers
Daniel Scheffler's avatar
Daniel Scheffler committed
893
            ransacInGDF = self.GDF[~self.GDF[self.new_cols].any(axis=1)].copy() \
894
                    if self.rs_exclude_previous_outliers else self.GDF
Daniel Scheffler's avatar
Daniel Scheffler committed
895

Daniel Scheffler's avatar