From dd01a2b70449503a0f0a02ce34b45f9c8a89f4ab Mon Sep 17 00:00:00 2001 From: John Walsh Date: Fri, 15 Aug 2025 16:44:35 -0400 Subject: [PATCH 1/6] switch to WHOI fork of phasepack containing optimized phase congruency function --- ifcb_features/phasecong.py | 393 ++----------------------------------- pyproject.toml | 3 +- requirements.txt | 3 +- 3 files changed, 23 insertions(+), 376 deletions(-) diff --git a/ifcb_features/phasecong.py b/ifcb_features/phasecong.py index 8570c69..dda186a 100644 --- a/ifcb_features/phasecong.py +++ b/ifcb_features/phasecong.py @@ -1,7 +1,5 @@ -# The following includes a modified version of the phasecong function from -# phasepack. - -# It skips several steps that are not used in IFCB segmentation +# The following uses the optimized phasecong function from phasepack v1.6.1+ +# with covariance_only=True for significant speedup in IFCB segmentation # Original license reproduced below @@ -26,370 +24,7 @@ # IFCB-specific optimizations and utilites by Joe Futrelle @ WHOI 2016 import numpy as np - -from scipy.fftpack import fftshift, ifftshift - -from phasepack.tools import rayleighmode as _rayleighmode -from phasepack.tools import lowpassfilter as _lowpassfilter - -# Try and use the faster Fourier transform functions from the pyfftw module if -# available -try: - from phasepack.tools import fft2, ifft2 -except ImportError: - # fall back to numpy (untested) - from numpy.fft import fft2, ifft2 - -def faster_phasecong(img, nscale=5, norient=6, minWaveLength=3, mult=2.1, - sigmaOnf=0.55, k=2., cutOff=0.5, g=10., noiseMethod=-1): - """ - Function for computing phase congruency on an image. This is a contrast- - invariant edge and corner detector. - Arguments: - ----------- - - img N/A The input image - nscale 5 Number of wavelet scales, try values 3-6 - norient 6 Number of filter orientations. - minWaveLength 3 Wavelength of smallest scale filter. - mult 2.1 Scaling factor between successive filters. - sigmaOnf 0.55 Ratio of the standard deviation of the Gaussian - describing the log Gabor filter's transfer function - in the frequency domain to the filter center - frequency. - k 2.0 No. of standard deviations of the noise energy - beyond the mean at which we set the noise threshold - point. You may want to vary this up to a value of - 10 or 20 for noisy images. - cutOff 0.5 The fractional measure of frequency spread below - which phase congruency values get penalized. - g 10 Controls the 'sharpness' of the transition in the - sigmoid function used to weight phase congruency - for frequency spread. - noiseMethod -1 Parameter specifies method used to determine - noise statistics. - -1 use median of smallest scale filter responses - -2 use mode of smallest scale filter responses - >=0 use this value as the fixed noise threshold - Returns: - --------- - M Maximum moment of phase congruency covariance, which can be used as - a measure of edge strength - m Minimum moment of phase congruency covariance, which can be used as - a measure of corner strength - - The convolutions are done via the FFT. Many of the parameters relate to - the specification of the filters in the frequency plane. The values do - not seem to be very critical and the defaults are usually fine. You may - want to experiment with the values of 'nscales' and 'k', the noise - compensation factor. - Notes on filter settings to obtain even coverage of the spectrum - sigmaOnf .85 mult 1.3 - sigmaOnf .75 mult 1.6 (filter bandwidth ~1 octave) - sigmaOnf .65 mult 2.1 - sigmaOnf .55 mult 3 (filter bandwidth ~2 octaves) - For maximum speed the input image should have dimensions that correspond - to powers of 2, but the code will operate on images of arbitrary size. - See also: phasecongmono, which uses monogenic filters for improved - speed, but does not return m, PC or EO. - References: - ------------ - Peter Kovesi, "Image Features From Phase Congruency". Videre: A Journal of - Computer Vision Research. MIT Press. Volume 1, Number 3, Summer 1999 - http://mitpress.mit.edu/e-journals/Videre/001/v13.html - Peter Kovesi, "Phase Congruency Detects Corners and Edges". Proceedings - DICTA 2003, Sydney Dec 10-12 - """ - - if img.dtype not in ['float32', 'float64']: - img = np.float64(img) - imgdtype = 'float64' - else: - imgdtype = img.dtype - - if img.ndim == 3: - img = img.mean(2) - - rows, cols = img.shape - - epsilon = 1E-4 # used to prevent /0. - IM = fft2(img) # Fourier transformed image - - zeromat = np.zeros((rows, cols), dtype=imgdtype) - - # matrices for covariance data - covx2 = zeromat.copy() - covy2 = zeromat.copy() - covxy = zeromat.copy() - - # Pre-compute some stuff to speed up filter construction - - # Set up X and Y matrices with ranges normalised to +/- 0.5 - if (cols % 2): - xvals = np.arange(-(cols - 1) / 2., - ((cols - 1) / 2.) + 1) / float(cols - 1) - else: - xvals = np.arange(-cols / 2., cols / 2.) / float(cols) - - if (rows % 2): - yvals = np.arange(-(rows - 1) / 2., - ((rows - 1) / 2.) + 1) / float(rows - 1) - else: - yvals = np.arange(-rows / 2., rows / 2.) / float(rows) - - x, y = np.meshgrid(xvals, yvals, sparse=True) - - # normalised distance from centre - radius = np.sqrt(x * x + y * y) - - # polar angle (-ve y gives +ve anti-clockwise angles) - theta = np.arctan2(-y, x) - - # Quadrant shift radius and theta so that filters are constructed with 0 - # frequency at the corners - radius = ifftshift(radius) # need to use ifftshift to bring 0 to (0,0) - theta = ifftshift(theta) - - # Get rid of the 0 radius value at the 0 frequency point (now at top-left - # corner) so that taking the log of the radius will not cause trouble. - radius[0, 0] = 1. - - sintheta = np.sin(theta) - costheta = np.cos(theta) - - del x, y, theta - - # Construct a bank of log-Gabor filters at different spatial scales - - # Filters are constructed in terms of two components. - # 1) The radial component, which controls the frequency band that the - # filter responds to - # 2) The angular component, which controls the orientation that the filter - # responds to. - # The two components are multiplied together to construct the overall - # filter. - - # Construct the radial filter components... First construct a low-pass - # filter that is as large as possible, yet falls away to zero at the - # boundaries. All log Gabor filters are multiplied by this to ensure no - # extra frequencies at the 'corners' of the FFT are incorporated as this - # seems to upset the normalisation process when calculating phase - # congrunecy. - - # Updated filter parameters 6/9/2013: radius .45, 'sharpness' 15 - lp = _lowpassfilter((rows, cols), .45, 15) - - logGaborDenom = 2. * np.log(sigmaOnf) ** 2. - logGabor = [] - - wavelengths = minWaveLength * mult ** np.arange(nscale) - - for wavelength in wavelengths: - # centre of frequency filter - fo = 1. / wavelength - - # log Gabor - logRadOverFo = np.log(radius / fo) - tmp = np.exp(-(logRadOverFo * logRadOverFo) / logGaborDenom) - - # apply low-pass filter - tmp = tmp * lp - - # set the value at the 0 frequency point of the filter back to - # zero (undo the radius fudge). - tmp[0, 0] = 0. - - logGabor.append(tmp) - - # MAIN LOOP - # for each orientation... - - # Construct the angular filter spread function - angls = np.arange(norient) * (np.pi / norient) - - for angl in angls: - # For each point in the filter matrix calculate the angular distance - # from the specified filter orientation. To overcome the angular wrap- - # around problem sine difference and cosine difference values are first - # computed and then the arctan2 function is used to determine angular - # distance. - - # difference in sine and cosine - sin_angl = np.sin(angl) - cos_angl = np.cos(angl) - - ds = sintheta * cos_angl - costheta * sin_angl - dc = costheta * cos_angl + sintheta * sin_angl - - # absolute angular difference - dtheta = np.abs(np.arctan2(ds, dc)) - - # Scale theta so that cosine spread function has the right wavelength - # and clamp to pi. - np.clip(dtheta * norient / 2., a_min=0, a_max=np.pi, out=dtheta) - - # The spread function is cos(dtheta) between -pi and pi. We add 1, and - # then divide by 2 so that the value ranges 0-1 - spread = (np.cos(dtheta) + 1.) / 2. - - # Initialize accumulators - sumE_ThisOrient = zeromat.copy() - sumO_ThisOrient = zeromat.copy() - sumAn_ThisOrient = zeromat.copy() - Energy = zeromat.copy() - - EOscale = [] - - # for each scale... - for ss in range(nscale): - - # Multiply radial and angular components to get filter - filt = logGabor[ss] * spread - - # Convolve image with even and odd filters - thisEO = ifft2(IM * filt) - - # Even + odd filter response amplitude - An = np.abs(thisEO) - - # Sum of amplitudes for even & odd filters across scales - sumAn_ThisOrient += An - - # Sum of even filter convolution results - sumE_ThisOrient += np.real(thisEO) - - # Sum of odd filter convolution results - sumO_ThisOrient += np.imag(thisEO) - - # At the smallest scale estimate noise characteristics from the - # distribution of the filter amplitude responses stored in sumAn. - # tau is the Rayleigh parameter that is used to describe the - # distribution. - if ss == 0: - # Use median to estimate noise statistics - if noiseMethod == -1: - tau = (np.median(sumAn_ThisOrient.ravel()) / - np.sqrt(np.log(4))) - - # Use the mode to estimate noise statistics - elif noiseMethod == -2: - tau = _rayleighmode(sumAn_ThisOrient.ravel()) - - maxAn = An - else: - # Record the maximum amplitude of components across scales to - # determine the frequency spread weighting - maxAn = np.maximum(maxAn, An) - - # append per-scale list - EOscale.append(thisEO) - - # Get weighted mean filter response vector, this gives the weighted - # mean phase angle. - XEnergy = np.sqrt(sumE_ThisOrient * sumE_ThisOrient + - sumO_ThisOrient * sumO_ThisOrient) + epsilon - MeanE = sumE_ThisOrient / XEnergy - MeanO = sumO_ThisOrient / XEnergy - - # Now calculate An(cos(phase_deviation)-| sin(phase_deviation))| by - # using dot and cross products between the weighted mean filter - # response vector and the individual filter response vectors at each - # scale. This quantity is phase congruency multiplied by An, which we - # call energy. - for ss in range(nscale): - E = np.real(EOscale[ss]) - O = np.imag(EOscale[ss]) - Energy += E * MeanE + O * MeanO - np.abs(E * MeanO - O * MeanE) - - # Automatically determine noise threshold - - # Assuming the noise is Gaussian the response of the filters to noise - # will form Rayleigh distribution. We use the filter responses at the - # smallest scale as a guide to the underlying noise level because the - # smallest scale filters spend most of their time responding to noise, - # and only occasionally responding to features. Either the median, or - # the mode, of the distribution of filter responses can be used as a - # robust statistic to estimate the distribution mean and standard - # deviation as these are related to the median or mode by fixed - # constants. The response of the larger scale filters to noise can then - # be estimated from the smallest scale filter response according to - # their relative bandwidths. - - # This code assumes that the expected reponse to noise on the phase - # congruency calculation is simply the sum of the expected noise - # responses of each of the filters. This is a simplistic overestimate, - # however these two quantities should be related by some constant that - # will depend on the filter bank being used. Appropriate tuning of the - # parameter 'k' will allow you to produce the desired output. - - # fixed noise threshold - if noiseMethod >= 0: - T = noiseMethod - - # Estimate the effect of noise on the sum of the filter responses as - # the sum of estimated individual responses (this is a simplistic - # overestimate). As the estimated noise response at succesive scales is - # scaled inversely proportional to bandwidth we have a simple geometric - # sum. - else: - totalTau = tau * (1. - (1. / mult) ** nscale) / (1. - (1. / mult)) - - # Calculate mean and std dev from tau using fixed relationship - # between these parameters and tau. See - # - EstNoiseEnergyMean = totalTau * np.sqrt(np.pi / 2.) - EstNoiseEnergySigma = totalTau * np.sqrt((4 - np.pi) / 2.) - - # Noise threshold, must be >= epsilon - T = np.maximum( - EstNoiseEnergyMean + k * EstNoiseEnergySigma, - epsilon) - - # Apply noise threshold, this is effectively wavelet denoising via soft - # thresholding. - Energy = np.maximum(Energy - T, 0) - - # Form weighting that penalizes frequency distributions that are - # particularly narrow. Calculate fractional 'width' of the frequencies - # present by taking the sum of the filter response amplitudes and - # dividing by the maximum amplitude at each point on the image. If - # there is only one non-zero component width takes on a value of 0, if - # all components are equal width is 1. - width = (sumAn_ThisOrient / (maxAn + epsilon) - 1.) / (nscale - 1) - - # Calculate the sigmoidal weighting function for this - # orientation - weight = 1. / (1. + np.exp(g * (cutOff - width))) - - # Apply weighting to energy, then calculate phase congruency - thisPC = weight * Energy / sumAn_ThisOrient - - # accumulate covariance data - covx = thisPC * cos_angl - covy = thisPC * sin_angl - covx2 += covx * covx - covy2 += covy * covy - covxy += covx * covy - - # Edge and Corner calculations - # The following is optimised code to calculate principal vectors of the - # phase congruency covariance data and to calculate the minimumum and - # maximum moments - these correspond to the singular values. - - # First normalise covariance values by the number of orientations/2 - covx2 /= norient / 2. - covy2 /= norient / 2. - covxy *= 4. / norient # This gives us 2*covxy/(norient/2) - denom = np.sqrt( - covxy * covxy + (covx2 - covy2) * (covx2 - covy2)) + epsilon - - # Maximum and minimum moments - M = (covx2 + covy2 + denom) / 2. - m = (covx2 + covy2 - denom) / 2. - - # don't compute anything else - - return M, m +from phasepack import phasecong # IFCB-specific utility @@ -403,9 +38,23 @@ def faster_phasecong(img, nscale=5, norient=6, minWaveLength=3, mult=2.1, PC_G=5 PC_NOISEMETHOD=-1 -def phasecong_Mm(roi, impl=faster_phasecong): - r = impl(roi,PC_NSCALE,PC_NORIENT,PC_MIN_WAVELENGTH,PC_MULT,PC_SIGMA_ONF,PC_K,PC_CUTOFF,PC_G,PC_NOISEMETHOD) - - M, m = r[0:2] +def phasecong_Mm(roi): + """ + IFCB-specific phase congruency function using optimized phasepack. + + Uses phasepack.phasecong with covariance_only=True and IFCB-optimized parameters. + Returns M + m (sum of maximum and minimum covariance moments). + """ + M, m = phasecong(roi, + nscale=PC_NSCALE, + norient=PC_NORIENT, + minWaveLength=PC_MIN_WAVELENGTH, + mult=PC_MULT, + sigmaOnf=PC_SIGMA_ONF, + k=PC_K, + cutOff=PC_CUTOFF, + g=PC_G, + noiseMethod=PC_NOISEMETHOD, + covariance_only=True) return M + m diff --git a/pyproject.toml b/pyproject.toml index 8510370..9fe2903 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -28,8 +28,7 @@ dependencies = [ "numpy", "scipy", "scikit-image", - "phasepack", - "pyfftw", + "phasepack @ git+https://github.com/WHOIGit/phasepack@v1.6.1", "scikit-learn", "pyifcb @ git+https://github.com/joefutrelle/pyifcb@v1.2.1" ] diff --git a/requirements.txt b/requirements.txt index 64ee293..0d1c87f 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,6 +1,5 @@ pandas==2.2.3 scipy==1.13.1 scikit-image==0.24.0 -phasepack==1.5 -pyfftw==0.15.0 +git+https://github.com/WHOIGit/phasepack@v1.6.1 scikit-learn==1.7.1 From 7f9fc613ec760e97560c62e52a9af717e268eb7a Mon Sep 17 00:00:00 2001 From: John Walsh Date: Mon, 18 Aug 2025 12:58:16 -0400 Subject: [PATCH 2/6] switch to uv to decrease image build time --- Dockerfile | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/Dockerfile b/Dockerfile index 3ae1249..cdec86c 100644 --- a/Dockerfile +++ b/Dockerfile @@ -1,4 +1,4 @@ -FROM python:3.12-slim +FROM ghcr.io/astral-sh/uv:python3.12-bookworm WORKDIR /app @@ -6,11 +6,8 @@ RUN apt-get update && apt-get install -y --no-install-recommends \ libgomp1 libopenblas0 libjpeg62-turbo tzdata git \ && rm -rf /var/lib/apt/lists/* -RUN pip install -U pip setuptools wheel - -COPY pyproject.toml ./ COPY . . -RUN pip install -v . +RUN uv pip install --system . ENTRYPOINT ["python", "extract_slim_features.py"] From e4d68c2fb7d24d9a0e51f7f7902a41debc165465 Mon Sep 17 00:00:00 2001 From: John Walsh Date: Mon, 18 Aug 2025 13:13:20 -0400 Subject: [PATCH 3/6] update github actions workflow to remove docker branch reference and push tagged images from prs --- .github/workflows/docker.yml | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/.github/workflows/docker.yml b/.github/workflows/docker.yml index a2b5c81..9ee8c92 100644 --- a/.github/workflows/docker.yml +++ b/.github/workflows/docker.yml @@ -2,7 +2,7 @@ name: Build and Push Docker Image on: push: - branches: [main, docker] + branches: [main] tags: ['v*'] pull_request: branches: [main] @@ -26,7 +26,6 @@ jobs: uses: docker/setup-buildx-action@v3 - name: Log in to Container Registry - if: github.event_name != 'pull_request' uses: docker/login-action@v3 with: registry: ${{ env.REGISTRY }} @@ -40,19 +39,18 @@ jobs: images: ${{ env.REGISTRY }}/${{ env.IMAGE_NAME }} tags: | type=ref,event=branch - type=ref,event=pr + type=ref,event=pr,prefix=pr- type=semver,pattern={{version}} type=semver,pattern={{major}}.{{minor}} type=sha type=raw,value=latest,enable={{is_default_branch}} - type=raw,value=docker,enable=${{ github.ref == 'refs/heads/docker' }} - name: Build and push Docker image uses: docker/build-push-action@v5 with: context: . platforms: linux/amd64,linux/arm64 - push: ${{ github.event_name != 'pull_request' }} + push: true tags: ${{ steps.meta.outputs.tags }} labels: ${{ steps.meta.outputs.labels }} cache-from: type=gha From 0392cb713c47e8bc3f75db4ecff6f84c457a20da Mon Sep 17 00:00:00 2001 From: John Walsh Date: Mon, 18 Aug 2025 14:08:32 -0400 Subject: [PATCH 4/6] remove system openblas install and use scipy openblas --- Dockerfile | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Dockerfile b/Dockerfile index cdec86c..5be2593 100644 --- a/Dockerfile +++ b/Dockerfile @@ -3,7 +3,7 @@ FROM ghcr.io/astral-sh/uv:python3.12-bookworm WORKDIR /app RUN apt-get update && apt-get install -y --no-install-recommends \ - libgomp1 libopenblas0 libjpeg62-turbo tzdata git \ + libjpeg62-turbo tzdata git \ && rm -rf /var/lib/apt/lists/* COPY . . From 0d2d026080747dcd087b6cf2a39c04b668a3a171 Mon Sep 17 00:00:00 2001 From: John Walsh Date: Mon, 18 Aug 2025 14:14:43 -0400 Subject: [PATCH 5/6] limit number of openblas threads --- Dockerfile | 2 ++ 1 file changed, 2 insertions(+) diff --git a/Dockerfile b/Dockerfile index 5be2593..c86510e 100644 --- a/Dockerfile +++ b/Dockerfile @@ -10,4 +10,6 @@ COPY . . RUN uv pip install --system . +ENV OPENBLAS_NUM_THREADS=64 OMP_NUM_THREADS=64 NUMEXPR_NUM_THREADS=64 + ENTRYPOINT ["python", "extract_slim_features.py"] From 752c2735cb6fd96e8273c6e25201471653bc942d Mon Sep 17 00:00:00 2001 From: John Walsh Date: Mon, 18 Aug 2025 14:54:59 -0400 Subject: [PATCH 6/6] remove libjpeg and tzdata installation --- Dockerfile | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/Dockerfile b/Dockerfile index c86510e..07278bd 100644 --- a/Dockerfile +++ b/Dockerfile @@ -2,8 +2,7 @@ FROM ghcr.io/astral-sh/uv:python3.12-bookworm WORKDIR /app -RUN apt-get update && apt-get install -y --no-install-recommends \ - libjpeg62-turbo tzdata git \ +RUN apt-get update && apt-get install -y --no-install-recommends git \ && rm -rf /var/lib/apt/lists/* COPY . .