Commit 1698d8fb authored by Niklas Bohn's avatar Niklas Bohn
Browse files

Added some functions to replace import of isofit library.

parent bdf820ba
...@@ -28,21 +28,139 @@ import numpy as np ...@@ -28,21 +28,139 @@ import numpy as np
from numpy import logical_and as aand from numpy import logical_and as aand
import dill import dill
from tqdm import tqdm from tqdm import tqdm
from os.path import split, abspath import os
from os.path import split, abspath, expandvars
import scipy as s import scipy as s
from scipy.interpolate import interp1d from scipy.interpolate import interp1d
from scipy.io import loadmat from scipy.io import loadmat
from scipy.linalg import norm from scipy.linalg import norm
from spectral.io import envi from spectral.io import envi
from sklearn.cluster import KMeans from sklearn.cluster import KMeans
import json
from isofit.core.common import json_load_ascii, expand_path, svd_inv
from sicor.Tools.EnMAP.LUT import read_lut_enmap_formatted, interpol_lut_c, get_data_file from sicor.Tools.EnMAP.LUT import read_lut_enmap_formatted, interpol_lut_c, get_data_file
from sicor.Tools.EnMAP.conversion import generate_filter from sicor.Tools.EnMAP.conversion import generate_filter
from sicor.Tools.EnMAP.metadata import varsol from sicor.Tools.EnMAP.metadata import varsol
def recursive_reencode(j, shell_replace=True):
"""
Recursively re-encode a mutable object (ascii->str).
:param j: object to reencode
:param shell_replace: boolean helper for recursive calls
:return: expanded, reencoded object
"""
if isinstance(j, dict):
for key, value in j.items():
j[key] = recursive_reencode(value)
return j
elif isinstance(j, list):
for i, k in enumerate(j):
j[i] = recursive_reencode(k)
return j
elif isinstance(j, tuple):
return tuple([recursive_reencode(k) for k in j])
else:
if shell_replace and isinstance(j, str):
try:
j = expandvars(j)
except IndexError:
pass
return j
def json_load_ascii(filename, shell_replace=True):
"""
Load a hierarchical structure, convert all unicode to ASCII and expand environment variables.
:param filename: json file to load from
:param shell_replace: boolean
:return: encoded dictionary
"""
with open(filename, 'r') as fin:
j = json.load(fin)
return recursive_reencode(j, shell_replace)
def expand_path(directory, subpath):
"""
Expand a path variable to an absolute path, if it is not one already.
:param directory: absolute location
:param subpath: path to expand
:return: expanded path
"""
if subpath.startswith('/'):
return subpath
return os.path.join(directory, subpath)
def svd_inv_sqrt(C, hashtable=None):
"""
Matrix inversion, based on decomposition. Built to be stable, and positive.
:param C: matrix to invert
hashtable: if used, the hashtable to store/retrieve results in/from
Return:
(np.array, np.array): inverse of C and square root of the inverse of C
"""
# If we have a hash table, look for the precalculated solution
h = None
if hashtable is not None:
# If arrays are in Fortran ordering, they are not hashable.
if not C.flags['C_CONTIGUOUS']:
C = C.copy(order='C')
h = xxhash.xxh64_digest(C)
if h in hashtable:
return hashtable[h]
D, P = np.linalg.eigh(C)
for count in range(3):
if np.any(D < 0) or np.any(np.isnan(D)):
inv_eps = 1e-6 * (count-1)*10
D, P = np.linalg.eigh(
C + np.diag(np.ones(C.shape[0]) * inv_eps))
else:
break
if count == 2:
raise ValueError('Matrix inversion contains negative values,' +
'even after adding {} to the diagonal.'.format(inv_eps))
Ds = np.diag(1/np.sqrt(D))
L = P@Ds
Cinv_sqrt = L@P.T
Cinv = L@L.T
# If there is a hash table, cache our solution. Bound the total cache
# size by removing any extra items in FIFO order.
if hashtable is not None:
hashtable[h] = (Cinv, Cinv_sqrt)
while len(hashtable) > max_table_size:
hashtable.popitem(last=False)
return Cinv, Cinv_sqrt
def svd_inv(C, hashtable=None):
"""
Matrix inversion, based on decomposition. Built to be stable, and positive.
:param C: matrix to invert
:param hashtable: if used, the hashtable to store/retrieve results in/from
:return: inverse of C
"""
return svd_inv_sqrt(C, hashtable)[0]
class MultiComponentSurface(object): class MultiComponentSurface(object):
""" """
A model of the surface based on a collection of multivariate Gaussians, with one or more equiprobable components A model of the surface based on a collection of multivariate Gaussians, with one or more equiprobable components
......
...@@ -29,7 +29,8 @@ from scipy.spatial import KDTree ...@@ -29,7 +29,8 @@ from scipy.spatial import KDTree
from scipy.stats import linregress from scipy.stats import linregress
from skimage.segmentation import slic from skimage.segmentation import slic
import numpy as np import numpy as np
from multiprocessing import Pool # from multiprocessing import Pool
from multiprocessing.pool import ThreadPool as Pool
def SLIC_segmentation(data_rad_all, n_pca, segs): def SLIC_segmentation(data_rad_all, n_pca, segs):
......
...@@ -66,7 +66,7 @@ ...@@ -66,7 +66,7 @@
} }
}, },
"surface": { "surface": {
"config": "/Users/niklasbohn/Desktop/SICOR_dev_ISOFIT_approach/sicor/sicor/options/surface_model_options.json" "config": "/Users/niklasbohn/Desktop/GFZ_Tasks/SICOR_dev_ISOFIT_approach/sicor/sicor/options/surface_model_options.json"
}, },
"unknowns": { "unknowns": {
"skyview": { "skyview": {
......
{ {
"output_model_file": "/Users/niklasbohn/Desktop/SICOR_dev_ISOFIT_approach/surface_model.mat", "output_model_file": "/Users/niklasbohn/Desktop/GFZ_Tasks/SICOR_dev_ISOFIT_approach/surface_model.mat",
"wavelength_file": "/Users/niklasbohn/Desktop/SICOR_dev_ISOFIT_approach/EnMAP_wavelengths_fwhm_DLR.txt", "wavelength_file": "/Users/niklasbohn/Desktop/GFZ_Tasks/SICOR_dev_ISOFIT_approach/EnMAP_wavelengths_fwhm_DLR.txt",
"normalize":"Euclidean", "normalize":"Euclidean",
"reference_windows":[[400,1300],[1450,1700],[2100,2450]], "reference_windows":[[400,1300],[1450,1700],[2100,2450]],
"sources": "sources":
...@@ -8,7 +8,7 @@ ...@@ -8,7 +8,7 @@
{ {
"input_spectrum_files": "input_spectrum_files":
[ [
"/Users/niklasbohn/Desktop/SICOR_dev_ISOFIT_approach/surface_model_ucsb" "/Users/niklasbohn/Desktop/GFZ_Tasks/SICOR_dev_ISOFIT_approach/surface_model_ucsb"
], ],
"n_components": 8, "n_components": 8,
"windows": [ "windows": [
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment