processing

organization

dhdt.processing.matching_tools_organization.estimate_precision(C, di, dj, method='gaussian')

given a similarity surface, estimate its matching dispersion. See for a full describtion [Al21]

Parameters

• C (numpy.ndarray, size=(m,n)) – cross correlation function

• di (float) – locational estimate of the cross-correlation peak, this can also be negative, see Notes below

• dj (float) – locational estimate of the cross-correlation peak, this can also be negative, see Notes below

• method ({'gaussian', 'hessian'}) –

Returns

• cov_ii, cov_jj (float) – co-variance estimate

• cov_ij (float) – off-diagonal co-variance estimate

Notes

It is important to know what type of coordinate systems exist, hence:

coordinate |              +--------> collumns
system 'ij'|              |
           |              |
           |       j      |  image frame
   --------+-------->     |
           |              |
           |              v
           | i            rows
           v

See also

estimate_subpixel, match_translation_of_two_subsets, hessian_spread, gauss_spread

References

Al21

Altena et al. “Correlation dispersion as a measure to better estimate uncertainty of remotely sensed glacier displacements” The cryosphere, vol.16(6) pp.2285-2300, 2021.

dhdt.processing.matching_tools_organization.estimate_subpixel(QC, subpix, m0=array([[0., 0.]]), **kwargs)

estimate the sub-pixel translational displacement, present in a cross-correlation function or spectra.

Parameters

• QC (numpy.ndarray, size=(m,n), dtype={complex,float}, ndim=2) – cross-correlation spectra or function

• subpix (string) – methodology to refine the to sub-pixel localization

• m0 (numpy.ndarray, size=(1,2), dtype=integer) – row, collumn location of the peak of the cross-correlation function

Returns

ddi,ddj – sub-pixel location of the peak/phase-plane

Return type

float

See also

match_translation_of_two_subsets, list_phase_estimators, list_peak_estimators

dhdt.processing.matching_tools_organization.estimate_translation_of_two_subsets(I1, I2, M1, M2, correlator='lucas_kan', **kwargs)

Parameters

• I1 (numpy.ndarray, size=(m,n), dtype={float,integer}) – grids with intensities

• I2 (numpy.ndarray, size=(m,n), dtype={float,integer}) – grids with intensities

• M1 (numpy.ndarray, size=(m,n), dtype=bool) – labelled array corresponding to I1,I2. Hence True denotes exclusion.

• M2 (numpy.ndarray, size=(m,n), dtype=bool) – labelled array corresponding to I1,I2. Hence True denotes exclusion.

• correlator (string, default='lucas_kanade') – methodology to use to correlate I1 and I2

Returns

• di, dj (float, unit=pixels) – relative displacement (translation) between I1 and I2

• score (float) – (dis)-similarity score

dhdt.processing.matching_tools_organization.list_differential_correlators()

” list the abbreviations of the different implemented differential correlators.

The following are implemented:

• lucas_kan

• lucas_aff

• hough_opt_flw

Returns

correlator_list

Return type

list of strings

See also

list_phase_estimators, list_peak_estimators

dhdt.processing.matching_tools_organization.list_frequency_correlators(unique=False)

list the abbreviations of the different implemented correlators

The following correlators are implemented:

• cosi_corr : cosicorr

• phas_only : phase only correlation

• symm_phas : symmetric phase correlation

• ampl_comp : amplitude compensation phase correlation

• orie_corr : orientation correlation

• grad_corr : gradient correlation

• grad_norm : normalize gradient correlation

• bina_phas : binary phase correlation

• wind_corr : windrose correlation

• gaus_phas : gaussian transformed phase correlation

• cros_corr : cross correlation

• robu_corr : robust phase correlation

• proj_phas : projected phase correlation

• phas_corr : phase correlation

Returns

correlator_list

Return type

list of strings

dhdt.processing.matching_tools_organization.list_peak_estimators()

list the abbreviations of the different implemented for the peak fitting of the similarity function.

The following methods are implemented:

• ‘gauss_1’ : 1D Gaussian fitting

• ‘parab_1’ : 1D parabolic fitting

• ‘moment’ : 2D moment of the peak

• ‘mass’ : 2D center of mass fitting

• ‘centroid’ : 2D estimate the centroid

• ‘blais’ : 1D estimate through forth order filter

• ‘ren’ : 1D parabolic fitting

• ‘birch’ : 1D weighted sum

• ‘eqang’ : 1D equiangular line fitting

• ‘trian’ : 1D triangular fitting

• ‘esinc’ : 1D exponential esinc function

• ‘gauss_2’ : 2D Gaussian fitting

• ‘parab_2’ : 2D parabolic fitting

• ‘optical_flow’ : optical flow refinement

Returns

subpix_list

Return type

list of strings

See also

list_phase_estimators, list_differential_correlators

dhdt.processing.matching_tools_organization.list_phase_estimators()

list the abbreviations of the different implemented phase plane estimation procedures

The following methods are implemented:

• ‘tpss’ : two point step size

• ‘svd’ : single value decomposition

• ‘radon’ : radon transform

• ‘hough’ : hough transform

• ‘ransac’ : random sampling and consensus

• ‘wpca’ : weighted principle component analysis

• ‘pca’ : principle component analysis

• ‘lsq’ : least squares estimation

• ‘diff’ : phase difference

Returns

subpix_list

Return type

list of strings

See also

list_peak_estimators, list_differential_correlators

dhdt.processing.matching_tools_organization.list_spatial_correlators(unique=False)

list the abbreviations of the different implemented correlators.

The following methods are implemented:

• norm_corr : normalized cross correlation

• sq_diff : sum of squared differences

• sad_diff : sum of absolute differences

• max_like : maximum likelihood

• wght_corr : weighted normalized cross correlation

Returns

correlator_list

Return type

list of strings

dhdt.processing.matching_tools_organization.match_translation_of_two_subsets(I1_sub, I2_sub, correlator, subpix, M1_sub=array([], dtype=float64), M2_sub=array([], dtype=float64))

Parameters

• I1_sub (numpy.ndarray, size={(m,n),(k,l)}, dtype={float,int}, ndim={2,3}) – grid with intensities, a.k.a. template to locate

• I2_sub (numpy.ndarray, size=(m,n), dtype={float,int}, ndim={2,3}) – grid with intensities, a.k.a. search space

• correlator (string) – methodology to use to correlate I1 and I2

• subpix (string) – methodology to refine the to sub-pixel localization

• M1_sub (numpy.ndarray, size=(m,n), dtype=bool) – labels corresponding to I1_sub,I2_sub. Hence True denotes exclusion.

• M2_sub (numpy.ndarray, size=(m,n), dtype=bool) – labels corresponding to I1_sub,I2_sub. Hence True denotes exclusion.

Returns

{Q,C} – cross-correlation spectra or function

Return type

numpy.ndarray, size={(m,n),(k,l)}, dtype={complex,float}, ndim=2

See also

estimate_subpixel, estimate_precision, list_frequency_correlators, list_spatial_correlators, list_phase_estimators, list_peak_estimators

frequency filters

dhdt.processing.matching_tools_frequency_filters.adaptive_masking(S, m=0.9)

mark significant intensities in spectrum, following [Le07].

Parameters

• S (numpy.array, size=(m,n), dtype=complex) – array with spectrum, i.e.: S = np.fft.fft2(I)

• m (float, default=.9) – cut-off intensity in respect to maximum

Returns

M – frequency mask

Return type

numpy.array, size=(m,n), dtype=bool

See also

tpss

References

Le07

Leprince, et.al. “Automatic and precise orthorectification, coregistration, and subpixel correlation of satellite images, application to ground deformation measurements”, IEEE Transactions on geoscience and remote sensing vol. 45.6 pp. 1529-1558, 2007.

dhdt.processing.matching_tools_frequency_filters.blackman_window(Z)

create two-dimensional Blackman filter

Parameters

Z (numpy.array, size=(m,n)) – array with intensities

Returns

W – weighting mask

Return type

numpy.array, size=(m,n), dtype=bool

See also

raised_cosine, cosine_bell, high_pass_circle, hamming_window, hanning_window

dhdt.processing.matching_tools_frequency_filters.construct_phase_plane(Z, di, dj, indexing='ij')

given a displacement, create what its phase plane in Fourier space

Parameters

• Z (np.array, size=(m,n)) – image domain

• di (float) – displacement along the vertical axis

• dj (float) – displacement along the horizantal axis

• indexing ({‘xy’, ‘ij’}) –

  • “xy” : using map coordinates

  • ”ij” : using local image coordinates

Returns

Q – array with phase angles

Return type

np.array, size=(m,n), complex

Notes

Two different coordinate system are used here:

indexing   |           indexing    ^ y
system 'ij'|           system 'xy' |
           |                       |
           |       i               |       x
   --------+-------->      --------+-------->
           |                       |
           |                       |
image      | j         map         |
based      v           based       |

dhdt.processing.matching_tools_frequency_filters.construct_phase_values(IJ, di, dj, indexing='ij', system='radians')

given a displacement, create what its phase plane in Fourier space

Parameters

• IJ (np.array, size=(_,2), dtype=float) – locations of phase values

• di (float) – displacment along the vertical axis

• dj (float) – displacment along the horizantal axis

• indexing ({‘radians’ (default), ‘unit’, 'normalized'}) –

  • “xy” : using map coordinates

  • ”ij” : using local image coordinates

• indexing –

  the extent of the cross-spectrum can span different ranges

  • ”radians” : -pi..+pi

  • ”unit” : -1…+1

  • ”normalized” : -0.5…+0.5

Returns

Q – array with phase angles

Return type

np.array, size=(_,1), complex

dhdt.processing.matching_tools_frequency_filters.cosine_bell(Z)

cosine bell filter

Parameters

Z (numpy.array, size=(m,n)) – array with intensities

Returns

W – weighting mask

Return type

numpy.array, size=(m,n), dtype=float

See also

raised_cosine

dhdt.processing.matching_tools_frequency_filters.cross_shading_filter(Q, az_1, az_2)

Parameters

• Q (numpy.array, size=(m,n)) – cross-power spectrum

• az_1 (float, unit=degrees, range=-180...+180) – illumination direction of template

• az_2 (float, unit=degrees, range=-180...+180) – illumination direction of search space

Returns

W – weigthing grid

Return type

numpy.array, size=(m,n), dtype=float

References

Wa15

Wan, “Phase correlation-based illumination-insensitive image matching for terrain-related applications” PhD dissertation at Imperical College London, 2015.

dhdt.processing.matching_tools_frequency_filters.cross_spectrum_to_coordinate_list(data, W=array([], dtype=float64))

if data is given in array for, then transform it to a coordinate list

Parameters

• data (np.array, size=(m,n), dtype=complex) – cross-spectrum

• W (np.array, size=(m,n), dtype=boolean) – weigthing matrix of the cross-spectrum

Returns

data_list – coordinate list with angles, in normalized ranges, i.e: -1 … +1

Return type

np.array, size=(m*n,3), dtype=float

dhdt.processing.matching_tools_frequency_filters.gaussian_mask(S)

mask significant intensities in spectrum, following [Ec08].

Parameters

S (numpy.array, size=(m,n), dtype=complex) – array with spectrum, i.e.: S = np.fft.fft2(I)

Returns

M – frequency mask

Return type

numpy.array, size=(m,n), dtype=bool

See also

tpss

References

Ec08

Eckstein et al. “Phase correlation processing for DPIV measurements”, Experiments in fluids, vol.45 pp.485-500, 2008.

Examples

>>> import numpy as np
>>> from dhdt.testing.matching_tools import create_sample_image_pair

>>> im1,im2,_,_,_ = create_sample_image_pair(d=2**4, max_range=1)
>>> spec1,spec2 = np.fft.fft2(im1), np.fft.fft2(im2)
>>> Q = spec1 * np.conjugate(spec2) # Fourier based image matching
>>> Qn = normalize_spectrum(Q)

>>> W = gaussian_mask(Q)
>>> C = np.fft.ifft2(W*Q)

dhdt.processing.matching_tools_frequency_filters.gradient_fourier(Z)

spectral derivative estimation

Parameters

Z (numpy.ndarray, size=(m,n), dtype=float) – intensity array

Returns

dZdx_1, dZdx_2 – first order derivative along both axis

Return type

numpy.ndarray, size=(m,n), dtype=float

dhdt.processing.matching_tools_frequency_filters.hamming_window(Z)

create two-dimensional Hamming filter

Parameters

Z (numpy.array, size=(m,n)) – array with intensities

Returns

W – weighting mask

Return type

numpy.array, size=(m,n), dtype=bool

See also

raised_cosine, cosine_bell, high_pass_circle, blackman_window, hamming_window

dhdt.processing.matching_tools_frequency_filters.hanning_window(Z)

create two-dimensional Hanning filter, also known as Cosine Bell

Parameters

Z (numpy.array, size=(m,n)) – array with intensities

Returns

W – weighting mask

Return type

numpy.array, size=(m,n), dtype=bool

See also

raised_cosine, cosine_bell, high_pass_circle, blackman_window, hamming_window

dhdt.processing.matching_tools_frequency_filters.high_pass_circle(Z, r=0.5)

create hard high-pass filter

Parameters

• Z (numpy.array, size=(m,n)) – array with intensities

• r (float, default=0.5) – radius of the circle, r=.5 is same as its width

Returns

W – weighting mask

Return type

numpy.array, size=(m,n), dtype=bool

See also

raised_cosine, cosine_bell, low_pass_circle

dhdt.processing.matching_tools_frequency_filters.kaiser_window(Z, beta=14.0)

create two dimensional Kaiser filter

Parameters

• Z (numpy.array, size=(m,n)) – array with intensities

• beta (float) – 0.0 - rectangular window 5.0 - similar to Hamming window 6.0 - similar to Hanning window 8.6 - similar to Blackman window

Returns

W – weighting mask

Return type

numpy.array, size=(m,n), dtype=bool

See also

raised_cosine, cosine_bell, high_pass_circle, hamming_window, hanning_window

dhdt.processing.matching_tools_frequency_filters.local_coherence(Q, ds=1)

estimate the local coherence of a spectrum

Parameters

• Q (numpy.ndarray, size=(m,n), dtype=complex) – array with cross-spectrum, with centered coordinate frame

• ds (integer, default=1) – kernel radius to describe the neighborhood

Returns

M – vector coherence from no to ideal, i.e.: 0…1

Return type

numpy.ndarray, size=(m,n), dtype=float

See also

thresh_masking

Examples

>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> from dhdt.testing.matching_tools import create_sample_image_pair

