hylite.hyimage

   1"""
   2Store and manipulate hyperspectral image data.
   3"""
   4
   5import os
   6import numpy as np
   7import hylite
   8from hylite.hydata import HyData
   9from hylite.hylibrary import HyLibrary
  10from hylite._deps import require
  11
  12
  13
  14class HyImage( HyData ):
  15    """
  16    A class for hyperspectral image data. These can be individual scenes or hyperspectral orthoimages.
  17    """
  18
  19    def __init__(self, data, **kwds):
  20        """
  21        Args:
  22            data (ndarray): a numpy array such that data[x][y][band] gives each pixel value.
  23            **kwds:
  24                wav = A numpy array containing band wavelengths for this image.
  25                affine = an affine transform of the format returned by GDAL.GetGeoTransform().
  26                projection = string defining the project. Default is None.
  27                sensor = sensor name. Default is "unknown".
  28                header = path to associated header file. Default is None.
  29        """
  30
  31        #call constructor for HyData
  32        super().__init__(data, **kwds)
  33
  34        # special case - if dataset only has oneband, slice it so it still has
  35        # the format data[x,y,b].
  36        if not self.data is None:
  37            if len(self.data.shape) == 1:
  38                self.data = self.data[None, None, :] # single pixel image
  39            if len(self.data.shape) == 2:
  40                self.data = self.data[:, :, None] # single band iamge
  41
  42        #load any additional project information (specific to images)
  43        self.set_projection(kwds.get("projection",None))
  44        self.affine = np.array( kwds.get("affine",[0,1,0,0,0,1]) )
  45        self.header['affine'] = kwds.get("affine",[0,1,0,0,0,1]) # also store this here
  46
  47        # wavelengths
  48        if 'wav' in kwds:
  49            self.set_wavelengths(kwds['wav'])
  50
  51        #special header formatting
  52        self.header['file type'] = 'ENVI Standard'
  53
  54    def copy(self,data=True):
  55        """
  56        Make a deep copy of this image instance.
  57
  58        Args:
  59            data (bool): True if a copy of the data should be made, otherwise only copy header.
  60
  61        Returns:
  62            a new `hylite.hyimage.HyImage` instance.
  63        """
  64        if not data:
  65            return HyImage(None, header=self.header.copy(), projection=self.projection, affine=self.affine)
  66        else:
  67            return HyImage( self.data.copy(), header=self.header.copy(), projection=self.projection, affine=self.affine)
  68
  69    def T(self):
  70        """
  71        Return a transposed view of the data matrix (corresponding with the [y,x] indexing used by matplotlib, opencv etc.
  72        """
  73        return np.transpose(self.data, (1,0,2))
  74
  75    def xdim(self):
  76        """
  77        Return number of pixels in x (first dimension of data array)
  78        """
  79        return self.data.shape[0]
  80
  81    def ydim(self):
  82        """
  83        Return number of pixels in y (second dimension of data array)
  84        """
  85        return self.data.shape[1]
  86
  87    def aspx(self):
  88        """
  89        Return the aspect ratio of this image (width/height).
  90        """
  91        return self.ydim() / self.xdim()
  92    
  93    #####################################
  94    ## GEOREFERENCING METHODS
  95    #####################################
  96
  97    def get_extent(self):
  98        """
  99        Returns the width and height of this image in world coordinates.
 100
 101        Returns:
 102            tuple with (width, height).
 103        """
 104        return self.xdim * self.pixel_size[0], self.ydim * self.pixel_size[1]
 105
 106    def set_projection(self,proj):
 107        """
 108        Set this project to an existing osgeo.osr.SpatialReference or GDAL georeference string.
 109
 110        Args:
 111            proj (str, osgeo.osr.SpatialReference): the project to use as osgeo.osr.SpatialReference or GDAL georeference string.
 112        """
 113        if proj is None:
 114            self.projection = None
 115        else:
 116            try:
 117                from osgeo.osr import SpatialReference
 118            except:
 119                assert False, "Error - GDAL must be installed to work with spatial projections in hylite."
 120            if isinstance(proj, SpatialReference):
 121                self.projection = proj
 122            elif isinstance(proj, str):
 123                self.projection = SpatialReference(proj)
 124            else:
 125                print("Invalid project %s" % proj)
 126                raise
 127
 128    def set_projection_EPSG(self,EPSG):
 129        """
 130        Sets this image project using an EPSG code.
 131
 132        Args:
 133            EPSG (str): string EPSG code that can be passed to SpatialReference.SetFromUserInput(...).
 134        """
 135
 136        try:
 137            from osgeo.osr import SpatialReference
 138        except:
 139            assert False, "Error - GDAL must be installed to work with spatial projections in hylite."
 140
 141        self.projection = SpatialReference()
 142        self.projection.SetFromUserInput(EPSG)
 143
 144    def get_projection_EPSG(self):
 145        """
 146        Gets a string describing this projections EPSG code (if it is an EPSG project).
 147
 148        Returns:
 149            an EPSG code string of the format "EPSG:XXXX".
 150        """
 151        if self.projection is None:
 152            return None
 153        else:
 154            return "%s:%s" % (self.projection.GetAttrValue("AUTHORITY",0),self.projection.GetAttrValue("AUTHORITY",1))
 155
 156    def pix_to_world(self, px, py, proj=None):
 157        """
 158        Take pixel coordinates and return world coordinates
 159
 160        Args:
 161            px (int): the pixel x-coord.
 162            py (int): the pixel y-coord.
 163            proj (str, osr.SpatialReference): the coordinate system to use. Default (None) uses the same system as this image. Otherwise
 164                   an osr.SpatialReference can be passed (`hylite.hyimage.HyImage`.project), or an EPSG string (e.g. get_projection_EPSG(...)).
 165        Returns:
 166            the world coordinates in the coordinate system defined by get_projection_EPSG(...).
 167        """
 168
 169        try:
 170            from osgeo import osr
 171            import osgeo.gdal as gdal
 172            from osgeo import ogr
 173        except:
 174            assert False, "Error - GDAL must be installed to work with spatial projections in hylite."
 175
 176        # parse project
 177        if proj is None:
 178            proj = self.projection
 179        elif isinstance(proj, str) or isinstance(proj, int):
 180            epsg = proj
 181            if isinstance(epsg, str):
 182                try:
 183                    epsg = int(str.split(':')[1])
 184                except:
 185                    assert False, "Error - %s is an invalid EPSG code." % proj
 186            proj = osr.SpatialReference()
 187            proj.ImportFromEPSG(epsg)
 188
 189        # check we have all the required info
 190        assert isinstance(proj, osr.SpatialReference), "Error - invalid spatial reference %s" % proj
 191        assert (not self.affine is None) and (
 192            not self.projection is None), "Error - project information is undefined."
 193
 194        #project to world coordinates in this images project/world coords
 195        x,y = gdal.ApplyGeoTransform(self.affine, px, py)
 196
 197        #project to target coords (if different)
 198        if not proj.IsSameGeogCS(self.projection):
 199            P = ogr.Geometry(ogr.wkbPoint)
 200            if proj.EPSGTreatsAsNorthingEasting():
 201                P.AddPoint(x, y)
 202            else:
 203                P.AddPoint(y, x)
 204            P.AssignSpatialReference(self.projection)  # tell the point what coordinates it's in
 205            P.TransformTo(proj)  # reproject it to the out spatial reference
 206            x, y = P.GetX(), P.GetY()
 207
 208            #do we need to transpose?
 209            if proj.EPSGTreatsAsLatLong():
 210                x,y=y,x #we want lon,lat not lat,lon
 211        return x, y
 212
 213    def world_to_pix(self, x, y, proj = None):
 214        """
 215        Take world coordinates and return pixel coordinates
 216
 217        Args:
 218            x (float): the world x-coord.
 219            y (float): the world y-coord.
 220            proj (str, osr.SpatialReference): the coordinate system of the input coordinates. Default (None) uses the same system as this image. Otherwise
 221                   an osr.SpatialReference can be passed (`hylite.hyimage.HyImage`.project), or an EPSG string (e.g. get_projection_EPSG(...)).
 222
 223        Returns:
 224            the pixel coordinates based on the affine transform stored in self.affine.
 225        """
 226
 227        try:
 228            from osgeo import osr
 229            import osgeo.gdal as gdal
 230            from osgeo import ogr
 231        except:
 232            assert False, "Error - GDAL must be installed to work with spatial projections in hylite."
 233
 234        # parse project
 235        if proj is None:
 236            proj = self.projection
 237        elif isinstance(proj, str) or isinstance(proj, int):
 238            epsg = proj
 239            if isinstance(epsg, str):
 240                try:
 241                    epsg = int(str.split(':')[1])
 242                except:
 243                    assert False, "Error - %s is an invalid EPSG code." % proj
 244            proj = osr.SpatialReference()
 245            proj.ImportFromEPSG(epsg)
 246
 247
 248        # check we have all the required info
 249        assert isinstance(proj, osr.SpatialReference), "Error - invalid spatial reference %s" % proj
 250        assert (not self.affine is None) and (not self.projection is None), "Error - project information is undefined."
 251
 252        # project to this images CS (if different)
 253        if not proj.IsSameGeogCS(self.projection):
 254            P = ogr.Geometry(ogr.wkbPoint)
 255            if proj.EPSGTreatsAsNorthingEasting():
 256                P.AddPoint(x, y)
 257            else:
 258                P.AddPoint(y, x)
 259            P.AssignSpatialReference(proj)  # tell the point what coordinates it's in
 260            P.AddPoint(x, y)
 261            P.TransformTo(self.projection)  # reproject it to the out spatial reference
 262            x, y = P.GetX(), P.GetY()
 263            if self.projection.EPSGTreatsAsLatLong(): # do we need to transpose?
 264                x, y = y, x  # we want lon,lat not lat,lon
 265
 266        inv = gdal.InvGeoTransform(self.affine)
 267        assert not inv is None, "Error - could not invert affine transform?"
 268
 269        #apply
 270        return gdal.ApplyGeoTransform(inv, x, y)
 271
 272    def crop(self, xmin, xmax, ymin, ymax, bands=None):
 273        """
 274        Return a cropped copy of this image.
 275
 276        Args:
 277            xmin, xmax (int): pixel bounds in x (rows)
 278            ymin, ymax (int): pixel bounds in y (columns)
 279            bands (None, list, tuple): optional band indices or (min,max) range
 280
 281        Returns:
 282            `hylite.hyimage.HyImage`: cropped image with updated affine transform
 283        """
 284
 285        # ---- validate bounds ----
 286        xmin = int(max(0, xmin))
 287        ymin = int(max(0, ymin))
 288        xmax = int(min(self.xdim(), xmax))
 289        ymax = int(min(self.ydim(), ymax))
 290
 291        assert xmin < xmax and ymin < ymax, "Invalid crop extent."
 292
 293        # ---- crop data ----
 294        if bands is None:
 295            data = self.data[xmin:xmax, ymin:ymax, :].copy()
 296            wav = self.get_wavelengths()
 297        else: # band selection
 298            if isinstance(bands, tuple):
 299                b0 = self.get_band_index(bands[0])
 300                b1 = self.get_band_index(bands[1])
 301                data = self.data[xmin:xmax, ymin:ymax, b0:b1].copy()
 302                wav = self.get_wavelengths()[b0:b1]
 303            else:
 304                idx = [self.get_band_index(b) for b in bands]
 305                data = self.data[xmin:xmax, ymin:ymax, idx].copy()
 306                wav = self.get_wavelengths()[idx]
 307
 308        # ---- update affine transform ----
 309        if self.affine is not None:
 310            a = list(self.affine)
 311            new_affine = a.copy()
 312
 313            # shift origin to new top-left pixel
 314            new_affine[0] = a[0] + xmin*a[1] + ymin*a[2]
 315            new_affine[3] = a[3] + xmin*a[4] + ymin*a[5]
 316        else:
 317            new_affine = None
 318
 319        # ---- construct output image ----
 320        out = HyImage(
 321            data,
 322            header=self.header.copy(),
 323            projection=self.projection,
 324            affine=new_affine,
 325            wav=wav
 326        )
 327
 328        return out
 329
 330    def resize(self, newdims: tuple, interpolation: int = 1):
 331        """
 332        Resize this image with opencv and update affine transform accordingly.
 333
 334        Args:
 335            newdims (tuple): the new image dimensions (xdim, ydim)
 336            interpolation (int): opencv interpolation method. Default is cv2.INTER_LINEAR.
 337        """
 338        import cv2  # avoid import issues if opencv is missing
 339
 340        old_x, old_y = self.xdim(), self.ydim()
 341        new_x, new_y = int(newdims[0]), int(newdims[1])
 342
 343        assert new_x > 0 and new_y > 0, "Invalid resize dimensions."
 344
 345        # resize data (opencv uses width, height = y, x)
 346        self.data = cv2.resize(
 347            self.data,
 348            (new_y, new_x),
 349            interpolation=interpolation
 350        )
 351
 352        # update affine transform
 353        if self.affine is not None:
 354            a = list(self.affine)
 355
 356            sx = old_x / new_x
 357            sy = old_y / new_y
 358
 359            self.affine = [
 360                a[0],          # x origin unchanged
 361                a[1] * sx,     # pixel width
 362                a[2] * sy,     # row rotation
 363                a[3],          # y origin unchanged
 364                a[4] * sx,     # column rotation
 365                a[5] * sy      # pixel height
 366            ]
 367
 368    def tile(self, tile_size):
 369        """
 370        Break image into tiles of given size and return a list of `hylite.hyimage.HyImage` tiles.
 371        Each tile has an updated affine transform reflecting its position in the original image.
 372
 373        Args:
 374            tile_size (tuple): (tile_x, tile_y) in pixels
 375        Returns:
 376            list of `hylite.hyimage.HyImage`
 377        """
 378        tiles = []
 379        tx, ty = tile_size
 380        nx, ny = self.xdim(), self.ydim()
 381        for i in range(0, nx, tx):
 382            for j in range(0, ny, ty):
 383                data_tile = self.data[i:min(i+tx, nx), j:min(j+ty, ny), :].copy() # slice data
 384
 385                # compute new affine
 386                # affine = [x0, dx, rotx, y0, roty, dy]
 387                x0, dx, rotx, y0, roty, dy = self.affine
 388                new_x0 = x0 + i*dx + j*rotx
 389                new_y0 = y0 + i*roty + j*dy
 390                new_affine = [new_x0, dx, rotx, new_y0, roty, dy]
 391
 392                # create tile
 393                tile_img = HyImage(
 394                    data_tile,
 395                    affine=new_affine,
 396                    projection=self.projection,
 397                    wav=self.get_wavelengths(),
 398                    header=self.header.copy()
 399                )
 400                tile_img.header['xleft'] = i
 401                tile_img.header['ytop'] = j
 402                tiles.append(tile_img)
 403
 404        return tiles
 405
 406    @staticmethod
 407    def mosaic(
 408        tiles,
 409        blend="mean",
 410        resampling="nearest",
 411        out_affine=None,
 412        out_shape=None,
 413    ):
 414        """
 415        Mosaic georeferenced `hylite.hyimage.HyImage` tiles using GDAL. Note that this assumes all tiles are in the same coordinate system.
 416
 417        Args:
 418            tiles (list[`hylite.hyimage.HyImage`])
 419            blend (str): 'first', 'min', 'max', 'mean', 'median'
 420            resampling (str): 'nearest', 'bilinear', 'cubic'
 421            out_affine (list): optional 6-element affine to define output grid. If None, the affine of the first tile is used.
 422            out_shape (tuple): optional (xdim, ydim) shape of the output grid. If None, the extent of all tiles will be used.
 423        Returns:
 424            HyImage
 425        """
 426        import numpy as np
 427        from hylite.project.align import resample_raster
 428        from osgeo import gdal, osr
 429
 430        assert len(tiles) > 0
 431        assert blend in ("first", "min", "max", "mean", "median")
 432
 433        # compute bounds in world coordinates
 434        # N.B. THIS ASSUMES ALL DATA ARE IN THE SAME CRS
 435        points = []
 436        for t in tiles:
 437            points.append( t.pix_to_world(0,0) )
 438            points.append( t.pix_to_world(t.xdim()+1,t.ydim()+1) )
 439        min_x, min_y = np.min(points, axis=0)
 440        max_x, max_y = np.max(points, axis=0)
 441        if out_shape is None:
 442            out_shape = (np.array(tiles[0].world_to_pix(max_x, max_y)) + np.array(tiles[0].world_to_pix(min_x, min_y))).round().astype(int)
 443        if out_affine is None:
 444            out_affine = list(tiles[0].affine)
 445            out_affine[0] = min_x
 446            out_affine[3] = max_y
 447
 448        if blend == "first": # fill output, first come, first served.
 449            out = np.full( tuple(out_shape) + (tiles[0].band_count(),), np.nan, dtype=np.float32)
 450            for t in tiles:
 451                r = resample_raster( t.data, t.affine, out_affine, out_shape )
 452                mask = np.isnan(out)
 453                out[mask] = r[mask]
 454        elif blend == "min": # keep minimum value in case of overlap
 455            out = np.nanmin( np.stack([ resample_raster( t.data, t.affine, out_affine, out_shape ) for t in tiles ], 
 456                                      axis=0), axis=0 )
 457        elif blend == "max": # keep maximum value in case of overlap
 458            out = np.nanmax( np.stack([ resample_raster( t.data, t.affine, out_affine, out_shape ) for t in tiles ], 
 459                                      axis=0), axis=0 )
 460        elif blend == "mean": # use average in case of overlap
 461            out = np.nanmean( np.stack([ resample_raster( t.data, t.affine, out_affine, out_shape ) for t in tiles ], 
 462                                       axis=0), axis=0 )
 463        elif blend == "median": # use mean in case of overlap
 464            out = np.nanmedian( np.stack([ resample_raster( t.data, t.affine, out_affine, out_shape ) for t in tiles ], 
 465                                         axis=0), axis=0 )
 466
 467        # Return HyImage
 468        out =  HyImage(
 469            out,
 470            affine=out_affine,
 471            projection=tiles[0].projection,
 472            wav=tiles[0].get_wavelengths(),
 473            header=tiles[0].header.copy()
 474        )
 475        if 'xleft' in out.header: del out.header['xleft'] # stored in tiles, but not meaningful here
 476        if 'ytop' in out.header: del out.header['ytop'] # stored in tiles, but not meaningful here
 477
 478        return out
 479
 480    #####################################
 481    ## BASIC TRANSFORMS
 482    #####################################
 483
 484    def flip(self, axis='x'):
 485        """
 486        Flip the image on the x or y axis. Note that this will remove any defined affine transform.
 487
 488        Args:
 489            axis (str): 'x' or 'y' or both 'xy'.
 490        """
 491
 492        if 'x' in axis.lower():
 493            self.data = np.flip(self.data,axis=0)
 494        if 'y' in axis.lower():
 495            self.data = np.flip(self.data,axis=1)
 496        self.affine = None
 497        if 'affine' in self.header: del self.header['affine']
 498        self.push_to_header() # update width and height info
 499
 500    def rot90(self):
 501        """
 502        Rotate this image by 90 degrees by transposing the underlying data array. Combine with flip('x') or flip('y')
 503        to achieve positive/negative rotations.
 504        """
 505        self.data = np.transpose( self.data, (1,0,2) )
 506        self.affine = None
 507        if 'affine' in self.header: del self.header['affine']
 508        self.push_to_header() # update width and height info
 509
 510    #####################################
 511    ##IMAGE FILTERING
 512    #####################################
 513    def fill_holes(self):
 514        """
 515        Replaces nan pixel with an average of their neighbours, thus removing 1-pixel large holes from an image. Note that
 516        for performance reasons this assumes that holes line up across bands. Note that this is not vectorized so very slow...
 517        """
 518
 519        # perform greyscale dilation
 520        ndimage = require("scipy").ndimage
 521        dilate = self.data.copy()
 522        mask = np.logical_not(np.isfinite(dilate))
 523        dilate[mask] = 0
 524        for b in range(self.band_count()):
 525            dilate[:, :, b] = ndimage.grey_dilation(dilate[:, :, b], size=(3, 3))
 526
 527        # map back to holes in dataset
 528        self.data[mask] = dilate[mask]
 529        #self.data[self.data == 0] = np.nan  # replace remaining 0's with nans
 530
 531    def blur(self, n=3):
 532        """
 533        Applies a gaussian kernel of size n to the image using OpenCV.
 534
 535        Args:
 536            n (int): the dimensions of the gaussian kernel to convolve. Default is 3. Increase for more blurry results.
 537        """
 538        import cv2 # import this here to avoid errors if opencv is not installed properly
 539
 540        nanmask = np.isnan(self.data)
 541        assert isinstance(n, int) and n >= 3, "Error - invalid kernel. N must be an integer > 3. "
 542        kernel = np.ones((n, n), np.float32) / (n ** 2)
 543        self.data = cv2.filter2D(self.data, -1, kernel)
 544        self.data[nanmask] = np.nan  # remove mask
 545
 546    def erode(self, size=3, iterations=1):
 547        """
 548        Apply an erode filter to this image to expand background (nan) pixels. Refer to open-cv's erode
 549        function for more details.
 550
 551        Args:
 552            size (int): the size of the erode filter. Default is a 3x3 kernel.
 553            iterations (int): the number of erode iterations. Default is 1.
 554        """
 555        import cv2 # import this here to avoid errors if opencv is not installed properly
 556
 557        # erode
 558        kernel = np.ones((size, size), np.uint8)
 559        if self.is_float():
 560            mask = np.isfinite(self.data).any(axis=-1)
 561            mask = cv2.erode(mask.astype(np.uint8), kernel, iterations=iterations)
 562            self.data[mask == 0, :] = np.nan
 563        else:
 564            mask = (self.data != 0).any( axis=-1 )
 565            mask = cv2.erode(mask.astype(np.uint8), kernel, iterations=iterations)
 566            self.data[mask == 0, :] = 0
 567
 568    def despeckle(self, size=5):
 569        """
 570        Despeckle each band of this image (independently) using a median filter.
 571
 572        Args:
 573            size (int): the size of the median filter kernel. Default is 5. Must be an odd number.
 574        """
 575
 576        assert (size % 2) == 1, "Error - size must be an odd integer"
 577        import cv2 # import this here to avoid errors if opencv is not installed properly
 578        if self.is_float():
 579            self.data = cv2.medianBlur( self.data.astype(np.float32), size )
 580        else:
 581            self.data = cv2.medianBlur( self.data, size )
 582
 583    #####################################
 584    ##FEATURES AND FEATURE MATCHING
 585    ######################################
 586    def get_keypoints(self, band, eq=False, mask=True, method='sift', cfac=0.0,bfac=0.0, **kwds):
 587        """
 588        Get feature descriptors from the specified band.
 589
 590        Args:
 591            band (int,float,str,tuple): the band index (int) or wavelength (float) to extract features from. Alternatively, a tuple can be passed
 592                    containing a range of bands (min : max) to average before feature matching.
 593            eq (bool): True if the image should be histogram equalized first. Default is False.
 594            mask (bool): True if 0 value pixels should be masked. Default is True.
 595            method (str): the feature detector to use. Options are 'SIFT' and 'ORB' (faster but less accurate). Default is 'SIFT'.
 596            cfac (float): contrast adjustment to apply to hyperspectral bands before matching. Default is 0.0.
 597            bfac (float): brightness adjustment to apply to hyperspectral bands before matching. Default is 0.0.
 598            **kwds: keyword arguments are passed to the opencv feature detector. For SIFT these are:
 599
 600                - contrastThreshold: default is 0.01.
 601                - edgeThreshold: default is 10.
 602                - sigma: default is 1.0
 603
 604                For ORB these are:
 605
 606                - nfeatures = the number of features to detect. Default is 5000.
 607
 608            Returns:
 609                Tuple containing
 610
 611                    - k (ndarray): the keypoints detected
 612                    - d (ndarray): corresponding feature descriptors
 613         """
 614        import cv2 # import this here to avoid errors if opencv is not installed properly
 615
 616        # get image
 617        if isinstance(band, int) or isinstance(band, float): #single band
 618            image = self.data[:, :, self.get_band_index(band)]
 619        elif isinstance(band,tuple): #range of bands (averaged)
 620            idx0 = self.get_band_index(band[0])
 621            idx1 = self.get_band_index(band[1])
 622
 623            #deal with out of range errors
 624            if idx0 is None:
 625                idx0 = 0
 626            if idx1 is None:
 627                idx1 = self.band_count()
 628
 629            #average bands
 630            image = np.nanmean(self.data[:,:,idx0:idx1],axis=2)
 631        else:
 632            assert False, "Error, unrecognised band %s" % band
 633
 634        #normalise image to range 0 - 1
 635        image -= np.nanmin(image)
 636        image = image / np.nanmax(image)
 637
 638        #apply brightness/contrast adjustment
 639        image = (1.0+cfac)*image + bfac
 640        image[image > 1.0] = 1.0
 641        image[image < 0.0] = 0.0
 642
 643        #convert image to uint8 for opencv
 644        image = np.uint8(255 * image)
 645        if eq:
 646            image = cv2.equalizeHist(image)
 647
 648        if mask:
 649            mask = np.zeros(image.shape, dtype=np.uint8)
 650            mask[image != 0] = 255  # include only non-zero pixels
 651        else:
 652            mask = None
 653
 654        if 'sift' in method.lower():  # SIFT
 655
 656            # setup default keywords
 657            kwds["contrastThreshold"] = kwds.get("contrastThreshold", 0.01)
 658            kwds["edgeThreshold"] = kwds.get("edgeThreshold", 10)
 659            kwds["sigma"] = kwds.get("sigma", 1.0)
 660
 661            # make feature detector
 662            #alg = cv2.xfeatures2d.SIFT_create(**kwds)
 663            alg = cv2.SIFT_create()
 664        elif 'orb' in method.lower():  # orb
 665            kwds['nfeatures'] = kwds.get('nfeatures', 5000)
 666            alg = cv2.ORB_create(scoreType=cv2.ORB_FAST_SCORE, **kwds)
 667        else:
 668            assert False, "Error - %s is not a recognised feature detector." % method
 669
 670        # detect keypoints
 671        kp = alg.detect(image, mask)
 672
 673        # extract and return feature vectors
 674        return alg.compute(image, kp)
 675
 676    @classmethod
 677    def match_keypoints(cls, kp1, kp2, d1, d2, method='SIFT', dist=0.7, tree = 5, check = 100, min_count=5):
 678        """
 679        Compares keypoint feature vectors from two images and returns matching pairs.
 680
 681        Args:
 682            kp1 (ndarray): keypoints from the first image
 683            kp2 (ndarray): keypoints from the second image
 684            d1 (ndarray): descriptors for the keypoints from the first image
 685            d2 (ndarray): descriptors for the keypoints from the second image
 686            method (str): the method used to calculate the feature descriptors. Should be 'sift' or 'orb'. Default is 'sift'.
 687            dist (float): minimum match distance (0 to 1), default is 0.7
 688            tree (int): not sure what this does? Default is 5. See open-cv docs.
 689            check (int): ditto. Default is 100.
 690            min_count (int): the minimum number of matches to consider a valid matching operation. If fewer matches are found,
 691                       then the function returns None, None. Default is 5.
 692        """
 693        import cv2 # import this here to avoid errors if opencv is not installed properly
 694        if 'sift' in method.lower():
 695            algorithm = cv2.NORM_INF
 696        elif 'orb' in method.lower():
 697            algorithm = cv2.NORM_HAMMING
 698        else:
 699            assert False, "Error - unknown matching algorithm %s" % method
 700
 701        #calculate flann matches
 702        index_params = dict(algorithm=algorithm, trees=tree)
 703        search_params = dict(checks=check)
 704        flann = cv2.FlannBasedMatcher(index_params, search_params)
 705        matches = flann.knnMatch(d1, d2, k=2)
 706
 707        # store all the good matches as per Lowe's ratio test.
 708        good = []
 709        for m, n in matches:
 710            if m.distance < dist * n.distance:
 711                good.append(m)
 712
 713        if len(good) < min_count:
 714            return None, None
 715        else:
 716            src_pts = np.float32([kp1[m.queryIdx].pt for m in good]).reshape(-1, 1, 2)
 717            dst_pts = np.float32([kp2[m.trainIdx].pt for m in good]).reshape(-1, 1, 2)
 718            return src_pts, dst_pts
 719
 720    ############################
 721    ## Visualisation methods
 722    ############################
 723    def quick_plot(self, bands=0, ax=None, bfac=0.0, cfac=0.0, samples=False, tscale=False, invert=False, rot=False, flipX=False, flipY=False,
 724                   **kwds):
 725        """
 726        Plot a band using matplotlib.imshow(...).
 727
 728        Args:
 729            bands (str,int,float,tuple): the band name (string), index (integer) or wavelength (float) to plot. Default is 0. If a tuple is passed then
 730                  each band in the tuple (string or index) will be mapped to rgb. Bands with negative wavelengths or indices will be inverted before plotting.
 731            ax: an axis object to plot to. If none, plt.imshow( ... ) is used.
 732            bfac (float): a brightness adjustment to apply to RGB mappings (-1 to 1)
 733            cfac (float): a contrast adjustment to apply to RGB mappings (-1 to 1)
 734            samples (bool): True if sample points (defined in the header file) should be plotted. Default is False. Otherwise, a list of
 735                     [ (x,y), ... ] points can be passed.
 736            tscale (bool): True if each band (for ternary images) should be scaled independently. Default is False.
 737                    When using scaling, vmin and vmax can be used to set the clipping percentiles (integers) or
 738                    (constant) values (float).
 739            invert (bool) : True if each band should be inverted before plotting. Only works for multiband (ternary) images.
 740            rot (bool): if True, the x and y axis will be flipped (90 degree rotation) before plotting. Default is False.
 741            flipX (bool): if True, the x axis will be flipped before plotting (after applying rotations).
 742            flipY (bool): if True, the y axis will be flippe before plotting (after applying rotations).
 743            **kwds: keywords are passed to matplotlib.imshow( ... ), except for the following:
 744
 745                 - mask = a 2D boolean mask containing true if pixels should be drawn and false otherwise.
 746                 - path = a file path to save the image too (at matching resolution; use fig.savefig(..) if you want to save the figure).
 747                 - ticks = True if x- and y- ticks should be plotted. Default is False.
 748                 - ps, pc = the size and color of sample points to plot. Can be constant or list.
 749                 - figsize = a figsize for the figure to create (if ax is None).
 750
 751        Returns:
 752            Tuple containing
 753
 754            - fig: matplotlib figure object
 755            - ax:  matplotlib axes object. If a colorbar is created, (band is an integer or a float), then this will be stored in ax.cbar.
 756        """
 757
 758        plt = require("matplotlib.pyplot")
 759
 760        #create new axes?
 761        if ax is None:
 762            fig, ax = plt.subplots(figsize=kwds.pop('figsize', (18,18*self.ydim()/self.xdim()) ))
 763
 764        # deal with ticks
 765        if not kwds.pop('ticks', False ):
 766            ax.set_xticks([])
 767            ax.set_yticks([])
 768
 769        #map individual band using colourmap
 770        if isinstance(bands, str) or isinstance(bands, int) or isinstance(bands, float):
 771            #get band
 772            if isinstance(bands, str):
 773                data = self.data[:, :, self.get_band_index(bands)]
 774            else:
 775                data = self.data[:, :, self.get_band_index(np.abs(bands))]
 776            if not isinstance(bands, str) and bands < 0:
 777                data = np.nanmax(data) - data # flip
 778
 779            # convert integer vmin and vmax values to percentiles
 780            if 'vmin' in kwds:
 781                if isinstance(kwds['vmin'], int):
 782                    kwds['vmin'] = np.nanpercentile( data, kwds['vmin'] )
 783            if 'vmax' in kwds:
 784                if isinstance(kwds['vmax'], int):
 785                    kwds['vmax'] = np.nanpercentile( data, kwds['vmax'] )
 786
 787            #mask nans (and apply custom mask)
 788            mask = np.isnan(data)
 789            if not np.isnan(self.header.get_data_ignore_value()):
 790                mask = mask + data == self.header.get_data_ignore_value()
 791            if 'mask' in kwds:
 792                mask = mask + kwds.get('mask')
 793                del kwds['mask']
 794            data = np.ma.array(data, mask = mask > 0 )
 795
 796            # apply rotations and flipping
 797            if rot:
 798                data = data.T
 799            if flipX:
 800                data = data[::-1, :]
 801            if flipY:
 802                data = data[:, ::-1]
 803
 804            # save?
 805            if 'path' in kwds:
 806                path = kwds.pop('path')
 807                imsave = require("matplotlib.pyplot").imsave
 808                if not os.path.exists(os.path.dirname(path)):
 809                    os.makedirs(os.path.dirname(path)) # ensure output directory exists
 810                imsave(path, data.T, **kwds)  # save the image
 811
 812            ax.cbar = ax.imshow(data.T, interpolation=kwds.pop('interpolation', 'none'), **kwds) # change default interpolation to None
 813
 814        #map 3 bands to RGB
 815        elif isinstance(bands, tuple) or isinstance(bands, list):
 816            #get band indices and range
 817            rgb = []
 818            for b in bands:
 819                if isinstance(b, str):
 820                    rgb.append(self.get_band_index(b))
 821                else:
 822                    rgb.append(self.get_band_index(np.abs(b)))
 823
 824            #slice image (as copy) and map to 0 - 1
 825            img = np.array(self.data[:, :, rgb]).copy()
 826            if np.isnan(img).all():
 827                print("Warning - image contains no data.")
 828                return ax.get_figure(), ax
 829
 830            # invert if needed
 831            if invert:
 832                bands = [-b for b in bands]
 833            for i,b in enumerate(bands):
 834                if not isinstance(b, str) and (b < 0):
 835                    img[..., i] = np.nanmax(img[..., i]) - img[..., i]
 836
 837            # do scaling
 838            if tscale: # scale bands independently
 839                for b in range(3):
 840                    mn = kwds.get("vmin", float(np.nanmin(img)))
 841                    mx = kwds.get("vmax", float(np.nanmax(img)))
 842                    if isinstance (mn, int):
 843                        assert mn >= 0 and mn <= 100, "Error - integer vmin values must be a percentile."
 844                        mn = float(np.nanpercentile(img[...,b], mn ))
 845                    if isinstance (mx, int):
 846                        assert mx >= 0 and mx <= 100, "Error - integer vmax values must be a percentile."
 847                        mx = float(np.nanpercentile(img[...,b], mx ))
 848                    img[...,b] = (img[..., b] - mn) / (mx - mn)
 849            else: # scale bands together
 850                mn = kwds.get("vmin", float(np.nanmin(img)))
 851                mx = kwds.get("vmax", float(np.nanmax(img)))
 852                if isinstance(mn, int):
 853                    assert mn >= 0 and mn <= 100, "Error - integer vmin values must be a percentile."
 854                    mn = float(np.nanpercentile(img, mn))
 855                if isinstance(mx, int):
 856                    assert mx >= 0 and mx <= 100, "Error - integer vmax values must be a percentile."
 857                    mx = float(np.nanpercentile(img, mx))
 858                img = (img - mn) / (mx - mn)
 859
 860            #apply brightness/contrast mapping
 861            img = np.clip((1.0 + cfac) * img + bfac, 0, 1.0 )
 862
 863            #apply masking so background is white
 864            img[np.logical_not( np.isfinite( img ) )] = 1.0
 865            if 'mask' in kwds:
 866                img[kwds.pop("mask"),:] = 1.0
 867
 868            # apply rotations and flipping
 869            if rot:
 870                img = np.transpose( img, (1,0,2) )
 871            if flipX:
 872                img = img[::-1, :, :]
 873            if flipY:
 874                img = img[:, ::-1, :]
 875
 876            # save?
 877            if 'path' in kwds:
 878                path = kwds.pop('path')
 879                imsave = require("matplotlib.pyplot").imsave
 880                if not os.path.exists(os.path.dirname(path)):
 881                    os.makedirs(os.path.dirname(path)) # ensure output directory exists
 882                imsave(path, np.transpose( np.clip( img*255, 0, 255).astype(np.uint8), (1, 0, 2)))  # save the image
 883
 884            # plot samples?
 885            ps = kwds.pop('ps', 5)
 886            pc = kwds.pop('pc', 'r')
 887            if samples:
 888                if isinstance(samples, list) or isinstance(samples, np.ndarray):
 889                    ax.scatter([s[0] for s in samples], [s[1] for s in samples], s=ps, c=pc)
 890                else:
 891                    for n in self.header.get_class_names():
 892                        points = np.array(self.header.get_sample_points(n))
 893                        ax.scatter(points[:, 0], points[:, 1], s=ps, c=pc)
 894
 895            #plot
 896            ax.imshow(np.transpose(img, (1,0,2)), interpolation=kwds.pop('interpolation', 'none'), **kwds)
 897            ax.cbar = None  # no colorbar
 898
 899        return ax.get_figure(), ax
 900
 901    ## masking
 902    def mask(self, mask=None, flag=np.nan, invert=False, crop=False, bands=None):
 903        """
 904         Apply a mask to an image, flagging masked pixels with the specified value. Note that this applies the mask to the
 905         image in-situ.
 906
 907         Args:
 908            flag (float): the value to use for masked pixels. Default is np.nan
 909            mask (ndarray): a numpy array defining the mask polygon of the format [[x1,y1],[x2,y2],...]. If None is passed then
 910                    pickPolygon( ... ) is used to interactively define a polygon. If a file path is passed then the polygon
 911                    will be loaded using np.load( ... ). Alternatively if mask.shape == image.shape[0,1] then it is treated as a
 912                    binary image mask (must be boolean) and True values will be masked across all bands. Default is None.
 913            invert (bool): if True, pixels within the polygon will be masked. If False, pixels outside the polygon are masked. Default is False.
 914            crop (bool): True if rows/columns containing only zeros should be removed. Default is False.
 915            bands (tuple): the bands of the image to plot if no mask is specified. If None, the middle band is used.
 916
 917         Returns:
 918            Tuple containing
 919
 920            - mask (ndarray): a boolean array with True where pixels are masked and False elsewhere.
 921            - poly (ndarray): the mask polygon array in the format described above. Useful if the polygon was interactively defined.
 922         """
 923
 924        if mask is None:  # pick mask interactively
 925            if bands is None:
 926                bands = int(self.band_count() / 2)
 927
 928            regions = self.pickPolygons(region_names=["mask"], bands=bands)
 929
 930            # the user bailed without picking a mask?
 931            if len(regions) == 0:
 932                print("Warning - no mask picked/applied.")
 933                return
 934
 935            # extract polygon mask
 936            mask = regions[0]
 937
 938        # convert polygon mask to binary mask
 939        if mask.shape[1] == 2:
 940
 941            # build meshgrid with pixel coords
 942            xx, yy = np.meshgrid(np.arange(self.xdim()), np.arange(self.ydim()))
 943            xx = xx.flatten()
 944            yy = yy.flatten()
 945            points = np.vstack([xx, yy]).T  # coordinates of each pixel
 946
 947            # calculate per-pixel mask
 948            MplPath = require("matplotlib.path").Path
 949            mask = MplPath(mask).contains_points(points)
 950            mask = mask.reshape((self.ydim(), self.xdim())).T
 951
 952            # flip as we want to mask (==True) outside points (unless invert is true)
 953            if not invert:
 954                mask = np.logical_not(mask)
 955
 956        # apply binary image mask
 957        assert mask.shape[0] == self.data.shape[0] and mask.shape[1] == self.data.shape[1], \
 958            "Error - mask shape %s does not match image shape %s" % (mask.shape, self.data.shape)
 959        for b in range(self.band_count()):
 960            self.data[:, :, b][mask] = flag
 961
 962        # crop image
 963        if crop:
 964            self.crop_to_data()
 965
 966        return mask
 967
 968    def crop_to_data(self):
 969        """
 970        Remove padding of nan or zero pixels from image. Note that this is performed in place.
 971        """
 972        valid = np.isfinite(self.data).any(axis=-1) & (self.data != 0).any(axis=-1)
 973
 974        # integrate along axes
 975        xdata = np.sum(valid, axis=1) > 0.0
 976        ydata = np.sum(valid, axis=0) > 0.0
 977
 978        # calculate domain containing valid pixels
 979        xmin = np.argmax(xdata)
 980        xmax = xdata.shape[0] - np.argmax(xdata[::-1])
 981        ymin = np.argmax(ydata)
 982        ymax = ydata.shape[0] - np.argmax(ydata[::-1])
 983
 984        # crop
 985        self.data = self.data[xmin:xmax, ymin:ymax, :]
 986        
 987        # shift affine origin to new top-left pixel
 988        if self.affine is not None: 
 989            a = self.affine # shorthand for affine
 990            new_affine = list(self.affine)
 991            new_affine[0] = a[0] + xmin*a[1] + ymin*a[2]
 992            new_affine[3] = a[3] + xmin*a[4] + ymin*a[5]
 993            self.affine = np.array(new_affine)
 994            self.header['affine'] = self.affine
 995            
 996    ##################################################
 997    ## Interactive tools for picking regions/pixels
 998    ##################################################
 999    def pickPolygons(self, region_names, bands=0):
