median not correct for even/odd arrays size

[PrometheusPi]

The current median implementation:

...
half = int(n / 2)
median = (sort_array[half] + sort_array[-half]) / 2.0
output_data[depth_global, col_global, row_global] = 0
...

will not produce the correct result. As an example, please have a look at this code:

myN = 6
myTest = np.arange(myN)
half = int(myN/2)

print("Length =", myN)
print("Numpy_median =", np.median(myTest))
print("Array =", myTest)
print("selection_left:", myTest[half] , " selecton right:", myTest[-half])
print("result =", (myTest[half] + myTest[-half])/2)

For an odd length array - the result is the center value, but your method selects to different array values:

Length = 5
Numpy_median = 2.0
Array = [0 1 2 3 4]
selection_left: 2  selecton right: 3
result = 2.5

For an even length array, you should compute the mean between the two array values around the center - but you select the larger value twice:

Length = 6
Numpy_median = 2.5
Array = [0 1 2 3 4 5]
selection_left: 3  selecton right: 3
result = 3.0

The misconception might arris from the fact that regular counting is c-like, but rrevese indexing starts with 1.

This can be precented by changing the negative indexed array value to ...[-half-1].

myN = 5
myTest = np.arange(myN)
half = int(myN/2)

print("Length =", myN)
print("Numpy_median =", np.median(myTest))
print("Array =", myTest)
print("selection_left:", myTest[half] , " selecton right:", myTest[-half-1])
print("result =", (myTest[half] + myTest[-half-1])/2)

This produces the correct result for both odd and even sized arrays.

[PrometheusPi]

Furthermore, the value is not yet added to the output array. (Currently it is set to zero.)

[chtietz]

Unfortunately we can't use that implementation. The sort_array is a local array and it's size has to be a simple constant expression (doc). In that case negative indexing is not purposeful. Pls review my suggestion: https://github.com/ComputationalRadiationPhysics/stencil_filter_on_GPU/commit/4633214ae2de4b2f82acd2a6f317f8e27bd8e2a0

[PrometheusPi]

I think your current implementation does not run into this error. You could close this issue, it you want.