>>> # create cross-spectrum with random displacement
>>> im1,im2,_,_,_ = create_sample_image_pair(d=2**4, max_range=1)
>>> spec1,spec2 = np.fft.fft2(im1), np.fft.fft2(im2)
>>> Q = spec1 * np.conjugate(spec2)
>>> Q = normalize_spectrum(Q)
>>> Q = np.fft.fftshift(Q) # transform to centered grid

>>> C = local_coherence(Q)

>>> plt.figure(), plt.imshow(C, cmap='OrRd'), plt.colorbar(), plt.show()
>>> plt.figure(), plt.imshow(np.angle(Q), cmap='twilight'),
>>> plt.colorbar(), plt.show()

dhdt.processing.matching_tools_frequency_filters.low_pass_bell(Z, r=0.5)

create low-pass two-dimensional filter with a bell shape, see also [Ta03].

Parameters

• Z (numpy.ndarray, size=(m,n)) – array with intensities

• r (float, default=0.5) – radius of the mother rectangle, r=.5 is same as its width

Returns

W – weighting mask

Return type

numpy.ndarray, size=(m,n), dtype=bool

See also

low_pass_rectancle, low_pass_circle, low_pass_pyramid

References

Ta03

Takita et al. “High-accuracy subpixel image registration based on phase-only correlation” IEICE transactions on fundamentals of electronics, communications and computer sciences, vol.86(8) pp.1925-1934, 2003.

dhdt.processing.matching_tools_frequency_filters.low_pass_circle(Z, r=0.5)

create hard two-dimensional low-pass filter

Parameters

• Z (numpy.array, size=(m,n)) – array with intensities

• r (float, default=0.5) – radius of the circle, r=.5 is same as its width

Returns

W – weighting mask

Return type

numpy.array, size=(m,n), dtype=bool

See also

raised_cosine, cosine_bell, high_pass_circle

dhdt.processing.matching_tools_frequency_filters.low_pass_ellipse(Z, r1=0.5, r2=0.5)

create hard low-pass filter

Parameters

• Z (numpy.array, size=(m,n)) – array with intensities

• r1 (float, default=0.5) – radius of the ellipse along the first axis, r=.5 is same as its width

• r2 (float, default=0.5) – radius of the ellipse along the second axis

Returns

W – weighting mask

Return type

numpy.array, size=(m,n), dtype=bool

See also

raised_cosine, cosine_bell, high_pass_circle, low_pass_circle

dhdt.processing.matching_tools_frequency_filters.low_pass_pyramid(Z, r=0.5)

create low-pass two-dimensional filter with pyramid shape, see also [Ta03].

Parameters

• Z (numpy.array, size=(m,n)) – array with intensities

• r (float, default=0.5) – radius of the mother rectangle, r=.5 is same as its width

Returns

W – weighting mask

Return type

numpy.array, size=(m,n), dtype=bool

See also

low_pass_rectancle, low_pass_circle, low_pass_bell

References

Ta03

Takita et al. “High-accuracy subpixel image registration based on phase-only correlation” IEICE transactions on fundamentals of electronics, communications and computer sciences, vol.86(8) pp.1925-1934, 2003.

dhdt.processing.matching_tools_frequency_filters.low_pass_rectancle(Z, r=0.5)

create hard two dimensional low-pass filter

Parameters

• Z (numpy.array, size=(m,n)) – array with intensities

• r (float, default=0.5) – radius of the rectangle, r=.5 is same as its width

Returns

W – weighting mask

Return type

numpy.array, size=(m,n), dtype=bool

See also

low_pass_circle, low_pass_pyramid, low_pass_bell

References

Ta03

Takita et al. “High-accuracy subpixel image registration based on phase-only correlation” IEICE transactions on fundamentals of electronics, communications and computer sciences, vol.86(8) pp.1925-1934, 2003.

dhdt.processing.matching_tools_frequency_filters.make_fourier_grid(Q, indexing='ij', system='radians', shift=True, axis='center')

The four quadrants of the coordinate system of the discrete Fourier transform are flipped. This function gives its coordinate system as it would be in a map (xy) or pixel based (ij) system.

Parameters

• Q (numpy.ndarray, size=(m,n), dtype=complex) – Fourier based (cross-)spectrum.

• indexing ({‘xy’, ‘ij’}) –

  • “xy” : using map coordinates

  • ”ij” : using local image coordinates

• system ({‘radians’, ‘unit’, 'normalized'}) –

  the extent of the cross-spectrum can span from

  • ”radians” : -pi..+pi (default)

  • ”unit” : -1…+1

  • ”normalized” : -0.5…+0.5

  • ”pixel” : -m/2…+m/2

• shift (boolean, default=True) – the coordinate frame of a Fourier transform is flipped

Returns

• F_1 (numpy.ndarray, size=(m,n), dtype=integer) – first coordinate index of the Fourier spectrum in a map system.

• F_2 (numpy.ndarray, size=(m,n), dtype=integer) – second coordinate index of the Fourier spectrum in a map system.

Notes

metric system:         Fourier-based flip
       y               ┼------><------┼
       ^               |              |
       |               |              |
       |               v              v
<------┼-------> x
       |               ^              ^
       |               |              |
       v               ┼------><------┼

It is important to know what type of coordinate systems exist, hence:

coordinate |           coordinate  ^ y
system 'ij'|           system 'xy' |
           |                       |
           |       j               |       x
   --------┼-------->      --------┼-------->
           |                       |
           |                       |
           | i                     |
           v                       |

dhdt.processing.matching_tools_frequency_filters.make_template_float(Z)

templates for frequency matching should be float and ideally have limited border issues.

Parameters

Z (numpy.ndarray, size={(m,n), (m,n,b)}) – grid with intensity values

Returns

Z – grid with intensity values

Return type

numpy.ndarray, size={(m,n), (m,n,b)}, dtype=float, range=0…1

dhdt.processing.matching_tools_frequency_filters.normalize_power_spectrum(Q)

transform spectrum to complex vectors with unit length

Parameters

Q (numpy.ndarray, size=(m,n), dtype=complex) – cross-spectrum

Returns

Qn – normalized cross-spectrum, that is elements with unit length

Return type

numpy.ndarray, size=(m,n), dtype=complex

Examples

>>> import numpy as np
>>> from dhdt.testing.matching_tools import create_sample_image_pair

>>> im1,im2,_,_,_ = create_sample_image_pair(d=2**4, max_range=1)
>>> spec1,spec2 = np.fft.fft2(im1), np.fft.fft2(im2)
>>> Q = spec1 * np.conjugate(spec2) # fourier based image matching
>>> Qn = normalize_spectrum(Q)

dhdt.processing.matching_tools_frequency_filters.perdecomp(img)

calculate the periodic and smooth components of an image, based upon [Mo11].

Parameters

img (numpy.ndarray, size=(m,n)) – array with intensities

Returns

• per (numpy.ndarray, size=(m,n)) – periodic component

• cor (numpy.ndarray, size=(m,n)) – smooth component

References

Mo11

Moisan, L. “Periodic plus smooth image decomposition”, Journal of mathematical imaging and vision vol. 39.2 pp. 161-179, 2011.

Examples

>>> import numpy as np
>>> from dhdt.testing.matching_tools import create_sample_image_pair

>>> im1,_,_,_,_ = create_sample_image_pair(d=2**7, max_range=1)
>>> per,cor = perdecomp(im1)

>>> spec1 = np.fft.fft2(per)

dhdt.processing.matching_tools_frequency_filters.raised_cosine(Z, beta=0.35)

raised cosine filter, based on [St01] and used by [Le07].

Parameters

• Z (numpy.array, size=(m,n)) – array with intensities

• beta (float, default=0.35) – roll-off factor

Returns

W – weighting mask

Return type

numpy.array, size=(m,n), dtype=float

See also

tpss

References

St01

Stone et al. “A fast direct Fourier-based algorithm for subpixel registration of images.” IEEE Transactions on geoscience and remote sensing. vol. 39(10) pp. 2235-2243, 2001.

Le07

Leprince, et.al. “Automatic and precise orthorectification, coregistration, and subpixel correlation of satellite images, application to ground deformation measurements”, IEEE Transactions on geoscience and remote sensing vol. 45.6 pp. 1529-1558, 2007.

Examples

>>> import numpy as np
>>> from dhdt.testing.matching_tools import create_sample_image_pair

>>> im1,im2,_,_,_ = create_sample_image_pair(d=2**4, max_range=1)
>>> spec1,spec2 = np.fft.fft2(im1), np.fft.fft2(im2)

>>> rc1 = raised_cosine(spec1, beta=0.35)
>>> rc2 = raised_cosine(spec2, beta=0.50)

>>> Q = (rc1*spec1) * np.conjugate((rc2*spec2)) # Fourier based image matching
>>> Qn = normalize_spectrum(Q)

dhdt.processing.matching_tools_frequency_filters.thresh_masking(S, m=0.0001, s=10)

Mask significant intensities in spectrum, following [St01] and [Le07].

Parameters

• S (numpy.array, size=(m,n), dtype=complex) – array with spectrum, i.e.: S = np.fft.fft2(I)

• m (float, default=1e-3) – cut-off intensity in respect to maximum

• s (integer, default=10) – kernel size of the median filter

Returns

M – frequency mask

Return type

numpy.array, size=(m,n), dtype=bool

See also

tpss

References

St01

Stone et al. “A fast direct Fourier-based algorithm for subpixel registration of images.” IEEE Transactions on geoscience and remote sensing vol. 39(10) pp. 2235-2243, 2001.

Le07

Leprince, et.al. “Automatic and precise orthorectification, coregistration, and subpixel correlation of satellite images, application to ground deformation measurements”, IEEE Transactions on geoscience and remote sensing vol. 45.6 pp. 1529-1558, 2007.

frequency correlators

dhdt.processing.matching_tools_frequency_correlators.amplitude_comp_corr(I1, I2, F_0=0.04)

match two imagery through amplitude compensated phase correlation

Parameters

• I1 (numpy.array, size=(m,n), ndim={2,3}) – array with intensities, if multiple bands are given the cross- correlation spectra are stacked, see [AL22].

• I2 (numpy.array, size=(m,n), ndim={2,3}) – array with intensities

• F_0 (float, default=4e-2) – cut-off intensity in respect to maximum

Returns

Q – cross-power spectrum

Return type

numpy.array, size=(m,n), dtype=complex

References

Mu88

Mu et al. “Amplitude-compensated matched filtering”, Applied optics, vol. 27(16) pp. 3461-3463, 1988.

AL22

Altena & Leinss, “Improved surface displacement estimation through stacking cross-correlation spectra from multi-channel imagery”, Science of Remote Sensing, vol.6 pp.100070, 2022.

Examples

>>> import numpy as np
>>> from dhdt.processing.matching_tools import get_integer_peak_location
>>> from dhdt.testing.matching_tools import create_sample_image_pair

>>> im1,im2,ti,tj,_ = create_sample_image_pair(d=2**4, max_range=1)
>>> Q = amplitude_comp_corr(im1, im2)
>>> C = np.fft.ifft2(Q)
>>> di,dj,_,_ = get_integer_peak_location(C)

>>> assert(np.isclose(ti, di, atol=1))
>>> assert(np.isclose(ti, di, atol=1))

dhdt.processing.matching_tools_frequency_correlators.binary_orientation_corr(I1, I2)

match two imagery through binary phase only correlation, see [KJ91].

Parameters

• I1 (numpy.array, size=(m,n), ndim={2,3}) – array with intensities, if multiple bands are given the cross- correlation spectra are stacked, see [AL22].

• I2 (numpy.array, size=(m,n), ndim={2,3}) – array with intensities

Returns

Q – cross-power spectrum

Return type

numpy.array, size=(m,n), dtype=complex

See also

dhdt.processing.matching_tools_frequency_correlators.orientation_corr, dhdt.processing.matching_tools_frequency_correlators.phase_only_corr

Notes

The matching equations are as follows:

\[\mathbf{S}_1, \mathbf{S}_2 = \mathcal{F}[\mathbf{I}_1],\quad \mathcal{F}[\mathbf{I}_2]\]

\[\mathbf{W} = sign(\mathfrak{Re} [ \mathbf{S}_2 ] )\]

\[\mathbf{Q}_{12} = \mathbf{S}_1 [\mathbf{W}\mathbf{S}_2]^{\star}\]

where \(\mathcal{F}\) denotes the Fourier transform and \(\star\) a complex conjugate operation

References

KJ91

Kumar & Juday, “Design of phase-only, binary phase-only, and complex ternary matched filters with increased signal-to-noise ratios for colored noise”, Optics letters, vol.16(13) pp.1025-1027, 1991.

AL22

Altena & Leinss, “Improved surface displacement estimation through stacking cross-correlation spectra from multi-channel imagery”, Science of Remote Sensing, vol.6 pp.100070, 2022.

Examples

>>> import numpy as np
>>> from dhdt.processing.matching_tools import get_integer_peak_location
>>> from dhdt.testing.matching_tools import create_sample_image_pair

>>> im1,im2,ti,tj,_ = create_sample_image_pair(d=2**4, max_range=1)
>>> Q = binary_orientation_corr(im1, im2)
>>> C = np.fft.ifft2(Q)
>>> di,dj,_,_ = get_integer_peak_location(C)

>>> assert(np.isclose(ti, di, atol=1))
>>> assert(np.isclose(ti, di, atol=1))

dhdt.processing.matching_tools_frequency_correlators.cosine_corr(I1, I2)

match two imagery through discrete cosine transformation, see [Li04].

Parameters

• I1 (numpy.array, size=(m,n), ndim={2,3}, dtype=float) – array with intensities, if multiple bands are given the cross- correlation spectra are stacked, see [AL22].

• I2 (numpy.array, size=(m,n), ndim={2,3}, dtype=float) – array with intensities

Returns

C – displacement surface

Return type

numpy.array, size=(m,n), dtype=real

See also

dhdt.processing.matching_tools_frequency_correlators.create_complex_DCT, dhdt.processing.matching_tools_frequency_correlators.sign_only_corr

References

Li04

Li, et al. “DCT-based phase correlation motion estimation”, IEEE international conference on image processing, vol. 1, 2004.

AL22

Altena & Leinss, “Improved surface displacement estimation through stacking cross-correlation spectra from multi-channel imagery”, Science of Remote Sensing, vol.6 pp.100070, 2022.

dhdt.processing.matching_tools_frequency_correlators.cross_corr(I1, I2)

match two imagery through cross correlation in FFT

Parameters

• I1 (numpy.array, size=(m,n), ndim={2,3}) – array with intensities, if multiple bands are given the cross- correlation spectra are stacked, see [AL22].

