review mlsharpen

This commit is contained in:
Ingo Weyrich
2020-01-17 14:16:54 +01:00
parent 906347ab2d
commit d5c60475e5
3 changed files with 225 additions and 236 deletions

View File

@@ -30,7 +30,7 @@
#include "settings.h"
#include "sleef.h"
//#define BENCHMARK
#define BENCHMARK
#include "StopWatch.h"
using namespace std;
@@ -385,241 +385,6 @@ BENCHFUN
}
// To the extent possible under law, Manuel Llorens <manuelllorens@gmail.com>
// has waived all copyright and related or neighboring rights to this work.
// This work is published from: Spain.
// Thanks to Manuel for this excellent job (Jacques Desmis JDC or frej83)
void ImProcFunctions::MLsharpen (LabImage* lab)
{
// JD: this algorithm maximize clarity of images; it does not play on accutance. It can remove (partially) the effects of the AA filter)
// I think we can use this algorithm alone in most cases, or first to clarify image and if you want a very little USM (unsharp mask sharpening) after...
if (!params->sharpenEdge.enabled) {
return;
}
MyTime t1e, t2e;
t1e.set();
int offset, c, i, j, p, width2;
int width = lab->W, height = lab->H;
float *L, lumH, lumV, lumD1, lumD2, v, contrast, s;
float difL, difR, difT, difB, difLT, difRB, difLB, difRT, wH, wV, wD1, wD2, chmax[3];
float f1, f2, f3, f4;
float templab;
int iii, kkk;
width2 = 2 * width;
const float epsil = 0.01f; //prevent divide by zero
const float eps2 = 0.001f; //prevent divide by zero
float amount;
amount = params->sharpenEdge.amount / 100.0f;
if (amount < 0.00001f) {
return;
}
if (settings->verbose) {
printf ("SharpenEdge amount %f\n", amount);
}
L = new float[width * height];
chmax[0] = 8.0f;
chmax[1] = 3.0f;
chmax[2] = 3.0f;
int channels;
if (params->sharpenEdge.threechannels) {
channels = 0;
} else {
channels = 2;
}
if (settings->verbose) {
printf ("SharpenEdge channels %d\n", channels);
}
int passes = params->sharpenEdge.passes;
if (settings->verbose) {
printf ("SharpenEdge passes %d\n", passes);
}
for (p = 0; p < passes; p++)
for (c = 0; c <= channels; c++) { // c=0 Luminance only
#ifdef _OPENMP
#pragma omp parallel for private(offset) shared(L)
#endif
for (offset = 0; offset < width * height; offset++) {
int ii = offset / width;
int kk = offset - ii * width;
if (c == 0) {
L[offset] = lab->L[ii][kk] / 327.68f; // adjust to RT and to 0..100
} else if (c == 1) {
L[offset] = lab->a[ii][kk] / 327.68f;
} else { /*if (c==2) */
L[offset] = lab->b[ii][kk] / 327.68f;
}
}
#ifdef _OPENMP
#pragma omp parallel for private(j,i,iii,kkk, templab,offset,wH,wV,wD1,wD2,s,lumH,lumV,lumD1,lumD2,v,contrast,f1,f2,f3,f4,difT,difB,difL,difR,difLT,difLB,difRT,difRB) shared(lab,L,amount)
#endif
for(j = 2; j < height - 2; j++)
for(i = 2, offset = j * width + i; i < width - 2; i++, offset++) {
// weight functions
wH = eps2 + fabs(L[offset + 1] - L[offset - 1]);
wV = eps2 + fabs(L[offset + width] - L[offset - width]);
s = 1.0f + fabs(wH - wV) / 2.0f;
wD1 = eps2 + fabs(L[offset + width + 1] - L[offset - width - 1]) / s;
wD2 = eps2 + fabs(L[offset + width - 1] - L[offset - width + 1]) / s;
s = wD1;
wD1 /= wD2;
wD2 /= s;
// initial values
int ii = offset / width;
int kk = offset - ii * width;
if (c == 0) {
lumH = lumV = lumD1 = lumD2 = v = lab->L[ii][kk] / 327.68f;
} else if (c == 1) {
lumH = lumV = lumD1 = lumD2 = v = lab->a[ii][kk] / 327.68f;
} else { /* if (c==2) */
lumH = lumV = lumD1 = lumD2 = v = lab->b[ii][kk] / 327.68f;
}
// contrast detection
contrast = sqrt(fabs(L[offset + 1] - L[offset - 1]) * fabs(L[offset + 1] - L[offset - 1]) + fabs(L[offset + width] - L[offset - width]) * fabs(L[offset + width] - L[offset - width])) / chmax[c];
if (contrast > 1.0f) {
contrast = 1.0f;
}
// new possible values
if (((L[offset] < L[offset - 1]) && (L[offset] > L[offset + 1])) || ((L[offset] > L[offset - 1]) && (L[offset] < L[offset + 1]))) {
f1 = fabs(L[offset - 2] - L[offset - 1]);
f2 = fabs(L[offset - 1] - L[offset]);
f3 = fabs(L[offset - 1] - L[offset - width]) * fabs(L[offset - 1] - L[offset + width]);
f4 = sqrt(fabs(L[offset - 1] - L[offset - width2]) * fabs(L[offset - 1] - L[offset + width2]));
difL = f1 * f2 * f2 * f3 * f3 * f4;
f1 = fabs(L[offset + 2] - L[offset + 1]);
f2 = fabs(L[offset + 1] - L[offset]);
f3 = fabs(L[offset + 1] - L[offset - width]) * fabs(L[offset + 1] - L[offset + width]);
f4 = sqrt(fabs(L[offset + 1] - L[offset - width2]) * fabs(L[offset + 1] - L[offset + width2]));
difR = f1 * f2 * f2 * f3 * f3 * f4;
if ((difR > epsil) && (difL > epsil)) {
lumH = (L[offset - 1] * difR + L[offset + 1] * difL) / (difL + difR);
lumH = v * (1.f - contrast) + lumH * contrast;
}
}
if (((L[offset] < L[offset - width]) && (L[offset] > L[offset + width])) || ((L[offset] > L[offset - width]) && (L[offset] < L[offset + width]))) {
f1 = fabs(L[offset - width2] - L[offset - width]);
f2 = fabs(L[offset - width] - L[offset]);
f3 = fabs(L[offset - width] - L[offset - 1]) * fabs(L[offset - width] - L[offset + 1]);
f4 = sqrt(fabs(L[offset - width] - L[offset - 2]) * fabs(L[offset - width] - L[offset + 2]));
difT = f1 * f2 * f2 * f3 * f3 * f4;
f1 = fabs(L[offset + width2] - L[offset + width]);
f2 = fabs(L[offset + width] - L[offset]);
f3 = fabs(L[offset + width] - L[offset - 1]) * fabs(L[offset + width] - L[offset + 1]);
f4 = sqrt(fabs(L[offset + width] - L[offset - 2]) * fabs(L[offset + width] - L[offset + 2]));
difB = f1 * f2 * f2 * f3 * f3 * f4;
if ((difB > epsil) && (difT > epsil)) {
lumV = (L[offset - width] * difB + L[offset + width] * difT) / (difT + difB);
lumV = v * (1.f - contrast) + lumV * contrast;
}
}
if (((L[offset] < L[offset - 1 - width]) && (L[offset] > L[offset + 1 + width])) || ((L[offset] > L[offset - 1 - width]) && (L[offset] < L[offset + 1 + width]))) {
f1 = fabs(L[offset - 2 - width2] - L[offset - 1 - width]);
f2 = fabs(L[offset - 1 - width] - L[offset]);
f3 = fabs(L[offset - 1 - width] - L[offset - width + 1]) * fabs(L[offset - 1 - width] - L[offset + width - 1]);
f4 = sqrt(fabs(L[offset - 1 - width] - L[offset - width2 + 2]) * fabs(L[offset - 1 - width] - L[offset + width2 - 2]));
difLT = f1 * f2 * f2 * f3 * f3 * f4;
f1 = fabs(L[offset + 2 + width2] - L[offset + 1 + width]);
f2 = fabs(L[offset + 1 + width] - L[offset]);
f3 = fabs(L[offset + 1 + width] - L[offset - width + 1]) * fabs(L[offset + 1 + width] - L[offset + width - 1]);
f4 = sqrt(fabs(L[offset + 1 + width] - L[offset - width2 + 2]) * fabs(L[offset + 1 + width] - L[offset + width2 - 2]));
difRB = f1 * f2 * f2 * f3 * f3 * f4;
if ((difLT > epsil) && (difRB > epsil)) {
lumD1 = (L[offset - 1 - width] * difRB + L[offset + 1 + width] * difLT) / (difLT + difRB);
lumD1 = v * (1.f - contrast) + lumD1 * contrast;
}
}
if (((L[offset] < L[offset + 1 - width]) && (L[offset] > L[offset - 1 + width])) || ((L[offset] > L[offset + 1 - width]) && (L[offset] < L[offset - 1 + width]))) {
f1 = fabs(L[offset - 2 + width2] - L[offset - 1 + width]);
f2 = fabs(L[offset - 1 + width] - L[offset]);
f3 = fabs(L[offset - 1 + width] - L[offset - width - 1]) * fabs(L[offset - 1 + width] - L[offset + width + 1]);
f4 = sqrt(fabs(L[offset - 1 + width] - L[offset - width2 - 2]) * fabs(L[offset - 1 + width] - L[offset + width2 + 2]));
difLB = f1 * f2 * f2 * f3 * f3 * f4;
f1 = fabs(L[offset + 2 - width2] - L[offset + 1 - width]);
f2 = fabs(L[offset + 1 - width] - L[offset]) * fabs(L[offset + 1 - width] - L[offset]);
f3 = fabs(L[offset + 1 - width] - L[offset + width + 1]) * fabs(L[offset + 1 - width] - L[offset - width - 1]);
f4 = sqrt(fabs(L[offset + 1 - width] - L[offset + width2 + 2]) * fabs(L[offset + 1 - width] - L[offset - width2 - 2]));
difRT = f1 * f2 * f2 * f3 * f3 * f4;
if ((difLB > epsil) && (difRT > epsil)) {
lumD2 = (L[offset + 1 - width] * difLB + L[offset - 1 + width] * difRT) / (difLB + difRT);
lumD2 = v * (1.f - contrast) + lumD2 * contrast;
}
}
s = amount;
// avoid sharpening diagonals too much
if (((fabs(wH / wV) < 0.45f) && (fabs(wH / wV) > 0.05f)) || ((fabs(wV / wH) < 0.45f) && (fabs(wV / wH) > 0.05f))) {
s = amount / 3.0f;
}
// final mix
if ((wH != 0.0f) && (wV != 0.0f) && (wD1 != 0.0f) && (wD2 != 0.0f)) {
iii = offset / width;
kkk = offset - iii * width;
float provL = lab->L[iii][kkk] / 327.68f;
if(c == 0) {
if(provL < 92.f) {
templab = v * (1.f - s) + (lumH * wH + lumV * wV + lumD1 * wD1 + lumD2 * wD2) / (wH + wV + wD1 + wD2) * s;
} else {
templab = provL;
}
} else {
templab = v * (1.f - s) + (lumH * wH + lumV * wV + lumD1 * wD1 + lumD2 * wD2) / (wH + wV + wD1 + wD2) * s;
}
if (c == 0) {
lab->L[iii][kkk] = fabs(327.68f * templab); // fabs because lab->L always >0
} else if (c == 1) {
lab->a[iii][kkk] = 327.68f * templab ;
} else if (c == 2) {
lab->b[iii][kkk] = 327.68f * templab ;
}
}
}
}
delete [] L;
t2e.set();
if (settings->verbose) {
printf("SharpenEdge gradient %d usec\n", t2e.etime(t1e));
}
}
// To the extent possible under law, Manuel Llorens <manuelllorens@gmail.com>
// has waived all copyright and related or neighboring rights to this work.
// This code is licensed under CC0 v1.0, see license information at