imfilter issues

[GunnarFarneback]

I tried to port some Matlab image processing code to Julia, but it turned out that imfilter in extras/image.jl is not in great shape.

1. First test.

julia> load("extras/image.jl")

julia> imfilter(ones(4,4), ones(3,3))
DomainError()
 in power_by_squaring at intfuncs.jl:77
 in imfilter at /home2/gunnar_no_backup/julia/extras/image.jl:686
 in imfilter at /home2/gunnar_no_backup/julia/extras/image.jl:730

This one was easy, replace 10^-7 with 1e-7 on line 686.
2. Try again.

julia> imfilter(ones(4,4), ones(3,3))
4x4 Float64 Array:
 9.0  9.0  9.0  9.0
 9.0  9.0  9.0  9.0
 9.0  9.0  9.0  9.0
 9.0  9.0  9.0  9.0

Correct.
3. Odd sized image.

julia> imfilter(ones(3,3), ones(3,3))
argument dimensions must match
 in promote_shape at operators.jl:168
 in .* at array.jl:927
 in imfilter at /home2/gunnar_no_backup/julia/extras/image.jl:700
 in imfilter at /home2/gunnar_no_backup/julia/extras/image.jl:730

On lines 697-699 something is bogus:

B[int((length(x))/2)+1:sa[1]+int((length(x))/2),int((length(y))/2)+1:sa[2]+int((length(y))/2)] = A
y = fft([zeros(T,int((m-length(y)-1)/2)); y; zeros(T,int((m-length(y)-1)/2))])
x = fft([zeros(T,int((m-length(x)-1)/2)); x; zeros(T,int((n-length(x)-1)/2))])

The lack of symmetry indicates that one m that should be an n, probably in the construction of x. But changing that is not enough, there's also some kind of mismatch by one in the sizes of B relative to x and y.
4. Try again but avoid the separable code path.

julia> imfilter(ones(3,3), [1 1 1;1 0.0 1;1 1 1])
3x3 Float64 Array:
 3.0  5.0  5.0
 5.0  8.0  8.0
 5.0  8.0  8.0

Wrong result, should be all eights. On line 717:

out = C[int(sc[1]/2-si[1]/2):int(sc[1]/2+si[1]/2)-1, int(sc[2]/2-si[2]/2):int(sc[2]/2+si[2]/2)-1]

Something is off by one in the odd size case.
5. In the border replicate code the corners of the padded area are filled with mirrored image values (lines 651-654). Where does that come from? It's not how Matlab does it and neither have I seen it elsewhere before.

[timholy]

Agreed wholeheartedly that it's not in great shape. I got started on a big refactoring early on, before realizing that (1) I didn't know enough Julia, and (2) that Julia didn't yet have enough tricks up her sleeve. I think we're there now, though, so image.jl could be rebuilt. In fact I'm doing a huge amount of image processing in Julia right now, it's just not in image.jl due to me not having enough time to release it (grants, ugh).

I think that in this case, however, this particular error predated my efforts. To me it sounds like you're running into boundary condition issues. Reflected boundary conditions are, in fact, one of the standard approaches, but I agree it is not what Matlab does by default.

I have a whole suite of tools for specifying different boundary conditions (fill with zeros, reflecting, periodic, NaN, etc.) for operations like interpolation, and those could be used for specifying what you want. Unfortunately I really don't have time to release this code yet, because every time I do that there's a whole explosion of issues that crop up, and dealing with those eats up a huge amount of time. So it will be mid-November or so. I have to encourage you to figure this one out on your own (sorry).

[pao]

@timholy agree with the tags? It sounds like there are definitely bugs here, but you won't be able to work on them for a while. (I'm about to add the same to the ODE bug since any issues have exceeded my depth.)

[timholy]

Yep.

[GunnarFarneback]

To me it sounds like you're running into boundary condition issues.

The problem is not with the boundary conditions as such but the padding and unpadding is not quite right.

Reflected boundary conditions are, in fact, one of the standard approaches, but I agree it is not what Matlab does by default.

But mixed reflection in the corners and replication on the sides is not so standard, is it?

I have to encourage you to figure this one out on your own (sorry).

No thanks, fft based image processing is not my cup of tea. I'll leave this to someone who is more familiar with the padding intricacies. If there's any interest in a spatial implementation of imfilter I might give it a shot though, although it probably would need to be a mixed Julia/C approach.

[StefanKarpinski]

I think the solution here is just to wait until Tim gets his grant proposals all done. Waiting until mid-November at this point is not bad. Good luck with the grants, Tim!

[timholy]

For small filters, spatial is often much faster. So I think that would be quite interesting, and I'd be very happy if you tackled it.

Out of curiosity, why would it need to be mixed Julia/C? Issue #1392? Aside from that, Julia's performance in such very low-level operations is identical to C (read that very long thread that I linked to in that issue).

[GunnarFarneback]

I was only thinking in terms of performance along the lines of half the speed of C. For a workhorse like imfilter that's often not good enough. But if a pure Julia implementation can give acceptable speed that's of course a more attractive solution. And yes, bounds-checking inside tight imfilter loops is not what you want to spend cycles on.

[timholy]

In 1d, the penalty from bounds checking in a tight loop seems to be about 3:4 (C execution time:Julia execution time). In higher dimensions it will get worse, of course. However, if you can do all of your indexing in terms of scalar offsets, then it should stay constant.

If you want to tackle this for multidimensions, then abstractarray.jl's make_arrayind_loop_nest might be a good way to go: it will autogenerate functions for you in arbitrary dimensions. (Functions in array.jl call this via gen_array_index_map, and hence might serve as examples of usage.) If you're not up for that and just want to handle the 2d case, then at least the comments above make_arrayind_loop_nest describe a very efficient mechanism for generating a scalar offset. There are also explicit 2d examples in ref and assign in array.jl for Matrix objects.