• I2 (numpy.array, size=(m,n), ndim={2,3}) – array with intensities

Returns

Q – cross-power spectrum

Return type

numpy.array, size=(m,n), dtype=complex

See also

dhdt.processing.matching_tools_frequency_correlators.phase_corr, dhdt.processing.matching_tools_frequency_correlators.upsampled_cross_corr

Notes

The matching equations are as follows:

\[\mathbf{S}_1, \mathbf{S}_2 = \mathcal{F}[\mathbf{I}_1],\quad \mathcal{F}[\mathbf{I}_2]\]

\[\mathbf{Q}_{12} = \mathbf{S}_1 {\mathbf{S}_2}^{\star}\]

where \(\mathcal{F}\) denotes the Fourier transform and \(\star\) a complex conjugate operation.

Schematically this looks like:

I1  ┌----┐ S1
 ---┤ F  ├-----┐
    └----┘     | Q12 ┌------┐  C12
               ×-----┤ F^-1 ├---
I2  ┌----┐ S2  |     └------┘
 ---┤ F  ├-----┘
    └----┘

References

HK12

Heid & Kääb. “Evaluation of existing image matching methods for deriving glacier surface displacements globally from optical satellite imagery”, Remote sensing of environment, vol.118 pp.339-355, 2012.

AL22

Altena & Leinss, “Improved surface displacement estimation through stacking cross-correlation spectra from multi-channel imagery”, Science of Remote Sensing, vol.6 pp.100070, 2022.

Examples

>>> import numpy as np
>>> from dhdt.processing.matching_tools import get_integer_peak_location
>>> from dhdt.testing.matching_tools import create_sample_image_pair

>>> im1,im2,ti,tj,_ = create_sample_image_pair(d=2**4, max_range=1)
>>> Q = cross_corr(im1, im2)
>>> C = np.fft.ifft2(Q)
>>> di,dj,_,_ = get_integer_peak_location(C)

>>> assert(np.isclose(ti, di, atol=1))
>>> assert(np.isclose(ti, di, atol=1))

dhdt.processing.matching_tools_frequency_correlators.gaussian_transformed_phase_corr(I1, I2)

match two imagery through Gaussian transformed phase correlation, see [Ec08].

Parameters

• I1 (numpy.array, size=(m,n), ndim={2,3}) – array with intensities, if multiple bands are given the cross- correlation spectra are stacked, see [AL22].

• I2 (numpy.array, size=(m,n), ndim={2,3}) – array with intensities

Returns

Q – cross-power spectrum

Return type

numpy.array, size=(m,n), dtype=complex

See also

dhdt.processing.matching_tools_frequency_correlators.phase_corr

Notes

The matching equations are as follows:

\[\mathbf{S}_1, \mathbf{S}_2 = \mathcal{F}[\mathbf{I}_1],\quad \mathcal{F}[\mathbf{I}_2]\]

\[\mathbf{Q}_{12} = \mathbf{S}_1 {\mathbf{S}_2}^{\star}\]

\[\mathbf{Q}_{12} / \| \mathbf{Q}_{12}\|\]

\[\mathbf{Q}_{12} = \mathbf{W} \mathbf{Q}_{12}\]

where \(\mathcal{F}\) denotes the Fourier transform and \(\star\) a complex conjugate operation.

Schematically this looks like:

I1  ┌----┐ S1                  ┌----┐
 ---┤ F  ├-----┐               | W  |
    └----┘     | Q12 ┌-------┐ └-┬--┘ ┌------┐ C12
               ×-----┤ 1/|x| ├---┴----┤ F^-1 ├----
I2  ┌----┐ S2  |     └-------┘        └------┘
 ---┤ F  ├-----┘
    └----┘

References

Ec08

Eckstein et al. “Phase correlation processing for DPIV measurements”, Experiments in fluids, vol.45 pp.485-500, 2008.

AL22

Altena & Leinss, “Improved surface displacement estimation through stacking cross-correlation spectra from multi-channel imagery”, Science of Remote Sensing, vol.6 pp.100070, 2022.

Examples

>>> import numpy as np
>>> from dhdt.processing.matching_tools import get_integer_peak_location
>>> from dhdt.testing.matching_tools import create_sample_image_pair

>>> im1,im2,ti,tj,_ = create_sample_image_pair(d=2**4, max_range=1)
>>> Q = gaussian_transformed_phase_corr(im1, im2)
>>> C = np.fft.ifft2(Q)
>>> di,dj,_,_ = get_integer_peak_location(C)

>>> assert(np.isclose(ti, di, atol=1))
>>> assert(np.isclose(ti, di, atol=1))

dhdt.processing.matching_tools_frequency_correlators.gradient_corr(I1, I2)

match two imagery through gradient correlation

Parameters

• I1 (numpy.array, size=(m,n), ndim={2,3}) – array with intensities, if multiple bands are given the cross- correlation spectra are stacked, see [AL22].

• I2 (numpy.array, size=(m,n), ndim={2,3}) – array with intensities

Returns

Q – cross-power spectrum

Return type

numpy.array, size=(m,n), dtype=complex

See also

dhdt.processing.matching_tools_frequency_correlators.normalized_gradient_corr, dhdt.processing.matching_tools_frequency_correlators.phase_corr, dhdt.processing.matching_tools_frequency_correlators.orientation_corr

Notes

The matching equations are as follows:

\[\mathbf{G}_1, \mathbf{G}_2 = \partial_{x}\mathbf{I}_1 + i \partial_{y}\mathbf{I}_1, \quad \partial_{x}\mathbf{I}_2 + i \partial_{y}\mathbf{I}_2\]

\[\mathbf{S}_1, \mathbf{S}_2 = \mathcal{F}[\mathbf{G}_1], \quad \mathcal{F}[\mathbf{G}_2]\]

\[\mathbf{Q}_{12} = \mathbf{S}_1 \mathbf{S}_2^{\star}\]

where \(\partial_{x}\) is the spatial derivate in horizontal direction, \(\mathcal{F}\) denotes the Fourier transform and \(\star\) a complex conjugate operation.

Schematically this looks like:

          ∂x I1
I1  ┌----┬-----┬----┐ S1
 ---┤ ∂  |     | F  ├------┐
    └----┴-----┴----┘      |
          ∂y I1            |
                           | Q12 ┌------┐  C12
          ∂x I2            ×-----┤ F^-1 ├---
I2  ┌----┬-----┬----┐ S2   |     └------┘
 ---┤ ∂  |     | F  ├------┘
    └----┴-----┴----┘
          ∂y I2

References

AV03

Argyriou & Vlachos. “Estimation of sub-pixel motion using gradient cross-correlation”, Electronics letters, vol.39(13) pp.980–982, 2003.

Tz10

Tzimiropoulos et al. “Subpixel registration with gradient correlation” IEEE transactions on image processing, vol.20(6) pp.1761–1767, 2010.

AL22

Altena & Leinss, “Improved surface displacement estimation through stacking cross-correlation spectra from multi-channel imagery”, Science of Remote Sensing, vol.6 pp.100070, 2022.

Examples

>>> import numpy as np
>>> from dhdt.processing.matching_tools import get_integer_peak_location
>>> from dhdt.testing.matching_tools import create_sample_image_pair

>>> im1,im2,ti,tj,_ = create_sample_image_pair(d=2**4, max_range=1)
>>> Q = gradient_corr(im1, im2)
>>> C = np.fft.ifft2(Q)
>>> di,dj,_,_ = get_integer_peak_location(C)

>>> assert(np.isclose(ti, di, atol=1))
>>> assert(np.isclose(ti, di, atol=1))

dhdt.processing.matching_tools_frequency_correlators.masked_corr(I1, I2, M1=array([], dtype=float64), M2=array([], dtype=float64))

match two imagery through masked normalized cross-correlation in FFT, see [Pa11].

Parameters

• I1 ({numpy.array, numpy.ma}, size=(m,n), ndim=2) – array with intensities

• I2 ({numpy.array, numpy.ma}, size=(m,n), ndim=2) – array with intensities

• M1 (numpy.array, size=(m,n)) – array with mask, hence False’s are neglected

• M2 (numpy.array, size=(m,n)) – array with mask, hence False’s are neglected

Returns

NCC – correlation surface

Return type

numpy.array, size=(m,n)

See also

dhdt.processing.matching_tools_spatial_correlators.normalized_cross_corr, dhdt.processing.matching_tools_spatial_correlators.weighted_normalized_cross_corr

References

Pa11

Padfield. “Masked object registration in the Fourier domain”, IEEE transactions on image processing, vol. 21(5) pp. 2706-2718, 2011.

Examples

>>> import numpy as np
>>> from dhdt.processing.matching_tools import get_integer_peak_location
>>> from dhdt.testing.matching_tools import create_sample_image_pair

>>> im1,im2,ti,tj,_ = create_sample_image_pair(d=2**4, max_range=1)
>>> msk1,msk2 = np.ones_like(im1), np.ones_like(im2)
>>> Q = masked_corr(im1, im2, msk1, msk2)
>>> C = np.fft.ifft2(Q)
>>> di,dj,_,_ = get_integer_peak_location(C)

>>> assert(np.isclose(ti, di, atol=1))
>>> assert(np.isclose(ti, di, atol=1))

dhdt.processing.matching_tools_frequency_correlators.masked_cosine_corr(I1, I2, M1, M2)

Parameters

• I1 ({numpy.array, numpy.ma}, size=(m,n), ndim=2) – array with intensities

• I2 ({numpy.array, numpy.ma}, size=(m,n), ndim=2) – array with intensities

• M1 (numpy.array, size=(m,n), ndim=2, dtype={bool,float}) – array with mask, hence False’s are neglected

• M2 (numpy.array, size=(m,n), ndim=2, dtype={bool,float}) – array with mask, hence False’s are neglected

dhdt.processing.matching_tools_frequency_correlators.normalized_gradient_corr(I1, I2)

match two imagery through gradient correlation

Parameters

• I1 (numpy.array, size=(m,n), ndim={2,3}) – array with intensities, if multiple bands are given the cross- correlation spectra are stacked, see [AL22].

• I2 (numpy.array, size=(m,n), ndim={2,3}) – array with intensities

Returns

Q – cross-power spectrum

Return type

numpy.array, size=(m,n), dtype=complex

See also

dhdt.processing.matching_tools_frequency_correlators.gradient_corr, dhdt.processing.matching_tools_frequency_correlators.phase_corr, dhdt.processing.matching_tools_frequency_correlators.windrose_corr

Notes

The matching equations are as follows:

\[\mathbf{G}_1, \mathbf{G}_2 = \partial_{x}\mathbf{I}_1 + i \partial_{y}\mathbf{I}_1, \quad \partial_{x}\mathbf{I}_2 + i \partial_{y}\mathbf{I}_2\]

\[\mathbf{S}_1, \mathbf{S}_2 = \mathcal{F}[\dot{\mathbf{G}}_1], \quad \mathcal{F}[\dot{\mathbf{G}}_2]\]

\[\mathbf{Q}_{12} = \mathbf{S}_1 \mathbf{S}_2^{\star}\]

where \(\partial_{x}\) is the spatial derivate in horizontal direction, \(\mathcal{F}\) denotes the Fourier transform and \(\star\) a complex conjugate operation and \(\dot{x}\) represents a statistical normalization through the sampled mean and standard deviation.

Schematically this looks like:

          ∂x I1
I1  ┌----┬-----┬----┐ S1 ┌------┐
 ---┤ ∂  |     | F  ├----┤ norm ├--┐
    └----┴-----┴----┘    └------┘  |
          ∂y I1                    |
                                   | Q12 ┌------┐  C12
          ∂x I2                    ×-----┤ F^-1 ├---
I2  ┌----┬-----┬----┐ S2 ┌------┐  |     └------┘
 ---┤ ∂  |     | F  ├----┤ norm ├--┘
    └----┴-----┴----┘    └------┘
          ∂y I2

References

Tz10

Tzimiropoulos et al. “Subpixel registration with gradient correlation” IEEE transactions on image processing, vol.20(6) pp.1761–1767, 2010.

AL22

Altena & Leinss, “Improved surface displacement estimation through stacking cross-correlation spectra from multi-channel imagery”, Science of Remote Sensing, vol.6 pp.100070, 2022.

Examples

>>> import numpy as np
>>> from dhdt.processing.matching_tools import get_integer_peak_location
>>> from dhdt.testing.matching_tools import create_sample_image_pair

>>> im1,im2,ti,tj,_ = create_sample_image_pair(d=2**4, max_range=1)
>>> Q = normalized_gradient_corr(im1, im2)
>>> C = np.fft.ifft2(Q)
>>> di,dj,_,_ = get_integer_peak_location(C)

>>> assert(np.isclose(ti, di, atol=1))
>>> assert(np.isclose(ti, di, atol=1))

dhdt.processing.matching_tools_frequency_correlators.orientation_corr(I1, I2)

match two imagery through orientation correlation, developed by [Fi02] while demonstrated its benefits over glaciated terrain by [HK12].

Parameters

• I1 (numpy.array, size=(m,n), ndim={2,3}) – array with intensities, if multiple bands are given the cross- correlation spectra are stacked, see [AL22].

• I2 (numpy.array, size=(m,n), ndim={2,3}) – array with intensities

Returns

Q – cross-power spectrum

Return type

numpy.array, size=(m,n), dtype=complex

See also

dhdt.processing.matching_tools_frequency_correlators.phase_corr, dhdt.processing.matching_tools_frequency_correlators.windrose_corr, dhdt.processing.matching_tools_spatial_correlators.cosine_similarity

Notes

The matching equations are as follows:

\[\mathbf{G}_1, \mathbf{G}_2 = \partial_{x}\mathbf{I}_1 + i \partial_{y}\mathbf{I}_1,\quad \partial_{x}\mathbf{I}_2 + i \partial_{y}\mathbf{I}_2\]

\[\hat{\mathbf{G}_1}, \hat{\mathbf{G}_2} = \mathbf{G}_1 / \| \mathbf{G}_1\|,\quad \mathbf{G}_2 / \| \mathbf{G}_2\|\]

\[\mathbf{S}_1, \mathbf{S}_2 = \mathcal{F}[\hat{\mathbf{G}_1}],\quad \mathcal{F}[\hat{\mathbf{G}_2}]\]

\[\mathbf{Q}_{12} = \mathbf{S}_1 \mathbf{S}_2^{\star}\]

where \(\partial_{x}\) is the spatial derivate in horizontal direction, \(\mathcal{F}\) denotes the Fourier transform and \(\star\) a complex conjugate operation and \(\hat{x}\) is a normalizing operator.

