Composition bug under ICE and Hugin

ICE and Hugin suffer from a compositional problem when they have to deal with columns composed of many images.

If the light is not perfectly homogeneous (which is often the case when using a microscope) and forms a slight gradient, then this gradient will be amplified on the final composition, which will give a disastrous result.

6502
ICE: Gradient correction error
© Boris Marmontel


To fix this problem, you have to make a manual gradient correction on each image.

To do this, we take a “plain” photograph of a white background, which will serve as a reference of the gradient.

6502


Now the idea is to code a program that will subtract this gradient.

Image im2("light_reference.png");
Image im1("input_image.png");

std::array<Layer64f, 3> data1;
std::array<Layer64f, 3> data2;
std::array<Channel, 3> channels = {
    Channel::Red,
    Channel::Green,
    Channel::Blue
};

for(unsigned i=0; i<3; ++i)
{
    Layer64f& l1 = data1[i];
    Layer64f& l2 = data2[i];
    
    l1 = extractChannel(im1, channels[i]);
    l2 = extractChannel(im2, channels[i]);
    
    double min1 = *std::min_element(l1.data.begin(), l1.data.end());
    double max1 = *std::max_element(l1.data.begin(), l1.data.end());
    double delta1 = max1 - min1;
    
    auto it1 = l1.data.begin();
    auto it2 = l2.data.begin();
    
    double factor = 0.3;
    while(it1 != l1.data.end() && it2 != l2.data.end())
    {
        *it1 -= *it2 * factor;
        
        ++it1;
        ++it2;
    }
    
    // normalize
    double min = *std::min_element(l1.data.begin(), l1.data.end());
    double max = *std::max_element(l1.data.begin(), l1.data.end());
    double delta = max - min;
    
    // [0-1]
    for(double& x: l1.data)
    {
        x = (x - min) / delta;
    }
    
    // restaure previous histogram state
    // [min1-max1]
    for(double& x: l1.data)
    {
        x = min1 + delta1*x;
    }
}


Layer64f is just an image represented with double-precision floating point numbers.

struct Layer64f
{
    std::vector<double> data;
    od::Sizeu size;
    
    Layer64f() {}
    Layer64f(od::Sizeu size):
        data(size.wid * size.hei), size(size) {}
    
    void set(od::Posu p, double c) { data[p.x + p.y * size.wid] = c; }
    double get(od::Posu p) const { return data[p.x + p.y * size.wid]; }
};


Once the layers have been calculated, the image is simply recomposed.

Here is a brightness correction with factor = 0.6 then factor = 0.3

6502
Left: input image, right: light corrected with factor = 0.6


6502
Left: input image, right: light corrected with factor = 0.3


As you can see, with 0.6 the correction is too strong and the panorama software will have a similar bug but this time with the too light part at the bottom. With the correction in 0.3 the panorama will be perfectly uniform.

6502     6502
0.6             0.3


The ideal value of the correction coefficient depends on the height in pixels of the column. Before making a complete panorama, it is therefore advisable to try different values on a single column to check that everything is fine.