hylite.hyimage

  1"""
  2Store and manipulate hyperspectral image data.
  3"""
  4
  5import os
  6import numpy as np
  7import matplotlib
  8import matplotlib.pyplot as plt
  9from matplotlib import path
 10from roipoly import MultiRoi
 11import imageio
 12import scipy as sp
 13import hylite
 14from hylite.hydata import HyData
 15from hylite.hylibrary import HyLibrary
 16
 17
 18
 19class HyImage( HyData ):
 20    """
 21    A class for hyperspectral image data. These can be individual scenes or hyperspectral orthoimages.
 22    """
 23
 24    def __init__(self, data, **kwds):
 25        """
 26        Args:
 27            data (ndarray): a numpy array such that data[x][y][band] gives each pixel value.
 28            **kwds:
 29                wav = A numpy array containing band wavelengths for this image.
 30                affine = an affine transform of the format returned by GDAL.GetGeoTransform().
 31                projection = string defining the project. Default is None.
 32                sensor = sensor name. Default is "unknown".
 33                header = path to associated header file. Default is None.
 34        """
 35
 36        #call constructor for HyData
 37        super().__init__(data, **kwds)
 38
 39        # special case - if dataset only has oneband, slice it so it still has
 40        # the format data[x,y,b].
 41        if not self.data is None:
 42            if len(self.data.shape) == 1:
 43                self.data = self.data[None, None, :] # single pixel image
 44            if len(self.data.shape) == 2:
 45                self.data = self.data[:, :, None] # single band iamge
 46
 47        #load any additional project information (specific to images)
 48        self.set_projection(kwds.get("projection",None))
 49        self.affine = kwds.get("affine",[0,1,0,0,0,1])
 50
 51        # wavelengths
 52        if 'wav' in kwds:
 53            self.set_wavelengths(kwds['wav'])
 54
 55        #special header formatting
 56        self.header['file type'] = 'ENVI Standard'
 57
 58    def copy(self,data=True):
 59        """
 60        Make a deep copy of this image instance.
 61
 62        Args:
 63            data (bool): True if a copy of the data should be made, otherwise only copy header.
 64
 65        Returns:
 66            a new HyImage instance.
 67        """
 68        if not data:
 69            return HyImage(None, header=self.header.copy(), projection=self.projection, affine=self.affine)
 70        else:
 71            return HyImage( self.data.copy(), header=self.header.copy(), projection=self.projection, affine=self.affine)
 72
 73    def T(self):
 74        """
 75        Return a transposed view of the data matrix (corresponding with the [y,x] indexing used by matplotlib, opencv etc.
 76        """
 77        return np.transpose(self.data, (1,0,2))
 78
 79    def xdim(self):
 80        """
 81        Return number of pixels in x (first dimension of data array)
 82        """
 83        return self.data.shape[0]
 84
 85    def ydim(self):
 86        """
 87        Return number of pixels in y (second dimension of data array)
 88        """
 89        return self.data.shape[1]
 90
 91    def aspx(self):
 92        """
 93        Return the aspect ratio of this image (width/height).
 94        """
 95        return self.ydim() / self.xdim()
 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 (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 (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 flip(self, axis='x'):
273        """
274        Flip the image on the x or y axis.
275
276        Args:
277            axis (str): 'x' or 'y' or both 'xy'.
278        """
279
280        if 'x' in axis.lower():
281            self.data = np.flip(self.data,axis=0)
282        if 'y' in axis.lower():
283            self.data = np.flip(self.data,axis=1)
284
285    def rot90(self):
286        """
287        Rotate this image by 90 degrees by transposing the underlying data array. Combine with flip('x') or flip('y')
288        to achieve positive/negative rotations.
289        """
290        self.data = np.transpose( self.data, (1,0,2) )
291        self.push_to_header()
292
293    #####################################
294    ##IMAGE FILTERING
295    #####################################
296    def fill_holes(self):
297        """
298        Replaces nan pixel with an average of their neighbours, thus removing 1-pixel large holes from an image. Note that
299        for performance reasons this assumes that holes line up across bands. Note that this is not vectorized so very slow...
300        """
301
302        # perform greyscale dilation
303        dilate = self.data.copy()
304        mask = np.logical_not(np.isfinite(dilate))
305        dilate[mask] = 0
306        for b in range(self.band_count()):
307            dilate[:, :, b] = sp.ndimage.grey_dilation(dilate[:, :, b], size=(3, 3))
308
309        # map back to holes in dataset
310        self.data[mask] = dilate[mask]
311        #self.data[self.data == 0] = np.nan  # replace remaining 0's with nans
312
313    def blur(self, n=3):
314        """
315        Applies a gaussian kernel of size n to the image using OpenCV.
316
317        Args:
318            n (int): the dimensions of the gaussian kernel to convolve. Default is 3. Increase for more blurry results.
319        """
320        import cv2 # import this here to avoid errors if opencv is not installed properly
321
322        nanmask = np.isnan(self.data)
323        assert isinstance(n, int) and n >= 3, "Error - invalid kernel. N must be an integer > 3. "
324        kernel = np.ones((n, n), np.float32) / (n ** 2)
325        self.data = cv2.filter2D(self.data, -1, kernel)
326        self.data[nanmask] = np.nan  # remove mask
327
328    def erode(self, size=3, iterations=1):
329        """
330        Apply an erode filter to this image to expand background (nan) pixels. Refer to open-cv's erode
331        function for more details.
332
333        Args:
334            size (int): the size of the erode filter. Default is a 3x3 kernel.
335            iterations (int): the number of erode iterations. Default is 1.
336        """
337        import cv2 # import this here to avoid errors if opencv is not installed properly
338
339        # erode
340        kernel = np.ones((size, size), np.uint8)
341        if self.is_float():
342            mask = np.isfinite(self.data).any(axis=-1)
343            mask = cv2.erode(mask.astype(np.uint8), kernel, iterations=iterations)
344            self.data[mask == 0, :] = np.nan
345        else:
346            mask = (self.data != 0).any( axis=-1 )
347            mask = cv2.erode(mask.astype(np.uint8), kernel, iterations=iterations)
348            self.data[mask == 0, :] = 0
349
350    def resize(self, newdims : tuple, interpolation : int = 1):
351        """
352        Resize this image with opencv.
353
354        Args:
355            newdims (tuple): the new image dimensions.
356            interpolation (int): opencv interpolation method. Default is cv2.INTER_LINEAR.
357        """
358        import cv2 # import this here to avoid errors if opencv is not installed properly
359        self.data = cv2.resize(self.data, (newdims[1],newdims[0]), interpolation=interpolation)
360
361    def despeckle(self, size=5):
362        """
363        Despeckle each band of this image (independently) using a median filter.
364
365        Args:
366            size (int): the size of the median filter kernel. Default is 5. Must be an odd number.
367        """
368
369        assert (size % 2) == 1, "Error - size must be an odd integer"
370        import cv2 # import this here to avoid errors if opencv is not installed properly
371        if self.is_float():
372            self.data = cv2.medianBlur( self.data.astype(np.float32), size )
373        else:
374            self.data = cv2.medianBlur( self.data, size )
375
376    #####################################
377    ##FEATURES AND FEATURE MATCHING
378    ######################################
379    def get_keypoints(self, band, eq=False, mask=True, method='sift', cfac=0.0,bfac=0.0, **kwds):
380        """
381        Get feature descriptors from the specified band.
382
383        Args:
384            band (int,float,str,tuple): the band index (int) or wavelength (float) to extract features from. Alternatively, a tuple can be passed
385                    containing a range of bands (min : max) to average before feature matching.
386            eq (bool): True if the image should be histogram equalized first. Default is False.
387            mask (bool): True if 0 value pixels should be masked. Default is True.
388            method (str): the feature detector to use. Options are 'SIFT' and 'ORB' (faster but less accurate). Default is 'SIFT'.
389            cfac (float): contrast adjustment to apply to hyperspectral bands before matching. Default is 0.0.
390            bfac (float): brightness adjustment to apply to hyperspectral bands before matching. Default is 0.0.
391            **kwds: keyword arguments are passed to the opencv feature detector. For SIFT these are:
392
393                - contrastThreshold: default is 0.01.
394                - edgeThreshold: default is 10.
395                - sigma: default is 1.0
396
397                For ORB these are:
398
399                - nfeatures = the number of features to detect. Default is 5000.
400
401            Returns:
402                Tuple containing
403
404                    - k (ndarray): the keypoints detected
405                    - d (ndarray): corresponding feature descriptors
406         """
407        import cv2 # import this here to avoid errors if opencv is not installed properly
408
409        # get image
410        if isinstance(band, int) or isinstance(band, float): #single band
411            image = self.data[:, :, self.get_band_index(band)]
412        elif isinstance(band,tuple): #range of bands (averaged)
413            idx0 = self.get_band_index(band[0])
414            idx1 = self.get_band_index(band[1])
415
416            #deal with out of range errors
417            if idx0 is None:
418                idx0 = 0
419            if idx1 is None:
420                idx1 = self.band_count()
421
422            #average bands
423            image = np.nanmean(self.data[:,:,idx0:idx1],axis=2)
424        else:
425            assert False, "Error, unrecognised band %s" % band
426
427        #normalise image to range 0 - 1
428        image -= np.nanmin(image)
429        image = image / np.nanmax(image)
430
431        #apply brightness/contrast adjustment
432        image = (1.0+cfac)*image + bfac
433        image[image > 1.0] = 1.0
434        image[image < 0.0] = 0.0
435
436        #convert image to uint8 for opencv
437        image = np.uint8(255 * image)
438        if eq:
439            image = cv2.equalizeHist(image)
440
441        if mask:
442            mask = np.zeros(image.shape, dtype=np.uint8)
443            mask[image != 0] = 255  # include only non-zero pixels
444        else:
445            mask = None
446
447        if 'sift' in method.lower():  # SIFT
448
449            # setup default keywords
450            kwds["contrastThreshold"] = kwds.get("contrastThreshold", 0.01)
451            kwds["edgeThreshold"] = kwds.get("edgeThreshold", 10)
452            kwds["sigma"] = kwds.get("sigma", 1.0)
453
454            # make feature detector
455            #alg = cv2.xfeatures2d.SIFT_create(**kwds)
456            alg = cv2.SIFT_create()
457        elif 'orb' in method.lower():  # orb
458            kwds['nfeatures'] = kwds.get('nfeatures', 5000)
459            alg = cv2.ORB_create(scoreType=cv2.ORB_FAST_SCORE, **kwds)
460        else:
461            assert False, "Error - %s is not a recognised feature detector." % method
462
463        # detect keypoints
464        kp = alg.detect(image, mask)
465
466        # extract and return feature vectors
467        return alg.compute(image, kp)
468
469    @classmethod
470    def match_keypoints(cls, kp1, kp2, d1, d2, method='SIFT', dist=0.7, tree = 5, check = 100, min_count=5):
471        """
472        Compares keypoint feature vectors from two images and returns matching pairs.
473
474        Args:
475            kp1 (ndarray): keypoints from the first image
476            kp2 (ndarray): keypoints from the second image
477            d1 (ndarray): descriptors for the keypoints from the first image
478            d2 (ndarray): descriptors for the keypoints from the second image
479            method (str): the method used to calculate the feature descriptors. Should be 'sift' or 'orb'. Default is 'sift'.
480            dist (float): minimum match distance (0 to 1), default is 0.7
481            tree (int): not sure what this does? Default is 5. See open-cv docs.
482            check (int): ditto. Default is 100.
483            min_count (int): the minimum number of matches to consider a valid matching operation. If fewer matches are found,
484                       then the function returns None, None. Default is 5.
485        """
486        import cv2 # import this here to avoid errors if opencv is not installed properly
487        if 'sift' in method.lower():
488            algorithm = cv2.NORM_INF
489        elif 'orb' in method.lower():
490            algorithm = cv2.NORM_HAMMING
491        else:
492            assert False, "Error - unknown matching algorithm %s" % method
493
494        #calculate flann matches
495        index_params = dict(algorithm=algorithm, trees=tree)
496        search_params = dict(checks=check)
497        flann = cv2.FlannBasedMatcher(index_params, search_params)
498        matches = flann.knnMatch(d1, d2, k=2)
499
500        # store all the good matches as per Lowe's ratio test.
501        good = []
502        for m, n in matches:
503            if m.distance < dist * n.distance:
504                good.append(m)
505
506        if len(good) < min_count:
507            return None, None
508        else:
509            src_pts = np.float32([kp1[m.queryIdx].pt for m in good]).reshape(-1, 1, 2)
510            dst_pts = np.float32([kp2[m.trainIdx].pt for m in good]).reshape(-1, 1, 2)
511            return src_pts, dst_pts
512
513    ############################
514    ## Visualisation methods
515    ############################
516    def quick_plot(self, band=0, ax=None, bfac=0.0, cfac=0.0, samples=False, tscale=False, invert=False, rot=False, flipX=False, flipY=False,
517                   **kwds):
518        """
519        Plot a band using matplotlib.imshow(...).
520
521        Args:
522            band (str,int,float,tuple): the band name (string), index (integer) or wavelength (float) to plot. Default is 0. If a tuple is passed then
523                  each band in the tuple (string or index) will be mapped to rgb. Bands with negative wavelengths or indices will be inverted before plotting.
524            ax: an axis object to plot to. If none, plt.imshow( ... ) is used.
525            bfac (float): a brightness adjustment to apply to RGB mappings (-1 to 1)
526            cfac (float): a contrast adjustment to apply to RGB mappings (-1 to 1)
527            samples (bool): True if sample points (defined in the header file) should be plotted. Default is False. Otherwise, a list of
528                     [ (x,y), ... ] points can be passed.
529            tscale (bool): True if each band (for ternary images) should be scaled independently. Default is False.
530                    When using scaling, vmin and vmax can be used to set the clipping percentiles (integers) or
531                    (constant) values (float).
532            invert (bool) : True if each band should be inverted before plotting. Only works for multiband (ternary) images.
533            rot (bool): if True, the x and y axis will be flipped (90 degree rotation) before plotting. Default is False.
534            flipX (bool): if True, the x axis will be flipped before plotting (after applying rotations).
535            flipY (bool): if True, the y axis will be flippe before plotting (after applying rotations).
536            **kwds: keywords are passed to matplotlib.imshow( ... ), except for the following:
537
538                 - mask = a 2D boolean mask containing true if pixels should be drawn and false otherwise.
539                 - path = a file path to save the image too (at matching resolution; use fig.savefig(..) if you want to save the figure).
540                 - ticks = True if x- and y- ticks should be plotted. Default is False.
541                 - ps, pc = the size and color of sample points to plot. Can be constant or list.
542                 - figsize = a figsize for the figure to create (if ax is None).
543
544        Returns:
545            Tuple containing
546
547            - fig: matplotlib figure object
548            - ax:  matplotlib axes object. If a colorbar is created, (band is an integer or a float), then this will be stored in ax.cbar.
549        """
550
551        #create new axes?
552        if ax is None:
553            fig, ax = plt.subplots(figsize=kwds.pop('figsize', (18,18*self.ydim()/self.xdim()) ))
554
555        # deal with ticks
556        if not kwds.pop('ticks', False ):
557            ax.set_xticks([])
558            ax.set_yticks([])
559
560        #map individual band using colourmap
561        if isinstance(band, str) or isinstance(band, int) or isinstance(band, float):
562            #get band
563            if isinstance(band, str):
564                data = self.data[:, :, self.get_band_index(band)]
565            else:
566                data = self.data[:, :, self.get_band_index(np.abs(band))]
567            if not isinstance(band, str) and band < 0:
568                data = np.nanmax(data) - data # flip
569
570            # convert integer vmin and vmax values to percentiles
571            if 'vmin' in kwds:
572                if isinstance(kwds['vmin'], int):
573                    kwds['vmin'] = np.nanpercentile( data, kwds['vmin'] )
574            if 'vmax' in kwds:
575                if isinstance(kwds['vmax'], int):
576                    kwds['vmax'] = np.nanpercentile( data, kwds['vmax'] )
577
578            #mask nans (and apply custom mask)
579            mask = np.isnan(data)
580            if not np.isnan(self.header.get_data_ignore_value()):
581                mask = mask + data == self.header.get_data_ignore_value()
582            if 'mask' in kwds:
583                mask = mask + kwds.get('mask')
584                del kwds['mask']
585            data = np.ma.array(data, mask = mask > 0 )
586
587            # apply rotations and flipping
588            if rot:
589                data = data.T
590            if flipX:
591                data = data[::-1, :]
592            if flipY:
593                data = data[:, ::-1]
594
595            # save?
596            if 'path' in kwds:
597                path = kwds.pop('path')
598                from matplotlib.pyplot import imsave
599                if not os.path.exists(os.path.dirname(path)):
600                    os.makedirs(os.path.dirname(path)) # ensure output directory exists
601                imsave(path, data.T, **kwds)  # save the image
602
603            ax.cbar = ax.imshow(data.T, interpolation=kwds.pop('interpolation', 'none'), **kwds) # change default interpolation to None
604
605        #map 3 bands to RGB
606        elif isinstance(band, tuple) or isinstance(band, list):
607            #get band indices and range
608            rgb = []
609            for b in band:
610                if isinstance(b, str):
611                    rgb.append(self.get_band_index(b))
612                else:
613                    rgb.append(self.get_band_index(np.abs(b)))
614
615            #slice image (as copy) and map to 0 - 1
616            img = np.array(self.data[:, :, rgb]).copy()
617            if np.isnan(img).all():
618                print("Warning - image contains no data.")
619                return ax.get_figure(), ax
620
621            # invert if needed
622            if invert:
623                band = [-b for b in band]
624            for i,b in enumerate(band):
625                if not isinstance(b, str) and (b < 0):
626                    img[..., i] = np.nanmax(img[..., i]) - img[..., i]
627
628            # do scaling
629            if tscale: # scale bands independently
630                for b in range(3):
631                    mn = kwds.get("vmin", float(np.nanmin(img)))
632                    mx = kwds.get("vmax", float(np.nanmax(img)))
633                    if isinstance (mn, int):
634                        assert mn >= 0 and mn <= 100, "Error - integer vmin values must be a percentile."
635                        mn = float(np.nanpercentile(img[...,b], mn ))
636                    if isinstance (mx, int):
637                        assert mx >= 0 and mx <= 100, "Error - integer vmax values must be a percentile."
638                        mx = float(np.nanpercentile(img[...,b], mx ))
639                    img[...,b] = (img[..., b] - mn) / (mx - mn)
640            else: # scale bands together
641                mn = kwds.get("vmin", float(np.nanmin(img)))
642                mx = kwds.get("vmax", float(np.nanmax(img)))
643                if isinstance(mn, int):
644                    assert mn >= 0 and mn <= 100, "Error - integer vmin values must be a percentile."
645                    mn = float(np.nanpercentile(img, mn))
646                if isinstance(mx, int):
647                    assert mx >= 0 and mx <= 100, "Error - integer vmax values must be a percentile."
648                    mx = float(np.nanpercentile(img, mx))
649                img = (img - mn) / (mx - mn)
650
651            #apply brightness/contrast mapping
652            img = np.clip((1.0 + cfac) * img + bfac, 0, 1.0 )
653
654            #apply masking so background is white
655            img[np.logical_not( np.isfinite( img ) )] = 1.0
656            if 'mask' in kwds:
657                img[kwds.pop("mask"),:] = 1.0
658
659            # apply rotations and flipping
660            if rot:
661                img = np.transpose( img, (1,0,2) )
662            if flipX:
663                img = img[::-1, :, :]
664            if flipY:
665                img = img[:, ::-1, :]
666
667            # save?
668            if 'path' in kwds:
669                path = kwds.pop('path')
670                from matplotlib.pyplot import imsave
671                if not os.path.exists(os.path.dirname(path)):
672                    os.makedirs(os.path.dirname(path)) # ensure output directory exists
673                imsave(path, np.transpose( np.clip( img*255, 0, 255).astype(np.uint8), (1, 0, 2)))  # save the image
674
675            # plot samples?
676            ps = kwds.pop('ps', 5)
677            pc = kwds.pop('pc', 'r')
678            if samples:
679                if isinstance(samples, list) or isinstance(samples, np.ndarray):
680                    ax.scatter([s[0] for s in samples], [s[1] for s in samples], s=ps, c=pc)
681                else:
682                    for n in self.header.get_class_names():
683                        points = np.array(self.header.get_sample_points(n))
684                        ax.scatter(points[:, 0], points[:, 1], s=ps, c=pc)
685
686            #plot
687            ax.imshow(np.transpose(img, (1,0,2)), interpolation=kwds.pop('interpolation', 'none'), **kwds)
688            ax.cbar = None  # no colorbar
689
690        return ax.get_figure(), ax
691
692    def createGIF(self, path, bands=None, figsize=(10,10), fps=10, **kwds):
693        """
694        Create and save an animated gif that loops through the bands of the image.
695
696        Args:
697            path (str): the path to save the .gif
698            bands (tuple): Tuple containing the range of band indices to draw. Default is the whole range.
699            figsize (tuple): the size of the image to draw. Default is (10,10).
700            fps (int): the framerate (frames per second) of the gif. Default is 10.
701            **kwds: keywords are passed directly to matplotlib.imshow. Use this to specify cmap etc.
702        """
703
704        frames = []
705        if bands is None:
706            bands = (0,self.band_count())
707        else:
708            assert 0 < bands[0] < self.band_count(), "Error - invalid range."
709            assert 0 < bands[1] < self.band_count(), "Error - invalid range."
710            assert bands[1] > bands[0], "Error - invalid range."
711
712        #plot frames
713        for i in range(bands[0],bands[1]):
714            fig, ax = plt.subplots(figsize=figsize)
715            ax.imshow(self.data[:, :, i], **kwds)
716            fig.canvas.draw()
717            frames.append(np.frombuffer(fig.canvas.tostring_rgb(), dtype='uint8'))
718            frames[-1] = np.reshape(frames[-1], (fig.canvas.get_width_height()[1], fig.canvas.get_width_height()[0], 3))
719            plt.close(fig)
720
721        #save gif
722        imageio.mimsave( os.path.splitext(path)[0] + ".gif", frames, fps=fps)
723
724    ## masking
725    def drop_bbl(self, drop=True):
726        """
727        Remove bad bands as stored in the 'bbl' key in the image header. Note that this operates in-place.
728
729        Args:
730            drop (bool): True if bad bands should be completely dropped. If False, these bands will be kept but replaced with nans.
731        """
732        assert 'bbl' in self.header, "Please specify a bad band list ('bbl') in the image header, as per the ENVI format definition."
733        mask = self.header.get_list('bbl') == 0
734        self.data[...,mask] = np.nan
735        if drop:
736            self.delete_nan_bands(inplace=True)
737    
738    def mask(self, mask=None, flag=np.nan, invert=False, crop=False, bands=None):
739        """
740         Apply a mask to an image, flagging masked pixels with the specified value. Note that this applies the mask to the
741         image in-situ.
742
743         Args:
744            flag (float): the value to use for masked pixels. Default is np.nan
745            mask (ndarray): a numpy array defining the mask polygon of the format [[x1,y1],[x2,y2],...]. If None is passed then
746                    pickPolygon( ... ) is used to interactively define a polygon. If a file path is passed then the polygon
747                    will be loaded using np.load( ... ). Alternatively if mask.shape == image.shape[0,1] then it is treated as a
748                    binary image mask (must be boolean) and True values will be masked across all bands. Default is None.
749            invert (bool): if True, pixels within the polygon will be masked. If False, pixels outside the polygon are masked. Default is False.
750            crop (bool): True if rows/columns containing only zeros should be removed. Default is False.
751            bands (tuple): the bands of the image to plot if no mask is specified. If None, the middle band is used.
752
753         Returns:
754            Tuple containing
755
756            - mask (ndarray): a boolean array with True where pixels are masked and False elsewhere.
757            - poly (ndarray): the mask polygon array in the format described above. Useful if the polygon was interactively defined.
758         """
759
760        if mask is None:  # pick mask interactively
761            if bands is None:
762                bands = int(self.band_count() / 2)
763
764            regions = self.pickPolygons(region_names=["mask"], bands=bands)
765
766            # the user bailed without picking a mask?
767            if len(regions) == 0:
768                print("Warning - no mask picked/applied.")
769                return
770
771            # extract polygon mask
772            mask = regions[0]
773
774        # convert polygon mask to binary mask
775        if mask.shape[1] == 2:
776
777            # build meshgrid with pixel coords
778            xx, yy = np.meshgrid(np.arange(self.xdim()), np.arange(self.ydim()))
779            xx = xx.flatten()
780            yy = yy.flatten()
781            points = np.vstack([xx, yy]).T  # coordinates of each pixel
782
783            # calculate per-pixel mask
784            mask = path.Path(mask).contains_points(points)
785            mask = mask.reshape((self.ydim(), self.xdim())).T
786
787            # flip as we want to mask (==True) outside points (unless invert is true)
788            if not invert:
789                mask = np.logical_not(mask)
790
791        # apply binary image mask
792        assert mask.shape[0] == self.data.shape[0] and mask.shape[1] == self.data.shape[1], \
793            "Error - mask shape %s does not match image shape %s" % (mask.shape, self.data.shape)
794        for b in range(self.band_count()):
795            self.data[:, :, b][mask] = flag
796
797        # crop image
798        if crop:
799            # calculate non-masked pixels
800            valid = np.logical_not(mask)
801
802            # integrate along axes
803            xdata = np.sum(valid, axis=1) > 0.0
804            ydata = np.sum(valid, axis=0) > 0.0
805
806            # calculate domain containing valid pixels
807            xmin = np.argmax(xdata)
808            xmax = xdata.shape[0] - np.argmax(xdata[::-1])
809            ymin = np.argmax(ydata)
810            ymax = ydata.shape[0] - np.argmax(ydata[::-1])
811
812            # crop
813            self.data = self.data[xmin:xmax, ymin:ymax, :]
814
815        return mask
816
817    def crop_to_data(self):
818        """
819        Remove padding of nan or zero pixels from image. Note that this is performed in place.
820        """
821
822        valid = np.isfinite(self.data).any(axis=-1) & (self.data != 0).any(axis=-1)
823        ymin, ymax = np.percentile(np.argwhere(np.sum(valid, axis=0) != 0), (0, 100))
824        xmin, xmax = np.percentile(np.argwhere(np.sum(valid, axis=1) != 0), (0, 100))
825        self.data = self.data[int(xmin):int(xmax), int(ymin):int(ymax), :]  # do clipping
826
827    ##################################################
828    ## Interactive tools for picking regions/pixels
829    ##################################################
830    def pickPolygons(self, region_names, bands=0):
831        """
832        Creates a matplotlib gui for selecting polygon regions in an image.
833
834        Args:
835            names (list, str): a list containing the names of the regions to pick. If a string is passed only one name is used.
836            bands (tuple): the bands of the image to plot.
837        """
838
839        if isinstance(region_names, str):
840            region_names = [region_names]
841
842        assert isinstance(region_names, list), "Error - names must be a list or a string."
843
844        # set matplotlib backend
845        backend = matplotlib.get_backend()
846        matplotlib.use('Qt5Agg')  # need this backend for ROIPoly to work
847
848        # plot image and extract roi's
849        fig, ax = self.quick_plot(bands)
850        roi = MultiRoi(roi_names=region_names)
851        plt.close(fig)  # close figure
852
853        # extract regions
854        regions = []
855        for name, r in roi.rois.items():
856            # store region
857            x = r.x
858            y = r.y
859            regions.append(np.vstack([x, y]).T)
860
861        # restore matplotlib backend (if possible)
862        try:
863            matplotlib.use(backend)
864        except:
865            print("Warning: could not reset matplotlib backend. Plots will remain interactive...")
866            pass
867
868        return regions
869
870    def pickPoints(self, n=-1, bands=hylite.RGB, integer=True, title="Pick Points", **kwds):
871        """
872        Creates a matplotlib gui for picking pixels from an image.
873
874        Args:
875            n (int): the number of pixels to pick, or -1 if the user can select as many as they wish. Default is -1.
876            bands (tuple): the bands of the image to plot. Default is HyImage.RGB
877            integer (bool): True if points coordinates should be cast to integers (for use as indices). Default is True.
878            title (str): The title of the point picking window.
879            **kwds: Keywords are passed to HyImage.quick_plot( ... ).
880
881        Returns:
882            A list containing the picked point coordinates [ (x1,y1), (x2,y2), ... ].
883        """
884
885        # set matplotlib backend
886        backend = matplotlib.get_backend()
887        matplotlib.use('Qt5Agg')  # need this backend for ROIPoly to work
888
889        # create figure
890        fig, ax = self.quick_plot( bands, **kwds )
891        ax.set_title(title)
892
893        # get points
894        points = fig.ginput( n )
895
896        if integer:
897            points = [ (int(p[0]), int(p[1])) for p in points ]
898
899        # restore matplotlib backend (if possible)
900        try:
901            matplotlib.use(backend)
902        except:
903            print("Warning: could not reset matplotlib backend. Plots will remain interactive...")
904            pass
905
906        return points
907
908    def pickSamples(self, names=None, store=True, **kwds):
909        """
910        Pick sample probe points and store these in the image header file.
911
912        Args:
913            names (str, list): the name of the sample to pick, or a list of names to pick multiple.
914            store (bool): True if sample should be stored in the image header file (for later access). Default is True.
915            **kwds: Keywords are passed to HyImage.quick_plot( ... )
916
917        Returns:
918            a list containing a list of points for each sample.
919        """
920
921        if isinstance(names, str):
922            names = [names]
923
924        # pick points
925        points = []
926        for s in names:
927            pnts = self.pickPoints(title="%s" % s, **kwds)
928            if store:
929                self.header['sample %s' % s] = pnts # store in header
930            points.append(pnts)
931        # add class to header file
932        if store:
933            cls_names = self.header.get_class_names()
934            if cls_names is None:
935                cls_names = []
936            self.header['class names'] = cls_names + names
937
938        return points