Schematically this looks like:

       ∂x I1┌-------┐
I1  ┌----┬--┤ 1/|x| ├--┬----┐ S1
 →--┤ ∂  |  ├-------┤  | F  ├--→-┐
    └----┴--┤ 1/|x| ├--┴----┘    |
       ∂y I1└-------┘            |
                                 | Q12 ┌------┐  C12
                                 ×--→--┤ F^-1 ├---
       ∂x I2┌-------┐            |     └------┘
I2  ┌----┬--┤ 1/|x| ├--┬----┐ S2 |
 →--┤ ∂  |  ├-------┤  | F  ├--→-┘
    └----┴--┤ 1/|x| ├--┴----┘
       ∂y I2└-------┘

References

Fi02

Fitch et al. “Orientation correlation”, Proceeding of the Britisch machine vison conference, pp. 1–10, 2002.

HK12

Heid & Kääb. “Evaluation of existing image matching methods for deriving glacier surface displacements globally from optical satellite imagery”, Remote sensing of environment, vol.118 pp.339-355, 2012.

AL22

Altena & Leinss, “Improved surface displacement estimation through stacking cross-correlation spectra from multi-channel imagery”, Science of Remote Sensing, vol.6 pp.100070, 2022.

Examples

>>> import numpy as np
>>> from dhdt.processing.matching_tools import get_integer_peak_location
>>> from dhdt.testing.matching_tools import create_sample_image_pair

>>> im1,im2,ti,tj,_ = create_sample_image_pair(d=2**4, max_range=1)
>>> Q = orientation_corr(im1, im2)
>>> C = np.fft.ifft2(Q)
>>> di,dj,_,_ = get_integer_peak_location(C)

>>> assert(np.isclose(ti, di, atol=1))
>>> assert(np.isclose(ti, di, atol=1))

dhdt.processing.matching_tools_frequency_correlators.phase_corr(I1, I2)

match two imagery through phase correlation

Parameters

• I1 (numpy.array, size=(m,n), ndim={2,3}) – array with intensities, if multiple bands are given the cross- correlation spectra are stacked, see [AL22].

• I2 (numpy.array, size=(m,n), ndim={2,3}) – array with intensities

Returns

Q – cross-power spectrum

Return type

numpy.array, size=(m,n), dtype=complex

See also

dhdt.processing.matching_tools_frequency_correlators.orientation_corr, dhdt.processing.matching_tools_frequency_correlators.cross_corr

Notes

The matching equations are as follows:

\[\mathbf{S}_1, \mathbf{S}_2 = \mathcal{F}[\mathbf{I}_1],\quad \mathcal{F}[\mathbf{I}_2]\]

\[\hat{\mathbf{S}}_1, \hat{\mathbf{S}}_2 = \mathbf{S}_1 / \| \mathbf{S}_1\|, \quad \mathbf{S}_2 / \| \mathbf{S}_2\|\]

\[\mathbf{Q}_{12} = \hat{\mathbf{S}}_1 {\hat{\mathbf{S}}_2}^{\star}\]

where \(\mathcal{F}\) denotes the Fourier transform and \(\star\) a complex conjugate operation and \(\hat{x}\) is a normalizing operator.

Schematically this looks like:

I1  ┌----┐ S1 ┌-------┐
 ---┤ F  ├----┤ 1/|x| ├--┐
    └----┘    └-------┘  | Q12 ┌------┐  C12
                         ×-----┤ F^-1 ├---
I2  ┌----┐ S2 ┌-------┐  |     └------┘
 ---┤ F  ├----┤ 1/|x| ├--┘
    └----┘    └-------┘

References

KH75

Kuglin & Hines. “The phase correlation image alignment method”, proceedings of the IEEE international conference on cybernetics and society, pp. 163-165, 1975.

AL22

Altena & Leinss, “Improved surface displacement estimation through stacking cross-correlation spectra from multi-channel imagery”, Science of Remote Sensing, vol.6 pp.100070, 2022.

Examples

>>> import numpy as np
>>> from dhdt.processing.matching_tools import get_integer_peak_location
>>> from dhdt.testing.matching_tools import create_sample_image_pair

>>> im1,im2,ti,tj,_ = create_sample_image_pair(d=2**4, max_range=1)
>>> Q = phase_corr(im1, im2)
>>> C = np.fft.ifft2(Q)
>>> di,dj,_,_ = get_integer_peak_location(C)

>>> assert(np.isclose(ti, di, atol=1))
>>> assert(np.isclose(ti, di, atol=1))

dhdt.processing.matching_tools_frequency_correlators.phase_only_corr(I1, I2)

match two imagery through phase only correlation, based upon [HG84] and [KJ91].

Parameters

• I1 (numpy.array, size=(m,n), ndim={2,3}) – array with intensities

• I2 (numpy.array, size=(m,n), ndim={2,3}) – array with intensities

Returns

Q – cross-power spectrum

Return type

numpy.array, size=(m,n), dtype=complex

See also

dhdt.processing.matching_tools_frequency_correlators.phase_corr, dhdt.processing.matching_tools_frequency_correlators.symmetric_phase_corr, dhdt.processing.matching_tools_frequency_correlators.amplitude_comp_corr

Notes

The matching equations are as follows:

\[\mathbf{S}_1, \mathbf{S}_2 = \mathcal{F}[\mathbf{I}_1],\quad \mathcal{F}[\mathbf{I}_2]\]

\[\mathbf{W} = 1 / \|\mathbf{S}_2 \|\]

\[\mathbf{Q}_{12} = \mathbf{S}_1 [\mathbf{W}\mathbf{S}_2]^{\star}\]

where \(\mathcal{F}\) denotes the Fourier transform and \(\star\) a complex conjugate operation

Schematically this looks like:

I2  ┌----┐ S2
 ---┤ F  ├--┬------------┐
    └----┘  | ┌-------┐W |
            └-┤ 1/|x| ├--×
              └-------┘  | Q12 ┌------┐  C12
                         ×-----┤ F^-1 ├---
I1  ┌----┐ S1            |     └------┘
 ---┤ F  ├---------------┘
    └----┘

If multiple bands are given, cross-spectral stacking is applied,see [AL22].

References

HG84

Horner & Gianino, “Phase-only matched filtering”, Applied optics, vol. 23(6) pp.812–816, 1984.

KJ91

Kumar & Juday, “Design of phase-only, binary phase-only, and complex ternary matched filters with increased signal-to-noise ratios for colored noise”, Optics letters, vol.16(13) pp.1025–1027, 1991.

AL22

Altena & Leinss, “Improved surface displacement estimation through stacking cross-correlation spectra from multi-channel imagery”, Science of Remote Sensing, vol.6 pp.100070, 2022.

Examples

>>> import numpy as np
>>> from dhdt.testing.matching_tools import create_sample_image_pair
>>> from dhdt.processing.matching_tools import get_integer_peak_location

>>> im1,im2,ti,tj,_ = create_sample_image_pair(d=2**4, max_range=1)
>>> Q = phase_only_corr(im1, im2)
>>> C = np.fft.ifft2(Q)
>>> di,dj,_,_ = get_integer_peak_location(C)

>>> assert(np.isclose(ti, di, atol=1))
>>> assert(np.isclose(ti, di, atol=1))

dhdt.processing.matching_tools_frequency_correlators.projected_phase_corr(I1, I2, M1=array([], dtype=float64), M2=array([], dtype=float64))

match two imagery through separated phase correlation

Parameters

• I1 ({numpy.array, numpy.ma}, size=(m,n), ndim=2) – array with intensities

• I2 ({numpy.array, numpy.ma}, size=(m,n), ndim=2) – array with intensities

• M1 (numpy.array, size=(m,n), ndim=2, dtype={bool,float}) – array with mask, hence False’s are neglected

• M2 (numpy.array, size=(m,n), ndim=2, dtype={bool,float}) – array with mask, hence False’s are neglected

Returns

C – displacement surface

Return type

numpy.array, size=(m,n), dtype=real

Notes

Schematically this looks like:

I1  ┌----┐ I1x  ┌----┐ S1x
 -┬-┤ Σx ├------┤ F  ├------┐ Q12x ┌------┐ C12x
  | └----┘      └----┘      ×------┤ F^-1 ├------┐
  | ┌----┐ I2x  ┌----┐ S2x  |      └------┘      |
  |┌┤ Σx ├------┤ F  ├------┘                    |  ┌---┐ C12
  ||└----┘      └----┘                           ×--┤ √ ├----
  ||┌----┐ I1y  ┌----┐ S1y                       |  └---┘
  └┼┤ Σy ├------┤ F  ├------┐ Q12y ┌------┐ C12y |
   |└----┘      └----┘      ×------┤ F^-1 ├------┘
I2 |┌----┐ I2y  ┌----┐ S2y  |      └------┘
 --┴┤ Σy ├------┤ F  ├------┘
    └----┘      └----┘

See also

dhdt.processing.matching_tools_frequency_correlators.phase_corr

References

Zh10

Zhang et al. “An efficient subpixel image registration based on the phase-only correlations of image projections”, IEEE proceedings of the 10th international symposium on communications and information technologies, 2010.

dhdt.processing.matching_tools_frequency_correlators.robust_corr(I1, I2)

match two imagery through fast robust correlation

Parameters

• I1 (numpy.array, size=(m,n), ndim=2) – array with intensities

• I2 (numpy.array, size=(m,n), ndim=2) – array with intensities

Returns

Q – cross-power spectrum

Return type

numpy.array, size=(m,n), dtype=complex

References

Fi05

Fitch et al. “Fast robust correlation”, IEEE transactions on image processing vol. 14(8) pp. 1063-1073, 2005.

Es07

Essannouni et al. “Adjustable SAD matching algorithm using frequency domain” Journal of real-time image processing, vol.1 pp.257-265, 2007.

Examples

>>> import numpy as np
>>> from dhdt.processing.matching_tools import get_integer_peak_location
>>> from dhdt.testing.matching_tools import create_sample_image_pair

>>> im1,im2,ti,tj,_ = create_sample_image_pair(d=2**4, max_range=1)
>>> Q = robust_corr(im1, im2)
>>> C = np.fft.ifft2(Q)
>>> di,dj,_,_ = get_integer_peak_location(C)

>>> assert(np.isclose(ti, di, atol=1))
>>> assert(np.isclose(ti, di, atol=1))

dhdt.processing.matching_tools_frequency_correlators.sign_only_corr(I1, I2)

match two imagery through phase only correlation

Parameters

• I1 (numpy.array, size=(m,n), ndim={2,3}) – array with intensities

• I2 (numpy.array, size=(m,n), ndim={2,3}) – array with intensities

Returns

C – displacement surface

Return type

numpy.array, size=(m,n), real

See also

dhdt.processing.matching_tools_frequency_correlators.cosine_corr

References

IK07

Ito & Kiya, “DCT sign-only correlation with application to image matching and the relationship with phase-only correlation”, IEEE international conference on acoustics, speech and signal processing, vol. 1, 2007.

dhdt.processing.matching_tools_frequency_correlators.symmetric_phase_corr(I1, I2)

match two imagery through symmetric phase only correlation (SPOF) also known as Smoothed Coherence Transform (SCOT)

Parameters

• I1 (numpy.array, size=(m,n), ndim={2,3}) – array with intensities, if multiple bands are given the cross- correlation spectra are stacked, see [3].

• I2 (numpy.array, size=(m,n), ndim={2,3}) – array with intensities

Returns

Q – cross-power spectrum

Return type

numpy.array, size=(m,n), dtype=complex

Notes

The matching equations are as follows:

\[\mathbf{S}_1, \mathbf{S}_2 = \mathcal{F}[\mathbf{I}_1], \quad \mathcal{F}[\mathbf{I}_2]\]

\[\mathbf{W} = 1 / \sqrt{||\mathbf{S}_1||||\mathbf{S}_2||}\]

\[\mathbf{Q}_{12} = \mathbf{S}_1 [\mathbf{W}\mathbf{S}_2]^{\star}\]

where \(\mathcal{F}\) denotes the Fourier transform and \(\star\) a complex conjugate operation

Schematically this looks like:

I2  ┌----┐ S2
 ---┤ F  ├---┬-------------------┐
    └----┘   | ┌--------------┐W |
             ├-┤ 1/(|x||x|)^½ ├--×
             | └--------------┘  | Q12 ┌------┐  C12
             |                   ×-----┤ F^-1 ├---
I1  ┌----┐ S1|                   |     └------┘
 ---┤ F  ├---┴-------------------┘
    └----┘

If multiple bands are given, cross-spectral stacking is applied, see [AL22].

References

NP93

Nikias & Petropoulou. “Higher order spectral analysis: a nonlinear signal processing framework”, Prentice hall. pp.313-322, 1993.

We05

Wernet. “Symmetric phase only filtering: a new paradigm for DPIV data processing”, Measurement science and technology, vol.16 pp.601-618, 2005.

AL22

Altena & Leinss, “Improved surface displacement estimation through stacking cross-correlation spectra from multi-channel imagery”, Science of Remote Sensing, vol.6 pp.100070, 2022.

Examples

>>> import numpy as np
>>> from dhdt.testing.matching_tools import create_sample_image_pair
>>> from dhdt.processing.matching_tools import get_integer_peak_location

>>> im1,im2,ti,tj,_ = create_sample_image_pair(d=2**4, max_range=1)
>>> Q = symmetric_phase_corr(im1, im2)
>>> C = np.fft.ifft2(Q)
>>> di,dj,_,_ = get_integer_peak_location(C)

>>> assert(np.isclose(ti, di, atol=1))
>>> assert(np.isclose(ti, di, atol=1))

dhdt.processing.matching_tools_frequency_correlators.upsampled_cross_corr(S1, S2, upsampling=2)

apply cros correlation, and upsample the correlation peak

Parameters

• S1 (numpy.array, size=(m,n), dtype=complex, ndim=2) – array with intensities

• S2 (numpy.array, size=(m,n), dtype=complex, ndim=2) – array with intensities

Returns

di,dj – sub-pixel displacement

Return type

float

See also

dhdt.processing.matching_tools_frequency_correlators.pad_dft, dhdt.processing.matching_tools_frequency_correlators.upsample_dft

References

GS08

Guizar-Sicairo, et al. “Efficient subpixel image registration algorithms”, Applied optics, vol. 33 pp.156–158, 2008.

Examples

>>> import numpy as np
>>> from dhdt.processing.matching_tools import get_integer_peak_location
>>> from dhdt.testing.matching_tools import create_sample_image_pair

