Commit 0d32f4ff authored by Daniel Scheffler's avatar Daniel Scheffler
Browse files

Moved calc_sam() and calc_sid() to top-level of the module.


Signed-off-by: Daniel Scheffler's avatarDaniel Scheffler <danschef@gfz-potsdam.de>
parent f2787388
Pipeline #3902 failed with stage
in 2 minutes and 10 seconds
......@@ -127,6 +127,7 @@ class _ImageClassifier(object):
# at nodata positions, the distances may have the initialzation value 1e6 (MinDist)
dists = distances[cmap[:] != cmap_nodataVal] if cmap_nodataVal is not None else distances
# noinspection PyTypeChecker
threshold = np.nanpercentile(dists, 100 - percent)
else:
raise ValueError(threshold)
......@@ -292,9 +293,7 @@ class SAM_Classifier(_ImageClassifier):
# loop over all training spectra and compute spectral angle for each pixel
for n_sample in range(self.n_samples):
train_spectrum = train_spectra_norm[n_sample, :].reshape(1, 1, self.n_features)
angles[:, :, n_sample] = self._calc_sam(tileimdata_norm,
train_spectrum,
axis=2)
angles[:, :, n_sample] = calc_sam(tileimdata_norm, train_spectrum, axis=2)
angles_min = np.min(angles, axis=2).astype(np.float32)
cmap = np.argmin(angles, axis=2).astype(np.int16)
......@@ -305,22 +304,6 @@ class SAM_Classifier(_ImageClassifier):
return tilepos, cmap.astype(np.int16), angles_min
@staticmethod
def _calc_sam(s1_norm, s2_norm, axis=0):
"""Compute spectral angle between two vectors or images (in radians)."""
upper = np.sum(s1_norm * s2_norm, axis=axis)
lower = np.sqrt(np.sum(s1_norm * s1_norm, axis=axis)) * np.sqrt(np.sum(s2_norm * s2_norm, axis=axis))
if lower.ndim > 1:
lower[lower == 0] = 1e-10
else:
lower = lower or 1e-10
quotient = upper / lower
quotient[np.isclose(quotient, 1)] = 1 # in case of pixels that are equal to the endmember
return np.arccos(quotient)
def label_unclassified_pixels(self, label_unclassified, threshold):
# type: (int, Union[str, int, float]) -> GeoArray
return self._label_unclassified_pixels(
......@@ -365,9 +348,7 @@ class SID_Classifier(_ImageClassifier):
# loop over all training spectra and compute spectral information divergence for each pixel
for n_sample in range(self.n_samples):
train_spectrum = train_spectra_norm[n_sample, :].reshape(1, 1, self.n_features)
sid[:, :, n_sample] = self._calc_sid(tileimdata_norm,
train_spectrum,
axis=2)
sid[:, :, n_sample] = calc_sid(tileimdata_norm, train_spectrum, axis=2)
sid_min = np.min(sid, axis=2).astype(np.float32)
cmap = np.argmin(sid, axis=2).astype(np.int16)
......@@ -378,24 +359,6 @@ class SID_Classifier(_ImageClassifier):
return tilepos, cmap.astype(np.int16), sid_min
@staticmethod
def _calc_sid(s1_norm, s2_norm, axis=0):
"""Compute the spectral information divergence between two vectors or images."""
def get_sum(x, axis=0):
s = np.sum(x, axis=axis)
s[s == 0] = 1e-10
return s
if s1_norm.ndim == 3 and s2_norm.ndim == 3:
p = (s1_norm / get_sum(s1_norm, axis=axis)[:, :, np.newaxis]) + np.spacing(1)
q = (s2_norm / get_sum(s1_norm, axis=axis)[:, :, np.newaxis]) + np.spacing(1)
else:
p = (s1_norm / get_sum(s1_norm, axis=axis)) + np.spacing(1)
q = (s2_norm / get_sum(s1_norm, axis=axis)) + np.spacing(1)
return np.sum(p * np.log(p / q) + q * np.log(q / p), axis=axis)
def label_unclassified_pixels(self, label_unclassified, threshold):
# type: (int, Union[str, int, float]) -> GeoArray
return self._label_unclassified_pixels(
......@@ -542,6 +505,40 @@ def normalize_endmembers_image(endmembers, image):
return em, im
def calc_sam(s1_norm, s2_norm, axis=0):
"""Compute spectral angle between two vectors or images (in radians)."""
upper = np.sum(s1_norm * s2_norm, axis=axis)
lower = np.sqrt(np.sum(s1_norm * s1_norm, axis=axis)) * np.sqrt(np.sum(s2_norm * s2_norm, axis=axis))
if lower.ndim > 1:
lower[lower == 0] = 1e-10
else:
lower = lower or 1e-10
quotient = upper / lower
quotient[np.isclose(quotient, 1)] = 1 # in case of pixels that are equal to the endmember
return np.arccos(quotient)
def calc_sid(s1_norm, s2_norm, axis=0):
"""Compute the spectral information divergence between two vectors or images."""
def get_sum(x, axis=0):
s = np.sum(x, axis=axis)
s[s == 0] = 1e-10
return s
if s1_norm.ndim == 3 and s2_norm.ndim == 3:
p = (s1_norm / get_sum(s1_norm, axis=axis)[:, :, np.newaxis]) + np.spacing(1)
q = (s2_norm / get_sum(s1_norm, axis=axis)[:, :, np.newaxis]) + np.spacing(1)
else:
p = (s1_norm / get_sum(s1_norm, axis=axis)) + np.spacing(1)
q = (s2_norm / get_sum(s1_norm, axis=axis)) + np.spacing(1)
return np.sum(p * np.log(p / q) + q * np.log(q / p), axis=axis)
def im2spectra(geoArr):
# type: (Union[GeoArray, np.ndarray]) -> np.ndarray
"""Convert 3D images to array of spectra samples (rows: samples; cols: spectral information)."""
......
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