diff --git a/.github/workflows/UnitTest.yml b/.github/workflows/UnitTest.yml
new file mode 100644
index 0000000..09dc64e
--- /dev/null
+++ b/.github/workflows/UnitTest.yml
@@ -0,0 +1,56 @@
+name: Unit test
+
+on:
+  create:
+    tags:
+  push:
+    branches:
+      - master
+  pull_request:
+
+jobs:
+  test:
+    runs-on: ${{ matrix.os }}
+    strategy:
+      fail-fast: false
+      matrix:
+        julia-version: ['1.3', '1', 'nightly']
+        os: [ubuntu-latest]
+        arch: [x64]
+        include:
+          - os: windows-latest
+            julia-version: '1'
+            arch: x64
+          - os: macOS-latest
+            julia-version: '1'
+            arch: x64
+          - os: ubuntu-latest
+            julia-version: '1'
+            arch: x86
+
+    steps:
+      - uses: actions/checkout@v2
+      - name: "Set up Julia"
+        uses: julia-actions/setup-julia@v1
+        with:
+          version: ${{ matrix.julia-version }}
+          arch: ${{ matrix.arch }}
+
+      - name: Cache artifacts
+        uses: actions/cache@v1
+        env:
+          cache-name: cache-artifacts
+        with:
+          path: ~/.julia/artifacts
+          key: ${{ runner.os }}-test-${{ env.cache-name }}-${{ hashFiles('**/Project.toml') }}
+          restore-keys: |
+            ${{ runner.os }}-test-${{ env.cache-name }}-
+            ${{ runner.os }}-test-
+            ${{ runner.os }}-
+      - name: "Unit Test"
+        uses: julia-actions/julia-runtest@master
+
+      - uses: julia-actions/julia-processcoverage@v1
+      - uses: codecov/codecov-action@v1
+        with:
+          file: lcov.info
diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml
new file mode 100644
index 0000000..b59bdbd
--- /dev/null
+++ b/.github/workflows/docs.yml
@@ -0,0 +1,41 @@
+name: Documentation
+
+on:
+  pull_request:
+  push:
+    branches:
+      - 'master'
+      - 'release-'
+    tags: '*'
+  release:
+    types: [published]
+
+jobs:
+  build:
+    runs-on: ${{ matrix.os }}
+    strategy:
+      matrix:
+        julia-version: [1]
+        os: [ubuntu-latest]
+    steps:
+      - uses: actions/checkout@v2
+      - uses: julia-actions/setup-julia@latest
+        with:
+          version: ${{ matrix.julia-version }}
+      - name: Cache artifacts
+        uses: actions/cache@v1
+        env:
+          cache-name: cache-artifacts
+        with:
+          path: ~/.julia/artifacts
+          key: ${{ runner.os }}-test-${{ env.cache-name }}-${{ hashFiles('**/Project.toml') }}
+          restore-keys: |
+            ${{ runner.os }}-test-${{ env.cache-name }}-
+            ${{ runner.os }}-test-
+            ${{ runner.os }}-
+      - name: Install dependencies
+        run: julia --project=docs/ -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()'
+      - name: Build and deploy
+        env:
+          GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
+        run: julia --project=docs/ docs/make.jl
diff --git a/.travis.yml b/.travis.yml
deleted file mode 100644
index a7cacc3..0000000
--- a/.travis.yml
+++ /dev/null
@@ -1,28 +0,0 @@
-language: julia
-os:
-  - linux
-  - osx
-julia:
-  - 1.3
-  - 1.4
-  - 1.5
-  - nightly
-notifications:
-  email: false
-git:
-  depth: 99999999
-matrix:
-  allow_failures:
-  - julia: nightly
-script:
-  - julia --color=yes -e 'using Pkg; Pkg.build()'
-  - julia --check-bounds=yes --color=yes -e 'using Pkg; Pkg.test(coverage=true)'
-jobs:
-  include:
-    - stage: "Documentation"
-      julia: 1.4
-      os: linux
-      script:
-        - julia --project=docs/ -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()'
-        - julia --project=docs/ docs/make.jl
-      after_success: skip
diff --git a/Project.toml b/Project.toml
index 5037ac9..7a27f67 100644
--- a/Project.toml
+++ b/Project.toml
@@ -4,28 +4,13 @@ authors = ["Peter Kovesi <peter.kovesi@gmail.com>"]
 version = "0.2.1"
 
 [deps]
-AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c"
-Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
-FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
-ImageMagick = "6218d12a-5da1-5696-b52f-db25d2ecc6d1"
-Images = "916415d5-f1e6-5110-898d-aaa5f9f070e0"
+ImageCore = "a09fc81d-aa75-5fe9-8630-4744c3626534"
+ImageMorphology = "787d08f9-d448-5407-9aad-5290dd7ab264"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
-Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
-TestImages = "5e47fb64-e119-507b-a336-dd2b206d9990"
 
 [compat]
-AbstractFFTs = "0.3, 0.4, 0.5, 1.0"
-Documenter = "0.24, 0.25, 0.26, 0.27"
 FFTW = "1"
-FileIO = "1"
-ImageMagick = "1"
-Images = "0.9, 0.10, 0.11, 0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 0.19, 0.20, 0.21, 0.22, 0.23, 0.24"
-TestImages = "1"
+ImageCore = "0.9"
+ImageMorphology = "0.3"
 julia = "1.3"
