From d00fd1bba7a4e527ca7959bb8a125595fe922b5a Mon Sep 17 00:00:00 2001 From: Sam Thiele Date: Mon, 30 Mar 2026 12:07:50 +0200 Subject: [PATCH 1/2] Merge dev with main (#18) * Added absorbance conversion function (#16) * added georeferencing aware crop, and tile / mosiac functions * bugfix matching header files with non-unix paths * minor fixes in projection workflow --- .DS_Store | Bin 10244 -> 14340 bytes .gitignore | 1 + hylite/.DS_Store | Bin 10244 -> 12292 bytes hylite/__init__.py | 7 - hylite/correct/detrend.py | 3 +- hylite/hyimage.py | 302 ++++++++++++++++++++++++++++++----- hylite/hyscene.py | 108 +++++++------ hylite/io/headers.py | 5 +- hylite/project/align.py | 96 +++++++++++ hylite/project/pmap.py | 4 - hylite/transform/__init__.py | 65 ++++++++ setup.py | 2 +- test_data/.DS_Store | Bin 6148 -> 6148 bytes tests/test_corrections.py | 37 +++++ tests/test_hyimage.py | 197 +++++++++++++++++++++-- 15 files changed, 704 insertions(+), 123 deletions(-) create mode 100644 tests/test_corrections.py diff --git a/.DS_Store b/.DS_Store index c1ae4f37befa008425431caa1ae218a42697e1f2..5a519821367c7d71b597b6b23976b2ac1edf4f62 100755 GIT binary patch literal 14340 zcmeI2YitzP6@bqf8$7#=VGPEGSsOOE;bDwfV`I##tY2XBGPMmhkJR1AV5Y2hoSj`8 zf>23cQPQN9()97CByC!$tyI!BeW?v;sQQRnRaGjN7AjRq`bebqN1NskQXW0$&g?R? zECi$yMPkO9d*A^;%M$x>A)c>o~H(sgh>pVu5^RDyKY z(Mi`TKne`Vf(~m)zVjiG3fq^{IS-D(NsoO3Oh`jBOj0`DpE48ilo=;YY6j~`o_$*! zbS28*DU!{uvLM-xi~BK1rs>x*pu2uY82X9t0pA0@2Ye6s9=PK?z;iDh@*G?GsqX>b z1HK0e9^n5EVLDkx9eKu<;$H_PyoBd^Uc!rmu4#cT%EFF}I`WJyg(=wxQD^iqkV&BNYMQk?H%wwfyYAjSWLo6ZwtKF;L#6=9S)1qF-Ss; zWUR=7JoHW$@o8()+0Xc4d@^NN*+}FfN@eA?6vJV(A~=$JV@l#OL`R#z-@m~+FX znTqGHWAOuqK0T&Zcbn;K%t-4QVPGV6QaLhWq?6`!mzhf^Enz!kSIB`tU`(w&b*gn^ zTT{5B?Q~Q4)W)_Jx@?b}J{^!3Hmuz|ICeZWkv1QsM?nWV|8Jqv&KJ^$=&FP(l2#%W zryFIx64fbds@v{&Sxvs*f6Dhe5S$Y_Fftk(Q|C`n=LROThMBew=oyRTm^#lQC-f)G z^uCPF^&<&0m77RgV`^0*mP+JOv8>*iN(~wJ>DFjgKau6DLzbDzb_>N(gKW8MOjUE1 zKB$v{F{SsAnN?lrmwkbObIQW{Wy@Db)^FO@)3cs4Ous{|H!KZ>4s%U+MyE9C7_3+s3JzP;Ej8vSHmwN;hmTVT zs21M>nYOG8hW4joaXqCb19m&@8-k&cDciT|E~jwwmSAW!VNOg?7p*F6a*Cs!q0k4p z<=rtWJ20uI>8_r4n(w8$L)NiaMjx1_D0E8usZ=~MDrckKq>(j6B{^gdb>txbTxT*V zy0NIHwlXvr8kW`H%k*%J^pru7PY>N50CGdL(S9t#vklRRoeR*KVKI&JDD=S*8uc`G z;bHg^oP}?}58wj)3SNLq@E3RqUWUKIzu+}^1Iw`z=b?&KsNo7+g^jox+c1iq*oC`r z9}eOW4&yx-$0X`_6iv+GG@ihd_)&ZaAI3-U)A$8^0-waM;@9vTp2wf!b9e!Njeo$4 z_(yyZU&H_6>-eTrA<0renlCMo7D_eJDydOwl{zH)f5nHCA6AF#d_I-B=&JO%asFRO z_m&zr+V&lN?~QSDNwG)Dg6f4!8dk1Z*S?< zwam^oYU=~CiguT;Gw_q=B$3E<*pdH5jID(^i2on^vM{yht%%Fu? zoWhUb84BEw;U_3;KZB3rXYq4*79Yni<5!#jej1;#gRc;hJM+P4OFUmXB&UWl>6Dq~ zDVu|NR$*GSil#$a0(hoFm!@kP{82uv-?eFQ5lw^bH21Z(w2OIgU8E`8*4kb)5Ar-m z$93mHUJki>ROA(f7CZdn;W#ZsQgk(a{e_RmN1{zX@jc*s!1sXf0pA0+%LBakUrv*f zcvJ43xS6fN{YJeauSf6da}Y1mzSZ~dyqCc5@rzrSA|B^=|4Gl=de-K>1h@8V_csGY zcfa25(P=3%MKe zYplyIrS(fR#ls*1y}4%D5>0dC3Q;Z> zFJGZ4SJ9hGRyQfiRrKbPNQ!w@sK=gdIYM$1o?PQ2{0r`wmzE>u+rswJb@_gkWSC-&z8bsn~h_19RJnjjoylY<% z%b|l9e3QnL3CmEyg;@Dq;2JC(mVCSHn&!^G{=&E@5R}vqR>bRa)v=9Woc$rn(J+%tdCCvc(~Dk*!yv0mgF9)OFYO=CM5gpsy3r zH;m^O;4-`jufpqa1!v$@la0XwpZ@Lh(yYii| zCB9736SGASTJ4@EXU)hN->I!oxf7a*O^e$?|#lYUFn5S(I+qyJw zI1%6OZE&Cey>aB{d=K~@@IByr;12P?EW4i(-UH~R?{=>9{r?@}4ZoYd2kuM{V0q6_ z55IWk_jsiai3-DXeu&QeLApYn60YWq_!N&a{3%{XUeBqEI`WJyg(=wxQh5hZX delta 368 zcmZoEXbF&DU|?W$DortDU{C-uIe-{M3-C-V6q~3gIoZI3MHI+q-~wVspg03V2}5d9 zd2zwUf1K=-6&UAjX6InxVC0(oK|z{v=VTtmLJi64Y9k{99R*!OW20Ifg=$L>$HLsS zww9AaR9W9TC_XzUH!r_?vaX~&W6$OzirQ?GjU)|U{Q94#b2*;PVna;L<@%_58+*#wzE&IkI6 z8%Vf<2pdrSOMP0dS zPGaK5jmE7K9>9%j6JNruF5UYo@nIzotmQ&nAhb<1j13}>9#Vt0saJn|(0Rw@Vq!*s1gCP}7rm{my5rUH zt$#Qi4mOFeh{G}6La1lI)sKj8n8!!o@}Kn8 zdFsKgMlN+}2lhU-rwl)i4Tqh~thT}_qk1t(_wnq}0Ag!syO52UNn7m|0{!I$wbjWz1oS=}@0`G#( zozV(Y_`JexB+#RpAVsA(>M^B|9rzNckNPwOJgZ3tQ81}MDMO8HynVqel}8i{vS?8o zSX24RpOXqNMzRQEI#6C$d*jnUITQU8VBoBaqBO99O4x#Aa0Ictw8}rt7rDzD<&w32 zrTnyHJ=s{lR=$`L7o-1{i6PJa;K5E#*T5z%?NEc5Fp_-!S zrqtjl6l}(a5Gc;@v2^G8sx3 z;(<_)psigQzmaw)ttOitz>ek zI{#!P#w(L~1nzH^72LqFS%j&LQIHvAG|)zFAmIwKb+aJTcjn3bI+l~|c_cWPAP&-) QoS-vvvb4YjpiW~Z0M~^kzyJUM diff --git a/hylite/__init__.py b/hylite/__init__.py index 0cb1998..6d94e16 100755 --- a/hylite/__init__.py +++ b/hylite/__init__.py @@ -67,13 +67,6 @@ os.environ["NUMEXPR_NUM_THREADS"] = "1" #disable annoying warnings -#np.warnings.filterwarnings('ignore') -#warnings.filterwarnings("ignore", category=DeprecationWarning) -# ignore all warnings -#def _warn(*args, **kwargs): -# pass -#warnings.warn = _warn - import warnings warnings.filterwarnings('ignore') warnings.filterwarnings("ignore", category=DeprecationWarning) diff --git a/hylite/correct/detrend.py b/hylite/correct/detrend.py index be3b1d7..3439111 100755 --- a/hylite/correct/detrend.py +++ b/hylite/correct/detrend.py @@ -6,6 +6,8 @@ import numpy as np from gfit.util import remove_hull +## TASNIM TO ADD DETRENDING FUNCTIONS HERE + def polynomial(data, degree = 1, method='div'): """ @@ -39,7 +41,6 @@ def polynomial(data, degree = 1, method='div'): return y.reshape(data.shape), t.reshape(data.shape) - def get_hull_corrected(data, band_range=None, method='div', hull='upper', vb=True): """ Apply a hull correction to an entire HyData instance (HyImage, HyCloud or HyLibrary). Returns a corrected copy of diff --git a/hylite/hyimage.py b/hylite/hyimage.py index a79f7f6..fac1214 100755 --- a/hylite/hyimage.py +++ b/hylite/hyimage.py @@ -46,7 +46,8 @@ def __init__(self, data, **kwds): #load any additional project information (specific to images) self.set_projection(kwds.get("projection",None)) - self.affine = kwds.get("affine",[0,1,0,0,0,1]) + self.affine = np.array( kwds.get("affine",[0,1,0,0,0,1]) ) + self.header['affine'] = kwds.get("affine",[0,1,0,0,0,1]) # also store this here # wavelengths if 'wav' in kwds: @@ -93,6 +94,10 @@ def aspx(self): Return the aspect ratio of this image (width/height). """ return self.ydim() / self.xdim() + + ##################################### + ## GEOREFERENCING METHODS + ##################################### def get_extent(self): """ @@ -269,9 +274,221 @@ def world_to_pix(self, x, y, proj = None): #apply return gdal.ApplyGeoTransform(inv, x, y) + def crop(self, xmin, xmax, ymin, ymax, bands=None): + """ + Return a cropped copy of this image. + + Args: + xmin, xmax (int): pixel bounds in x (rows) + ymin, ymax (int): pixel bounds in y (columns) + bands (None, list, tuple): optional band indices or (min,max) range + + Returns: + HyImage: cropped image with updated affine transform + """ + + # ---- validate bounds ---- + xmin = int(max(0, xmin)) + ymin = int(max(0, ymin)) + xmax = int(min(self.xdim(), xmax)) + ymax = int(min(self.ydim(), ymax)) + + assert xmin < xmax and ymin < ymax, "Invalid crop extent." + + # ---- crop data ---- + if bands is None: + data = self.data[xmin:xmax, ymin:ymax, :].copy() + wav = self.get_wavelengths() + else: # band selection + if isinstance(bands, tuple): + b0 = self.get_band_index(bands[0]) + b1 = self.get_band_index(bands[1]) + data = self.data[xmin:xmax, ymin:ymax, b0:b1].copy() + wav = self.get_wavelengths()[b0:b1] + else: + idx = [self.get_band_index(b) for b in bands] + data = self.data[xmin:xmax, ymin:ymax, idx].copy() + wav = self.get_wavelengths()[idx] + + # ---- update affine transform ---- + if self.affine is not None: + a = list(self.affine) + new_affine = a.copy() + + # shift origin to new top-left pixel + new_affine[0] = a[0] + xmin*a[1] + ymin*a[2] + new_affine[3] = a[3] + xmin*a[4] + ymin*a[5] + else: + new_affine = None + + # ---- construct output image ---- + out = HyImage( + data, + header=self.header.copy(), + projection=self.projection, + affine=new_affine, + wav=wav + ) + + return out + + def resize(self, newdims: tuple, interpolation: int = 1): + """ + Resize this image with opencv and update affine transform accordingly. + + Args: + newdims (tuple): the new image dimensions (xdim, ydim) + interpolation (int): opencv interpolation method. Default is cv2.INTER_LINEAR. + """ + import cv2 # avoid import issues if opencv is missing + + old_x, old_y = self.xdim(), self.ydim() + new_x, new_y = int(newdims[0]), int(newdims[1]) + + assert new_x > 0 and new_y > 0, "Invalid resize dimensions." + + # resize data (opencv uses width, height = y, x) + self.data = cv2.resize( + self.data, + (new_y, new_x), + interpolation=interpolation + ) + + # update affine transform + if self.affine is not None: + a = list(self.affine) + + sx = old_x / new_x + sy = old_y / new_y + + self.affine = [ + a[0], # x origin unchanged + a[1] * sx, # pixel width + a[2] * sy, # row rotation + a[3], # y origin unchanged + a[4] * sx, # column rotation + a[5] * sy # pixel height + ] + + def tile(self, tile_size): + """ + Break image into tiles of given size and return a list of HyImage tiles. + Each tile has an updated affine transform reflecting its position in the original image. + + Args: + tile_size (tuple): (tile_x, tile_y) in pixels + Returns: + list of HyImage + """ + tiles = [] + tx, ty = tile_size + nx, ny = self.xdim(), self.ydim() + for i in range(0, nx, tx): + for j in range(0, ny, ty): + data_tile = self.data[i:min(i+tx, nx), j:min(j+ty, ny), :].copy() # slice data + + # compute new affine + # affine = [x0, dx, rotx, y0, roty, dy] + x0, dx, rotx, y0, roty, dy = self.affine + new_x0 = x0 + i*dx + j*rotx + new_y0 = y0 + i*roty + j*dy + new_affine = [new_x0, dx, rotx, new_y0, roty, dy] + + # create tile + tile_img = HyImage( + data_tile, + affine=new_affine, + projection=self.projection, + wav=self.get_wavelengths(), + header=self.header.copy() + ) + tile_img.header['xleft'] = i + tile_img.header['ytop'] = j + tiles.append(tile_img) + + return tiles + + @staticmethod + def mosaic( + tiles, + blend="mean", + resampling="nearest", + out_affine=None, + out_shape=None, + ): + """ + Mosaic georeferenced HyImage tiles using GDAL. Note that this assumes all tiles are in the same coordinate system. + + Args: + tiles (list[HyImage]) + blend (str): 'first', 'min', 'max', 'mean', 'median' + resampling (str): 'nearest', 'bilinear', 'cubic' + out_affine (list): optional 6-element affine to define output grid. If None, the affine of the first tile is used. + out_shape (tuple): optional (xdim, ydim) shape of the output grid. If None, the extent of all tiles will be used. + Returns: + HyImage + """ + import numpy as np + from hylite.project.align import resample_raster + from osgeo import gdal, osr + + assert len(tiles) > 0 + assert blend in ("first", "min", "max", "mean", "median") + + # compute bounds in world coordinates + # N.B. THIS ASSUMES ALL DATA ARE IN THE SAME CRS + points = [] + for t in tiles: + points.append( t.pix_to_world(0,0) ) + points.append( t.pix_to_world(t.xdim()+1,t.ydim()+1) ) + min_x, min_y = np.min(points, axis=0) + max_x, max_y = np.max(points, axis=0) + if out_shape is None: + 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) + if out_affine is None: + out_affine = list(tiles[0].affine) + out_affine[0] = min_x + out_affine[3] = max_y + + if blend == "first": # fill output, first come, first served. + out = np.full( tuple(out_shape) + (tiles[0].band_count(),), np.nan, dtype=np.float32) + for t in tiles: + r = resample_raster( t.data, t.affine, out_affine, out_shape ) + mask = np.isnan(out) + out[mask] = r[mask] + elif blend == "min": # keep minimum value in case of overlap + out = np.nanmin( np.stack([ resample_raster( t.data, t.affine, out_affine, out_shape ) for t in tiles ], + axis=0), axis=0 ) + elif blend == "max": # keep maximum value in case of overlap + out = np.nanmax( np.stack([ resample_raster( t.data, t.affine, out_affine, out_shape ) for t in tiles ], + axis=0), axis=0 ) + elif blend == "mean": # use average in case of overlap + out = np.nanmean( np.stack([ resample_raster( t.data, t.affine, out_affine, out_shape ) for t in tiles ], + axis=0), axis=0 ) + elif blend == "median": # use mean in case of overlap + out = np.nanmedian( np.stack([ resample_raster( t.data, t.affine, out_affine, out_shape ) for t in tiles ], + axis=0), axis=0 ) + + # Return HyImage + out = HyImage( + out, + affine=out_affine, + projection=tiles[0].projection, + wav=tiles[0].get_wavelengths(), + header=tiles[0].header.copy() + ) + if 'xleft' in out.header: del out.header['xleft'] # stored in tiles, but not meaningful here + if 'ytop' in out.header: del out.header['ytop'] # stored in tiles, but not meaningful here + + return out + + ##################################### + ## BASIC TRANSFORMS + ##################################### + def flip(self, axis='x'): """ - Flip the image on the x or y axis. + Flip the image on the x or y axis. Note that this will remove any defined affine transform. Args: axis (str): 'x' or 'y' or both 'xy'. @@ -281,6 +498,9 @@ def flip(self, axis='x'): self.data = np.flip(self.data,axis=0) if 'y' in axis.lower(): self.data = np.flip(self.data,axis=1) + self.affine = None + if 'affine' in self.header: del self.header['affine'] + self.push_to_header() # update width and height info def rot90(self): """ @@ -288,7 +508,9 @@ def rot90(self): to achieve positive/negative rotations. """ self.data = np.transpose( self.data, (1,0,2) ) - self.push_to_header() + self.affine = None + if 'affine' in self.header: del self.header['affine'] + self.push_to_header() # update width and height info ##################################### ##IMAGE FILTERING @@ -347,17 +569,6 @@ def erode(self, size=3, iterations=1): mask = cv2.erode(mask.astype(np.uint8), kernel, iterations=iterations) self.data[mask == 0, :] = 0 - def resize(self, newdims : tuple, interpolation : int = 1): - """ - Resize this image with opencv. - - Args: - newdims (tuple): the new image dimensions. - interpolation (int): opencv interpolation method. Default is cv2.INTER_LINEAR. - """ - import cv2 # import this here to avoid errors if opencv is not installed properly - self.data = cv2.resize(self.data, (newdims[1],newdims[0]), interpolation=interpolation) - def despeckle(self, size=5): """ Despeckle each band of this image (independently) using a median filter. @@ -513,13 +724,13 @@ def match_keypoints(cls, kp1, kp2, d1, d2, method='SIFT', dist=0.7, tree = 5, ch ############################ ## Visualisation methods ############################ - 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, + 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, **kwds): """ Plot a band using matplotlib.imshow(...). Args: - 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 + 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 each band in the tuple (string or index) will be mapped to rgb. Bands with negative wavelengths or indices will be inverted before plotting. ax: an axis object to plot to. If none, plt.imshow( ... ) is used. bfac (float): a brightness adjustment to apply to RGB mappings (-1 to 1) @@ -558,13 +769,13 @@ def quick_plot(self, band=0, ax=None, bfac=0.0, cfac=0.0, samples=False, tscale= ax.set_yticks([]) #map individual band using colourmap - if isinstance(band, str) or isinstance(band, int) or isinstance(band, float): + if isinstance(bands, str) or isinstance(bands, int) or isinstance(bands, float): #get band - if isinstance(band, str): - data = self.data[:, :, self.get_band_index(band)] + if isinstance(bands, str): + data = self.data[:, :, self.get_band_index(bands)] else: - data = self.data[:, :, self.get_band_index(np.abs(band))] - if not isinstance(band, str) and band < 0: + data = self.data[:, :, self.get_band_index(np.abs(bands))] + if not isinstance(bands, str) and bands < 0: data = np.nanmax(data) - data # flip # convert integer vmin and vmax values to percentiles @@ -603,10 +814,10 @@ def quick_plot(self, band=0, ax=None, bfac=0.0, cfac=0.0, samples=False, tscale= ax.cbar = ax.imshow(data.T, interpolation=kwds.pop('interpolation', 'none'), **kwds) # change default interpolation to None #map 3 bands to RGB - elif isinstance(band, tuple) or isinstance(band, list): + elif isinstance(bands, tuple) or isinstance(bands, list): #get band indices and range rgb = [] - for b in band: + for b in bands: if isinstance(b, str): rgb.append(self.get_band_index(b)) else: @@ -620,8 +831,8 @@ def quick_plot(self, band=0, ax=None, bfac=0.0, cfac=0.0, samples=False, tscale= # invert if needed if invert: - band = [-b for b in band] - for i,b in enumerate(band): + bands = [-b for b in bands] + for i,b in enumerate(bands): if not isinstance(b, str) and (b < 0): img[..., i] = np.nanmax(img[..., i]) - img[..., i] @@ -796,21 +1007,7 @@ def mask(self, mask=None, flag=np.nan, invert=False, crop=False, bands=None): # crop image if crop: - # calculate non-masked pixels - valid = np.logical_not(mask) - - # integrate along axes - xdata = np.sum(valid, axis=1) > 0.0 - ydata = np.sum(valid, axis=0) > 0.0 - - # calculate domain containing valid pixels - xmin = np.argmax(xdata) - xmax = xdata.shape[0] - np.argmax(xdata[::-1]) - ymin = np.argmax(ydata) - ymax = ydata.shape[0] - np.argmax(ydata[::-1]) - - # crop - self.data = self.data[xmin:xmax, ymin:ymax, :] + self.crop_to_data() return mask @@ -818,12 +1015,29 @@ def crop_to_data(self): """ Remove padding of nan or zero pixels from image. Note that this is performed in place. """ - valid = np.isfinite(self.data).any(axis=-1) & (self.data != 0).any(axis=-1) - ymin, ymax = np.percentile(np.argwhere(np.sum(valid, axis=0) != 0), (0, 100)) - xmin, xmax = np.percentile(np.argwhere(np.sum(valid, axis=1) != 0), (0, 100)) - self.data = self.data[int(xmin):int(xmax), int(ymin):int(ymax), :] # do clipping + # integrate along axes + xdata = np.sum(valid, axis=1) > 0.0 + ydata = np.sum(valid, axis=0) > 0.0 + + # calculate domain containing valid pixels + xmin = np.argmax(xdata) + xmax = xdata.shape[0] - np.argmax(xdata[::-1]) + ymin = np.argmax(ydata) + ymax = ydata.shape[0] - np.argmax(ydata[::-1]) + + # crop + self.data = self.data[xmin:xmax, ymin:ymax, :] + + # shift affine origin to new top-left pixel + if self.affine is not None: + new_affine = list(self.affine) + new_affine[0] = a[0] + xmin*a[1] + ymin*a[2] + new_affine[3] = a[3] + xmin*a[4] + ymin*a[5] + self.affine = np.array(new_affine) + self.header['affine'] = self.affine + ################################################## ## Interactive tools for picking regions/pixels ################################################## diff --git a/hylite/hyscene.py b/hylite/hyscene.py index aa8c87d..1a61105 100755 --- a/hylite/hyscene.py +++ b/hylite/hyscene.py @@ -8,6 +8,7 @@ from tqdm import tqdm from scipy.ndimage import grey_dilation import numpy as np +from hylite.correct.equalize import hist_eq, norm_eq class HyScene( HyCollection ): """ @@ -98,6 +99,9 @@ def construct(self, image, cloud, camera, s=1, occ_tol = 10, maxf=0, bf=True, vb prg.update(1) self.pmap.filter_footprint( self.maxf ) + # filter nans + self.pmap.remove_nan_pixels() + # do projections if vb: prg.set_description("Converting geometry") @@ -287,58 +291,58 @@ def push_to_image(self, bands, method='closest', image=None, cloud=None): return push_to_image( self.pmap, bands, method, image=image, cloud=cloud ) - # def match_colour_to(self, reference, uniform=True, method='norm', inplace=True): - # - # """ - # Identifies matching pixels between two hyperspectral scenes and uses them to minimise - # colour differences using a linear model (aka by adjusting the brightness/contrast of this scene - # to match the brightness/contrast of the reference scene). WARNING: by default this modifies this scene's - # image IN PLACE. - # - # Args: - # - reference = the scene to match colours to. - # - uniform = True if a single brightness contrast adjustment is applied to all bands (avoids introducing spectral - # artefacts). If False, different corrections area applied to each band - use with CARE! Default is True. - # - method = The colour matching method to use. Current options are: - # - 'norm' = centre-means and scale to match standard deviation. Only compares points known to match. - # - 'hist' = histogram equalisation. Applies to all pixels in scene - use with care! - # Default is 'norm'. - # - inplace = True if the correction should be applied to self.image in-place. If False, no correction is - # applied, and the correction weights (cfac and mfac) returned for future use. Default is True. - # Returns: - # - The corrected image as a HyImage object. If inplace=True (default) then this will be the same as self.image. - # """ - # - # image = self.image - # if not inplace: - # image = image.copy() - # - # if 'norm' in method.lower(): - # # get matching pixels - # px1, px2 = self.intersect_pixels(reference) - # assert px1.shape[0] > 0, "Error - no overlap between images." - # if px1.shape[0] < 1000: - # print("Warning: images only have %d overlapping pixels," - # " which may result in poor colour matching." % px1.shape[0]) - # - # # extract data to create vector of matching values - # px1 = image.data[px1[:, 0], px1[:, 1], :] - # px2 = reference.image.data[px2[:, 0], px2[:, 1], :] - # - # # apply correction - # image.data = norm_eq( image.data, px1, px2, per_band=not uniform, inplace=True) - # - # elif 'hist' in method.lower(): - # if uniform: # apply to whole dataset - # image.data = hist_eq(image.data, reference.image.data) - # else: # apply per band - # for b in range(self.image.band_count()): - # image.data[:, :, b] = hist_eq(image.data[:, :, b], reference.image.data[:, :, b]) - # else: - # assert False, "Error - %s is an unrecognised colour correction method." % method - # - # return image - # + def match_colour_to(self, reference, uniform=True, method='norm', inplace=True): + + """ + Identifies matching pixels between two hyperspectral scenes and uses them to minimise + colour differences using a linear model (aka by adjusting the brightness/contrast of this scene + to match the brightness/contrast of the reference scene). WARNING: by default this modifies this scene's + image IN PLACE. + + Args: + - reference = the scene to match colours to. + - uniform = True if a single brightness contrast adjustment is applied to all bands (avoids introducing spectral + artefacts). If False, different corrections area applied to each band - use with CARE! Default is True. + - method = The colour matching method to use. Current options are: + - 'norm' = centre-means and scale to match standard deviation. Only compares points known to match. + - 'hist' = histogram equalisation. Applies to all pixels in scene - use with care! + Default is 'norm'. + - inplace = True if the correction should be applied to self.image in-place. If False, no correction is + applied, and the correction weights (cfac and mfac) returned for future use. Default is True. + Returns: + - The corrected image as a HyImage object. If inplace=True (default) then this will be the same as self.image. + """ + + image = self.image + if not inplace: + image = image.copy() + + if 'norm' in method.lower(): + # get matching pixels + px1, px2 = self.pmap.intersect_pixels(reference.pmap) + assert px1.shape[0] > 0, "Error - no overlap between images." + if px1.shape[0] < 1000: + print("Warning: images only have %d overlapping pixels," + " which may result in poor colour matching." % px1.shape[0]) + + # extract data to create vector of matching values + px1 = image.data[px1[:, 0], px1[:, 1], :] + px2 = reference.image.data[px2[:, 0], px2[:, 1], :] + + # apply correction + image.data = norm_eq( image.data, px1, px2, per_band=not uniform, inplace=True) + + elif 'hist' in method.lower(): + if uniform: # apply to whole dataset + image.data = hist_eq(image.data, reference.image.data) + else: # apply per band + for b in range(self.image.band_count()): + image.data[:, :, b] = hist_eq(image.data[:, :, b], reference.image.data[:, :, b]) + else: + assert False, "Error - %s is an unrecognised colour correction method." % method + + return image + # ################################### # ##PLOTTING AND EXPORT FUNCTIONS # ################################### diff --git a/hylite/io/headers.py b/hylite/io/headers.py index d325570..0dbb85b 100755 --- a/hylite/io/headers.py +++ b/hylite/io/headers.py @@ -169,13 +169,14 @@ def matchHeader(path): - header = file path to the associated .hdr or .HDR file (if found, otherwise None) - image = file path to the associated image data (if found, otherwise None). """ + path = os.path.normpath(path) # ensure sensibly formatted path # find files with the same name but different extensions path, ext = os.path.splitext(path) header = None image = None - match = glob.glob(path + "*") - + match = [os.path.normpath(p) for p in glob.glob(path + "*")] + assert (path + ext) in match, "Error - file not found (%s)" % (path + ext) match.remove(path + ext) # remove self-match diff --git a/hylite/project/align.py b/hylite/project/align.py index d15b3f4..e683470 100755 --- a/hylite/project/align.py +++ b/hylite/project/align.py @@ -96,6 +96,102 @@ def align_to_cloud_manual( cloud, cam, points, pixels, **kwds ): return est, err +def resample_raster(data, src_gt, dst_gt, dst_shape, order=1): + """ + Resample a raster using GDAL Warp. + + Parameters + ---------- + data : np.ndarray | hylite.HyImage + Source raster array, shape (x, y, band) or (x, y) + src_gt : tuple + Source GDAL affine (x0, px_w, rot_x, y0, rot_y, px_h) + dst_gt : tuple + Destination GDAL affine (x0, px_w, rot_x, y0, rot_y, px_h) + dst_shape : tuple + Output shape (height, width) = (y, x) + order : int + Resampling order (0=nearest, 1=bilinear, 3=cubic) + + Returns + ------- + np.ndarray + Resampled array, same shape as input ([x, y, band] or [x, y]) + """ + + try: + from osgeo import gdal, osr + except: + assert False, "Please install GDAL to use this functionality." + + header = None + if isinstance(data, HyImage): + header = data.header.copy() + data = data.data # extract numpy array + + # check for bad dims + single_band = False + if data.ndim == 2: + data = data[:, :, np.newaxis] + single_band = True + + # get arguments and ensure they are integers + width, height, bands = data.shape + dst_width = int(dst_shape[1]) + dst_height = int(dst_shape[0]) + bands = int(bands) + width = int(width) + height = int(height) + + # Determine resampling type + resample_map = { + 0: gdal.GRA_NearestNeighbour, + 1: gdal.GRA_Bilinear, + 3: gdal.GRA_Cubic, + } + resample_alg = resample_map.get(order, gdal.GRA_Bilinear) + + # reate source MEM dataset + mem_driver = gdal.GetDriverByName("MEM") + src_ds = mem_driver.Create("", width, height, bands, gdal.GDT_Float32) + + src_ds.SetGeoTransform(src_gt) + + # Write data band-by-band + for b in range(bands): + band = src_ds.GetRasterBand(b + 1) + band.WriteArray(data[:, :, b].T) + band.SetNoDataValue(np.nan) + + # Create destination MEM dataset for output + dst_ds = mem_driver.Create("", dst_width, dst_height, bands, gdal.GDT_Float32) + dst_ds.SetGeoTransform(dst_gt) + for b in range(bands): + band = dst_ds.GetRasterBand(b + 1) + band.SetNoDataValue(np.nan) + band.WriteArray(np.full((dst_height, dst_width), np.nan, dtype=np.float32)) # Explicitly initialize to NaN + + # Perform warp (resampling) + gdal.Warp( + destNameOrDestDS=dst_ds, + srcDSOrSrcDSTab=src_ds, + resampleAlg=resample_alg, + srcNodata=np.nan, + dstNodata=np.nan, + warpOptions=["INIT_DEST=NO_DATA"] + ) + + # Read back into NumPy + out = np.full((dst_width, dst_height, bands), np.nan, dtype=np.float32) + for b in range(bands): + out[:, :, b] = dst_ds.GetRasterBand(b + 1).ReadAsArray().T + + if single_band: + return out[:, :, 0] + + if header is not None: + return hylite.HyImage( out, header=header ) + return out def refine_alignment(image, cloud, cam, bands=hylite.RGB, s=2, maxdist=25, histeq=True, method='sift', diff --git a/hylite/project/pmap.py b/hylite/project/pmap.py index 8762d43..b099a35 100755 --- a/hylite/project/pmap.py +++ b/hylite/project/pmap.py @@ -583,7 +583,6 @@ def push_to_cloud(pmap, bands=(0, -1), method='best', image=None, cloud=None ): return out - def push_to_image(pmap, bands='xyz', method='closest', image=None, cloud=None): """ Project the specified data from a point cloud onto an image using a (precalculated) PMap instance. If multiple points map @@ -661,7 +660,6 @@ def push_to_image(pmap, bands='xyz', method='closest', image=None, cloud=None): return out # functions for blending scenes together - def push_geomattr(scene, method='best'): """ Return a HyCloud instance containing the geometric attributes of a given projection. This will have @@ -698,7 +696,6 @@ def push_geomattr(scene, method='best'): attr.set_band_names(bands) return attr - def get_blend_weights(scenes, method, ascloud=True): """ Compute weights for doing scene blending based on geometric attributes computed using @@ -741,7 +738,6 @@ def get_blend_weights(scenes, method, ascloud=True): weights.set_band_names([s.name for s in scenes]) return weights - def blend_scenes(scenes, weights, bands=(0, -1), chunksize=25, trim=False, ooc=True, vb=True): """ Blend together collections of scenes that reference the same underlying cloud to create a fused diff --git a/hylite/transform/__init__.py b/hylite/transform/__init__.py index da3c2b4..c64aeee 100644 --- a/hylite/transform/__init__.py +++ b/hylite/transform/__init__.py @@ -13,6 +13,71 @@ from sklearn.covariance import EmpiricalCovariance import hylite import matplotlib.pyplot as plt +from hylite import HyData, HyImage +import warnings + +def convertToAbsorbance( data : HyData, method: str = 'kubelka-munk', band_range = None) -> HyImage: + """ + Converts Reflectance to Pseudo-Absorbance using methods such as Kubelka-Munk. + Applies to an entire HyData instance (HyImage, HyCloud or HyLibrary). + + Kubelka-Munk transformation: + Based on Escobedo-Morales et al. 2019, 'Automated method for the determination of + the band gap energy of pure and mixed powder samples using diffuse reflectance spectroscopy', + DOI: 10.1016/j.heliyon.2019.e01505 + + This function is developed by Andréa de Lima Ribeiro (Orchid Id: 0000-0003-0096-3627) + + Args: + data: a numpy array or HyData instance to detrend. + method: string specifying the conversion method. Currently only 'kubelka-munk' is supported. + band_range: Tuple containing the (min,max) band indices or wavelengths to run the correction between. If None + (default) then the correction is run of the entire range. Only works if data is a HyData instance. + + Returns: + HyImage: Pseudo-absorbance dataset + """ + + if isinstance(data, HyData): + # create copy containing the bands of interest + if band_range is None: + band_range = (0, -1) + else: + band_range = (data.get_band_index(band_range[0]), data.get_band_index(band_range[1])) + corrected = data.export_bands(band_range) + + # selected range check + arr = corrected.data + valid = arr[np.isfinite(arr)] + valid = valid[ valid >= 0.0 ] # make the negative values as nan + min_val = valid.min() + max_val = valid.max() + + # check if reflectance values are within 0–1 or 0–100 + in_01 = (min_val >= 0.0) and (max_val <= 1.0) + in_0100 = (min_val >= 0.0) and (max_val <= 100.0) + + if not (in_01 or in_0100): + warnings.warn(f"Reflectance values out of expected range: " + f"min={min_val:.3f}, max={max_val:.3f}. Expected 0–1 or 0–100.") + raise ValueError("Invalid reflectance range. Function stopped.") + + # if reflectance values between 0–100, convert to 0–1 for the formula + if in_0100 and not in_01: + arr = arr / 100.0 + + # Kubelka–Munk transformation + if method.lower() == 'kubelka-munk': + eps = 1e-12 #factor added to avoid division by 0 or by negative numbers + arr = np.clip(arr, eps, 1.0) # reflectance should be (0,1] + km = ((1.0 - arr) ** 2) / (2.0 * arr) # Kubelka-Munk transformation + + corrected.data = km + return corrected + + else: + raise TypeError("data must be a HyData instance (HyImage/HyCloud/HyLibrary).") + class NoiseWhitener(BaseEstimator, TransformerMixin): """ diff --git a/setup.py b/setup.py index a4ecca9..a718184 100755 --- a/setup.py +++ b/setup.py @@ -5,7 +5,7 @@ setuptools.setup( name="hylite", - version="1.36", + version="1.37", author="Helmholtz Institute Freiberg", author_email="s.thiele@hzdr.de", description="Open-source toolbox for hyperspectral geology.", diff --git a/test_data/.DS_Store b/test_data/.DS_Store index b46a92244d202e1fa3d231f003aeaf73059b81cb..74adfc331e2f1711900ab6c08550d58818170455 100755 GIT binary patch delta 232 zcmZoMXfc=|#>B`mu~2NHo)!;7CPPwEd2vBfPJYtF!sU|XHE8B!R0kR=%y z7QB)qu~2NHo)%A1d2vBfPJR*t1H;M9j7-ZJCmV>cxUevkFr+dRGn61p zCl%ym7MB+Tu_PmDvL~ZY8c;>N zfJAk*v9YC&f~8TdjzYDefuV(tg1Lc7Z7nBM8XF_geRdr2m-Ate(fPj$^LNo9~ zX&5yN$bh@7EVw8yCqFM8D8>kMI1}q;b`E|HpvQom@640=MJzdh=75db93ZlV8358c BH}wDj diff --git a/tests/test_corrections.py b/tests/test_corrections.py new file mode 100644 index 0000000..b1cb497 --- /dev/null +++ b/tests/test_corrections.py @@ -0,0 +1,37 @@ +import unittest +import os +from pathlib import Path +from hylite import io +import numpy as np + +# TASNIM TO IMPLEMENT THIS +class MyTestCase(unittest.TestCase): + def test_baseline(self): + #image = io.load(os.path.join(os.path.join(str(Path(__file__).parent.parent), "test_data"), "image.hdr")) + #cloud = io.load(os.path.join(os.path.join(str(Path(__file__).parent.parent), "test_data"), "hypercloud.hdr")) + + #self.assertEqual(pca.band_count(), 10) + #self.assertLess( np.nanmax( np.abs( pca.data - pca2.data ) ), 1e-4 ) + + #self.assertLess( np.nanmin(hc.data), 1 ) # assert some values are less than one + #self.assertGreaterEqual( np.nanmin(hc.data), 0 ) # no negatives! + + pass + + def test_absorbance(self): + from hylite.transform import convertToAbsorbance + image = io.load(os.path.join(os.path.join(str(Path(__file__).parent.parent), "test_data"), "image.hdr")) + cloud = io.load(os.path.join(os.path.join(str(Path(__file__).parent.parent), "test_data"), "hypercloud.hdr")) + lib = io.load(os.path.join(os.path.join(str(Path(__file__).parent.parent), "test_data"), "library.csv")) + + for D in [image, cloud, lib]: + absorbance = convertToAbsorbance(D, method='kubelka-munk') + self.assertIsNotNone(absorbance) + self.assertEqual(absorbance.band_count(), D.band_count()) + self.assertEqual(absorbance.data.shape[:-1], D.data.shape[:-1]) + self.assertGreaterEqual(np.nanmin(absorbance.data), 0) + self.assertGreater(np.sum(np.isfinite(absorbance.data)), 0) + + +if __name__ == '__main__': + unittest.main() diff --git a/tests/test_hyimage.py b/tests/test_hyimage.py index f0587ef..e0980a4 100755 --- a/tests/test_hyimage.py +++ b/tests/test_hyimage.py @@ -13,31 +13,118 @@ def test_image(self): # create test image image = genImage(dimx = 1464, dimy=401, nbands=10) + # check basics + self.assertEqual(image.xdim(), 1464) + self.assertEqual(image.ydim(), 401) + self.assertEqual(image.band_count(), 10) + self.assertEqual(image.aspx(), 401 / 1464) + # run plotting functions image.quick_plot( (0,1,2), vmin=2, vmax=98 ) image.quick_plot( 0 ) - self.assertEqual(image.xdim(), 1464) - self.assertEqual(image.ydim(), 401) + # ------------------------------------------------------------------ + # Set a known affine transform for georeferencing tests + # GDAL format: [x0, px_w, rot_x, y0, rot_y, px_h] + # ------------------------------------------------------------------ + + # do we have GDAL? + gdal = True + try: + from osgeo import gdal + gdal = True + except: + gdal = False + + image.affine = np.array( [1000.0, 2.0, 0.0, 2000.0, 0.0, -2.0] ) + if gdal: + image.set_projection_EPSG('EPSG:32633') # set an EPSG code as this is needed for pix to world tests + + # ------------------------------------------------------------------ + # Resize and test affine update + # ------------------------------------------------------------------ + old_affine = image.affine.copy() + old_x, old_y = image.xdim(), image.ydim() + + nx, ny = int(old_x / 2), int(old_y / 2) + image.resize(newdims=(nx, ny)) + + self.assertEqual(image.xdim(), nx) + self.assertEqual(image.ydim(), ny) self.assertEqual(image.band_count(), 10) - self.assertEqual(image.aspx(), 401 / 1464) - # todo - add test code for georeferencing code - # get_extent, set_projection, set_projection_EPSG, get_projection_EPSG, pix_to_world, world_to_pix + # pixel scaling factors + sx = old_x / nx + sy = old_y / ny + + # origin should be unchanged + self.assertAlmostEqual(image.affine[0], old_affine[0]) + self.assertAlmostEqual(image.affine[3], old_affine[3]) + + # pixel size should scale + self.assertAlmostEqual(image.affine[1], old_affine[1] * sx) + self.assertAlmostEqual(image.affine[5], old_affine[5] * sy) + + # rotation terms (if any) should scale consistently + self.assertAlmostEqual(image.affine[2], old_affine[2] * sy) + self.assertAlmostEqual(image.affine[4], old_affine[4] * sx) + + # world-space invariant check (strongest test) + if gdal: + x0, y0 = image.pix_to_world(0, 0) + x1, y1 = image.pix_to_world(nx, ny) + + ox0, oy0 = hylite.HyImage( + np.zeros((old_x, old_y, 1)), + affine=old_affine, + projection=image.projection + ).pix_to_world(0, 0) + + ox1, oy1 = hylite.HyImage( + np.zeros((old_x, old_y, 1)), + affine=old_affine, + projection=image.projection + ).pix_to_world(old_x, old_y) + + self.assertAlmostEqual(x0, ox0) + self.assertAlmostEqual(y0, oy0) + self.assertAlmostEqual(x1, ox1) + self.assertAlmostEqual(y1, oy1) + + # Check crop invariance + xmin, xmax = 100, 300 + ymin, ymax = 50, 200 + + cropped = image.crop(xmin, xmax, ymin, ymax) + + # shape check + self.assertEqual(cropped.xdim(), xmax - xmin) + self.assertEqual(cropped.ydim(), ymax - ymin) + # cropped origin must match original pixel location + cx0, cy0 = cropped.pix_to_world(0, 0) + oxc, oyc = image.pix_to_world(xmin, ymin) + + self.assertAlmostEqual(cx0, oxc) + self.assertAlmostEqual(cy0, oyc) + + # cropped far corner must align too + cx1, cy1 = cropped.pix_to_world( + cropped.xdim(), + cropped.ydim() + ) + ox1, oy1 = image.pix_to_world(xmax, ymax) + + self.assertAlmostEqual(cx1, ox1) + self.assertAlmostEqual(cy1, oy1) + + # test some image manipulations image.flip(axis='y') image.data[10,10,:] = np.nan image.fill_holes() self.assertEqual( np.isfinite( image.data ).all(), True ) image.blur() - # resize - nx,ny = int(1464/2), int(401/2) - image.resize(newdims=(nx, ny)) - self.assertEqual(image.xdim(), nx) - self.assertEqual(image.ydim(), ny) - self.assertEqual(image.band_count(), 10) - # extract features k, d = image.get_keypoints( band=0 ) src, dst = image.match_keypoints(k,k,d,d) @@ -46,6 +133,92 @@ def test_image(self): # masking image.mask( np.sum(image.data,axis=2) > 0.75 ) self.assertEqual(np.isfinite(image.data).all(), False) + + def test_tile_and_mosaic_affine_correctness(self): + import numpy as np + import hylite + + try: + from osgeo import gdal + except ImportError: + self.skipTest("GDAL not available") + + # Create a rotated / skewed test image with structured signal + nx, ny, nb = 128, 128, 5 + + x = np.linspace(0, 2 * np.pi, nx) + y = np.linspace(0, 2 * np.pi, ny) + xx, yy = np.meshgrid(x, y, indexing="ij") + + data = np.zeros((nx, ny, nb), dtype=np.float32) + + for b in range(nb): + # long-wavelength spatial signal + band-specific phase + data[..., b] = ( + np.sin(xx * 0.5 + b * 0.3) + + np.cos(yy * 0.5 - b * 0.2) + ) + + data[::32, ::32, :] = 1.0 # add some spikes + # normalize to [0,1] for numerical stability + data -= data.min() + data /= data.max() + + affine = [ + 1000.0, # x origin + 2.0, # pixel width + 0.5, # x skew + 2000.0, # y origin + -0.3, # y skew + -2.0 # pixel height (north-up) + ] + + img = hylite.HyImage( + data, + affine=affine, + wav=np.arange(nb) + ) + img.set_projection_EPSG('EPSG:32633') # UTM zone 33N + + # Tile image + tile_size = (32, 32) + tiles = img.tile(tile_size) + + self.assertGreater(len(tiles), 0) + + # 1. Verify tile affine correctness + for t in tiles: + # Pick tile origin pixel in original image + # Compute expected world coordinate + tx0, ty0 = t.pix_to_world(0, 0) + px0, py0 = img.pix_to_world(t.header['xleft'], t.header['ytop']) + self.assertAlmostEqual( tx0, px0, places=5 ) + self.assertAlmostEqual( ty0, py0, places=5 ) + + # Recover pixel offset by inverse mapping + px, py = img.world_to_pix(tx0, ty0) + self.assertAlmostEqual(px, t.header['xleft'], places=5) + self.assertAlmostEqual(py, t.header['ytop'], places=5) + + # World-space invariant + ox, oy = img.pix_to_world(int(round(px)), int(round(py))) + self.assertAlmostEqual(tx0, ox, places=6) + self.assertAlmostEqual(ty0, oy, places=6) + + # 2. Mosaic tiles back into a NEW grid + for b in ['first', 'mean', 'max', 'min', 'median']: + mosaic = hylite.HyImage.mosaic(tiles, blend=b, out_shape=[128,128]) + + self.assertAlmostEqual(mosaic.affine[0], img.affine[0]) + self.assertAlmostEqual(mosaic.affine[1], img.affine[1]) + self.assertAlmostEqual(mosaic.affine[2], img.affine[2]) + self.assertAlmostEqual(mosaic.affine[3], img.affine[3]) + self.assertAlmostEqual(mosaic.affine[4], img.affine[4]) + self.assertAlmostEqual(mosaic.affine[5], img.affine[5]) + + diff = np.abs(mosaic.data - data) + mask = np.isfinite(diff) + self.assertLess(np.nanmean(diff[mask]), 1e-3) if __name__ == '__main__': unittest.main() From 3122048ba28bff308213ee4080f894e681d52d71 Mon Sep 17 00:00:00 2001 From: Cevedal <69959754+sandralorenz268@users.noreply.github.com> Date: Wed, 22 Apr 2026 09:19:20 +0200 Subject: [PATCH 2/2] added functions to mask and plot spectra in hyperclouds from rendering (#19) --- hylite/hycloud.py | 326 +++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 325 insertions(+), 1 deletion(-) diff --git a/hylite/hycloud.py b/hylite/hycloud.py index 784efec..7d28e47 100755 --- a/hylite/hycloud.py +++ b/hylite/hycloud.py @@ -13,6 +13,10 @@ from hylite.hyimage import HyImage from hylite.project import proj_persp, proj_pano, rasterize, Camera +import matplotlib.pyplot as plt +from matplotlib import path +from roipoly import MultiRoi + class HyCloud( HyData ): @@ -592,7 +596,327 @@ def quick_plot(self, band='rgb', cam='ortho', s=1, step=1, fill_holes=False, blu else: return img.quick_plot((0, 1, 2), **kwds) else: - return img.quick_plot(0, **kwds) + return img.quick_plot(0, **kwds) + + def mask(self, mask=None, flag=np.nan, invert=False, bands=['rgb'], cam='ortho', **kwds): + """ + Apply a mask to a point cloud by rendering and polygon selection. + Masks ALL points projecting into the polygon, even if multiple points + project to the same pixel. + + Args: + mask (ndarray, str, None): Polygon [[x1,y1],[x2,y2],...] or None for interactive + flag (float): Value for masked points (default np.nan) or 'delete' to remove + invert (bool): If True, mask inside polygon; if False, mask outside + bands (list, str, tuple): Display bands for rendering + cam (str or Camera): Camera for rendering ('ortho' default) + **kwds: Additional render args (s, fill_holes, blur, etc.) + + Returns: + Tuple of (point_mask, polygon) + """ + + from hylite.project import proj_persp, proj_pano + + if mask is None: + # Interactive polygon picking + print("Rendering cloud for visualization...") + rendered_img = self.render(cam=cam, bands=bands, **kwds) + + n_display_bands = rendered_img.band_count() + if n_display_bands >= 3: + display_bands = list(range(n_display_bands))[:3] + else: + display_bands = 0 + + print("\nSelect masking polygon. Close window when done.") + print("Draw polygon by clicking points, right-click to finalize.") + + backend = plt.matplotlib.get_backend() + plt.matplotlib.use('Qt5Agg') + + try: + fig, ax = rendered_img.quick_plot(bands=display_bands, vmax=98, vmin=2) + ax.set_title("Draw polygon for masking (right-click to finish)") + + roi = MultiRoi(roi_names=["mask"]) + plt.close(fig) + + if len(roi.rois) == 0: + print("Warning - no mask picked/applied.") + return None, None + + regions = [] + for name, r in roi.rois.items(): + x = r.x + y = r.y + regions.append(np.vstack([x, y]).T) + + mask = regions[0] + + finally: + try: + plt.matplotlib.use(backend) + except: + print("Warning: could not reset matplotlib backend.") + + elif isinstance(mask, str): + mask = np.load(mask) + + assert isinstance(mask, np.ndarray) and mask.shape[1] == 2, \ + "Error - mask must be (N,2) polygon vertices" + + print("Projecting all points to determine masking...") + + # Get camera and projection parameters from the visualization render + if isinstance(cam, str) and 'ortho' in cam: + # Orthographic projection - use render to get bounds + dummy = self.render(cam=cam, bands=['rgb'], **{k:v for k,v in kwds.items() if k != 'step'}) + + # Get world bounds from render + cloudxmin = np.amin(self.xyz[:, 0]) + cloudxmax = np.amax(self.xyz[:, 0]) + cloudymin = np.amin(self.xyz[:, 1]) + cloudymax = np.amax(self.xyz[:, 1]) + + res = kwds.get("res", None) + if res is None: + res = max(cloudxmax - cloudxmin, cloudymax - cloudymin) / 1000 + + # Project ALL points to pixel coordinates + C = np.array([cloudxmin, cloudymax, 0]) + pp = np.abs(self.xyz - C[None, :]) / res + vis = np.ones(self.xyz.shape[0], dtype=bool) + dims = (int((cloudxmax - cloudxmin) / res + 1), int((cloudymax - cloudymin) / res + 1)) + + else: + # Perspective or panoramic projection + if isinstance(cam, str): + raise ValueError("Only 'ortho' camera is supported for string cameras") + + if 'persp' in cam.proj.lower(): + pp, vis = proj_persp(self.xyz, C=cam.pos, a=cam.ori, fov=cam.fov, dims=cam.dims) + elif 'pano' in cam.proj.lower(): + pp, vis = proj_pano(self.xyz, C=cam.pos, a=cam.ori, fov=cam.fov, dims=cam.dims, step=cam.step) + else: + raise ValueError("Unknown camera projection") + + # Extract pixel coordinates (x, y only, ignore depth) + points_2d = pp[:, :2] # Shape: (n_points, 2) + + # Check which projected points fall within the polygon + polygon_mask = path.Path(mask).contains_points(points_2d) + + # Apply inversion if requested + if not invert: + polygon_mask = np.logical_not(polygon_mask) + + # Only consider visible points + point_mask = polygon_mask & vis + + print(f"Masking {np.sum(point_mask)} points out of {self.point_count()}") + + if flag == 'delete': + self.filter_points(band=0, val=point_mask[:, None], trim=False) + print(f"Deleted {np.sum(point_mask)} points. Remaining: {self.point_count()}") + else: + if self.has_bands(): + for b in range(self.band_count()): + self.data[point_mask, b] = flag + + return point_mask, mask + + def plot_from_render(self, click=None, bands='rgb', band_range=None, cam='ortho', radius=5, **kwds): + """ + Render the point cloud, let the user click a pixel, and show the spectral + data next to the rendered image in a single Qt popup window. + Wavelengths are read automatically from self.get_wavelengths(). + + Args: + click (array-like, None): Pixel coordinate [x, y], or None for interactive + bands (list, str): Display bands for rendering (default 'rgb') + band_range (tuple, None): Optional subset of bands to plot. + - (int, int) → interpreted as band index range, e.g. band_range=(10, 50) + - (float, float) → interpreted as wavelength range in nm, e.g. band_range=(400.0, 700.0) + If None, the full spectrum is plotted. + cam (str or Camera): Camera for rendering ('ortho' default) + radius (int): Pixel radius around the clicked point (default 5) + **kwds: Additional render args passed to self.render() + + Returns: + Tuple of (selected_data, click_coord) + selected_data : ndarray of shape (n_selected_points, n_bands_plotted) + click_coord : ndarray [x, y] of the clicked pixel + """ + + from hylite.project import proj_persp, proj_pano + + # ------------------------------------------------------------------ # + # 1. Render image once # + # ------------------------------------------------------------------ # + print("Rendering cloud for visualization...") + rendered_img = self.render(cam=cam, bands=bands, **kwds) + + n_display_bands = rendered_img.band_count() + display_bands = list(range(n_display_bands))[:3] if n_display_bands >= 3 else 0 + + # ------------------------------------------------------------------ # + # 2. Project all points to pixel space # + # ------------------------------------------------------------------ # + print("Projecting points...") + + if isinstance(cam, str) and 'ortho' in cam: + cloudxmin = np.amin(self.xyz[:, 0]) + cloudxmax = np.amax(self.xyz[:, 0]) + cloudymin = np.amin(self.xyz[:, 1]) + cloudymax = np.amax(self.xyz[:, 1]) + + res = kwds.get("res", None) + if res is None: + res = max(cloudxmax - cloudxmin, cloudymax - cloudymin) / 1000 + + C = np.array([cloudxmin, cloudymax, 0]) + pp = np.abs(self.xyz - C[None, :]) / res + vis = np.ones(self.xyz.shape[0], dtype=bool) + + else: + if isinstance(cam, str): + raise ValueError("Only 'ortho' camera is supported as a string; pass a Camera object.") + + if 'persp' in cam.proj.lower(): + pp, vis = proj_persp(self.xyz, C=cam.pos, a=cam.ori, fov=cam.fov, dims=cam.dims) + elif 'pano' in cam.proj.lower(): + pp, vis = proj_pano(self.xyz, C=cam.pos, a=cam.ori, fov=cam.fov, dims=cam.dims, step=cam.step) + else: + raise ValueError("Unknown camera projection type.") + + points_2d = pp[:, :2] + + # ------------------------------------------------------------------ # + # 3. Resolve wavelengths and optional band subset # + # ------------------------------------------------------------------ # + wavelengths = self.get_wavelengths() # None if not set + + n_bands = self.data.shape[1] + wavelengths = self.get_wavelengths() # None if not set + + if wavelengths is not None: + wavelengths = np.asarray(wavelengths) + + if band_range is not None: + b_min, b_max = band_range + if isinstance(b_min, int) and isinstance(b_max, int): + # Integer → treat as band index + band_mask = np.zeros(n_bands, dtype=bool) + band_mask[b_min:b_max + 1] = True + x = wavelengths[band_mask] if wavelengths is not None else np.where(band_mask)[0] + xlabel = "Wavelength (nm)" if wavelengths is not None else "Band index" + else: + # Float → treat as wavelength range + if wavelengths is None: + raise ValueError("band_range given as float (wavelength range) but no wavelengths found in self.get_wavelengths().") + band_mask = (wavelengths >= b_min) & (wavelengths <= b_max) + x = wavelengths[band_mask] + xlabel = "Wavelength (nm)" + else: + band_mask = np.ones(n_bands, dtype=bool) + x = wavelengths if wavelengths is not None else np.arange(n_bands) + xlabel = "Wavelength (nm)" if wavelengths is not None else "Band index" + + # ------------------------------------------------------------------ # + # 4. Setup Qt figure with image + empty spectrum side by side # + # ------------------------------------------------------------------ # + backend = plt.matplotlib.get_backend() + plt.matplotlib.use('Qt5Agg') + + fig, (ax_img, ax_spec) = plt.subplots(1, 2, figsize=(14, 5)) + fig.suptitle("Click on a pixel – spectrum appears on the right", fontsize=11) + + rendered_img.quick_plot(bands=display_bands, vmax=98, vmin=2, ax=ax_img) + ax_img.set_title("Rendered Image") + + ax_spec.set_title("Spectrum") + ax_spec.set_xlabel(xlabel) + ax_spec.set_ylabel("Intensity") + ax_spec.text(0.5, 0.5, "← Click a pixel", + ha='center', va='center', transform=ax_spec.transAxes, + fontsize=12, color='gray') + + plt.tight_layout() + plt.ion() + plt.show() + + # ------------------------------------------------------------------ # + # 5. Click loop – update spectrum on every click, close on Enter # + # ------------------------------------------------------------------ # + last_data, last_click = None, None + + print("Click a pixel to display its spectrum. Press Enter to finish.") + + def on_click(event): + nonlocal last_data, last_click + + if event.inaxes is not ax_img or event.button != 1: + return + + click = np.array([event.xdata, event.ydata]) + + distances = np.linalg.norm(points_2d - click[None, :], axis=1) + in_radius = (distances <= radius) & vis + n_selected = int(np.sum(in_radius)) + + for artist in ax_img.lines + ax_img.patches: + artist.remove() + circle = plt.Circle(click, radius, color='red', fill=False, + linewidth=1.5, linestyle='--') + ax_img.add_patch(circle) + ax_img.plot(*click, 'r+', markersize=10, markeredgewidth=2) + + ax_spec.cla() + ax_spec.set_xlabel(xlabel) + ax_spec.set_ylabel("Intensity") + + if n_selected == 0: + ax_spec.text(0.5, 0.5, f"No points found (r={radius}px)", + ha='center', va='center', transform=ax_spec.transAxes, + fontsize=11, color='gray') + ax_spec.set_title("Spectrum – no points found") + else: + selected_data = self.data[in_radius][:, band_mask] + mean_spectrum = np.nanmean(selected_data, axis=0) + + for row in selected_data: + ax_spec.plot(x, row, color='steelblue', alpha=0.3, linewidth=0.8) + ax_spec.plot(x, mean_spectrum, color='navy', linewidth=2, + label=f'Mean (n={n_selected})') + ax_spec.set_title(f"Pixel ({click[0]:.1f}, {click[1]:.1f}), r={radius}px") + ax_spec.legend() + + last_data, last_click = selected_data, click + print(f"{n_selected} points at pixel ({click[0]:.1f}, {click[1]:.1f})") + + fig.canvas.draw_idle() + + def on_key(event): + if event.key in ('enter', 'escape'): + plt.close(fig) + + fig.canvas.mpl_connect('button_press_event', on_click) + fig.canvas.mpl_connect('key_press_event', on_key) + + plt.ioff() + plt.show(block=True) + + try: + plt.matplotlib.use(backend) + except: + print("Warning: could not reset matplotlib backend.") + + if last_click is None: + print("No pixel selected.") + return None, None + + return last_data, last_click def colourise(self, bands, stretch=(1, 99)):