jcelaya / hdrmerge

HDR exposure merging
http://jcelaya.github.io/hdrmerge/
Other
371 stars 78 forks source link

Wrong computed `satThreshold` #126

Closed Floessie closed 6 years ago

Floessie commented 6 years ago

Hi,

with my new camera (Canon 6D) I sometimes have an issue with HDRMerge: The mask is solid green (brightest image). This seems to be due to the value fed to Image::isSaturated() via Image::isSaturatedAround() either being too low or the computed satThreshold too high. The max value fed to Image::isSaturated() was 13314, the computed satThreshold 13556.

If I set the custom WB to 13314, everything is fine, but I can only get this value by adding printf() code to Image::isSaturated(), which is impractical at least. Why do those values differ?

Here's the output of -vv with the additional patch below:

Using LibRaw 0.17.2-Release
Number of frames : Number of frames : 11

Number of frames : 1
20180320-183828-1.cr2: 5496x3670 (5568x3708+72+38, by Canon EOS 6D, 400ISO 1/50sec f3.2 EV:-7
b4b4b4b4 RGBG, sat 15490, black 2047, flip 5, wb: 2305 1024 1458 1024, cblack: 2047 2047 2048 2047
20180320-183828.cr2: 5496x3670 (5568x3708+72+38, by Canon EOS 6D, 400ISO 1/3200sec f3.2 EV:-13
b4b4b4b4 RGBG, sat 15490, black 2047, flip 5, wb: 2305 1024 1458 1024, cblack: 2047 2047 2048 2048
20180320-183826.cr2: 5496x3670 (5568x3708+72+38, by Canon EOS 6D, 400ISO 1/400sec f3.2 EV:-10
b4b4b4b4 RGBG, sat 15490, black 2047, flip 5, wb: 2305 1024 1458 1024, cblack: 2047 2047 2048 2047
Load files: 0.710646 seconds
13693 14336 13859 14335 
params.max = 15490, satThreshold = 15490
A: satThreshold = 13693
B: satThreshold = 13556
Compute response functions: 0.187964 seconds
Generate mask: 0.0791689 seconds
diff --git a/src/ImageStack.cpp b/src/ImageStack.cpp
index 9097bb4..ab918e5 100644
--- a/src/ImageStack.cpp
+++ b/src/ImageStack.cpp
@@ -87,17 +87,26 @@ void ImageStack::calculateSaturationLevel(const RawParameters & params, bool use
         }
     }

+    for (auto m : maxPerColor) {
+        printf("%u ", m);
+    }
+    printf("\n");
+
     satThreshold = params.max == 0 ? maxPerColor[0] : params.max;
+    printf("params.max = %u, satThreshold = %u\n", params.max, satThreshold);
     for (int c = 0; c < 4; ++c) {
         if (maxPerColor[c] > 0 && maxPerColor[c] < satThreshold) {
             satThreshold = maxPerColor[c];
         }
     }
+    printf("A: satThreshold = %u\n", satThreshold);
     if(!useCustomWl) // only scale when no custom white level was specified
         satThreshold *= 0.99;
     else
         Log::debug( "Using custom white level ", params.max );

+    printf("B: satThreshold = %u\n", satThreshold);
+
     for (auto & i : images) {
         i.setSaturationThreshold(satThreshold);
     }

Best, Flössie

heckflosse commented 6 years ago

@Floessie Does 13443 (15490 - 2047) also work?

Floessie commented 6 years ago

@heckflosse No. Why should it?

You can find the files here. Those are crap handheld shots, but expose the problem well.

Thanks for looking into it.

Best, Flössie

Floessie commented 6 years ago

Ingo, I've come to the conclusion, that ImageStack::calculateSaturationLevel() computes the right value. There seem to be some outliers (single digit number in the histogram) on every channel (stuck pixels?) which are taken into account, while the majority clips at 13314 or 13313. The latter you can also see in "No demosaic" mode in RT. Is there something wrong with my camera?

Here's a patch that looks for the highest index with more than one per mill pixels in the histogram instead of reducing the computed satThreshold to 99%. Maybe this is rubbish, and surely the problem lies somewhere else, but it solves the problem with this stack.

diff --git a/src/Image.hpp b/src/Image.hpp
index dc0fded..972fdeb 100644
--- a/src/Image.hpp
+++ b/src/Image.hpp
@@ -74,6 +74,10 @@ public:
         return brightness > r.brightness;
     }
     void setSaturationThreshold(uint16_t sat);
+    uint16_t getMax() const
+    {
+        return max;
+    }

 private:
     struct ResponseFunction {
diff --git a/src/ImageStack.cpp b/src/ImageStack.cpp
index 9097bb4..b271236 100644
--- a/src/ImageStack.cpp
+++ b/src/ImageStack.cpp
@@ -50,39 +50,38 @@ int ImageStack::addImage(Image && i) {

 void ImageStack::calculateSaturationLevel(const RawParameters & params, bool useCustomWl) {
     // Calculate max value of brightest image and assume it is saturated
-    uint16_t maxPerColor[4] = { 0, 0, 0, 0 };
     Image & brightest = images.front();
-    #pragma omp parallel
-    {
-        uint16_t maxPerColorThr[4] = { 0, 0, 0, 0 };
-        #pragma omp for schedule(dynamic,16) nowait
-        for (size_t y = 0; y < height; ++y) {
-            // get the color codes from x = 0 to 5, works for bayer and xtrans
-            uint16_t fcrow[6];
-            for(size_t i = 0; i < 6; ++i) {
-                fcrow[i] = params.FC(i, y);
-            }
-            size_t x = 0;
-            for (; x < width - 5; x+=6) {
-                for(size_t j = 0; j < 6; ++j) {
-                    uint16_t v = brightest(x + j, y);
-                    if (v > maxPerColorThr[fcrow[j]]) {
-                        maxPerColorThr[fcrow[j]] = v;
-                    }
-                }
+
+    std::vector<std::vector<size_t>> histograms(4, std::vector<size_t>(brightest.getMax() + 1));
+
+    for (size_t y = 0; y < height; ++y) {
+        // get the color codes from x = 0 to 5, works for bayer and xtrans
+        uint16_t fcrow[6];
+        for(size_t i = 0; i < 6; ++i) {
+            fcrow[i] = params.FC(i, y);
+        }
+        size_t x = 0;
+        for (; x < width - 5; x+=6) {
+            for(size_t j = 0; j < 6; ++j) {
+                uint16_t v = brightest(x + j, y);
+                ++histograms[fcrow[j]][v];
             }
-            // remaining pixels
-            for (size_t j = 0; x < width; ++x, ++j) {
-                 uint16_t v = brightest(x, y);
-                if (v > maxPerColorThr[fcrow[j]]) {
-                    maxPerColorThr[fcrow[j]] = v;
-                 }
-             }
         }
-        #pragma omp critical
-        {
-            for(int c = 0; c < 4; ++c) {
-                maxPerColor[c] = std::max(maxPerColorThr[c], maxPerColor[c]);
+        // remaining pixels
+        for (size_t j = 0; x < width; ++x, ++j) {
+            uint16_t v = brightest(x, y);
+            ++histograms[fcrow[j]][v];
+        }
+    }
+
+    uint16_t maxPerColor[4] = { 0, 0, 0, 0 };
+
+    for (int c = 0; c < 4; ++c) {
+        for (int i = histograms[c].size() - 1; i >= 0; --i) {
+            const size_t v = histograms[c][i];
+            if (v > width * height / 1000) {
+                maxPerColor[c] = i;
+                break;
             }
         }
     }
@@ -93,9 +92,8 @@ void ImageStack::calculateSaturationLevel(const RawParameters & params, bool use
             satThreshold = maxPerColor[c];
         }
     }
-    if(!useCustomWl) // only scale when no custom white level was specified
-        satThreshold *= 0.99;
-    else
+
+    if(useCustomWl)
         Log::debug( "Using custom white level ", params.max );

     for (auto & i : images) {

Sorry, I killed your OMP on the way...

Best, Flössie

heckflosse commented 6 years ago

@Floessie According to your measures the whitelevel doesn't fit to the white level in rt camconst.json https://github.com/Beep6581/RawTherapee/blob/dev/rtengine/camconst.json#L517

For ISO 400 it should be 15180 while the one you measured is 13314 which matches the wl (13100) in camconst.json reported for 160, 320, 640, 1250, 2500 much better.

??? Maybe ping Ilias?

Floessie commented 6 years ago

@heckflosse I'm still puzzled as to why this relates to camconst.json. Does HDRMerge use that data, too? When using Amaze in RT and the Neutral profile, the sky is perfectly clipped (RGB 255). Without demosaicing it shows 13314 for one of RGB.

@iliasg I'm lost here. I suspect, that some setting on my new camera is wrong, but it could as well be something more serious. Could you give me a clue what's wrong with the 20180320-183828-1.cr2 in this filebin? The vast majority of the clipped values are at 13313/13314, but when creating a histogram from the data provided by libraw in HDRMerge there are a few points with values above that. Here's a table directly generated within the histogram patch given above:

Color 0: [13693] = 1
Color 0: [13398] = 1
Color 0: [13375] = 1
Color 0: [13362] = 123
Color 0: [13361] = 1
Color 0: [13347] = 1
Color 0: [13332] = 1
Color 0: [13330] = 1
Color 0: [13316] = 1
Color 0: [13314] = 1560014
Color 1: [14336] = 3
Color 1: [14122] = 1
Color 1: [14040] = 2
Color 1: [13865] = 1
Color 1: [13859] = 1
Color 1: [13852] = 1
Color 1: [13831] = 1
Color 1: [13756] = 1
Color 1: [13735] = 1
Color 1: [13730] = 1
Color 1: [13715] = 1
Color 1: [13704] = 1
Color 1: [13586] = 1
Color 1: [13572] = 1
Color 1: [13560] = 1
Color 1: [13518] = 1
Color 1: [13424] = 1
Color 1: [13375] = 1
Color 1: [13362] = 131
Color 1: [13314] = 1718826
Color 2: [13859] = 1
Color 2: [13783] = 1
Color 2: [13466] = 1
Color 2: [13440] = 1
Color 2: [13385] = 2
Color 2: [13362] = 1
Color 2: [13361] = 116
Color 2: [13348] = 1
Color 2: [13341] = 1
Color 2: [13324] = 1
Color 2: [13319] = 1
Color 2: [13313] = 1646195
Color 3: [14335] = 8
Color 3: [14320] = 1
Color 3: [14306] = 2
Color 3: [14182] = 1
Color 3: [14162] = 1
Color 3: [14144] = 1
Color 3: [14091] = 1
Color 3: [14040] = 1
Color 3: [13918] = 1
Color 3: [13804] = 1
Color 3: [13757] = 1
Color 3: [13688] = 1
Color 3: [13682] = 1
Color 3: [13674] = 2
Color 3: [13659] = 1
Color 3: [13654] = 1
Color 3: [13633] = 1
Color 3: [13620] = 1
Color 3: [13566] = 1
Color 3: [13538] = 1
Color 3: [13483] = 2
Color 3: [13380] = 1
Color 3: [13361] = 118
Color 3: [13326] = 1
Color 3: [13313] = 1718105

Thanks in advance, Flössie

Beep6581 commented 6 years ago

Here are histograms from RawDigger for 20180320-183828-1.cr2 with and without black level subtraction: 20180320-183828-1-rawdata with black subtraction.csv.txt 20180320-183828-1-rawdata no black subtraction.csv.txt

Maybe relevant: https://github.com/LibRaw/LibRaw/issues/82#issuecomment-271242569

iliasg commented 6 years ago

@Floessie I'll digg in in a while ..

Floessie commented 6 years ago

@iliasg Ilias, thanks. Until then, I can live with my patch. Morgan's LibRaw link seems indeed relevant.

iliasg commented 6 years ago

@Floessie Your 6D looks perfectly normal .. The ISO400 WhiteLevel at normal situations is 15283 (before BlackLevel subtraction) but with wide apperture settings, Canon operates a digital scaling so for f3.2 it's 15361 (see comment in camconst.json) https://github.com/Beep6581/RawTherapee/blob/806e1b32bb2e108abd32b448f664b341184f884f/rtengine/camconst.json#L511

It's common that these WhiteLevel values are violated by some outlier pixels (which normally are out of spec pixels) so here you have some sporadic higher than WL values and among them a group of some 120 pixels at each channel, at 15409 (13361/13362 after black level subtraction).

20180320-183828-1-full-5496x3670 See the gaps of the digital scaling ..

RT uses

I think HDRmerge should use the same .. ;)

Beep6581 commented 6 years ago

@iliasg you wrote,

It's common that these WhiteLevel values are violated by some outlier pixels

While @LibRaw wrote,

So to avoid magenta cast you need to set data cut limit based on histogram or on real data maximum

Since the brightest pixel value cannot be trusted to be a safe data maximum, I suppose this means the white level detection should somehow know how to find the "real, safe" white, to ignore the outliers. Is that correct? If so, how? Do you just discard the top 1% or is there a better or standard way? Short of having humans do the measurement and come up with a camconst.json file, that is.

iliasg commented 6 years ago

@Beep6581 As you saw .. I used the histogram ;) .. I don't know any other reliable way .. Here we had less than 0.01% outliers .. but I guess if we know that we have a shot with alot of clipping the real clipping value is the one with the most population in the histogram ..

The safety of this detection for abnormal cases (LENR or some raws where there is no clippin) is variable

I guess you have your A7II in mind .. no ? ;)

Beep6581 commented 6 years ago

@iliasg well I had RT in mind considering we need to adopt LibrRaw or RawSpeed, but also just the general problem of detecting a safe white level automatically on a per-image basis, considering some images could have clipping, some not, some could have special types of pixels (e.g. Sony), some not.

heckflosse commented 6 years ago

@Floessie @Beep6581 @iliasg

integrating hdrmerge into rt would solve this issue for us. That would also allow raw ca correction and so in prior to merge.

Should we try to do this for rt 5.6?

Beep6581 commented 6 years ago

That would be nice, but I would prioritize the pipeline unification, adopting Exiv2 and supplementing dcraw with one of the other options (in my mind, all three are critical). Besides, integrating HDRMerge in RT would require being able to paint a mask, and I don't think that will happen before the pipeline unification.

LibRaw commented 6 years ago

Starting from LibRaw 0.18, LibRaw library provides imgdata.color.linear_max[4] values based on camera metadata. It is not available for all vendors, but many vendors provide this data.

Floessie commented 6 years ago

@iliasg

Your 6D looks perfectly normal ..

Phew! That's to me the most important outcome of this discussion. :relieved:

I think HDRmerge should use the same .. ;)

That's not something we can do overnight (except for Alberto, perhaps).

@heckflosse

integrating hdrmerge into rt would solve this issue for us.

I use HDRMerge to archive my bracketed shots. Furthermore the batch mode is handy. How would one integrate that into RT?

@iliasg

As you saw .. I used the histogram ;) .. I don't know any other reliable way .. Here we had less than 0.01% outliers .. but I guess if we know that we have a shot with alot of clipping the real clipping value is the one with the most population in the histogram ..

With my patch above I search the histogram from top to bottom for the first value with more than one per mill of pixels. One could also take let's say the top tenth (relative to Image::max) of the histogram to look for the maximum...

Anyway, thanks for all the pondering.

Best, Flössie

Floessie commented 6 years ago

Fixed in #134. Closing.

heckflosse commented 6 years ago

@Floessie Reopened because I found a bug: Merging the raw files from https://discuss.pixls.us/t/play-raw-warning-baroque-can-harm-your-eyes/7922 only the darkest will be used.

Edit: In this case the blue values even in brightest image are very dark. Your method leads to a satThreshold of 3968 for red and green and 108 for blue. Because the minimum of the values is used, 108 will be used. That leads to the usage of darkest raw...

heckflosse commented 6 years ago

@Floessie Using the max instead of the min fixes the issue. Please test also with your stack.

diff --git a/src/ImageStack.cpp b/src/ImageStack.cpp
index 78c6888..6deb4ad 100644
--- a/src/ImageStack.cpp
+++ b/src/ImageStack.cpp
@@ -103,13 +103,13 @@ void ImageStack::calculateSaturationLevel(const RawParameters & params, bool use
         }
     }

-    satThreshold = params.max == 0 ? maxPerColor[0] : params.max;
-    for (int c = 0; c < 4; ++c) {
-        if (maxPerColor[c] > 0 && maxPerColor[c] < satThreshold) {
-            satThreshold = maxPerColor[c];
-        }
-    }

+    uint16_t maxPerColors = std::max(maxPerColor[0], std::max(maxPerColor[1],std::max(maxPerColor[2], maxPerColor[3])));
+    satThreshold = params.max == 0 ? maxPerColors : params.max;
+
+    if(maxPerColors > 0) {
+        satThreshold = std::min(satThreshold, maxPerColors);
+    }
     if(useCustomWl) {
         Log::debug( "Using custom white level ", params.max );
     }
heckflosse commented 6 years ago

@Floessie I already commited my patch. Please test nevertheless with your stack

Floessie commented 6 years ago

@heckflosse I'll test it later this week. Was this a bug, I accidentally introduced? If so: Sorry!

heckflosse commented 6 years ago

@Floessie The bug was introduced by your pr. But no need to apologize. I also did not detect it in my tests ;-)

In a good contribution, there are at least two contributors responsible. The one who introduced the bug and the one who did not detect it ;-) :+1:

Floessie commented 6 years ago

@heckflosse

The bug was introduced by your pr.

Ingo, I must be missing something. The minimum of maxPerColor[] has always been picked, and my PR didn't change that.

Best, Flössie

heckflosse commented 6 years ago

@Floessie Yes, that's correct. But now the maxPerColor[] values are different. If for example for one colour the 1% is reached at very dark areas because only few bright areas with this colour exist, the minimum of maxPerColor is very low. In the example above your method leads to a satThreshold of 3968 for red and green and 108 for blue and in consequence to a minimum of 108 which. Then only the darkest raw will be used.

Floessie commented 6 years ago

@heckflosse Aha, I see: The old values were better balanced, thus the minimum was fine. Now we pick lower values and don't want to settle on the lowest. :+1:

heckflosse commented 6 years ago

@Floessie Can we close the issue?

Floessie commented 6 years ago

@heckflosse Yes We Can!

heckflosse commented 6 years ago

Reopened because it's still buggy...

Beep6581 commented 6 years ago

Still buggy = #151

heckflosse commented 6 years ago

@Floessie I would like to revert this, because it introcuded more bugs than it fixed. We can let the issue open and start from scratch. What do you think?

Floessie commented 6 years ago

@heckflosse I don't mind. It fixes the problems I had, so I'll just keep on using my branch.

There must have been something odd with the calculation of satThreshold in the first place. My change might have been too simple to be true...

Floessie commented 6 years ago

@heckflosse Here's a fix:

diff --git a/src/ImageStack.cpp b/src/ImageStack.cpp
index 94677be..7d76f4d 100644
--- a/src/ImageStack.cpp
+++ b/src/ImageStack.cpp
@@ -111,6 +111,10 @@ void ImageStack::calculateSaturationLevel(const RawParameters & params, bool use
         satThreshold = std::min(satThreshold, maxPerColors);
     }

+    if (!useCustomWl) { // only scale when no custom white level was specified
+        satThreshold *= 0.99;
+    }
+
     Log::debug( "Using white level ", satThreshold );

     for (auto& i : images) {

Explanation: I originally removed the safety margin (*= .99), 'cause I deemed it no longer important. Now after bisecting and printing the values I figured out, the only difference in the resulting satThreshold was this multiplication.

HTH, Flössie

Beep6581 commented 6 years ago

@Floessie using your patch I can now successfully generate HDR DNGs from the Pentax K10D shots as well as from my Sony ILCE-7M2 shots without needing to manually set a custom white level.

heckflosse commented 6 years ago

@Floessie This is the log from merging 3 files from @Beep6581 using your latest patch. Because white level is very low (403) only the darkest image is used.

Invalid aperture: 0 replaced by aperture: f8
2018-05-02_soderasen_40.dng: 6024x4024 (6048x4024+0+0, by Sony ILCE-7M2, 200ISO 1/1000sec f8 EV:-14.9658
b4b4b4b4 RGBG, sat 16300, black 512, flip 5, wb: 2.46875 1 1.61328 0, cblack: 512 512 512 512
Invalid aperture: 0 replaced by aperture: f8
2018-05-02_soderasen_41.dng: 6024x4024 (6048x4024+0+0, by Sony ILCE-7M2, 200ISO 1/250sec f8 EV:-12.9658
b4b4b4b4 RGBG, sat 16300, black 512, flip 5, wb: 2.46875 1 1.61328 0, cblack: 512 512 512 512
Invalid aperture: 0 replaced by aperture: f8
2018-05-02_soderasen_42.dng: 6024x4024 (6048x4024+0+0, by Sony ILCE-7M2, 200ISO 1/60sec f8 EV:-10.9069
b4b4b4b4 RGBG, sat 16300, black 512, flip 5, wb: 2.46875 1 1.61328 0, cblack: 512 512 512 512
Load files: 2.8663 seconds
Using white level 403

Using this patch

diff --git a/src/ImageStack.cpp b/src/ImageStack.cpp
index 94677be..6074321 100644
--- a/src/ImageStack.cpp
+++ b/src/ImageStack.cpp
@@ -89,7 +89,7 @@ void ImageStack::calculateSaturationLevel(const RawParameters & params, bool use
         }
     }

-    const size_t threshold = width * height / 1000;
+    const size_t threshold = width * height / 10000;

     uint16_t maxPerColor[4] = {0, 0, 0, 0};

@@ -111,6 +111,10 @@ void ImageStack::calculateSaturationLevel(const RawParameters & params, bool use
         satThreshold = std::min(satThreshold, maxPerColors);
     }

+    if (!useCustomWl) { // only scale when no custom white level was specified
+        satThreshold *= 0.99;
+    }
+
     Log::debug( "Using white level ", satThreshold );

     for (auto& i : images) {

it works fine:

Invalid aperture: 0 replaced by aperture: f8
2018-05-02_soderasen_40.dng: 6024x4024 (6048x4024+0+0, by Sony ILCE-7M2, 200ISO 1/1000sec f8 EV:-14.9658
b4b4b4b4 RGBG, sat 16300, black 512, flip 5, wb: 2.46875 1 1.61328 0, cblack: 512 512 512 512
Invalid aperture: 0 replaced by aperture: f8
2018-05-02_soderasen_41.dng: 6024x4024 (6048x4024+0+0, by Sony ILCE-7M2, 200ISO 1/250sec f8 EV:-12.9658
b4b4b4b4 RGBG, sat 16300, black 512, flip 5, wb: 2.46875 1 1.61328 0, cblack: 512 512 512 512
Invalid aperture: 0 replaced by aperture: f8
2018-05-02_soderasen_42.dng: 6024x4024 (6048x4024+0+0, by Sony ILCE-7M2, 200ISO 1/60sec f8 EV:-10.9069
b4b4b4b4 RGBG, sat 16300, black 512, flip 5, wb: 2.46875 1 1.61328 0, cblack: 512 512 512 512
Load files: 2.87216 seconds
Using white level 15705
Floessie commented 6 years ago

@heckflosse Using 0.01% instead of 0.1% as a threshold should still be fine with my initial stack. :+1:

Beep6581 commented 6 years ago

Because white level is very low (403) only the darkest image is used.

:open_mouth:

The mask looked good, it used all source images, but I didn't examine the HDR DNG in RT to see whether the noise matched the mask...

heckflosse commented 6 years ago

@Beep6581

This is the mask with 0.1% threshold. Looks not good imho

grafik

And now with 0.01% threshold

grafik

Beep6581 commented 6 years ago

That first screenshot certainly looks bad, but it doesn't look like my image set. Using only @Floessie 's patch the mask in my images looks good: screenshot_20180620_142105

I tested again using 148c8e0 and it still looks good for both PEF and ARW: screenshot_20180620_142351 screenshot_20180620_142546

ARW set: https://filebin.net/kflpexzvnuf7gbkn

PEF set: http://rawtherapee.com/shared/test_images/amsterdam_moving_boat_1.pef http://rawtherapee.com/shared/test_images/amsterdam_moving_boat_2.pef http://rawtherapee.com/shared/test_images/amsterdam_moving_boat_3.pef