implement convex hull algorithm

[Andreas Janz]

The simple, but slow, solution would be to use the Scipy implementation:
https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.ConvexHull.html

A mutch more performant solution would be to implement it using Numba. Problem is: some users had problems installing Numba on there systems.

Maybe we implement both optiuons, and let the user choose which one to use.

[Andreas Janz]

Stéphane Guillaso already provided the Numba code:

‌

from numba import jit, float32, void, intc, int64, int32
import numpy as np
#   This routine calculate the continuum line to be removed from the original spectrum.
#   It is a recursive routine.
#
#   Algo: 
#   1) calculate the line equation between the extreme points of the spectrum y
#     we obtain 2 coefficients: a and b
#   2) calculate the continumm line h between extreme points
#     h = a * x + b
#   3) calculate distance from each spectrum position to the line dlm
#   4) we choose the lowest value corresponding the the biggest distance between
#      line and a point above the line (this is due to the calculation line-point)
#   5) Test this distance
#     - if dlm >= 0 then we exit the routine, meaning that we don't have any point
#       above the line
#   6) We split the spectrum into 2 parts at the position of the selected point
#      (dlp)
#   7) We call the same routine over the 2 sub lines and we continue.
#
#   This routines is based on the algorithm described here:
#   https://en.wikipedia.org/wiki/Quickhull
#   the main difference is that we calculate only the upper part (as we are not 
#   interesting to the lower part of the spectrum)
# 
#   At the end, the variable h (having the same size than x and y) is returning
#   the continuum line.
@jit(void(float32[:], float32[:], float32[:], intc, intc), nopython=True, fastmath=True)
def get_continuum_line(x, y, h, beg, end):

    a = (y[end] - y[beg]) / (x[end] - x[beg])
    b = y[beg] - a * x[beg]
    dlv = 0.
    dlm = 10000.
    dlp = -1
    for k in range(beg, end+1, 1):
        h[k] = a * x[k] + b
        dlv = (a * x[k] - y[k] + b) / np.sqrt(a*a + 1)
        if (dlv<dlm):
            dlm = dlv
            dlp = k

    if (dlm == 0): return
    if (dlp > beg and dlp < end):
        get_continuum_line(x, y, h, beg, dlp)
        get_continuum_line(x, y, h, dlp, end)

Where x is wavelength value, y is spectrum value and h is the hull function.

Be aware that the function is only calculating the upper part of the hull and not the total (which can easy be implemented by following the code given in https://en.wikipedia.org/wiki/Quickhull

I hope its help

‌

Best

Stéphane

‌

[Andreas Janz]

Hi @Agustin Lobo , I would allow to calculate the hull and the continuum-removed outputs at the same time. Both of them could be skipped. Also, the user can decide to use band number or nanomters as X values.

Any further suggestions from your side?

‌

[Agustin Lobo]

I understand out_contrem = input - out_convhull, correct?

Could I test it?

‌

[Andreas Janz]

Not yet implemented

Isn’t it out_contrem = input / out_convhull ? Like in ENVI: https://www.l3harrisgeospatial.com/docs/continuumremoval.html

[Agustin Lobo]

See Hecker, Christoph, Frank J.A. van Ruitenbeek, Harald M.A. van der Werff, Wim H. Bakker, Robert D. Hewson, and Freek D. van der Meer. 2019. “Spectral Absorption Feature Analysis for Finding Ore: A Tutorial on Using the Method in Geological Remote Sensing.” IEEE Geoscience and Remote Sensing Magazine 7 (2): 51–71. https://doi.org/10.1109/MGRS.2019.2899193.

“In the past, several methods to approximate the continuum
have been applied: straight line, second-order polynomial,
cubic spline, and convex hull [21], [41]–[43]Currently, most
continuum removal methods in image processing seem
to be using the convex hull method (see, for instance, the
ENVI software package). Once the continuum is calculated,
it can be removed by dividing it into the original spectrum,
assuming that the continuum has a multiplicative effect on
the spectrum. This process results in a continuum-removed
spectrum, or hull quotient, because it is the ratio between
the spectrum and the hull [40]. In this way, the absorption
features are always between zero and one and can read-
ily be subjected to further analysis. Alternatively, the hull
could also be subtracted from the spectra, which would be
relevant when using apparent absorbance spectra [21].”

and its Fig. 2:

‌

[Andreas Janz]

Interesting, so we could have three outputs: i) the hull, ii) the hull quotient and iii) the hull difference.

Just to be sure: is this the correct formular for the hull difference?
out_contrem = input - out_convhull

Could you check the IDL code that produced the screenshot above?

[Agustin Lobo]

so we could have three outputs: i) the hull, ii) the hull quotient and iii) the hull difference

Yes, but please make selecting the outputs optional (ratio and/or difference) so that the user can save disk space and writing time.

Could you check the IDL code that produced the screenshot above?

Not easily, but will try.

‌

‌

[Agustin Lobo]

Note that HypPy also has the option and there you have the python code:

‌

[Andreas Janz]

Yes, but please make selecting the outputs optional (ratio and/or difference) so that the user can save disk space and writing time.

Yes, all 3 outputs can be skipped by the user:

‌

[Andreas Janz]

Note that HypPy also has the option and there you have the python code:

I’ll check it out

‌

[Agustin Lobo]

Regarding IDL DISPEC, the authors provide the binary *.sav file only:

https://www.itc.nl/about-itc/organization/scientific-departments/applied-earth-sciences/geological-remote-sensing/

Perhaps you could contact them to offer porting the IDL source to python.

[Andreas Janz]

I’ll do that, in case I didn’t find it in HypPy source.

[Andreas Janz]

BTW, I think we should have this option …

‌

… available here:

‌

[Andreas Janz]

ConvexHull and ContinuumRemoved-HullQuotient done:

‌

[Andreas Janz]

@Agustin Lobo in the HypPy source I found the formular for Hull_Difference as:

‌

[Andreas Janz]

resolves #926

→ <<cset 9f39d83dc649>>