Add mask method

rtengine/shmap.cc

@@ -64,6 +64,20 @@ void SHMap::fillLuminance( Imagefloat * img, float **luminance, double lumi[3] )
 
 }
 
+void SHMap::fillLuminanceL( float ** L, float **luminance) // fill with luminance
+{
+
+#ifdef _OPENMP
+    #pragma omp parallel for
+#endif
+
+    for (int i = 0; i < H; i++)
+        for (int j = 0; j < W; j++) {
+            luminance[i][j] = std::max(L[i][j], 0.f) ;//we can put here some enhancements Gamma, compression data,...
+        }
+
+}
+
 void SHMap::update (Imagefloat* img, double radius, double lumi[3], bool hq, int skip)
 {
 
@@ -210,6 +224,154 @@ void SHMap::update (Imagefloat* img, double radius, double lumi[3], bool hq, int
 
 }
 
+void SHMap::updateL (float** L, double radius, bool hq, int skip)
+{
+
+    if (!hq) {
+        fillLuminanceL( L, map);
+#ifdef _OPENMP
+        #pragma omp parallel
+#endif
+        {
+            gaussianBlur<float> (map, map, W, H, radius);
+        }
+    }
+
+    else
+
+    {
+        //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+        //experimental dirpyr shmap
+
+        float thresh = (100.f * radius); //1000;
+
+        // set up range function
+        // calculate size of Lookup table. That's possible because from a value k for all i>=k rangefn[i] will be exp(-10)
+        // So we use this fact and the automatic clip of lut to reduce the size of lut and the number of calculations to fill the lut
+        // In past this lut had only integer precision with rangefn[i] = 0 for all i>=k
+        // We set the last element to a small epsilon 1e-15 instead of zero to avoid divisions by zero
+        const int lutSize = thresh * sqrtf(10.f) + 1;
+        thresh *= thresh;
+        LUTf rangefn(lutSize);
+
+        for (int i = 0; i < lutSize - 1; i++) {
+            rangefn[i] = xexpf(-min(10.f, (static_cast<float>(i) * i) / thresh )); //*intfactor;
+        }
+
+        rangefn[lutSize - 1] = 1e-15f;
+
+        // We need one temporary buffer
+        float ** buffer = allocArray<float> (W, H);
+
+        // the final result has to be in map
+        // for an even number of levels that means: map => buffer, buffer => map
+        // for an odd number of levels that means: buffer => map, map => buffer, buffer => map
+        // so let's calculate the number of levels first
+        // There are at least two levels
+        int numLevels = 2;
+        int scale = 2;
+
+        while (skip * scale < 16) {
+            scale *= 2;
+            numLevels++;
+        }
+        //printf("numlev=%d\n",numLevels);
+        float ** dirpyrlo[2];
+
+        if(numLevels & 1) { // odd number of levels, start with buffer
+            dirpyrlo[0] = buffer;
+            dirpyrlo[1] = map;
+        } else { // even number of levels, start with map
+            dirpyrlo[0] = map;
+            dirpyrlo[1] = buffer;
+        }
+
+        fillLuminanceL( L, dirpyrlo[0]);
+
+        scale = 1;
+        int level = 0;
+        int indx = 0;
+        dirpyr_shmap(dirpyrlo[indx], dirpyrlo[1 - indx], W, H, rangefn, level, scale );
+        scale *= 2;
+        level ++;
+        indx = 1 - indx;
+
+        while (skip * scale < 16) {
+            dirpyr_shmap(dirpyrlo[indx], dirpyrlo[1 - indx], W, H, rangefn, level, scale );
+            scale *= 2;
+            level ++;
+            indx = 1 - indx;
+        }
+
+        dirpyr_shmap(dirpyrlo[indx], dirpyrlo[1 - indx], W, H, rangefn, level, scale );
+
+        freeArray<float>(buffer, H);
+
+        //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+        /*
+                // anti-alias filtering the result
+        #ifdef _OPENMP
+        #pragma omp for
+        #endif
+                for (int i=0; i<H; i++)
+                    for (int j=0; j<W; j++)
+                        if (i>0 && j>0 && i<H-1 && j<W-1)
+                            map[i][j] = (buffer[i-1][j-1]+buffer[i-1][j]+buffer[i-1][j+1]+buffer[i][j-1]+buffer[i][j]+buffer[i][j+1]+buffer[i+1][j-1]+buffer[i+1][j]+buffer[i+1][j+1])/9;
+                        else
+                            map[i][j] = buffer[i][j];
+        */
+
+    }
+
+    // update average, minimum, maximum
+    double _avg = 0.0f; // use double precision to gain precision especially at systems with few cores and big pictures (error for 36 MPixel on single core was about 8% with float)
+    min_f = 65535;
+    max_f = 0;
+#ifdef _OPENMP
+    #pragma omp parallel
+#endif
+    {
+        float _min_f = 65535.0f;
+        float _max_f = 0.0f;
+        float _val;
+#ifdef _OPENMP
+        #pragma omp for reduction(+:_avg) schedule(dynamic,16) nowait
+#endif
+
+        for (int i = 0; i < H; i++)
+            for (int j = 0; j < W; j++) {
+                _val = map[i][j];
+
+                if (_val < _min_f) {
+                    _min_f = _val;
+                }
+
+                if (_val > _max_f) {
+                    _max_f = _val;
+                }
+
+                _avg += _val;
+            }
+
+#ifdef _OPENMP
+        #pragma omp critical
+#endif
+        {
+            if(_min_f < min_f ) {
+                min_f = _min_f;
+            }
+
+            if(_max_f > max_f ) {
+                max_f = _max_f;
+            }
+        }
+    }
+    _avg /= ((H) * (W));
+    avg = _avg;
+
+}
+
+
 void SHMap::forceStat (float max_, float min_, float avg_)
 {