liz/rawTherapee
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Filmenegative core cleanup

Browse Source
This commit is contained in:
Flössie
2019-06-17 08:03:46 +02:00
parent 80f2b6a002
commit 54cc02eea9
6 changed files with 422 additions and 344 deletions
Download Patch File Download Diff File Expand all files Collapse all files

338
rtengine/filmnegativeproc.cc
View File

@@ -1,7 +1,7 @@
/*
* This file is part of RawTherapee.
*
* Copyright (c) 2004-2010 Gabor Horvath <hgabor@rawtherapee.com>
* Copyright (c) 2019 rom9
*
* RawTherapee is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
@@ -19,15 +19,17 @@
#include <cmath>
#include <iostream>
#include "rtengine.h"
#include "rawimagesource.h"
#include "mytime.h"
#include "procparams.h"
#ifdef _OPENMP
#include <omp.h>
#endif
#include "rawimagesource.h"
#include "mytime.h"
#include "opthelper.h"
#include "procparams.h"
#include "rt_algo.h"
#include "rtengine.h"
//#define BENCHMARK
#include "StopWatch.h"
@@ -37,168 +39,221 @@ namespace rtengine
extern const Settings* settings;
bool RawImageSource::channelsAvg(Coord spotPos, int spotSize, float avgs[3], const FilmNegativeParams &params)
}
namespace
{
avgs[0] = avgs[1] = avgs[2] = 0.f; // Channel averages
if(ri->getSensorType() != ST_BAYER && ri->getSensorType() != ST_FUJI_XTRANS)
bool channelsAvg(
const rtengine::RawImage* ri,
int width,
int height,
array2D<float>& rawData,
rtengine::Coord spotPos,
int spotSize,
const rtengine::procparams::FilmNegativeParams& params,
std::array<float, 3>& avgs
)
{
avgs = {}; // Channel averages
if (ri->getSensorType() != rtengine::ST_BAYER && ri->getSensorType() != rtengine::ST_FUJI_XTRANS) {
return false;
}
if (settings->verbose)
if (rtengine::settings->verbose) {
printf("Spot coord: x=%d y=%d\n", spotPos.x, spotPos.y);
}
int x1 = spotPos.x - spotSize / 2;
int x2 = spotPos.x + spotSize / 2;
int y1 = spotPos.y - spotSize / 2;
int y2 = spotPos.y + spotSize / 2;
const int half_spot_size = spotSize / 2;
if(x1<0 || x2>W || y1<0 || y2>H)
const int& x1 = spotPos.x - half_spot_size;
const int& x2 = spotPos.x + half_spot_size;
const int& y1 = spotPos.y - half_spot_size;
const int& y2 = spotPos.y + half_spot_size;
if (x1 < 0 || x2 > width || y1 < 0 || y2 > height) {
return false; // Spot goes outside bounds, bail out.
}
int pxCount[3] = {0}; // Per-channel sample counts
for(int c=spotPos.x-spotSize; c<spotPos.x+spotSize; c++) {
for(int r=spotPos.y-spotSize; r<spotPos.y+spotSize; r++) {
std::array<int, 3> pxCount = {}; // Per-channel sample counts
for (int c = spotPos.x - spotSize; c < spotPos.x + spotSize; ++c) {
for (int r = spotPos.y - spotSize; r < spotPos.y + spotSize; ++r) {
const int ch = ri->getSensorType() == rtengine::ST_BAYER ? ri->FC(r,c) : ri->XTRANSFC(r,c);
int ch = (ri->getSensorType() == ST_BAYER) ? FC(r,c) : ri->XTRANSFC(r,c);
++pxCount[ch];
pxCount[ch]++;
// If film negative is currently enabled, undo the effect by elevating to 1/exp,
// in order to sample the original, linear value
if(params.enabled)
avgs[ch] += powf(rawData[r][c], -1 / (ch==0 ? params.redExp : ch==1 ? params.greenExp : params.blueExp));
else
if (params.enabled) {
avgs[ch] += powf(rawData[r][c], -1.f / (ch == 0 ? params.redExp : ch == 1 ? params.greenExp : params.blueExp));
} else {
avgs[ch] += rawData[r][c];
}
}
}
for(int ch=0; ch<3; ch++)
avgs[ch] = avgs[ch] / (pxCount[ch]);
for (int ch = 0; ch < 3; ++ch) {
avgs[ch] /= pxCount[ch];
}
return true;
}
// Calculate logarithms in arbitrary base
float logBase(float base, float num) {
return log(num) / log(base);
}
bool RawImageSource::getFilmNegativeExponents(Coord2D spotA, Coord2D spotB, int tran, const FilmNegativeParams &currentParams, std::array<float, 3>& newExps)
bool rtengine::RawImageSource::getFilmNegativeExponents(Coord2D spotA, Coord2D spotB, int tran, const FilmNegativeParams &currentParams, std::array<float, 3>& newExps)
{
float clearVals[3], denseVals[3];
newExps = {
static_cast<float>(currentParams.redExp),
static_cast<float>(currentParams.greenExp),
static_cast<float>(currentParams.blueExp)
};
newExps[0] = currentParams.redExp;
newExps[1] = currentParams.greenExp;
newExps[2] = currentParams.blueExp;
constexpr int spotSize = 32; // TODO: Make this configurable?
int spotSize = 32; // TODO : make this confugurable ?
Coord spot;
std::array<float, 3> clearVals;
std::array<float, 3> denseVals;
// Sample first spot
transformPosition (spotA.x, spotA.y, tran, spot.x, spot.y);
if(!channelsAvg(spot, spotSize, clearVals, currentParams))
return false;
// Sample second spot
transformPosition (spotB.x, spotB.y, tran, spot.x, spot.y);
if(!channelsAvg(spot, spotSize, denseVals, currentParams))
return false;
// Detect which one is the dense spot, based on green channel
if(clearVals[1] < denseVals[1])
std::swap(clearVals, denseVals);
if (settings->verbose) {
printf("Clear film values: R=%g G=%g B=%g\n", clearVals[0], clearVals[1], clearVals[2]);
printf("Dense film values: R=%g G=%g B=%g\n", denseVals[0], denseVals[1], denseVals[2]);
transformPosition(spotA.x, spotA.y, tran, spot.x, spot.y);
if (!channelsAvg(ri, W, H, rawData, spot, spotSize, currentParams, clearVals)) {
return false;
}
float denseGreenRatio = clearVals[1] / denseVals[1];
// Sample second spot
transformPosition(spotB.x, spotB.y, tran, spot.x, spot.y);
if (!channelsAvg(ri, W, H, rawData, spot, spotSize, currentParams, denseVals)) {
return false;
}
// Detect which one is the dense spot, based on green channel
if (clearVals[1] < denseVals[1]) {
std::swap(clearVals, denseVals);
}
if (settings->verbose) {
printf("Clear film values: R=%g G=%g B=%g\n", clearVals[0], clearVals[1], clearVals[2]);
printf("Dense film values: R=%g G=%g B=%g\n", denseVals[0], denseVals[1], denseVals[2]);
}
const float denseGreenRatio = clearVals[1] / denseVals[1];
// Calculate logarithms in arbitrary base
const auto logBase =
[](float base, float num) -> float
{
return std::log(num) / std::log(base);
};
// Calculate exponents for each channel, based on the ratio between the bright and dark values,
// compared to the ratio in the reference channel (green)
for(int ch=0; ch<3; ch++)
if(ch==1)
for (int ch = 0; ch < 3; ++ch) {
if (ch == 1) {
newExps[ch] = 2.f; // Green is the reference channel
else
} else {
newExps[ch] = CLAMP(2.f * logBase(clearVals[ch] / denseVals[ch], denseGreenRatio), 0.3f, 6.f);
}
}
if (settings->verbose)
if (settings->verbose) {
printf("New exponents: R=%g G=%g B=%g\n", newExps[0], newExps[1], newExps[2]);
}
return true;
}
void RawImageSource::filmNegativeProcess(const procparams::FilmNegativeParams &params)
void rtengine::RawImageSource::filmNegativeProcess(const procparams::FilmNegativeParams &params)
{
// BENCHFUNMICRO
if(!params.enabled)
if (!params.enabled) {
return;
}
float exps[3] = { (float)params.redExp, (float)params.greenExp, (float)params.blueExp };
const std::array<float, 3> exps = {
static_cast<float>(params.redExp),
static_cast<float>(params.greenExp),
static_cast<float>(params.blueExp)
};
MyTime t1, t2, t3,t4, t5, t6;
MyTime t1, t2, t3,t4, t5;
t1.set();
// Channel vectors to calculate medians
std::vector<float> cvs[3] = {
std::vector<float>(),
std::vector<float>(),
std::vector<float>()
};
std::array<std::vector<float>, 3> cvs;
// Sample one every 5 pixels, and push the value in the appropriate channel vector.
// Chose an odd step, not multiple of the CFA size, to get a chance to visit each channel.
if(ri->getSensorType() == ST_BAYER) {
for (int row = 0; row < H; row+=5) {
for (int col = 0; col < W; col+=5) {
int c = FC(row, col); // three colors, 0=R, 1=G, 2=B
// Choose an odd step, not a multiple of the CFA size, to get a chance to visit each channel.
if (ri->getSensorType() == ST_BAYER) {
for (int row = 0; row < H; row += 5) {
for (int col = 0; col < W; col += 5) {
const int c = FC(row, col); // three colors: 0=R, 1=G, 2=B
cvs[c].push_back(rawData[row][col]);
}
}
} else if(ri->getSensorType() == ST_FUJI_XTRANS) {
for (int row = 0; row < H; row+=5) {
for (int col = 0; col < W; col+=5) {
int c = ri->XTRANSFC(row, col); // three colors, 0=R, 1=G, 2=B
}
else if (ri->getSensorType() == ST_FUJI_XTRANS) {
for (int row = 0; row < H; row += 5) {
for (int col = 0; col < W; col += 5) {
const int c = ri->XTRANSFC(row, col); // three colors: 0=R, 1=G, 2=B
cvs[c].push_back(rawData[row][col]);
}
}
}
const float MAX_OUT_VALUE = 65000.f;
constexpr float MAX_OUT_VALUE = 65000.f;
t2.set();
if (settings->verbose)
if (settings->verbose) {
printf("Median vector fill loop time us: %d\n", t2.etime(t1));
}
float medians[3]; // Channel median values
float mults[3] = { 1.f }; // Channel normalization multipliers
t2.set();
for (int c=0; c<3; c++) {
std::array<float, 3> medians; // Channel median values
std::array<float, 3> mults = {
1.f,
1.f,
1.f
}; // Channel normalization multipliers
for (int c = 0; c < 3; ++c) {
// Find median values for each channel
if(cvs[c].size() > 0) {
findMinMaxPercentile(&cvs[c][0], cvs[c].size(), 0.5f, medians[c], 0.5f, medians[c], true);
medians[c] = pow_F(max(medians[c], 1.f), -exps[c]);
if (!cvs[c].empty()) {
findMinMaxPercentile(cvs[c].data(), cvs[c].size(), 0.5f, medians[c], 0.5f, medians[c], true);
medians[c] = pow_F(rtengine::max(medians[c], 1.f), -exps[c]);
// Determine the channel multipler so that N times the median becomes 65k. This clips away
// the values in the dark border surrounding the negative (due to the film holder, for example),
// the reciprocal of which have blown up to stellar values.
mults[c] = MAX_OUT_VALUE / (medians[c] * 24);
mults[c] = MAX_OUT_VALUE / (medians[c] * 24.f);
}
}
t3.set();
if (settings->verbose) {
printf("Sample count : %lu, %lu, %lu\n", cvs[0].size(), cvs[1].size(), cvs[2].size());
printf("Medians : %g %g %g\n", medians[0], medians[1], medians[2] );
printf("Computed multipliers : %g %g %g\n", mults[0], mults[1], mults[2] );
printf("Sample count: %zu, %zu, %zu\n", cvs[0].size(), cvs[1].size(), cvs[2].size());
printf("Medians: %g %g %g\n", medians[0], medians[1], medians[2] );
printf("Computed multipliers: %g %g %g\n", mults[0], mults[1], mults[2] );
printf("Median calc time us: %d\n", t3.etime(t2));
}
t3.set();
if (ri->getSensorType() == ST_BAYER) {
#ifdef __SSE2__
const vfloat onev = F2V(1.f);
const vfloat c65535v = F2V(65535.f);
#endif
if(ri->getSensorType() == ST_BAYER) {
#ifdef _OPENMP
#pragma omp parallel for schedule(dynamic, 16)
#endif
for (int row = 0; row < H; row ++) {
for (int row = 0; row < H; ++row) {
int col = 0;
// Exponents are expressed as positive in the parameters, so negate them in order
// to get the reciprocals. Avoid trouble with zeroes, minimum pixel value is 1.
@@ -209,21 +264,24 @@ void RawImageSource::filmNegativeProcess(const procparams::FilmNegativeParams &p
#ifdef __SSE2__
const vfloat expsv = _mm_setr_ps(exps0, exps1, exps0, exps1);
const vfloat multsv = _mm_setr_ps(mult0, mult1, mult0, mult1);
const vfloat onev = F2V(1.f);
const vfloat c65535v = F2V(65535.f);
for (; col < W - 3; col+=4) {
for (; col < W - 3; col += 4) {
STVFU(rawData[row][col], vminf(multsv * pow_F(vmaxf(LVFU(rawData[row][col]), onev), expsv), c65535v));
}
#endif // __SSE2__
for (; col < W - 1; col+=2) {
rawData[row][col] = rtengine::min(mult0 * pow_F(max(rawData[row][col], 1.f), exps0), 65535.f);
rawData[row][col + 1] = rtengine::min(mult1 * pow_F(max(rawData[row][col + 1], 1.f), exps1), 65535.f);
for (; col < W - 1; col += 2) {
rawData[row][col] = rtengine::min(mult0 * pow_F(rtengine::max(rawData[row][col], 1.f), exps0), 65535.f);
rawData[row][col + 1] = rtengine::min(mult1 * pow_F(rtengine::max(rawData[row][col + 1], 1.f), exps1), 65535.f);
}
if (col < W) {
rawData[row][col] = rtengine::min(mult0 * pow_F(max(rawData[row][col], 1.f), exps0), 65535.f);
rawData[row][col] = rtengine::min(mult0 * pow_F(rtengine::max(rawData[row][col], 1.f), exps0), 65535.f);
}
}
} else if(ri->getSensorType() == ST_FUJI_XTRANS) {
} else if (ri->getSensorType() == ST_FUJI_XTRANS) {
#ifdef __SSE2__
const vfloat onev = F2V(1.f);
const vfloat c65535v = F2V(65535.f);
#endif
#ifdef _OPENMP
#pragma omp parallel for schedule(dynamic, 16)
#endif
@@ -231,8 +289,22 @@ void RawImageSource::filmNegativeProcess(const procparams::FilmNegativeParams &p
int col = 0;
// Exponents are expressed as positive in the parameters, so negate them in order
// to get the reciprocals. Avoid trouble with zeroes, minimum pixel value is 1.
const float expsc[6] = {-exps[ri->XTRANSFC(row, 0)], -exps[ri->XTRANSFC(row, 1)], -exps[ri->XTRANSFC(row, 2)], -exps[ri->XTRANSFC(row, 3)], -exps[ri->XTRANSFC(row, 4)], -exps[ri->XTRANSFC(row, 5)]};
const float multsc[6] = {mults[ri->XTRANSFC(row, 0)], mults[ri->XTRANSFC(row, 1)], mults[ri->XTRANSFC(row, 2)], mults[ri->XTRANSFC(row, 3)], mults[ri->XTRANSFC(row, 4)], mults[ri->XTRANSFC(row, 5)]};
const std::array<float, 6> expsc = {
-exps[ri->XTRANSFC(row, 0)],
-exps[ri->XTRANSFC(row, 1)],
-exps[ri->XTRANSFC(row, 2)],
-exps[ri->XTRANSFC(row, 3)],
-exps[ri->XTRANSFC(row, 4)],
-exps[ri->XTRANSFC(row, 5)]
};
const std::array<float, 6> multsc = {
mults[ri->XTRANSFC(row, 0)],
mults[ri->XTRANSFC(row, 1)],
mults[ri->XTRANSFC(row, 2)],
mults[ri->XTRANSFC(row, 3)],
mults[ri->XTRANSFC(row, 4)],
mults[ri->XTRANSFC(row, 5)]
};
#ifdef __SSE2__
const vfloat expsv0 = _mm_setr_ps(expsc[0], expsc[1], expsc[2], expsc[3]);
const vfloat expsv1 = _mm_setr_ps(expsc[4], expsc[5], expsc[0], expsc[1]);
@@ -240,81 +312,75 @@ void RawImageSource::filmNegativeProcess(const procparams::FilmNegativeParams &p
const vfloat multsv0 = _mm_setr_ps(multsc[0], multsc[1], multsc[2], multsc[3]);
const vfloat multsv1 = _mm_setr_ps(multsc[4], multsc[5], multsc[0], multsc[1]);
const vfloat multsv2 = _mm_setr_ps(multsc[2], multsc[3], multsc[4], multsc[5]);
const vfloat onev = F2V(1.f);
const vfloat c65535v = F2V(65535.f);
for (; col < W - 11; col+=12) {
for (; col < W - 11; col += 12) {
STVFU(rawData[row][col], vminf(multsv0 * pow_F(vmaxf(LVFU(rawData[row][col]), onev), expsv0), c65535v));
STVFU(rawData[row][col + 4], vminf(multsv1 * pow_F(vmaxf(LVFU(rawData[row][col + 4]), onev), expsv1), c65535v));
STVFU(rawData[row][col + 8], vminf(multsv2 * pow_F(vmaxf(LVFU(rawData[row][col + 8]), onev), expsv2), c65535v));
}
#endif // __SSE2__
for (; col < W - 5; col+=6) {
for (; col < W - 5; col += 6) {
for (int c = 0; c < 6; ++c) {
rawData[row][col + c] = rtengine::min(multsc[c] * pow_F(max(rawData[row][col + c], 1.f), expsc[c]), 65535.f);
rawData[row][col + c] = rtengine::min(multsc[c] * pow_F(rtengine::max(rawData[row][col + c], 1.f), expsc[c]), 65535.f);
}
}
for (int c = 0; col < W; col++, c++) {
rawData[row][col + c] = rtengine::min(multsc[c] * pow_F(max(rawData[row][col + c], 1.f), expsc[c]), 65535.f);
rawData[row][col + c] = rtengine::min(multsc[c] * pow_F(rtengine::max(rawData[row][col + c], 1.f), expsc[c]), 65535.f);
}
}
}
t4.set();
if (settings->verbose) {
printf("Pow loop time us: %d\n", t4.etime(t3));
}
t4.set();
if (settings->verbose)
printf("Pow loop time us: %d\n", t4.etime(t3));
t5.set();
PixelsMap bitmapBads(W, H);
int totBP = 0; // Hold count of bad pixels to correct
if(ri->getSensorType() == ST_BAYER) {
if (ri->getSensorType() == ST_BAYER) {
#ifdef _OPENMP
#pragma omp parallel for reduction(+:totBP) schedule(dynamic,16)
#endif
for(int i = 0; i < H; i++)
for(int j = 0; j < W; j++) {
for (int i = 0; i < H; ++i) {
for (int j = 0; j < W; ++j) {
if (rawData[i][j] >= MAX_OUT_VALUE) {
bitmapBads.set(j, i);
totBP++;
++totBP;
}
}
if (totBP > 0) {
interpolateBadPixelsBayer( bitmapBads, rawData );
}
} else if(ri->getSensorType() == ST_FUJI_XTRANS) {
#ifdef _OPENMP
#pragma omp parallel for reduction(+:totBP) schedule(dynamic,16)
#endif
for(int i = 0; i < H; i++)
for(int j = 0; j < W; j++) {
if (rawData[i][j] >= MAX_OUT_VALUE) {
bitmapBads.set(j, i);
totBP++;
}
}
if (totBP > 0) {
interpolateBadPixelsXtrans( bitmapBads );
interpolateBadPixelsBayer(bitmapBads, rawData);
}
}
else if (ri->getSensorType() == ST_FUJI_XTRANS) {
#ifdef _OPENMP
#pragma omp parallel for reduction(+:totBP) schedule(dynamic,16)
#endif
for (int i = 0; i < H; ++i) {
for (int j = 0; j < W; ++j) {
if (rawData[i][j] >= MAX_OUT_VALUE) {
bitmapBads.set(j, i);
totBP++;
}
}
}
if (totBP > 0) {
interpolateBadPixelsXtrans(bitmapBads);
}
}
t5.set();
t6.set();
if (settings->verbose) {
printf("Bad pixels count: %d\n", totBP);
printf("Bad pixels interpolation time us: %d\n", t6.etime(t5));
printf("Bad pixels interpolation time us: %d\n", t5.etime(t4));
}
}
}