[JeffBezanson]

I have high hopes that llvm optimizations will be able to hoist redundant bounds checks and offset computations if I can make the right code gen improvements. I can also easily add a switch to disable bounds checks.

[StefanKarpinski]

In general, if pure Julia code is not fast enough for something, our approach is to work on optimizations until it is fast enough, not to code it in C. That's how we've gotten where we are :-)

[ViralBShah]

Could we do bounds check disabling through a macro selectively, rather than with a global switch?

@noboundscheck for ....

end

[timholy]

Last I knew, that was that plan.

@GunnarFarneback , if Jeff disables bounds-checking (or if the cost goes to zero through other optimizations), then for tight loops our previous measurements suggest that C and Julia get identical performance. So I'd definitely urge you to aim for a pure-Julia solution.

[JeffBezanson]

BTW @GunnarFarneback thanks for the useful report.

[JeffBezanson]

I think I fixed some of the off-by-one errors, but I'm not sure what to do with the border replicate case.

[GunnarFarneback]

The replicate and mirror padding code is still broken for non-quadratic filters, e.g.

julia> imfilter(ones(4,4),ones(1,3),"replicate")
argument dimensions must match
 in assign at array.jl:531
 in imfilter at /home2/gunnar_no_backup/julia/extras/image.jl:657
 in imfilter at /home2/gunnar_no_backup/julia/extras/image.jl:747

Unless someone needs the exact padding behaviour of the current code I'd suggest dropping it and calling this Matlab/Octave compatible padarray:

function padarray{T,n}(img::Array{T,n}, prepad::Vector{Int}, postpad::Vector{Int}, border::String, value::T)
    I = Array(Vector{Int}, n)
    for d = 1:n
        M = size(img, d)
        I[d] = [(1 - prepad[d]):(M + postpad[d])]
        if border == "value"
            I[d] = int((I[d] .>= 1) & (I[d] .<= M))
        elseif border == "replicate"
            I[d] = min(max(I[d], 1), M)
        elseif border == "circular"
            I[d] = 1 + mod(I[d] - 1, M)
        elseif border == "symmetric"
            I[d] = [1:M, M:-1:1][1 + mod(I[d] - 1, 2 * M)]
        elseif border == "reflect"
            I[d] = [1:M, M-1:-1:2][1 + mod(I[d] - 1, 2 * M - 2)]
        else
            error("unknown border condition")
        end
    end

    if border == "value"
        A = Array(T, map(length, I)...)
        fill!(A, value)
        A[map(x->map(bool, x), I)...] = img
    else
        A = img[I...]
    end

    return A
end

For more convenience and better Matlab/Octave call compatibility:

padarray{T,n}(img::Array{T,n}, padding::Vector{Int}, border::String) = padarray(img, padding, padding, border, zero(T))

padarray{T,n}(img::Array{T,n}, padding::Vector{Int}, value::T) = padarray(img, padding, padding, "value", value)

function padarray{T,n}(img::Array{T,n}, padding::Vector{Int}, border::String, direction::String)
    if direction == "both"
        return padarray(img, padding, padding, border, zero(T))
    elseif direction == "pre"
        return padarray(img, padding, 0 * padding, border, zero(T))
    elseif direction == "post"
        return padarray(img, 0 * padding, padding, border, zero(T))
    end
end

function padarray{T,n}(img::Array{T,n}, padding::Vector{Int}, value::T, direction::String)
    if direction == "both"
        return padarray(img, padding, padding, "value", value)
    elseif direction == "pre"
        return padarray(img, padding, 0 * padding, "value", value)
    elseif direction == "post"
        return padarray(img, 0 * padding, padding, "value", value)
    end
end

[JeffBezanson]

Thanks, looks great. Please pull request it and I will merge it.

[vtjnash]

@JeffBezanson is it generally recommend to use strings, symbols, or enum values for options like this?

[JeffBezanson]

Symbols are the most efficient, but I haven't been fussy about it.

[StefanKarpinski]

Symbols do seem like the best choice and they're also the shortest. A speculative proposal would be having key: sym as a shorthand for key=:sym in function arguments, like so:

foo(arg1, arg2, key: sym)

[JeffBezanson]

Probably not, since that would make x:y and x: y utterly different.

[StefanKarpinski]

Oh, right. Poor overloaded : — so many meanings! The space was only for clarity, not to imply that you couldn't write x:y, which we obviously can't change the pre-existing meaning of. So forget that idea.

[GunnarFarneback]

I have made some patches for imfilter.

1. GunnarFarneback@62e58c6 is a bugfix for a broken 180 degree rotation of the filter kernel.
2. GunnarFarneback@4711b93 adds the padarray implementation shown some comments above and replaces the imfilter border handling with a call to padarray.
3. GunnarFarneback@51d8c86 implements a spatial filter for arbitrary image dimensionality, temporarily called imfilterN.

Is a pull request still desired in these times of packages or should it wait for someone to packagize extras/image.jl?

[pao]

I'd pull req at least GunnarFarneback@62e58c6 since it's a bugfix. I've been doing the same sorts of fixes for extras/strpack.jl, even though I'm working towards a package which has turned into a major overhaul.

[timholy]

While I haven't packaged it officially yet, there's https://github.com/timholy/Image.jl. More interesting than the README is the code in src/; it's in a half-baked state, but to a capable programmer such as yourself it may be sufficiently fleshed out for you to see where it's going. I'd be interested in your views, and/or happy with you taking things in other directions.

[JeffBezanson]

@GunnarFarneback could you make a pull request with your patches against Tim's package, https://github.com/timholy/Images.jl ?