-
-[extras]
-Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
-
-[targets]
-test = ["Test"]
diff --git a/docs/Project.toml b/docs/Project.toml
new file mode 100644
index 0000000..4e34057
--- /dev/null
+++ b/docs/Project.toml
@@ -0,0 +1,12 @@
+[deps]
+DemoCards = "311a05b2-6137-4a5a-b473-18580a3d38b5"
+Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+ImageContrastAdjustment = "f332f351-ec65-5f6a-b3d1-319c6670881a"
+ImageIO = "82e4d734-157c-48bb-816b-45c225c6df19"
+ImageMagick = "6218d12a-5da1-5696-b52f-db25d2ecc6d1"
+ImagePhaseCongruency = "10e51d30-6ba1-539a-b97e-c69c597142c4"
+Images = "916415d5-f1e6-5110-898d-aaa5f9f070e0"
+
+[compat]
+Documenter = "0.27"
+DemoCards = "0.4"
diff --git a/docs/examples/config.json b/docs/examples/config.json
new file mode 100644
index 0000000..18fa980
--- /dev/null
+++ b/docs/examples/config.json
@@ -0,0 +1,11 @@
+{
+    "theme": "grid",
+    "order": [
+        "phase_congruency",
+        "test_images",
+        "misc"
+    ],
+    "properties": {
+        "notebook": "false"
+    }
+}
diff --git a/docs/examples/misc/perfft2.jl b/docs/examples/misc/perfft2.jl
new file mode 100644
index 0000000..f1797d4
--- /dev/null
+++ b/docs/examples/misc/perfft2.jl
@@ -0,0 +1,54 @@
+# ---
+# title: Fourier transform of Moisan periodic image component
+# id: demo_perfft2
+# cover: assets/perfft2.png
+# author: Peter Kovesi
+# date: 2018-10-26
+# ---
+
+# The function `perfft2()` implements Moisan's "Periodic plus Smooth Image
+# Decomposition" which decomposes an image into two components
+#
+#         img = p + s
+#
+# where `s` is the 'smooth' component with mean 0 and `p` is the 'periodic' component
+# which has no sharp discontinuities when one moves cyclically across the image
+# boundaries.
+#
+# This decomposition is very useful when one wants to obtain an FFT of an image
+# with minimal artifacts introduced from the boundary discontinuities. The image
+# `p` gathers most of the image information but avoids periodization artifacts.
+#
+# Reference:
+# L. Moisan, "Periodic plus Smooth Image Decomposition", Journal of
+# Mathematical Imaging and Vision, vol 39:2, pp. 161-179, 2011.
+
+using Images
+using FFTW
+using ImagePhaseCongruency
+using ImageContrastAdjustment
+using TestImages
+
+img = Float64.(Gray.(testimage("lena")))
+
+IMG = fft(img)               # 'Standard' fft
+(P, S, p, s) = perfft2(img)  # 'Periodic' fft
+
+mosaic(
+    adjust_histogram(Gray.(p), LinearStretching()),
+    adjust_histogram(s, LinearStretching()),
+    ## Note the vertical and horizontal cross in
+    ## the spectrum induced by the non-periodic edges.
+    adjust_histogram(log.(abs.(fftshift(IMG)) .+ 1), LinearStretching()),
+    ## Note the clean spectrum because p is periodic.
+    adjust_histogram(log.(abs.(fftshift(P)) .+ 1), LinearStretching());
+    nrow=2, rowmajor=true
+)
+# Top 1) left: periodic component 2) right: smooth component
+#
+# Bottom 3) left: spectrum of standard FFT 4) right: spectrum of periodic component
+
+# save cover image #src
+isdir("assets") || mkdir("assets") #src
+cover = Gray.(adjust_histogram(log.(abs.(fftshift(P)) .+ 1), LinearStretching())) #src
+save(joinpath("assets", "perfft2.png"), cover) #src
diff --git a/docs/examples/phase_congruency/phasecong3.jl b/docs/examples/phase_congruency/phasecong3.jl
new file mode 100644
index 0000000..067a0ee
--- /dev/null
+++ b/docs/examples/phase_congruency/phasecong3.jl
@@ -0,0 +1,30 @@
+# ---
+# title: Log-Gabor filters v3
+# id: demo_phasecong3
+# cover: assets/phasecong3.png
+# author: Peter Kovesi
+# date: 2018-10-26
+# ---
+
+# Use of the function `phasecong3()` allows corner points to be detected as well. These
+# corner points are a subset of the edge image and, unlike other corner detectors, their
+# location is precise and stable over different scales.
+
+using TestImages
+using Images
+using ImagePhaseCongruency
+
+img = restrict(testimage("mandril_gray"))
+(edges, corners) = phasecong3(img)
+
+mosaic(
+    img,
+    adjust_histogram(Gray.(edges), LinearStretching()),
+    adjust_histogram(corners, LinearStretching()),
+    nrow=1
+)
+# Images from top to right: 1) original image 2) edges 3) corners
+
+# save cover image #src
+isdir("assets") || mkdir("assets") #src
+save(joinpath("assets", "phasecong3.png"), adjust_histogram(Gray.(edges), LinearStretching())) #src
diff --git a/docs/examples/phase_congruency/phasecongmono.jl b/docs/examples/phase_congruency/phasecongmono.jl
new file mode 100644
index 0000000..23c517b
--- /dev/null
+++ b/docs/examples/phase_congruency/phasecongmono.jl
@@ -0,0 +1,35 @@
+# ---
+# title: Monogenic filters
+# id: demo_phasecongmono
+# cover: assets/phasecongmono.png
+# author: Peter Kovesi
+# date: 2018-10-26
+# ---
+
+# Phase congruency marks all classes of features from steps to lines and is a dimensionless
+# quantity that ranges from 0 to 1. This allows fixed thresholds to be used over wide
+# classes of images.
+
+using TestImages
+using Images
+using ImagePhaseCongruency
+
+img = restrict(testimage("mandril_gray"))
+
+(pc, or, ft, T) = phasecongmono(img)
+nonmax = Images.thin_edges(pc, or)
+
+mosaic(
+    img,
+    adjust_histogram(pc, LinearStretching()),
+    nonmax,
+    hysthresh(nonmax, 0.1, 0.2);
+    nrow=2, rowmajor=true
+)
+
+# Images: 1) top left: original image 2) top right: phase congruency 3) bottom left:
+# non-maximal suppression 4) bottom right: Hystersis thresholded
+
+# save cover image #src
+isdir("assets") || mkdir("assets") #src
+save(joinpath("assets", "phasecongmono.png"), adjust_histogram(Gray.(pc), LinearStretching())) #src
diff --git a/docs/examples/phase_congruency/phasesymmetry.jl b/docs/examples/phase_congruency/phasesymmetry.jl
new file mode 100644
index 0000000..6e38a00
--- /dev/null
+++ b/docs/examples/phase_congruency/phasesymmetry.jl
@@ -0,0 +1,30 @@
+# ---
+# title: Symmetric monogenic filters
+# id: demo_phasesymmono
+# cover: assets/phasesymmono.gif
+# author: Peter Kovesi
+# date: 2018-10-26
+# ---
+
+# Phase symmetry responds well to line like features and circular objects.  The number of
+# filter scales will affect the scale of features that are marked. Phase symmetry marks
+# features independently of contrast (a bright circle is not more symmetric than a grey
+# circle) and is a dimensionless quantity between 0 and 1.  However this may not be what one
+# desires in which case the symmetry energy may be of greater interest.
+
+using TestImages
+using Images
+using ImagePhaseCongruency
+
+img = Gray.(testimage("blobs"))
+## Detect regions of bright symmetry (polarity = 1)
+phase_bright, = phasesymmono(img; nscale=5, polarity=1)
+
+## Detect regions of dark symmetry (polarity = -1)
+phase_dark, = phasesymmono(img; nscale=5, polarity=-1)
+
+mosaic(img, phase_bright, phase_dark; nrow=1)
+
+# save cover image #src
+isdir("assets") || mkdir("assets") #src
+save(joinpath("assets", "phasesymmono.gif"), Images.gif([phase_bright, phase_dark]); fps=1) #src
diff --git a/docs/examples/phase_congruency/ppdenoise.jl b/docs/examples/phase_congruency/ppdenoise.jl
new file mode 100644
index 0000000..0a7cf89
--- /dev/null
+++ b/docs/examples/phase_congruency/ppdenoise.jl
@@ -0,0 +1,33 @@
+# ---
+# title: Denoise
+# id: demo_ppdenoise
+# cover: assets/ppdenoise.png
+# author: Peter Kovesi
+# date: 2018-10-26
+# ---
+
+using TestImages
+using Images
+using ImageContrastAdjustment
+using ImagePhaseCongruency
+using Random #hide
+Random.seed!(1234) #hide
+
+## Values in the range 0 to 1
+img = centered(Gray.(restrict(testimage("lighthouse"))))[-127:128, -127:128]
+
+## Add noise with standard deviation of 0.25
+img .+= 0.25 * randn(size(img))
+
+cleanimg = ppdenoise(img; nscale=6, norient=6, mult=2.5, minwavelength=2, sigmaonf=0.55, dthetaonsigma=1.0, k=3, softness=1.0)
+
+mosaic(
+    adjust_histogram(img, LinearStretching()),
+    adjust_histogram(cleanimg, LinearStretching());
+    nrow=1
+)
+
+# save cover image #src
+isdir("assets") || mkdir("assets") #src
+cover = adjust_histogram(Gray.(cleanimg), LinearStretching()) #src
+save(joinpath("assets", "ppdenoise.png"), cover) #src
diff --git a/docs/examples/phase_congruency/ppdrc.jl b/docs/examples/phase_congruency/ppdrc.jl
new file mode 100644
index 0000000..b195387
--- /dev/null
+++ b/docs/examples/phase_congruency/ppdrc.jl
@@ -0,0 +1,42 @@
+# ---
+# title: Dynamic Range Compression
+# id: demo_ppdrc
+# cover: assets/ppdrc.png
+# author: Peter Kovesi
+# date: 2018-10-26
+# ---
+
+# An example using the 16 bit M51 image.  Phase preserving dynamic range compression allows
+# the scale of analysis to be controlled.  Here we process the image at wavelengths up to 50
+# pixels and up to 200 pixels.  Longer wavelengths allow larger structures to be seen. Small
+# wavelengths allow fine structures to be seen.  Note the image size is (510, 320).
+
+using TestImages
+using Images
+using ImageContrastAdjustment
+using ImagePhaseCongruency
+
+img = float64.(testimage("m51"))
+
+## Histogram equalization for reference (with a very large number of bins!)
+img_histeq = histeq(img, 100_000)
+
+## Phase presserving dynamic range compression at cutoff wavelengths of 50 and
+## 200 pixels.  Note we scale the image because its raw values are between 0 and
+## 1, see the help information for ppdrc() for details.
+scale = 1e4
+img_ppdrc1 = ppdrc(img*scale, 50)
+img_ppdrc2 = ppdrc(img*scale, 200)
+
+mosaic(
+    adjust_histogram(img, LinearStretching()),
+    adjust_histogram(img_histeq, LinearStretching()),
+    adjust_histogram(img_ppdrc1, LinearStretching()),
+    adjust_histogram(img_ppdrc2, LinearStretching()),
+    nrow=1
+)
+
+# save cover image #src
+isdir("assets") || mkdir("assets") #src
+cropped_cover = adjust_histogram(centered(img_ppdrc1)[-128:127, -128:127], LinearStretching()) #src
+save(joinpath("assets", "ppdrc.png"), cropped_cover) #src
diff --git a/docs/examples/phase_congruency/quantizephase.jl b/docs/examples/phase_congruency/quantizephase.jl
new file mode 100644
index 0000000..0a7556e
--- /dev/null
+++ b/docs/examples/phase_congruency/quantizephase.jl
@@ -0,0 +1,27 @@
+# ---
+# title: Phase Quantization
+# id: demo_quantizephase
+# cover: assets/quantizephase.gif
+# author: Peter Kovesi
+# date: 2018-10-26
+# ---
+
+# Phase values in an image are important.  However, despite this, phase can be quantized
+# very heavily with little perceptual loss.  It can be quantized to a few as four levels, or
+# even three.  Quantizing to two levels still gives an image that can be interpreted.
+
+using TestImages
+using Images
+using ImagePhaseCongruency
+
+img = Float64.(restrict(testimage("mandril_gray")))
+
+results = map((8, 4, 3, 2)) do nlevels
+    out = quantizephase(img, nlevels)
+    clamp01!(Gray.(out))
+end
+mosaic(results; nrow=1)
+
+# save cover image #src
+isdir("assets") || mkdir("assets") #src
+save(joinpath("assets", "quantizephase.gif"), Images.gif([results...]); fps=1) #src
diff --git a/docs/examples/phase_congruency/swapphase.jl b/docs/examples/phase_congruency/swapphase.jl
new file mode 100644
index 0000000..9f371d1
--- /dev/null
+++ b/docs/examples/phase_congruency/swapphase.jl
@@ -0,0 +1,32 @@
+# ---
+# title: Amplitude swapping
+# id: demo_swapphase
+# cover: assets/swapphase.gif
+# author: Peter Kovesi
+# date: 2018-10-26
+# ---
+
+# A demonstration of the importance of phase information in images. Given two
+# images`swapphase()` takes their Fourier transforms and constructs two new, synthetic,
+# images formed from the swapped phase and amplitude imformation.  In general it is the
+# phase information that dominates.  However, for textures where the amplitude spectra can
+# be concentrated in a limited set of locations, the reverse can apply.
+
+# See [Oppenheim and Lim's paper "The importance of phase in signals". Proceedings of the
+# IEEE. Volume: 69 , Issue: 5 , May 1981](https://ieeexplore.ieee.org/document/1456290)
+
+using TestImages
+using Images
+using ImagePhaseCongruency
+
+img1 = centered(Float64.(Gray.(restrict(testimage("lighthouse")))))[-127:128, -127:128]
+img2 = restrict(Float64.(testimage("mandril_gray")))[1:256, 1:256]
+
+(newimg1, newimg2) = swapphase(img1, img2)
+
+mosaic(Gray.(img1), newimg1, img2, newimg2; nrow=2)
+# Bottom 1) left: phase of lighthouse, amplitude of Mandrill 2) right: amplitude of lighthouse, phase of Mandrill
+
+# save cover image #src
+isdir("assets") || mkdir("assets") #src
+save(joinpath("assets", "swapphase.gif"), Images.gif([img1, newimg1]); fps=1) #src
diff --git a/docs/examples/test_images/circsine.jl b/docs/examples/test_images/circsine.jl
new file mode 100644
index 0000000..8035532
--- /dev/null
+++ b/docs/examples/test_images/circsine.jl
@@ -0,0 +1,19 @@
+# ---
+# title: circsine
+# id: demo_circsine
+# cover: assets/circsine.png
+# author: Peter Kovesi
+# date: 2018-10-26
+# ---
+
+using Images
+using ImagePhaseCongruency
+
+## Circular features at a phase congruent angle of pi/4 and
+## an amplitude decay exponent of 1.5
+img = circsine(offset = pi/4, ampexponent = -1.5)
+adjust_histogram(Gray.(img), LinearStretching())
+
+# save cover image #src
+isdir("assets") || mkdir("assets") #src
+save(joinpath("assets", "circsine.png"), adjust_histogram(Gray.(img), LinearStretching())) #src
diff --git a/docs/examples/test_images/noiseonf.jl b/docs/examples/test_images/noiseonf.jl
new file mode 100644
index 0000000..d0d3fee
--- /dev/null
+++ b/docs/examples/test_images/noiseonf.jl
@@ -0,0 +1,25 @@
+# ---
+# title: noiseonf
+# id: demo_noiseonf
+# cover: assets/noiseonf.png
+# author: Peter Kovesi
+# date: 2018-10-26
+# ---
+
+using Images
+using ImageContrastAdjustment
+using ImagePhaseCongruency
+
+## Noise images with amplitude decay exponents of 1.5 and 2.5
+img1 = noiseonf(512, 1.5)
+img2 = noiseonf(512, 2.5)
+
+mosaic(
+    adjust_histogram(Gray.(img1), LinearStretching()),
+    adjust_histogram(img2, LinearStretching());
+    nrow=1
+)
+
+# save cover image #src
+isdir("assets") || mkdir("assets") #src
+save(joinpath("assets", "noiseonf.png"), adjust_histogram(Gray.(img1), LinearStretching())) #src
diff --git a/docs/examples/test_images/starsine.jl b/docs/examples/test_images/starsine.jl
new file mode 100644
index 0000000..2e9778a
--- /dev/null
+++ b/docs/examples/test_images/starsine.jl
@@ -0,0 +1,19 @@
+# ---
+# title: starsine
+# id: demo_starsine
+# cover: assets/starsine.png
+# author: Peter Kovesi
+# date: 2018-10-26
+# ---
+
+using Images
+using ImagePhaseCongruency
+
+## Circular features at a phase congruent angle of pi/2 and
+## an amplitude decay exponent of 2
+img = starsine(offset = pi/4, ampexponent = -2)
+adjust_histogram(Gray.(img), LinearStretching())
+
+# save cover image #src
+isdir("assets") || mkdir("assets") #src
+save(joinpath("assets", "starsine.png"), adjust_histogram(Gray.(img), LinearStretching())) #src
diff --git a/docs/examples/test_images/step2line.jl b/docs/examples/test_images/step2line.jl
new file mode 100644
index 0000000..34a57a8
--- /dev/null
+++ b/docs/examples/test_images/step2line.jl
@@ -0,0 +1,43 @@
+# ---
+# title: step2line
+# id: demo_step2line
+# cover: assets/step2line.png
+# author: Peter Kovesi
+# date: 2018-10-26
+# ---
+
+# The `step2line()` function generates a phase congruent test image where angle at which the
+# congruency occurs is interpolated from 0 at the top of the image to pi/2 at the bottom.
+# This produces an interpolation of feature type from step edge to line.  The point being
+# that phase congruency at any angle produces a feature and the angle at which the
+# congruency occurs defines the feature type. Gradient based edge detectors will only
+# correctly mark the step-like feature towards the top of the image and incorrectly mark two
+# features towards the bottom of the image whereas phase congruency will correctly mark a
+# single feature from top to bottom.  In general, natural images contain a roughly uniform
+# distribution of the full continuum of feature types from step to line.
+
+using Images
+using ImagePhaseCongruency
+
+img1 = step2line(ampexponent=-1)
+## note the softer features
+img2 = step2line(ampexponent=-1.5)
+
+## Compute phase congruency on the `step2line` image using default parameters
+(pc,) = phasecongmono(step2line(ampexponent = -1))
+
+fimg = imfilter(step2line(ampexponent = -1), KernelFactors.gaussian((2, 2)))
+(gx, gy) = imgradients(fimg, KernelFactors.ando3)
+∇img = sqrt.(gx.^2 + gy.^2)
+
+mosaicview(
+    adjust_histogram(Gray.(img1), LinearStretching()),
+    adjust_histogram(img2, LinearStretching()),
+    adjust_histogram(pc, LinearStretching()),
+    adjust_histogram(∇img, LinearStretching()),
+    nrow=2, rowmajor=true
+)
+
+# save cover image #src
+isdir("assets") || mkdir("assets") #src
+save(joinpath("assets", "step2line.png"), adjust_histogram(Gray.(img1), LinearStretching())) #src
diff --git a/docs/make.jl b/docs/make.jl
index 8b9398c..0dd8f04 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -1,16 +1,29 @@
-#push!(LOAD_PATH, "../src/")
-
 using Documenter, ImagePhaseCongruency
