OpenPIV / openpiv-python

OpenPIV is an open source Particle Image Velocimetry analysis software written in Python and Cython
http://www.openpiv.net
GNU General Public License v3.0
237 stars 141 forks source link

Clipping values below zero in normalize_intensity function #211

Closed ErichZimmer closed 2 years ago

ErichZimmer commented 3 years ago

Is your feature request related to a problem? Please describe. Some image data is lost when pixels less than zero are clipped. This could lead to an increase in RMS and bias errors.

Describe the solution you'd like Find a normalization method that keeps as much data as possible to extract as much information as possible.

Describe alternatives you've considered Don't do anything as it performs pretty good as is, especially for the zero-mean cross correlation algorithm implemented already.

Additional context Some random tests I did: bad_quality good_quality Nomenclature: corr = correlation int = image pixel intensities abs = absolute value

ErichZimmer commented 3 years ago

Using this function for normalized cross correlation,

def norm_corr(image_a, image_b, correlation_method = 'circular'):
    s1 = np.array(image_a.shape[-2:])
    s2 = np.array(image_b.shape[-2:])
    norm = (image_a.std(axis = (-2,-1))[:, np.newaxis, np.newaxis] * s2[0] *
            image_b.std(axis = (-2,-1))[:, np.newaxis, np.newaxis] * s2[1])
    if correlation_method == 'circular':
        f2a = rfft2(image_a - image_a.mean(axis = (-2,-1))[:,np.newaxis, np.newaxis]).conj()
        f2b = rfft2(image_b - image_b.mean(axis = (-2,-1))[:,np.newaxis, np.newaxis])
        corr = fftshift(irfft2(f2a * f2b).real, axes=(-2, -1)) / norm
    else:
        size = s1 + s2 - 1
        fsize = 2 ** np.ceil(np.log2(size)).astype(int)
        fslice = (slice(0, image_a.shape[0]),
                  slice((fsize[0]-s1[0])//2, (fsize[0]+s1[0])//2),
                  slice((fsize[1]-s1[1])//2, (fsize[1]+s1[1])//2))
        f2a = rfft2(image_a - image_a.mean(axis = (-2,-1))[:,np.newaxis, np.newaxis], fsize, axes=(-2, -1)).conj()
        f2b = rfft2(image_b - image_b.mean(axis = (-2,-1))[:,np.newaxis, np.newaxis], fsize, axes=(-2, -1))
        corr = fftshift(irfft2(f2a * f2b).real, axes=(-2, -1))[fslice] / norm
    return corr

the results seem a lot better. To deal with negative correlation values, the correlation was offset like such;

corr_min = corr2.min(axis = (-2,-1)) # avoid negative values for gaussian method
corr_min[corr_min < 0] = 0
u2, v2 = vectorized_correlation_to_displacements(corr2 - corr_min[:, np.newaxis, np.newaxis], n_rows, n_cols, 'gaussian')

The new normalized correlation algorithm seems to perform better. circular_good_quality circular_bad_quality linear_good_quality linear_bad_quality Settings: image size: 256x256 particle count: 10,000 window size: 32 overlap: 50%

alexlib commented 3 years ago

Looks a reasonable step to include it, @ErichZimmer please submit a pull request, starting from an updated forked master.

ErichZimmer commented 2 years ago

To avoid adding a new correlation function(s), I would suggest modifying the normalize_intensity function in a way that a user can choose to clip intensity values less than zero.

alexlib commented 2 years ago

To avoid adding a new correlation function(s), I would suggest modifying the normalize_intensity function in a way that a user can choose to clip intensity values less than zero.

Hi @ErichZimmer - please refresh my memory - what we're talking about ? is there a new pull request that modifies the correlation function? do we have a comparative test?

thanks, I'm just not sure that I know to recover the algorithm from 2021

ErichZimmer commented 2 years ago

@alexlib The current normalize_intensity function clips values less than zero. While this is good for noisy images, there is a slight loss of information on low noise images and an increase in outliers. For the normalized correlation function stated earlier in this issue, it would be unecessary to add it as it would produce the same results as fft_correlate_images when there is no clipping of negative pixel intensities during the normalization process.

alexlib commented 2 years ago

@ErichZimmer you mean this line https://github.com/OpenPIV/openpiv-python/blob/455ef8eaff225cd3d68edcdb2a995b63a9842005/openpiv/pyprocess.py#L729

to be optional with a flag in the Settings(), e.g. settings.clip_to_zero_intensity = False ?

ErichZimmer commented 2 years ago

Yes, that is correct. However, I am not sure if it goes against literature due to negative correlation values. To combat negative correlation values, I normally offset the correlation matrix by its lowest value and then perform subpixel estimation.

ErichZimmer commented 2 years ago

After doing a quick synthetic test on a new set of 1000 image pairs, I find that the clipping of negative pixel values makes the cross-correlation algorithm more stable, especially to varying particle intensities and diameters. The vector field may not look as "nice", but it is more accurate.