pub fn gradients<P, F, Q>(
    image: &Image<P>,
    horizontal_kernel: Kernel<i32>,
    vertical_kernel: Kernel<i32>,
    f: F,
) -> Image<Q>
where
    P: Pixel<Subpixel = u8> + WithChannel<u16> + WithChannel<i16>,
    Q: Pixel,
    ChannelMap<P, u16>: HasBlack,
    F: Fn(ChannelMap<P, u16>) -> Q,
{
    let horizontal = filter_clamped::<_, _, i16>(image, horizontal_kernel);
    let vertical = filter_clamped::<_, _, i16>(image, vertical_kernel);

    let (width, height) = image.dimensions();
    let mut out = Image::<Q>::new(width, height);

    // This would be more concise using itertools::multizip over image pixels, but that increased runtime by around 20%
    for y in 0..height {
        for x in 0..width {
            // JUSTIFICATION
            //  Benefit
            //      Using checked indexing here makes this gradients 1.1x slower,
            //      as measured by bench_sobel_gradients
            //  Correctness
            //      x and y are in bounds for image by construction,
            //      vertical and horizontal are the result of calling filter_clamped on image,
            //      and filter_clamped returns an image of the same size as its input
            let (h, v) = unsafe {
                (
                    horizontal.unsafe_get_pixel(x, y),
                    vertical.unsafe_get_pixel(x, y),
                )
            };
            let mut p = ChannelMap::<P, u16>::black();

            for (h, v, p) in multizip((h.channels(), v.channels(), p.channels_mut())) {
                *p = gradient_magnitude(*h as f32, *v as f32);
            }

            // JUSTIFICATION
            //  Benefit
            //      Using checked indexing here makes this sobel_gradients 1.1x slower,
            //      as measured by bench_sobel_gradients
            //  Correctness
            //      x and y are in bounds for image by construction,
            //      and out has the same dimensions
            unsafe {
                out.unsafe_put_pixel(x, y, f(p));
            }
        }
    }

    out
}
Computes per-channel gradients using the given horizontal and vertical kernels and calls `f`
to compute each output pixel.

# Examples
```
# extern crate image;
# #[macro_use]
# extern crate imageproc;
# fn main() {
use imageproc::gradients::gradients;
use imageproc::kernel::Kernel;
use imageproc::kernel;
use image::Luma;
use std::cmp;

// A shallow horizontal gradient in the red
// channel, a steeper vertical gradient in the
// blue channel, constant in the green channel.
let input = rgb_image!(
[0, 5, 0], [1, 5, 0], [2, 5, 0];
[0, 5, 2], [1, 5, 2], [2, 5, 2];
[0, 5, 4], [1, 5, 4], [2, 5, 4]
);

// Computing independent per-channel gradients.
let channel_gradient = rgb_image!(type: u16,
[ 4,  0,  8], [ 8,  0,  8], [ 4,  0,  8];
[ 4,  0, 16], [ 8,  0, 16], [ 4,  0, 16];
[ 4,  0,  8], [ 8,  0,  8], [ 4,  0,  8]
);

let horizontal_kernel = kernel::SOBEL_HORIZONTAL_3X3;
let vertical_kernel = kernel::SOBEL_VERTICAL_3X3;

assert_pixels_eq!(
gradients(&input, horizontal_kernel, vertical_kernel, |p| p),
channel_gradient
);

// Defining the gradient of an RGB image to be the
// mean of its per-channel gradients.
let mean_gradient = gray_image!(type: u16,
4, 5, 4;
6, 8, 6;
4, 5, 4
);

assert_pixels_eq!(
gradients(&input, horizontal_kernel, vertical_kernel, |p| {
let mean = (p[0] + p[1] + p[2]) / 3;
Luma([mean])
}),
mean_gradient
);

// Defining the gradient of an RGB image to be the pixelwise
// maximum of its per-channel gradients.
let max_gradient = gray_image!(type: u16,
8,  8,  8;
16, 16, 16;
8,  8,  8
);

assert_pixels_eq!(
gradients(&input, horizontal_kernel, vertical_kernel, |p| {
let max = cmp::max(cmp::max(p[0], p[1]), p[2]);
Luma([max])
}),
max_gradient
);
# }
src/gradients.rs