+using TestImages
+using DemoCards
+
+testimage("cameraman") # used to trigger artifact downloading
 
+# generate
+demopage, postprocess_cb, demo_assets = makedemos("examples") # this is the relative path to docs/
+assets = []
+isnothing(demo_assets) || (push!(assets, demo_assets))
+
+format = Documenter.HTML(edit_link = "master",
+                         prettyurls = get(ENV, "CI", nothing) == "true",
+                         assets = assets)
 makedocs(
+    format=format,
     sitename = "ImagePhaseCongruency",
     pages = [
         "index.md",
-        "examples.md",
+        demopage,
         "functions.md"
     ]
 )
 
+postprocess_cb()
+
 deploydocs(
     repo = "github.com/peterkovesi/ImagePhaseCongruency.jl.git",
 )
diff --git a/docs/src/examples.md b/docs/src/examples.md
deleted file mode 100644
index 7b27113..0000000
--- a/docs/src/examples.md
+++ /dev/null
@@ -1,381 +0,0 @@
-# Examples
-
-Note that these examples use PyPlot for the output.  However I have had difficulty getting the automated documentation building process to handle PyPlot, accordingly all calls have been commented out.  If you want to execute these examples simple reinstate the PyPlot calls and all should be well.
-
-## [Phase Congruency](@id PhaseCongruencyExample)
-
-Phase congruency marks all classes of features from steps to lines and is a
-dimensionless quantity that ranges from 0 to 1.  This allows fixed
-thresholds to be used over wide classes of images.
-
-```@setup phasecong
-using TestImages, Images
-img = Float64.(Gray.(testimage("lighthouse")))[1:512, 1:512];
-```
-
-```@example phasecong
-using ImagePhaseCongruency, Images, TestImages #, PyPlot
-
-img = Float64.(Gray.(testimage("lighthouse")))[1:512, 1:512];
-#set_cmap(PyPlot.ColorMap("gray"))
-#imshow(img); axis("off")
-save("testimg.png", img) # hide
-
-(pc, or, ft, T) =
-         phasecongmono(img; nscale=4, minwavelength=3, mult=2,
-                        sigmaonf=0.55, k=3, cutoff=0.5, g=10,
-                        deviationgain=1.5, noisemethod=-1)
-
-#imshow(pc); axis("off")
-save("testimg_pc.png", imadjustintensity(pc)) # hide
-nonmax = thin_edges_nonmaxsup(pc, or)
-#imshow(nonmax); axis("off")
-save("testimg_nm.png", imadjustintensity(nonmax)) # hide
-# Hysteresis threshold between Phase Congruency of 0.1 and 0.2
-bw = hysthresh(nonmax, 0.1, 0.2)
-#imshow(bw); axis("off")
-save("testimg_bw.png", bw) # hide
-```
-
-|  Test image               | Phase Congruency          |
-|---------------------------|---------------------------|
-|![](testimg.png)           |![](testimg_pc.png)        |
-|**Non-maximal suppression**|**Hysteresis thresholded** |
-|![](testimg_nm.png)        |![](testimg_bw.png)        |
-
-
-Use of the function `phasecong3()` allows corner points to be detected as well.
-These corner points are a subset of the edge image and, unlike other corner
-detectors, their location is precise and stable over different scales.
-
-```@setup phasecong3
-using TestImages, Images
-img = Float64.(Gray.(testimage("lighthouse")))[1:512, 1:512];
-```
-
-```@example phasecong3
-using ImagePhaseCongruency, Images, TestImages #, PyPlot
-
-img = Float64.(Gray.(testimage("lighthouse")))[1:512, 1:512];
-#set_cmap(PyPlot.ColorMap("gray"))
-(M, m) = phasecong3(img)
-#imshow(M); axis("off")  # Edge image
-save("testimg_Me.png", imadjustintensity(M)) # hide
-#imshow(m); axis("off")  # 'Corner' image
-save("testimg_mc.png", imadjustintensity(m)) # hide
-```
-
-|  Test image edges         | Test image corners        |
-|---------------------------|---------------------------|
-|![](testimg_Me.png)        |![](testimg_mc.png)        |
-
-
-## [Phase Symmetry](@id PhaseSymmetryExample)
-
-Phase symmetry responds well to line like features and circular objects.  The
-number of filter scales will affect the scale of features that are marked.
-Phase symmetry marks features independently of contrast (a bright circle is not
-more symmetric than a grey circle) and is a dimensionless quantity between 0 and 1.  However this may not be what one desires in which case the symmetry energy
-may be of greater interest.
-
-```@setup phasesym
-using TestImages
-img = testimage("blobs")
-```
-
-```@example phasesym
-using ImagePhaseCongruency, Images, TestImages #, PyPlot
-
-img = Float64.(Gray.(testimage("blobs")))
-
-#set_cmap(PyPlot.ColorMap("gray"))
-#imshow(img); axis("off")
-save("blobs.png", img)  #hide
-
-# Detect regions of bright symmetry (polarity = 1)
-(phaseSym, symmetryEnergy, T) = phasesymmono(img; nscale=5, polarity=1)
-#imshow(phaseSym); axis("off")
-save("blobs_sym1.png", phaseSym) #hide
-
-# Detect regions of dark symmetry (polarity = -1)
-(phaseSym, symmetryEnergy, T) = phasesymmono(img; nscale=5, polarity=-1)
-#imshow(phaseSym); axis("off")
-save("blobs_sym-1.png", phaseSym) #hide
-
-```
-
-|  Blobs        | .                      |
-|--------------------|------------------------|
-|![](blobs.png) |                        |
-|**Bright symmetry** | **Dark Symmetry**      |
-|![](blobs_sym1.png) |![](blobs_sym-1.png) |
-
-
-## [Phase Preserving Dynamic Range Compression](@id ppdrcExample)
-
-An example using the 16 bit M51 image.  Phase preserving dynamic range
-compression allows the scale of analysis to be controlled.  Here we process the
-image at wavelengths up to 50 pixels and up to 200 pixels.  Longer wavelengths
-allow larger structures to be seen. Small wavelengths allow fine structures to
-be seen.  Note the image size is (510, 320).
-
-```@setup ppdrc
-using TestImages
-img = testimage("m51")
-```
-
-```@example ppdrc
-using ImagePhaseCongruency, TestImages, Images #, PyPlot
-
-#set_cmap(PyPlot.ColorMap("gray"))
-
-img = Float64.(testimage("m51"))
-#imshow(img)
-save("m51.png", imadjustintensity(img)) #hide
-
-# Histogram equalization for reference (with a very large number of bins!)
-#imshow(histeq(img, 100000))
-save("m51histeq.png", histeq(img, 100000)) #hide
-
-# Phase presserving dynamic range compression at cutoff wavelengths of 50 and
-# 200 pixels.  Note we scale the image because its raw values are between 0 and
-# 1, see the help information for ppdrc() for details.
-scale = 1e4
-#imshow(ppdrc(img*scale, 50))
-save("m51ppdrc50.png", imadjustintensity(ppdrc(img*scale, 50))) #hide
-#imshow(ppdrc(img*scale, 200))
-save("m51ppdrc200.png", imadjustintensity(ppdrc(img*scale, 200))) #hide
-```
-
-| M51                | Histogram equalized (100000 bins)|
-|--------------------|--------------------|
-|![](m51.png)        |![](m51histeq.png)  |
-|**ppdrc: wavelength cutoff 50 pixels**|**ppdrc: wavelength cutoff 200 pixels**|
-|![](m51ppdrc50.png) |![](m51ppdrc200.png)|
-
-
-## [Phase Preserving Denoising](@id ppdenoiseExample)
-
-```@setup denoise
-using TestImages, Images
-img = Float64.(Gray.(testimage("lighthouse")))[1:512, 1:512]
-```
-
-```@example denoise
-using ImagePhaseCongruency, TestImages, Images #, PyPlot
-
-img = Float64.(Gray.(testimage("lighthouse")))[1:512, 1:512]; # Values in the range 0 to 1
-img .+= 0.25 * randn(size(img))             # Add noise with standard deviation of 0.25
-
-cleanimg = ppdenoise(img,  nscale = 6, norient = 6, mult = 2.5, minwavelength = 2,
-                     sigmaonf = 0.55, dthetaonsigma = 1.0, k = 3, softness = 1.0)
-
-#set_cmap(PyPlot.ColorMap("gray"))
-#imshow(img)
-save("testimgplusnoise.png", imadjustintensity(img)) #hide
-#imshow(cleanimg)
-save("testimgdenoised.png", imadjustintensity(cleanimg)) #hide
-```
-
-| Test image + noise        | Test image denoised      |
-|---------------------------|--------------------------|
-| ![](testimgplusnoise.png) | ![](testimgdenoised.png) |
-
-## Phase-Amplitude Swapping
-
-A demonstration of the importance of phase information in images. Given two
-images`swapphase()` takes their Fourier transforms and constructs two new,
-synthetic, images formed from the swapped phase and amplitude imformation.  In
-general it is the phase information that dominates.  However, for textures where
-the amplitude spectra can be concentrated in a limited set of locations, the
-reverse can apply.
-
-See [Oppenheim and Lim's paper "The importance of phase in signals". Proceedings of the IEEE. Volume: 69 , Issue: 5 , May 1981](https://ieeexplore.ieee.org/document/1456290)
-
-```@setup phaseswap
-using TestImages, Images
-img = Float64.(Gray.(testimage("lighthouse")))[1:512, 1:512]
-img = testimage("mandril_gray")
-```
-
-```@example phaseswap
-using ImagePhaseCongruency, Images, TestImages #, PyPlot
-
-img1 = Float64.(Gray.(testimage("lighthouse")))[1:512, 1:512]
-img2 = Float64.(testimage("mandril_gray"))
-
-(newimg1, newimg2) = swapphase(img1, img2)
-
-#set_cmap(PyPlot.ColorMap("gray"))
-#imshow(newimg1)
-save("phaselighthouseampmandril.png", imadjustintensity(newimg1)) #hide
-#imshow(newimg2)
-save("amplighthousephasemandril.png", imadjustintensity(newimg2)) #hide
-```
-
-| Phase of lighthouse, amplitude of Mandrill | Amplitude of lighthouse, phase of Mandrill |
-|-------------------------------------------|-------------------------------------------|
-| ![](phaselighthouseampmandril.png)        | ![](amplighthousephasemandril.png)        |
-
-
-## Phase Quantization
-
-Phase values in an image are important.  However, despite this, phase can be
-quantized very heavily with little perceptual loss.  It can be quantized to a
-few as four levels, or even three.  Quantizing to two levels still gives an
-image that can be interpreted.
-
-```@setup phasequant
-using TestImages
-img = testimage("mandril_gray")
-```
-
-```@example phasequant
-using ImagePhaseCongruency, Images, TestImages #, PyPlot
-
-img = Float64.(testimage("mandril_gray"))
-
-#set_cmap(PyPlot.ColorMap("gray"))
-
-#imshow(quantizephase(img,8))
-save("testimg8.png", imadjustintensity(quantizephase(img,8))) #hide
-#imshow(quantizephase(img,4))
-save("testimg4.png", imadjustintensity(quantizephase(img,4))) #hide
-#imshow(quantizephase(img,3))
-save("testimg3.png", imadjustintensity(quantizephase(img,3))) #hide
-#imshow(quantizephase(img,2))
-save("testimg2.png", imadjustintensity(quantizephase(img,2))) #hide
-
-```
-| Mandrill: 8 phase values | Mandrill: 4 phase values |
-|:------------------------:|:------------------------:|
-|![](testimg8.png)         | ![](testimg4.png)        |
-| **Mandrill: 3 phase values** | **Mandrill: 2 phase values** |
-|![](testimg3.png)         | ![](testimg2.png)        |
-
-
-## Test Images
-
-The `step2line()` function generates a phase congruent test image where angle at
-which the congruency occurs is interpolated from 0 at the top of the image to
-pi/2 at the bottom. This produces an interpolation of feature type from step
-edge to line.  The point being that phase congruency at any angle produces a
-feature and the angle at which the congruency occurs defines the feature
-type. Gradient based edge detectors will only correctly mark the step-like
-feature towards the top of the image and incorrectly mark two features towards
-the bottom of the image whereas phase congruency will correctly mark a single
-feature from top to bottom.  In general, natural images contain a roughly
-uniform distribution of the full continuum of feature types from step to line.
-
-
-```@example
-using ImagePhaseCongruency, Images #, PyPlot
-
-#set_cmap(PyPlot.ColorMap("gray"))
-
-#imshow(step2line(ampexponent = -1))
-save("step2line-1.png", imadjustintensity(step2line(ampexponent = -1))) #hide
-#imshow(step2line(ampexponent = -1.5))  # Note the softer features
-save("step2line-15.png", imadjustintensity(step2line(ampexponent = -1.5))) #hide
-
-# Compute phase congruency on the step2line image using default parameters
-(pc,) = phasecongmono(step2line(ampexponent = -1))
-#imshow(pc)
-save("step2line_pc.png", imadjustintensity(pc)) #hide
-
-# Compute gradient magnitude of the step2line image
-fimg = imfilter(step2line(ampexponent = -1), KernelFactors.gaussian((2, 2)))
-(gx, gy) = imgradients(fimg, KernelFactors.ando3)
-#imshow(sqrt.(gx.^2 + gy.^2)) # Note the doubled responses at the bottom on the image.
-save("step2line_gr.png", imadjustintensity(sqrt.(gx.^2 + gy.^2))) #hide
-```
-
-| step2line ampexponent = -1 | step2line ampexponent = -1.5 |
-|--------------------|---------------------|
-|![](step2line-1.png)|![](step2line-15.png)|
-| **Phase Congruency on step2line ampexp=-1** | **Gradient magnitude of step2line: ampexp=-1** |
-|![](step2line_pc.png) | ![](step2line_gr.png) |   |
-
-
-```@example
-using ImagePhaseCongruency, Images #, PyPlot
-
-# Circular features at a phase congruent angle of pi/4
-# and an amplitude decay exponent of 1.5
-#imshow(circsine(offset = pi/4, ampexponent = -1.5))
-save("circsine.png", imadjustintensity(circsine(offset = pi/4, ampexponent = -1.5))) #hide
-
-# Radial features at a phase congruent angle of pi/2
-# and an amplitude decay exponent of 2
-#imshow(starsine(offset = pi/2, ampexponent = -2))
-save("starsine.png", imadjustintensity(starsine(offset = pi/2, ampexponent = -2))) #hide
-
-# Noise images with amplitude decay exponents of 1.5 and 2.5
-#imshow(noiseonf(512, 1.5))
-save("noiseonf_15.png", imadjustintensity(noiseonf(512, 1.5))) #hide
-#imshow(noiseonf(512, 2.5))
-save("noiseonf_25.png", imadjustintensity(noiseonf(512, 2.5))) #hide
-```
-
-| circsine           | starsine            |
-|--------------------|---------------------|
-|![](circsine.png)   |![](starsine.png)    |
-|**noiseonf: p=1.5** |**noiseonf: p=2.5**  |
-|![](noiseonf_15.png)|![](noiseonf_25.png) |
-
-
-## Fourier transform of Moisan's periodic image component
-
-The function `perfft2()` implements Moisan's "Periodic plus Smooth Image
-Decomposition" which decomposes an image into two components
-
-      img = p + s
-
-where `s` is the 'smooth' component with mean 0 and `p` is the 'periodic' component
-which has no sharp discontinuities when one moves cyclically across the image
-boundaries.
-
-This decomposition is very useful when one wants to obtain an FFT of an image
-with minimal artifacts introduced from the boundary discontinuities. The image
-`p` gathers most of the image information but avoids periodization artifacts.
-
-Reference:
-L. Moisan, "Periodic plus Smooth Image Decomposition", Journal of
-Mathematical Imaging and Vision, vol 39:2, pp. 161-179, 2011.
-
-```@setup perfft
-using TestImages, Images
-img = Float64.(Gray.(testimage("lighthouse")))[1:512, 1:512]
-```
-
-```@example perfft
-using ImagePhaseCongruency, TestImages, Images, FFTW #,PyPlot
-
-img = Float64.(Gray.(testimage("lighthouse")))[1:512, 1:512]
-
-IMG = fft(img)               # 'Standard' fft
-(P, S, p, s) = perfft2(img)  # 'Periodic' fft
-
-#set_cmap(PyPlot.ColorMap("gray"))
-
-#imshow(img)  # (img = p + s)
-#imshow(p)    # The periodic component
-save("testimg_p.png", imadjustintensity(p)) # hide
-#imshow(s)    # The smooth component
-save("testimg_s.png", imadjustintensity(s)) # hide
-
-#imshow(log.(abs.(fftshift(IMG)))) # Note the vertical and horizontal cross in
-                                  # the spectrum induced by the non-periodic edges.
-save("testimg_fft.png", imadjustintensity(log.(abs.(fftshift(IMG))))) # hide
-#imshow(log.(abs.(fftshift(P))))   # Note the clean spectrum because p is periodic.
-save("testimg_fft_p.png", imadjustintensity(log.(abs.(fftshift(P))))) # hide
-
-```
-
-|  Test image                  |                            |
-|------------------------------|----------------------------|
-| ![](testimg.png)             |                            |
-| **Periodic component**       | **Smooth component**       |
-| ![](testimg_p.png)           | ![](testimg_s.png)         |
-| **Spectrum of periodic component**| **Spectrum of standard FFT** |
-| ![](testimg_fft_p.png)       | ![](testimg_fft.png)          |
\ No newline at end of file
diff --git a/docs/src/index.md b/docs/src/index.md
index e67151e..2743ec6 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -1,6 +1,22 @@
 
 # ImagePhaseCongruency
 
+```@setup genimages
+using TestImages, Images, ImageContrastAdjustment
+using Random
+Random.seed!(1234)
+
+save("testimg.png", restrict(testimage("mandril_gray")))
+save("blobs.png", Gray.(testimage("blobs")))
+
+img = centered(Gray.(restrict(testimage("lighthouse"))))[-127:128, -127:128]
+img .+= 0.25 * randn(size(img))
+save("testimgplusnoise.png", clamp01!(img))
+
+img = testimage("m51")
+img = adjust_histogram(centered(img)[-128:127, -128:127], LinearStretching())
+save("m51.png", img)
+```
 
 This package provides a collection of image processing functions that exploit
 the importance of phase information in our perception of images. The functions
@@ -20,7 +36,7 @@ form two main groups:
 
 |.  |   |
 |---|---|
-|![](testimg.png)              |![](testimg_bw.png)           |
+|![](testimg.png)              |![](examples/covers/phasecongmono.png)           |
 
 Rather than assume a feature is a point of maximal intensity gradient, the Local
 Energy Model postulates that features are perceived at points in an image where
@@ -39,13 +55,13 @@ values to be used over large classes of images.
 
 * [`phasecongmono()`](@ref) Phase congruency of an image using monogenic filters.
 * [`phasecong3()`](@ref) Computes edge and corner phase congruency in an image via log-Gabor filters.
-* [Example](@ref PhaseCongruencyExample) of using `phasecongmono()` and `phasecong3()`.
+* [Example](@ref Phase-congruency) of using `phasecongmono()` and `phasecong3()`.
 
 ## Phase symmetry
 
 |.  |   |
 |---|---|
-|![](blobs.png) |![](blobs_sym-1.png) |
+|![](blobs.png) |![](examples/covers/phasesymmetry.gif) |
 
 A point of local symmetry in an image corresponds to a point where the local
 frequency components are at either the minimum or maximum points in their
@@ -55,14 +71,14 @@ quantity.
 
 * [`phasesym()`](@ref) Compute phase symmetry on an image via log-Gabor filters.
 * [`phasesymmono()`](@ref) Phase symmetry of an image using monogenic filters.
-* [Example](@ref PhaseSymmetryExample) of using `phasesymmono()`.
+* [Example](@ref demo_phasesymmono) of using `phasesymmono()`.
 
 
 ## Phase preserving denoising
 
 |.  |   |
 |---|---|
-| ![](testimgplusnoise.png) | ![](testimgdenoised.png) |
+| ![](testimgplusnoise.png) | ![](examples/covers/ppdenoise.png) |
 
 This is a wavelet denoising scheme that uses non-orthogonal, complex valued,
 log-Gabor wavelets, rather than the more usual orthogonal or bi-orthogonal
@@ -72,14 +88,14 @@ corrupted. It is also allows threshold values can be determined automatically
 from the statistics of the wavelet responses to the image.
 
 * [`ppdenoise()`](@ref) Phase preserving wavelet image denoising.
-* [Example](@ref ppdenoiseExample) of using `ppdenoise()`.
+* [Example](@ref demo_ppdenoise) of using `ppdenoise()`.
 
 
 ## Phase preserving dynamic range compression
 
 |.  |   |
 |---|---|
-|![](m51.png) |![](m51ppdrc200.png)|
+|![](m51.png) |![](examples/covers/ppdrc.png)|
 
 A common method for displaying images with a high dynamic range is to use some
 variant of histogram equalization.  The problem with histogram equalization is
@@ -92,7 +108,7 @@ information is preserved and the contrast amplification of structures in the
 signal is purely a function of their amplitude.
 
 * [`ppdrc()`](@ref) Phase Preserving Dynamic Range Compression.
-* [Example](@ref ppdrcExample) of using `ppdrc()`.
+* [Example](@ref demo_ppdrc) of using `ppdrc()`.
 
 
 ## Supporting filtering functions
@@ -150,4 +166,4 @@ Peter Kovesi, "Edges Are Not Just Steps". Proceedings of ACCV2002 The Fifth Asia
 
 Peter Kovesi, "Phase Preserving Denoising of Images". The Australian Pattern Recognition Society Conference: DICTA'99. December 1999. Perth WA. pp 212-217. [preprint](https://www.peterkovesi.com/papers/denoise.pdf)
 
-Peter Kovesi, "Phase Preserving Tone Mapping of Non-Photographic High Dynamic Range Images". Proceedings: The Australian Pattern Recognition Society Conference: Digital Image Computing: Techniques and Applications DICTA 2012. [preprint](https://www.peterkovesi.com/papers/DICTA2012-tonemapping.pdf)
\ No newline at end of file
+Peter Kovesi, "Phase Preserving Tone Mapping of Non-Photographic High Dynamic Range Images". Proceedings: The Australian Pattern Recognition Society Conference: Digital Image Computing: Techniques and Applications DICTA 2012. [preprint](https://www.peterkovesi.com/papers/DICTA2012-tonemapping.pdf)
diff --git a/src/phasecongruency.jl b/src/phasecongruency.jl
index c77b9b5..d0d8688 100755
--- a/src/phasecongruency.jl
+++ b/src/phasecongruency.jl
@@ -33,7 +33,8 @@ August 2015   Original conversion from MATLAB to Julia
 November 2017  Julia 0.6
 October  2018  Julia 0.7/1.0
 ---------------------------------------------------------------------=#
-using Images, FFTW, Statistics
+using FFTW, Statistics
+using ImageCore
 
 export phasecongmono, phasesymmono, ppdrc
 export highpassmonogenic, bandpassmonogenic
@@ -645,12 +646,10 @@ http://mathworld.wolfram.com/RayleighDistribution.html
 http://en.wikipedia.org/wiki/Rayleigh_distribution
 ```
 """
-function rayleighmode(data, nbins::Integer= 50)
-
-    edges, count = Images.imhist(data, nbins, 0, maximum(data))
-    ind = indmax(count)   # Find the index of maximum in histogram
-
-    return rmode = (edges[ind]+edges[ind+1])/2
+function rayleighmode(X, nbins::Integer= 50)
+    edges, counts = build_histogram(X, nbins=nbins)
+    ind = argmax(counts)
+    return (edges[ind]+edges[ind+1])/2
 end
 
 #-------------------------------------------------------------------------
diff --git a/src/utilities.jl b/src/utilities.jl
index e65a5a7..b3feca1 100644
--- a/src/utilities.jl
+++ b/src/utilities.jl
@@ -18,8 +18,7 @@ PK August 2015
    October 2018 Julia 0.7/1.0
 
 ---------------------------------------------------------------------=#
-import Images
-
+import ImageMorphology: label_components, component_indices
 export replacenan, fillnan
 export hysthresh, imgnormalize
 
@@ -126,8 +125,8 @@ function hysthresh(img::AbstractArray{T0,2}, T1::Real, T2::Real) where T0 <: Rea
 
     # Form 8-connected components of pixels with a value above the
     # lower threshold and get the indices of pixels in each component.
-    label = Images.label_components(img .>= T2, trues(3,3))
-    pix = Images.component_indices(label)
+    label = label_components(img .>= T2, trues(3,3))
+    pix = component_indices(label)
 
     # For each list of pixels in pix test to see if there are any
     # image values above T1.  If so, set these pixels in the output
diff --git a/test/Project.toml b/test/Project.toml
new file mode 100644
index 0000000..82a804d
--- /dev/null
+++ b/test/Project.toml
@@ -0,0 +1,9 @@
+[deps]
+FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
+ImageCore = "a09fc81d-aa75-5fe9-8630-4744c3626534"
+ImageIO = "82e4d734-157c-48bb-816b-45c225c6df19"
+ImageMagick = "6218d12a-5da1-5696-b52f-db25d2ecc6d1"
+ImageMorphology = "787d08f9-d448-5407-9aad-5290dd7ab264"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+TestImages = "5e47fb64-e119-507b-a336-dd2b206d9990"
diff --git a/test/test_frequencyfilt.jl b/test/test_frequencyfilt.jl
index d67a367..dd79b64 100644
--- a/test/test_frequencyfilt.jl
+++ b/test/test_frequencyfilt.jl
@@ -9,7 +9,7 @@ Set disp = true to display images for visual verification
 
 disp = false
 
-using ImagePhaseCongruency, Images, Test, TestImages, AbstractFFTs
+using ImagePhaseCongruency, ImageCore, Test, TestImages, FFTW
 
 if disp
     using PyPlot
diff --git a/test/test_phasecongruency.jl b/test/test_phasecongruency.jl
index 495b4c9..420b0d8 100644
--- a/test/test_phasecongruency.jl
+++ b/test/test_phasecongruency.jl
@@ -8,7 +8,7 @@ Set the variable 'disp' to true to display the processed images
 
 disp = false
 
-using Test, Images, ImagePhaseCongruency, TestImages
+using Test, ImageCore, ImagePhaseCongruency, TestImages
 
 if disp
     using PyPlot
diff --git a/test/test_syntheticimages.jl b/test/test_syntheticimages.jl
index a857074..83a489c 100644
--- a/test/test_syntheticimages.jl
+++ b/test/test_syntheticimages.jl
@@ -9,7 +9,7 @@ Set 'disp' = true to display the images
 
 disp = false
 
-using ImagePhaseCongruency, Images, TestImages
+using ImagePhaseCongruency, ImageCore, TestImages
 
 if disp
     using PyPlot