>>> im1,im2,ti,tj,_ = create_sample_image_pair(d=2**4, max_range=1)
>>> di,dj = upsampled_cross_corr(im1, im2)

>>> assert(np.isclose(ti, di, atol=1))
>>> assert(np.isclose(ti, di, atol=1))

dhdt.processing.matching_tools_frequency_correlators.windrose_corr(I1, I2)

match two imagery through windrose phase correlation, see [KJ91].

Parameters

• I1 (numpy.array, size=(m,n), ndim={2,3}) – array with intensities, if multiple bands are given the cross- correlation spectra are stacked, see [AL22].

• I2 (numpy.array, size=(m,n), ndim={2,3}) – array with intensities

Returns

Q – cross-power spectrum

Return type

numpy.array, size=(m,n), dtype=complex

See also

dhdt.processing.matching_tools_frequency_correlators.binary_orientation_corr, dhdt.processing.matching_tools_frequency_correlators.orientation_corr, dhdt.processing.matching_tools_frequency_correlators.phase_only_corr

Notes

The matching equations are as follows:

\[\mathbf{S}_1, \mathbf{S}_2 = \mathcal{F}[\mathbf{I}_1],\quad \mathcal{F}[\mathbf{I}_2]\]

\[\mathbf{Q}_{12} = sign(\mathbf{S}_1) sign(\mathbf{S}_2)^{\star}\]

where \(\mathcal{F}\) denotes the Fourier transform and \(\star\) a complex conjugate operation

References

KJ91

Kumar & Juday, “Design of phase-only, binary phase-only, and complex ternary matched filters with increased signal-to-noise ratios for colored noise”, Optics letters, vol. 16(13) pp. 1025–1027, 1991.

AL22

Altena & Leinss, “Improved surface displacement estimation through stacking cross-correlation spectra from multi-channel imagery” Science of Remote Sensing, vol.6 pp.100070, 2022.

Examples

>>> import numpy as np
>>> from dhdt.processing.matching_tools import get_integer_peak_location
>>> from dhdt.testing.matching_tools import create_sample_image_pair

>>> im1,im2,ti,tj,_ = create_sample_image_pair(d=2**4, max_range=1)
>>> Q = windrose_corr(im1, im2)
>>> C = np.fft.ifft2(Q)
>>> di,dj,_,_ = get_integer_peak_location(C)

>>> assert(np.isclose(ti, di, atol=1))
>>> assert(np.isclose(ti, di, atol=1))

spatial correlators

dhdt.processing.matching_tools_spatial_correlators.cosine_similarity(I1, I2)

estimate the similarity of orientation imagery, via their dot product. Going along via different names, such as Gradient Inner Product (GIP)[Gl08]_ or Cosine Similarity [De21].

Parameters

• I1 (numpy.array, size=(m,n), dtype=complex, {x ∈ ℂ | ║x║ = 1}) – image with intensities (template), given as complex orientation image

• I2 (numpy.array, size=(k,l), dtype=complex, {x ∈ ℂ | ║x║ = 1}) – image with intensities (search space), as complex orientation image

Returns

score

Return type

numpy.array, dtype=float, range=0…+1

See also

dhdt.processing.matching_tools_frequency_correlators.orientation_corr

Examples

>>> import numpy as np
>>> from scipy import ndimage
>>> from dhdt.testing.matching_tools import create_sample_image_pair
>>> from dhdt.generic.handler_im import get_grad_filters
>>> from dhdt.processing.matching_tools import get_integer_peak_location
>>> from dhdt.processing.matching_tools_frequency_filters import             normalize_power_spectrum

>>> _,im2,ti,tj,im1 = create_sample_image_pair(d=2**4, max_range=1)

complex orientation imagery are used here, so these need to be transformed

>>> H_x = get_grad_filters(ftype='kroon', tsize=3, order=1)[0]
>>> H_y = np.transpose(H_x)
>>> G1 = ndimage.convolve(im2, H_x) +              1j*ndimage.convolve(im2, H_y)
>>> G2 = ndimage.convolve(im1_big, H_x) +              1j*ndimage.convolve(im1_big, H_y)
>>> G1, G2 = normalize_power_spectrum(G1), normalize_power_spectrum(G2)

now processing that actual processing can be done

>>> C = normalized_cross_corr(G1,G2)
>>> di,dj,_,_ = get_integer_peak_location(C)

>>> assert(np.isclose(ti, di, atol=1))
>>> assert(np.isclose(ti, di, atol=1))

References

Gl08

Glocker et al. “Optical flow estimation with uncertainties through dynamic MRFs” IEEE conference on computer vision and pattern recognition, 2008.

De21

Dematteis et al., “Comparison of digital image correlation methods and the impact of noise in geoscience applications” Remote sensing, vol.13(2), pp.327, 2021.

dhdt.processing.matching_tools_spatial_correlators.cumulative_cross_corr(I1, I2)

doing normalized cross correlation on distance imagery

Parameters

• I1 (numpy.array, ndim=2, size=(m,n)) – grid with intensities (template)

• I2 (numpy.array, ndim=2, size=(k,l), {k>=m, l>=n}) – grid with intensities (search space)

Returns

ccs – similarity surface

Return type

numpy.array

Notes

The following nomenclature is used:

ccs: cross correlation surface

dhdt.processing.matching_tools_spatial_correlators.maximum_likelihood(I1, I2)

maximum likelihood texture tracking

Parameters

• I1 (numpy.array, ndim=2, size=(m,n)) – grid with intensities (template)

• I2 (numpy.array, ndim=2, size=(k,l), {k>=m, l>=n}) – grid with intensities (search space)

Returns

C – texture correspondence surface

Return type

float

References

Er09

Erten et al. “Glacier velocity monitoring by maximum likelihood texture tracking” IEEE transactions on geosciences and remote sensing, vol.47(2) pp.394–405, 2009.

dhdt.processing.matching_tools_spatial_correlators.normalized_cross_corr(I1, I2)

simple normalized cross correlation

Parameters

• I1 (numpy.ndarray, ndim=2, size=(m,n)) – grid with intensities (template)

• I2 (numpy.ndarray, ndim=2, size=(k,l), {k>=m, l>=n}) – grid with intensities (search space)

Returns

ccs – similarity surface,

Return type

numpy.array

Examples

>>> import numpy as np
>>> from dhdt.processing.matching_tools import get_integer_peak_location
>>> from dhdt.testing.matching_tools import create_sample_image_pair

>>> _,im2,ti,tj,im1 = create_sample_image_pair(d=2**4, max_range=1)
>>> C = normalized_cross_corr(im1,im2)
>>> di,dj,_,_ = get_integer_peak_location(C)

>>> assert(np.isclose(ti, di, atol=1))
>>> assert(np.isclose(ti, di, atol=1))

Notes

The following nomenclature is used:

ccs: cross correlation surface

dhdt.processing.matching_tools_spatial_correlators.sum_sad_diff(I1, I2)

sum of absolute difference correlation

Parameters

• I1 (numpy.array, ndim=2, size=(m,n)) – grid with intensities (template)

• I2 (numpy.array, ndim=2, size=(k,l), {k>=m, l>=n}) – grid with intensities (search space)

Returns

sad – dissimilarity surface

Return type

numpy.array, ndim=2

Notes

The following nomenclature is used:

sad: sum of absolute difference

dhdt.processing.matching_tools_spatial_correlators.sum_sq_diff(I1, I2)

sum of squared difference correlation

Parameters

• I1 (numpy.array, ndim=2, size=(m,n)) – grid with intensities (template)

• I2 (numpy.array, ndim=2, size=(k,l), {k>=m, l>=n}) – grid with intensities (search space)

Returns

ssd – dissimilarity surface

Return type

numpy.array, ndim=2

Notes

The following nomenclature is used:

ssd: sum of squared differences

dhdt.processing.matching_tools_spatial_correlators.total_gradients(I1, I2)

calculate total gradients and its kurtosis, inspired by [Ch17]

Parameters

• I1 (numpy.ndarray, ndim=2, size=(m,n)) – grid with intensities (template)

• I2 (numpy.ndarray, ndim=2, size=(k,l), {k>=m, l>=n}) – grid with intensities (search space)

Returns

tg – similarity score surface

Return type

numpy.array, ndim=2

References

Ch17

Chen et al. “Normalized total gradients: a new measure for multispectral image registration” IEEE transactions on image processing, vol.27(3) pp.1297–1310, 2017.

dhdt.processing.matching_tools_spatial_correlators.weighted_normalized_cross_correlation(I1, I2, W1=None, W2=None)

Parameters

• I1 (numpy.array, ndim=2, size=(m,n)) – grid with intensities (template)

• I2 (numpy.array, ndim=2, size=(k,l), {k>=m, l>=n}) – grid with intensities (search space)

• W1 (numpy.array, size=(m,n), dtype={boolean,float}) – weighting matrix for the template

• W2 (numpy.array, size=(m,n), dtype={boolean,float}) – weighting matrix for the search range

Returns

• wncc (numpy.array) – weighted normalized cross-correlation

• support (numpy.array, range=0…1) – amount of overlap present and used in both imagery

References

AK21

Altena & Kääb “Quantifying river ice movement through a combination of European satellite monitoring services” International journal of applied Earth observation and geoinformation, vol.98 pp.102315, 2021.

phase angle estimators

class dhdt.processing.matching_tools_frequency_subpixel.PlaneModel

Least squares estimator for phase plane.

Vectors/lines are parameterized using polar coordinates as functional model:

z = x * dx + y * dy

This estimator minimizes the squared distances from all points to the plane, independent of distance.

estimate(data)

Estimate plane from data using least squares. :param data: N points with (x, y) coordinates of vector, respectively. :type data: (N, 3) array

Returns

success – True, if model estimation succeeds.

Return type

bool

predict_xy(xy, params=None)

Predict vector using the estimated heading. :param xy: x,y-coordinates. :type xy: numpy.array :param params: Optional custom parameter set. :type params: (2, ) array, optional

Returns

Q_hat – Predicted plane height at x,y-coordinates.

Return type

numpy.array

residuals(data)

Determine residuals of data to model For each point the shortest distance to the plane is returned.

Parameters

data ((N, 3) array) – N points with x, y coordinates and z values, respectively.

Returns

residuals – Residual for each data point.

Return type

(N, ) array

class dhdt.processing.matching_tools_frequency_subpixel.SawtoothModel

Least squares estimator for phase sawtooth.

estimate(data, params_bound=0)

Estimate plane from data using least squares. :param data: N points with (x, y) coordinates of vector, respectively. :type data: (N, 3) array :param params_bound: bound of the parameter space :type params_bound: float

Returns

success – True, if model estimation succeeds.

Return type

bool

predict_xy(xy, params=None)

Predict vector using the estimated heading. :param xy: x,y-coordinates. :type xy: numpy.array :param params: Optional custom parameter set. :type params: (2, ) array, optional

Returns

Q_hat – Predicted plane height at x,y-coordinates.

Return type

numpy.array

residuals(data)

Determine residuals of data to model For each point the shortest distance to the plane is returned.

Parameters

data ((N, 3) array) – N points with x, y coordinates and z values, respectively.

Returns

residuals – Residual for each data point.

Return type

(N, ) array

dhdt.processing.matching_tools_frequency_subpixel.phase_difference(Q, W=array([], dtype=float64))

get displacement from phase plane through neighbouring vector difference

find slope of the phase plane through local difference of the phase angles

Parameters

• Q (numpy.array, size=(m,n), dtype=complex) – normalized cross spectrum

• W (numpy.array, size=(m,n), dtype=boolean) – weigthing matrix

Returns

di,dj – sub-pixel displacement

Return type

float

See also

phase_tpss, phase_svd, phase_difference

References

Ka89

Kay, S. “A fast and accurate frequency estimator”, IEEE transactions on acoustics, speech and signal processing, vol.37(12) pp.1987-1990, 1989.

Examples

>>> import numpy as np
>>> from dhdt.testing.matching_tools import create_sample_image_pair

>>> im1,im2,ti,tj,_ = create_sample_image_pair(d=2**4, max_range=1)
>>> Q = phase_corr(im1, im2)
>>> di,dj,_,_ = phase_svd(Q)

>>> assert(np.isclose(ti, di, atol=.2))
>>> assert(np.isclose(tj, dj, atol=.2))

dhdt.processing.matching_tools_frequency_subpixel.phase_difference_1d(Q, W=array([], dtype=float64), axis=0)

get displacement from phase plane along one axis through differencing

find slope of the phase plane through local difference of the phase angles, see [Ka89]

Parameters

• Q (numpy.array, size=(m,n), dtype=complex) – normalized cross spectrum

• W (numpy.array, size=(m,n), dtype=boolean) – weigthing matrix

Returns

dj – sub-pixel displacement

Return type

float

See also

phase_tpss, phase_svd, phase_difference

References

Ka89

Kay, S. “A fast and accurate frequency estimator”, IEEE transactions on acoustics, speech and signal processing, vol.37(12) pp.1987-1990, 1989.

dhdt.processing.matching_tools_frequency_subpixel.phase_gradient_descend(data, W=array([], dtype=float64), x_0=array([0., 0.]), learning_rate=1, n_iters=50)

get phase plane of cross-spectrum through principle component analysis

find slope of the phase plane through principle component analysis

Parameters

• data (numpy.array, size=(m*n,3), dtype=complex) – normalized cross spectrum

• data – coordinate list with complex cross-sprectum at last collumn

• W (numpy.array, size=(m*n,1), dtype=boolean) – index of data that is correct

• W – list with classification of correct data

Returns

di,dj – sub-pixel displacement

Return type

float

See also

phase_lsq

Examples

>>> import numpy as np
>>> from dhdt.testing.matching_tools import create_sample_image_pair

>>> im1,im2,ti,tj,_ = create_sample_image_pair(d=2**5, max_range=1)
>>> Q = phase_corr(im1, im2)
>>> di,dj,_,_ = phase_gradient_descend(Q)

>>> assert(np.isclose(ti, di, atol=.2))
>>> assert(np.isclose(tj, dj, atol=.2))

dhdt.processing.matching_tools_frequency_subpixel.phase_jac(Q, m, W=array([], dtype=float64), F1=array([], dtype=float64), F2=array([], dtype=float64), rank=2)

Parameters

• Q (numpy.array, size=(_,_), dtype=complex) – cross spectrum

• m (numpy.array, size=(2,1), dtype=float) – displacement estimate, in pixel coordinate system

• W (numpy.array, size=(m,n), dtype=float | boolean) – weigthing matrix, in a range of 0…1

