Multidimensional gradients

[juliohm]

Would you be interested in having a multidimensional gradients implementation?

I finished one for my package, but I think it could be useful for others:

https://github.com/juliohm/ImageQuilting.jl/blob/master/src/ndgradients.jl

It assumes separable kernels and right now it only computes the gradients at specified points of the image. I could add an option to use the imfilter in the entire input image when no points are specified.

The imgradients() available in Images.jl is 2D.

[timholy]

:+1:

[juliohm]

Hi Tim, where should this function be put? Can I create a new file in the source called ndgradients.jl or do you prefer to have it dumped in edge.jl? The code is kind of verbose, I would say it deserves a separate file, but let me know of your preferences. I can submit a PR soon you confirm the location.

[juliohm]

This is how the implementation looks like:

"""
G = ndgradients(img, [points], [method], [border])

Performs edge detection filtering in the N-dimensional array `img`.
Gradients are computed at specified `points` (or indexes) in the
array or everywhere. Available methods: `"sobel"` and `"ando3"`.
Border options:`"replicate"`, `"circular"`, `"reflect"`, `"symmetric"`.

Returns a 2D array `G` with the gradients as rows. The number of rows
is the number of points at which the gradient was computed and the
number of columns is the dimensionality of the array.

Author: Júlio Hoffimann Mendes (juliohm@stanford.edu)
"""
function ndgradients(img::AbstractArray; points::AbstractVector=[],
                     method::AbstractString="ando3", border::AbstractString="replicate")
  extent = size(img)
  ndirs = length(extent)
  npoints = isempty(points) ? prod(extent) : length(points)

  # smoothing weights
  weights = (method == "sobel" ? [1,2,1] :
             method == "ando3" ? [.112737,.274526,.112737] :
             error("Unknown gradient method: $method"))

  # border settings
  border ∉ ["replicate", "circular", "reflect", "symmetric"] && error("border option not supported")

  # pad input image
  imgpad = padarray(img, ones(Int, ndirs), ones(Int, ndirs), border)

  # gradient matrix
  G = zeros(npoints, ndirs)

  for i in 1:ndirs
    # kernel = centered difference + perpendicular smoothing
    if extent[i] > 1
      # centered difference
      idx = ones(Int, ndirs); idx[i] = 3
      kern = reshape([-1,0,1], idx...)
      # perpendicular smoothing
      for j in setdiff(1:ndirs, i)
        if extent[j] > 1
          idx = ones(Int, ndirs); idx[j] = 3
          kern = broadcast(*, kern, reshape(weights, idx...))
        end
      end

      # compute gradient at specified points
      if !isempty(points)
        A = zeros(kern)
        shape = size(kern)
        for (k, p) in enumerate(points)
          icenter = CartesianIndex(ind2sub(extent, p))
          i1 = CartesianIndex(tuple(ones(Int, ndirs)...))
          for ii in CartesianRange(shape)
            A[ii] = imgpad[ii + icenter - i1]
          end

          G[k,i] = sum(kern .* A)
        end
      else
        # simple convolution
        G[:,i] = imfilter(imgpad, kern, "inner")[:]
      end
    end
  end

  G
end

Any suggestions on how to make it more Julian?

[Evizero]

weights = (method == "sobel" ? [1,2,1] :
             method == "ando3" ? [.112737,.274526,.112737] :
             error("Unknown gradient method: $method"))

This does not look type stable. Add a dot to the vector in the first line ([1.,2,1]) to make it so. Given that weights is used within the loop it should have some influence on performance

[juliohm]

Good catch @Evizero , fixed it.

[Evizero]

      # compute gradient at specified points
      if !isempty(points)
        A = zeros(kern)
        shape = size(kern)
        for (k, p) in enumerate(points)
          icenter = CartesianIndex(ind2sub(extent, p))
          i1 = CartesianIndex(tuple(ones(Int, ndirs)...))

does i1 have to be computed in the inner loop? As far as I know the splatting (the ... operation) has some performance penalty, so if you can pre-compute it further up the nested loops you may gain some performance.

[juliohm]

@Evizero I suspect the performance gain is not that big, but I can profile it. It is more natural to have variables defined near their usage. I think the JIT compiler can figure out this optimization, am I wrong?

[Evizero]

I do not know. I would be interested to know if the JIT compiler helps out in such instances.

It is just something that caught my eye, since I enjoy thinking about these kinds of things. I could very well see there not being any empirical performance difference whatsoever.

[timholy]

This is very reasonable to have in Images. For the "everywhere" option, I would split it into a separate function and make the output a d-by-sz array, where d is the dimensionality of img and sz = size(img). For the "specific points" option indeed a 2D array seems reasonable.

We should also just call the function imgradients.

I definitely recommend profiling.

[juliohm]

Sure Tim, I fully agree, will submit a PR with your suggestions soon.

[juliohm]

This is part of the new Images somewhere I think. Closing the issue. :+1: