Update

MSRCP: Add option chroma_protect to attenuate chroma adjustment.
MSRCP: Intensity channel is calculated from average of R, G, B value
now, instead of weighted average.
Delete Specification.h.
Update README.md

Author: mawen1250

README.md

@@ -14,7 +14,7 @@ The Retinex theory and algorithm mainly aims at simulating the color constancy f
 
 The light perceived by visual receptors can be separated into illuminance and reflectance. Retinex estimates the illuminance and derive the reflectance from the light, the filtered result of which is an image represents the reflectance characteristics of the scene, regardless of the illuminance.
 
-Retinex is a very powerful filter in dynamic range compression, contrast enhancement, color constancy, etc.
+Retinex is a very powerful filter in dynamic range compression, contrast enhancement, color constancy, de-fog, etc.
 
 ## MSRCP
 
@@ -28,10 +28,12 @@ As MSRCP preserves chromaticity, it is excellent for dynamic range compression a
 
 This function accept 8-16bit integer Gray/YUV/RGB/YCoCg input. Sub-sampled format is not supported. If you want to process YUV420/YUV422 clip, convert it to YUV444 or RGB first.
 
+For processing in YUV444 and RGB, the filtering results are different. MSR is applied to intesity channel, which is Y for YUV444 input and (R+G+B)/3 for RGB input. Since Y is a weighted average of R, G, B, processing in YUV444 may produce imbalanced chromaticity preservation. Also when chroma_protect is larger than 1 (default 1.2), the saturation of YUV444 processing result will be different from that of RGB processing result.
+
 ### Usage
 
 ```python
-retinex.MSRCP(clip input, float[] sigmaS=[25,80,250], float lower_thr=0, float upper_thr=0, bool fulls, bool fulld=fulls)
+retinex.MSRCP(clip input, float[] sigmaS=[25,80,250], float lower_thr=0, float upper_thr=0, bool fulls, bool fulld=fulls, float chroma_protect=1.2)
 ```
 
 - input:<br />
@@ -58,23 +60,35 @@ retinex.MSRCP(clip input, float[] sigmaS=[25,80,250], float lower_thr=0, float u
 
 - fulld: (Default: fulls)<br />
     Determine the value range of output clip. True means full range/PC range, and False means limited range/TV range.<br />
-    Set different value for fulls and fulld will result in range conversion.<br />
-    When fulls and fulld are the same, it is safe to assign either True or False for any kinds of input clip except YUV. When fulls=False it will automatically determine the Floor and Ceil of input image, which may produce a stronger filter result under some circumstances.
+    Set different value for fulls and fulld will result in range conversion.
+
+- chroma_protect: (Default: 1.2)<br />
+    The base of log function to attenuate chroma adjustment. It could avoid extreme chroma amplifying, while the saturation of the result is changed.<br />
+    Available range is [1, +inf), 1 means no attenuation.<br />
+    It is only available for YUV/YCoCg input.
 
 ### Example
 
-TV range YUV420P8 input, filter in TV range YUV444P16
+TV range YUV420P8 input, filtered in TV range YUV444P16 with chroma protect, output TV range YUV444P16
 
 ```python
 v = core.fmtc.resample(v, csp=vs.YUV444P16)
-v = core.retinex.MSRCP(v)
+v = core.retinex.MSRCP(v, chroma_protect=1.2)
 ```
 
-JPEG image(PC range YUV420P8 with MPEG-1 chroma placement) input, filter in PC range RGB48
+JPEG image(PC range YUV420P8 with MPEG-1 chroma placement) input, filtered in PC range YUV444P16 without chroma protect, output PC range RGB48
 
 ```python
 i = core.lsmas.LWLibavSource(r'Image.jpg')
 i = core.fmtc.resample(i, csp=vs.YUV444P16, fulls=True, cplace="MPEG1")
+i = core.retinex.MSRCP(i, fulls=True, chroma_protect=1)
 i = core.fmtc.matrix(i, mat="601", fulls=True, csp=vs.RGB48)
+```
+
+PNG image(PC range RGB24) input, filtered in PC range RGB48, output PC range RGB48
+
+```python
+i = core.lsmas.LWLibavSource(r'Image.png')
+i = core.fmtc.bitdepth(i, bits=16)
 i = core.retinex.MSRCP(i)
 ```

include/MSR.h

@@ -37,8 +37,8 @@ class MSRData
     std::vector<double> sigma;
     double lower_thr = 0;
     double upper_thr = 0;
-    int fulls = 1;
-    int fulld = fulls;
+    bool fulls = true;
+    bool fulld = fulls;
 
     int process[3];
 
@@ -60,9 +60,9 @@ class MSRData
     void fulls_select()
     {
         if (vi->format->colorFamily == cmGray || vi->format->colorFamily == cmYUV)
-            fulls = 0;
+            fulls = false;
         else if (vi->format->colorFamily == cmRGB || vi->format->colorFamily == cmYCoCg)
-            fulls = 1;
+            fulls = true;
     }
 };
 