• F1 (np,array, size=(m,n), dtype=integer) – coordinate of the first axis from the Fourier spectrum.

• F2 (np,array, size=(m,n), dtype=integer) – coordinate of the second axis from the Fourier spectrum

• rank (TYPE, optional) – DESCRIPTION. The default is 2.

Returns

dQdm – Jacobian of phase estimate

Return type

numpy.array, size=(m,n)

dhdt.processing.matching_tools_frequency_subpixel.phase_lsq(data, W=array([], dtype=float64))

get phase plane of cross-spectrum through least squares plane fitting

find slope of the phase plane through principle component analysis

Parameters

• data (numpy.array, size=(m,n), dtype=complex) – normalized cross spectrum

• numpy.array (or) – coordinate list with complex cross-sprectum at last

• size=(m*n – coordinate list with complex cross-sprectum at last

• 3) – coordinate list with complex cross-sprectum at last

• dtype=complex – coordinate list with complex cross-sprectum at last

• W (numpy.array, size=(m,n), dtype=boolean) – index of data that is correct

• numpy.array – list with classification of correct data

• size=(m*n – list with classification of correct data

• 1) – list with classification of correct data

• dtype=boolean – list with classification of correct data

Returns

di,dj – sub-pixel displacement

Return type

float

See also

phase_pca, phase_ransac, phase_hough

Examples

>>> import numpy as np
>>> from dhdt.testing.matching_tools import create_sample_image_pair

>>> im1,im2,ti,tj,_ = create_sample_image_pair(d=2**5, max_range=1)
>>> Q = phase_corr(im1, im2)
>>> di,dj,_,_ = phase_lsq(Q)

>>> assert(np.isclose(ti, di, atol=.2))
>>> assert(np.isclose(tj, dj, atol=.2))

dhdt.processing.matching_tools_frequency_subpixel.phase_pca(data, W=array([], dtype=float64))

get phase plane of cross-spectrum through principle component analysis

find slope of the phase plane through principle component analysis

Parameters

• data (numpy.array, size=(m,n), dtype=complex) – normalized cross spectrum

• numpy.array (or) – coordinate list with complex cross-sprectum at last

• size=(m*n – coordinate list with complex cross-sprectum at last

• 3) – coordinate list with complex cross-sprectum at last

• dtype=complex – coordinate list with complex cross-sprectum at last

• W (numpy.array, size=(m,n), dtype=boolean) – index of data that is correct

• numpy.array – list with classification of correct data

• size=(m*n – list with classification of correct data

• 1) – list with classification of correct data

• dtype=boolean – list with classification of correct data

Returns

di,dj – sub-pixel displacement

Return type

float

See also

phase_lsq

Examples

>>> import numpy as np
>>> from dhdt.testing.matching_tools import create_sample_image_pair

>>> im1,im2,ti,tj,_ = create_sample_image_pair(d=2**5, max_range=1)
>>> Q = phase_corr(im1, im2)
>>> di,dj = phase_pca(Q)

>>> assert(np.isclose(ti, di, atol=.2))
>>> assert(np.isclose(tj, dj, atol=.2))

dhdt.processing.matching_tools_frequency_subpixel.phase_radon(Q, coord_system='ij')

get direction and magnitude from phase plane through Radon transform

find slope of the phase plane through single value decomposition

Parameters

Q (numpy.array, size=(m,n), dtype=complex) – cross spectrum

Returns

• * di,dj (float, default) – cartesian displacement

• * θ,ρ (float) – magnitude and direction of displacement

See also

phase_tpss, phase_svd, phase_difference

References

BF06

Balci & Foroosh. “Subpixel registration directly from the phase difference” EURASIP journal on advances in signal processing, pp.1-11, 2006.

Examples

>>> import numpy as np
>>> from dhdt.testing.matching_tools import create_sample_image_pair

>>> im1,im2,ti,tj,_ = create_sample_image_pair(d=2**5, max_range=1)
>>> Q = phase_corr(im1, im2)
>>> di,dj,_,_ = phase_radon(Q)

>>> assert(np.isclose(ti, di, atol=.2))
>>> assert(np.isclose(tj, dj, atol=.2))

dhdt.processing.matching_tools_frequency_subpixel.phase_ransac(data, max_displacement=1, precision_threshold=0.05)

robustly fit plane using RANSAC algorithm

find slope of the phase plane through random sampling and consensus

Parameters

• data (numpy.array, size=(m,n), dtype=complex) – normalized cross spectrum

• numpy.array (or) – coordinate list with complex cross-sprectum at last

• size=(m*n – coordinate list with complex cross-sprectum at last

• 3) – coordinate list with complex cross-sprectum at last

• dtype=complex – coordinate list with complex cross-sprectum at last

Returns

• di (float) – sub-pixel displacement along vertical axis

• dj (float) – sub-pixel displacement along horizontal axis

See also

phase_lsq, phase_svd, phase_hough, phase_pca

References

FB81

Fischler & Bolles. “Random sample consensus: a paradigm for model fitting with applications to image analysis and automated cartography” Communications of the ACM vol.24(6) pp.381-395, 1981.

To15

Tong et al. “A novel subpixel phase correlation method using singular value decomposition and unified random sample consensus” IEEE transactions on geoscience and remote sensing vol.53(8) pp.4143-4156, 2015.

Examples

>>> import numpy as np
>>> from dhdt.testing.matching_tools import create_sample_image_pair

>>> im1,im2,ti,tj,_ = create_sample_image_pair(d=2**5, max_range=1)
>>> Q = phase_corr(im1, im2)
>>> di,dj,_,_ = phase_ransac(Q)

>>> assert(np.isclose(ti, di, atol=.2))
>>> assert(np.isclose(tj, dj, atol=.2))

dhdt.processing.matching_tools_frequency_subpixel.phase_secant(data, W=array([], dtype=float64), x_0=array([0., 0.]))

get phase plane of cross-spectrum through secant

find slope of the phase plane through secant method (or Newton’s method) in multiple dimensions it is known as the Broyden’s method.

Parameters

• data (numpy.array, size=(m,n), dtype=complex) – normalized cross spectrum

• numpy.array (or) – coordinate list with complex cross-sprectum at last collumn

• size=(m*n – coordinate list with complex cross-sprectum at last collumn

• 3) – coordinate list with complex cross-sprectum at last collumn

• dtype=complex – coordinate list with complex cross-sprectum at last collumn

• W (numpy.array, size=(m,n), dtype=boolean) – index of data that is correct

• numpy.array – list with classification of correct data

• size=(m*n – list with classification of correct data

• 1) – list with classification of correct data

• dtype=boolean – list with classification of correct data

Returns

di,dj – sub-pixel displacement

Return type

float

See also

phase_gradient_descend

References

Br65

Broyden, “A class of methods for solving nonlinear simultaneous equations” Mathematics and computation. vol.19(92) pp.577–593, 1965.

Examples

>>> import numpy as np
>>> from dhdt.testing.matching_tools import create_sample_image_pair

>>> im1,im2,ti,tj,_ = create_sample_image_pair(d=2**5, max_range=1)
>>> Q = phase_corr(im1, im2)
>>> di,dj,_,_ = phase_secant(Q)

>>> assert(np.isclose(ti, di, atol=.2))
>>> assert(np.isclose(tj, dj, atol=.2))

dhdt.processing.matching_tools_frequency_subpixel.phase_slope_1d(t, rad=0.1)

estimate the slope and intercept for one-dimensional signal

Parameters

• t (numpy.array, size=(m,1), dtype=complex) – angle values.

• rad (float, range=(0.0,0.5)) – radial inclusion, seen from the center

Returns

x_hat – estimated slope and intercept.

Return type

numpy.array, size=(2,1)

See also

phase_svd

dhdt.processing.matching_tools_frequency_subpixel.phase_svd(Q, W, rad=0.1)

get phase plane of cross-spectrum through single value decomposition

find slope of the phase plane through single value decomposition, see [Ho03]

Parameters

• Q (numpy.array, size=(m,n), dtype=complex) – cross spectrum

• W (numpy.array, size=(m,n), dtype=float) – weigthing matrix

• rad (float, range=(0.0,0.5)) – radial inclusion, seen from the center

Returns

di,dj – sub-pixel displacement

Return type

float

See also

phase_tpss, phase_radon, phase_difference

References

Ho03

Hoge, W.S. “A subspace identification extension to the phase correlation method”, IEEE transactions on medical imaging, vol.22(2) pp.277-280, 2003.

Examples

>>> import numpy as np
>>> from dhdt.testing.matching_tools import create_sample_image_pair

>>> im1,im2,ti,tj,_ = create_sample_image_pair(d=2**5, max_range=1)
>>> Q = phase_corr(im1, im2)
>>> di,dj,_,_ = phase_svd(Q)

>>> assert(np.isclose(ti, di, atol=.2))
>>> assert(np.isclose(tj, dj, atol=.2))

dhdt.processing.matching_tools_frequency_subpixel.phase_tpss(Q, W, m, p=0.0001, k=4, j=5, n=3)

get phase plane of cross-spectrum through two point step size iteration

find slope of the phase plane through two point step size for phase correlation minimization [BB88] which is implemented by [Le07]

Parameters

• Q (numpy.ndarray, size=(_,_), dtype=complex) – cross spectrum

• m0 (numpy.ndarray, size=(2,1)) – initial displacement estimate

• p (float, default=1e4) – closing error threshold

• k (integer, default=4) – number of refinements in iteration

• j (integer, default=5) – number of sub routines during an estimation

• n (integer, default=3) – mask convergence factor

Returns

• di,dj (float, size=(2,1)) – sub-pixel displacement

• snr (float) – signal-to-noise ratio

See also

phase_svd, phase_radon, phase_difference, phase_jac

References

BB88

Barzilai & Borwein. “Two-point step size gradient methods”, IMA journal of numerical analysis. vol.8 pp.141–148, 1988.

Le07

Leprince, et.al. “Automatic and precise orthorectification, coregistration, and subpixel correlation of satellite images, application to ground deformation measurements”, IEEE Transactions on geoscience and remote sensing vol. 45.6 pp. 1529-1558, 2007.

dhdt.processing.matching_tools_frequency_subpixel.phase_weighted_pca(Q, W)

get phase plane of cross-spectrum through principle component analysis

find slope of the phase plane through principle component analysis

Parameters

• data (numpy.array, size=(m,n), dtype=complex) – normalized cross spectrum

• numpy.array (or) – coordinate list with complex cross-sprectum at last

• size=(m*n – coordinate list with complex cross-sprectum at last

• 3) – coordinate list with complex cross-sprectum at last

• dtype=complex – coordinate list with complex cross-sprectum at last

• W (numpy.array, size=(m,n), dtype=boolean) – index of data that is correct

• numpy.array – list with classification of correct data

• size=(m*n – list with classification of correct data

• 1) – list with classification of correct data

• dtype=boolean – list with classification of correct data

Returns

di,dj – sub-pixel displacement

Return type

float

See also

phase_lsq

Examples

>>> import numpy as np
>>> from dhdt.testing.matching_tools import create_sample_image_pair

>>> im1,im2,ti,tj,_ = create_sample_image_pair(d=2**5, max_range=1)
>>> Q = phase_corr(im1, im2)
>>> W = gaussian_mask(Q)
>>> di,dj,_,_ = phase_weighted_pca(Q, W)

>>> assert(np.isclose(ti, di, atol=.2))
>>> assert(np.isclose(tj, dj, atol=.2))

dhdt.processing.matching_tools_frequency_subpixel.ransac(data, model_class, min_samples, residual_threshold, is_data_valid=None, is_model_valid=None, max_trials=100, stop_sample_num=inf, stop_residuals_sum=0, stop_probability=1, random_state=None, initial_inliers=None, params_bounds=0)

Fit a model to data with the RANSAC (random sample consensus) algorithm. RANSAC is an iterative algorithm for the robust estimation of parameters from a subset of inliers from the complete data set. Each iteration performs the following tasks: 1. Select min_samples random samples from the original data and check

whether the set of data is valid (see is_data_valid).

