From ee4e56ce0094ef329552cad7f47f59919818e282 Mon Sep 17 00:00:00 2001 From: Ingo Date: Wed, 30 Sep 2015 17:45:36 +0200 Subject: [PATCH 01/16] Enable DCP options regardless of DCP source. Neutral profile disables DCP curve and look table. Fixes #2922 --- rtengine/dcp.cc | 10 +-- rtengine/dcp.h | 4 +- rtengine/procparams.cc | 4 +- rtgui/icmpanel.cc | 159 ++++++++++++++++++++++------------------- rtgui/icmpanel.h | 1 + 5 files changed, 94 insertions(+), 84 deletions(-) diff --git a/rtengine/dcp.cc b/rtengine/dcp.cc index 965ed47b2..9b3d852fc 100644 --- a/rtengine/dcp.cc +++ b/rtengine/dcp.cc @@ -814,7 +814,7 @@ const DCPProfile::HSBModify* DCPProfile::MakeHueSatMap(ColorTemp &wb, int prefer return aDeltas; } -DCPProfile::DCPProfile(Glib::ustring fname, bool isRTProfile) +DCPProfile::DCPProfile(Glib::ustring fname) { const int TIFFFloatSize = 4; const int TagColorMatrix1 = 50721, TagColorMatrix2 = 50722, TagProfileHueSatMapDims = 50937; @@ -989,7 +989,7 @@ DCPProfile::DCPProfile(Glib::ustring fname, bool isRTProfile) // Read tone curve points, if any, but disable to RTs own profiles tag = tagDir->getTag(TagProfileToneCurve); - if (tag != NULL && !isRTProfile) { + if (tag != NULL) { std::vector cPoints; cPoints.push_back(double(DCT_Spline)); // The first value is the curve type @@ -1760,7 +1760,7 @@ void DCPStore::init (Glib::ustring rtProfileDir) } } -DCPProfile* DCPStore::getProfile (Glib::ustring filename, bool isRTProfile) +DCPProfile* DCPStore::getProfile (Glib::ustring filename) { MyMutex::MyLock lock(mtx); @@ -1771,7 +1771,7 @@ DCPProfile* DCPStore::getProfile (Glib::ustring filename, bool isRTProfile) } // Add profile - profileCache[filename] = new DCPProfile(filename, isRTProfile); + profileCache[filename] = new DCPProfile(filename); return profileCache[filename]; } @@ -1783,7 +1783,7 @@ DCPProfile* DCPStore::getStdProfile(Glib::ustring camShortName) // Warning: do NOT use map.find(), since it does not seem to work reliably here for (std::map::iterator i = fileStdProfiles.begin(); i != fileStdProfiles.end(); i++) if (name2 == (*i).first) { - return getProfile((*i).second, true); + return getProfile((*i).second); } return NULL; diff --git a/rtengine/dcp.h b/rtengine/dcp.h index 0d3bcc029..7a960edc9 100644 --- a/rtengine/dcp.h +++ b/rtengine/dcp.h @@ -76,7 +76,7 @@ class DCPProfile void HSDApply(const HSDTableInfo &ti, const HSBModify *tableBase, float &h, float &s, float &v) const; public: - DCPProfile(Glib::ustring fname, bool isRTProfile); + DCPProfile(Glib::ustring fname); ~DCPProfile(); bool getHasToneCurve() @@ -122,7 +122,7 @@ public: bool isValidDCPFileName(Glib::ustring filename) const; - DCPProfile* getProfile(Glib::ustring filename, bool isRTProfile = false); + DCPProfile* getProfile(Glib::ustring filename); DCPProfile* getStdProfile(Glib::ustring camShortName); static DCPStore* getInstance(); diff --git a/rtengine/procparams.cc b/rtengine/procparams.cc index 66cb76b3b..32506cb0c 100644 --- a/rtengine/procparams.cc +++ b/rtengine/procparams.cc @@ -831,8 +831,8 @@ void ColorManagementParams::setDefaults() { input = "(cameraICC)"; blendCMSMatrix = false; - toneCurve = true; - applyLookTable = true; + toneCurve = false; + applyLookTable = false; applyBaselineExposureOffset = true; applyHueSatMap = true; dcpIlluminant = 0; diff --git a/rtgui/icmpanel.cc b/rtgui/icmpanel.cc index c6f11035b..e9e4e05ff 100644 --- a/rtgui/icmpanel.cc +++ b/rtgui/icmpanel.cc @@ -30,7 +30,7 @@ using namespace rtengine::procparams; extern Options options; -ICMPanel::ICMPanel () : FoldableToolPanel(this, "icm", M("TP_ICM_LABEL")), iunchanged(NULL), icmplistener(NULL), lastRefFilename("") +ICMPanel::ICMPanel () : FoldableToolPanel(this, "icm", M("TP_ICM_LABEL")), iunchanged(NULL), icmplistener(NULL), lastRefFilename(""), camName("") { isBatchMode = lastToneCurve = lastApplyLookTable = lastApplyBaselineExposureOffset = lastApplyHueSatMap = lastBlendCMSMatrix = lastgamfree = false; @@ -348,72 +348,76 @@ void ICMPanel::updateDCP (int dcpIlluminant, Glib::ustring dcp_name) dcpIll->set_sensitive (false); dcpFrame->set_sensitive(false); - if (ifromfile->get_active() && dcpStore->isValidDCPFileName(dcp_name)) { - DCPProfile* dcp = dcpStore->getProfile(dcp_name, false); + DCPProfile* dcp = NULL; - if (dcp) { - dcpFrame->set_sensitive(true); + if(dcp_name == "(cameraICC)") { + dcp = dcpStore->getStdProfile(camName); + } else if (ifromfile->get_active() && dcpStore->isValidDCPFileName(dcp_name)) { + dcp = dcpStore->getProfile(dcp_name); + } - if (dcp->getHasToneCurve()) { - ckbToneCurve->set_sensitive (true); + if (dcp) { + dcpFrame->set_sensitive(true); + + if (dcp->getHasToneCurve()) { + ckbToneCurve->set_sensitive (true); + } + + if (dcp->getHasLookTable()) { + ckbApplyLookTable->set_sensitive (true); + } + + if (dcp->getHasBaselineExposureOffset()) { + ckbApplyBaselineExposureOffset->set_sensitive (true); + } + + if (dcp->getHasHueSatMap()) { + ckbApplyHueSatMap->set_sensitive (true); + } + + int i1, i2; + double temp1, temp2; + bool willInterpolate; + dcp->getIlluminants(i1, temp1, i2, temp2, willInterpolate); + + if (willInterpolate) { + if (dcpTemperatures[0] != temp1 || dcpTemperatures[1] != temp2) { + char tempstr1[64], tempstr2[64]; + sprintf(tempstr1, "%.0fK", temp1); + sprintf(tempstr2, "%.0fK", temp2); + int curr_active = dcpIll->get_active_row_number(); + ignoreDcpSignal = true; + dcpIll->clear_items (); + dcpIll->append_text (M("TP_ICM_DCPILLUMINANT_INTERPOLATED")); + dcpIll->append_text (tempstr1); + dcpIll->append_text (tempstr2); + dcpTemperatures[0] = temp1; + dcpTemperatures[1] = temp2; + dcpIll->set_active (curr_active); + ignoreDcpSignal = false; } - if (dcp->getHasLookTable()) { - ckbApplyLookTable->set_sensitive (true); + if (dcpIlluminant > 2) { + dcpIlluminant = 0; } - if (dcp->getHasBaselineExposureOffset()) { - ckbApplyBaselineExposureOffset->set_sensitive (true); + if (dcpIll->get_active_row_number() == -1 && dcpIlluminant == -1) { + ignoreDcpSignal = true; + dcpIll->set_active(0); + ignoreDcpSignal = false; + } else if (dcpIlluminant >= 0 && dcpIlluminant != dcpIll->get_active_row_number()) { + ignoreDcpSignal = true; + dcpIll->set_active(dcpIlluminant); + ignoreDcpSignal = false; } - if (dcp->getHasHueSatMap()) { - ckbApplyHueSatMap->set_sensitive (true); - } - - int i1, i2; - double temp1, temp2; - bool willInterpolate; - dcp->getIlluminants(i1, temp1, i2, temp2, willInterpolate); - - if (willInterpolate) { - if (dcpTemperatures[0] != temp1 || dcpTemperatures[1] != temp2) { - char tempstr1[64], tempstr2[64]; - sprintf(tempstr1, "%.0fK", temp1); - sprintf(tempstr2, "%.0fK", temp2); - int curr_active = dcpIll->get_active_row_number(); - ignoreDcpSignal = true; - dcpIll->clear_items (); - dcpIll->append_text (M("TP_ICM_DCPILLUMINANT_INTERPOLATED")); - dcpIll->append_text (tempstr1); - dcpIll->append_text (tempstr2); - dcpTemperatures[0] = temp1; - dcpTemperatures[1] = temp2; - dcpIll->set_active (curr_active); - ignoreDcpSignal = false; - } - - if (dcpIlluminant > 2) { - dcpIlluminant = 0; - } - - if (dcpIll->get_active_row_number() == -1 && dcpIlluminant == -1) { - ignoreDcpSignal = true; - dcpIll->set_active(0); - ignoreDcpSignal = false; - } else if (dcpIlluminant >= 0 && dcpIlluminant != dcpIll->get_active_row_number()) { - ignoreDcpSignal = true; - dcpIll->set_active(dcpIlluminant); - ignoreDcpSignal = false; - } - - dcpIll->set_sensitive (true); - dcpIllLabel->set_sensitive (true); - } else { - if (dcpIll->get_active_row_number() != -1) { - ignoreDcpSignal = true; - dcpIll->set_active(-1); - ignoreDcpSignal = false; - } + dcpIll->set_sensitive (true); + dcpIllLabel->set_sensitive (true); + } else { + if (dcpIll->get_active_row_number() != -1) { + ignoreDcpSignal = true; + dcpIll->set_active(-1); + ignoreDcpSignal = false; } } } @@ -467,7 +471,7 @@ void ICMPanel::read (const ProcParams* pp, const ParamsEdited* pedited) } else if ((pp->icm.input == "(cameraICC)") && icameraICC->get_state() != Gtk::STATE_INSENSITIVE) { icameraICC->set_active (true); ckbBlendCMSMatrix->set_sensitive (true); - updateDCP(pp->icm.dcpIlluminant, ""); + updateDCP(pp->icm.dcpIlluminant, "(cameraICC)"); } else if ((pp->icm.input == "(cameraICC)") && icamera->get_state() != Gtk::STATE_INSENSITIVE && icameraICC->get_state() == Gtk::STATE_INSENSITIVE) { // this is the case when (cameraICC) is instructed by packaged profiles, but ICC file is not found // therefore falling back UI to explicitly reflect the (camera) option @@ -478,7 +482,7 @@ void ICMPanel::read (const ProcParams* pp, const ParamsEdited* pedited) // If neither (camera) nor (cameraICC) are available, as is the case when loading a non-raw, activate (embedded). iembedded->set_active (true); ckbBlendCMSMatrix->set_sensitive (false); - updateDCP(pp->icm.dcpIlluminant, ""); + updateDCP(pp->icm.dcpIlluminant, "(cameraICC)"); } else if ((pp->icm.input == "(camera)" || pp->icm.input == "") && icamera->get_state() != Gtk::STATE_INSENSITIVE) { icamera->set_active (true); ckbBlendCMSMatrix->set_sensitive (false); @@ -603,25 +607,29 @@ void ICMPanel::write (ProcParams* pp, ParamsEdited* pedited) pp->icm.freegamma = freegamma->get_active(); + DCPProfile* dcp = NULL; + if (ifromfile->get_active() && pp->icm.input.substr(0, 5) == "file:" && dcpStore->isValidDCPFileName(pp->icm.input.substr(5))) { - DCPProfile* dcp = dcpStore->getProfile(pp->icm.input.substr(5), false); + dcp = dcpStore->getProfile(pp->icm.input.substr(5)); + } else if(icameraICC->get_active()) { + dcp = dcpStore->getStdProfile(camName); + } - if (dcp) { - if (dcp->getHasToneCurve()) { - pp->icm.toneCurve = ckbToneCurve->get_active (); - } + if (dcp) { + if (dcp->getHasToneCurve()) { + pp->icm.toneCurve = ckbToneCurve->get_active (); + } - if (dcp->getHasLookTable()) { - pp->icm.applyLookTable = ckbApplyLookTable->get_active (); - } + if (dcp->getHasLookTable()) { + pp->icm.applyLookTable = ckbApplyLookTable->get_active (); + } - if (dcp->getHasBaselineExposureOffset()) { - pp->icm.applyBaselineExposureOffset = ckbApplyBaselineExposureOffset->get_active (); - } + if (dcp->getHasBaselineExposureOffset()) { + pp->icm.applyBaselineExposureOffset = ckbApplyBaselineExposureOffset->get_active (); + } - if (dcp->getHasHueSatMap()) { - pp->icm.applyHueSatMap = ckbApplyHueSatMap->get_active (); - } + if (dcp->getHasHueSatMap()) { + pp->icm.applyHueSatMap = ckbApplyHueSatMap->get_active (); } } @@ -880,6 +888,7 @@ void ICMPanel::setRawMeta (bool raw, const rtengine::ImageData* pMeta) icamera->set_active (raw); iembedded->set_active (!raw); icamera->set_sensitive (raw); + camName = pMeta->getCamera(); icameraICC->set_sensitive (raw && (iccStore->getStdProfile(pMeta->getCamera()) != NULL || dcpStore->getStdProfile(pMeta->getCamera()) != NULL)); iembedded->set_sensitive (!raw); diff --git a/rtgui/icmpanel.h b/rtgui/icmpanel.h index 856558dff..93828f5fd 100644 --- a/rtgui/icmpanel.h +++ b/rtgui/icmpanel.h @@ -93,6 +93,7 @@ private: double dcpTemperatures[2]; bool enableLastICCWorkDirChange; Glib::ustring lastRefFilename; + Glib::ustring camName; void updateDCP(int dcpIlluminant, Glib::ustring dcp_name); public: ICMPanel (); From aaf21544a3e85eda247c7a9dd1eac9621c3aec9d Mon Sep 17 00:00:00 2001 From: heckflosse Date: Thu, 1 Oct 2015 21:16:47 +0200 Subject: [PATCH 02/16] Speeduo for Hot/Dead pixel detection --- rtengine/rawimagesource.cc | 118 ++++++++++++++++++++++--------------- rtengine/sleefsseavx.c | 6 ++ 2 files changed, 77 insertions(+), 47 deletions(-) diff --git a/rtengine/rawimagesource.cc b/rtengine/rawimagesource.cc index b425a9f88..6197dba17 100644 --- a/rtengine/rawimagesource.cc +++ b/rtengine/rawimagesource.cc @@ -40,6 +40,7 @@ #include #endif #include "opthelper.h" +#include "StopWatch.h" namespace rtengine { @@ -946,89 +947,112 @@ int RawImageSource::interpolateBadPixelsXtrans( PixelsMap &bitmapBads ) */ int RawImageSource::findHotDeadPixels( PixelsMap &bpMap, float thresh, bool findHotPixels, bool findDeadPixels ) { - float varthresh = (20.0 * (thresh / 100.0) + 1.0 ); - - // counter for dead or hot pixels - int counter = 0; + StopWatch Stop1("findHotDeadPixels"); + float varthresh = (20.0 * (thresh / 100.0) + 1.0 ) / 24.f; // allocate temporary buffer float (*cfablur); cfablur = (float (*)) malloc (H * W * sizeof * cfablur); + // counter for dead or hot pixels + int counter = 0; + #pragma omp parallel { - #pragma omp for + #pragma omp for schedule(dynamic,16) nowait - for (int i = 0; i < H; i++) { - int iprev, inext, jprev, jnext; + for (int i = 2; i < H - 2; i++) { float p[9], temp; - if (i < 2) { - iprev = i + 2; - } else { - iprev = i - 2; - } + for (int j = 2; j < W - 2; j++) { - if (i > H - 3) { - inext = i - 2; - } else { - inext = i + 2; - } - - for (int j = 0; j < W; j++) { - if (j < 2) { - jprev = j + 2; - } else { - jprev = j - 2; - } - - if (j > W - 3) { - jnext = j - 2; - } else { - jnext = j + 2; - } - - med3x3(rawData[iprev][jprev], rawData[iprev][j], rawData[iprev][jnext], - rawData[i][jprev], rawData[i][j], rawData[i][jnext], - rawData[inext][jprev], rawData[inext][j], rawData[inext][jnext], temp); + med3x3(rawData[i - 2][j - 2], rawData[i - 2][j], rawData[i - 2][j + 2], + rawData[i][j - 2], rawData[i][j], rawData[i][j + 2], + rawData[i + 2][j - 2], rawData[i + 2][j], rawData[i + 2][j + 2], temp); cfablur[i * W + j] = rawData[i][j] - temp; } } + // process borders. Former version calculated the median using mirrored border which does not make sense because the original pixel loses weight + // Setting the difference between pixel and median for border pixels to zero should do the job not worse then former version + #pragma omp critical + { + for(int i = 0; i < 2; i++) { + for(int j = 0; j < W; j++) { + cfablur[i * W + j] = 0.f; + } + } + + for(int i = 2; i < H - 2; i++) { + for(int j = 0; j < 2; j++) { + cfablur[i * W + j] = 0.f; + } + + for(int j = W - 2; j < W; j++) { + cfablur[i * W + j] = 0.f; + } + } + + for(int i = H - 2; i < H; i++) { + for(int j = 0; j < W; j++) { + cfablur[i * W + j] = 0.f; + } + } + } + #pragma omp barrier #pragma omp for reduction(+:counter) schedule (dynamic,16) //cfa pixel heat/death evaluation - for (int rr = 0; rr < H; rr++) { - int top = max(0, rr - 2); - int bottom = min(H - 1, rr + 2); - int rrmWpcc = rr * W; + for (int rr = 2; rr < H - 2; rr++) { + int rrmWpcc = rr * W + 2; - for (int cc = 0; cc < W; cc++, rrmWpcc++) { + for (int cc = 2; cc < W - 2; cc++, rrmWpcc++) { //evaluate pixel for heat/death float pixdev = cfablur[rrmWpcc]; - if((!findDeadPixels) && pixdev <= 0) { + if(pixdev == 0.f) { continue; } - if((!findHotPixels) && pixdev >= 0) { + if((!findDeadPixels) && pixdev < 0) { + continue; + } + + if((!findHotPixels) && pixdev > 0) { continue; } pixdev = fabsf(pixdev); float hfnbrave = -pixdev; - int left = max(0, cc - 2); - int right = min(W - 1, cc + 2); - for (int mm = top; mm <= bottom; mm++) { - int mmmWpnn = mm * W + left; - for (int nn = left; nn <= right; nn++, mmmWpnn++) { +#ifdef __SSE2__ + // sum up 5*4 = 20 values using SSE + vfloat sum = vabsf(LVFU(cfablur[(rr - 2) * W + cc - 2])) + vabsf(LVFU(cfablur[(rr - 1) * W + cc - 2])); + sum += vabsf(LVFU(cfablur[(rr) * W + cc - 2])); + sum += vabsf(LVFU(cfablur[(rr + 1) * W + cc - 2])); + sum += vabsf(LVFU(cfablur[(rr + 2) * W + cc - 2])); + // horizontally add the values and add the result to hfnbrave + hfnbrave += vhadd(sum); + + // add remaining 5 values of last column + for (int mm = rr - 2; mm <= rr + 2; mm++) { + hfnbrave += fabsf(cfablur[mm * W + cc + 2]); + } + +#else + + for (int mm = rr - 2; mm <= rr + 2; mm++) { + int mmmWpnn = mm * W + cc - 2; + + for (int nn = cc - 2; nn <= cc + 2; nn++, mmmWpnn++) { hfnbrave += fabsf(cfablur[mmmWpnn]); } } - if (pixdev * ((bottom - top + 1) * (right - left + 1) - 1) > varthresh * hfnbrave) { +#endif + + if (pixdev > varthresh * hfnbrave) { // mark the pixel as "bad" bpMap.set(cc, rr); counter++; diff --git a/rtengine/sleefsseavx.c b/rtengine/sleefsseavx.c index 6b83c859a..c2d1278f4 100644 --- a/rtengine/sleefsseavx.c +++ b/rtengine/sleefsseavx.c @@ -1322,5 +1322,11 @@ static inline void vswap( vmask condition, vfloat &a, vfloat &b) { b = temp; } +static inline float vhadd( vfloat a ) +{ + a += _mm_movehl_ps(a, a); + return _mm_cvtss_f32(_mm_add_ss(a, _mm_shuffle_ps(a, a, 1))); +} + #endif // __SSE2__ #endif // SLEEFSSEAVX From ad4e96b9ca4ef672b9c162ec45108e6e4c307dc2 Mon Sep 17 00:00:00 2001 From: heckflosse Date: Fri, 2 Oct 2015 23:00:10 +0200 Subject: [PATCH 03/16] Fixed comments and added some #ifdef _OPENMP --- rtengine/rawimagesource.cc | 228 ++++++++++++++++++------------------- 1 file changed, 111 insertions(+), 117 deletions(-) diff --git a/rtengine/rawimagesource.cc b/rtengine/rawimagesource.cc index 6197dba17..2421b33cf 100644 --- a/rtengine/rawimagesource.cc +++ b/rtengine/rawimagesource.cc @@ -31,7 +31,6 @@ #include "curves.h" #include "dfmanager.h" #include "ffmanager.h" -#include "slicer.h" #include "../rtgui/options.h" #include "dcp.h" #include "rt_math.h" @@ -551,7 +550,9 @@ int RawImageSource::interpolateBadPixelsBayer( PixelsMap &bitmapBads ) { static const float eps = 1.f; int counter = 0; +#ifdef _OPENMP #pragma omp parallel for reduction(+:counter) schedule(dynamic,16) +#endif for( int row = 2; row < H - 2; row++ ) { for(int col = 2; col < W - 2; col++ ) { @@ -666,7 +667,9 @@ int RawImageSource::interpolateBadPixelsNColours( PixelsMap &bitmapBads, const i { static const float eps = 1.f; int counter = 0; +#ifdef _OPENMP #pragma omp parallel for reduction(+:counter) schedule(dynamic,16) +#endif for( int row = 2; row < H - 2; row++ ) { for(int col = 2; col < W - 2; col++ ) { @@ -718,7 +721,7 @@ int RawImageSource::interpolateBadPixelsNColours( PixelsMap &bitmapBads, const i } } - if (LIKELY(norm[0] > 0.f)) { // This means, we found at least one pair of valid pixels in the steps above, likelyhood of this case is about 99.999% + if (LIKELY(norm[0] > 0.f)) { // This means, we found at least one pair of valid pixels in the steps above, likelihood of this case is about 99.999% for(int c = 0; c < colours; c++) { rawData[row][col * colours + c] = wtdsum[c] / (2.f * norm[c]); //gradient weighted average, Factor of 2.f is an optimization to avoid multiplications in former steps } @@ -760,13 +763,15 @@ int RawImageSource::interpolateBadPixelsNColours( PixelsMap &bitmapBads, const i return counter; // Number of interpolated pixels. } /* interpolateBadPixelsXtrans: correct raw pixels looking at the bitmap - * takes into consideration if there are multiple bad pixels in the neighborhood + * takes into consideration if there are multiple bad pixels in the neighbourhood */ int RawImageSource::interpolateBadPixelsXtrans( PixelsMap &bitmapBads ) { static const float eps = 1.f; int counter = 0; +#ifdef _OPENMP #pragma omp parallel for reduction(+:counter) schedule(dynamic,16) +#endif for( int row = 2; row < H - 2; row++ ) { for(int col = 2; col < W - 2; col++ ) { @@ -788,10 +793,10 @@ int RawImageSource::interpolateBadPixelsXtrans( PixelsMap &bitmapBads ) if(pixelColor == 1) { // green channel. A green pixel can either be a solitary green pixel or a member of a 2x2 square of green pixels if(ri->XTRANSFC(row, col - 1) == ri->XTRANSFC(row, col + 1)) { - // If left and right neighbour have same color, then this is a solitary green pixel - // For these the following pixels will be used for interpolation. Pixel to be interpolated is in center and marked with a P. + // If left and right neighbour have same colour, then this is a solitary green pixel + // For these the following pixels will be used for interpolation. Pixel to be interpolated is in centre and marked with a P. // Pairs of pixels used in this step are numbered. A pair will be used if none of the pixels of the pair is marked bad - // 0 means, the pixel has a different color and will not be used + // 0 means, the pixel has a different colour and will not be used // 0 1 0 2 0 // 3 5 0 6 4 // 0 0 P 0 0 @@ -830,12 +835,14 @@ int RawImageSource::interpolateBadPixelsXtrans( PixelsMap &bitmapBads ) // this is a member of a 2x2 square of green pixels // For these the following pixels will be used for interpolation. Pixel to be interpolated is at position P in the example. // Pairs of pixels used in this step are numbered. A pair will be used if none of the pixels of the pair is marked bad - // 0 means, the pixel has a different color and will not be used + // 0 means, the pixel has a different colour and will not be used // 1 0 0 3 // 0 P 2 0 // 0 2 1 0 // 3 0 0 0 - int offset1 = ri->XTRANSFC(row - 1, col - 1) == ri->XTRANSFC(row + 1, col + 1) ? 1 : -1; // pixels marked 1 in above example. Distance to P is sqrt(2) => weighting is 0.70710678f + + // pixels marked 1 in above example. Distance to P is sqrt(2) => weighting is 0.70710678f + int offset1 = ri->XTRANSFC(row - 1, col - 1) == ri->XTRANSFC(row + 1, col + 1) ? 1 : -1; if( !(bitmapBads.get(col - offset1, row - 1) || bitmapBads.get(col + offset1, row + 1))) { float dirwt = 0.70710678f / ( fabsf( rawData[row - 1][col - offset1] - rawData[row + 1][col + offset1]) + eps); @@ -843,6 +850,7 @@ int RawImageSource::interpolateBadPixelsXtrans( PixelsMap &bitmapBads ) norm += dirwt; } + // pixels marked 2 in above example. Distance to P is 1 => weighting is 1.f int offsety = (ri->XTRANSFC(row - 1, col) != 1 ? 1 : -1); int offsetx = offset1 * offsety; @@ -866,9 +874,9 @@ int RawImageSource::interpolateBadPixelsXtrans( PixelsMap &bitmapBads ) } } else { // red and blue channel. - // Each red or blue pixel has exactly one neigbour of same color in distance 2 and four neighbours of same color which can be reached by a move of a knight in chess. + // Each red or blue pixel has exactly one neighbour of same colour in distance 2 and four neighbours of same colour which can be reached by a move of a knight in chess. // For the distance 2 pixel (marked with an X) we generate a virtual counterpart (marked with a V) - // For red and blue channel following pixels will be used for interpolation. Pixel to be interpolated is in center and marked with a P. + // For red and blue channel following pixels will be used for interpolation. Pixel to be interpolated is in centre and marked with a P. // Pairs of pixels used in this step are numbered except for distance 2 pixels which are marked X and V. A pair will be used if none of the pixels of the pair is marked bad // 0 1 0 0 0 0 0 X 0 0 remaining cases are symmetric // 0 0 0 0 2 1 0 0 0 2 @@ -876,12 +884,12 @@ int RawImageSource::interpolateBadPixelsXtrans( PixelsMap &bitmapBads ) // 0 0 0 0 1 0 0 0 0 0 // 0 2 0 0 0 0 2 V 1 0 - // Find two knight moves landing on a pixel of same color as the pixel to be interpolated. + // Find two knight moves landing on a pixel of same colour as the pixel to be interpolated. // If we look at first and last row of 5x5 square, we will find exactly two knight pixels. - // Additionally we know that the column of this pixel has 1 or -1 horizontal distance to the center pixel + // Additionally we know that the column of this pixel has 1 or -1 horizontal distance to the centre pixel // When we find a knight pixel, we get its counterpart, which has distance (+-3,+-3), where the signs of distance depend on the corner of the found knight pixel. // These pixels are marked 1 or 2 in above examples. Distance to P is sqrt(5) => weighting is 0.44721359f - // The following loop simply scans the four possible places. To keep things simple, it doesn't stop after finding two knight pixels, because it will not find more than two + // The following loop simply scans the four possible places. To keep things simple, it does not stop after finding two knight pixels, because it will not find more than two for(int d1 = -2, offsety = 3; d1 <= 2; d1 += 4, offsety -= 6) { for(int d2 = -1, offsetx = 3; d2 < 1; d2 += 2, offsetx -= 6) { if(ri->XTRANSFC(row + d1, col + d2) == pixelColor) { @@ -927,10 +935,8 @@ int RawImageSource::interpolateBadPixelsXtrans( PixelsMap &bitmapBads ) norm += dirwt; } - if (LIKELY(norm > 0.f)) { // This means, we found at least one pair of valid pixels in the steps above, likelyhood of this case is about 99.999% + if (LIKELY(norm > 0.f)) { // This means, we found at least one pair of valid pixels in the steps above, likelihood of this case is about 99.999% rawData[row][col] = wtdsum / (2.f * norm); //gradient weighted average, Factor of 2.f is an optimization to avoid multiplications in former steps -//#pragma omp critical -// printf("%s Pixel at (col/row) : (%4d/%4d) : Original : %f, interpolated: %f\n",pixelColor == 0 ? "Red " : pixelColor==1 ? "Green" : "Blue ", col-7,row-7,oldval, rawData[row][col]); counter++; } } @@ -941,9 +947,9 @@ int RawImageSource::interpolateBadPixelsXtrans( PixelsMap &bitmapBads ) //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% /* Search for hot or dead pixels in the image and update the map - * For each pixel compare its value to the average of similar color surrounding + * For each pixel compare its value to the average of similar colour surrounding * (Taken from Emil Martinec idea) - * (Optimized by Ingo Weyrich 2013) + * (Optimized by Ingo Weyrich 2013 and 2015) */ int RawImageSource::findHotDeadPixels( PixelsMap &bpMap, float thresh, bool findHotPixels, bool findDeadPixels ) { @@ -957,15 +963,18 @@ int RawImageSource::findHotDeadPixels( PixelsMap &bpMap, float thresh, bool find // counter for dead or hot pixels int counter = 0; +#ifdef _OPENMP #pragma omp parallel +#endif { +#ifdef _OPENMP #pragma omp for schedule(dynamic,16) nowait +#endif for (int i = 2; i < H - 2; i++) { - float p[9], temp; + float p[9], temp; // we need this for med3x3 macro for (int j = 2; j < W - 2; j++) { - med3x3(rawData[i - 2][j - 2], rawData[i - 2][j], rawData[i - 2][j + 2], rawData[i][j - 2], rawData[i][j], rawData[i][j + 2], rawData[i + 2][j - 2], rawData[i + 2][j], rawData[i + 2][j + 2], temp); @@ -975,7 +984,9 @@ int RawImageSource::findHotDeadPixels( PixelsMap &bpMap, float thresh, bool find // process borders. Former version calculated the median using mirrored border which does not make sense because the original pixel loses weight // Setting the difference between pixel and median for border pixels to zero should do the job not worse then former version +#ifdef _OPENMP #pragma omp critical +#endif { for(int i = 0; i < 2; i++) { for(int j = 0; j < W; j++) { @@ -999,8 +1010,11 @@ int RawImageSource::findHotDeadPixels( PixelsMap &bpMap, float thresh, bool find } } } - #pragma omp barrier - #pragma omp for reduction(+:counter) schedule (dynamic,16) +#ifdef _OPENMP + #pragma omp barrier // barrier because of nowait clause above + + #pragma omp for reduction(+:counter) schedule(dynamic,16) +#endif //cfa pixel heat/death evaluation for (int rr = 2; rr < H - 2; rr++) { @@ -1025,7 +1039,6 @@ int RawImageSource::findHotDeadPixels( PixelsMap &bpMap, float thresh, bool find pixdev = fabsf(pixdev); float hfnbrave = -pixdev; - #ifdef __SSE2__ // sum up 5*4 = 20 values using SSE vfloat sum = vabsf(LVFU(cfablur[(rr - 2) * W + cc - 2])) + vabsf(LVFU(cfablur[(rr - 1) * W + cc - 2])); @@ -1043,10 +1056,8 @@ int RawImageSource::findHotDeadPixels( PixelsMap &bpMap, float thresh, bool find #else for (int mm = rr - 2; mm <= rr + 2; mm++) { - int mmmWpnn = mm * W + cc - 2; - - for (int nn = cc - 2; nn <= cc + 2; nn++, mmmWpnn++) { - hfnbrave += fabsf(cfablur[mmmWpnn]); + for (int nn = cc - 2; nn <= cc + 2; nn++) { + hfnbrave += fabsf(cfablur[mm * W + nn]); } } @@ -1563,7 +1574,10 @@ void RawImageSource::preprocess (const RAWParams &raw, const LensProfParams &le int totBP = 0; // Hold count of bad pixels to correct if(ri->zeroIsBad()) { // mark all pixels with value zero as bad, has to be called before FF and DF. dcraw sets this flag only for some cameras (mainly Panasonic and Leica) +#ifdef _OPENMP #pragma omp parallel for reduction(+:totBP) +#endif + for(int i = 0; i < H; i++) for(int j = 0; j < W; j++) { if(ri->data[i][j] == 0.f) { @@ -1638,7 +1652,9 @@ void RawImageSource::preprocess (const RAWParams &raw, const LensProfParams &le if (pLCPProf && idata->getFocalLen() > 0.f) { LCPMapper map(pLCPProf, idata->getFocalLen(), idata->getFocalLen35mm(), idata->getFocusDist(), idata->getFNumber(), true, false, W, H, coarse, -1); +#ifdef _OPENMP #pragma omp parallel for +#endif for (int y = 0; y < H; y++) { for (int x = 0; x < W; x++) { @@ -1672,7 +1688,9 @@ void RawImageSource::preprocess (const RAWParams &raw, const LensProfParams &le int ng1 = 0, ng2 = 0, i = 0; double avgg1 = 0., avgg2 = 0.; +#ifdef _OPENMP #pragma omp parallel for default(shared) private(i) reduction(+: ng1, ng2, avgg1, avgg2) +#endif for (i = border; i < H - border; i++) for (int j = border; j < W - border; j++) @@ -1689,7 +1707,9 @@ void RawImageSource::preprocess (const RAWParams &raw, const LensProfParams &le double corrg1 = ((double)avgg1 / ng1 + (double)avgg2 / ng2) / 2.0 / ((double)avgg1 / ng1); double corrg2 = ((double)avgg1 / ng1 + (double)avgg2 / ng2) / 2.0 / ((double)avgg2 / ng2); +#ifdef _OPENMP #pragma omp parallel for default(shared) +#endif for (int i = border; i < H - border; i++) for (int j = border; j < W - border; j++) @@ -1909,10 +1929,14 @@ void RawImageSource::processFlatField(const RAWParams &raw, RawImage *riFlatFile float maxval = 0.f; int c = FC(m, n); int c4 = ( c == 1 && !(m & 1) ) ? 3 : c; +#ifdef _OPENMP #pragma omp parallel +#endif { float maxvalthr = 0.f; +#ifdef _OPENMP #pragma omp for +#endif for (int row = 0; row < H - m; row += 2) { for (int col = 0; col < W - n; col += 2) { @@ -1924,7 +1948,9 @@ void RawImageSource::processFlatField(const RAWParams &raw, RawImage *riFlatFile } } +#ifdef _OPENMP #pragma omp critical +#endif { if(maxvalthr > maxval) { @@ -1954,11 +1980,15 @@ void RawImageSource::processFlatField(const RAWParams &raw, RawImage *riFlatFile for (int m = 0; m < 2; m++) for (int n = 0; n < 2; n++) { +#ifdef _OPENMP #pragma omp parallel +#endif { int c = FC(m, n); int c4 = ( c == 1 && !(m & 1) ) ? 3 : c; +#ifdef _OPENMP #pragma omp for +#endif for (int row = 0; row < H - m; row += 2) { @@ -1994,10 +2024,14 @@ void RawImageSource::processFlatField(const RAWParams &raw, RawImage *riFlatFile int clipControlGui = 0; float maxval = 0.f; // xtrans files have only one black level actually, so we can simplify the code a bit +#ifdef _OPENMP #pragma omp parallel +#endif { float maxvalthr = 0.f; +#ifdef _OPENMP #pragma omp for schedule(dynamic,16) nowait +#endif for (int row = 0; row < H; row++) { for (int col = 0; col < W; col++) { @@ -2009,7 +2043,9 @@ void RawImageSource::processFlatField(const RAWParams &raw, RawImage *riFlatFile } } +#ifdef _OPENMP #pragma omp critical +#endif { if(maxvalthr > maxval) { maxval = maxvalthr; @@ -2031,7 +2067,9 @@ void RawImageSource::processFlatField(const RAWParams &raw, RawImage *riFlatFile refcolor[c] *= limitFactor; } +#ifdef _OPENMP #pragma omp parallel for +#endif for (int row = 0; row < H; row++) { for (int col = 0; col < W; col++) { @@ -2054,7 +2092,9 @@ void RawImageSource::processFlatField(const RAWParams &raw, RawImage *riFlatFile if(ri->getSensorType() == ST_BAYER) { for (int m = 0; m < 2; m++) for (int n = 0; n < 2; n++) { +#ifdef _OPENMP #pragma omp parallel for +#endif for (int row = 0; row < H - m; row += 2) { int c = FC(row, 0); @@ -2068,7 +2108,9 @@ void RawImageSource::processFlatField(const RAWParams &raw, RawImage *riFlatFile } } } else if(ri->getSensorType() == ST_FUJI_XTRANS) { +#ifdef _OPENMP #pragma omp parallel for +#endif for (int row = 0; row < H; row++) { for (int col = 0; col < W; col++) { @@ -2177,9 +2219,13 @@ SSEFUNCTION void RawImageSource::cfaboxblur(RawImage *riFlatFile, float* cfablur cfatmp = (float (*)) calloc (H * W, sizeof * cfatmp); // const float hotdeadthresh = 0.5; +#ifdef _OPENMP #pragma omp parallel +#endif { +#ifdef _OPENMP #pragma omp for +#endif for (int i = 0; i < H; i++) { int iprev, inext, jprev, jnext; @@ -2228,7 +2274,9 @@ SSEFUNCTION void RawImageSource::cfaboxblur(RawImage *riFlatFile, float* cfablur //box blur cfa image; box size = BS //horizontal blur +#ifdef _OPENMP #pragma omp for +#endif for (int row = 0; row < H; row++) { int len = boxW / 2 + 1; @@ -2267,7 +2315,9 @@ SSEFUNCTION void RawImageSource::cfaboxblur(RawImage *riFlatFile, float* cfablur __m128 onev = _mm_set1_ps( 1.0f ); __m128 temp1v, temp2v, lenv, lenp1v, lenm1v; int row; +#ifdef _OPENMP #pragma omp for +#endif for (int col = 0; col < W - 3; col += 4) { lenv = leninitv; @@ -2356,7 +2406,9 @@ SSEFUNCTION void RawImageSource::cfaboxblur(RawImage *riFlatFile, float* cfablur } #else +#ifdef _OPENMP #pragma omp for +#endif for (int col = 0; col < W; col++) { int len = boxH / 2 + 1; @@ -2438,12 +2490,15 @@ void RawImageSource::scaleColors(int winx, int winy, int winw, int winh, const R // scale image colors if( ri->getSensorType() == ST_BAYER) { +#ifdef _OPENMP #pragma omp parallel +#endif { float tmpchmax[3]; tmpchmax[0] = tmpchmax[1] = tmpchmax[2] = 0.0f; - +#ifdef _OPENMP #pragma omp for nowait +#endif for (int row = winy; row < winy + winh; row ++) { @@ -2458,7 +2513,9 @@ void RawImageSource::scaleColors(int winx, int winy, int winw, int winh, const R } } +#ifdef _OPENMP #pragma omp critical +#endif { chmax[0] = max(tmpchmax[0], chmax[0]); chmax[1] = max(tmpchmax[1], chmax[1]); @@ -2466,11 +2523,14 @@ void RawImageSource::scaleColors(int winx, int winy, int winw, int winh, const R } } } else if ( ri->get_colors() == 1 ) { +#ifdef _OPENMP #pragma omp parallel +#endif { float tmpchmax = 0.0f; - +#ifdef _OPENMP #pragma omp for nowait +#endif for (int row = winy; row < winy + winh; row ++) { @@ -2483,18 +2543,23 @@ void RawImageSource::scaleColors(int winx, int winy, int winw, int winh, const R } } +#ifdef _OPENMP #pragma omp critical +#endif { chmax[0] = chmax[1] = chmax[2] = chmax[3] = max(tmpchmax, chmax[0]); } } } else if(ri->getSensorType() == ST_FUJI_XTRANS) { +#ifdef _OPENMP #pragma omp parallel +#endif { float tmpchmax[3]; tmpchmax[0] = tmpchmax[1] = tmpchmax[2] = 0.0f; - +#ifdef _OPENMP #pragma omp for nowait +#endif for (int row = winy; row < winy + winh; row ++) { @@ -2509,7 +2574,9 @@ void RawImageSource::scaleColors(int winx, int winy, int winw, int winh, const R } } +#ifdef _OPENMP #pragma omp critical +#endif { chmax[0] = max(tmpchmax[0], chmax[0]); chmax[1] = max(tmpchmax[1], chmax[1]); @@ -2517,12 +2584,15 @@ void RawImageSource::scaleColors(int winx, int winy, int winw, int winh, const R } } } else { +#ifdef _OPENMP #pragma omp parallel +#endif { float tmpchmax[3]; tmpchmax[0] = tmpchmax[1] = tmpchmax[2] = 0.0f; - +#ifdef _OPENMP #pragma omp for nowait +#endif for (int row = winy; row < winy + winh; row ++) { @@ -2537,7 +2607,9 @@ void RawImageSource::scaleColors(int winx, int winy, int winw, int winh, const R } } +#ifdef _OPENMP #pragma omp critical +#endif { chmax[0] = max(tmpchmax[0], chmax[0]); chmax[1] = max(tmpchmax[1], chmax[1]); @@ -2850,8 +2922,9 @@ void RawImageSource::colorSpaceConversion_ (Imagefloat* im, ColorManagementParam mat[i][j] += work[i][k] * camMatrix[k][j]; // rgb_xyz * imatrices.xyz_cam } - +#ifdef _OPENMP #pragma omp parallel for +#endif for (int i = 0; i < im->height; i++) for (int j = 0; j < im->width; j++) { @@ -3174,91 +3247,6 @@ void RawImageSource::colorSpaceConversion_ (Imagefloat* im, ColorManagementParam // printf ("ICM TIME: %d usec\n", t3.etime(t1)); } -//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -// Converts raw image including ICC input profile to working space - 16bit int version -/*void RawImageSource::colorSpaceConversion16 (Image16* im, ColorManagementParams cmp, cmsHPROFILE embedded, cmsHPROFILE camprofile, double camMatrix[3][3], std::string camName) { - cmsHPROFILE in; - DCPProfile *dcpProf; - - if (!findInputProfile(cmp.input, embedded, camName, &dcpProf, in)) return; - - if (dcpProf!=NULL) { - dcpProf->Apply(im, (DCPLightType)cmp.preferredProfile, cmp.working, cmp.toneCurve); - } else { - if (in==NULL) { - // Take camprofile from DCRAW - // in this case we avoid using the slllllooooooowwww lcms - TMatrix work = iccStore->workingSpaceInverseMatrix (cmp.working); - double mat[3][3] = {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}}; - for (int i=0; i<3; i++) - for (int j=0; j<3; j++) - for (int k=0; k<3; k++) - mat[i][j] += work[i][k] * camMatrix[k][j]; // rgb_xyz * imatrices.xyz_cam - -#pragma omp parallel for - for (int i=0; iheight; i++) - for (int j=0; jwidth; j++) { - - float newr = mat[0][0]*im->r(i,j) + mat[0][1]*im->g(i,j) + mat[0][2]*im->b(i,j); - float newg = mat[1][0]*im->r(i,j) + mat[1][1]*im->g(i,j) + mat[1][2]*im->b(i,j); - float newb = mat[2][0]*im->r(i,j) + mat[2][1]*im->g(i,j) + mat[2][2]*im->b(i,j); - - im->r(i,j) = CLIP((int)newr); - im->g(i,j) = CLIP((int)newg); - im->b(i,j) = CLIP((int)newb); - } - } else { - // Gamma preprocessing - float gammaFac, lineFac, lineSum; - getProfilePreprocParams(in, gammaFac, lineFac, lineSum); - - if (gammaFac>0) { -#pragma omp parallel for - for ( int h = 0; h < im->height; ++h ) - for ( int w = 0; w < im->width; ++w ) { - im->r(h,w)= (int) (pow ((double)(im->r(h,w) / 65535.0), (double)gammaFac) * 65535.0); - im->g(h,w)= (int) (pow ((double)(im->g(h,w) / 65535.0), (double)gammaFac) * 65535.0); - im->b(h,w)= (int) (pow ((double)(im->b(h,w) / 65535.0), (double)gammaFac) * 65535.0); - } - } - - cmsHPROFILE out = iccStore->workingSpace (cmp.working); - //out = iccStore->workingSpaceGamma (wProfile); - lcmsMutex->lock (); - cmsHTRANSFORM hTransform = cmsCreateTransform (in, TYPE_RGB_16, out, TYPE_RGB_16, settings->colorimetricIntent, - cmsFLAGS_NOCACHE); // NOCACHE is important for thread safety - lcmsMutex->unlock (); - - if (hTransform) { - im->ExecCMSTransform(hTransform); - - // There might be Nikon postprocessings - if (lineSum>0) { -#pragma omp parallel for - for ( int h = 0; h < im->height; ++h ) - for ( int w = 0; w < im->width; ++w ) { - im->r(h,w) *= im->r(h,w) * lineFac / 65535.0 + lineSum; - im->g(h,w) *= im->g(h,w) * lineFac / 65535.0 + lineSum; - im->b(h,w) *= im->b(h,w) * lineFac / 65535.0 + lineSum; - } - } - } - else { - lcmsMutex->lock (); - hTransform = cmsCreateTransform (camprofile, TYPE_RGB_16, out, TYPE_RGB_16, - settings->colorimetricIntent, cmsFLAGS_NOCACHE); - lcmsMutex->unlock (); - - im->ExecCMSTransform(hTransform); - } - - cmsDeleteTransform(hTransform); - } - } - //t3.set (); - //printf ("ICM TIME: %d\n", t3.etime(t1)); -}*/ // Determine RAW input and output profiles. Returns TRUE on success bool RawImageSource::findInputProfile(Glib::ustring inProfile, cmsHPROFILE embedded, std::string camName, DCPProfile **dcpProf, cmsHPROFILE& in) @@ -3557,11 +3545,15 @@ void RawImageSource::getAutoExpHistogram (LUTu & histogram, int& histcompr) histogram(65536 >> histcompr); histogram.clear(); +#ifdef _OPENMP #pragma omp parallel +#endif { LUTu tmphistogram(65536 >> histcompr); tmphistogram.clear(); +#ifdef _OPENMP #pragma omp for nowait +#endif for (int i = border; i < H - border; i++) { int start, end; @@ -3600,7 +3592,9 @@ void RawImageSource::getAutoExpHistogram (LUTu & histogram, int& histcompr) } } +#ifdef _OPENMP #pragma omp critical +#endif { for(int i = 0; i < (65536 >> histcompr); i++) { histogram[i] += tmphistogram[i]; From 643a9e7d54d6c3266aa4d20fbf9e791c3038753f Mon Sep 17 00:00:00 2001 From: heckflosse Date: Fri, 2 Oct 2015 23:01:56 +0200 Subject: [PATCH 04/16] added a file I forgot at last commit --- rtengine/stdimagesource.h | 1 - 1 file changed, 1 deletion(-) diff --git a/rtengine/stdimagesource.h b/rtengine/stdimagesource.h index 3848fe5e9..6017fae74 100644 --- a/rtengine/stdimagesource.h +++ b/rtengine/stdimagesource.h @@ -92,7 +92,6 @@ public: void convertColorSpace(Imagefloat* image, ColorManagementParams cmp, ColorTemp &wb);// RAWParams raw will not be used for non-raw files (see imagesource.h) static void colorSpaceConversion (Imagefloat* im, ColorManagementParams cmp, cmsHPROFILE embedded, IIOSampleFormat sampleFormat); - //static void colorSpaceConversion16 (Image16* im, ColorManagementParams cmp, cmsHPROFILE embedded); static inline double intpow (double a, int b) { From d05754fbf8bbb3943a72beea4ee4dcbbcb96cb54 Mon Sep 17 00:00:00 2001 From: heckflosse Date: Fri, 2 Oct 2015 23:03:32 +0200 Subject: [PATCH 05/16] Reduced updates to progress bar for amaze demosaic --- rtengine/amaze_demosaic_RT.cc | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/rtengine/amaze_demosaic_RT.cc b/rtengine/amaze_demosaic_RT.cc index 10a3687e6..289bfa46e 100644 --- a/rtengine/amaze_demosaic_RT.cc +++ b/rtengine/amaze_demosaic_RT.cc @@ -31,13 +31,14 @@ #include "procparams.h" #include "sleef.c" #include "opthelper.h" +#include "StopWatch.h" namespace rtengine { SSEFUNCTION void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, int winh) { - +StopWatch Stop1("amaze_demosaic_RT"); #define HCLIP(x) x //is this still necessary??? //min(clip_pt,x) @@ -1615,10 +1616,10 @@ SSEFUNCTION void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, if(plistener) { progresscounter++; - if(progresscounter % 4 == 0) { + if(progresscounter % 16 == 0) { #pragma omp critical { - progress += (double)4 * ((TS - 32) * (TS - 32)) / (height * width); + progress += (double)16 * ((TS - 32) * (TS - 32)) / (height * width); if (progress > 1.0) { From 929c0ad1b3aecf18fce6245e0b75809990d8c364 Mon Sep 17 00:00:00 2001 From: heckflosse Date: Sat, 3 Oct 2015 00:17:18 +0200 Subject: [PATCH 06/16] Removed Stopwatches --- rtengine/amaze_demosaic_RT.cc | 2 -- rtengine/rawimagesource.cc | 2 -- 2 files changed, 4 deletions(-) diff --git a/rtengine/amaze_demosaic_RT.cc b/rtengine/amaze_demosaic_RT.cc index 289bfa46e..dc2d4316e 100644 --- a/rtengine/amaze_demosaic_RT.cc +++ b/rtengine/amaze_demosaic_RT.cc @@ -31,14 +31,12 @@ #include "procparams.h" #include "sleef.c" #include "opthelper.h" -#include "StopWatch.h" namespace rtengine { SSEFUNCTION void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, int winh) { -StopWatch Stop1("amaze_demosaic_RT"); #define HCLIP(x) x //is this still necessary??? //min(clip_pt,x) diff --git a/rtengine/rawimagesource.cc b/rtengine/rawimagesource.cc index 2421b33cf..6381c45b6 100644 --- a/rtengine/rawimagesource.cc +++ b/rtengine/rawimagesource.cc @@ -39,7 +39,6 @@ #include #endif #include "opthelper.h" -#include "StopWatch.h" namespace rtengine { @@ -953,7 +952,6 @@ int RawImageSource::interpolateBadPixelsXtrans( PixelsMap &bitmapBads ) */ int RawImageSource::findHotDeadPixels( PixelsMap &bpMap, float thresh, bool findHotPixels, bool findDeadPixels ) { - StopWatch Stop1("findHotDeadPixels"); float varthresh = (20.0 * (thresh / 100.0) + 1.0 ) / 24.f; // allocate temporary buffer From 20839d3a4a0de88a217a3709a089749268bc6bfc Mon Sep 17 00:00:00 2001 From: heckflosse Date: Sat, 3 Oct 2015 13:33:52 +0200 Subject: [PATCH 07/16] Color propagation highlight reconstruction causes crash on some files, fixes #2928 --- rtengine/hilite_recon.cc | 7 ++++--- rtengine/rawimagesource.cc | 2 +- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/rtengine/hilite_recon.cc b/rtengine/hilite_recon.cc index 411bb0823..bc5e973e2 100644 --- a/rtengine/hilite_recon.cc +++ b/rtengine/hilite_recon.cc @@ -847,21 +847,22 @@ void RawImageSource :: HLRecovery_inpaint (float** red, float** green, float** b green[i][j] = clipfix[1] * factor; blue[i][j] = clipfix[2] * factor; } else {//some channels clipped + const float eps = 0.0001f; int notclipped[3] = {pixel[0] < max_f[0] ? 1 : 0, pixel[1] < max_f[1] ? 1 : 0, pixel[2] < max_f[2] ? 1 : 0}; if (notclipped[0] == 0) { //red clipped red[i][j] = max(red[i][j], (clipfix[0] * ((notclipped[1] * pixel[1] + notclipped[2] * pixel[2]) / - (notclipped[1] * clipfix[1] + notclipped[2] * clipfix[2])))); + (notclipped[1] * clipfix[1] + notclipped[2] * clipfix[2] + eps)))); } if (notclipped[1] == 0) { //green clipped green[i][j] = max(green[i][j], (clipfix[1] * ((notclipped[2] * pixel[2] + notclipped[0] * pixel[0]) / - (notclipped[2] * clipfix[2] + notclipped[0] * clipfix[0])))); + (notclipped[2] * clipfix[2] + notclipped[0] * clipfix[0] + eps)))); } if (notclipped[2] == 0) { //blue clipped blue[i][j] = max(blue[i][j], (clipfix[2] * ((notclipped[0] * pixel[0] + notclipped[1] * pixel[1]) / - (notclipped[0] * clipfix[0] + notclipped[1] * clipfix[1])))); + (notclipped[0] * clipfix[0] + notclipped[1] * clipfix[1] + eps)))); } } diff --git a/rtengine/rawimagesource.cc b/rtengine/rawimagesource.cc index 6381c45b6..3194e2f32 100644 --- a/rtengine/rawimagesource.cc +++ b/rtengine/rawimagesource.cc @@ -950,7 +950,7 @@ int RawImageSource::interpolateBadPixelsXtrans( PixelsMap &bitmapBads ) * (Taken from Emil Martinec idea) * (Optimized by Ingo Weyrich 2013 and 2015) */ -int RawImageSource::findHotDeadPixels( PixelsMap &bpMap, float thresh, bool findHotPixels, bool findDeadPixels ) +SSEFUNCTION int RawImageSource::findHotDeadPixels( PixelsMap &bpMap, float thresh, bool findHotPixels, bool findDeadPixels ) { float varthresh = (20.0 * (thresh / 100.0) + 1.0 ) / 24.f; From 9dc786f0ea236bc793cd77b365fae1fbf9772017 Mon Sep 17 00:00:00 2001 From: heckflosse Date: Tue, 6 Oct 2015 22:29:44 +0200 Subject: [PATCH 08/16] Segfault when curve black point is moved below white point, fixes #2923 --- rtengine/diagonalcurves.cc | 6 +++++- rtengine/rawimagesource.cc | 2 +- 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/rtengine/diagonalcurves.cc b/rtengine/diagonalcurves.cc index f46cdb8cb..f09eb3459 100644 --- a/rtengine/diagonalcurves.cc +++ b/rtengine/diagonalcurves.cc @@ -72,7 +72,11 @@ DiagonalCurve::DiagonalCurve (const std::vector& p, int poly_pn) if(x[0] == 0.f && x[1] == 0.f) // Avoid crash when first two points are at x = 0 (git Issue 2888) - x[1] = 0.5f; + x[1] = 0.01f; + + if(x[0] == 1.f && x[1] == 1.f) + // Avoid crash when first two points are at x = 1 (100 in gui) (git Issue 2923) + x[0] = 0.99f; if (!identity) { if (kind == DCT_Spline && N > 2) { diff --git a/rtengine/rawimagesource.cc b/rtengine/rawimagesource.cc index 3194e2f32..a85e7a7d3 100644 --- a/rtengine/rawimagesource.cc +++ b/rtengine/rawimagesource.cc @@ -983,7 +983,7 @@ SSEFUNCTION int RawImageSource::findHotDeadPixels( PixelsMap &bpMap, float thres // process borders. Former version calculated the median using mirrored border which does not make sense because the original pixel loses weight // Setting the difference between pixel and median for border pixels to zero should do the job not worse then former version #ifdef _OPENMP - #pragma omp critical + #pragma omp single #endif { for(int i = 0; i < 2; i++) { From b1dd9dd59a4895f5d6928406c381072b68a205d5 Mon Sep 17 00:00:00 2001 From: heckflosse Date: Sat, 10 Oct 2015 01:05:08 +0200 Subject: [PATCH 09/16] Speedup and reduced memory usage for Colour Propagation --- rtengine/array2D.h | 37 +- rtengine/hilite_recon.cc | 760 +++++++++++++++++++++++--------------- rtengine/rawimagesource.h | 4 +- 3 files changed, 477 insertions(+), 324 deletions(-) diff --git a/rtengine/array2D.h b/rtengine/array2D.h index 93734a421..f9d92ac36 100644 --- a/rtengine/array2D.h +++ b/rtengine/array2D.h @@ -23,7 +23,7 @@ * Usage: * * array2D name (X-size,Y-size); - * array2D name (X-size,Y-size type ** data); + * array2D name (X-size,Y-size,type ** data); * * creates an array which is valid within the normal C/C++ scope "{ ... }" * @@ -75,7 +75,7 @@ private: T ** ptr; T * data; bool lock; // useful lock to ensure data is not changed anymore. - void ar_realloc(int w, int h) + void ar_realloc(int w, int h, int offset = 0) { if ((ptr) && ((h > y) || (4 * h < y))) { delete[] ptr; @@ -92,14 +92,14 @@ private: } if (data == NULL) { - data = new T[h * w]; + data = new T[h * w + offset]; } x = w; y = h; for (int i = 0; i < h; i++) { - ptr[i] = data + w * i; + ptr[i] = data + offset + w * i; } owner = 1; @@ -184,6 +184,19 @@ public: } } + void free() + { + if ((owner) && (data)) { + delete[] data; + data = NULL; + } + + if (ptr) { + delete [] ptr; + ptr = NULL; + } + } + // use with indices T * operator[](int index) { @@ -207,7 +220,7 @@ public: // useful within init of parent object // or use as resize of 2D array - void operator()(int w, int h, unsigned int flgs = 0) + void operator()(int w, int h, unsigned int flgs = 0, int offset = 0) { flags = flgs; @@ -223,10 +236,10 @@ public: lock = flags & ARRAY2D_LOCK_DATA; - ar_realloc(w, h); + ar_realloc(w, h, offset); if (flags & ARRAY2D_CLEAR_DATA) { - memset(data, 0, w * h * sizeof(T)); + memset(data + offset, 0, w * h * sizeof(T)); } } @@ -298,10 +311,10 @@ private: array2D list[num]; public: - multi_array2D(int x, int y, int flags = 0) + multi_array2D(int x, int y, int flags = 0, int offset = 0) { for (size_t i = 0; i < num; i++) { - list[i](x, y, flags); + list[i](x, y, flags, (i + 1) * offset); } } @@ -312,11 +325,7 @@ public: array2D & operator[](int index) { - if (static_cast(index) >= num) { - printf("index %0u is out of range[0..%0lu]", index, num - 1); - raise( SIGSEGV); - } - + assert(static_cast(index) < num); return list[index]; } }; diff --git a/rtengine/hilite_recon.cc b/rtengine/hilite_recon.cc index bc5e973e2..46c666c80 100644 --- a/rtengine/hilite_recon.cc +++ b/rtengine/hilite_recon.cc @@ -29,25 +29,15 @@ #include "rt_math.h" #include "opthelper.h" -#ifdef _OPENMP -#include -#endif - - #define FOREACHCOLOR for (int c=0; c < ColorCount; c++) -//#include "RGBdefringe.cc" - namespace rtengine { - //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -SSEFUNCTION void RawImageSource::boxblur2(float** src, float** dst, int H, int W, int box ) +SSEFUNCTION void RawImageSource::boxblur2(float** src, float** dst, float** temp, int H, int W, int box ) { - array2D temp(W, H); - //box blur image channel; box size = 2*box+1 //horizontal blur #ifdef _OPENMP @@ -84,10 +74,10 @@ SSEFUNCTION void RawImageSource::boxblur2(float** src, float** dst, int H, int W #endif { float len = box + 1; - __m128 lenv = _mm_set1_ps( len ); - __m128 lenp1v = _mm_set1_ps( len + 1.0f ); - __m128 onev = _mm_set1_ps( 1.0f ); - __m128 tempv, temp2v; + vfloat lenv = F2V( len ); + vfloat lenp1v = F2V( len + 1.0f ); + vfloat onev = F2V( 1.0f ); + vfloat tempv, temp2v; #ifdef _OPENMP #pragma omp for nowait #endif @@ -130,7 +120,9 @@ SSEFUNCTION void RawImageSource::boxblur2(float** src, float** dst, int H, int W } } +#ifdef _OPENMP #pragma omp single +#endif { for (int col = W - (W % 8); col < W - 3; col += 4) { tempv = LVFU(temp[0][col]) / lenv; @@ -160,29 +152,28 @@ SSEFUNCTION void RawImageSource::boxblur2(float** src, float** dst, int H, int W _mm_storeu_ps( &dst[row][col], tempv ); } } - } - } + for (int col = W - (W % 4); col < W; col++) { + int len = box + 1; + dst[0][col] = temp[0][col] / len; - for (int col = W - (W % 4); col < W; col++) { - int len = box + 1; - dst[0][col] = temp[0][col] / len; + for (int i = 1; i <= box; i++) { + dst[0][col] += temp[i][col] / len; + } - for (int i = 1; i <= box; i++) { - dst[0][col] += temp[i][col] / len; - } + for (int row = 1; row <= box; row++) { + dst[row][col] = (dst[(row - 1)][col] * len + temp[(row + box)][col]) / (len + 1); + len ++; + } - for (int row = 1; row <= box; row++) { - dst[row][col] = (dst[(row - 1)][col] * len + temp[(row + box)][col]) / (len + 1); - len ++; - } + for (int row = box + 1; row < H - box; row++) { + dst[row][col] = dst[(row - 1)][col] + (temp[(row + box)][col] - temp[(row - box - 1)][col]) / len; + } - for (int row = box + 1; row < H - box; row++) { - dst[row][col] = dst[(row - 1)][col] + (temp[(row + box)][col] - temp[(row - box - 1)][col]) / len; - } - - for (int row = H - box; row < H; row++) { - dst[row][col] = (dst[(row - 1)][col] * len - temp[(row - box - 1)][col]) / (len - 1); - len --; + for (int row = H - box; row < H; row++) { + dst[row][col] = (dst[(row - 1)][col] * len - temp[(row - box - 1)][col]) / (len - 1); + len --; + } + } } } @@ -219,11 +210,9 @@ SSEFUNCTION void RawImageSource::boxblur2(float** src, float** dst, int H, int W } -void RawImageSource::boxblur_resamp(float **src, float **dst, int H, int W, int box, int samp ) +void RawImageSource::boxblur_resamp(float **src, float **dst, float ** temp, int H, int W, int box, int samp ) { - array2D temp((W / samp) + ((W % samp) == 0 ? 0 : 1), H); - #ifdef _OPENMP #pragma omp parallel #endif @@ -235,7 +224,8 @@ void RawImageSource::boxblur_resamp(float **src, float **dst, int H, int W, int //box blur image channel; box size = 2*box+1 //horizontal blur - for (int row = 0; row < H; row++) { + for (int row = 0; row < H; row++) + { int len = box + 1; tempval = src[row][0] / len; @@ -255,8 +245,10 @@ void RawImageSource::boxblur_resamp(float **src, float **dst, int H, int W, int len ++; } + float oneByLen = 1.f / (float)len; + for (int col = box + 1; col < W - box; col++) { - tempval = tempval + (src[row][col + box] - src[row][col - box - 1]) / len; + tempval = tempval + (src[row][col + box] - src[row][col - box - 1]) * oneByLen; if(col % samp == 0) { temp[row][col / samp] = tempval; @@ -275,63 +267,125 @@ void RawImageSource::boxblur_resamp(float **src, float **dst, int H, int W, int } } + static const int numCols = 8; // process numCols columns at once for better L1 CPU cache usage #ifdef _OPENMP #pragma omp parallel #endif { - float tempval; + float tempvalN[numCols] ALIGNED16; #ifdef _OPENMP - #pragma omp for + #pragma omp for nowait #endif //vertical blur - for (int col = 0; col < W / samp; col++) { + for (int col = 0; col < (W / samp) - (numCols - 1); col += numCols) { int len = box + 1; - tempval = temp[0][col] / len; - for (int i = 1; i <= box; i++) { - tempval += temp[i][col] / len; + for(int n = 0; n < numCols; n++) { + tempvalN[n] = temp[0][col + n] / len; } - dst[0][col] = tempval; + for (int i = 1; i <= box; i++) { + for(int n = 0; n < numCols; n++) { + tempvalN[n] += temp[i][col + n] / len; + } + } + + for(int n = 0; n < numCols; n++) { + dst[0][col + n] = tempvalN[n]; + } for (int row = 1; row <= box; row++) { - tempval = (tempval * len + temp[(row + box)][col]) / (len + 1); + for(int n = 0; n < numCols; n++) { + tempvalN[n] = (tempvalN[n] * len + temp[(row + box)][col + n]) / (len + 1); + } if(row % samp == 0) { - dst[row / samp][col] = tempval; + for(int n = 0; n < numCols; n++) { + dst[row / samp][col + n] = tempvalN[n]; + } } len ++; } for (int row = box + 1; row < H - box; row++) { - tempval = tempval + (temp[(row + box)][col] - temp[(row - box - 1)][col]) / len; + for(int n = 0; n < numCols; n++) { + tempvalN[n] = tempvalN[n] + (temp[(row + box)][col + n] - temp[(row - box - 1)][col + n]) / len; + } if(row % samp == 0) { - dst[row / samp][col] = tempval; + for(int n = 0; n < numCols; n++) { + dst[row / samp][col + n] = tempvalN[n]; + } } } for (int row = H - box; row < H; row++) { - tempval = (tempval * len - temp[(row - box - 1)][col]) / (len - 1); + for(int n = 0; n < numCols; n++) { + tempvalN[n] = (tempvalN[n] * len - temp[(row - box - 1)][col + n]) / (len - 1); + } if(row % samp == 0) { - dst[row / samp][col] = tempval; + for(int n = 0; n < numCols; n++) { + dst[row / samp][col + n] = tempvalN[n]; + } } len --; } } + + // process remaining columns + #pragma omp single + { + float tempval; + + //vertical blur + for (int col = (W / samp) - ((W / samp) % numCols); col < W / samp; col++) { + int len = box + 1; + tempval = temp[0][col] / len; + + for (int i = 1; i <= box; i++) { + tempval += temp[i][col] / len; + } + + dst[0][col] = tempval; + + for (int row = 1; row <= box; row++) { + tempval = (tempval * len + temp[(row + box)][col]) / (len + 1); + + if(row % samp == 0) { + dst[row / samp][col] = tempval; + } + + len ++; + } + + for (int row = box + 1; row < H - box; row++) { + tempval = tempval + (temp[(row + box)][col] - temp[(row - box - 1)][col]) / len; + + if(row % samp == 0) { + dst[row / samp][col] = tempval; + } + } + + for (int row = H - box; row < H; row++) { + tempval = (tempval * len - temp[(row - box - 1)][col]) / (len - 1); + + if(row % samp == 0) { + dst[row / samp][col] = tempval; + } + + len --; + } + } + } } } -//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - void RawImageSource :: HLRecovery_inpaint (float** red, float** green, float** blue) { double progress = 0.0; @@ -344,12 +398,8 @@ void RawImageSource :: HLRecovery_inpaint (float** red, float** green, float** b int height = H; int width = W; - const int range = 2; - const int pitch = 4; - - - int hfh = (height - (height % pitch)) / pitch; - int hfw = (width - (width % pitch)) / pitch; + static const int range = 2; + static const int pitch = 4; static const int numdirs = 4; @@ -360,7 +410,7 @@ void RawImageSource :: HLRecovery_inpaint (float** red, float** green, float** b //%%%%%%%%%%%%%%%%%%%% //for blend algorithm: static const float blendthresh = 1.0; - const int ColorCount = 3; + static const int ColorCount = 3; // Transform matrixes rgb>lab and back static const float trans[2][ColorCount][ColorCount] = { { { 1, 1, 1 }, { 1.7320508, -1.7320508, 0 }, { -1, -1, 2 } }, @@ -370,49 +420,12 @@ void RawImageSource :: HLRecovery_inpaint (float** red, float** green, float** b { { 1, 0.8660254, -0.5 }, { 1, -0.8660254, -0.5 }, { 1, 0, 1 } }, { { 1, 1, 1 }, { 1, -1, 1 }, { 1, 1, -1 } } }; - //%%%%%%%%%%%%%%%%%%%% - - float max_f[3], thresh[3], fixthresh[3], norm[3]; - - //float red1, green1, blue1;//diagnostic -// float chmaxalt[4]={0,0,0,0};//diagnostic - - //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - //halfsize demosaic - - /* - multi_array2D hfsize (hfw+1,hfh+1,ARRAY2D_CLEAR_DATA); - - boxblur_resamp(red,hfsize[0],chmaxalt[0],height,width,range,pitch); - if(plistener){ - progress += 0.05; - plistener->setProgress(progress); - } - boxblur_resamp(green,hfsize[1],chmaxalt[1],height,width,range,pitch); - if(plistener){ - progress += 0.05; - plistener->setProgress(progress); - } - boxblur_resamp(blue,hfsize[2],chmaxalt[2],height,width,range,pitch); - if(plistener){ - progress += 0.05; - plistener->setProgress(progress); - } - - //blur image - //for (int m=0; m<3; m++) - // boxblur2(hfsize[m],hfsizeblur[m],hfh,hfw,3); - - */ - //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + float max_f[3], thresh[3]; for (int c = 0; c < 3; c++) { thresh[c] = chmax[c] * threshpct; - fixthresh[c] = chmax[c] * fixthreshpct; max_f[c] = chmax[c] * maxpct; //min(chmax[0],chmax[1],chmax[2])*maxpct; - norm[c] = 1.0 / (max_f[c] - fixthresh[c]); } float whitept = max(max_f[0], max_f[1], max_f[2]); @@ -420,43 +433,57 @@ void RawImageSource :: HLRecovery_inpaint (float** red, float** green, float** b float medpt = max_f[0] + max_f[1] + max_f[2] - whitept - clippt; float blendpt = blendthresh * clippt; - float camwb[4]; - - for (int c = 0; c < 4; c++) { - camwb[c] = ri->get_cam_mul(c); - } - - multi_array2D channelblur(width, height, ARRAY2D_CLEAR_DATA); - multi_array2D hilite_full(width, height, ARRAY2D_CLEAR_DATA); - - if(plistener) { - progress += 0.05; - plistener->setProgress(progress); - } + multi_array2D channelblur(width, height, 0, 48); + array2D temp(width, height); // allocate temporary buffer // blur RGB channels - boxblur2(red , channelblur[0], height, width, 4); + + boxblur2(red , channelblur[0], temp, height, width, 4); + if(plistener) { + progress += 0.05; + plistener->setProgress(progress); + } + + boxblur2(green, channelblur[1], temp, height, width, 4); + if(plistener) { + progress += 0.05; + plistener->setProgress(progress); + } + + boxblur2(blue , channelblur[2], temp, height, width, 4); + if(plistener) { + progress += 0.05; + plistener->setProgress(progress); + } + + // reduce channel blur to one array +#ifdef _OPENMP + #pragma omp parallel for +#endif + + for(int i = 0; i < height; i++) + for(int j = 0; j < width; j++) { + channelblur[0][i][j] = fabsf(channelblur[0][i][j] - red[i][j]) + fabsf(channelblur[1][i][j] - green[i][j]) + fabsf(channelblur[2][i][j] - blue[i][j]); + } + + for (int c = 1; c < 3; c++) { + channelblur[c].free(); //free up some memory + } if(plistener) { progress += 0.05; plistener->setProgress(progress); } - boxblur2(green, channelblur[1], height, width, 4); - + multi_array2D hilite_full(width, height, ARRAY2D_CLEAR_DATA, 32); if(plistener) { - progress += 0.05; + progress += 0.10; plistener->setProgress(progress); } - boxblur2(blue , channelblur[2], height, width, 4); - if(plistener) { - progress += 0.05; - plistener->setProgress(progress); - } - - float hipass_sum = 0, hipass_norm = 0.00; + double hipass_sum = 0.f; + int hipass_norm = 0; // set up which pixels are clipped or near clipping #ifdef _OPENMP @@ -465,40 +492,34 @@ void RawImageSource :: HLRecovery_inpaint (float** red, float** green, float** b for (int i = 0; i < height; i++) { for (int j = 0; j < width; j++) { - //if one or more channels is highlight but none are blown, add to highlight accumulator - if ((red[i][j] > thresh[0] || green[i][j] > thresh[1] || blue[i][j] > thresh[2]) && (red[i][j] < max_f[0] && green[i][j] < max_f[1] && blue[i][j] < max_f[2])) { - hipass_sum += fabs(channelblur[0][i][j] - red[i][j]) + fabs(channelblur[1][i][j] - green[i][j]) + fabs(channelblur[2][i][j] - blue[i][j]); - hipass_norm++; + hipass_sum += channelblur[0][i][j]; + hipass_norm ++; hilite_full[0][i][j] = red[i][j]; hilite_full[1][i][j] = green[i][j]; hilite_full[2][i][j] = blue[i][j]; - hilite_full[3][i][j] = 1; - hilite_full[4][i][j] = 1; + hilite_full[3][i][j] = 1.f; } - - //if (i%100==0 && j%100==0) - // printf("row=%d col=%d r=%f g=%f b=%f hilite=%f \n",i,j,hilite_full[0][i][j],hilite_full[1][i][j],hilite_full[2][i][j],hilite_full[3][i][j]); } }//end of filling highlight array - hipass_norm += 0.01; - - float hipass_ave = (hipass_sum / hipass_norm); - + float hipass_ave = 2.f * hipass_sum / (hipass_norm + 0.01); if(plistener) { progress += 0.05; plistener->setProgress(progress); } + array2D hilite_full4(width, height); //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% //blur highlight data - boxblur2(hilite_full[4], hilite_full[4], height, width, 1); + boxblur2(hilite_full[3], hilite_full4, temp, height, width, 1); + + temp.free(); // free temporary buffer if(plistener) { progress += 0.05; @@ -511,32 +532,34 @@ void RawImageSource :: HLRecovery_inpaint (float** red, float** green, float** b for (int i = 0; i < height; i++) { for (int j = 0; j < width; j++) { - - float hipass = fabs(channelblur[0][i][j] - red[i][j]) + fabs(channelblur[1][i][j] - green[i][j]) + fabs(channelblur[2][i][j] - blue[i][j]); - - if (hipass > 2 * hipass_ave) { + if (channelblur[0][i][j] > hipass_ave) { //too much variation - hilite_full[0][i][j] = hilite_full[1][i][j] = hilite_full[2][i][j] = hilite_full[3][i][j] = 0; + hilite_full[0][i][j] = hilite_full[1][i][j] = hilite_full[2][i][j] = hilite_full[3][i][j] = 0.f; continue; } - if (hilite_full[4][i][j] > 0.00001 && hilite_full[4][i][j] < 0.95) { + if (hilite_full4[i][j] > 0.00001 && hilite_full4[i][j] < 0.95) { //too near an edge, could risk using CA affected pixels, therefore omit - hilite_full[0][i][j] = hilite_full[1][i][j] = hilite_full[2][i][j] = hilite_full[3][i][j] = 0; + hilite_full[0][i][j] = hilite_full[1][i][j] = hilite_full[2][i][j] = hilite_full[3][i][j] = 0.f; } } } - for (int c = 0; c < 3; c++) { - channelblur[c](1, 1); //free up some memory - } + channelblur[0].free(); //free up some memory + hilite_full4.free(); //free up some memory - multi_array2D hilite(hfw + 1, hfh + 1, ARRAY2D_CLEAR_DATA); + int hfh = (height - (height % pitch)) / pitch; + int hfw = (width - (width % pitch)) / pitch; + + multi_array2D hilite(hfw + 1, hfh + 1, ARRAY2D_CLEAR_DATA, 48); //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% // blur and resample highlight data; range=size of blur, pitch=sample spacing + + array2D temp2((width / pitch) + ((width % pitch) == 0 ? 0 : 1), height); + for (int m = 0; m < 4; m++) { - boxblur_resamp(hilite_full[m], hilite[m], height, width, range, pitch); + boxblur_resamp(hilite_full[m], hilite[m], temp2, height, width, range, pitch); if(plistener) { progress += 0.05; @@ -544,101 +567,235 @@ void RawImageSource :: HLRecovery_inpaint (float** red, float** green, float** b } } - for (int c = 0; c < 5; c++) { - hilite_full[c](1, 1); //free up some memory + temp2.free(); + + for (int c = 0; c < 4; c++) { + hilite_full[c].free(); //free up some memory } - multi_array2D hilite_dir(hfw, hfh, ARRAY2D_CLEAR_DATA); - + multi_array2D hilite_dir(hfw, hfh, ARRAY2D_CLEAR_DATA, 64); + // for faster processing we create two buffers using (height,width) instead of (width,height) + multi_array2D hilite_dir0(hfh, hfw, ARRAY2D_CLEAR_DATA, 64); + multi_array2D hilite_dir4(hfh, hfw, ARRAY2D_CLEAR_DATA, 64); //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - //blur highlights - //for (int m=0; m<4; m++) - // boxblur2(hilite[m],hilite[m],hfh,hfw,4); - - //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% if(plistener) { progress += 0.05; plistener->setProgress(progress); } - LUTf invfn(0x10000); - - for (int i = 0; i < 0x10000; i++) { - invfn[i] = 1.0 / (1.0 + i); - } - //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% //fill gaps in highlight map by directional extension //raster scan from four corners - for (int j = 1; j < hfw - 1; j++) + for (int j = 1; j < hfw - 1; j++) { for (int i = 2; i < hfh - 2; i++) { //from left if (hilite[3][i][j] > 0.01) { - for (int c = 0; c < 4; c++) { - hilite_dir[c][i][j] = hilite[c][i][j] / hilite[3][i][j]; - } + hilite_dir0[3][j][i] = 1.f; } else { - for (int c = 0; c < 4; c++) { - hilite_dir[c][i][j] = 0.1 * ((hilite_dir[0 + c][i - 2][j - 1] + hilite_dir[0 + c][i - 1][j - 1] + hilite_dir[0 + c][i][j - 1] + hilite_dir[0 + c][i + 1][j - 1] + hilite_dir[0 + c][i + 2][j - 1]) / - (hilite_dir[0 + 3][i - 2][j - 1] + hilite_dir[0 + 3][i - 1][j - 1] + hilite_dir[0 + 3][i][j - 1] + hilite_dir[0 + 3][i + 1][j - 1] + hilite_dir[0 + 3][i + 2][j - 1] + 0.00001)); - hilite_dir[4 + c][i][j + 1] += hilite_dir[c][i][j]; - hilite_dir[8 + c][i - 2][j] += hilite_dir[c][i][j]; - hilite_dir[12 + c][i + 2][j] += hilite_dir[c][i][j]; - } + hilite_dir0[3][j][i] = (hilite_dir0[0 + 3][j - 1][i - 2] + hilite_dir0[0 + 3][j - 1][i - 1] + hilite_dir0[0 + 3][j - 1][i] + hilite_dir0[0 + 3][j - 1][i + 1] + hilite_dir0[0 + 3][j - 1][i + 2]) == 0.f ? 0.f : 0.1f; } } + if(hilite[3][2][j] <= 0.01) { + hilite_dir[0 + 3][0][j] = hilite_dir0[3][j][2]; + } + + if(hilite[3][3][j] <= 0.01) { + hilite_dir[0 + 3][1][j] = hilite_dir0[3][j][3]; + } + + if(hilite[3][hfh - 3][j] <= 0.01) { + hilite_dir[4 + 3][hfh - 1][j] = hilite_dir0[3][j][hfh - 3]; + } + + if(hilite[3][hfh - 4][j] <= 0.01) { + hilite_dir[4 + 3][hfh - 2][j] = hilite_dir0[3][j][hfh - 4]; + } + } + + for (int i = 2; i < hfh - 2; i++) { + if(hilite[3][i][hfw - 2] <= 0.01) { + hilite_dir4[3][hfw - 1][i] = hilite_dir0[3][hfw - 2][i]; + } + } + + +#ifdef _OPENMP + #pragma omp parallel for +#endif + + for (int c = 0; c < 3; c++) { + for (int j = 1; j < hfw - 1; j++) { + for (int i = 2; i < hfh - 2; i++) { + //from left + if (hilite[3][i][j] > 0.01f) { + hilite_dir0[c][j][i] = hilite[c][i][j] / hilite[3][i][j]; + } else { + hilite_dir0[c][j][i] = 0.1f * ((hilite_dir0[0 + c][j - 1][i - 2] + hilite_dir0[0 + c][j - 1][i - 1] + hilite_dir0[0 + c][j - 1][i] + hilite_dir0[0 + c][j - 1][i + 1] + hilite_dir0[0 + c][j - 1][i + 2]) / + (hilite_dir0[0 + 3][j - 1][i - 2] + hilite_dir0[0 + 3][j - 1][i - 1] + hilite_dir0[0 + 3][j - 1][i] + hilite_dir0[0 + 3][j - 1][i + 1] + hilite_dir0[0 + 3][j - 1][i + 2] + 0.00001f)); + } + } + + if(hilite[3][2][j] <= 0.01f) { + hilite_dir[0 + c][0][j] = hilite_dir0[c][j][2]; + } + + if(hilite[3][3][j] <= 0.01f) { + hilite_dir[0 + c][1][j] = hilite_dir0[c][j][3]; + } + + if(hilite[3][hfh - 3][j] <= 0.01f) { + hilite_dir[4 + c][hfh - 1][j] = hilite_dir0[c][j][hfh - 3]; + } + + if(hilite[3][hfh - 4][j] <= 0.01f) { + hilite_dir[4 + c][hfh - 2][j] = hilite_dir0[c][j][hfh - 4]; + } + } + + for (int i = 2; i < hfh - 2; i++) { + if(hilite[3][i][hfw - 2] <= 0.01f) { + hilite_dir4[c][hfw - 1][i] = hilite_dir0[c][hfw - 2][i]; + } + } + } + if(plistener) { progress += 0.05; plistener->setProgress(progress); } - for (int j = hfw - 2; j > 0; j--) + for (int j = hfw - 2; j > 0; j--) { for (int i = 2; i < hfh - 2; i++) { //from right if (hilite[3][i][j] > 0.01) { - for (int c = 0; c < 4; c++) { - hilite_dir[4 + c][i][j] = hilite[c][i][j] / hilite[3][i][j]; - } + hilite_dir4[3][j][i] = 1.f; } else { - for (int c = 0; c < 4; c++) { - hilite_dir[4 + c][i][j] = 0.1 * ((hilite_dir[4 + c][(i - 2)][(j + 1)] + hilite_dir[4 + c][(i - 1)][(j + 1)] + hilite_dir[4 + c][(i)][(j + 1)] + hilite_dir[4 + c][(i + 1)][(j + 1)] + hilite_dir[4 + c][(i + 2)][(j + 1)]) / - (hilite_dir[4 + 3][(i - 2)][(j + 1)] + hilite_dir[4 + 3][(i - 1)][(j + 1)] + hilite_dir[4 + 3][(i)][(j + 1)] + hilite_dir[4 + 3][(i + 1)][(j + 1)] + hilite_dir[4 + 3][(i + 2)][(j + 1)] + 0.00001)); - hilite_dir[8 + c][i - 2][j] += hilite_dir[4 + c][i][j]; - hilite_dir[12 + c][i + 2][j] += hilite_dir[4 + c][i][j]; - } + hilite_dir4[3][j][i] = (hilite_dir4[3][(j + 1)][(i - 2)] + hilite_dir4[3][(j + 1)][(i - 1)] + hilite_dir4[3][(j + 1)][(i)] + hilite_dir4[3][(j + 1)][(i + 1)] + hilite_dir4[3][(j + 1)][(i + 2)]) == 0.f ? 0.f : 0.1f; } } + if(hilite[3][2][j] <= 0.01) { + hilite_dir[0 + 3][0][j] += hilite_dir4[3][j][2]; + } + + if(hilite[3][hfh - 3][j] <= 0.01) { + hilite_dir[4 + 3][hfh - 1][j] += hilite_dir4[3][j][hfh - 3]; + } + } + + for (int i = 2; i < hfh - 2; i++) { + if(hilite[3][i][0] <= 0.01) { + hilite_dir[0 + 3][i - 2][0] += hilite_dir4[3][0][i]; + hilite_dir[4 + 3][i + 2][0] += hilite_dir4[3][0][i]; + } + + if(hilite[3][i][1] <= 0.01) { + hilite_dir[0 + 3][i - 2][1] += hilite_dir4[3][1][i]; + hilite_dir[4 + 3][i + 2][1] += hilite_dir4[3][1][i]; + } + + if(hilite[3][i][hfw - 2] <= 0.01) { + hilite_dir[0 + 3][i - 2][hfw - 2] += hilite_dir4[3][hfw - 2][i]; + hilite_dir[4 + 3][i + 2][hfw - 2] += hilite_dir4[3][hfw - 2][i]; + } + } + +#ifdef _OPENMP + #pragma omp parallel for +#endif + + for (int c = 0; c < 3; c++) { + for (int j = hfw - 2; j > 0; j--) { + for (int i = 2; i < hfh - 2; i++) { + //from right + if (hilite[3][i][j] > 0.01) { + hilite_dir4[c][j][i] = hilite[c][i][j] / hilite[3][i][j]; + } else { + hilite_dir4[c][j][i] = 0.1 * ((hilite_dir4[c][(j + 1)][(i - 2)] + hilite_dir4[c][(j + 1)][(i - 1)] + hilite_dir4[c][(j + 1)][(i)] + hilite_dir4[c][(j + 1)][(i + 1)] + hilite_dir4[c][(j + 1)][(i + 2)]) / + (hilite_dir4[3][(j + 1)][(i - 2)] + hilite_dir4[3][(j + 1)][(i - 1)] + hilite_dir4[3][(j + 1)][(i)] + hilite_dir4[3][(j + 1)][(i + 1)] + hilite_dir4[3][(j + 1)][(i + 2)] + 0.00001)); + } + } + + if(hilite[3][2][j] <= 0.01) { + hilite_dir[0 + c][0][j] += hilite_dir4[c][j][2]; + } + + if(hilite[3][hfh - 3][j] <= 0.01) { + hilite_dir[4 + c][hfh - 1][j] += hilite_dir4[c][j][hfh - 3]; + } + } + + for (int i = 2; i < hfh - 2; i++) { + if(hilite[3][i][0] <= 0.01) { + hilite_dir[0 + c][i - 2][0] += hilite_dir4[c][0][i]; + hilite_dir[4 + c][i + 2][0] += hilite_dir4[c][0][i]; + } + + if(hilite[3][i][1] <= 0.01) { + hilite_dir[0 + c][i - 2][1] += hilite_dir4[c][1][i]; + hilite_dir[4 + c][i + 2][1] += hilite_dir4[c][1][i]; + } + + if(hilite[3][i][hfw - 2] <= 0.01) { + hilite_dir[0 + c][i - 2][hfw - 2] += hilite_dir4[c][hfw - 2][i]; + hilite_dir[4 + c][i + 2][hfw - 2] += hilite_dir4[c][hfw - 2][i]; + } + } + } + if(plistener) { progress += 0.05; plistener->setProgress(progress); } + for (int i = 1; i < hfh - 1; i++) for (int j = 2; j < hfw - 2; j++) { - //if (i%100==0 && j%100==0) - // printf("row=%d col=%d r=%f g=%f b=%f hilite=%f \n",i,j,hilite[0][i][j],hilite[1][i][j],hilite[2][i][j],hilite[3][i][j]); - //from top if (hilite[3][i][j] > 0.01) { - for (int c = 0; c < 4; c++) { - hilite_dir[8 + c][i][j] = hilite[c][i][j] / hilite[3][i][j]; - } + hilite_dir[0 + 3][i][j] = 1.f; } else { - for (int c = 0; c < 4; c++) { - hilite_dir[8 + c][i][j] = 0.1 * ((hilite_dir[8 + c][i - 1][j - 2] + hilite_dir[8 + c][i - 1][j - 1] + hilite_dir[8 + c][i - 1][j] + hilite_dir[8 + c][i - 1][j + 1] + hilite_dir[8 + c][i - 1][j + 2]) / - (hilite_dir[8 + 3][i - 1][j - 2] + hilite_dir[8 + 3][i - 1][j - 1] + hilite_dir[8 + 3][i - 1][j] + hilite_dir[8 + 3][i - 1][j + 1] + hilite_dir[8 + 3][i - 1][j + 2] + 0.00001)); - hilite_dir[12 + c][i + 1][j] += hilite_dir[8 + c][i][j]; + hilite_dir[0 + 3][i][j] = (hilite_dir[0 + 3][i - 1][j - 2] + hilite_dir[0 + 3][i - 1][j - 1] + hilite_dir[0 + 3][i - 1][j] + hilite_dir[0 + 3][i - 1][j + 1] + hilite_dir[0 + 3][i - 1][j + 2]) == 0.f ? 0.f : 0.1f; + } + } + + for (int j = 2; j < hfw - 2; j++) { + if(hilite[3][hfh - 2][j] <= 0.01) { + hilite_dir[4 + 3][hfh - 1][j] += hilite_dir[0 + 3][hfh - 2][j]; + } + } + +#ifdef _OPENMP + #pragma omp parallel for +#endif + + for (int c = 0; c < 3; c++) { + for (int i = 1; i < hfh - 1; i++) { + for (int j = 2; j < hfw - 2; j++) { + //from top + if (hilite[3][i][j] > 0.01) { + hilite_dir[0 + c][i][j] = hilite[c][i][j] / hilite[3][i][j]; + } else { + hilite_dir[0 + c][i][j] = 0.1 * ((hilite_dir[0 + c][i - 1][j - 2] + hilite_dir[0 + c][i - 1][j - 1] + hilite_dir[0 + c][i - 1][j] + hilite_dir[0 + c][i - 1][j + 1] + hilite_dir[0 + c][i - 1][j + 2]) / + (hilite_dir[0 + 3][i - 1][j - 2] + hilite_dir[0 + 3][i - 1][j - 1] + hilite_dir[0 + 3][i - 1][j] + hilite_dir[0 + 3][i - 1][j + 1] + hilite_dir[0 + 3][i - 1][j + 2] + 0.00001)); } } } + for (int j = 2; j < hfw - 2; j++) { + if(hilite[3][hfh - 2][j] <= 0.01) { + hilite_dir[4 + c][hfh - 1][j] += hilite_dir[0 + c][hfh - 2][j]; + } + } + + } + if(plistener) { progress += 0.05; plistener->setProgress(progress); @@ -648,17 +805,29 @@ void RawImageSource :: HLRecovery_inpaint (float** red, float** green, float** b for (int j = 2; j < hfw - 2; j++) { //from bottom if (hilite[3][i][j] > 0.01) { - for (int c = 0; c < 4; c++) { - hilite_dir[12 + c][i][j] = hilite[c][i][j] / hilite[3][i][j]; - } + hilite_dir[4 + 3][i][j] = 1.f; } else { - for (int c = 0; c < 4; c++) { - hilite_dir[12 + c][i][j] = 0.1 * ((hilite_dir[12 + c][(i + 1)][(j - 2)] + hilite_dir[12 + c][(i + 1)][(j - 1)] + hilite_dir[12 + c][(i + 1)][(j)] + hilite_dir[12 + c][(i + 1)][(j + 1)] + hilite_dir[12 + c][(i + 1)][(j + 2)]) / - (hilite_dir[12 + 3][(i + 1)][(j - 2)] + hilite_dir[12 + 3][(i + 1)][(j - 1)] + hilite_dir[12 + 3][(i + 1)][(j)] + hilite_dir[12 + 3][(i + 1)][(j + 1)] + hilite_dir[12 + 3][(i + 1)][(j + 2)] + 0.00001)); + hilite_dir[4 + 3][i][j] = (hilite_dir[4 + 3][(i + 1)][(j - 2)] + hilite_dir[4 + 3][(i + 1)][(j - 1)] + hilite_dir[4 + 3][(i + 1)][(j)] + hilite_dir[4 + 3][(i + 1)][(j + 1)] + hilite_dir[4 + 3][(i + 1)][(j + 2)]) == 0.f ? 0.f : 0.1f; + } + } + +#ifdef _OPENMP + #pragma omp parallel for +#endif + + for (int c = 0; c < 4; c++) { + for (int i = hfh - 2; i > 0; i--) { + for (int j = 2; j < hfw - 2; j++) { + //from bottom + if (hilite[3][i][j] > 0.01) { + hilite_dir[4 + c][i][j] = hilite[c][i][j] / hilite[3][i][j]; + } else { + hilite_dir[4 + c][i][j] = 0.1 * ((hilite_dir[4 + c][(i + 1)][(j - 2)] + hilite_dir[4 + c][(i + 1)][(j - 1)] + hilite_dir[4 + c][(i + 1)][(j)] + hilite_dir[4 + c][(i + 1)][(j + 1)] + hilite_dir[4 + c][(i + 1)][(j + 2)]) / + (hilite_dir[4 + 3][(i + 1)][(j - 2)] + hilite_dir[4 + 3][(i + 1)][(j - 1)] + hilite_dir[4 + 3][(i + 1)][(j)] + hilite_dir[4 + 3][(i + 1)][(j + 1)] + hilite_dir[4 + 3][(i + 1)][(j + 2)] + 0.00001)); } } - } + } if(plistener) { progress += 0.05; @@ -666,7 +835,7 @@ void RawImageSource :: HLRecovery_inpaint (float** red, float** green, float** b } //fill in edges - for (int dir = 0; dir < numdirs; dir++) { + for (int dir = 0; dir < 2; dir++) { for (int i = 1; i < hfh - 1; i++) for (int c = 0; c < 4; c++) { hilite_dir[dir * 4 + c][i][0] = hilite_dir[dir * 4 + c][i][1]; @@ -687,27 +856,63 @@ void RawImageSource :: HLRecovery_inpaint (float** red, float** green, float** b } } + for (int i = 1; i < hfh - 1; i++) + for (int c = 0; c < 4; c++) { + hilite_dir0[c][0][i] = hilite_dir0[c][1][i]; + hilite_dir0[c][hfw - 1][i] = hilite_dir0[c][hfw - 2][i]; + } + + for (int j = 1; j < hfw - 1; j++) + for (int c = 0; c < 4; c++) { + hilite_dir0[c][j][0] = hilite_dir0[c][j][1]; + hilite_dir0[c][j][hfh - 1] = hilite_dir0[c][j][hfh - 2]; + } + + for (int c = 0; c < 4; c++) { + hilite_dir0[c][0][0] = hilite_dir0[c][0][1] = hilite_dir0[c][1][0] = hilite_dir0[c][1][1] = hilite_dir0[c][2][2]; + hilite_dir0[c][hfw - 1][0] = hilite_dir0[c][hfw - 1][1] = hilite_dir0[c][hfw - 2][0] = hilite_dir0[c][hfw - 2][1] = hilite_dir0[c][hfw - 3][2]; + hilite_dir0[c][0][hfh - 1] = hilite_dir0[c][0][hfh - 2] = hilite_dir0[c][1][hfh - 1] = hilite_dir0[c][1][hfh - 2] = hilite_dir0[c][2][hfh - 3]; + hilite_dir0[c][hfw - 1][hfh - 1] = hilite_dir0[c][hfw - 1][hfh - 2] = hilite_dir0[c][hfw - 2][hfh - 1] = hilite_dir0[c][hfw - 2][hfh - 2] = hilite_dir0[c][hfw - 3][hfh - 3]; + } + + for (int i = 1; i < hfh - 1; i++) + for (int c = 0; c < 4; c++) { + hilite_dir4[c][0][i] = hilite_dir4[c][1][i]; + hilite_dir4[c][hfw - 1][i] = hilite_dir4[c][hfw - 2][i]; + } + + for (int j = 1; j < hfw - 1; j++) + for (int c = 0; c < 4; c++) { + hilite_dir4[c][j][0] = hilite_dir4[c][j][1]; + hilite_dir4[c][j][hfh - 1] = hilite_dir4[c][j][hfh - 2]; + } + + for (int c = 0; c < 4; c++) { + hilite_dir4[c][0][0] = hilite_dir4[c][0][1] = hilite_dir4[c][1][0] = hilite_dir4[c][1][1] = hilite_dir4[c][2][2]; + hilite_dir4[c][hfw - 1][0] = hilite_dir4[c][hfw - 1][1] = hilite_dir4[c][hfw - 2][0] = hilite_dir4[c][hfw - 2][1] = hilite_dir4[c][hfw - 3][2]; + hilite_dir4[c][0][hfh - 1] = hilite_dir4[c][0][hfh - 2] = hilite_dir4[c][1][hfh - 1] = hilite_dir4[c][1][hfh - 2] = hilite_dir4[c][2][hfh - 3]; + hilite_dir4[c][hfw - 1][hfh - 1] = hilite_dir4[c][hfw - 1][hfh - 2] = hilite_dir4[c][hfw - 2][hfh - 1] = hilite_dir4[c][hfw - 2][hfh - 2] = hilite_dir4[c][hfw - 3][hfh - 3]; + } + + + if(plistener) { progress += 0.05; plistener->setProgress(progress); } - //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - /*for (int m=0; m<4*numdirs; m++) { - boxblur2(hilite_dir[m],hilite_dir[m],hfh,hfw,4); - }*/ - - //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - //now we have highlight data extended in various directions - //next step is to build back local data by averaging - - //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% // now reconstruct clipped channels using color ratios - //const float Yclip = 0.3333*(max[0] + max[1] + max[2]); - //float sumwt=0, counts=0; + for(int c = 0; c < 4; c++) { + hilite[c].free(); + } + + LUTf invfn(0x10000); + + for (int i = 0; i < 0x10000; i++) { + invfn[i] = 1.0 / (1.0 + i); + } #ifdef _OPENMP #pragma omp parallel for schedule(dynamic,16) @@ -725,8 +930,6 @@ void RawImageSource :: HLRecovery_inpaint (float** red, float** green, float** b continue; //pixel not clipped } - //if (pixel[0] 0.5) { + dirwt = invfn[65535 * (SQR(rgb_blend[0] / Y - hilite_dir0[0][j1][i1] / Yhi) + + SQR(rgb_blend[1] / Y - hilite_dir0[1][j1][i1] / Yhi) + + SQR(rgb_blend[2] / Y - hilite_dir0[2][j1][i1] / Yhi))]; + totwt += dirwt; + clipfix[0] += dirwt * hilite_dir0[0][j1][i1] / (hilite_dir0[3][j1][i1] + 0.00001); + clipfix[1] += dirwt * hilite_dir0[1][j1][i1] / (hilite_dir0[3][j1][i1] + 0.00001); + clipfix[2] += dirwt * hilite_dir0[2][j1][i1] / (hilite_dir0[3][j1][i1] + 0.00001); + } + + for (int dir = 0; dir < 2; dir++) { float Yhi = 0.001 + (hilite_dir[dir * 4 + 0][i1][j1] + hilite_dir[dir * 4 + 1][i1][j1] + hilite_dir[dir * 4 + 2][i1][j1]); float Y = 0.001 + (rgb_blend[0] + rgb_blend[1] + rgb_blend[2]); @@ -827,6 +1043,19 @@ void RawImageSource :: HLRecovery_inpaint (float** red, float** green, float** b } } + Yhi = 0.001 + (hilite_dir4[0][j1][i1] + hilite_dir4[1][j1][i1] + hilite_dir4[2][j1][i1]); + Y = 0.001 + (rgb_blend[0] + rgb_blend[1] + rgb_blend[2]); + + if (hilite_dir4[0][j1][i1] + hilite_dir4[1][j1][i1] + hilite_dir4[2][j1][i1] > 0.5) { + dirwt = invfn[65535 * (SQR(rgb_blend[0] / Y - hilite_dir4[0][j1][i1] / Yhi) + + SQR(rgb_blend[1] / Y - hilite_dir4[1][j1][i1] / Yhi) + + SQR(rgb_blend[2] / Y - hilite_dir4[2][j1][i1] / Yhi))]; + totwt += dirwt; + clipfix[0] += dirwt * hilite_dir4[0][j1][i1] / (hilite_dir4[3][j1][i1] + 0.00001); + clipfix[1] += dirwt * hilite_dir4[1][j1][i1] / (hilite_dir4[3][j1][i1] + 0.00001); + clipfix[2] += dirwt * hilite_dir4[2][j1][i1] / (hilite_dir4[3][j1][i1] + 0.00001); + } + if(totwt == 0.f) { continue; } @@ -841,7 +1070,7 @@ void RawImageSource :: HLRecovery_inpaint (float** red, float** green, float** b if (pixel[0] > max_f[0] && pixel[1] > max_f[1] && pixel[2] > max_f[2]) { //all channels clipped float Y = (0.299 * clipfix[0] + 0.587 * clipfix[1] + 0.114 * clipfix[2]); - //float Y = (clipfix[0] + clipfix[1] + clipfix[2]); + factor = whitept / Y; red[i][j] = clipfix[0] * factor; green[i][j] = clipfix[1] * factor; @@ -866,28 +1095,11 @@ void RawImageSource :: HLRecovery_inpaint (float** red, float** green, float** b } } - /*if (hilite[3][i1][j1]>0.01) { - red[i][j] = (red[i][j] + hilite[0][i1][j1])/(1+hilite[3][i1][j1]); - green[i][j] = (green[i][j]+ hilite[1][i1][j1])/(1+hilite[3][i1][j1]); - blue[i][j] = (blue[i][j] + hilite[2][i1][j1])/(1+hilite[3][i1][j1]); - }*/ - Y = (0.299 * red[i][j] + 0.587 * green[i][j] + 0.114 * blue[i][j]); if (Y > whitept) { factor = whitept / Y; - /*I = (0.596 * red[i][j] - 0.275 * green[i][j] - 0.321 * blue[i][j]); - Q = (0.212 * red[i][j] - 0.523 * green[i][j] + 0.311 * blue[i][j]); - - Y *= factor; - I *= factor;//max(0,min(1,(whitept-Y)/(whitept-clippt))); - Q *= factor;//max(0,min(1,(whitept-Y)/(whitept-clippt))); - - red[i][j] = Y + 0.956*I + 0.621*Q; - green[i][j] = Y - 0.272*I - 0.647*Q; - blue[i][j] = Y - 1.105*I + 1.702*Q;*/ - red[i][j] *= factor; green[i][j] *= factor; blue[i][j] *= factor; @@ -899,76 +1111,8 @@ void RawImageSource :: HLRecovery_inpaint (float** red, float** green, float** b plistener->setProgress(1.00); } - //printf("ave wt=%f\n",sumwt/counts); - - - // diagnostic output - /*for (int i=0; i & hfsize, multi_array2D & hilite, int box ); From 6eaf68e8500efd1672508c1026490c6888c49692 Mon Sep 17 00:00:00 2001 From: heckflosse Date: Sun, 11 Oct 2015 19:17:09 +0200 Subject: [PATCH 10/16] Improved quality of Colour Propagation for some cases --- rtengine/hilite_recon.cc | 39 +++++++++++++++++++++++++++++++++++++-- 1 file changed, 37 insertions(+), 2 deletions(-) diff --git a/rtengine/hilite_recon.cc b/rtengine/hilite_recon.cc index 46c666c80..0c8b40fd9 100644 --- a/rtengine/hilite_recon.cc +++ b/rtengine/hilite_recon.cc @@ -421,11 +421,46 @@ void RawImageSource :: HLRecovery_inpaint (float** red, float** green, float** b { { 1, 1, 1 }, { 1, -1, 1 }, { 1, 1, -1 } } }; + + for(int c=0;c<3;c++) + printf("chmax[%d] : %f\tclmax[%d] : %f\tratio[%d] : %f\n",c,chmax[c],c,clmax[c],c,chmax[c]/clmax[c]); + + float factor[3]; + for(int i=0;i<3;i++) + factor[i] = chmax[i]/clmax[i]; + + float minFactor = min(factor[0],factor[1],factor[2]); + if(minFactor > 1.f) { // all 3 channels clipped + // calculate clip factor per channel + for(int c=0;c<3;c++) + factor[c] /= minFactor; + // get max clip factor + int maxpos = 0; + float maxValNew = 0.f; + for(int c=0;c<3;c++) + if(chmax[c]/factor[c] > maxValNew) { + maxValNew = chmax[c]/factor[c]; + maxpos = c; + } + float clipFactor = (clmax[maxpos]) / maxValNew; + if(clipFactor < maxpct) + // if max clipfactor < clippct (0.95) adjust per channel factors + for(int c=0;c<3;c++) { + factor[c] *= (maxpct / clipFactor); + } + } else { + factor[0] = factor[1] = factor[2] = 1.f; + } + + for(int c=0;c<3;c++) { + printf("correction factor [%d] : %f\n",c,factor[c]); + } + float max_f[3], thresh[3]; for (int c = 0; c < 3; c++) { - thresh[c] = chmax[c] * threshpct; - max_f[c] = chmax[c] * maxpct; //min(chmax[0],chmax[1],chmax[2])*maxpct; + thresh[c] = chmax[c] * threshpct / factor[c]; + max_f[c] = chmax[c] * maxpct / factor[c]; } float whitept = max(max_f[0], max_f[1], max_f[2]); From 33182f743cfa80066be24edcb2284f65b29e4cb3 Mon Sep 17 00:00:00 2001 From: heckflosse Date: Mon, 12 Oct 2015 17:57:38 +0200 Subject: [PATCH 11/16] sometimes getimage is not called after demosaic, fixes #2932 --- rtengine/improccoordinator.cc | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/rtengine/improccoordinator.cc b/rtengine/improccoordinator.cc index fa1fae502..933373920 100644 --- a/rtengine/improccoordinator.cc +++ b/rtengine/improccoordinator.cc @@ -220,13 +220,11 @@ void ImProcCoordinator::updatePreviewImage (int todo, Crop* cropCall) } imgsrc->demosaic( rp ); + // if a demosaic happened we should also call getimage later, so we need to set the M_INIT flag + todo |= M_INIT; if (highDetailNeeded) { highDetailRawComputed = true; - - if (params.toneCurve.hrenabled && params.toneCurve.method == "Color") { - todo |= M_INIT; - } } else { highDetailRawComputed = false; } From 9061b2549d4a43f4d7aa51567b54d18c9da4f2ad Mon Sep 17 00:00:00 2001 From: heckflosse Date: Tue, 13 Oct 2015 14:36:19 +0200 Subject: [PATCH 12/16] Fix inconsistent behaviour in Colour Propagation --- rtengine/hilite_recon.cc | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/rtengine/hilite_recon.cc b/rtengine/hilite_recon.cc index 0c8b40fd9..c35eb9f63 100644 --- a/rtengine/hilite_recon.cc +++ b/rtengine/hilite_recon.cc @@ -467,6 +467,10 @@ void RawImageSource :: HLRecovery_inpaint (float** red, float** green, float** b float clippt = min(max_f[0], max_f[1], max_f[2]); float medpt = max_f[0] + max_f[1] + max_f[2] - whitept - clippt; float blendpt = blendthresh * clippt; + float medFactor[3]; + for(int c=0;c<3;c++) { + medFactor[c] = max(1.0f, max_f[c] / medpt) / (-blendpt); + } multi_array2D channelblur(width, height, 0, 48); array2D temp(width, height); // allocate temporary buffer @@ -1023,9 +1027,9 @@ void RawImageSource :: HLRecovery_inpaint (float** red, float** green, float** b FOREACHCOLOR rgb[c] = cam[0][c] / ColorCount; // Copy converted pixel back - float rfrac = min(1.0f, max(1.0f, max_f[0] / medpt) * (pixel[0] - blendpt) / (hlmax[0] - blendpt)); - float gfrac = min(1.0f, max(1.0f, max_f[1] / medpt) * (pixel[1] - blendpt) / (hlmax[1] - blendpt)); - float bfrac = min(1.0f, max(1.0f, max_f[2] / medpt) * (pixel[2] - blendpt) / (hlmax[2] - blendpt)); + float rfrac = max(0.f,min(1.0f, medFactor[0] * (pixel[0] - blendpt))); + float gfrac = max(0.f,min(1.0f, medFactor[1] * (pixel[1] - blendpt))); + float bfrac = max(0.f,min(1.0f, medFactor[2] * (pixel[2] - blendpt))); if (pixel[0] > blendpt) { rgb_blend[0] = rfrac * rgb[0] + (1 - rfrac) * pixel[0]; From 17eca17e61924e551b6bb433f23b9474b8059356 Mon Sep 17 00:00:00 2001 From: Marcin Bajor Date: Wed, 14 Oct 2015 10:43:33 +0200 Subject: [PATCH 13/16] Fix for build with sigc++ >= 2.5.2 --- rtgui/adjuster.cc | 2 +- rtgui/preferences.cc | 2 +- rtgui/thresholdadjuster.cc | 2 +- rtgui/tonecurve.cc | 2 +- 4 files changed, 4 insertions(+), 4 deletions(-) diff --git a/rtgui/adjuster.cc b/rtgui/adjuster.cc index bbd6e76cf..f0e53fbbf 100644 --- a/rtgui/adjuster.cc +++ b/rtgui/adjuster.cc @@ -17,7 +17,7 @@ * along with RawTherapee. If not, see . */ #include "adjuster.h" -#include +#include #include #include "multilangmgr.h" #include "../rtengine/rtengine.h" diff --git a/rtgui/preferences.cc b/rtgui/preferences.cc index 9c0918189..141ad35d3 100644 --- a/rtgui/preferences.cc +++ b/rtgui/preferences.cc @@ -16,7 +16,7 @@ * You should have received a copy of the GNU General Public License * along with RawTherapee. If not, see . */ -#include +#include #include "preferences.h" #include "multilangmgr.h" #include "splash.h" diff --git a/rtgui/thresholdadjuster.cc b/rtgui/thresholdadjuster.cc index a7e21a26d..800f6edf5 100644 --- a/rtgui/thresholdadjuster.cc +++ b/rtgui/thresholdadjuster.cc @@ -17,7 +17,7 @@ * along with RawTherapee. If not, see . */ #include "thresholdadjuster.h" -#include +#include #include #include "multilangmgr.h" #include "../rtengine/rtengine.h" diff --git a/rtgui/tonecurve.cc b/rtgui/tonecurve.cc index 770652430..799513c2a 100644 --- a/rtgui/tonecurve.cc +++ b/rtgui/tonecurve.cc @@ -18,7 +18,7 @@ */ #include "tonecurve.h" #include "adjuster.h" -#include +#include #include #include "ppversion.h" #include "edit.h" From c6620b415db72cbdf8833476a21af5190e15fa22 Mon Sep 17 00:00:00 2001 From: Marcin Bajor Date: Wed, 14 Oct 2015 10:51:27 +0200 Subject: [PATCH 14/16] Fix for build with C++11 flags --- rtengine/dcraw.cc | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/rtengine/dcraw.cc b/rtengine/dcraw.cc index add5b525d..8e0a36df2 100644 --- a/rtengine/dcraw.cc +++ b/rtengine/dcraw.cc @@ -136,10 +136,10 @@ const float d65_white[3] = { 0.950456, 1, 1.088754 }; #define SQR(x) rtengine::SQR(x) #define ABS(x) (((int)(x) ^ ((int)(x) >> 31)) - ((int)(x) >> 31)) -#define MIN(a,b) rtengine::min(a,static_cast(b)) -#define MAX(a,b) rtengine::max(a,static_cast(b)) -#define LIM(x,min,max) rtengine::LIM(x,static_cast(min),static_cast(max)) -#define ULIM(x,y,z) rtengine::ULIM(x,static_cast(y),static_cast(z)) +#define MIN(a,b) rtengine::min(a,static_cast<__typeof__(a)>(b)) +#define MAX(a,b) rtengine::max(a,static_cast<__typeof__(a)>(b)) +#define LIM(x,min,max) rtengine::LIM(x,static_cast<__typeof__(x)>(min),static_cast<__typeof__(x)>(max)) +#define ULIM(x,y,z) rtengine::ULIM(x,static_cast<__typeof__(x)>(y),static_cast<__typeof__(x)>(z)) #define CLIP(x) rtengine::CLIP(x) #define SWAP(a,b) { a=a+b; b=a-b; a=a-b; } From aef0fc7ec8847b75650a86fd63b31bc0dc63ab49 Mon Sep 17 00:00:00 2001 From: heckflosse Date: Wed, 14 Oct 2015 16:37:54 +0200 Subject: [PATCH 15/16] Updated dcraw.patch file --- rtengine/dcraw.patch | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/rtengine/dcraw.patch b/rtengine/dcraw.patch index e1603864a..db24e7fe8 100644 --- a/rtengine/dcraw.patch +++ b/rtengine/dcraw.patch @@ -1,5 +1,5 @@ ---- dcraw.c 2015-08-15 15:35:27 +0000 -+++ dcraw.cc 2015-08-16 13:46:33 +0000 +--- dcraw.c 2015-09-21 10:08:04 +0000 ++++ dcraw.cc 2015-10-14 14:29:55 +0000 @@ -1,3 +1,15 @@ +/*RT*/#include +/*RT*/#include @@ -148,10 +148,10 @@ -#define LIM(x,min,max) MAX(min,MIN(x,max)) -#define ULIM(x,y,z) ((y) < (z) ? LIM(x,y,z) : LIM(x,z,y)) -#define CLIP(x) LIM(x,0,65535) -+#define MIN(a,b) rtengine::min(a,static_cast(b)) -+#define MAX(a,b) rtengine::max(a,static_cast(b)) -+#define LIM(x,min,max) rtengine::LIM(x,static_cast(min),static_cast(max)) -+#define ULIM(x,y,z) rtengine::ULIM(x,static_cast(y),static_cast(z)) ++#define MIN(a,b) rtengine::min(a,static_cast<__typeof__(a)>(b)) ++#define MAX(a,b) rtengine::max(a,static_cast<__typeof__(a)>(b)) ++#define LIM(x,min,max) rtengine::LIM(x,static_cast<__typeof__(x)>(min),static_cast<__typeof__(x)>(max)) ++#define ULIM(x,y,z) rtengine::ULIM(x,static_cast<__typeof__(x)>(y),static_cast<__typeof__(x)>(z)) +#define CLIP(x) rtengine::CLIP(x) #define SWAP(a,b) { a=a+b; b=a-b; a=a-b; } From a0dabf8067ef34fa5fedb0ef142bf82a28f11c89 Mon Sep 17 00:00:00 2001 From: heckflosse Date: Wed, 14 Oct 2015 23:15:01 +0200 Subject: [PATCH 16/16] Colour Propagation, cleanup and a small speedup --- rtengine/hilite_recon.cc | 312 +++++++++++++++++++-------------------- 1 file changed, 156 insertions(+), 156 deletions(-) diff --git a/rtengine/hilite_recon.cc b/rtengine/hilite_recon.cc index c35eb9f63..9f618ace5 100644 --- a/rtengine/hilite_recon.cc +++ b/rtengine/hilite_recon.cc @@ -28,12 +28,11 @@ #include "rawimagesource.h" #include "rt_math.h" #include "opthelper.h" - -#define FOREACHCOLOR for (int c=0; c < ColorCount; c++) - namespace rtengine { +extern const Settings* settings; + //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% SSEFUNCTION void RawImageSource::boxblur2(float** src, float** dst, float** temp, int H, int W, int box ) @@ -152,6 +151,7 @@ SSEFUNCTION void RawImageSource::boxblur2(float** src, float** dst, float** temp _mm_storeu_ps( &dst[row][col], tempv ); } } + for (int col = W - (W % 4); col < W; col++) { int len = box + 1; dst[0][col] = temp[0][col] / len; @@ -403,62 +403,70 @@ void RawImageSource :: HLRecovery_inpaint (float** red, float** green, float** b static const int numdirs = 4; - static const float threshpct = 0.25; - static const float fixthreshpct = 0.7; - static const float maxpct = 0.95; - + static const float threshpct = 0.25f; + static const float fixthreshpct = 0.7f; + static const float maxpct = 0.95f; + static const float epsilon = 0.00001f; //%%%%%%%%%%%%%%%%%%%% //for blend algorithm: static const float blendthresh = 1.0; static const int ColorCount = 3; // Transform matrixes rgb>lab and back - static const float trans[2][ColorCount][ColorCount] = { - { { 1, 1, 1 }, { 1.7320508, -1.7320508, 0 }, { -1, -1, 2 } }, - { { 1, 1, 1 }, { 1, -1, 1 }, { 1, 1, -1 } } - }; - static const float itrans[2][ColorCount][ColorCount] = { - { { 1, 0.8660254, -0.5 }, { 1, -0.8660254, -0.5 }, { 1, 0, 1 } }, - { { 1, 1, 1 }, { 1, -1, 1 }, { 1, 1, -1 } } - }; + static const float trans[ColorCount][ColorCount] = + { { 1.f, 1.f, 1.f }, { 1.7320508f, -1.7320508f, 0.f }, { -1.f, -1.f, 2.f } }; + static const float itrans[ColorCount][ColorCount] = + { { 1.f, 0.8660254f, -0.5f }, { 1.f, -0.8660254f, -0.5f }, { 1.f, 0.f, 1.f } }; - - for(int c=0;c<3;c++) - printf("chmax[%d] : %f\tclmax[%d] : %f\tratio[%d] : %f\n",c,chmax[c],c,clmax[c],c,chmax[c]/clmax[c]); + if(settings->verbose) + for(int c = 0; c < 3; c++) { + printf("chmax[%d] : %f\tclmax[%d] : %f\tratio[%d] : %f\n", c, chmax[c], c, clmax[c], c, chmax[c] / clmax[c]); + } float factor[3]; - for(int i=0;i<3;i++) - factor[i] = chmax[i]/clmax[i]; - float minFactor = min(factor[0],factor[1],factor[2]); + for(int c = 0; c < ColorCount; c++) { + factor[c] = chmax[c] / clmax[c]; + } + + float minFactor = min(factor[0], factor[1], factor[2]); + if(minFactor > 1.f) { // all 3 channels clipped // calculate clip factor per channel - for(int c=0;c<3;c++) + for (int c = 0; c < ColorCount; c++) { factor[c] /= minFactor; + } + // get max clip factor int maxpos = 0; float maxValNew = 0.f; - for(int c=0;c<3;c++) - if(chmax[c]/factor[c] > maxValNew) { - maxValNew = chmax[c]/factor[c]; + + for (int c = 0; c < ColorCount; c++) { + if(chmax[c] / factor[c] > maxValNew) { + maxValNew = chmax[c] / factor[c]; maxpos = c; } - float clipFactor = (clmax[maxpos]) / maxValNew; + } + + float clipFactor = clmax[maxpos] / maxValNew; + if(clipFactor < maxpct) - // if max clipfactor < clippct (0.95) adjust per channel factors - for(int c=0;c<3;c++) { + + // if max clipFactor < maxpct (0.95) adjust per channel factors + for (int c = 0; c < ColorCount; c++) { factor[c] *= (maxpct / clipFactor); } } else { factor[0] = factor[1] = factor[2] = 1.f; } - for(int c=0;c<3;c++) { - printf("correction factor [%d] : %f\n",c,factor[c]); - } + if(settings->verbose) + for (int c = 0; c < ColorCount; c++) { + printf("correction factor[%d] : %f\n", c, factor[c]); + } float max_f[3], thresh[3]; - for (int c = 0; c < 3; c++) { + for (int c = 0; c < ColorCount; c++) { thresh[c] = chmax[c] * threshpct / factor[c]; max_f[c] = chmax[c] * maxpct / factor[c]; } @@ -468,7 +476,8 @@ void RawImageSource :: HLRecovery_inpaint (float** red, float** green, float** b float medpt = max_f[0] + max_f[1] + max_f[2] - whitept - clippt; float blendpt = blendthresh * clippt; float medFactor[3]; - for(int c=0;c<3;c++) { + + for (int c = 0; c < ColorCount; c++) { medFactor[c] = max(1.0f, max_f[c] / medpt) / (-blendpt); } @@ -477,19 +486,22 @@ void RawImageSource :: HLRecovery_inpaint (float** red, float** green, float** b // blur RGB channels - boxblur2(red , channelblur[0], temp, height, width, 4); + boxblur2(red, channelblur[0], temp, height, width, 4); + if(plistener) { progress += 0.05; plistener->setProgress(progress); } boxblur2(green, channelblur[1], temp, height, width, 4); + if(plistener) { progress += 0.05; plistener->setProgress(progress); } - boxblur2(blue , channelblur[2], temp, height, width, 4); + boxblur2(blue, channelblur[2], temp, height, width, 4); + if(plistener) { progress += 0.05; plistener->setProgress(progress); @@ -515,12 +527,12 @@ void RawImageSource :: HLRecovery_inpaint (float** red, float** green, float** b } multi_array2D hilite_full(width, height, ARRAY2D_CLEAR_DATA, 32); + if(plistener) { progress += 0.10; plistener->setProgress(progress); } - double hipass_sum = 0.f; int hipass_norm = 0; @@ -547,7 +559,8 @@ void RawImageSource :: HLRecovery_inpaint (float** red, float** green, float** b } }//end of filling highlight array - float hipass_ave = 2.f * hipass_sum / (hipass_norm + 0.01); + float hipass_ave = 2.f * hipass_sum / (hipass_norm + epsilon); + if(plistener) { progress += 0.05; plistener->setProgress(progress); @@ -577,7 +590,7 @@ void RawImageSource :: HLRecovery_inpaint (float** red, float** green, float** b continue; } - if (hilite_full4[i][j] > 0.00001 && hilite_full4[i][j] < 0.95) { + if (hilite_full4[i][j] > epsilon && hilite_full4[i][j] < 0.95f) { //too near an edge, could risk using CA affected pixels, therefore omit hilite_full[0][i][j] = hilite_full[1][i][j] = hilite_full[2][i][j] = hilite_full[3][i][j] = 0.f; } @@ -616,49 +629,43 @@ void RawImageSource :: HLRecovery_inpaint (float** red, float** green, float** b // for faster processing we create two buffers using (height,width) instead of (width,height) multi_array2D hilite_dir0(hfh, hfw, ARRAY2D_CLEAR_DATA, 64); multi_array2D hilite_dir4(hfh, hfw, ARRAY2D_CLEAR_DATA, 64); - //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - if(plistener) { progress += 0.05; plistener->setProgress(progress); } - - //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - //fill gaps in highlight map by directional extension //raster scan from four corners for (int j = 1; j < hfw - 1; j++) { for (int i = 2; i < hfh - 2; i++) { //from left - if (hilite[3][i][j] > 0.01) { + if (hilite[3][i][j] > epsilon) { hilite_dir0[3][j][i] = 1.f; } else { hilite_dir0[3][j][i] = (hilite_dir0[0 + 3][j - 1][i - 2] + hilite_dir0[0 + 3][j - 1][i - 1] + hilite_dir0[0 + 3][j - 1][i] + hilite_dir0[0 + 3][j - 1][i + 1] + hilite_dir0[0 + 3][j - 1][i + 2]) == 0.f ? 0.f : 0.1f; } } - if(hilite[3][2][j] <= 0.01) { + if(hilite[3][2][j] <= epsilon) { hilite_dir[0 + 3][0][j] = hilite_dir0[3][j][2]; } - if(hilite[3][3][j] <= 0.01) { + if(hilite[3][3][j] <= epsilon) { hilite_dir[0 + 3][1][j] = hilite_dir0[3][j][3]; } - if(hilite[3][hfh - 3][j] <= 0.01) { + if(hilite[3][hfh - 3][j] <= epsilon) { hilite_dir[4 + 3][hfh - 1][j] = hilite_dir0[3][j][hfh - 3]; } - if(hilite[3][hfh - 4][j] <= 0.01) { + if(hilite[3][hfh - 4][j] <= epsilon) { hilite_dir[4 + 3][hfh - 2][j] = hilite_dir0[3][j][hfh - 4]; } } for (int i = 2; i < hfh - 2; i++) { - if(hilite[3][i][hfw - 2] <= 0.01) { + if(hilite[3][i][hfw - 2] <= epsilon) { hilite_dir4[3][hfw - 1][i] = hilite_dir0[3][hfw - 2][i]; } } @@ -672,33 +679,33 @@ void RawImageSource :: HLRecovery_inpaint (float** red, float** green, float** b for (int j = 1; j < hfw - 1; j++) { for (int i = 2; i < hfh - 2; i++) { //from left - if (hilite[3][i][j] > 0.01f) { + if (hilite[3][i][j] > epsilon) { hilite_dir0[c][j][i] = hilite[c][i][j] / hilite[3][i][j]; } else { hilite_dir0[c][j][i] = 0.1f * ((hilite_dir0[0 + c][j - 1][i - 2] + hilite_dir0[0 + c][j - 1][i - 1] + hilite_dir0[0 + c][j - 1][i] + hilite_dir0[0 + c][j - 1][i + 1] + hilite_dir0[0 + c][j - 1][i + 2]) / - (hilite_dir0[0 + 3][j - 1][i - 2] + hilite_dir0[0 + 3][j - 1][i - 1] + hilite_dir0[0 + 3][j - 1][i] + hilite_dir0[0 + 3][j - 1][i + 1] + hilite_dir0[0 + 3][j - 1][i + 2] + 0.00001f)); + (hilite_dir0[0 + 3][j - 1][i - 2] + hilite_dir0[0 + 3][j - 1][i - 1] + hilite_dir0[0 + 3][j - 1][i] + hilite_dir0[0 + 3][j - 1][i + 1] + hilite_dir0[0 + 3][j - 1][i + 2] + epsilon)); } } - if(hilite[3][2][j] <= 0.01f) { + if(hilite[3][2][j] <= epsilon) { hilite_dir[0 + c][0][j] = hilite_dir0[c][j][2]; } - if(hilite[3][3][j] <= 0.01f) { + if(hilite[3][3][j] <= epsilon) { hilite_dir[0 + c][1][j] = hilite_dir0[c][j][3]; } - if(hilite[3][hfh - 3][j] <= 0.01f) { + if(hilite[3][hfh - 3][j] <= epsilon) { hilite_dir[4 + c][hfh - 1][j] = hilite_dir0[c][j][hfh - 3]; } - if(hilite[3][hfh - 4][j] <= 0.01f) { + if(hilite[3][hfh - 4][j] <= epsilon) { hilite_dir[4 + c][hfh - 2][j] = hilite_dir0[c][j][hfh - 4]; } } for (int i = 2; i < hfh - 2; i++) { - if(hilite[3][i][hfw - 2] <= 0.01f) { + if(hilite[3][i][hfw - 2] <= epsilon) { hilite_dir4[c][hfw - 1][i] = hilite_dir0[c][hfw - 2][i]; } } @@ -712,34 +719,34 @@ void RawImageSource :: HLRecovery_inpaint (float** red, float** green, float** b for (int j = hfw - 2; j > 0; j--) { for (int i = 2; i < hfh - 2; i++) { //from right - if (hilite[3][i][j] > 0.01) { + if (hilite[3][i][j] > epsilon) { hilite_dir4[3][j][i] = 1.f; } else { hilite_dir4[3][j][i] = (hilite_dir4[3][(j + 1)][(i - 2)] + hilite_dir4[3][(j + 1)][(i - 1)] + hilite_dir4[3][(j + 1)][(i)] + hilite_dir4[3][(j + 1)][(i + 1)] + hilite_dir4[3][(j + 1)][(i + 2)]) == 0.f ? 0.f : 0.1f; } } - if(hilite[3][2][j] <= 0.01) { + if(hilite[3][2][j] <= epsilon) { hilite_dir[0 + 3][0][j] += hilite_dir4[3][j][2]; } - if(hilite[3][hfh - 3][j] <= 0.01) { + if(hilite[3][hfh - 3][j] <= epsilon) { hilite_dir[4 + 3][hfh - 1][j] += hilite_dir4[3][j][hfh - 3]; } } for (int i = 2; i < hfh - 2; i++) { - if(hilite[3][i][0] <= 0.01) { + if(hilite[3][i][0] <= epsilon) { hilite_dir[0 + 3][i - 2][0] += hilite_dir4[3][0][i]; hilite_dir[4 + 3][i + 2][0] += hilite_dir4[3][0][i]; } - if(hilite[3][i][1] <= 0.01) { + if(hilite[3][i][1] <= epsilon) { hilite_dir[0 + 3][i - 2][1] += hilite_dir4[3][1][i]; hilite_dir[4 + 3][i + 2][1] += hilite_dir4[3][1][i]; } - if(hilite[3][i][hfw - 2] <= 0.01) { + if(hilite[3][i][hfw - 2] <= epsilon) { hilite_dir[0 + 3][i - 2][hfw - 2] += hilite_dir4[3][hfw - 2][i]; hilite_dir[4 + 3][i + 2][hfw - 2] += hilite_dir4[3][hfw - 2][i]; } @@ -753,35 +760,35 @@ void RawImageSource :: HLRecovery_inpaint (float** red, float** green, float** b for (int j = hfw - 2; j > 0; j--) { for (int i = 2; i < hfh - 2; i++) { //from right - if (hilite[3][i][j] > 0.01) { + if (hilite[3][i][j] > epsilon) { hilite_dir4[c][j][i] = hilite[c][i][j] / hilite[3][i][j]; } else { hilite_dir4[c][j][i] = 0.1 * ((hilite_dir4[c][(j + 1)][(i - 2)] + hilite_dir4[c][(j + 1)][(i - 1)] + hilite_dir4[c][(j + 1)][(i)] + hilite_dir4[c][(j + 1)][(i + 1)] + hilite_dir4[c][(j + 1)][(i + 2)]) / - (hilite_dir4[3][(j + 1)][(i - 2)] + hilite_dir4[3][(j + 1)][(i - 1)] + hilite_dir4[3][(j + 1)][(i)] + hilite_dir4[3][(j + 1)][(i + 1)] + hilite_dir4[3][(j + 1)][(i + 2)] + 0.00001)); + (hilite_dir4[3][(j + 1)][(i - 2)] + hilite_dir4[3][(j + 1)][(i - 1)] + hilite_dir4[3][(j + 1)][(i)] + hilite_dir4[3][(j + 1)][(i + 1)] + hilite_dir4[3][(j + 1)][(i + 2)] + epsilon)); } } - if(hilite[3][2][j] <= 0.01) { + if(hilite[3][2][j] <= epsilon) { hilite_dir[0 + c][0][j] += hilite_dir4[c][j][2]; } - if(hilite[3][hfh - 3][j] <= 0.01) { + if(hilite[3][hfh - 3][j] <= epsilon) { hilite_dir[4 + c][hfh - 1][j] += hilite_dir4[c][j][hfh - 3]; } } for (int i = 2; i < hfh - 2; i++) { - if(hilite[3][i][0] <= 0.01) { + if(hilite[3][i][0] <= epsilon) { hilite_dir[0 + c][i - 2][0] += hilite_dir4[c][0][i]; hilite_dir[4 + c][i + 2][0] += hilite_dir4[c][0][i]; } - if(hilite[3][i][1] <= 0.01) { + if(hilite[3][i][1] <= epsilon) { hilite_dir[0 + c][i - 2][1] += hilite_dir4[c][1][i]; hilite_dir[4 + c][i + 2][1] += hilite_dir4[c][1][i]; } - if(hilite[3][i][hfw - 2] <= 0.01) { + if(hilite[3][i][hfw - 2] <= epsilon) { hilite_dir[0 + c][i - 2][hfw - 2] += hilite_dir4[c][hfw - 2][i]; hilite_dir[4 + c][i + 2][hfw - 2] += hilite_dir4[c][hfw - 2][i]; } @@ -797,7 +804,7 @@ void RawImageSource :: HLRecovery_inpaint (float** red, float** green, float** b for (int i = 1; i < hfh - 1; i++) for (int j = 2; j < hfw - 2; j++) { //from top - if (hilite[3][i][j] > 0.01) { + if (hilite[3][i][j] > epsilon) { hilite_dir[0 + 3][i][j] = 1.f; } else { hilite_dir[0 + 3][i][j] = (hilite_dir[0 + 3][i - 1][j - 2] + hilite_dir[0 + 3][i - 1][j - 1] + hilite_dir[0 + 3][i - 1][j] + hilite_dir[0 + 3][i - 1][j + 1] + hilite_dir[0 + 3][i - 1][j + 2]) == 0.f ? 0.f : 0.1f; @@ -805,7 +812,7 @@ void RawImageSource :: HLRecovery_inpaint (float** red, float** green, float** b } for (int j = 2; j < hfw - 2; j++) { - if(hilite[3][hfh - 2][j] <= 0.01) { + if(hilite[3][hfh - 2][j] <= epsilon) { hilite_dir[4 + 3][hfh - 1][j] += hilite_dir[0 + 3][hfh - 2][j]; } } @@ -818,21 +825,20 @@ void RawImageSource :: HLRecovery_inpaint (float** red, float** green, float** b for (int i = 1; i < hfh - 1; i++) { for (int j = 2; j < hfw - 2; j++) { //from top - if (hilite[3][i][j] > 0.01) { + if (hilite[3][i][j] > epsilon) { hilite_dir[0 + c][i][j] = hilite[c][i][j] / hilite[3][i][j]; } else { hilite_dir[0 + c][i][j] = 0.1 * ((hilite_dir[0 + c][i - 1][j - 2] + hilite_dir[0 + c][i - 1][j - 1] + hilite_dir[0 + c][i - 1][j] + hilite_dir[0 + c][i - 1][j + 1] + hilite_dir[0 + c][i - 1][j + 2]) / - (hilite_dir[0 + 3][i - 1][j - 2] + hilite_dir[0 + 3][i - 1][j - 1] + hilite_dir[0 + 3][i - 1][j] + hilite_dir[0 + 3][i - 1][j + 1] + hilite_dir[0 + 3][i - 1][j + 2] + 0.00001)); + (hilite_dir[0 + 3][i - 1][j - 2] + hilite_dir[0 + 3][i - 1][j - 1] + hilite_dir[0 + 3][i - 1][j] + hilite_dir[0 + 3][i - 1][j + 1] + hilite_dir[0 + 3][i - 1][j + 2] + epsilon)); } } } for (int j = 2; j < hfw - 2; j++) { - if(hilite[3][hfh - 2][j] <= 0.01) { + if(hilite[3][hfh - 2][j] <= epsilon) { hilite_dir[4 + c][hfh - 1][j] += hilite_dir[0 + c][hfh - 2][j]; } } - } if(plistener) { @@ -843,7 +849,7 @@ void RawImageSource :: HLRecovery_inpaint (float** red, float** green, float** b for (int i = hfh - 2; i > 0; i--) for (int j = 2; j < hfw - 2; j++) { //from bottom - if (hilite[3][i][j] > 0.01) { + if (hilite[3][i][j] > epsilon) { hilite_dir[4 + 3][i][j] = 1.f; } else { hilite_dir[4 + 3][i][j] = (hilite_dir[4 + 3][(i + 1)][(j - 2)] + hilite_dir[4 + 3][(i + 1)][(j - 1)] + hilite_dir[4 + 3][(i + 1)][(j)] + hilite_dir[4 + 3][(i + 1)][(j + 1)] + hilite_dir[4 + 3][(i + 1)][(j + 2)]) == 0.f ? 0.f : 0.1f; @@ -858,11 +864,11 @@ void RawImageSource :: HLRecovery_inpaint (float** red, float** green, float** b for (int i = hfh - 2; i > 0; i--) { for (int j = 2; j < hfw - 2; j++) { //from bottom - if (hilite[3][i][j] > 0.01) { + if (hilite[3][i][j] > epsilon) { hilite_dir[4 + c][i][j] = hilite[c][i][j] / hilite[3][i][j]; } else { hilite_dir[4 + c][i][j] = 0.1 * ((hilite_dir[4 + c][(i + 1)][(j - 2)] + hilite_dir[4 + c][(i + 1)][(j - 1)] + hilite_dir[4 + c][(i + 1)][(j)] + hilite_dir[4 + c][(i + 1)][(j + 1)] + hilite_dir[4 + c][(i + 1)][(j + 2)]) / - (hilite_dir[4 + 3][(i + 1)][(j - 2)] + hilite_dir[4 + 3][(i + 1)][(j - 1)] + hilite_dir[4 + 3][(i + 1)][(j)] + hilite_dir[4 + 3][(i + 1)][(j + 1)] + hilite_dir[4 + 3][(i + 1)][(j + 2)] + 0.00001)); + (hilite_dir[4 + 3][(i + 1)][(j - 2)] + hilite_dir[4 + 3][(i + 1)][(j - 1)] + hilite_dir[4 + 3][(i + 1)][(j)] + hilite_dir[4 + 3][(i + 1)][(j + 1)] + hilite_dir[4 + 3][(i + 1)][(j + 2)] + epsilon)); } } } @@ -933,25 +939,18 @@ void RawImageSource :: HLRecovery_inpaint (float** red, float** green, float** b hilite_dir4[c][hfw - 1][hfh - 1] = hilite_dir4[c][hfw - 1][hfh - 2] = hilite_dir4[c][hfw - 2][hfh - 1] = hilite_dir4[c][hfw - 2][hfh - 2] = hilite_dir4[c][hfw - 3][hfh - 3]; } - - if(plistener) { progress += 0.05; plistener->setProgress(progress); } - //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - // now reconstruct clipped channels using color ratios - + //free up some memory for(int c = 0; c < 4; c++) { hilite[c].free(); } - LUTf invfn(0x10000); - - for (int i = 0; i < 0x10000; i++) { - invfn[i] = 1.0 / (1.0 + i); - } + //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + // now reconstruct clipped channels using color ratios #ifdef _OPENMP #pragma omp parallel for schedule(dynamic,16) @@ -961,7 +960,6 @@ void RawImageSource :: HLRecovery_inpaint (float** red, float** green, float** b int i1 = min((i - (i % pitch)) / pitch, hfh - 1); for (int j = 0; j < width; j++) { - int j1 = min((j - (j % pitch)) / pitch, hfw - 1); float pixel[3] = {red[i][j], green[i][j], blue[i][j]}; @@ -969,9 +967,9 @@ void RawImageSource :: HLRecovery_inpaint (float** red, float** green, float** b continue; //pixel not clipped } - //%%%%%%%%%%%%%%%%%%%%%%% - //estimate recovered values using modified HLRecovery_blend algorithm + int j1 = min((j - (j % pitch)) / pitch, hfw - 1); + //estimate recovered values using modified HLRecovery_blend algorithm float rgb[ColorCount], rgb_blend[ColorCount] = {}, cam[2][ColorCount], lab[2][ColorCount], sum[2], chratio; // Copy input pixel to rgb so it's easier to access in loops @@ -980,23 +978,22 @@ void RawImageSource :: HLRecovery_inpaint (float** red, float** green, float** b rgb[2] = pixel[2]; // Initialize cam with raw input [0] and potentially clipped input [1] - FOREACHCOLOR { + for (int c = 0; c < ColorCount; c++) { cam[0][c] = rgb[c]; cam[1][c] = min(cam[0][c], clippt); } // Calculate the lightness correction ratio (chratio) for (int i2 = 0; i2 < 2; i2++) { - FOREACHCOLOR { + for (int c = 0; c < ColorCount; c++) { lab[i2][c] = 0; - for (int j = 0; j < ColorCount; j++) - { - lab[i2][c] += trans[ColorCount - 3][c][j] * cam[i2][j]; + for (int j = 0; j < ColorCount; j++) { + lab[i2][c] += trans[c][j] * cam[i2][j]; } } - sum[i2] = 0; + sum[i2] = 0.f; for (int c = 1; c < ColorCount; c++) { sum[i2] += SQR(lab[i2][c]); @@ -1004,10 +1001,10 @@ void RawImageSource :: HLRecovery_inpaint (float** red, float** green, float** b } if(sum[0] == 0.f) { // avoid division by zero - sum[0] = 0.0001f; + sum[0] = epsilon; } - chratio = sqrt(sum[1] / sum[0]); + chratio = sqrtf(sum[1] / sum[0]); // Apply ratio to lightness in lab space @@ -1016,83 +1013,89 @@ void RawImageSource :: HLRecovery_inpaint (float** red, float** green, float** b } // Transform back from lab to RGB - FOREACHCOLOR { + for (int c = 0; c < ColorCount; c++) { cam[0][c] = 0; - for (int j = 0; j < ColorCount; j++) - { - cam[0][c] += itrans[ColorCount - 3][c][j] * lab[0][j]; + for (int j = 0; j < ColorCount; j++) { + cam[0][c] += itrans[c][j] * lab[0][j]; } } - FOREACHCOLOR rgb[c] = cam[0][c] / ColorCount; + + for (int c = 0; c < ColorCount; c++) { + rgb[c] = cam[0][c] / ColorCount; + } // Copy converted pixel back - float rfrac = max(0.f,min(1.0f, medFactor[0] * (pixel[0] - blendpt))); - float gfrac = max(0.f,min(1.0f, medFactor[1] * (pixel[1] - blendpt))); - float bfrac = max(0.f,min(1.0f, medFactor[2] * (pixel[2] - blendpt))); + float rfrac = max(0.f, min(1.f, medFactor[0] * (pixel[0] - blendpt))); + float gfrac = max(0.f, min(1.f, medFactor[1] * (pixel[1] - blendpt))); + float bfrac = max(0.f, min(1.f, medFactor[2] * (pixel[2] - blendpt))); if (pixel[0] > blendpt) { - rgb_blend[0] = rfrac * rgb[0] + (1 - rfrac) * pixel[0]; + rgb_blend[0] = rfrac * rgb[0] + (1.f - rfrac) * pixel[0]; } if (pixel[1] > blendpt) { - rgb_blend[1] = gfrac * rgb[1] + (1 - gfrac) * pixel[1]; + rgb_blend[1] = gfrac * rgb[1] + (1.f - gfrac) * pixel[1]; } if (pixel[2] > blendpt) { - rgb_blend[2] = bfrac * rgb[2] + (1 - bfrac) * pixel[2]; + rgb_blend[2] = bfrac * rgb[2] + (1.f - bfrac) * pixel[2]; } //end of HLRecovery_blend estimation //%%%%%%%%%%%%%%%%%%%%%%% - //float pixref[3]={min(Yclip,hfsize[0][i1][j1]),min(Yclip,hfsize[1][i1][j1]),min(Yclip,hfsize[2][i1][j1])}; - //there are clipped highlights //first, determine weighted average of unclipped extensions (weighting is by 'hue' proximity) - float dirwt, factor; - float totwt = 0; //0.0003; - float clipfix[3] = {0, 0, 0}; //={totwt*rgb_blend[0],totwt*rgb_blend[1],totwt*rgb_blend[2]}; + float totwt = 0.f; + float clipfix[3] = {0.f, 0.f, 0.f}; - float Yhi = 0.001 + (hilite_dir0[0][j1][i1] + hilite_dir0[1][j1][i1] + hilite_dir0[2][j1][i1]); - float Y = 0.001 + (rgb_blend[0] + rgb_blend[1] + rgb_blend[2]); + float Y = epsilon + rgb_blend[0] + rgb_blend[1] + rgb_blend[2]; - if (hilite_dir0[0][j1][i1] + hilite_dir0[1][j1][i1] + hilite_dir0[2][j1][i1] > 0.5) { - dirwt = invfn[65535 * (SQR(rgb_blend[0] / Y - hilite_dir0[0][j1][i1] / Yhi) + - SQR(rgb_blend[1] / Y - hilite_dir0[1][j1][i1] / Yhi) + - SQR(rgb_blend[2] / Y - hilite_dir0[2][j1][i1] / Yhi))]; - totwt += dirwt; - clipfix[0] += dirwt * hilite_dir0[0][j1][i1] / (hilite_dir0[3][j1][i1] + 0.00001); - clipfix[1] += dirwt * hilite_dir0[1][j1][i1] / (hilite_dir0[3][j1][i1] + 0.00001); - clipfix[2] += dirwt * hilite_dir0[2][j1][i1] / (hilite_dir0[3][j1][i1] + 0.00001); + for (int c = 0; c < ColorCount; c++) { + rgb_blend[c] /= Y; + } + + float Yhi = 1.f / (hilite_dir0[0][j1][i1] + hilite_dir0[1][j1][i1] + hilite_dir0[2][j1][i1]); + + if (Yhi < 2.f) { + float dirwt = 1.f / (1.f + 65535.f * (SQR(rgb_blend[0] - hilite_dir0[0][j1][i1] * Yhi) + + SQR(rgb_blend[1] - hilite_dir0[1][j1][i1] * Yhi) + + SQR(rgb_blend[2] - hilite_dir0[2][j1][i1] * Yhi))); + totwt = dirwt; + dirwt /= (hilite_dir0[3][j1][i1] + epsilon); + clipfix[0] = dirwt * hilite_dir0[0][j1][i1]; + clipfix[1] = dirwt * hilite_dir0[1][j1][i1]; + clipfix[2] = dirwt * hilite_dir0[2][j1][i1]; } for (int dir = 0; dir < 2; dir++) { - float Yhi = 0.001 + (hilite_dir[dir * 4 + 0][i1][j1] + hilite_dir[dir * 4 + 1][i1][j1] + hilite_dir[dir * 4 + 2][i1][j1]); - float Y = 0.001 + (rgb_blend[0] + rgb_blend[1] + rgb_blend[2]); + float Yhi = 1.f / ( hilite_dir[dir * 4 + 0][i1][j1] + hilite_dir[dir * 4 + 1][i1][j1] + hilite_dir[dir * 4 + 2][i1][j1]); - if (hilite_dir[dir * 4 + 0][i1][j1] + hilite_dir[dir * 4 + 1][i1][j1] + hilite_dir[dir * 4 + 2][i1][j1] > 0.5) { - dirwt = invfn[65535 * (SQR(rgb_blend[0] / Y - hilite_dir[dir * 4 + 0][i1][j1] / Yhi) + - SQR(rgb_blend[1] / Y - hilite_dir[dir * 4 + 1][i1][j1] / Yhi) + - SQR(rgb_blend[2] / Y - hilite_dir[dir * 4 + 2][i1][j1] / Yhi))]; + if (Yhi < 2.f) { + float dirwt = 1.f / (1.f + 65535.f * (SQR(rgb_blend[0] - hilite_dir[dir * 4 + 0][i1][j1] * Yhi) + + SQR(rgb_blend[1] - hilite_dir[dir * 4 + 1][i1][j1] * Yhi) + + SQR(rgb_blend[2] - hilite_dir[dir * 4 + 2][i1][j1] * Yhi))); totwt += dirwt; - clipfix[0] += dirwt * hilite_dir[dir * 4 + 0][i1][j1] / (hilite_dir[dir * 4 + 3][i1][j1] + 0.00001); - clipfix[1] += dirwt * hilite_dir[dir * 4 + 1][i1][j1] / (hilite_dir[dir * 4 + 3][i1][j1] + 0.00001); - clipfix[2] += dirwt * hilite_dir[dir * 4 + 2][i1][j1] / (hilite_dir[dir * 4 + 3][i1][j1] + 0.00001); + dirwt /= (hilite_dir[dir * 4 + 3][i1][j1] + epsilon); + clipfix[0] += dirwt * hilite_dir[dir * 4 + 0][i1][j1]; + clipfix[1] += dirwt * hilite_dir[dir * 4 + 1][i1][j1]; + clipfix[2] += dirwt * hilite_dir[dir * 4 + 2][i1][j1]; } } - Yhi = 0.001 + (hilite_dir4[0][j1][i1] + hilite_dir4[1][j1][i1] + hilite_dir4[2][j1][i1]); - Y = 0.001 + (rgb_blend[0] + rgb_blend[1] + rgb_blend[2]); - if (hilite_dir4[0][j1][i1] + hilite_dir4[1][j1][i1] + hilite_dir4[2][j1][i1] > 0.5) { - dirwt = invfn[65535 * (SQR(rgb_blend[0] / Y - hilite_dir4[0][j1][i1] / Yhi) + - SQR(rgb_blend[1] / Y - hilite_dir4[1][j1][i1] / Yhi) + - SQR(rgb_blend[2] / Y - hilite_dir4[2][j1][i1] / Yhi))]; + Yhi = 1.f / (hilite_dir4[0][j1][i1] + hilite_dir4[1][j1][i1] + hilite_dir4[2][j1][i1]); + + if (Yhi < 2.f) { + float dirwt = 1.f / (1.f + 65535.f * (SQR(rgb_blend[0] - hilite_dir4[0][j1][i1] * Yhi) + + SQR(rgb_blend[1] - hilite_dir4[1][j1][i1] * Yhi) + + SQR(rgb_blend[2] - hilite_dir4[2][j1][i1] * Yhi))); totwt += dirwt; - clipfix[0] += dirwt * hilite_dir4[0][j1][i1] / (hilite_dir4[3][j1][i1] + 0.00001); - clipfix[1] += dirwt * hilite_dir4[1][j1][i1] / (hilite_dir4[3][j1][i1] + 0.00001); - clipfix[2] += dirwt * hilite_dir4[2][j1][i1] / (hilite_dir4[3][j1][i1] + 0.00001); + dirwt /= (hilite_dir4[3][j1][i1] + epsilon); + clipfix[0] += dirwt * hilite_dir4[0][j1][i1]; + clipfix[1] += dirwt * hilite_dir4[1][j1][i1]; + clipfix[2] += dirwt * hilite_dir4[2][j1][i1]; } if(totwt == 0.f) { @@ -1102,42 +1105,39 @@ void RawImageSource :: HLRecovery_inpaint (float** red, float** green, float** b clipfix[0] /= totwt; clipfix[1] /= totwt; clipfix[2] /= totwt; - //sumwt += totwt; - //counts ++; //now correct clipped channels if (pixel[0] > max_f[0] && pixel[1] > max_f[1] && pixel[2] > max_f[2]) { //all channels clipped float Y = (0.299 * clipfix[0] + 0.587 * clipfix[1] + 0.114 * clipfix[2]); - factor = whitept / Y; + float factor = whitept / Y; red[i][j] = clipfix[0] * factor; green[i][j] = clipfix[1] * factor; blue[i][j] = clipfix[2] * factor; } else {//some channels clipped - const float eps = 0.0001f; - int notclipped[3] = {pixel[0] < max_f[0] ? 1 : 0, pixel[1] < max_f[1] ? 1 : 0, pixel[2] < max_f[2] ? 1 : 0}; + float notclipped[3] = {pixel[0] <= max_f[0] ? 1.f : 0.f, pixel[1] <= max_f[1] ? 1.f : 0.f, pixel[2] <= max_f[2] ? 1.f : 0.f}; - if (notclipped[0] == 0) { //red clipped + if (notclipped[0] == 0.f) { //red clipped red[i][j] = max(red[i][j], (clipfix[0] * ((notclipped[1] * pixel[1] + notclipped[2] * pixel[2]) / - (notclipped[1] * clipfix[1] + notclipped[2] * clipfix[2] + eps)))); + (notclipped[1] * clipfix[1] + notclipped[2] * clipfix[2] + epsilon)))); } - if (notclipped[1] == 0) { //green clipped + if (notclipped[1] == 0.f) { //green clipped green[i][j] = max(green[i][j], (clipfix[1] * ((notclipped[2] * pixel[2] + notclipped[0] * pixel[0]) / - (notclipped[2] * clipfix[2] + notclipped[0] * clipfix[0] + eps)))); + (notclipped[2] * clipfix[2] + notclipped[0] * clipfix[0] + epsilon)))); } - if (notclipped[2] == 0) { //blue clipped + if (notclipped[2] == 0.f) { //blue clipped blue[i][j] = max(blue[i][j], (clipfix[2] * ((notclipped[0] * pixel[0] + notclipped[1] * pixel[1]) / - (notclipped[0] * clipfix[0] + notclipped[1] * clipfix[1] + eps)))); + (notclipped[0] * clipfix[0] + notclipped[1] * clipfix[1] + epsilon)))); } } Y = (0.299 * red[i][j] + 0.587 * green[i][j] + 0.114 * blue[i][j]); if (Y > whitept) { - factor = whitept / Y; + float factor = whitept / Y; red[i][j] *= factor; green[i][j] *= factor;