include/MSRCP.h

@@ -22,7 +22,6 @@
 
 
 #include "Helper.h"
-#include "Specification.h"
 #include "MSR.h"
 
 
@@ -33,19 +32,13 @@ class MSRCPData
     : public MSRData
 {
 public:
-    ColorMatrix ColorMatrix_ = ColorMatrix::Unspecified;
+    double chroma_protect = 1.2;
 
 public:
     MSRCPData(const VSAPI *_vsapi = nullptr)
         : MSRData(_vsapi) {}
 
     ~MSRCPData() {}
-
-    void ColorMatrix_select()
-    {
-        if (ColorMatrix_ == ColorMatrix::Unspecified)
-            ColorMatrix_ = ColorMatrix_Default(vi->width, vi->height);
-    }
 };
 
 
@@ -71,14 +64,14 @@ void Retinex_MSRCP(VSFrameRef * dst, const VSFrameRef * src, const VSAPI * vsapi
     int sNeutral = 128 << (bps - 8);
     T sCeil = (1 << bps) - 1;
     //T sCeilC = (1 << bps) - 1;
-    T sRange = sCeil - sFloor;
+    T sRange = d.fulls ? (1 << bps) - 1 : 219 << (bps - 8);
     T sRangeC = d.fulls ? (1 << bps) - 1 : 224 << (bps - 8);
     T dFloor = d.fulld ? 0 : 16 << (bps - 8);
     //T dFloorC = d.fulld ? 0 : 16 << (bps - 8);
     int dNeutral = 128 << (bps - 8);
     T dCeil = d.fulld ? (1 << bps) - 1 : 235 << (bps - 8);
     //T dCeilC = d.fulld ? (1 << bps) - 1 : 240 << (bps - 8);
-    T dRange = dCeil - dFloor;
+    T dRange = d.fulld ? (1 << bps) - 1 : 219 << (bps - 8);
     T dRangeC = d.fulld ? (1 << bps) - 1 : 224 << (bps - 8);
     FLType sFloorFL = static_cast<FLType>(sFloor);
     //FLType sFloorCFL = static_cast<FLType>(sFloorC);
@@ -140,13 +133,11 @@ void Retinex_MSRCP(VSFrameRef * dst, const VSFrameRef * src, const VSAPI * vsapi
             {
                 sFloor = min;
                 sCeil = max;
-                sRange = sCeil - sFloor;
                 sFloorFL = static_cast<FLType>(sFloor);
                 //sCeilFL = static_cast<FLType>(sCeil);
-                sRangeFL = static_cast<FLType>(sRange);
             }
 
-            gain = 1 / sRangeFL;
+            gain = 1 / static_cast<FLType>(sCeil - sFloor);
             for (j = 0; j < height; j++)
             {
                 i = stride * j;
@@ -157,25 +148,12 @@ void Retinex_MSRCP(VSFrameRef * dst, const VSFrameRef * src, const VSAPI * vsapi
 
         Retinex_MSR(odata, idata, d, height, width, stride);
 
-        if (d.fulld)
-        {
-            offset = FLType(0.5);
-            for (j = 0; j < height; j++)
-            {
-                i = stride * j;
-                for (upper = i + width; i < upper; i++)
-                    Ydstp[i] = static_cast<T>(odata[i] * dRangeFL + offset);
-            }
-        }
-        else
+        offset = dFloorFL + FLType(0.5);
+        for (j = 0; j < height; j++)
         {
-            offset = dFloorFL + FLType(0.5);
-            for (j = 0; j < height; j++)
-            {
-                i = stride * j;
-                for (upper = i + width; i < upper; i++)
-                    Ydstp[i] = static_cast<T>(odata[i] * dRangeFL + offset);
-            }
+            i = stride * j;
+            for (upper = i + width; i < upper; i++)
+                Ydstp[i] = static_cast<T>(odata[i] * dRangeFL + offset);
         }
     }
     else if (fi->colorFamily == cmRGB)
@@ -190,17 +168,14 @@ void Retinex_MSRCP(VSFrameRef * dst, const VSFrameRef * src, const VSAPI * vsapi
         Bsrcp = reinterpret_cast<const T *>(vsapi->getReadPtr(src, 2));
         Bdstp = reinterpret_cast<T *>(vsapi->getWritePtr(dst, 2));
 
-        FLType Kr, Kg, Kb;
-        ColorMatrix_Parameter(d.ColorMatrix_, Kr, Kg, Kb);
-
         if (d.fulls)
         {
-            gain = 1 / sRangeFL;
+            gain = 1 / (sRangeFL * 3);
             for (j = 0; j < height; j++)
             {
                 i = stride * j;
                 for (upper = i + width; i < upper; i++)
-                    idata[i] = (Kr*Rsrcp[i] + Kg*Gsrcp[i] + Kb*Bsrcp[i]) * gain;
+                    idata[i] = (Rsrcp[i] + Gsrcp[i] + Bsrcp[i]) * gain;
             }
         }
         else
@@ -222,18 +197,17 @@ void Retinex_MSRCP(VSFrameRef * dst, const VSFrameRef * src, const VSAPI * vsapi
             {
                 sFloor = min;
                 sCeil = max;
-                sRange = sCeil - sFloor;
                 sFloorFL = static_cast<FLType>(sFloor);
                 //sCeilFL = static_cast<FLType>(sCeil);
-                sRangeFL = static_cast<FLType>(sRange);
             }
 
-            gain = 1 / sRangeFL;
+            offset = sFloorFL * -3;
+            gain = 1 / (static_cast<FLType>(sCeil - sFloor) * 3);
             for (j = 0; j < height; j++)
             {
                 i = stride * j;
                 for (upper = i + width; i < upper; i++)
-                    idata[i] = (Kr*Rsrcp[i] + Kg*Gsrcp[i] + Kb*Bsrcp[i] - sFloorFL) * gain;
+                    idata[i] = (Rsrcp[i] + Gsrcp[i] + Bsrcp[i] + offset) * gain;
             }
         }
 
@@ -247,7 +221,8 @@ void Retinex_MSRCP(VSFrameRef * dst, const VSFrameRef * src, const VSAPI * vsapi
                 i = stride * j;
                 for (upper = i + width; i < upper; i++)
                 {
-                    gain = Min(sRangeFL / Max(Rsrcp[i], Max(Gsrcp[i], Bsrcp[i])), idata[i] <= 0 ? 1 : odata[i] / idata[i]);
+                    gain = idata[i] <= 0 ? 1 : odata[i] / idata[i];
+                    gain = Min(sRangeFL / Max(Rsrcp[i], Max(Gsrcp[i], Bsrcp[i])), gain);
                     Rdstp[i] = static_cast<T>(Rsrcp[i] * gain + offset);
                     Gdstp[i] = static_cast<T>(Gsrcp[i] * gain + offset);
                     Bdstp[i] = static_cast<T>(Bsrcp[i] * gain + offset);
@@ -263,7 +238,8 @@ void Retinex_MSRCP(VSFrameRef * dst, const VSFrameRef * src, const VSAPI * vsapi
                 i = stride * j;
                 for (upper = i + width; i < upper; i++)
                 {
-                    gain = Min(sRangeFL / Max(Rsrcp[i], Max(Gsrcp[i], Bsrcp[i])), idata[i] <= 0 ? 1 : odata[i] / idata[i]) * scale;
+                    gain = idata[i] <= 0 ? 1 : odata[i] / idata[i];
+                    gain = Min(sRangeFL / Max(Rsrcp[i], Max(Gsrcp[i], Bsrcp[i])), gain) * scale;
                     Rdstp[i] = static_cast<T>((Rsrcp[i] - sFloor) * gain + offset);
                     Gdstp[i] = static_cast<T>((Gsrcp[i] - sFloor) * gain + offset);
                     Bdstp[i] = static_cast<T>((Bsrcp[i] - sFloor) * gain + offset);
@@ -312,13 +288,11 @@ void Retinex_MSRCP(VSFrameRef * dst, const VSFrameRef * src, const VSAPI * vsapi
             {
                 sFloor = min;
                 sCeil = max;
-                sRange = sCeil - sFloor;
                 sFloorFL = static_cast<FLType>(sFloor);
                 //sCeilFL = static_cast<FLType>(sCeil);
-                sRangeFL = static_cast<FLType>(sRange);
             }
 
-            gain = 1 / sRangeFL;
+            gain = 1 / static_cast<FLType>(sCeil - sFloor);
             for (j = 0; j < height; j++)
             {
                 i = stride * j;
@@ -329,35 +303,34 @@ void Retinex_MSRCP(VSFrameRef * dst, const VSFrameRef * src, const VSAPI * vsapi
 
         Retinex_MSR(odata, idata, d, height, width, stride);
 
-        if (dRangeCFL == sRangeCFL)
-        {
-            offset = dNeutralFL + FLType(0.5);
-            for (j = 0; j < height; j++)
-            {
-                i = stride * j;
-                for (upper = i + width; i < upper; i++)
-                {
-                    gain = Min(sRangeC2FL / Max(Abs(Usrcp[i] - sNeutral), Abs(Vsrcp[i] - sNeutral)), idata[i] <= 0 ? 1 : odata[i] / idata[i]);
-                    Ydstp[i] = static_cast<T>(odata[i] * dRangeFL + dFloorFL + FLType(0.5));
-                    Udstp[i] = static_cast<T>((Usrcp[i] - sNeutral) * gain + offset);
-                    Vdstp[i] = static_cast<T>((Vsrcp[i] - sNeutral) * gain + offset);
-                }
-            }
-        }
+        FLType chroma_protect_mul1 = static_cast<FLType>(d.chroma_protect - 1);
+        FLType chroma_protect_mul2 = static_cast<FLType>(1 / log(d.chroma_protect));
+
+        int Uval, Vval;
+        scale = dRangeCFL / sRangeCFL;
+        if (d.fulld)
+            offset = dNeutralFL + FLType(0.499999);
         else
-        {
-            scale = dRangeCFL / sRangeCFL;
             offset = dNeutralFL + FLType(0.5);
-            for (j = 0; j < height; j++)
+
+        for (j = 0; j < height; j++)
+        {
+            i = stride * j;
+            for (upper = i + width; i < upper; i++)
             {
-                i = stride * j;
-                for (upper = i + width; i < upper; i++)
-                {
-                    gain = Min(sRangeC2FL / Max(Abs(Usrcp[i] - sNeutral), Abs(Vsrcp[i] - sNeutral)), idata[i] <= 0 ? 1 : odata[i] / idata[i]) * scale;
-                    Ydstp[i] = static_cast<T>(odata[i] * dRangeFL + dFloorFL + FLType(0.5));
-                    Udstp[i] = static_cast<T>((Usrcp[i] - sNeutral) * gain + offset);
-                    Vdstp[i] = static_cast<T>((Vsrcp[i] - sNeutral) * gain + offset);
-                }
+                if (d.chroma_protect > 1)
+                    gain = idata[i] <= 0 ? 1 : log(odata[i] / idata[i] * chroma_protect_mul1 + 1) * chroma_protect_mul2;
+                else
+                    gain = idata[i] <= 0 ? 1 : odata[i] / idata[i];
+                Uval = Usrcp[i] - sNeutral;
+                Vval = Vsrcp[i] - sNeutral;
+                if (dRangeCFL == sRangeCFL)
+                    gain = Min(sRangeC2FL / Max(Abs(Uval), Abs(Vval)), gain);
+                else
+                    gain = Min(sRangeC2FL / Max(Abs(Uval), Abs(Vval)), gain) * scale;
+                Ydstp[i] = static_cast<T>(odata[i] * dRangeFL + dFloorFL + FLType(0.5));
+                Udstp[i] = static_cast<T>(Uval * gain + offset);
+                Vdstp[i] = static_cast<T>(Vval * gain + offset);
             }
         }
     }