2. Estimate a model to the random subset (model_cls.estimate(*data[random_subset]) and check whether the estimated model is valid (see is_model_valid).

3. Classify all data as inliers or outliers by calculating the residuals to the estimated model (model_cls.residuals(*data)) - all data samples with residuals smaller than the residual_threshold are considered as inliers.

4. Save estimated model as best model if number of inlier samples is maximal. In case the current estimated model has the same number of inliers, it is only considered as the best model if it has less sum of residuals.

These steps are performed either a maximum number of times or until one of the special stop criteria are met. The final model is estimated using all inlier samples of the previously determined best model. :param data: Data set to which the model is fitted, where N is the number of data

points and the remaining dimension are depending on model requirements. If the model class requires multiple input data arrays (e.g. source and destination coordinates of skimage.transform.AffineTransform), they can be optionally passed as tuple or list. Note, that in this case the functions estimate(*data), residuals(*data), is_model_valid(model, *random_data) and is_data_valid(*random_data) must all take each data array as separate arguments.

Parameters

• model_class (object) –

  Object with the following object methods:

  • success = estimate(*data)

  • residuals(*data)

  where success indicates whether the model estimation succeeded (True or None for success, False for failure).

• min_samples (int in range (0, N)) – The minimum number of data points to fit a model to.

• residual_threshold (float larger than 0) – Maximum distance for a data point to be classified as an inlier.

• is_data_valid (function, optional) – This function is called with the randomly selected data before the model is fitted to it: is_data_valid(*random_data).

• is_model_valid (function, optional) – This function is called with the estimated model and the randomly selected data: is_model_valid(model, *random_data), .

• max_trials (int, optional) – Maximum number of iterations for random sample selection.

• stop_sample_num (int, optional) – Stop iteration if at least this number of inliers are found.

• stop_residuals_sum (float, optional) – Stop iteration if sum of residuals is less than or equal to this threshold.

• stop_probability (float in range [0, 1], optional) –

  RANSAC iteration stops if at least one outlier-free set of the training data is sampled with probability >= stop_probability, depending on the current best model’s inlier ratio and the number of trials. This requires to generate at least N samples (trials):

  N >= log(1 - probability) / log(1 - e**m)

  where the probability (confidence) is typically set to a high value such as 0.99, e is the current fraction of inliers w.r.t. the total number of samples, and m is the min_samples value.

• random_state ({None, int, numpy.random.Generator}, optional) – If random_state is None the numpy.random.Generator singleton is used. If random_state is an int, a new Generator instance is used, seeded with random_state. If random_state is already a Generator instance then that instance is used.

• initial_inliers (array-like of bool, shape (N,), optional) – Initial samples selection for model estimation

Returns

• model (object) – Best model with largest consensus set.

• inliers ((N, ) array) – Boolean mask of inliers classified as True.

References

1

“RANSAC”, Wikipedia, https://en.wikipedia.org/wiki/RANSAC

Examples

Generate ellipse data without tilt and add noise:

>>> t = np.linspace(0, 2 * np.pi, 50)
>>> xc, yc = 20, 30
>>> a, b = 5, 10
>>> x = xc + a * np.cos(t)
>>> y = yc + b * np.sin(t)
>>> data = np.column_stack([x, y])
>>> rng = np.random.default_rng(203560)  # do not copy this value
>>> data += rng.normal(size=data.shape)

Add some faulty data:

>>> data[0] = (100, 100)
>>> data[1] = (110, 120)
>>> data[2] = (120, 130)
>>> data[3] = (140, 130)

Estimate ellipse model using all available data:

>>> model = EllipseModel()
>>> model.estimate(data)
True
>>> np.round(model.params)  
array([ 72.,  75.,  77.,  14.,   1.])

Estimate ellipse model using RANSAC:

>>> ransac_model, inliers = ransac(
...     data, EllipseModel, 20, 3, max_trials=50
... )
>>> abs(np.round(ransac_model.params))
array([20., 30., 10.,  6.,  2.])

>>> inliers  
array([False, False, False, False,  True,  True,  True,  True,  True,
        True,  True,  True,  True,  True,  True,  True,  True,  True,
        True,  True,  True,  True,  True,  True,  True,  True,  True,
        True,  True,  True,  True,  True,  True,  True,  True,  True,
        True,  True,  True,  True,  True,  True,  True,  True,  True,
        True,  True,  True,  True,  True], dtype=bool)

>>> sum(inliers) > 40
True

RANSAC can be used to robustly estimate a geometric transformation. In this section, we also show how to use a proportion of the total samples, rather than an absolute number.

>>> from skimage.transform import SimilarityTransform
>>> rng = np.random.default_rng()
>>> src = 100 * rng.random((50, 2))
>>> model0 = SimilarityTransform(scale=0.5, rotation=1,
...                              translation=(10, 20))
>>> dst = model0(src)
>>> dst[0] = (10000, 10000)
>>> dst[1] = (-100, 100)
>>> dst[2] = (50, 50)
>>> ratio = 0.5  # use half of the samples
>>> min_samples = int(ratio * len(src))
>>> model, inliers = ransac((src, dst), SimilarityTransform, min_samples,
...                         10,
...                         initial_inliers=np.ones(len(src), dtype=bool))

>>> inliers  
array([False, False, False,  True,  True,  True,  True,  True,  True,
        True,  True,  True,  True,  True,  True,  True,  True,  True,
        True,  True,  True,  True,  True,  True,  True,  True,  True,
        True,  True,  True,  True,  True,  True,  True,  True,  True,
        True,  True,  True,  True,  True,  True,  True,  True,  True,
        True,  True,  True,  True,  True])

subpixel localization

dhdt.processing.matching_tools_spatial_subpixel.get_top_2d_gaussian(C, top=array([], dtype=float64))

find location of highest score using 2D Gaussian

Parameters

• C (numpy.array, size=(_,_)) – similarity surface

• top (numpy.array, size=(1,2)) – location of the maximum score

Returns

• ddi (float) – estimated subpixel location on the vertical axis of the peak

• ddj (float) – estimated subpixel location on the horizontal axis of the peak

• i_int (integer) – location of highest score on the vertical axis

• j_int (integer) – location of highest score on the horizontal axis

Notes

NH05

Nobach & Honkanen, “Two-dimensional Gaussian regression for sub-pixel displacement estimation in particle image velocimetry or particle position estimation in particle tracking velocimetry”, Experiments in fluids, vol.38 pp.511-515, 2005

dhdt.processing.matching_tools_spatial_subpixel.get_top_birchfield(C, top=array([], dtype=float64))

find location of highest score along each axis

Parameters

• C (numpy.array, size=(_,_)) – similarity surface

• top (numpy.array, size=(1,2)) – location of the maximum score

Returns

• ddi (float) – estimated subpixel location on the vertical axis of the peak

• ddj (float) – estimated subpixel location on the horizontal axis of the peak

• i_int (integer) – location of highest score on the vertical axis

• j_int (integer) – location of highest score on the horizontal axis

Notes

BT99

Birchfield & Tomasi. “Depth discontinuities by pixel-to-pixel stereo” International journal of computer vision, vol. 35(3) pp.269-293, 1999.

dhdt.processing.matching_tools_spatial_subpixel.get_top_blais(C, top=array([], dtype=float64))

find location of highest score through forth order filter

Parameters

• C (numpy.array, size=(_,_)) – similarity surface

• top (numpy.array, size=(1,2)) – location of the maximum score

Returns

• ddi (float) – estimated subpixel location on the vertical axis of the peak

• ddj (float) – estimated subpixel location on the horizontal axis of the peak

• i_int (integer) – location of highest score on the vertical axis

• j_int (integer) – location of highest score on the horizontal axis

Notes

BR86

Blais & Rioux, “Real-time numerical peak detector” Signal processing vol.11 pp.145-155, 1986.

dhdt.processing.matching_tools_spatial_subpixel.get_top_blue(C, ds=1, top=array([], dtype=float64))

find top of correlation peak through best linear unbiased estimation

References

1

dhdt.processing.matching_tools_spatial_subpixel.get_top_centroid(C, top=array([], dtype=float64))

find location of highest score through 1D centorid fit

Parameters

• C (numpy.ndarray, size=(_,_)) – similarity surface

• top (numpy.ndarray, size=(1,2)) – location of the maximum score

Returns

• ddi (float) – estimated subpixel location on the vertical axis of the peak

• ddj (float) – estimated subpixel location on the horizontal axis of the peak ddi : float estimated subpixel location on the vertical axis of the peak

• i_int (integer) – location of highest score on the vertical axis

• j_int (integer) – location of highest score on the horizontal axis

Notes

Ra18

Raffel et al. “Particle Image Velocimetry” Ch.6 pp.184, 2018.

dhdt.processing.matching_tools_spatial_subpixel.get_top_equiangular(C, top=array([], dtype=float64))

find location of highest score along each axis by equiangular line

Parameters

• C (numpy.array, size=(_,_)) – similarity surface

• top (numpy.array, size=(1,2)) – location of the crossing

Returns

• ddi (float) – estimated subpixel location on the vertical axis of the peak

• ddj (float) – estimated subpixel location on the horizontal axis of the peak

• i_int (integer) – location of highest score on the vertical axis

• j_int (integer) – location of highest score on the horizontal axis

Notes

SO05

Shimizu & Okutomi. “Sub-pixel estimation error cancellation on area-based matching” International journal of computer vision, vol.63(3), pp.207–224, 2005.

dhdt.processing.matching_tools_spatial_subpixel.get_top_esinc(C, ds=1, top=array([], dtype=float64))

find location of highest score using exponential esinc function

Parameters

• C (numpy.array, size=(_,_)) – similarity surface

• top (numpy.array, size=(1,2)) – location of the maximum score

Returns

• ddi (float) – estimated subpixel location on the vertical axis of the peak

• ddj (float) – estimated subpixel location on the horizontal axis of the peak

• i_int (integer) – location of highest score on the vertical axis

• j_int (integer) – location of highest score on the horizontal axis

Notes

AV06

Argyriou & Vlachos, “A study of sub-pixel motion estimation using phase correlation”, proceedings of the British machine vision conference, 2006

dhdt.processing.matching_tools_spatial_subpixel.get_top_gaussian(C, top=array([], dtype=float64))

find location of highest score through 1D gaussian fit

Parameters

• C (numpy.ndarray, size=(_,_)) – similarity surface

• top (numpy.ndarray, size=(1,2)) – location of the maximum score

Returns

• ddi (float) – estimated subpixel location on the vertical axis of the peak

• ddj (float) – estimated subpixel location on the horizontal axis of the peak ddi : float estimated subpixel location on the vertical axis of the peak

• i_int (integer) – location of highest score on the vertical axis

• j_int (integer) – location of highest score on the horizontal axis

Notes

WG91

Willert & Gharib, “Digital particle image velocimetry”, Experiments in fluids, vol.10 pp.181-193, 1991.

AV06

Argyriou & Vlachos, “A Study of sub-pixel motion estimation using phase correlation” Proceeding of the British machine vision conference, pp.387-396, 2006.

Ra18

Raffel et al. “Particle Image Velocimetry” Ch.6 pp.184 2018.

dhdt.processing.matching_tools_spatial_subpixel.get_top_mass(C, top=array([], dtype=float64))

find location of highest score through 1D center of mass

Parameters

• C (numpy.array, size=(_,_)) – similarity surface

• top (numpy.array, size=(1,2)) – location of the maximum score

Returns

• ddi (float) – estimated subpixel location on the vertical axis of the peak

• ddj (float) – estimated subpixel location on the horizontal axis of the peak

• i_int (integer) – location of highest score on the vertical axis

• j_int (integer) – location of highest score on the horizontal axis

Notes

FN96

Fisher & Naidu, “A Comparison of algorithms for subpixel peak detection” in Image Technology - Advances in image processing, multimedia and machine vision pp.385-404, 1996.

dhdt.processing.matching_tools_spatial_subpixel.get_top_moment(C, ds=1, top=array([], dtype=float64))

find location of highest score through bicubic fitting

Parameters

• C (numpy.ndarray, size=(_,_)) – similarity surface

• ds (integer, default=1) – size of the radius to use neighboring information

• top (numpy.ndarray, size=(1,2)) – location of the maximum score

Returns

• ddi (float) – estimated subpixel location on the vertical axis of the peak

• ddj (float) – estimated subpixel location on the horizontal axis of the peak

• i_int (integer) – location of highest score on the vertical axis

• j_int (integer) – location of highest score on the horizontal axis

References

Fe12

Feng et al. “A subpixel registration algorithm for low PSNR images” IEEE international conference on advanced computational intelligence, pp.626-630, 2012.

MG15

Messerli & Grinstad, “Image georectification and feature tracking toolbox: ImGRAFT” Geoscientific instrumentation, methods and data systems, vol.4(1) pp.23-34, 2015.

dhdt.processing.matching_tools_spatial_subpixel.get_top_parabolic(C, top=array([], dtype=float64), ds=1)

find location of highest score through 1D parabolic fit

Parameters

• C (numpy.array, size=(_,_)) – similarity surface

• top (numpy.array, size=(1,2)) – location of the maximum score

Returns

• ddi (float) – estimated subpixel location on the vertical axis of the peak

• ddj (float) – estimated subpixel location on the horizontal axis of the peak

• i_int (integer) – location of highest score on the vertical axis

• j_int (integer) – location of highest score on the horizontal axis

Notes

AV06

Argyriou & Vlachos, “A Study of sub-pixel motion estimation using phase correlation” Proceeding of the British machine vision conference, pp.387-396), 2006.

Ra18

Raffel et al. “Particle Image Velocimetry” Ch.6 pp.184 2018.

Br21

Bradley, “Sub-pixel registration for low cost, low dosage, X-ray phase contrast imaging”, 2021.

dhdt.processing.matching_tools_spatial_subpixel.get_top_paraboloid(C, top=array([], dtype=float64))

find location of highest score using paraboloid

Parameters

• C (numpy.array, size=(_,_)) – similarity surface

• top (numpy.array, size=(1,2)) – location of the maximum score

Returns

• ddi (float) – estimated subpixel location on the vertical axis of the peak

• ddj (float) – estimated subpixel location on the horizontal axis of the peak

• i_int (integer) – location of highest score on the vertical axis

• j_int (integer) – location of highest score on the horizontal axis

Notes

Pa20

Pallotta et al. “Subpixel SAR image registration through parabolic interpolation of the 2-D cross correlation”, IEEE transactions on geoscience and remote sensing, vol.58(6) pp.4132–4144, 2020.

dhdt.processing.matching_tools_spatial_subpixel.get_top_ren(C, top=array([], dtype=float64))

find location of highest score

Parameters

• C (numpy.array, size=(_,_)) – similarity surface

• top (numpy.array, size=(1,2)) – location of the maximum score

Returns

• ddi (float) – estimated subpixel location on the vertical axis of the peak

• ddj (float) – estimated subpixel location on the horizontal axis of the peak

• i_int (integer) – location of highest score on the vertical axis

• j_int (integer) – location of highest score on the horizontal axis

Notes

Re10

Ren et al. “High-accuracy sub-pixel motion estimation from noisy images in Fourier domain.” IEEE transactions on image processing, vol.19(5) pp.1379-1384, 2010.

dhdt.processing.matching_tools_spatial_subpixel.get_top_triangular(C, top=array([], dtype=float64))

find location of highest score through triangular fit

Parameters

• C (numpy.array, size=(_,_)) – similarity surface

• top (numpy.array, size=(1,2)) – location of the maximum score

Returns

• ddi (float) – estimated subpixel location on the vertical axis of the peak

• ddj (float) – estimated subpixel location on the horizontal axis of the peak

• i_int (integer) – location of highest score on the vertical axis

• j_int (integer) – location of highest score on the horizontal axis

Notes

OC91

Olsen & Coombs, “Real-time vergence control for binocular robots” International journal of computer vision, vol.7(1), pp.67-89, 1991.

frequency metrics

dhdt.processing.matching_tools_frequency_metrics.get_phase_metric(Q, di, dj, metric='phase_fit')

redistribution function, to get a metric from a matching score surface

Parameters

• Q (numpy.ndarray, size=(m,n), type=complex) – phase angles of the cross-correlation spectra

• di ({float,integer}) – estimated displacement

• dj ({float,integer}) – estimated displacement

• metric (string) – abbreviation for the metric type to be calculated, for the options see list_pahse_metrics() for the options

Returns

score – metric of the matching phase angle in relation to its surface

Return type

float

See also

list_phase_metrics, dhdt.processing.matching_tools_spatial_metrics.list_matching_metrics

dhdt.processing.matching_tools_frequency_metrics.list_phase_metrics()

list the abbreviations of the different implemented phase metrics, there are:

• 'phase_fit' : the absolute score

• 'phase_sup' : the primary peak ratio i.r.t. the second peak

See also

get_phase_metric

dhdt.processing.matching_tools_frequency_metrics.phase_fitness(Q, di, dj, norm=2)

estimate the goodness of fit of the phase diffence

Parameters

• Q (numpy.ndarray, size={(m,n), (k,3)}) – cross-spectrum

