Error in match_histogram_mut()

[JakobZahn]

The function match_histogram_mut does not compute the optimal image for the following example:

let mut input_image = gray_image!(type: u8, 0, 1, 2, 3, 4, 5, 6, 7);
let target_image = gray_image!(type: u8, 30, 130, 30, 130);

match_histogram_mut(&mut input_image, &target_image);

for p in input_image.iter_mut() {
    println!("{p}");
}

Output:

29
30
30
30
129
130
130
130

However, the two histograms of the images would match better, if the 29 were a 30 and the 129 were a 130.

The function match_histogram() is based on the function histogram_lut():

pub fn match_histogram_mut(image: &mut GrayImage, target: &GrayImage) {
    let image_histc = cumulative_histogram(image).channels[0];
    let target_histc = cumulative_histogram(target).channels[0];
    let lut = histogram_lut(&image_histc, &target_histc);

    for p in image.iter_mut() {
        *p = lut[*p as usize] as u8;
    }
}

I.e. for each possible old gray value the array lut (result of histogram_lut()) stores an updated gray value.

And histogram_lut() has a unit test that expects this kind of suboptimal result:

    [test]
    fn test_histogram_lut_gradient_to_step_contrast() {
        let mut grad_histc = [0u32; 256];
        for i in 0..grad_histc.len() {
            grad_histc[i] = i as u32;
        }

        let mut step_histc = [0u32; 256];
        for i in 30..130 {
            step_histc[i] = 100;
        }
        for i in 130..256 {
            step_histc[i] = 200;
        }

        let lut = histogram_lut(&grad_histc, &step_histc);
        let mut expected = [0usize; 256];

        // No black pixels in either image
        expected[0] = 0;

        for i in 1..64 {
            expected[i] = 29;
        }
        for i in 64..128 {
            expected[i] = 30;
        }
        for i in 128..192 {
            expected[i] = 129;
        }
        for i in 192..256 {
            expected[i] = 130;
        }

        assert_eq!(&lut[0..256], &expected[0..256]);
    }

The comment above the definition of histogram_lut() states:
l = histogram_lut(s, t) is chosen so that target_histc[l[i]] / sum(target_histc) is as close as possible to source_histc[i] / sum(source_histc)

I interpret this as:
Find the function l so that for every entry i in the histogram of the source image the ratio target_histc[l[i]] / sum(target_histc) is most similar to source_histc[i] / sum(source_histc)
where target_histc and source_histc are the respective cumulative histograms.

In my little example at the top this function is chosen so that
l(1) == 29
l(2) == 30
l(3) == 30
l(4) == 30
l(5) == 129
l(6) == 130
l(7) == 130
l(8) == 130

target_histc should be an array with 256 elements: 30 times a 0, then 100 times a 2 and the remaining entries a 4
source_histc should be an array with 256 elements: [0, 1, 2, 3, 4, 5, 6, 7, 8, 8, 8, 8, ...]

But if I insert the values i==1 to i==8 in the equation I do not get reasonable results, so I might be wrong somewhere or the comment is misleading.

If this is indeed a bug, I could look into possibilities to fix it.