@pavanky I wish I saw the PR before it got merged...
If I'm not mistaken with the merged implementation, at least for the CPU version, the missing value case (i.e. when none of the masks are set) the result will be set to numeric_limits<T>::max() instead of numeric_limits<T>::infinity(). For example for float the former is 3.40282e+38, while the latter is inf. This might be problematic for a few reasons

1. since the value will be a real number, subsequent operations on it will not trigger any special inf-arithmetic, so it will be very hard to detect that something went wrong after morph() was called.
2. I understand that ArrayFire doesn't need to replicate Matlab behaviour, but I think the fact that the result will be different than inf, which is what Matlab does, will cause more confusion, especially because it's not really obvious why 3.40282e+38 was returned.
3. in many situations, after dilate() is called, people do something like this

out = dilate(in, filter);
out(isInf(out)) = 0;

If we use max instead of inf, this type of code will be messier to write.

So, I was thinking, ideally we should let the user pass an optional missingValue argument to these functions. By default, if the missingValue is not given, it should use the special initial value (but still with inf for floats). I'm not sure if c++ allows passing template dependent default arguments to a function. If that's not possible, we might need to split these functions into two:

template<typename T, bool IsDilation>
void morph(Array<T> out, Array<T> const in, Array<T> const mask, T missingValue)
{
}

template<typename T, bool IsDilation>
void morph(Array<T> out, Array<T> const in, Array<T> const mask)
{
    T missingValue = ... // -inf/inf for floats, min/max for integrals
    morph(out, int, mask, missingValue);
}

Note that, for the first function, since the passed missingValue still cannot be used as the initial value, we would have to keep a boolean when we do the inner loop:

T filterResult = init;
bool missing = true;
for (wj) {
    for (wi) {
        if (...) {
            missing = false;
            // ...
        }
    }
}
outData[ getIdx(ostrides, i, j) ] = missing ? missingValue : filterResult;

The nice thing about having this default missingValue, it allows the user not to write the messy and expensive code as in 3 above, instead it just becomes

out = dilate(in, filter, 0); // no need to worry about replacing inf or max anymore

BTW, these "missing values" are not really some strange corner case. They can happen very often near the boundaries of the input array, because near the boundaries, part of the filter falls outside of the array, and the rest could easily have all zeros.

This is an easy fix. I can change the init value to be +/- infinity for float types in our Binary<T, af_max_t> and Binary<T, af_min_t> functors.

Passing an additional parameter now will either break the API or will need implementing a new function.

Also changing the value from max to infinity does not break other functions where we are using this functor. This will be the proper and consistent fix.

@vakopian Thanks for pointing out the flaw btw. If you agree that using inifnity is a fair solution, I'll go ahead and do that.

@pavanky of course, let's start by changing it to infinity. We can perhaps add the additional parameter/extra methods in a later release. BTW, is there a chance this fix will be in a 3.3.3 release?

@vakopian There isn't going to be a 3.3.3. The next release will be 3.4 and it should be out very soon.

@pavanky thanks for the information. While we're on the subject of morph, I'm seeing very strange performance results for this function. The following code takes less than 1 second on matlab, it takes about 4 seconds in GNU octave, but takes about 110 seconds with ArrayFire (timing the dilate function only). All are CPU only versions:

in = randn(3000, 3000);
filter = constant(1, 20, 20);
dilate(in, filter);

I've looked at the CPU implementation of morph() in ArrayFire, and I can't see anything obviously inefficient about it. I've played with some trivial optimizations there, such as moving the evaluation of the maskValue inside the if condition, or defining the loop boundaries outside of the loops to avoid their repeated evaluations, but they did not result in any significant improvement. I was wandering if you might have an idea how Matlab and octave achieve so much better performance.

@vakopian The CPU backend is not at all optimized. It runs on a single core and is not vectorized. It exists as a fallback solution. This will change in the future. For now, if you want a better performing CPU version try out with the intel opencl driver installed.

@pavanky: sure, I understand that we can get better performance with vectorization and such, but I was thinking there might be purely algorithmic ways to improve this. After some digging around, I think we could achieve a big performance improvement by first pre-calculating the linear offsets for the given neighborhood (i.e. the non-zeros of the filter), then use that information directly, rather than testing the entire neighborhood for every pixel.
I'm thinking to try to implement this idea for ArrayFire (at least CPU version). But first I'd like to know if you would be interested in such a contribution, even if it's CPU only, and also if there is already code in the library that does such neighborhood-to-offsets translations, so I don't have to reinvent it.

@vakopian We do not have any translations like that as of now. If you can send in a PR for such a contribution, we could consider generalizing it for other window operations.

Can you perhaps create a new issue to discuss this further?