Piszę sobie mały wykrywacz rozmytych obrazów z użyciem filtra LoG (później pewnie pójdą i inne). No i mam problem z wygenerowaną macierzą. Efekt raczej nie spełnia moich oczekiwań.

Kod:

fn log(x : f32, y : f32, sigma : f32) -> f32 {
    let sigma_sq = sigma * sigma;
    let m = (x * x + y * y) / (2. * sigma_sq);

    -1. / (consts::PI * sigma_sq * sigma_sq) * (1. - m) * consts::E.powf(-m)
}

fn generate_kernel(size : usize) -> Matrix<f32> {
    let val = (size / 2) as f32;
    let mut matrix = Matrix::new(size, size);

    for i in (0..size) {
        for j in (0..size) {
            matrix[(i, j)] = log(-val + (i as f32), -val + (j as f32), 1.4);
        }
    }

    matrix
}

Macierz, która mi wychodzi:

    0.000169    0.000757    0.002068    0.003616    0.004310    0.003616    0.002068    0.000757    0.000169
    0.000757    0.003016    0.006964    0.010024    0.010810    0.010024    0.006964    0.003016    0.000757
    0.002068    0.006964    0.011205    0.006376    0.000610    0.006376    0.011205    0.006964    0.002068
    0.003616    0.010024    0.006376   -0.024365   -0.047824   -0.024365    0.006376    0.010024    0.003616
    0.004310    0.010810    0.000610   -0.047824   -0.082859   -0.047824    0.000610    0.010810    0.004310
    0.003616    0.010024    0.006376   -0.024365   -0.047824   -0.024365    0.006376    0.010024    0.003616
    0.002068    0.006964    0.011205    0.006376    0.000610    0.006376    0.011205    0.006964    0.002068
    0.000757    0.003016    0.006964    0.010024    0.010810    0.010024    0.006964    0.003016    0.000757
    0.000169    0.000757    0.002068    0.003616    0.004310    0.003616    0.002068    0.000757    0.000169

Obraz wyjściowy:

out.png

Obraz oczekiwany (wygenerowany w Octave):

expected.png

Ma ktoś jakieś pomysły? Zakładam, że błąd jest idiotyczny, ale ja już nie mam pomysłu co to mogło by być.

EDIT:

Po głębszych testach wyszło, że błąd jest gdzieś w aplikowaniu filtra:

fn clamp<N: PartialOrd>(val: N, min: N, max: N) -> N {
    if val > max { max }
    else if val < min { min }
    else { val }
}

fn filter<P: Primitive + 'static, T: Pixel<P> + 'static, I: GenericImage<T>>(
    image: &I,
    kernel: &Matrix<f32>) -> ImageBuffer<Vec<P>, P, T> {
        let xcenter = kernel.rows() / 2;
        let ycenter = kernel.cols() / 2;
        let (width, height) = image.dimensions();

        let sum = match kernel.sum() {
            0.0 => 1.0,
            v   => v
        };
        let sumx4 = f32x4(sum, sum, sum, sum);
        let max : f32 = Primitive::max_value();

        let mut out = ImageBuffer::new(width, height);

        for i in (xcenter..width - xcenter) {
            for j in (ycenter..height - ycenter) {
                let mut sum = f32x4(0.0, 0.0, 0.0, 0.0);

                for x in (0..kernel.rows()) {
                    for y in (0..kernel.cols()) {
                        let diffx = x as isize - xcenter as isize;
                        let diffy = y as isize - ycenter as isize;
                        let x0 = i as isize + diffx;
                        let y0 = j as isize + diffy;
                        let k = kernel[(x, y)];
                        let kx4 = f32x4(k, k, k, k);
                        let pixel = image.get_pixel(x0 as u32, y0 as u32);

                        let (c1, c2, c3, c4) = pixel.channels4();
                        let vec = f32x4(
                            cast(c1).unwrap(),
                            cast(c2).unwrap(),
                            cast(c3).unwrap(),
                            cast(c4).unwrap());

                        sum += vec * kx4;
                    }
                }

                let f32x4(s1, s2, s3, s4) = sum / sumx4;
                let t: T = Pixel::from_channels(
                    cast(clamp(s1, 0.0, max)).unwrap(),
                    cast(clamp(s2, 0.0, max)).unwrap(),
                    cast(clamp(s3, 0.0, max)).unwrap(),
                    cast(clamp(s4, 0.0, max)).unwrap());

                out.put_pixel(i, j, t);
            }
        }

        out
    }