1000        """
1001        Creates a matplotlib gui for selecting polygon regions in an image.
1002
1003        Args:
1004            names (list, str): a list containing the names of the regions to pick. If a string is passed only one name is used.
1005            bands (tuple): the bands of the image to plot.
1006        """
1007
1008        if isinstance(region_names, str):
1009            region_names = [region_names]
1010
1011        assert isinstance(region_names, list), "Error - names must be a list or a string."
1012
1013        matplotlib = require("matplotlib")
1014        plt = require("matplotlib.pyplot")
1015        MultiRoi = require("roipoly").MultiRoi
1016
1017        # set matplotlib backend
1018        backend = matplotlib.get_backend()
1019        matplotlib.use('Qt5Agg')  # need this backend for ROIPoly to work
1020
1021        # plot image and extract roi's
1022        fig, ax = self.quick_plot(bands)
1023        roi = MultiRoi(roi_names=region_names)
1024        plt.close(fig)  # close figure
1025
1026        # extract regions
1027        regions = []
1028        for name, r in roi.rois.items():
1029            # store region
1030            x = r.x
1031            y = r.y
1032            regions.append(np.vstack([x, y]).T)
1033
1034        # restore matplotlib backend (if possible)
1035        try:
1036            matplotlib.use(backend)
1037        except:
1038            print("Warning: could not reset matplotlib backend. Plots will remain interactive...")
1039            pass
1040
1041        return regions
1042
1043    def pickPoints(self, n=-1, bands=hylite.RGB, integer=True, title="Pick Points", **kwds):
1044        """
1045        Creates a matplotlib gui for picking pixels from an image.
1046
1047        Args:
1048            n (int): the number of pixels to pick, or -1 if the user can select as many as they wish. Default is -1.
1049            bands (tuple): the bands of the image to plot. Default is `hylite.hyimage.HyImage`.RGB
1050            integer (bool): True if points coordinates should be cast to integers (for use as indices). Default is True.
1051            title (str): The title of the point picking window.
1052            **kwds: Keywords are passed to `hylite.hyimage.HyImage`.quick_plot( ... ).
1053
1054        Returns:
1055            A list containing the picked point coordinates [ (x1,y1), (x2,y2), ... ].
1056        """
1057
1058        matplotlib = require("matplotlib")
1059        plt = require("matplotlib.pyplot")
1060
1061        # set matplotlib backend
1062        backend = matplotlib.get_backend()
1063        matplotlib.use('Qt5Agg')  # need this backend for ROIPoly to work
1064
1065        # create figure
1066        fig, ax = self.quick_plot( bands, **kwds )
1067        ax.set_title(title)
1068
1069        # get points
1070        points = fig.ginput( n )
1071
1072        if integer:
1073            points = [ (int(p[0]), int(p[1])) for p in points ]
1074
1075        # restore matplotlib backend (if possible)
1076        try:
1077            matplotlib.use(backend)
1078        except:
1079            print("Warning: could not reset matplotlib backend. Plots will remain interactive...")
1080            pass
1081
1082        return points
1083
1084    def pickSamples(self, names=None, store=True, **kwds):
1085        """
1086        Pick sample probe points and store these in the image header file.
1087
1088        Args:
1089            names (str, list): the name of the sample to pick, or a list of names to pick multiple.
1090            store (bool): True if sample should be stored in the image header file (for later access). Default is True.
1091            **kwds: Keywords are passed to `hylite.hyimage.HyImage`.quick_plot( ... )
1092
1093        Returns:
1094            a list containing a list of points for each sample.
1095        """
1096
1097        if isinstance(names, str):
1098            names = [names]
1099
1100        # pick points
1101        points = []
1102        for s in names:
1103            pnts = self.pickPoints(title="%s" % s, **kwds)
1104            if store:
1105                self.header['sample %s' % s] = pnts # store in header
1106            points.append(pnts)
1107        # add class to header file
1108        if store:
1109            cls_names = self.header.get_class_names()
1110            if cls_names is None:
1111                cls_names = []
1112            self.header['class names'] = cls_names + names
1113
1114        return points