add all scripts/tools

5 files changed, 309 insertions, 0 deletions

diff --git a/julia/anisodiff.jl b/julia/anisodiff.jl
new file mode 100644
index 0000000..e01f6db
--- /dev/null
+++ b/julia/anisodiff.jl
@@ -0,0 +1,101 @@
+# -*- coding: utf-8 -*-
+#
+# ANISODIFF - Anisotropic diffusion
+#
+# Usage:
+#   diff = anisodiff(img, niter, kappa, lambda, option)
+#
+# Arguments:
+# | img    - input image (2D grayscale)
+# | niter  - number of iterations
+# | kappa  - conduction coefficient (gradient modulus threshold)
+# |          This parameter controls conduction as a function of gradient.
+# |          If kappa is low, small intensity gradients are able to block
+# |          conduction and hence diffusion across step edges.  A large value
+# |          reduces the influence of intensity gradients on conduction.
+# | lambda - integration constant for stability (0 <= lambda <= 0.25)
+# |          This parameter controls the diffusion speed, and you
+# |          usually want it at the maximum value of 0.25.
+# |          default value: 0.25
+# | option - conduction coefficient functions proposed by Perona & Malik:
+# |          1: c(x, y, t) = exp(-(nablaI/kappa).^2)
+# |             privileges high-contrast edges over low-contrast ones
+# |          2: c(x, y, t) = 1 ./ (1 + (nablaI/kappa).^2)
+# |             privileges wide regions over smaller ones
+# |          default value: 1
+#
+# Returns:
+# | diff   - anisotropic diffused image
+#
+# Reference:
+# [1] P. Perona and J. Malik.
+#     Scale-space and edge detection using ansotropic diffusion.
+#     IEEE Transactions on Pattern Analysis and Machine Intelligence,
+#     12(7):629-639, July 1990.
+#     https://dx.doi.org/10.1109%2F34.56205
+#
+# Credits:
+# [1] Peter Kovesi
+#     pk@peterkovesi.com
+#     MATLAB and Octave Functions for Computer Vision and Image Processing
+#     http://www.peterkovesi.com/matlabfns/Spatial/anisodiff.m
+#     --
+#     June 2000  original version
+#     March 2002 corrected diffusion eqn No 2.
+# [2] Daniel Lopes
+#     Anisotropic Diffusion (Perona & Malik)
+#     http://www.mathworks.com/matlabcentral/fileexchange/14995-anisotropic-diffusion--perona---malik-
+#
+#
+# Aaron LI <aaronly.me@gmail.com>
+# 2015/07/17
+#
+
+include("calc_k_percentile.jl");
+
+function anisodiff(img, niter, k=calc_k_percentile, lambda=0.25, option=1)
+    diff = float(img)
+    rows, cols = size(diff)
+
+    for i = 1:niter
+        println("anisodiff - iteration: ", i)
+
+        # Construct diffl which is the same as diff but
+        # has an extra padding of zeros around it.
+        diffl = zeros(rows+2, cols+2)
+        diffl[2:rows+1, 2:cols+1] = diff
+
+        # North, South, East and West differences
+        deltaN = diffl[1:rows,   2:cols+1] - diff
+        deltaS = diffl[3:rows+2, 2:cols+1] - diff
+        deltaE = diffl[2:rows+1, 3:cols+2] - diff
+        deltaW = diffl[2:rows+1, 1:cols]   - diff
+
+        # Calculate the kappa
+        if isa(k, Function)
+            kappa = k(diff)
+        else
+            kappa = k
+        end
+
+        println("  kappa: ", kappa)
+
+        # Conduction
+        if option == 1
+            cN = exp(-(deltaN/kappa).^2)
+            cS = exp(-(deltaS/kappa).^2)
+            cE = exp(-(deltaE/kappa).^2)
+            cW = exp(-(deltaW/kappa).^2)
+        elseif option == 2
+            cN = 1 ./ (1 + (deltaN/kappa).^2)
+            cS = 1 ./ (1 + (deltaS/kappa).^2)
+            cE = 1 ./ (1 + (deltaE/kappa).^2)
+            cW = 1 ./ (1 + (deltaW/kappa).^2)
+        end
+
+        diff += lambda * (cN.*deltaN + cS.*deltaS + cE.*deltaE + cW.*deltaW)
+    end
+
+    return diff
+end
+
diff --git a/julia/calc_k_percentile.jl b/julia/calc_k_percentile.jl
new file mode 100644
index 0000000..36b15e1
--- /dev/null
+++ b/julia/calc_k_percentile.jl
@@ -0,0 +1,32 @@
+# -*- coding: utf-8 -*-
+#
+# Calculate the percentile of the gradient image, which
+# used as the 'kappa' parameter of the anisotropic diffusion.
+#
+# Credits:
+# [1] KAZE: nldiffusion_functions.cpp / compute_k_percentile()
+#
+# Aaron LI
+# 2015/07/20
+#
+
+include("scharr.jl");
+
+function calc_k_percentile(img, percent=0.7, nbins=300)
+    rows, cols = size(img)
+    # derivatives of the image
+    img_gx = scharr(img, 1, 0)
+    img_gy = scharr(img, 0, 1)
+    img_modg = sqrt(img_gx.^2 + img_gy.^2)
+    # histogram
+    hmax = maximum(img_modg)
+    hist_e, hist_counts = hist(reshape(img_modg, length(img_modg)), nbins)
+    hist_cum = cumsum(hist_counts)
+    # find the percent of the histogram percentile
+    npoints = sum(img_modg .> 0.0)
+    nthreshold = npoints * percent
+    k = sum(hist_cum .<= nthreshold)
+    kperc = (k == length(hist_cum)) ? 0.03 : (hmax * k / nbins)
+    return kperc
+end
+
diff --git a/julia/forcefield.jl b/julia/forcefield.jl
new file mode 100644
index 0000000..bf2c236
--- /dev/null
+++ b/julia/forcefield.jl
@@ -0,0 +1,87 @@
+# -*- coding: utf-8 -*-
+#
+# Force field transform with specified size of mask.
+#
+# Aaron LI
+# 2015/07/19
+#
+
+# Make the specified sized force field mask.
+# NOTE: the number of rows and cols must be odd.
+function ff_mask(rows=5, cols=5)
+    rows % 2 == 1 || error("rows must be odd number")
+    cols % 2 == 1 || error("cols must be odd number")
+    mask = complex(zeros(rows, cols))
+    for r = range(-div(rows, 2), rows)
+        for c = range(-div(cols, 2), cols)
+            i, j = r + div(rows+1, 2), c + div(cols+1, 2)
+            #@printf("(r,c) = (%d,%d); (i,j) = (%d,%d)\n", r, c, i, j)
+            d = c + r*im
+            if abs(d) < 1e-8
+                mask[i, j] = 0.0
+            else
+                mask[i, j] = d / abs(d)^3
+            end
+        end
+    end
+    return mask / sum(abs(mask))
+end
+
+
+# Padding image by specified number of rows and cols.
+# Default padding mode: mirror
+function pad_image(img, pad_rows, pad_cols, mode="mirror")
+    rows, cols = size(img)
+    rows_new, cols_new = rows + 2*pad_rows, cols + 2*pad_cols
+    img_pad = zeros(rows_new, cols_new)
+    img_pad[(pad_rows+1):(pad_rows+rows), (pad_cols+1):(pad_cols+cols)] = img
+    for r = 1:rows_new
+        for c = 1:cols_new
+            if mode == "mirror"
+                if r <= pad_rows
+                    r_mirror = 2*(pad_rows+1) - r
+                elseif r <= pad_rows+rows
+                    r_mirror = r
+                else
+                    r_mirror = 2*(pad_rows+rows) - r
+                end
+                if c <= pad_cols
+                    c_mirror = 2*(pad_cols+1) - c
+                elseif c <= pad_cols+cols
+                    c_mirror = c
+                else
+                    c_mirror = 2*(pad_cols+cols) - c
+                end
+                if (r_mirror, c_mirror) != (r, c)
+                    #@printf("(%d,%d) <= (%d,%d)\n", r, c, r_mirror, c_mirror)
+                    img_pad[r, c] = img_pad[r_mirror, c_mirror]
+                end
+            else
+                error("mode not supported")
+            end
+        end
+    end
+    return img_pad
+end
+
+
+# Perform force field transform for the image.
+function ff_transform(img, mask, mode="mirror")
+    rows, cols = size(img)
+    mask_rows, mask_cols = size(mask)
+    pad_rows, pad_cols = div(mask_rows, 2), div(mask_cols, 2)
+    img_pad = pad_image(img, pad_rows, pad_cols)
+    # result images
+    ff_amplitudes = zeros(rows, cols)
+    ff_angles = zeros(rows, cols)
+    # calculate transformed values
+    for r = (pad_rows+1):(pad_rows+rows)
+        for c = (pad_cols+1):(pad_cols+cols)
+            force = sum(img_pad[r, c] * img_pad[(r-pad_rows):(r+pad_rows), (c-pad_cols):(c+pad_cols)] .* mask)
+            ff_amplitudes[r-pad_rows, c-pad_cols] = abs(force)
+            ff_angles[r-pad_rows, c-pad_cols] = angle(force)
+        end
+    end
+    return ff_amplitudes, ff_angles
+end
+
diff --git a/julia/ndgrid.jl b/julia/ndgrid.jl
new file mode 100644
index 0000000..688a246
--- /dev/null
+++ b/julia/ndgrid.jl
@@ -0,0 +1,52 @@
+# This file is a part of Julia. License is MIT: http://julialang.org/license
+
+ndgrid(v::AbstractVector) = copy(v)
+
+function ndgrid{T}(v1::AbstractVector{T}, v2::AbstractVector{T})
+    m, n = length(v1), length(v2)
+    v1 = reshape(v1, m, 1)
+    v2 = reshape(v2, 1, n)
+    (repmat(v1, 1, n), repmat(v2, m, 1))
+end
+
+function ndgrid_fill(a, v, s, snext)
+    for j = 1:length(a)
+        a[j] = v[div(rem(j-1, snext), s)+1]
+    end
+end
+
+function ndgrid{T}(vs::AbstractVector{T}...)
+    n = length(vs)
+    sz = map(length, vs)
+    out = ntuple(i->Array(T, sz), n)
+    s = 1
+    for i=1:n
+        a = out[i]::Array
+        v = vs[i]
+        snext = s*size(a,i)
+        ndgrid_fill(a, v, s, snext)
+        s = snext
+    end
+    out
+end
+
+meshgrid(v::AbstractVector) = meshgrid(v, v)
+
+function meshgrid{T}(vx::AbstractVector{T}, vy::AbstractVector{T})
+    m, n = length(vy), length(vx)
+    vx = reshape(vx, 1, n)
+    vy = reshape(vy, m, 1)
+    (repmat(vx, m, 1), repmat(vy, 1, n))
+end
+
+function meshgrid{T}(vx::AbstractVector{T}, vy::AbstractVector{T},
+                     vz::AbstractVector{T})
+    m, n, o = length(vy), length(vx), length(vz)
+    vx = reshape(vx, 1, n, 1)
+    vy = reshape(vy, m, 1, 1)
+    vz = reshape(vz, 1, 1, o)
+    om = ones(Int, m)
+    on = ones(Int, n)
+    oo = ones(Int, o)
+    (vx[om, :, oo], vy[:, on, oo], vz[om, on, :])
+end
diff --git a/julia/scharr.jl b/julia/scharr.jl
new file mode 100644
index 0000000..02daeb6
--- /dev/null
+++ b/julia/scharr.jl
@@ -0,0 +1,37 @@
+# -*- coding: utf-8 -*-
+#
+# Calculate the derivatives of an image using the Scharr operator
+# of kernal size 3x3.
+#
+# References:
+# [1] https://en.wikipedia.org/wiki/Sobel_operator
+# [2] http://docs.opencv.org/doc/tutorials/imgproc/imgtrans/sobel_derivatives/sobel_derivatives.html
+#
+# Aaron LI
+# 2015/07/20
+#
+
+# Calculate the derivatives of the image using the Scharr operator
+# img - input image
+# dx  - order of the derivative x
+# dy  - order of the derivative y
+function scharr(img, dx, dy)
+    rows, cols = size(img)
+    img_d = float(img)
+    (isa(dx, Int) && dx >= 0) || error("dx should be an integer >= 0")
+    (isa(dy, Int) && dy >= 0) || error("dy should be an integer >= 0")
+    # Scharr operator
+    Gy = [-3.0 -10.0 -3.0; 0.0 0.0 0.0; 3.0 10.0 3.0];
+    Gx = Gy'
+    # calculate the derivatives using convolution
+    for i = 1:dx
+        img_d = conv2(img_d, Gx)
+    end
+    for i = 1:dy
+        img_d = conv2(img_d, Gy)
+    end
+    # FIXME: 'conv2' will increase the image size
+    rows_d, cols_d = size(img_d)
+    return img_d[(div(rows_d-rows, 2)+1):(div(rows_d-rows, 2)+rows), (div(cols_d-cols, 2)+1):(div(cols_d-cols, 2)+cols)]
+end
+