• di ({float,integer}) – estimated displacement

• dj ({float,integer}) – estimated displacement

• norm (integer) – norm for the difference

Returns

fitness – goodness of fit, see also [1]

Return type

float, range=0…1

See also

signal_to_noise, phase_support

References

Le07

Leprince, et.al. “Automatic and precise orthorectification, coregistration, and subpixel correlation of satellite images, application to ground deformation measurements”, IEEE Transactions on geoscience and remote sensing vol. 45.6 pp. 1529-1558, 2007.

dhdt.processing.matching_tools_frequency_metrics.phase_support(Q, di, dj, thres=1.4826)

estimate the support of the fit for the phase differences, following the implementation by [AL22].

Parameters

• Q (numpy.ndarray, size={(m,n), (k,3)}) – cross-spectrum

• di ({float,integer}) – estimated displacement

• dj ({float,integer}) – estimated displacement

• norm (integer) – norm for the difference

Returns

support – ammount of support, see also [1]

Return type

float, range=0…1

See also

phase_fitness

References

AL22

Altena & Leinss, “Improved surface displacement estimation through stacking cross-correlation spectra from multi-channel imagery”, Science of Remote Sensing, vol.6 pp.100070, 2022.

dhdt.processing.matching_tools_frequency_metrics.signal_to_noise(Q, C, norm=2)

calculate the signal to noise from a theoretical and an experimental cross-spectrum, see [Le07].

Parameters

• Q (numpy.ndarray, size=(m,n), dtype=complex) – cross-spectrum

• C (numpy.ndarray, size=(m,n), dtype=complex) – phase plane

• norm (integer) – norm for the difference

Returns

snr – signal to noise ratio, see also [1]

Return type

float, range=0…1

See also

phase_fitness

References

Le07

Leprince, et.al. “Automatic and precise orthorectification, coregistration, and subpixel correlation of satellite images, application to ground deformation measurements”, IEEE Transactions on geoscience and remote sensing vol. 45.6 pp. 1529-1558, 2007.

correlation metrics

dhdt.processing.matching_tools_spatial_metrics.entropy_corr(C)

metric for uniqueness of the correlation estimate

Parameters

C (numpy.ndarray, size=(m,n), dtype=float) – correlation or similarity surface

Returns

disorder – entropy of the correlation surface

Return type

float

See also

dhdt.processing.matching_tools_spatial_metrics.peak_confidence, dhdt.processing.matching_tools_spatial_metrics.peak_to_noise, dhdt.processing.matching_tools_spatial_metrics.peak_corr_energy, dhdt.processing.matching_tools_spatial_metrics.peak_rms_ratio, dhdt.processing.matching_tools_spatial_metrics.num_of_peaks, dhdt.processing.matching_tools_spatial_metrics.peak_winner_margin, dhdt.processing.matching_tools_spatial_metrics.primary_peak_margin, dhdt.processing.matching_tools_spatial_metrics.primary_peak_ratio

References

SS98

Scharstein & Szeliski, “Stereo matching with nonlinear diffusion” International journal of computer vision, vol.28(2) pp.155-174, 1998.

Xu14

Xue et al. “Particle image velocimetry correlation signal-to-noise ratio metrics and measurement uncertainty quantification” Measurement science and technology, vol.25 pp.115301, 2014.

St26

Sturges, “The choice of a class interval”. Journal of the american statistical association. vol.21(153) pp.65–66, 1926.

dhdt.processing.matching_tools_spatial_metrics.gauss_spread(C, intI, intJ, dI, dJ, est='dist')

estimate an oriented gaussian function through the vicinity of the correlation function. See for a full describtion [Al21].

Parameters

• C (numpy.ndarray, size=(m,n), dtype=float) – array with correlation values

• intI (integer, {x ∈ ℕ | x ≥ 0}) – location in rows of highest value in search space

• intJ (integer, {x ∈ ℕ | x ≥ 0}) – location in collumns of highest value in search space

• dI (float) – sub-pixel bias of top location along the row axis

• dJ (float) – sub-pixel bias of top location along the collumn axis

• est (string, default='dist') –

  • ‘dist’ do weighted least sqaures dependent on the distance from the top location

  otherwise : ordinary least squares

Returns

• cov_ii, cov_jj (float) – standard deviation of the vertical and horizontal axis

• ρ (float) – orientation of the Gaussian

• hess (numpy.ndarray, size=(4)) – estimate of the least squares computation

• frac – scaling of the correlation peak

Notes

The image frame is used in this function, hence the difference is:

coordinate |              +--------> collumns
system 'ij'|              |
           |              |
           |       j      |  image frame
   --------+-------->     |
           |              |
           |              v
           | i            rows
           v

References

Al21

Altena et al. “Correlation dispersion as a measure to better estimate uncertainty of remotely sensed glacier displacements” The cryosphere, vol.16(6) pp.2285-2300, 2021.

dhdt.processing.matching_tools_spatial_metrics.get_correlation_metric(C, metric='peak_abs')

redistribution function, to get a metric from a matching score surface

Parameters

• C (numpy.ndarray, size=(m,n)) – grid with correlation scores

• metric ({'peak_ratio', 'peak_rms', 'peak_ener', 'peak_nois', 'peak_conf',) – ‘peak_entr’, ‘peak_abs’, ‘peak_marg’, ‘peak_win’, ‘peak_num’} abbreviation for the metric type to be calculated, for the options see list_matching_metrics() for the options

Returns

score – metric of the matching peak in relation to its surface

Return type

float

See also

list_matching_metrics, dhdt.processing.matching_tools_frequency_metrics.list_phase_metrics

dhdt.processing.matching_tools_spatial_metrics.hessian_spread(C, intI, intJ)

Parameters

• C (numpy.ndarray, size=(m,n), dtype=float) – array with correlation values

• intI (integer, {x ∈ ℕ | x ≥ 0}) – location in rows of highest value in search space

• intJ (integer, {x ∈ ℕ | x ≥ 0}) – location in collumns of highest value in search space

Returns

• cov_ii (float, real, unit=pixel) – first diagonal element of the co-variance matrix

• cov_jj (float, real, unit=pixel) – second diagonal element of the co-variance matrix

• cov_ij (float, real, unit=pixel) – off-diagonal element of the co-variance matrix

Notes

The image frame is used in this function, hence the difference is:

coordinate |              +--------> collumns
system 'ij'|              |
           |              |
           |       j      |  image frame
   --------+-------->     |
           |              |
           |              v
           | i            rows
           v

References

Ro04

Rosen et al. “Updated repeat orbit interferometry package released” EOS, vol.85(5) pp.47, 2004.

Ca11

Casu et al. “Deformation time-series generation in areas characterized by large displacement dynamics: The SAR amplitude pixel-offset SBAS Technique” IEEE transactions in geoscience and remote sensing vol.49(7) pp.2752–2763, 2011.

dhdt.processing.matching_tools_spatial_metrics.intensity_disparity(I1, I2)

Parameters

• I1 (numpy.ndarray, size=(m,n), dtype=float) – template with intensities

• I2 (numpy.ndarray, size=(m,n), dtype=float) – template with intensities

Returns

sigma_dI – temporal intensity variation

Return type

float

References

Sc13

Sciacchitano et al. “PIV uncertainty quantification by image matching” Measurement science and technology, vol.24 pp.045302, 2013.

dhdt.processing.matching_tools_spatial_metrics.list_matching_metrics()

list the abbreviations of the different implemented correlation metrics

The following metrics are implemented:

• ‘peak_abs’ : the absolute score

• 'peak_ratio' : the primary peak ratio i.r.t. the second peak

• 'peak_rms' : the peak ratio i.r.t. the root mean square error

• 'peak_ener' : the peaks’ energy

• 'peak_noise' : the peak score i.r.t. to the noise level

• 'peak_conf' : the peak confidence

• 'peak_entr' : the peaks’ entropy

• 'peak_num' : the number of correlation peaks’

• 'peak_marg' : the dominance of the peak

See also

get_correlation_metric

dhdt.processing.matching_tools_spatial_metrics.num_of_peaks(C, filtering=True)

metric for the uniqueness of a match

Parameters

• C (numpy.ndarray, size=(m,n), dtype=float) – correlation or similarity surface

• filtering (boolean) – apply low-pass filtering, otherwise noise is also seen as a peak, see [2] for more motivation

Returns

ppm – margin between primary and secondary peak

Return type

float

See also

entropy_corr, peak_confidence, peak_to_noise, peak_corr_energy, peak_rms_ratio, peak_winner_margin, primary_peak_margin, primary_peak_ratio

References

Le07

Lefebvre et al., “A colour correlation-based stereo matching using 1D windows” IEEE conference on signal-image technologies and internet-based system, pp.702-710, 2007.

HM12

Hu & Mordohai, “A quantitative evaluation of confidence measures for stereo vision” IEEE transactions on pattern analysis and machine intelligence, vol.34(11) pp.2121-2133, 2012.

dhdt.processing.matching_tools_spatial_metrics.peak_confidence(C, radius=1)

metric for steepness of a correlation peak

Parameters

• C (numpy.ndarray, size=(m,n), dtype=float) – correlation or similarity surface

• radius (integer, {x ∈ ℕ | x ≥ 1}, default=1) – how far away should the sampling set be taken

Returns

q – confidence measure

Return type

float

See also

entropy_corr, peak_to_noise, peak_corr_energy, peak_rms_ratio, num_of_peaks, peak_winner_margin, primary_peak_margin, primary_peak_ratio

References

Er09

Erten et al. “Glacier velocity monitoring by maximum likelihood texture tracking” IEEE transactions on geosciences and remote sensing, vol.47(2) pp.394–405, 2009.

Er13

Erten “Glacier velocity estimation by means of polarimetric similarity measure” IEEE transactions on geosciences and remote sensing, vol.51 pp.3319–3327, 2013.

dhdt.processing.matching_tools_spatial_metrics.peak_corr_energy(C)

metric for uniqueness of the correlation estimate, also known as, “signal to noise”

Parameters

C (numpy.ndarray, size=(m,n), dtype=float) – correlation or similarity surface

Returns

prmsr – peak to correlation energy

Return type

float

See also

entropy_corr, peak_confidence, peak_to_noise, peak_rms_ratio, num_of_peaks, peak_winner_margin, primary_peak_margin, primary_peak_ratio

References

Xu14

Xue et al. “Particle image velocimetry correlation signal-to-noise ratio metrics and measurement uncertainty quantification” Measurement science and technology, vol.25 pp.115301, 2014.

dhdt.processing.matching_tools_spatial_metrics.peak_rms_ratio(C)

metric for uniqueness of the correlation estimate this is a similar metric as implemented in ROI_PAC [RoXX]

Parameters

C (numpy.ndarray, size=(m,n), dtype=float) – correlation or similarity surface

Returns

prmsr – peak to root mean square ratio

Return type

float

See also

entropy_corr, peak_confidence, peak_to_noise, peak_corr_energy, num_of_peaks, peak_winner_margin, primary_peak_margin, primary_peak_ratio

References

Xu14

Xue et al. “Particle image velocimetry correlation signal-to-noise ratio metrics and measurement uncertainty quantification” Measurement science and technology, vol.25 pp.115301, 2014.

Ro04

Rosen et al. “Updated repeat orbit interferometry package released” EOS, vol.85(5) pp.47, 2004.

dhdt.processing.matching_tools_spatial_metrics.peak_to_noise(C)

metric for uniqueness of the correlation estimate

Parameters

C (numpy.ndarray, size=(m,n), dtype=float) – correlation or similarity surface

Returns

snr – signal to noise ratio

Return type

float

See also

entropy_corr, peak_confidence, peak_corr_energy, peak_rms_ratio, num_of_peaks, peak_winner_margin, primary_peak_margin, primary_peak_ratio

References

dL07

de Lange et al. “Improvement of satellite radar feature tracking for ice velocity derivation by spatial frequency filtering” IEEE transactions on geosciences and remote sensing, vol.45(7) pp.2309–2318, 2007.

HK12

Heid & Kääb “Evaluation of existing image matching methods for deriving glacier surface displacements globally from optical satellite imagery.” Remote sensing of environment vol.118 pp.339-355, 2012.

dhdt.processing.matching_tools_spatial_metrics.peak_winner_margin(C)

metric for dominance of the correlation estimate in relation to other candidates and their surrounding scores. See also [SS98].

Parameters

C (numpy.ndarray, size=(m,n), dtype=float) – correlation or similarity surface

Returns

ppm – margin between primary and secondary peak

Return type

float

See also

entropy_corr, peak_confidence, peak_to_noise, peak_corr_energy, peak_rms_ratio, num_of_peaks, primary_peak_margin, primary_peak_ratio

References

SS98

Scharstein & Szeliski, “Stereo matching with nonlinear diffusion” International journal of computer vision, vol.28(2) pp.155-174, 1998.

dhdt.processing.matching_tools_spatial_metrics.primary_peak_margin(C)

metric for dominance of the correlation estimate in relation to other candidates

Parameters

C (numpy.ndarray, size=(m,n), dtype=float) – correlation or similarity surface

Returns

ppm – margin between primary and secondary peak

Return type

float

See also

entropy_corr, peak_confidence, peak_to_noise, peak_corr_energy, peak_rms_ratio, num_of_peaks, peak_winner_margin, primary_peak_ratio

References

HM12

Hu & Mordohai, “A quantitative evaluation of confidence measures for stereo vision” IEEE transactions on pattern analysis and machine intelligence, vol.34(11) pp.2121-2133, 2012.

dhdt.processing.matching_tools_spatial_metrics.primary_peak_ratio(C)

metric for uniqueness of the correlation estimate also known as, “Cross correlation peak ratio”

Parameters

C (numpy.ndarray, size=(m,n), dtype=float) – correlation or similarity surface

Returns

ppr – peak ratio between primary and secondary peak

Return type

float, range={0..1}

See also

entropy_corr, peak_confidence, peak_to_noise, peak_corr_energy, peak_rms_ratio, num_of_peaks, peak_winner_margin, primary_peak_margin

References

KA90

Keane & Adrian, “Optimization of particle image velocimeters. I. Double pulsed systems” Measurement science and technology. vol.1 pp.1202, 1990.

CV13

Charonk & Vlachos “Estimation of uncertainty bounds for individual particle image velocimetry measurements from cross-correlation peak ratio” Measurement science and technology. vol.24 pp.065301, 2013.

Xu14

Xue et al. “Particle image velocimetry correlation signal-to-noise ratio metrics and measurement uncertainty quantification” Measurement science and technology, vol.25 pp.115301, 2014.