diff --git a/README.md b/README.md index baa566ca0..d7ab40a5c 100644 --- a/README.md +++ b/README.md @@ -2,6 +2,19 @@ RawTherapee is a powerful, cross-platform raw photo processing program, released under the GNU General Public License Version 3. It is written in C++ using a GTK+ front-end and a patched version of dcraw for reading raw files. It is notable for the advanced control it gives the user over the demosaicing and developing process. +## Target audience and feature set + +Rawtherapee is a free software (as in free beer) aimed at developing raw files from a broad range of camera (all files than can be decoded by *dcraw* from Dave Coffin), but can also edit standard image files (jpeg, png, tiff) and some HDR files. The targeted user base ranges from the enthusiast photographer with some basic technical knowledge in color science (it's not mandatory, but it's better if you e.g. know what a color profile is) to semi-professional users. + +Of course, professionals can use RawTherapee (for free) too, but will probably lack some peripheral features like Digital Assets Management, printing, uploading, etc. RawTherapee is not and probably will not be an all-in-one software in the foreseeable future. But the Open Source community is sufficiently developed now to offer you all that missing peripheral features. + +The general feature set includes : +* basic file handling (move, delete, review, ranking, tagging) through thumnails +* image editing through 32 bits floating points filters (called *Tools*) that affect the whole image uniformly (by opposition to local editing, which is planned but still on the TODO list) +* export to standard 8 bits image file, with possibility to automatically load the resulting image in you favorite image editor for local adjustment + +## Links + Website: http://rawtherapee.com/ diff --git a/rtengine/LUT.h b/rtengine/LUT.h index b6eb1a05c..965180c92 100644 --- a/rtengine/LUT.h +++ b/rtengine/LUT.h @@ -71,9 +71,7 @@ #include #include #endif -#ifdef __SSE2__ -#include "sleefsseavx.c" -#endif +#include "opthelper.h" #include #include "rt_math.h" @@ -91,10 +89,9 @@ protected: private: unsigned int owner; #if defined( __SSE2__ ) && defined( __x86_64__ ) - __m128 maxsv __attribute__ ((aligned (16))); - __m128 sizev __attribute__ ((aligned (16))); - __m128i maxsiv __attribute__ ((aligned (16))); - __m128i sizeiv __attribute__ ((aligned (16))); + vfloat maxsv __attribute__ ((aligned (16))); + vfloat sizev __attribute__ ((aligned (16))); + vint sizeiv __attribute__ ((aligned (16))); #endif public: /// convenience flag! If one doesn't want to delete the buffer but want to flag it to be recomputed... @@ -120,10 +117,9 @@ public: maxs = size - 2; maxsf = (float)maxs; #if defined( __SSE2__ ) && defined( __x86_64__ ) - maxsv = _mm_set1_ps( maxs ); - maxsiv = _mm_cvttps_epi32( maxsv ); + maxsv = F2V( maxs ); sizeiv = _mm_set1_epi32( (int)(size - 1) ); - sizev = _mm_set1_ps( size - 1 ); + sizev = F2V( size - 1 ); #endif } void operator ()(int s, int flags = 0xfffffff) @@ -150,10 +146,9 @@ public: maxs = size - 2; maxsf = (float)maxs; #if defined( __SSE2__ ) && defined( __x86_64__ ) - maxsv = _mm_set1_ps( maxs ); - maxsiv = _mm_cvttps_epi32( maxsv ); + maxsv = F2V( maxs ); sizeiv = _mm_set1_epi32( (int)(size - 1) ); - sizev = _mm_set1_ps( size - 1 ); + sizev = F2V( size - 1 ); #endif } @@ -167,11 +162,11 @@ public: assert (s > 0); - if (source == NULL) { + if (!source) { printf("source is NULL!\n"); } - assert (source != NULL); + assert (source != nullptr); #endif dirty = false; // Assumption clip = flags; @@ -182,10 +177,9 @@ public: maxs = size - 2; maxsf = (float)maxs; #if defined( __SSE2__ ) && defined( __x86_64__ ) - maxsv = _mm_set1_ps( size - 2); - maxsiv = _mm_cvttps_epi32( maxsv ); + maxsv = F2V( size - 2); sizeiv = _mm_set1_epi32( (int)(size - 1) ); - sizev = _mm_set1_ps( size - 1 ); + sizev = F2V( size - 1 ); #endif for (int i = 0; i < s; i++) { @@ -195,7 +189,7 @@ public: LUT() { - data = NULL; + data = nullptr; reset(); } @@ -237,10 +231,10 @@ public: if (this != &rhs) { if (rhs.size > this->size) { delete [] this->data; - this->data = NULL; + this->data = nullptr; } - if (this->data == NULL) { + if (this->data == nullptr) { this->data = new T[rhs.size]; } @@ -252,10 +246,9 @@ public: this->maxs = this->size - 2; this->maxsf = (float)this->maxs; #if defined( __SSE2__ ) && defined( __x86_64__ ) - this->maxsv = _mm_set1_ps( this->size - 2); - this->maxsiv = _mm_cvttps_epi32( this->maxsv ); + this->maxsv = F2V( this->size - 2); this->sizeiv = _mm_set1_epi32( (int)(this->size - 1) ); - this->sizev = _mm_set1_ps( this->size - 1 ); + this->sizev = F2V( this->size - 1 ); #endif } @@ -268,14 +261,15 @@ public: } #if defined( __SSE2__ ) && defined( __x86_64__ ) - __m128 operator[](__m128 indexv ) const +/* + vfloat operator[](vfloat indexv ) const { // printf("don't use this operator. It's not ready for production"); return _mm_setzero_ps(); // convert floats to ints - __m128i idxv = _mm_cvttps_epi32( indexv ); - __m128 tempv, resultv, p1v, p2v; + vint idxv = _mm_cvttps_epi32( indexv ); + vfloat tempv, resultv, p1v, p2v; vmask maxmask = vmaskf_gt(indexv, maxsv); idxv = _mm_castps_si128(vself(maxmask, maxsv, _mm_castsi128_ps(idxv))); vmask minmask = vmaskf_lt(indexv, _mm_setzero_ps()); @@ -327,15 +321,55 @@ public: p2v = _mm_move_ss( p2v, tempv); // now p2v is 3 2 1 0 - __m128 diffv = indexv - _mm_cvtepi32_ps ( idxv ); + vfloat diffv = indexv - _mm_cvtepi32_ps ( idxv ); diffv = vself(vorm(maxmask, minmask), _mm_setzero_ps(), diffv); resultv = p1v + p2v * diffv; return resultv ; } - - __m128 operator[](__m128i idxv ) const +*/ +#ifdef __SSE4_1__ + vfloat operator[](vint idxv ) const { - __m128 tempv, p1v; + vfloat tempv, p1v; + idxv = _mm_max_epi32( _mm_setzero_si128(), _mm_min_epi32(idxv, sizeiv)); + // access the LUT 4 times and shuffle the values into p1v + + int idx; + + // get 4th value + idx = _mm_extract_epi32(idxv, 3); + tempv = _mm_load_ss(&data[idx]); + p1v = PERMUTEPS(tempv, _MM_SHUFFLE(0, 0, 0, 0)); + // now p1v is 3 3 3 3 + + // get 3rd value + idx = _mm_extract_epi32(idxv, 2); + tempv = _mm_load_ss(&data[idx]); + p1v = _mm_move_ss( p1v, tempv); + // now p1v is 3 3 3 2 + + // get 2nd value + idx = _mm_extract_epi32(idxv, 1); + tempv = _mm_load_ss(&data[idx]); + p1v = PERMUTEPS( p1v, _MM_SHUFFLE(1, 0, 1, 0)); + // now p1v is 3 2 3 2 + p1v = _mm_move_ss( p1v, tempv ); + // now p1v is 3 2 3 1 + + // get 1st value + idx = _mm_cvtsi128_si32(idxv); + tempv = _mm_load_ss(&data[idx]); + p1v = PERMUTEPS( p1v, _MM_SHUFFLE(3, 2, 0, 0)); + // now p1v is 3 2 1 1 + p1v = _mm_move_ss( p1v, tempv ); + // now p1v is 3 2 1 0 + + return p1v; + } +#else + vfloat operator[](vint idxv ) const + { + vfloat tempv, p1v; tempv = _mm_cvtepi32_ps(idxv); tempv = _mm_min_ps( tempv, sizev ); idxv = _mm_cvttps_epi32(_mm_max_ps( tempv, _mm_setzero_ps( ) )); @@ -346,7 +380,7 @@ public: // get 4th value idx = _mm_cvtsi128_si32 (_mm_shuffle_epi32(idxv, _MM_SHUFFLE(3, 3, 3, 3))); tempv = _mm_load_ss(&data[idx]); - p1v = _mm_shuffle_ps(tempv, tempv, _MM_SHUFFLE(0, 0, 0, 0)); + p1v = PERMUTEPS(tempv, _MM_SHUFFLE(0, 0, 0, 0)); // now p1v is 3 3 3 3 // get 3rd value @@ -358,7 +392,7 @@ public: // get 2nd value idx = _mm_cvtsi128_si32 (_mm_shuffle_epi32(idxv, _MM_SHUFFLE(1, 1, 1, 1))); tempv = _mm_load_ss(&data[idx]); - p1v = _mm_shuffle_ps( p1v, p1v, _MM_SHUFFLE(1, 0, 1, 0)); + p1v = PERMUTEPS( p1v, _MM_SHUFFLE(1, 0, 1, 0)); // now p1v is 3 2 3 2 p1v = _mm_move_ss( p1v, tempv ); // now p1v is 3 2 3 1 @@ -366,13 +400,14 @@ public: // get 1st value idx = _mm_cvtsi128_si32 (idxv); tempv = _mm_load_ss(&data[idx]); - p1v = _mm_shuffle_ps( p1v, p1v, _MM_SHUFFLE(3, 2, 0, 0)); + p1v = PERMUTEPS( p1v, _MM_SHUFFLE(3, 2, 0, 0)); // now p1v is 3 2 1 1 p1v = _mm_move_ss( p1v, tempv ); // now p1v is 3 2 1 0 return p1v; } +#endif #endif // use with float indices @@ -465,7 +500,7 @@ public: } dirty = true; - data = NULL; + data = nullptr; owner = 1; size = 0; upperBound = 0; @@ -484,7 +519,7 @@ class HueLUT : public LUTf { public: HueLUT() : LUTf() {} - HueLUT(bool createArray) : LUTf() + explicit HueLUT(bool createArray) : LUTf() { if (createArray) { this->operator () (501, LUT_CLIP_BELOW | LUT_CLIP_ABOVE); diff --git a/rtengine/camconst.json b/rtengine/camconst.json index 1cfbfddf6..05e21576c 100644 --- a/rtengine/camconst.json +++ b/rtengine/camconst.json @@ -927,6 +927,13 @@ Quality X: unknown, ie we knowing to little about the camera properties to know "ranges": { "white": 4040 } }, + { // Quality B + "make_model": "FUJIFILM X-PRO2", + "dcraw_matrix": [ 11434,-4948,-1210,-3746,12042,1903,-666,1479,5235 ], // DNG_v9.4 D65 + "raw_crop": [ 0, 0, 6032, 4032 ], // full raw 6160,4032, Usable 6032,4032 - experimental crop + "ranges": { "white": 16100 } + }, + { // Quality B, Matrix from Adobe's dcp D65 instead of the internal in Leica's DNG "make_model": [ "LEICA SL (Typ 601)", "LEICA Q (Typ 116)" ], "dcraw_matrix": [ 10068,-4043,-1068,-5319,14268,1044,-765,1701,6522 ], // DCP D65 diff --git a/rtengine/dcraw.cc b/rtengine/dcraw.cc index 09cdb8d0f..c0dbd816d 100644 --- a/rtengine/dcraw.cc +++ b/rtengine/dcraw.cc @@ -6682,7 +6682,7 @@ void CLASS parse_fuji (int offset) } else if (tag == 0xc000) { c = order; order = 0x4949; - if ((tag = get4()) > 10000) tag = get4(); + while ((tag = get4()) > 10000); width = tag; height = get4(); order = c; diff --git a/rtengine/dcraw.h b/rtengine/dcraw.h index e79b821b1..856da84d1 100644 --- a/rtengine/dcraw.h +++ b/rtengine/dcraw.h @@ -7,7 +7,7 @@ * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. - * + * * RawTherapee is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the @@ -59,7 +59,7 @@ public: ,RT_blacklevel_from_constant(0) ,RT_matrix_from_constant(0) ,getbithuff(this,ifp,zero_after_ff) - ,ph1_bithuff(this,ifp,order) + ,ph1_bithuff(this,ifp,order) ,pana_bits(ifp,load_flags) { memset(&hbd, 0, sizeof(hbd)); @@ -100,7 +100,7 @@ protected: double gamm[6]; dcrawImage_t image; float bright, threshold, user_mul[4]; - + int half_size, four_color_rgb, document_mode, highlight; int verbose, use_auto_wb, use_camera_wb, use_camera_matrix; int output_color, output_bps, output_tiff, med_passes; @@ -111,7 +111,7 @@ protected: int RT_matrix_from_constant; float cam_mul[4], pre_mul[4], cmatrix[3][4], rgb_cam[3][4]; - + int histogram[4][0x2000]; void (DCraw::*write_thumb)(), (DCraw::*write_fun)(); void (DCraw::*load_raw)(), (DCraw::*thumb_load_raw)(); @@ -163,7 +163,7 @@ protected: int rat[10]; unsigned gps[26]; char desc[512], make[64], model[64], soft[32], date[20], artist[64]; - }; + }; protected: int fcol (int row, int col); @@ -362,7 +362,7 @@ void crop_masked_pixels(); void tiff_get (unsigned base, unsigned *tag, unsigned *type, unsigned *len, unsigned *save); void parse_thumb_note (int base, unsigned toff, unsigned tlen); -int parse_tiff_ifd (int base); +int parse_tiff_ifd (int base); void parse_makernote (int base, int uptag); void get_timestamp (int reversed); void parse_exif (int base); @@ -397,6 +397,19 @@ void identify(); void apply_profile (const char *input, const char *output); void jpeg_thumb() {} // not needed bool dcraw_coeff_overrides(const char make[], const char model[], int iso_speed, short trans[12], int *black_level, int *white_level); +void shiftXtransMatrix( const int offsy, const int offsx) { + char xtransTemp[6][6]; + for(int row = 0;row < 6; row++) { + for(int col = 0;col < 6; col++) { + xtransTemp[row][col] = xtrans[(row+offsy)%6][(col+offsx)%6]; + } + } + for(int row = 0;row < 6; row++) { + for(int col = 0;col < 6; col++) { + xtrans[row][col] = xtransTemp[row][col]; + } + } +} }; diff --git a/rtengine/dcraw.patch b/rtengine/dcraw.patch index ce2d186d8..88a56068c 100644 --- a/rtengine/dcraw.patch +++ b/rtengine/dcraw.patch @@ -1,5 +1,5 @@ ---- dcraw.c 2015-09-08 08:08:11.000000000 +0200 -+++ dcraw.cc 2016-01-08 15:37:02.884467080 +0100 +--- dcraw.c 2016-02-11 22:56:58.043957200 +0100 ++++ dcraw.cc 2016-02-11 23:13:28.708268000 +0100 @@ -1,3 +1,15 @@ +/*RT*/#include +/*RT*/#include @@ -1348,6 +1348,15 @@ strcpy (make, "Phase One"); if (model[0]) return; switch (raw_height) { +@@ -6658,7 +6682,7 @@ + } else if (tag == 0xc000) { + c = order; + order = 0x4949; +- if ((tag = get4()) > 10000) tag = get4(); ++ while ((tag = get4()) > 10000); + width = tag; + height = get4(); + order = c; @@ -6688,7 +6712,11 @@ order = get2(); hlen = get4(); diff --git a/rtengine/demosaic_algos.cc b/rtengine/demosaic_algos.cc index aa1b90278..5f5e07139 100644 --- a/rtengine/demosaic_algos.cc +++ b/rtengine/demosaic_algos.cc @@ -21,6 +21,7 @@ #include "rawimagesource.h" #include "rawimagesource_i.h" +#include "jaggedarray.h" #include "median.h" #include "rawimage.h" #include "mytime.h" @@ -94,40 +95,12 @@ void RawImageSource::eahd_demosaic () // end of cielab preparation - float* rh[3]; - float* gh[4]; - float* bh[3]; - float* rv[3]; - float* gv[4]; - float* bv[3]; - float* lLh[3]; - float* lah[3]; - float* lbh[3]; - float* lLv[3]; - float* lav[3]; - float* lbv[3]; - float* homh[3]; - float* homv[3]; - - for (int i = 0; i < 4; i++) { - gh[i] = new float[W]; - gv[i] = new float[W]; - } - - for (int i = 0; i < 3; i++) { - rh[i] = new float[W]; - bh[i] = new float[W]; - rv[i] = new float[W]; - bv[i] = new float[W]; - lLh[i] = new float[W]; - lah[i] = new float[W]; - lbh[i] = new float[W]; - lLv[i] = new float[W]; - lav[i] = new float[W]; - lbv[i] = new float[W]; - homh[i] = new float[W]; - homv[i] = new float[W]; - } + const JaggedArray + rh (W, 3), gh (W, 4), bh (W, 3), + rv (W, 3), gv (W, 4), bv (W, 3), + lLh (W, 3), lah (W, 3), lbh (W, 3), + lLv (W, 3), lav (W, 3), lbv (W, 3), + homh (W, 3), homv (W, 3); // interpolate first two lines interpolate_row_g (gh[0], gv[0], 0); @@ -311,27 +284,9 @@ void RawImageSource::eahd_demosaic () } } - freeArray2(rh, 3); - freeArray2(gh, 4); - freeArray2(bh, 3); - freeArray2(rv, 3); - freeArray2(gv, 4); - freeArray2(bv, 3); - freeArray2(lLh, 3); - freeArray2(lah, 3); - freeArray2(lbh, 3); - freeArray2(homh, 3); - freeArray2(lLv, 3); - freeArray2(lav, 3); - freeArray2(lbv, 3); - freeArray2(homv, 3); - // Interpolate R and B for (int i = 0; i < H; i++) { - if (i == 0) - // rm, gm, bm must be recovered - //interpolate_row_rb_mul_pp (red, blue, NULL, green[i], green[i+1], i, rm, gm, bm, 0, W, 1); - { + if (i == 0) { interpolate_row_rb_mul_pp (red[i], blue[i], NULL, green[i], green[i + 1], i, 1.0, 1.0, 1.0, 0, W, 1); } else if (i == H - 1) { interpolate_row_rb_mul_pp (red[i], blue[i], green[i - 1], green[i], NULL, i, 1.0, 1.0, 1.0, 0, W, 1); @@ -545,7 +500,7 @@ void RawImageSource::hphd_demosaic () plistener->setProgress (0.0); } - float** hpmap = allocArray< float >(W, H, true); + const JaggedArray hpmap (W, H, true); #ifdef _OPENMP #pragma omp parallel @@ -590,17 +545,13 @@ void RawImageSource::hphd_demosaic () #endif hphd_green (hpmap); - freeArray(hpmap, H);//TODO: seems to cause sigabrt ??? why??? if (plistener) { plistener->setProgress (0.66); } for (int i = 0; i < H; i++) { - if (i == 0) - // rm, gm, bm must be recovered - //interpolate_row_rb_mul_pp (red, blue, NULL, green[i], green[i+1], i, rm, gm, bm, 0, W, 1); - { + if (i == 0) { interpolate_row_rb_mul_pp (red[i], blue[i], NULL, green[i], green[i + 1], i, 1.0, 1.0, 1.0, 0, W, 1); } else if (i == H - 1) { interpolate_row_rb_mul_pp (red[i], blue[i], green[i - 1], green[i], NULL, i, 1.0, 1.0, 1.0, 0, W, 1); diff --git a/rtengine/helpersse2.h b/rtengine/helpersse2.h index 3f2bf6299..3e4365e99 100644 --- a/rtengine/helpersse2.h +++ b/rtengine/helpersse2.h @@ -39,6 +39,11 @@ typedef __m128i vint2; #define STVFU(x,y) _mm_storeu_ps(&x,y) #endif +#if defined(__x86_64__) && defined(__AVX__) +#define PERMUTEPS(a,mask) _mm_permute_ps(a,mask) +#else +#define PERMUTEPS(a,mask) _mm_shuffle_ps(a,a,mask) +#endif static INLINE vfloat LC2VFU(float &a) { diff --git a/rtengine/hlmultipliers.cc b/rtengine/hlmultipliers.cc index 8f7c5ef85..35e1d3e87 100644 --- a/rtengine/hlmultipliers.cc +++ b/rtengine/hlmultipliers.cc @@ -20,6 +20,7 @@ #include #include "rawimagesource.h" #include "rawimagesource_i.h" +#include "jaggedarray.h" #include "../rtgui/options.h" namespace rtengine @@ -321,7 +322,7 @@ void RawImageSource::updateHLRecoveryMap_ColorPropagation () int** rec[3]; for (int i = 0; i < 3; i++) { - rec[i] = allocArray (dw, dh); + rec[i] = allocJaggedArray (dw, dh); } float* reds[HR_SCALE]; @@ -333,10 +334,10 @@ void RawImageSource::updateHLRecoveryMap_ColorPropagation () } if (needhr) { - freeArray(needhr, H); + freeJaggedArray(needhr); } - needhr = allocArray (W, H); + needhr = allocJaggedArray (W, H); for (int i = 0; i < dh; i++) { for (int j = 0; j < HR_SCALE; j++) { @@ -413,14 +414,14 @@ void RawImageSource::updateHLRecoveryMap_ColorPropagation () hlmultipliers (rec, max_3, dh, dw); if (hrmap[0] != NULL) { - freeArray (hrmap[0], dh); - freeArray (hrmap[1], dh); - freeArray (hrmap[2], dh); + freeJaggedArray (hrmap[0]); + freeJaggedArray (hrmap[1]); + freeJaggedArray (hrmap[2]); } - hrmap[0] = allocArray (dw, dh); - hrmap[1] = allocArray (dw, dh); - hrmap[2] = allocArray (dw, dh); + hrmap[0] = allocJaggedArray (dw, dh); + hrmap[1] = allocJaggedArray (dw, dh); + hrmap[2] = allocJaggedArray (dw, dh); for (int i = 0; i < dh; i++) for (int j = 0; j < dw; j++) { @@ -431,9 +432,9 @@ void RawImageSource::updateHLRecoveryMap_ColorPropagation () delete ds; - freeArray (rec[0], dh); - freeArray (rec[1], dh); - freeArray (rec[2], dh); + freeJaggedArray (rec[0]); + freeJaggedArray (rec[1]); + freeJaggedArray (rec[2]); } } diff --git a/rtengine/improcfun.h b/rtengine/improcfun.h index 7a2e28441..996f91c58 100644 --- a/rtengine/improcfun.h +++ b/rtengine/improcfun.h @@ -286,14 +286,13 @@ public: void dirpyrdenoise (LabImage* src);//Emil's pyramid denoise void dirpyrequalizer (LabImage* lab, int scale);//Emil's wavelet - void EPDToneMapResid(float * WavCoeffs_L0, unsigned int Iterates, int skip, struct cont_params cp, int W_L, int H_L, float max0, float min0); + void EPDToneMapResid(float * WavCoeffs_L0, unsigned int Iterates, int skip, struct cont_params& cp, int W_L, int H_L, float max0, float min0); float *CompressDR(float *Source, int skip, struct cont_params &cp, int W_L, int H_L, float Compression, float DetailBoost, float max0, float min0, float ave, float ah, float bh, float al, float bl, float factorx, float *Compressed); void ContrastResid(float * WavCoeffs_L0, unsigned int Iterates, int skip, struct cont_params &cp, int W_L, int H_L, float max0, float min0, float ave, float ah, float bh, float al, float bl, float factorx); float *ContrastDR(float *Source, int skip, struct cont_params &cp, int W_L, int H_L, float Compression, float DetailBoost, float max0, float min0, float ave, float ah, float bh, float al, float bl, float factorx, float *Contrast = NULL); void EPDToneMap(LabImage *lab, unsigned int Iterates = 0, int skip = 1); void EPDToneMapCIE(CieImage *ncie, float a_w, float c_, float w_h, int Wid, int Hei, int begh, int endh, float minQ, float maxQ, unsigned int Iterates = 0, int skip = 1); - // void CAT02 (Imagefloat* baseImg, const ProcParams* params); // pyramid denoise procparams::DirPyrDenoiseParams dnparams; @@ -303,13 +302,6 @@ public: void idirpyr (LabImage* data_coarse, LabImage* data_fine, int level, LUTf &rangefn_L, LUTf & nrwt_l, LUTf & nrwt_ab, int pitch, int scale, const int luma, const int chroma/*, LUTf & Lcurve, LUTf & abcurve*/ ); - // FT denoise - //void RGB_InputTransf(Imagefloat * src, LabImage * dst, const procparams::DirPyrDenoiseParams & dnparams, const procparams::DefringeParams & defringe); - //void RGB_OutputTransf(LabImage * src, Imagefloat * dst, const procparams::DirPyrDenoiseParams & dnparams); - //void output_tile_row (float *Lbloxrow, float ** Lhipassdn, float ** tilemask, int height, int width, int top, int blkrad ); - -// void MSR(LabImage* lab, int width, int height, int skip); - void Tile_calc (int tilesize, int overlap, int kall, int imwidth, int imheight, int &numtiles_W, int &numtiles_H, int &tilewidth, int &tileheight, int &tileWskip, int &tileHskip); void ip_wavelet(LabImage * lab, LabImage * dst, int kall, const procparams::WaveletParams & waparams, const WavCurve & wavCLVCcurve, const WavOpacityCurveRG & waOpacityCurveRG, const WavOpacityCurveBY & waOpacityCurveBY, const WavOpacityCurveW & waOpacityCurveW, const WavOpacityCurveWL & waOpacityCurveWL, LUTf &wavclCurve, bool wavcontlutili, int skip); @@ -327,13 +319,13 @@ public: void ContAllAB (LabImage * lab, int maxlvl, float **varhue, float **varchrom, float ** WavCoeffs_a, float * WavCoeffs_a0, int level, int dir, const WavOpacityCurveW & waOpacityCurveW, struct cont_params &cp, int W_ab, int H_ab, const bool useChannelA); void Evaluate2(wavelet_decomposition &WaveletCoeffs_L, - struct cont_params cp, int ind, float *mean, float *meanN, float *sigma, float *sigmaN, float *MaxP, float *MaxN, float madL[8][3]); - void Eval2 (float ** WavCoeffs_L, int level, struct cont_params cp, + const struct cont_params& cp, int ind, float *mean, float *meanN, float *sigma, float *sigmaN, float *MaxP, float *MaxN, float madL[8][3]); + void Eval2 (float ** WavCoeffs_L, int level, const struct cont_params& cp, int W_L, int H_L, int skip_L, int ind, float *mean, float *meanN, float *sigma, float *sigmaN, float *MaxP, float *MaxN, float *madL); void Aver(float * HH_Coeffs, int datalen, float &averagePlus, float &averageNeg, float &max, float &min); void Sigma(float * HH_Coeffs, int datalen, float averagePlus, float averageNeg, float &sigmaPlus, float &sigmaNeg); - void calckoe(float ** WavCoeffs_LL, struct cont_params cp, float ** koeLi, int level, int dir, int W_L, int H_L, float edd, float *maxkoeLi, float **tmC = NULL); + void calckoe(float ** WavCoeffs_LL, const struct cont_params& cp, float ** koeLi, int level, int dir, int W_L, int H_L, float edd, float *maxkoeLi, float **tmC = NULL); diff --git a/rtengine/iplab2rgb.cc b/rtengine/iplab2rgb.cc index c36e185d0..89deac08f 100644 --- a/rtengine/iplab2rgb.cc +++ b/rtengine/iplab2rgb.cc @@ -30,19 +30,10 @@ namespace rtengine { -#define CLIP01(a) ((a)>0?((a)<1?(a):1):0) - extern const Settings* settings; -const double (*wprof[])[3] = {xyz_sRGB, xyz_adobe, xyz_prophoto, xyz_widegamut, xyz_bruce, xyz_beta, xyz_best}; -const double (*iwprof[])[3] = {sRGB_xyz, adobe_xyz, prophoto_xyz, widegamut_xyz, bruce_xyz, beta_xyz, best_xyz}; -const char* wprofnames[] = {"sRGB", "Adobe RGB", "ProPhoto", "WideGamut", "BruceRGB", "Beta RGB", "BestRGB"}; -const int numprof = 7; - void ImProcFunctions::lab2monitorRgb (LabImage* lab, Image8* image) { - //gamutmap(lab); - if (monitorTransform) { int W = lab->W; @@ -211,18 +202,7 @@ Image8* ImProcFunctions::lab2rgb (LabImage* lab, int cx, int cy, int cw, int ch, } } else { - double rgb_xyz[3][3]; - - for (int i = 0; i < numprof; i++) { - if (profile == wprofnames[i]) { - for (int m = 0; m < 3; m++) - for (int n = 0; n < 3; n++) { - rgb_xyz[m][n] = iwprof[i][m][n]; - } - - break; - } - } + const auto rgb_xyz = iccStore->workingSpaceMatrix (profile); #ifdef _OPENMP #pragma omp parallel for schedule(dynamic,16) if (multiThread) @@ -386,7 +366,7 @@ Image16* ImProcFunctions::lab2rgb16b (LabImage* lab, int cx, int cy, int cw, int Image16* image = new Image16 (cw, ch); float p1, p2, p3, p4, p5, p6; //primaries - //double ga0,ga1,ga2,ga3,ga4,ga5=0.0,ga6=0.0;//gamma parameters + double g_a0, g_a1, g_a2, g_a3, g_a4, g_a5; //gamma parameters double pwr; double ts; @@ -401,15 +381,7 @@ Image16* ImProcFunctions::lab2rgb16b (LabImage* lab, int cx, int cy, int cw, int //primaries for 7 working profiles ==> output profiles // eventually to adapt primaries if RT used special profiles ! - if(profi == "ProPhoto") { - p1 = 0.7347; //Prophoto primaries - p2 = 0.2653; - p3 = 0.1596; - p4 = 0.8404; - p5 = 0.0366; - p6 = 0.0001; - select_temp = 1; - } else if (profi == "WideGamut") { + if (profi == "WideGamut") { p1 = 0.7350; //Widegamut primaries p2 = 0.2650; p3 = 0.1150; @@ -457,6 +429,14 @@ Image16* ImProcFunctions::lab2rgb16b (LabImage* lab, int cx, int cy, int cw, int p5 = 0.1300; p6 = 0.0350; select_temp = 1; + } else { + p1 = 0.7347; //ProPhoto and default primaries + p2 = 0.2653; + p3 = 0.1596; + p4 = 0.8404; + p5 = 0.0366; + p6 = 0.0001; + select_temp = 1; } if (!freegamma) {//if Free gamma not selected diff --git a/rtengine/ipretinex.cc b/rtengine/ipretinex.cc index 085e2683a..090841a9e 100644 --- a/rtengine/ipretinex.cc +++ b/rtengine/ipretinex.cc @@ -45,6 +45,7 @@ #include "rawimagesource.h" #include "improcfun.h" #include "opthelper.h" +#define BENCHMARK #include "StopWatch.h" #define MAX_RETINEX_SCALES 8 @@ -208,6 +209,8 @@ void mean_stddv( float **dst, float &mean, float &stddv, int W_L, int H_L, const void RawImageSource::MSR(float** luminance, float** originalLuminance, float **exLuminance, LUTf & mapcurve, bool &mapcontlutili, int width, int height, RetinexParams deh, const RetinextransmissionCurve & dehatransmissionCurve, float &minCD, float &maxCD, float &mini, float &maxi, float &Tmean, float &Tsigma, float &Tmin, float &Tmax) { + BENCHFUN + if (deh.enabled) {//enabled float mean, stddv, maxtr, mintr; //float mini, delta, maxi; @@ -427,7 +430,7 @@ void RawImageSource::MSR(float** luminance, float** originalLuminance, float **e mapmet = 4; } - double shradius = (double) deh.radius; + const double shradius = mapmet == 4 ? (double) deh.radius : 40.; if(deh.viewMethod == "mask") { viewmet = 1; @@ -477,131 +480,128 @@ void RawImageSource::MSR(float** luminance, float** originalLuminance, float **e pond /= log(elogt); } - auto shmap = mapmet > 1 ? new SHMap (W_L, H_L, true) : nullptr; - - + auto shmap = ((mapmet == 2 || mapmet == 3 || mapmet == 4) && it == 1) ? new SHMap (W_L, H_L, true) : nullptr; float *buffer = new float[W_L * H_L];; + + for ( int scale = scal - 1; scale >= 0; scale-- ) { #ifdef _OPENMP - #pragma omp parallel + #pragma omp parallel #endif - { - for ( int scale = scal - 1; scale >= 0; scale-- ) { - if(scale == scal - 1) { + { + if(scale == scal - 1) + { gaussianBlur (src, out, W_L, H_L, RetinexScales[scale], buffer); } else { // reuse result of last iteration + // out was modified in last iteration => restore it + if((((mapmet == 2 && scale > 1) || mapmet == 3 || mapmet == 4) || (mapmet > 0 && mapcontlutili)) && it == 1) + { +#ifdef _OPENMP + #pragma omp for +#endif + + for (int i = 0; i < H_L; i++) { + for (int j = 0; j < W_L; j++) { + out[i][j] = buffer[i * W_L + j]; + } + } + } + gaussianBlur (out, out, W_L, H_L, sqrtf(SQR(RetinexScales[scale]) - SQR(RetinexScales[scale + 1])), buffer); } + if((((mapmet == 2 && scale > 2) || mapmet == 3 || mapmet == 4) || (mapmet > 0 && mapcontlutili)) && it == 1 && scale > 0) + { + // out will be modified => store it for use in next iteration. We even don't need a new buffer because 'buffer' is free after gaussianBlur :) +#ifdef _OPENMP + #pragma omp for +#endif - if(mapmet == 4) { - shradius /= 1.; - } else { - shradius = 40.; + for (int i = 0; i < H_L; i++) { + for (int j = 0; j < W_L; j++) { + buffer[i * W_L + j] = out[i][j]; + } + } } + } - //if(shHighlights > 0 || shShadows > 0) { - if(mapmet == 3) if(it == 1) { - shmap->updateL (out, shradius, true, 1); //wav Total - } - - if(mapmet == 2 && scale > 2) if(it == 1) { - shmap->updateL (out, shradius, true, 1); //wav partial - } - - if(mapmet == 4) if(it == 1) { - shmap->updateL (out, shradius, false, 1); //gauss - } - - //} - if (shmap) { - h_th = shmap->max_f - deh.htonalwidth * (shmap->max_f - shmap->avg) / 100; - s_th = deh.stonalwidth * (shmap->avg - shmap->min_f) / 100; - } + if(((mapmet == 2 && scale > 2) || mapmet == 3 || mapmet == 4) && it == 1) { + shmap->updateL (out, shradius, true, 1); + h_th = shmap->max_f - deh.htonalwidth * (shmap->max_f - shmap->avg) / 100; + s_th = deh.stonalwidth * (shmap->avg - shmap->min_f) / 100; + } #ifdef __SSE2__ - vfloat pondv = F2V(pond); - vfloat limMinv = F2V(ilimdx); - vfloat limMaxv = F2V(limdx); + vfloat pondv = F2V(pond); + vfloat limMinv = F2V(ilimdx); + vfloat limMaxv = F2V(limdx); #endif - if(mapmet > 0) { + if(mapmet > 0 && mapcontlutili && it == 1) { #ifdef _OPENMP - #pragma omp for -#endif - - for (int i = 0; i < H_L; i++) { - if(mapcontlutili) { - int j = 0; - - for (; j < W_L; j++) { - if(it == 1) { - out[i][j] = mapcurve[2.f * out[i][j]] / 2.f; - } - } - } - } - - } - - //if(shHighlights > 0 || shShadows > 0) { - if(((mapmet == 2 && scale > 2) || mapmet == 3 || mapmet == 4) && it == 1) { - - -#ifdef _OPENMP - #pragma omp for -#endif - - for (int i = 0; i < H_L; i++) { - int j = 0; - - for (; j < W_L; j++) { - double mapval = 1.0 + shmap->map[i][j]; - double factor = 1.0; - - if (mapval > h_th) { - factor = (h_th + (100.0 - shHighlights) * (mapval - h_th) / 100.0) / mapval; - } else if (mapval < s_th) { - factor = (s_th - (100.0 - shShadows) * (s_th - mapval) / 100.0) / mapval; - } - - out[i][j] *= factor; - - } - } - } - - //} - -#ifdef _OPENMP - #pragma omp for + #pragma omp parallel for #endif for (int i = 0; i < H_L; i++) { - int j = 0; + for (int j = 0; j < W_L; j++) { + out[i][j] = mapcurve[2.f * out[i][j]] / 2.f; + } + } + + } + + if(((mapmet == 2 && scale > 2) || mapmet == 3 || mapmet == 4) && it == 1) { + + +#ifdef _OPENMP + #pragma omp parallel for +#endif + + for (int i = 0; i < H_L; i++) { + for (int j = 0; j < W_L; j++) { + double mapval = 1.0 + shmap->map[i][j]; + double factor = 1.0; + + if (mapval > h_th) { + factor = (h_th + (100.0 - shHighlights) * (mapval - h_th) / 100.0) / mapval; + } else if (mapval < s_th) { + factor = (s_th - (100.0 - shShadows) * (s_th - mapval) / 100.0) / mapval; + } + + out[i][j] *= factor; + + } + } + } + +#ifdef _OPENMP + #pragma omp parallel for +#endif + + for (int i = 0; i < H_L; i++) { + int j = 0; #ifdef __SSE2__ - if(useHslLin) { - for (; j < W_L - 3; j += 4) { - _mm_storeu_ps(&luminance[i][j], LVFU(luminance[i][j]) + pondv * (LIMV(LVFU(src[i][j]) / LVFU(out[i][j]), limMinv, limMaxv) )); - } - } else { - for (; j < W_L - 3; j += 4) { - _mm_storeu_ps(&luminance[i][j], LVFU(luminance[i][j]) + pondv * xlogf(LIMV(LVFU(src[i][j]) / LVFU(out[i][j]), limMinv, limMaxv) )); - } + if(useHslLin) { + for (; j < W_L - 3; j += 4) { + _mm_storeu_ps(&luminance[i][j], LVFU(luminance[i][j]) + pondv * (LIMV(LVFU(src[i][j]) / LVFU(out[i][j]), limMinv, limMaxv) )); } + } else { + for (; j < W_L - 3; j += 4) { + _mm_storeu_ps(&luminance[i][j], LVFU(luminance[i][j]) + pondv * xlogf(LIMV(LVFU(src[i][j]) / LVFU(out[i][j]), limMinv, limMaxv) )); + } + } #endif - if(useHslLin) { - for (; j < W_L; j++) { - luminance[i][j] += pond * (LIM(src[i][j] / out[i][j], ilimdx, limdx)); - } - } else { - for (; j < W_L; j++) { - luminance[i][j] += pond * xlogf(LIM(src[i][j] / out[i][j], ilimdx, limdx)); // /logt ? - } + if(useHslLin) { + for (; j < W_L; j++) { + luminance[i][j] += pond * (LIM(src[i][j] / out[i][j], ilimdx, limdx)); + } + } else { + for (; j < W_L; j++) { + luminance[i][j] += pond * xlogf(LIM(src[i][j] / out[i][j], ilimdx, limdx)); // /logt ? } } } diff --git a/rtengine/ipwavelet.cc b/rtengine/ipwavelet.cc index 8680a17dc..4e433b453 100644 --- a/rtengine/ipwavelet.cc +++ b/rtengine/ipwavelet.cc @@ -23,9 +23,9 @@ // //////////////////////////////////////////////////////////////// +#include +#include - -#include #include "../rtgui/threadutils.h" #include "rtengine.h" @@ -163,9 +163,6 @@ SSEFUNCTION void ImProcFunctions::ip_wavelet(LabImage * lab, LabImage * dst, int { - MyTime t1e, t2e ; - t1e.set(); - #ifdef _DEBUG // init variables to display Munsell corrections MunsellDebugInfo* MunsDebugInfo = new MunsellDebugInfo(); @@ -219,10 +216,6 @@ SSEFUNCTION void ImProcFunctions::ip_wavelet(LabImage * lab, LabImage * dst, int cp.tonemap = true; } - /*else if(params->wavelet.TMmethod=="std") {cp.TMmeth=1;cp.tonemap=true;} - else if(params->wavelet.TMmethod=="high") {cp.TMmeth=2;cp.tonemap=true;} - else if(params->wavelet.TMmethod=="lowhigh") {cp.TMmeth=3;cp.tonemap=true;} - */ if(params->wavelet.TMmethod == "cont") { cp.contmet = 1; } else if(params->wavelet.TMmethod == "tm") { @@ -1314,15 +1307,7 @@ SSEFUNCTION void ImProcFunctions::ip_wavelet(LabImage * lab, LabImage * dst, int delete dsttmp; } - if (settings->verbose) { - t2e.set(); - printf("Wavelet performed in %d usec:\n", t2e.etime(t1e)); - } - -}//end o - - - +} #undef TS #undef fTS @@ -1427,7 +1412,7 @@ void ImProcFunctions::Sigma( float * RESTRICT DataList, int datalen, float aver } void ImProcFunctions::Evaluate2(wavelet_decomposition &WaveletCoeffs_L, - struct cont_params cp, int ind, float *mean, float *meanN, float *sigma, float *sigmaN, float *MaxP, float *MaxN, float madL[8][3]) + const struct cont_params& cp, int ind, float *mean, float *meanN, float *sigma, float *sigmaN, float *MaxP, float *MaxN, float madL[8][3]) { //StopWatch Stop1("Evaluate2"); int maxlvl = WaveletCoeffs_L.maxlevel(); @@ -1445,7 +1430,7 @@ void ImProcFunctions::Evaluate2(wavelet_decomposition &WaveletCoeffs_L, } } -void ImProcFunctions::Eval2 (float ** WavCoeffs_L, int level, struct cont_params cp, +void ImProcFunctions::Eval2 (float ** WavCoeffs_L, int level, const struct cont_params& cp, int W_L, int H_L, int skip_L, int ind, float *mean, float *meanN, float *sigma, float *sigmaN, float *MaxP, float *MaxN, float *madL) { const float eps = 0.01f; @@ -1460,7 +1445,6 @@ void ImProcFunctions::Eval2 (float ** WavCoeffs_L, int level, struct cont_param for (int dir = 1; dir < 4; dir++) { Aver(WavCoeffs_L[dir], W_L * H_L, avLP[dir], avLN[dir], maxL[dir], minL[dir]); Sigma(WavCoeffs_L[dir], W_L * H_L, avLP[dir], avLN[dir], sigP[dir], sigN[dir]); - // printf("dir=%d level=%d avLP=%f max=%f avLN=%f min=%f sigP=%f sigN=%f\n",dir,level,avLP[dir] ,maxL[dir], avLN[dir] ,minL[dir], sigP[dir], sigN[dir]); } AvL = 0.f; @@ -1501,9 +1485,6 @@ void ImProcFunctions::Eval2 (float ** WavCoeffs_L, int level, struct cont_param sigmaN[level] = SN; MaxP[level] = maxLP; MaxN[level] = maxLN; - // if(params->wavelet.CLmethod!="all") {//display only if user choose different from all - // printf("Ind=%d Level=%d MadL=%f AvL=%f AvN=%f SL=%f SN=%f maxP=%f maxN=%f\n",ind, level,MADL,mean[level],meanN[level],sigma[level],sigmaN[level],MaxP[level],MaxN[level]); - // } } float *ImProcFunctions::ContrastDR(float *Source, int skip, struct cont_params &cp, int W_L, int H_L, float Compression, float DetailBoost, float max0, float min0, float ave, float ah, float bh, float al, float bl, float factorx, float *Contrast) @@ -1520,27 +1501,7 @@ float *ImProcFunctions::ContrastDR(float *Source, int skip, struct cont_params & #endif for (int i = 0; i < W_L * H_L; i++) { //contrast - /* - //source between 0 and 1 - if(Source[i] < 1.f) { - float prov; - if( 32768.f*Source[i]> ave) { - float kh = ah*(Source[i]*100.f)+bh; - prov=32768.f*Source[i]; - Contrast[i]=ave+kh*(Source[i]*32768.f-ave); - } else { - float kl = al*(Source[i]*100.f)+bl; - prov=32768.f*Source[i]; - Contrast[i]=ave-kl*(ave-Source[i]*32768.f); - } - float diflc=Contrast[i]-prov; - diflc*=factorx; - Contrast[i] = (prov + diflc)/32768.f; - //contrast between 0 and 1 - } - */ Contrast[i] = Source[i] ; - } return Contrast; @@ -1608,9 +1569,6 @@ SSEFUNCTION float *ImProcFunctions::CompressDR(float *Source, int skip, struct c temp = (Compression - 1.0f) / 20.f; } -// printf("temp=%f \n", temp); - - #ifdef __SSE2__ #ifdef _RT_NESTED_OPENMP #pragma omp parallel @@ -1708,7 +1666,7 @@ void ImProcFunctions::ContrastResid(float * WavCoeffs_L0, unsigned int Iterates -void ImProcFunctions::EPDToneMapResid(float * WavCoeffs_L0, unsigned int Iterates, int skip, struct cont_params cp, int W_L, int H_L, float max0, float min0) +void ImProcFunctions::EPDToneMapResid(float * WavCoeffs_L0, unsigned int Iterates, int skip, struct cont_params& cp, int W_L, int H_L, float max0, float min0) { @@ -2031,16 +1989,10 @@ void ImProcFunctions::WaveletcontAllL(LabImage * labco, float ** varhue, float * float interm = 0.f; if(cp.lip3 && cp.lipp) { - // comparaison betwwen pixel and neighbours - float kneigh, somm; - - if(cp.neigh == 0) { - kneigh = 38.f; - somm = 50.f; - } else if(cp.neigh == 1) { - kneigh = 28.f; - somm = 40.f; - } + // comparaison between pixel and neighbours + const auto neigh = cp.neigh == 1; + const auto kneigh = neigh ? 28.f : 38.f; + const auto somm = neigh ? 40.f : 50.f; for (int dir = 1; dir < 4; dir++) { //neighbours proxi koeLi[lvl * 3 + dir - 1][i * W_L + j] = (kneigh * koeLi[lvl * 3 + dir - 1][i * W_L + j] + 2.f * koeLi[lvl * 3 + dir - 1][(i - 1) * W_L + j] + 2.f * koeLi[lvl * 3 + dir - 1][(i + 1) * W_L + j] @@ -2372,7 +2324,7 @@ void ImProcFunctions::WaveletcontAllAB(LabImage * labco, float ** varhue, float //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -void ImProcFunctions::calckoe(float ** WavCoeffs_LL, struct cont_params cp, float *koeLi[12], int level, int dir, int W_L, int H_L, float edd, float *maxkoeLi, float **tmC) +void ImProcFunctions::calckoe(float ** WavCoeffs_LL, const struct cont_params& cp, float *koeLi[12], int level, int dir, int W_L, int H_L, float edd, float *maxkoeLi, float **tmC) { int borderL = 2; @@ -2688,6 +2640,8 @@ void ImProcFunctions::finalContAllL (float ** WavCoeffs_L, float * WavCoeffs_L0, void ImProcFunctions::ContAllL (float *koeLi[12], float *maxkoeLi, bool lipschitz, int maxlvl, LabImage * labco, float ** varhue, float **varchrom, float ** WavCoeffs_L, float * WavCoeffs_L0, int level, int dir, struct cont_params &cp, int W_L, int H_L, int skip, float *mean, float *meanN, float *sigma, float *sigmaN, float *MaxP, float *MaxN, const WavCurve & wavCLVCcurve, const WavOpacityCurveW & waOpacityCurveW, FlatCurve* ChCurve, bool Chutili) { + assert (level >= 0); + assert (maxlvl > level); static const float scales[10] = {1.f, 2.f, 4.f, 8.f, 16.f, 32.f, 64.f, 128.f, 256.f, 512.f}; float scaleskip[10]; @@ -2700,17 +2654,11 @@ void ImProcFunctions::ContAllL (float *koeLi[12], float *maxkoeLi, bool lipschit float t_r = settings->top_right; float t_l = settings->top_left; float b_r = settings->bot_right; - float b_l = settings->bot_left; float edd = settings->ed_detec; - float eddlow = settings->ed_low; float eddstrength = settings->ed_detecStr; float aedstr = (eddstrength - 1.f) / 90.f; float bedstr = 1.f - 10.f * aedstr; - bool refi = false; - - // if(cp.lev0s > 0.f || cp.lev1s > 0.f || cp.lev2s > 0.f) refi=true; - // if(cp.val > 0 || refi) {//edge if(cp.val > 0 && cp.edgeena) { float * koe; float maxkoe = 0.f; @@ -2724,7 +2672,7 @@ void ImProcFunctions::ContAllL (float *koeLi[12], float *maxkoeLi, bool lipschit maxkoe = 0.f; - if(cp.detectedge) {// + if(cp.detectedge) { float** tmC; int borderL = 1; tmC = new float*[H_L]; @@ -2785,8 +2733,6 @@ void ImProcFunctions::ContAllL (float *koeLi[12], float *maxkoeLi, bool lipschit } } - //printf("maxkoe=%f \n",maxkoe); - for (int i = 0; i < H_L; i++) { delete [] tmC[i]; } @@ -2846,9 +2792,6 @@ void ImProcFunctions::ContAllL (float *koeLi[12], float *maxkoeLi, bool lipschit } } - float coefsd = 0.85f; //arbitray value to reduce effect after sigma in all case - float coefmean = 0.85f; //arbitray value to reduce effect after sigma in all case -// edge = 1.f + value * exp (expkoef);//estimate edge "pseudo variance" //take into account local contrast float refin = value * exp (expkoef); @@ -2874,7 +2817,6 @@ void ImProcFunctions::ContAllL (float *koeLi[12], float *maxkoeLi, bool lipschit float edgePrecalc = 1.f + refin; //estimate edge "pseudo variance" - //bool exa=false; if(cp.EDmet == 2) { //curve // if(exa) {//curve float insigma = 0.666f; //SD @@ -2889,7 +2831,6 @@ void ImProcFunctions::ContAllL (float *koeLi[12], float *maxkoeLi, bool lipschit float absciss; float kinterm; float kmul; - // for (int i=0; iL[ii * 2][jj * 2]; @@ -3209,18 +3150,6 @@ void ImProcFunctions::ContAllL (float *koeLi[12], float *maxkoeLi, bool lipschit if(useChromAndHue) { - /* - int ii=i/W_L; - int jj=i-ii*W_L; - float LL = labco->L[ii*2][jj*2]; - LL100=LL100init=LL/327.68f; - LL100res=WavCoeffs_L0[i]/327.68f; - float delta=fabs(LL100init-LL100res)/(maxlvl/2); - for(int ml=0;ml + * Copyright (c) 2016 Adam Reichold + * + * RawTherapee is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * RawTherapee is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with RawTherapee. If not, see . + */ +#ifndef JAGGEDARRAY_H +#define JAGGEDARRAY_H + +namespace rtengine +{ + +// These emulate a jagged array, but use only 2 allocations instead of 1 + H. + +template +inline T** const allocJaggedArray (const int W, const int H, const bool initZero = false) +{ + T** const a = new T*[H]; + a[0] = new T[H * W]; + + for (int i = 1; i < H; ++i) { + a[i] = a[i - 1] + W; + } + + if (initZero) { + std::memset(a[0], 0, sizeof(T) * W * H); + } + + return a; +} + +template +inline void freeJaggedArray (T** const a) +{ + delete [] a[0]; + delete [] a; +} + +template +class JaggedArray +{ +public: + JaggedArray (const int W, const int H, const bool initZero = false) + { + a = allocJaggedArray (W, H, initZero); + } + ~JaggedArray () + { + if (a) { + freeJaggedArray (a); + a = nullptr; + } + } + + JaggedArray (const JaggedArray&) = delete; + JaggedArray& operator= (const JaggedArray&) = delete; + +public: + operator T** const () const + { + return a; + } + +private: + T** a; + +}; + +// Declared but not defined to prevent +// explicitly freeing a JaggedArray implicitly cast to T**. +template +void freeJaggedArray (JaggedArray&); + +} // rtengine + +#endif // JAGGEDARRAY_H diff --git a/rtengine/procparams.cc b/rtengine/procparams.cc index 90eafe375..f79f51c91 100644 --- a/rtengine/procparams.cc +++ b/rtengine/procparams.cc @@ -2609,17 +2609,21 @@ int ProcParams::save (Glib::ustring fname, Glib::ustring fname2, bool fnameAbsol if (!pedited || pedited->icm.outputIntent) { Glib::ustring intent; + switch (icm.outputIntent) { default: case RI_PERCEPTUAL: intent = "Perceptual"; break; + case RI_RELATIVE: intent = "Relative"; break; + case RI_SATURATION: intent = "Saturation"; break; + case RI_ABSOLUTE: intent = "Absolute"; break; @@ -4110,7 +4114,7 @@ int ProcParams::load (Glib::ustring fname, ParamsEdited* pedited) if (keyFile.has_key ("Retinex", "Radius")) { - sh.radius = keyFile.get_integer ("Retinex", "Radius"); + retinex.radius = keyFile.get_integer ("Retinex", "Radius"); if (pedited) { pedited->retinex.radius = true; @@ -5845,6 +5849,7 @@ int ProcParams::load (Glib::ustring fname, ParamsEdited* pedited) if (keyFile.has_key ("Color Management", "OutputProfileIntent")) { Glib::ustring intent = keyFile.get_string ("Color Management", "OutputProfileIntent"); + if (intent == "Perceptual") { icm.outputIntent = RI_PERCEPTUAL; } else if (intent == "Relative") { diff --git a/rtengine/rawimage.cc b/rtengine/rawimage.cc index 3b491eab5..89731992a 100644 --- a/rtengine/rawimage.cc +++ b/rtengine/rawimage.cc @@ -495,11 +495,13 @@ int RawImage::loadRaw (bool loadData, bool closeFile, ProgressListener *plistene if (cc && cc->has_rawCrop()) { int lm, tm, w, h; cc->get_rawCrop(lm, tm, w, h); - - if(((int)top_margin - tm) & 1) { // we have an odd border difference - filters = (filters << 4) | (filters >> 28); // left rotate filters by 4 bits + if(isXtrans()) { + shiftXtransMatrix(6 - ((top_margin - tm)%6), 6 - ((left_margin - lm)%6)); + } else { + if(((int)top_margin - tm) & 1) { // we have an odd border difference + filters = (filters << 4) | (filters >> 28); // left rotate filters by 4 bits + } } - left_margin = lm; top_margin = tm; diff --git a/rtengine/rawimagesource.cc b/rtengine/rawimagesource.cc index 2a5f42d56..d1d4da54f 100644 --- a/rtengine/rawimagesource.cc +++ b/rtengine/rawimagesource.cc @@ -22,16 +22,14 @@ #include "rtengine.h" #include "rawimagesource.h" #include "rawimagesource_i.h" +#include "jaggedarray.h" #include "median.h" #include "rawimage.h" #include "mytime.h" -#include "iccmatrices.h" #include "iccstore.h" -#include "image8.h" #include "curves.h" #include "dfmanager.h" #include "ffmanager.h" -#include "../rtgui/options.h" #include "dcp.h" #include "rt_math.h" #include "improcfun.h" @@ -39,11 +37,378 @@ #include #endif #include "opthelper.h" -#include "StopWatch.h" #define clipretinex( val, minv, maxv ) (( val = (val < minv ? minv : val ) ) > maxv ? maxv : val ) #undef CLIPD #define CLIPD(a) ((a)>0.0f?((a)<1.0f?(a):1.0f):0.0f) +namespace +{ + +void rotateLine (const float* const line, rtengine::PlanarPtr &channel, const int tran, const int i, const int w, const int h) +{ + switch(tran & TR_ROT) { + case TR_R180: + for (int j = 0; j < w; j++) { + channel(h - 1 - i, w - 1 - j) = line[j]; + } + + break; + + case TR_R90: + for (int j = 0; j < w; j++) { + channel(j, h - 1 - i) = line[j]; + } + + break; + + case TR_R270: + for (int j = 0; j < w; j++) { + channel(w - 1 - j, i) = line[j]; + } + + break; + + case TR_NONE: + default: + for (int j = 0; j < w; j++) { + channel(i, j) = line[j]; + } + } +} + +void transLineStandard (const float* const red, const float* const green, const float* const blue, const int i, rtengine::Imagefloat* const image, const int tran, const int imwidth, const int imheight) +{ + // conventional CCD coarse rotation + rotateLine (red, image->r, tran, i, imwidth, imheight); + rotateLine (green, image->g, tran, i, imwidth, imheight); + rotateLine (blue, image->b, tran, i, imwidth, imheight); +} + +void transLineFuji (const float* const red, const float* const green, const float* const blue, const int i, rtengine::Imagefloat* const image, const int tran, const int imwidth, const int imheight, const int fw) +{ + + // Fuji SuperCCD rotation + coarse rotation + int start = ABS(fw - i); + int w = fw * 2 + 1; + int h = (imheight - fw) * 2 + 1; + int end = min(h + fw - i, w - fw + i); + + switch(tran & TR_ROT) { + case TR_R180: + for (int j = start; j < end; j++) { + int y = i + j - fw; + int x = fw - i + j; + + if (x >= 0 && y < image->height && y >= 0 && x < image->width) { + image->r(image->height - 1 - y, image->width - 1 - x) = red[j]; + image->g(image->height - 1 - y, image->width - 1 - x) = green[j]; + image->b(image->height - 1 - y, image->width - 1 - x) = blue[j]; + } + } + + break; + + case TR_R270: + for (int j = start; j < end; j++) { + int y = i + j - fw; + int x = fw - i + j; + + if (x >= 0 && x < image->height && y >= 0 && y < image->width) { + image->r(image->height - 1 - x, y) = red[j]; + image->g(image->height - 1 - x, y) = green[j]; + image->b(image->height - 1 - x, y) = blue[j]; + } + } + + break; + + case TR_R90: + for (int j = start; j < end; j++) { + int y = i + j - fw; + int x = fw - i + j; + + if (x >= 0 && y < image->width && y >= 0 && x < image->height) { + image->r(x, image->width - 1 - y) = red[j]; + image->g(x, image->width - 1 - y) = green[j]; + image->b(x, image->width - 1 - y) = blue[j]; + } + } + + break; + + case TR_NONE: + default: + for (int j = start; j < end; j++) { + int y = i + j - fw; + int x = fw - i + j; + + if (x >= 0 && y < image->height && y >= 0 && x < image->width) { + image->r(y, x) = red[j]; + image->g(y, x) = green[j]; + image->b(y, x) = blue[j]; + } + } + } +} + +void transLineD1x (const float* const red, const float* const green, const float* const blue, const int i, rtengine::Imagefloat* const image, const int tran, const int imwidth, const int imheight, const bool oddHeight, const bool clip) +{ + // Nikon D1X has an uncommon sensor with 4028 x 1324 sensels. + // Vertical sensel size is 2x horizontal sensel size + // We have to do vertical interpolation for the 'missing' rows + // We do that in combination with coarse rotation + + switch(tran & TR_ROT) { + case TR_R180: // rotate 180 degree + for (int j = 0; j < imwidth; j++) { + image->r(2 * (imheight - 1 - i), imwidth - 1 - j) = red[j]; + image->g(2 * (imheight - 1 - i), imwidth - 1 - j) = green[j]; + image->b(2 * (imheight - 1 - i), imwidth - 1 - j) = blue[j]; + } + + if (i == 0) { + for (int j = 0; j < imwidth; j++) { + image->r(2 * imheight - 1, imwidth - 1 - j) = red[j]; + image->g(2 * imheight - 1, imwidth - 1 - j) = green[j]; + image->b(2 * imheight - 1, imwidth - 1 - j) = blue[j]; + } + } + + if (i == 1 || i == 2) { // linear interpolation + int row = 2 * imheight - 1 - 2 * i; + + for (int j = 0; j < imwidth; j++) { + int col = imwidth - 1 - j; + image->r(row, col) = (red[j] + image->r(row + 1, col)) / 2; + image->g(row, col) = (green[j] + image->g(row + 1, col)) / 2; + image->b(row, col) = (blue[j] + image->b(row + 1, col)) / 2; + } + + if(i == 2 && oddHeight) { + int row = 2 * imheight; + + for (int j = 0; j < imwidth; j++) { + int col = imwidth - 1 - j; + image->r(row, col) = (red[j] + image->r(row - 2, col)) / 2; + image->g(row, col) = (green[j] + image->g(row - 2, col)) / 2; + image->b(row, col) = (blue[j] + image->b(row - 2, col)) / 2; + } + } + } else if (i == imheight - 1 || i == imheight - 2) { + int row = 2 * imheight - 1 - 2 * i; + + for (int j = 0; j < imwidth; j++) { + int col = imwidth - 1 - j; + image->r(row, col) = (red[j] + image->r(row + 1, col)) / 2; + image->g(row, col) = (green[j] + image->g(row + 1, col)) / 2; + image->b(row, col) = (blue[j] + image->b(row + 1, col)) / 2; + } + + row = 2 * imheight - 1 - 2 * i + 2; + + for (int j = 0; j < imwidth; j++) { + int col = imwidth - 1 - j; + image->r(row, col) = (red[j] + image->r(row + 1, col)) / 2; + image->g(row, col) = (green[j] + image->g(row + 1, col)) / 2; + image->b(row, col) = (blue[j] + image->b(row + 1, col)) / 2; + } + } else if (i > 2 && i < imheight - 1) { // vertical bicubic interpolation + int row = 2 * imheight - 1 - 2 * i + 2; + + for (int j = 0; j < imwidth; j++) { + int col = imwidth - 1 - j; + image->r(row, col) = MAX(0.f, -0.0625f * (red[j] + image->r(row + 3, col)) + 0.5625f * (image->r(row - 1, col) + image->r(row + 1, col))); + image->g(row, col) = MAX(0.f, -0.0625f * (green[j] + image->g(row + 3, col)) + 0.5625f * (image->g(row - 1, col) + image->g(row + 1, col))); + image->b(row, col) = MAX(0.f, -0.0625f * (blue[j] + image->b(row + 3, col)) + 0.5625f * (image->b(row - 1, col) + image->b(row + 1, col))); + + if(clip) { + image->r(row, col) = MIN(image->r(row, col), rtengine::MAXVALF); + image->g(row, col) = MIN(image->g(row, col), rtengine::MAXVALF); + image->b(row, col) = MIN(image->b(row, col), rtengine::MAXVALF); + } + } + } + + break; + + case TR_R90: // rotate right + if( i == 0) { + for (int j = 0; j < imwidth; j++) { + image->r(j, 2 * imheight - 1) = red[j]; + image->g(j, 2 * imheight - 1) = green[j]; + image->b(j, 2 * imheight - 1) = blue[j]; + } + } + + for (int j = 0; j < imwidth; j++) { + image->r(j, 2 * (imheight - 1 - i)) = red[j]; + image->g(j, 2 * (imheight - 1 - i)) = green[j]; + image->b(j, 2 * (imheight - 1 - i)) = blue[j]; + } + + if (i == 1 || i == 2) { // linear interpolation + int col = 2 * imheight - 1 - 2 * i; + + for (int j = 0; j < imwidth; j++) { + image->r(j, col) = (red[j] + image->r(j, col + 1)) / 2; + image->g(j, col) = (green[j] + image->g(j, col + 1)) / 2; + image->b(j, col) = (blue[j] + image->b(j, col + 1)) / 2; + + if(oddHeight && i == 2) { + image->r(j, 2 * imheight) = (red[j] + image->r(j, 2 * imheight - 2)) / 2; + image->g(j, 2 * imheight) = (green[j] + image->g(j, 2 * imheight - 2)) / 2; + image->b(j, 2 * imheight) = (blue[j] + image->b(j, 2 * imheight - 2)) / 2; + } + } + } else if (i == imheight - 1) { + int col = 2 * imheight - 1 - 2 * i; + + for (int j = 0; j < imwidth; j++) { + image->r(j, col) = (red[j] + image->r(j, col + 1)) / 2; + image->g(j, col) = (green[j] + image->g(j, col + 1)) / 2; + image->b(j, col) = (blue[j] + image->b(j, col + 1)) / 2; + } + + col = 2 * imheight - 1 - 2 * i + 2; + + for (int j = 0; j < imwidth; j++) { + image->r(j, col) = (red[j] + image->r(j, col + 1)) / 2; + image->g(j, col) = (green[j] + image->g(j, col + 1)) / 2; + image->b(j, col) = (blue[j] + image->b(j, col + 1)) / 2; + } + } else if (i > 2 && i < imheight - 1) { // vertical bicubic interpolation + int col = 2 * imheight - 1 - 2 * i + 2; + + for (int j = 0; j < imwidth; j++) { + image->r(j, col) = MAX(0.f, -0.0625f * (red[j] + image->r(j, col + 3)) + 0.5625f * (image->r(j, col - 1) + image->r(j, col + 1))); + image->g(j, col) = MAX(0.f, -0.0625f * (green[j] + image->g(j, col + 3)) + 0.5625f * (image->g(j, col - 1) + image->g(j, col + 1))); + image->b(j, col) = MAX(0.f, -0.0625f * (blue[j] + image->b(j, col + 3)) + 0.5625f * (image->b(j, col - 1) + image->b(j, col + 1))); + + if(clip) { + image->r(j, col) = MIN(image->r(j, col), rtengine::MAXVALF); + image->g(j, col) = MIN(image->g(j, col), rtengine::MAXVALF); + image->b(j, col) = MIN(image->b(j, col), rtengine::MAXVALF); + } + } + } + + break; + + case TR_R270: // rotate left + if (i == 0) { + for (int j = imwidth - 1, row = 0; j >= 0; j--, row++) { + image->r(row, 2 * i) = red[j]; + image->g(row, 2 * i) = green[j]; + image->b(row, 2 * i) = blue[j]; + } + } else if (i == 1 || i == 2) { // linear interpolation + for (int j = imwidth - 1, row = 0; j >= 0; j--, row++) { + image->r(row, 2 * i) = red[j]; + image->g(row, 2 * i) = green[j]; + image->b(row, 2 * i) = blue[j]; + image->r(row, 2 * i - 1) = (red[j] + image->r(row, 2 * i - 2)) * 0.5f; + image->g(row, 2 * i - 1) = (green[j] + image->g(row, 2 * i - 2)) * 0.5f; + image->b(row, 2 * i - 1) = (blue[j] + image->b(row, 2 * i - 2)) * 0.5f; + } + } else if (i > 0 && i < imheight) { // vertical bicubic interpolation + for (int j = imwidth - 1, row = 0; j >= 0; j--, row++) { + image->r(row, 2 * i - 3) = MAX(0.f, -0.0625f * (red[j] + image->r(row, 2 * i - 6)) + 0.5625f * (image->r(row, 2 * i - 2) + image->r(row, 2 * i - 4))); + image->g(row, 2 * i - 3) = MAX(0.f, -0.0625f * (green[j] + image->g(row, 2 * i - 6)) + 0.5625f * (image->g(row, 2 * i - 2) + image->g(row, 2 * i - 4))); + image->b(row, 2 * i - 3) = MAX(0.f, -0.0625f * (blue[j] + image->b(row, 2 * i - 6)) + 0.5625f * (image->b(row, 2 * i - 2) + image->b(row, 2 * i - 4))); + + if(clip) { + image->r(row, 2 * i - 3) = MIN(image->r(row, 2 * i - 3), rtengine::MAXVALF); + image->g(row, 2 * i - 3) = MIN(image->g(row, 2 * i - 3), rtengine::MAXVALF); + image->b(row, 2 * i - 3) = MIN(image->b(row, 2 * i - 3), rtengine::MAXVALF); + } + + image->r(row, 2 * i) = red[j]; + image->g(row, 2 * i) = green[j]; + image->b(row, 2 * i) = blue[j]; + } + } + + if (i == imheight - 1) { + for (int j = imwidth - 1, row = 0; j >= 0; j--, row++) { + image->r(row, 2 * i - 1) = MAX(0.f, -0.0625f * (red[j] + image->r(row, 2 * i - 4)) + 0.5625f * (image->r(row, 2 * i) + image->r(row, 2 * i - 2))); + image->g(row, 2 * i - 1) = MAX(0.f, -0.0625f * (green[j] + image->g(row, 2 * i - 4)) + 0.5625f * (image->g(row, 2 * i) + image->g(row, 2 * i - 2))); + image->b(row, 2 * i - 1) = MAX(0.f, -0.0625f * (blue[j] + image->b(row, 2 * i - 4)) + 0.5625f * (image->b(row, 2 * i) + image->b(row, 2 * i - 2))); + + if(clip) { + image->r(j, 2 * i - 1) = MIN(image->r(j, 2 * i - 1), rtengine::MAXVALF); + image->g(j, 2 * i - 1) = MIN(image->g(j, 2 * i - 1), rtengine::MAXVALF); + image->b(j, 2 * i - 1) = MIN(image->b(j, 2 * i - 1), rtengine::MAXVALF); + } + + image->r(row, 2 * i + 1) = (red[j] + image->r(row, 2 * i - 1)) / 2; + image->g(row, 2 * i + 1) = (green[j] + image->g(row, 2 * i - 1)) / 2; + image->b(row, 2 * i + 1) = (blue[j] + image->b(row, 2 * i - 1)) / 2; + + if (oddHeight) { + image->r(row, 2 * i + 2) = (red[j] + image->r(row, 2 * i - 2)) / 2; + image->g(row, 2 * i + 2) = (green[j] + image->g(row, 2 * i - 2)) / 2; + image->b(row, 2 * i + 2) = (blue[j] + image->b(row, 2 * i - 2)) / 2; + } + } + } + + break; + + case TR_NONE: // no coarse rotation + default: + rotateLine (red, image->r, tran, 2 * i, imwidth, imheight); + rotateLine (green, image->g, tran, 2 * i, imwidth, imheight); + rotateLine (blue, image->b, tran, 2 * i, imwidth, imheight); + + if (i == 1 || i == 2) { // linear interpolation + for (int j = 0; j < imwidth; j++) { + image->r(2 * i - 1, j) = (red[j] + image->r(2 * i - 2, j)) / 2; + image->g(2 * i - 1, j) = (green[j] + image->g(2 * i - 2, j)) / 2; + image->b(2 * i - 1, j) = (blue[j] + image->b(2 * i - 2, j)) / 2; + } + } else if (i > 2 && i < imheight) { // vertical bicubic interpolation + for (int j = 0; j < imwidth; j++) { + image->r(2 * i - 3, j) = MAX(0.f, -0.0625f * (red[j] + image->r(2 * i - 6, j)) + 0.5625f * (image->r(2 * i - 2, j) + image->r(2 * i - 4, j))); + image->g(2 * i - 3, j) = MAX(0.f, -0.0625f * (green[j] + image->g(2 * i - 6, j)) + 0.5625f * (image->g(2 * i - 2, j) + image->g(2 * i - 4, j))); + image->b(2 * i - 3, j) = MAX(0.f, -0.0625f * (blue[j] + image->b(2 * i - 6, j)) + 0.5625f * (image->b(2 * i - 2, j) + image->b(2 * i - 4, j))); + + if(clip) { + image->r(2 * i - 3, j) = MIN(image->r(2 * i - 3, j), rtengine::MAXVALF); + image->g(2 * i - 3, j) = MIN(image->g(2 * i - 3, j), rtengine::MAXVALF); + image->b(2 * i - 3, j) = MIN(image->b(2 * i - 3, j), rtengine::MAXVALF); + } + } + } + + if (i == imheight - 1) { + for (int j = 0; j < imwidth; j++) { + image->r(2 * i - 1, j) = MAX(0.f, -0.0625f * (red[j] + image->r(2 * i - 4, j)) + 0.5625f * (image->r(2 * i, j) + image->r(2 * i - 2, j))); + image->g(2 * i - 1, j) = MAX(0.f, -0.0625f * (green[j] + image->g(2 * i - 4, j)) + 0.5625f * (image->g(2 * i, j) + image->g(2 * i - 2, j))); + image->b(2 * i - 1, j) = MAX(0.f, -0.0625f * (blue[j] + image->b(2 * i - 4, j)) + 0.5625f * (image->b(2 * i, j) + image->b(2 * i - 2, j))); + + if(clip) { + image->r(2 * i - 1, j) = MIN(image->r(2 * i - 1, j), rtengine::MAXVALF); + image->g(2 * i - 1, j) = MIN(image->g(2 * i - 1, j), rtengine::MAXVALF); + image->b(2 * i - 1, j) = MIN(image->b(2 * i - 1, j), rtengine::MAXVALF); + } + + image->r(2 * i + 1, j) = (red[j] + image->r(2 * i - 1, j)) / 2; + image->g(2 * i + 1, j) = (green[j] + image->g(2 * i - 1, j)) / 2; + image->b(2 * i + 1, j) = (blue[j] + image->b(2 * i - 1, j)) / 2; + + if (oddHeight) { + image->r(2 * i + 2, j) = (red[j] + image->r(2 * i - 2, j)) / 2; + image->g(2 * i + 2, j) = (green[j] + image->g(2 * i - 2, j)) / 2; + image->b(2 * i + 2, j) = (blue[j] + image->b(2 * i - 2, j)) / 2; + } + } + } + } +} + +} + + namespace rtengine { @@ -116,15 +481,11 @@ RawImageSource::~RawImageSource () if (hrmap[0] != NULL) { int dh = H / HR_SCALE; - freeArray(hrmap[0], dh); - freeArray(hrmap[1], dh); - freeArray(hrmap[2], dh); + freeJaggedArray(hrmap[0]); + freeJaggedArray(hrmap[1]); + freeJaggedArray(hrmap[2]); } - //if (needhr) - // freeArray(needhr, H); - //if (hpmap) - // freeArray(hpmap, H); if (camProfile) { cmsCloseProfile (camProfile); } @@ -308,7 +669,7 @@ void RawImageSource::getImage (ColorTemp ctemp, int tran, Imagefloat* image, Pre defGain = 0.0; // compute image area to render in order to provide the requested part of the image - int sx1, sy1, imwidth, imheight, fw; + int sx1, sy1, imwidth, imheight, fw, d1xHeightOdd; transformRect (pp, tran, sx1, sy1, imwidth, imheight, fw); // check possible overflows @@ -323,7 +684,12 @@ void RawImageSource::getImage (ColorTemp ctemp, int tran, Imagefloat* image, Pre } if (d1x) { + // D1X has only half of the required rows + // we interpolate the missing ones later to get correct aspect ratio + // if the height is odd we also have to add an additional row to avoid a black line + d1xHeightOdd = maximheight & 1; maximheight /= 2; + imheight = maximheight; } // correct if overflow (very rare), but not fuji because it is corrected in transline @@ -341,9 +707,9 @@ void RawImageSource::getImage (ColorTemp ctemp, int tran, Imagefloat* image, Pre hlmax[0] = clmax[0] * rm; hlmax[1] = clmax[1] * gm; hlmax[2] = clmax[2] * bm; - const bool has_clipping = (chmax[0] >= clmax[0] || chmax[1] >= clmax[1] || chmax[2] >= clmax[2]); - //if (sx1+skip*imwidth>maxx) imwidth --; // very hard to fix this situation without an 'if' in the loop. + const bool doClip = (chmax[0] >= clmax[0] || chmax[1] >= clmax[1] || chmax[2] >= clmax[2]) && !hrp.hrenabled; + float area = skip * skip; rm /= area; gm /= area; @@ -354,11 +720,9 @@ void RawImageSource::getImage (ColorTemp ctemp, int tran, Imagefloat* image, Pre { #endif // render the requested image part - float* line_red = new float[imwidth]; - float* line_grn = new float[imwidth]; - float* line_blue = new float[imwidth]; - //printf("clip[0]=%f clip[1]=%f clip[2]=%f\n",hlmax[0],hlmax[1],hlmax[2]); - + float line_red[imwidth] ALIGNED16; + float line_grn[imwidth] ALIGNED16; + float line_blue[imwidth] ALIGNED16; #ifdef _OPENMP #pragma omp for @@ -373,12 +737,9 @@ void RawImageSource::getImage (ColorTemp ctemp, int tran, Imagefloat* image, Pre if (ri->getSensorType() == ST_BAYER || ri->getSensorType() == ST_FUJI_XTRANS || ri->get_colors() == 1) { for (int j = 0, jx = sx1; j < imwidth; j++, jx += skip) { - if (jx >= maxx - skip) { - jx = maxx - skip - 1; // avoid trouble - } + jx = jx >= (maxx - skip) ? jx = maxx - skip - 1 : jx; // avoid trouble - float rtot, gtot, btot; - rtot = gtot = btot = 0; + float rtot = 0.f, gtot = 0.f, btot = 0.f; for (int m = 0; m < skip; m++) for (int n = 0; n < skip; n++) { @@ -391,12 +752,12 @@ void RawImageSource::getImage (ColorTemp ctemp, int tran, Imagefloat* image, Pre gtot *= gm; btot *= bm; - if (!hrp.hrenabled && has_clipping) { + if (doClip) { // note: as hlmax[] can be larger than CLIP and we can later apply negative // exposure this means that we can clip away local highlights which actually // are not clipped. We have to do that though as we only check pixel by pixel // and don't know if this will transition into a clipped area, if so we need - // to clip also surrounding to make a good color transition + // to clip also surrounding to make a good colour transition rtot = CLIP(rtot); gtot = CLIP(gtot); btot = CLIP(btot); @@ -426,7 +787,7 @@ void RawImageSource::getImage (ColorTemp ctemp, int tran, Imagefloat* image, Pre gtot *= gm; btot *= bm; - if (!hrp.hrenabled && has_clipping) { + if (doClip) { rtot = CLIP(rtot); gtot = CLIP(gtot); btot = CLIP(btot); @@ -444,13 +805,16 @@ void RawImageSource::getImage (ColorTemp ctemp, int tran, Imagefloat* image, Pre hlRecovery (hrp.method, line_red, line_grn, line_blue, i, sx1, imwidth, skip, raw, hlmax); } - transLine (line_red, line_grn, line_blue, ix, image, tran, imwidth, imheight, fw); + if(d1x) { + transLineD1x (line_red, line_grn, line_blue, ix, image, tran, imwidth, imheight, d1xHeightOdd, doClip); + } else if(fuji) { + transLineFuji (line_red, line_grn, line_blue, ix, image, tran, imwidth, imheight, fw); + } else { + transLineStandard (line_red, line_grn, line_blue, ix, image, tran, imwidth, imheight); + } } - delete [] line_red; - delete [] line_grn; - delete [] line_blue; #ifdef _OPENMP } #endif @@ -503,8 +867,6 @@ void RawImageSource::getImage (ColorTemp ctemp, int tran, Imagefloat* image, Pre } } - - // Flip if needed if (tran & TR_HFLIP) { hflip (image); @@ -514,11 +876,14 @@ void RawImageSource::getImage (ColorTemp ctemp, int tran, Imagefloat* image, Pre vflip (image); } - // Color correction (only when running on full resolution) - if (ri->getSensorType() != ST_NONE && pp.skip == 1) { - if (ri->getSensorType() == ST_BAYER) { + // Colour correction (only when running on full resolution) + if(pp.skip == 1) { + switch(ri->getSensorType()) { + case ST_BAYER: processFalseColorCorrection (image, raw.bayersensor.ccSteps); - } else if (ri->getSensorType() == ST_FUJI_XTRANS) { + break; + + case ST_FUJI_XTRANS: processFalseColorCorrection (image, raw.xtranssensor.ccSteps); } } @@ -547,7 +912,7 @@ void RawImageSource::convertColorSpace(Imagefloat* image, ColorManagementParams //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% /* interpolateBadPixelsBayer: 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::interpolateBadPixelsBayer( PixelsMap &bitmapBads ) { @@ -574,8 +939,8 @@ int RawImageSource::interpolateBadPixelsBayer( PixelsMap &bitmapBads ) // diagonal interpolation if(FC(row, col) == 1) { - // green channel. We can use closer pixels than for red or blue channel. Distance to center pixel is sqrt(2) => weighting is 0.70710678 - // For green channel following pixels will be used for interpolation. Pixel to be interpolated is in center. + // green channel. We can use closer pixels than for red or blue channel. Distance to centre pixel is sqrt(2) => weighting is 0.70710678 + // For green channel following pixels will be used for interpolation. Pixel to be interpolated is in centre. // 1 means that pixel is used in this step, if itself and his counterpart are not marked bad // 0 0 0 0 0 // 0 1 0 1 0 @@ -592,8 +957,8 @@ int RawImageSource::interpolateBadPixelsBayer( PixelsMap &bitmapBads ) norm += dirwt; } } else { - // red and blue channel. Distance to center pixel is sqrt(8) => weighting is 0.35355339 - // For red and blue channel following pixels will be used for interpolation. Pixel to be interpolated is in center. + // red and blue channel. Distance to centre pixel is sqrt(8) => weighting is 0.35355339 + // For red and blue channel following pixels will be used for interpolation. Pixel to be interpolated is in centre. // 1 means that pixel is used in this step, if itself and his counterpart are not marked bad // 1 0 0 0 1 // 0 0 0 0 0 @@ -611,8 +976,8 @@ int RawImageSource::interpolateBadPixelsBayer( PixelsMap &bitmapBads ) } } - // channel independent. Distance to center pixel is 2 => weighting is 0.5 - // Additionally for all channel following pixels will be used for interpolation. Pixel to be interpolated is in center. + // channel independent. Distance to centre pixel is 2 => weighting is 0.5 + // Additionally for all channel following pixels will be used for interpolation. Pixel to be interpolated is in centre. // 1 means that pixel is used in this step, if itself and his counterpart are not marked bad // 0 0 1 0 0 // 0 0 0 0 0 @@ -634,7 +999,7 @@ int RawImageSource::interpolateBadPixelsBayer( 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 counter++; } else { //backup plan -- simple average. Same method for all channels. We could improve this, but it's really unlikely that this case happens @@ -664,7 +1029,7 @@ int RawImageSource::interpolateBadPixelsBayer( PixelsMap &bitmapBads ) } /* interpolateBadPixels3Colours: 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::interpolateBadPixelsNColours( PixelsMap &bitmapBads, const int colours ) { @@ -1043,6 +1408,7 @@ SSEFUNCTION int RawImageSource::findHotDeadPixels( PixelsMap &bpMap, float thres #ifdef __SSE2__ // sum up 5*4 = 20 values using SSE + // 10 fabs function calls and float 10 additions with 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])); @@ -1056,7 +1422,7 @@ SSEFUNCTION int RawImageSource::findHotDeadPixels( PixelsMap &bpMap, float thres } #else - + // 25 fabs function calls and 25 float additions without SSE for (int mm = rr - 2; mm <= rr + 2; mm++) { for (int nn = cc - 2; nn <= cc + 2; nn++) { hfnbrave += fabsf(cfablur[mm * W + nn]); @@ -1077,250 +1443,6 @@ SSEFUNCTION int RawImageSource::findHotDeadPixels( PixelsMap &bpMap, float thres return counter; } -void RawImageSource::rotateLine (float* line, PlanarPtr &channel, int tran, int i, int w, int h) -{ - - if ((tran & TR_ROT) == TR_R180) - for (int j = 0; j < w; j++) { - channel(h - 1 - i, w - 1 - j) = line[j]; - } - - else if ((tran & TR_ROT) == TR_R90) - for (int j = 0; j < w; j++) { - channel(j, h - 1 - i) = line[j]; - } - - else if ((tran & TR_ROT) == TR_R270) - for (int j = 0; j < w; j++) { - channel(w - 1 - j, i) = line[j]; - } - else - for (int j = 0; j < w; j++) { - channel(i, j) = line[j]; - } -} - -//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -void RawImageSource::transLine (float* red, float* green, float* blue, int i, Imagefloat* image, int tran, int imwidth, int imheight, int fw) -{ - - // Fuji SuperCCD rotation + coarse rotation - if (fuji) { - int start = ABS(fw - i); - int w = fw * 2 + 1; - int h = (imheight - fw) * 2 + 1; - - if ((tran & TR_ROT) == TR_R180) { - int end = min(h + fw - i, w - fw + i); - - for (int j = start; j < end; j++) { - int y = i + j - fw; - int x = fw - i + j; - - if (x >= 0 && y < image->height && y >= 0 && x < image->width) { - image->r(image->height - 1 - y, image->width - 1 - x) = red[j]; - image->g(image->height - 1 - y, image->width - 1 - x) = green[j]; - image->b(image->height - 1 - y, image->width - 1 - x) = blue[j]; - } - } - } else if ((tran & TR_ROT) == TR_R270) { - int end = min(h + fw - i, w - fw + i); - - for (int j = start; j < end; j++) { - int y = i + j - fw; - int x = fw - i + j; - - if (x >= 0 && x < image->height && y >= 0 && y < image->width) { - image->r(image->height - 1 - x, y) = red[j]; - image->g(image->height - 1 - x, y) = green[j]; - image->b(image->height - 1 - x, y) = blue[j]; - } - } - } else if ((tran & TR_ROT) == TR_R90) { - int end = min(h + fw - i, w - fw + i); - - for (int j = start; j < end; j++) { - int y = i + j - fw; - int x = fw - i + j; - - if (x >= 0 && y < image->width && y >= 0 && x < image->height) { - image->r(x, image->width - 1 - y) = red[j]; - image->g(x, image->width - 1 - y) = green[j]; - image->b(x, image->width - 1 - y) = blue[j]; - } - } - } else { - int end = min(h + fw - i, w - fw + i); - - for (int j = start; j < end; j++) { - int y = i + j - fw; - int x = fw - i + j; - - if (x >= 0 && y < image->height && y >= 0 && x < image->width) { - image->r(y, x) = red[j]; - image->g(y, x) = green[j]; - image->b(y, x) = blue[j]; - } - } - } - } - // Nikon D1X vertical interpolation + coarse rotation - else if (d1x) { - // copy new pixels - if ((tran & TR_ROT) == TR_R180) { - for (int j = 0; j < imwidth; j++) { - image->r(2 * imheight - 2 - 2 * i, imwidth - 1 - j) = red[j]; - image->g(2 * imheight - 2 - 2 * i, imwidth - 1 - j) = green[j]; - image->b(2 * imheight - 2 - 2 * i, imwidth - 1 - j) = blue[j]; - } - - if (i == 1 || i == 2) { // linear interpolation - int row = 2 * imheight - 1 - 2 * i; - - for (int j = 0; j < imwidth; j++) { - int col = imwidth - 1 - j; - image->r(row, col) = (red[j] + image->r(row + 1, col)) / 2; - image->g(row, col) = (green[j] + image->g(row + 1, col)) / 2; - image->b(row, col) = (blue[j] + image->b(row + 1, col)) / 2; - } - } else if (i == imheight - 1) { - int row = 2 * imheight - 1 - 2 * i; - - for (int j = 0; j < imwidth; j++) { - int col = imwidth - 1 - j; - image->r(row, col) = (red[j] + image->r(row + 1, col)) / 2; - image->g(row, col) = (green[j] + image->g(row + 1, col)) / 2; - image->b(row, col) = (blue[j] + image->b(row + 1, col)) / 2; - } - - row = 2 * imheight - 1 - 2 * i + 2; - - for (int j = 0; j < imwidth; j++) { - int col = imwidth - 1 - j; - image->r(row, col) = (red[j] + image->r(row + 1, col)) / 2; - image->g(row, col) = (green[j] + image->g(row + 1, col)) / 2; - image->b(row, col) = (blue[j] + image->b(row + 1, col)) / 2; - } - } else if (i > 2 && i < imheight - 1) { // vertical bicubic interpolation - int row = 2 * imheight - 1 - 2 * i + 2; - - for (int j = 0; j < imwidth; j++) { - int col = imwidth - 1 - j; - image->r(row, col) = CLIP((int)(-0.0625 * red[j] + 0.5625 * image->r(row - 1, col) + 0.5625 * image->r(row + 1, col) - 0.0625 * image->r(row + 3, col))); - image->g(row, col) = CLIP((int)(-0.0625 * green[j] + 0.5625 * image->g(row - 1, col) + 0.5625 * image->g(row + 1, col) - 0.0625 * image->g(row + 3, col))); - image->b(row, col) = CLIP((int)(-0.0625 * blue[j] + 0.5625 * image->b(row - 1, col) + 0.5625 * image->b(row + 1, col) - 0.0625 * image->b(row + 3, col))); - } - } - } else if ((tran & TR_ROT) == TR_R90) { - for (int j = 0; j < imwidth; j++) { - image->r(j, 2 * imheight - 2 - 2 * i) = red[j]; - image->g(j, 2 * imheight - 2 - 2 * i) = green[j]; - image->b(j, 2 * imheight - 2 - 2 * i) = blue[j]; - } - - if (i == 1 || i == 2) { // linear interpolation - int col = 2 * imheight - 1 - 2 * i; - - for (int j = 0; j < imwidth; j++) { - image->r(j, col) = (red[j] + image->r(j, col + 1)) / 2; - image->g(j, col) = (green[j] + image->g(j, col + 1)) / 2; - image->b(j, col) = (blue[j] + image->b(j, col + 1)) / 2; - } - } else if (i == imheight - 1) { - int col = 2 * imheight - 1 - 2 * i; - - for (int j = 0; j < imwidth; j++) { - image->r(j, col) = (red[j] + image->r(j, col + 1)) / 2; - image->g(j, col) = (green[j] + image->g(j, col + 1)) / 2; - image->b(j, col) = (blue[j] + image->b(j, col + 1)) / 2; - } - - col = 2 * imheight - 1 - 2 * i + 2; - - for (int j = 0; j < imwidth; j++) { - image->r(j, col) = (red[j] + image->r(j, col + 1)) / 2; - image->g(j, col) = (green[j] + image->g(j, col + 1)) / 2; - image->b(j, col) = (blue[j] + image->b(j, col + 1)) / 2; - } - } else if (i > 2 && i < imheight - 1) { // vertical bicubic interpolation - int col = 2 * imheight - 1 - 2 * i + 2; - - for (int j = 0; j < imwidth; j++) { - image->r(j, col) = CLIP((int)(-0.0625 * red[j] + 0.5625 * image->r(j, col - 1) + 0.5625 * image->r(j, col + 1) - 0.0625 * image->r(j, col + 3))); - image->g(j, col) = CLIP((int)(-0.0625 * green[j] + 0.5625 * image->g(j, col - 1) + 0.5625 * image->g(j, col + 1) - 0.0625 * image->g(j, col + 3))); - image->b(j, col) = CLIP((int)(-0.0625 * blue[j] + 0.5625 * image->b(j, col - 1) + 0.5625 * image->b(j, col + 1) - 0.0625 * image->b(j, col + 3))); - } - } - } else if ((tran & TR_ROT) == TR_R270) { - for (int j = 0; j < imwidth; j++) { - image->r(imwidth - 1 - j, 2 * i) = red[j]; - image->g(imwidth - 1 - j, 2 * i) = green[j]; - image->b(imwidth - 1 - j, 2 * i) = blue[j]; - } - - if (i == 1 || i == 2) { // linear interpolation - for (int j = 0; j < imwidth; j++) { - int row = imwidth - 1 - j; - image->r(row, 2 * i - 1) = (red[j] + image->r(row, 2 * i - 2)) * 0.5f; - image->g(row, 2 * i - 1) = (green[j] + image->g(row, 2 * i - 2)) * 0.5f; - image->b(row, 2 * i - 1) = (blue[j] + image->b(row, 2 * i - 2)) * 0.5f; - } - } else if (i == imheight - 1) { - for (int j = 0; j < imwidth; j++) { - int row = imwidth - 1 - j; - image->r(row, 2 * i - 1) = (red[j] + image->r(row, 2 * i - 2)) * 0.5f; - image->g(row, 2 * i - 1) = (green[j] + image->g(row, 2 * i - 2)) * 0.5f; - image->b(row, 2 * i - 1) = (blue[j] + image->b(row, 2 * i - 2)) * 0.5f; - image->r(row, 2 * i - 3) = (image->r(row, 2 * i - 2) + image->r(row, 2 * i - 4)) * 0.5f; - image->g(row, 2 * i - 3) = (image->g(row, 2 * i - 2) + image->g(row, 2 * i - 4)) * 0.5f; - image->b(row, 2 * i - 3) = (image->b(row, 2 * i - 2) + image->b(row, 2 * i - 4)) * 0.5f; - } - } else if (i > 0 && i < imheight - 1) { // vertical bicubic interpolationi - for (int j = 0; j < imwidth; j++) { - int row = imwidth - 1 - j; - image->r(row, 2 * i - 3) = CLIP((int)(-0.0625 * red[j] + 0.5625 * image->r(row, 2 * i - 2) + 0.5625 * image->r(row, 2 * i - 4) - 0.0625 * image->r(row, 2 * i - 6))); - image->g(row, 2 * i - 3) = CLIP((int)(-0.0625 * green[j] + 0.5625 * image->g(row, 2 * i - 2) + 0.5625 * image->g(row, 2 * i - 4) - 0.0625 * image->g(row, 2 * i - 6))); - image->b(row, 2 * i - 3) = CLIP((int)(-0.0625 * blue[j] + 0.5625 * image->b(row, 2 * i - 2) + 0.5625 * image->b(row, 2 * i - 4) - 0.0625 * image->b(row, 2 * i - 6))); - } - } - } else { - rotateLine (red, image->r, tran, 2 * i, imwidth, imheight); - rotateLine (green, image->g, tran, 2 * i, imwidth, imheight); - rotateLine (blue, image->b, tran, 2 * i, imwidth, imheight); - - if (i == 1 || i == 2) { // linear interpolation - for (int j = 0; j < imwidth; j++) { - image->r(2 * i - 1, j) = (red[j] + image->r(2 * i - 2, j)) / 2; - image->g(2 * i - 1, j) = (green[j] + image->g(2 * i - 2, j)) / 2; - image->b(2 * i - 1, j) = (blue[j] + image->b(2 * i - 2, j)) / 2; - } - } else if (i == imheight - 1) { - for (int j = 0; j < imwidth; j++) { - image->r(2 * i - 3, j) = (image->r(2 * i - 4, j) + image->r(2 * i - 2, j)) / 2; - image->g(2 * i - 3, j) = (image->g(2 * i - 4, j) + image->g(2 * i - 2, j)) / 2; - image->b(2 * i - 3, j) = (image->b(2 * i - 4, j) + image->b(2 * i - 2, j)) / 2; - image->r(2 * i - 1, j) = (red[j] + image->r(2 * i - 2, j)) / 2; - image->g(2 * i - 1, j) = (green[j] + image->g(2 * i - 2, j)) / 2; - image->b(2 * i - 1, j) = (blue[j] + image->b(2 * i - 2, j)) / 2; - } - } else if (i > 2 && i < imheight - 1) { // vertical bicubic interpolationi - for (int j = 0; j < imwidth; j++) { - image->r(2 * i - 3, j) = CLIP((int)(-0.0625 * red[j] + 0.5625 * image->r(2 * i - 2, j) + 0.5625 * image->r(2 * i - 4, j) - 0.0625 * image->r(2 * i - 6, j))); - image->g(2 * i - 3, j) = CLIP((int)(-0.0625 * green[j] + 0.5625 * image->g(2 * i - 2, j) + 0.5625 * image->g(2 * i - 4, j) - 0.0625 * image->g(2 * i - 6, j))); - image->b(2 * i - 3, j) = CLIP((int)(-0.0625 * blue[j] + 0.5625 * image->b(2 * i - 2, j) + 0.5625 * image->b(2 * i - 4, j) - 0.0625 * image->b(2 * i - 6, j))); - } - } - } - } // if nikon dx1 - // other (conventional) CCD coarse rotation - else { - rotateLine (red, image->r, tran, i, imwidth, imheight); - rotateLine (green, image->g, tran, i, imwidth, imheight); - rotateLine (blue, image->b, tran, i, imwidth, imheight); - } -} - //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% void RawImageSource::getFullSize (int& w, int& h, int tr) @@ -1333,7 +1455,7 @@ void RawImageSource::getFullSize (int& w, int& h, int tr) h = (H - ri->get_FujiWidth()) * 2 + 1; } else if (d1x) { w = W; - h = 2 * H - 1; + h = 2 * H; } else { w = W; h = H; @@ -1355,17 +1477,8 @@ void RawImageSource::getSize (int tran, PreviewProps pp, int& w, int& h) { tran = defTransform (tran); - -// if (fuji) { -// return; -// } -// else if (d1x) { -// return; -// } -// else { w = pp.w / pp.skip + (pp.w % pp.skip > 0); h = pp.h / pp.skip + (pp.h % pp.skip > 0); -// } } //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -1426,10 +1539,6 @@ int RawImageSource::load (Glib::ustring fname, bool batch) d1x = ! ri->get_model().compare("D1X"); - if (d1x) { - border = 8; - } - if(ri->getSensorType() == ST_FUJI_XTRANS) { border = 7; } else if(ri->getSensorType() == ST_FOVEON) { @@ -1991,11 +2100,11 @@ void RawImageSource::retinexPrepareBuffers(ColorManagementParams cmp, RetinexPar vfloat H, S, L; int pos; Color::rgb2hsl(LVFU(red[i][j]), LVFU(green[i][j]), LVFU(blue[i][j]), H, S, L); - _mm_storeu_ps(&conversionBuffer[0][i - border][j - border], H); - _mm_storeu_ps(&conversionBuffer[1][i - border][j - border], S); + STVFU(conversionBuffer[0][i - border][j - border], H); + STVFU(conversionBuffer[1][i - border][j - border], S); L *= c32768; - _mm_storeu_ps(&conversionBuffer[2][i - border][j - border], L); - _mm_storeu_ps(&conversionBuffer[3][i - border][j - border], H); + STVFU(conversionBuffer[2][i - border][j - border], L); + STVFU(conversionBuffer[3][i - border][j - border], H); if(lhist16RETI) { for(int p = 0; p < 4; p++) { @@ -2116,6 +2225,7 @@ void RawImageSource::retinexPrepareCurves(RetinexParams retinexParams, LUTf &cdc } else { CurveFactory::curveDehaContL (retinexcontlutili, retinexParams.cdcurve, cdcurve, 1, lhist16RETI, histLRETI); } + CurveFactory::mapcurve (mapcontlutili, retinexParams.mapcurve, mapcurve, 1, lhist16RETI, histLRETI); retinexParams.getCurves(retinextransmissionCurve); @@ -2314,9 +2424,9 @@ void RawImageSource::retinex(ColorManagementParams cmp, RetinexParams deh, ToneC vfloat R, G, B; Color::hsl2rgb(LVFU(conversionBuffer[0][i - border][j - border]), LVFU(conversionBuffer[1][i - border][j - border]), LVFU(LBuffer[i - border][j - border]) / c32768, R, G, B); - _mm_storeu_ps(&red[i][j], R); - _mm_storeu_ps(&green[i][j], G); - _mm_storeu_ps(&blue[i][j], B); + STVFU(red[i][j], R); + STVFU(green[i][j], G); + STVFU(blue[i][j], B); } #endif @@ -2401,9 +2511,9 @@ void RawImageSource::retinex(ColorManagementParams cmp, RetinexParams deh, ToneC Color::Lab2XYZ(LVFU(LBuffer[i - border][j - border]), LVFU(conversionBuffer[0][i - border][j - border]), LVFU(conversionBuffer[1][i - border][j - border]), x_, y_, z_) ; Color::xyz2rgb(x_, y_, z_, R, G, B, wipv); - _mm_storeu_ps(&red[i][j], R); - _mm_storeu_ps(&green[i][j], G); - _mm_storeu_ps(&blue[i][j], B); + STVFU(red[i][j], R); + STVFU(green[i][j], G); + STVFU(blue[i][j], B); } @@ -2442,6 +2552,7 @@ void RawImageSource::retinex(ColorManagementParams cmp, RetinexParams deh, ToneC } } } + rgbSourceModified = false; // tricky handling for Color propagation t5.set(); @@ -2748,8 +2859,9 @@ void RawImageSource::copyOriginalPixels(const RAWParams &raw, RawImage *src, Raw { // TODO: Change type of black[] to float to avoid conversions unsigned short black[4] = { - (unsigned short)ri->get_cblack(0), (unsigned short)ri->get_cblack(1), - (unsigned short)ri->get_cblack(2), (unsigned short)ri->get_cblack(3)}; + (unsigned short)ri->get_cblack(0), (unsigned short)ri->get_cblack(1), + (unsigned short)ri->get_cblack(2), (unsigned short)ri->get_cblack(3) + }; if (ri->getSensorType() == ST_BAYER || ri->getSensorType() == ST_FUJI_XTRANS) { if (!rawData) { @@ -2923,9 +3035,9 @@ SSEFUNCTION void RawImageSource::cfaboxblur(RawImage *riFlatFile, float* cfablur //vertical blur #ifdef __SSE2__ - __m128 leninitv = _mm_set1_ps( (float)((int)(boxH / 2 + 1))); - __m128 onev = _mm_set1_ps( 1.0f ); - __m128 temp1v, temp2v, lenv, lenp1v, lenm1v; + vfloat leninitv = F2V(boxH / 2 + 1); + vfloat onev = F2V( 1.0f ); + vfloat temp1v, temp2v, lenv, lenp1v, lenm1v; int row; #ifdef _OPENMP #pragma omp for @@ -2941,29 +3053,29 @@ SSEFUNCTION void RawImageSource::cfaboxblur(RawImage *riFlatFile, float* cfablur temp2v += LVFU(cfatmp[(i + 1) * W + col]) / lenv; } - _mm_storeu_ps(&cfablur[0 * W + col], temp1v); - _mm_storeu_ps(&cfablur[1 * W + col], temp2v); + STVFU(cfablur[0 * W + col], temp1v); + STVFU(cfablur[1 * W + col], temp2v); for (row = 2; row < boxH + 2; row += 2) { lenp1v = lenv + onev; temp1v = (temp1v * lenv + LVFU(cfatmp[(row + boxH) * W + col])) / lenp1v; temp2v = (temp2v * lenv + LVFU(cfatmp[(row + boxH + 1) * W + col])) / lenp1v; - _mm_storeu_ps( &cfablur[row * W + col], temp1v); - _mm_storeu_ps( &cfablur[(row + 1)*W + col], temp2v); + STVFU(cfablur[row * W + col], temp1v); + STVFU(cfablur[(row + 1)*W + col], temp2v); lenv = lenp1v; } for (; row < H - boxH - 1; row += 2) { temp1v = temp1v + (LVFU(cfatmp[(row + boxH) * W + col]) - LVFU(cfatmp[(row - boxH - 2) * W + col])) / lenv; temp2v = temp2v + (LVFU(cfatmp[(row + 1 + boxH) * W + col]) - LVFU(cfatmp[(row + 1 - boxH - 2) * W + col])) / lenv; - _mm_storeu_ps(&cfablur[row * W + col], temp1v); - _mm_storeu_ps(&cfablur[(row + 1)*W + col], temp2v); + STVFU(cfablur[row * W + col], temp1v); + STVFU(cfablur[(row + 1)*W + col], temp2v); } for(; row < H - boxH; row++) { temp1v = temp1v + (LVFU(cfatmp[(row + boxH) * W + col]) - LVFU(cfatmp[(row - boxH - 2) * W + col])) / lenv; - _mm_storeu_ps(&cfablur[row * W + col], temp1v); - __m128 swapv = temp1v; + STVFU(cfablur[row * W + col], temp1v); + vfloat swapv = temp1v; temp1v = temp2v; temp2v = swapv; @@ -2973,15 +3085,15 @@ SSEFUNCTION void RawImageSource::cfaboxblur(RawImage *riFlatFile, float* cfablur lenm1v = lenv - onev; temp1v = (temp1v * lenv - LVFU(cfatmp[(row - boxH - 2) * W + col])) / lenm1v; temp2v = (temp2v * lenv - LVFU(cfatmp[(row - boxH - 1) * W + col])) / lenm1v; - _mm_storeu_ps(&cfablur[row * W + col], temp1v); - _mm_storeu_ps(&cfablur[(row + 1)*W + col], temp2v); + STVFU(cfablur[row * W + col], temp1v); + STVFU(cfablur[(row + 1)*W + col], temp2v); lenv = lenm1v; } for(; row < H; row++) { lenm1v = lenv - onev; temp1v = (temp1v * lenv - LVFU(cfatmp[(row - boxH - 2) * W + col])) / lenm1v; - _mm_storeu_ps(&cfablur[(row)*W + col], temp1v); + STVFU(cfablur[(row)*W + col], temp1v); } } @@ -4225,7 +4337,8 @@ void RawImageSource::getRAWHistogram (LUTu & histRedRaw, LUTu & histGreenRaw, LU const float mult[4] = { 65535.0f / ri->get_white(0), 65535.0f / ri->get_white(1), 65535.0f / ri->get_white(2), - 65535.0f / ri->get_white(3) }; + 65535.0f / ri->get_white(3) + }; #ifdef _OPENMP int numThreads; diff --git a/rtengine/rawimagesource.h b/rtengine/rawimagesource.h index 10057d58a..bef4ad4e2 100644 --- a/rtengine/rawimagesource.h +++ b/rtengine/rawimagesource.h @@ -31,38 +31,6 @@ namespace rtengine { -// these two functions "simulate" and jagged array, but just use two allocs -template T** allocArray (int W, int H, bool initZero = false) -{ - - T** t = new T*[H]; - t[0] = new T[H * W]; - - if (initZero) { - memset(t[0], 0, sizeof(T)*W * H); - } - - for (int i = 1; i < H; i++) { - t[i] = t[i - 1] + W; - } - - return t; -} - -template void freeArray (T** a, int H) -{ - - delete [] a[0]; - delete [] a; -} - - -template void freeArray2 (T** a, int H) -{ - //for (int i=0; i &channel, int tran, int i, int w, int h); void transformRect (PreviewProps pp, int tran, int &sx1, int &sy1, int &width, int &height, int &fw); void transformPosition (int x, int y, int tran, int& tx, int& ty); @@ -151,7 +118,6 @@ public: int load (Glib::ustring fname, bool batch = false); void preprocess (const RAWParams &raw, const LensProfParams &lensProf, const CoarseTransformParams& coarse); void demosaic (const RAWParams &raw); -// void retinex (RAWParams raw, ColorManagementParams cmp, RetinexParams lcur, LUTf & cdcurve, bool dehacontlutili); void retinex (ColorManagementParams cmp, RetinexParams deh, ToneCurveParams Tc, LUTf & cdcurve, LUTf & mapcurve, const RetinextransmissionCurve & dehatransmissionCurve, multi_array2D &conversionBuffer, bool dehacontlutili, bool mapcontlutili, bool useHsl, float &minCD, float &maxCD, float &mini, float &maxi, float &Tmean, float &Tsigma, float &Tmin, float &Tmax, LUTu &histLRETI); void retinexPrepareCurves (RetinexParams retinexParams, LUTf &cdcurve, LUTf &mapcurve, RetinextransmissionCurve &retinextransmissionCurve, bool &retinexcontlutili, bool &mapcontlutili, bool &useHsl, LUTu & lhist16RETI, LUTu & histLRETI); void retinexPrepareBuffers (ColorManagementParams cmp, RetinexParams retinexParams, multi_array2D &conversionBuffer, LUTu &lhist16RETI); @@ -231,13 +197,7 @@ public: void boxblur2(float** src, float** dst, float** temp, int H, int W, int box ); void boxblur_resamp(float **src, float **dst, float** temp, int H, int W, int box, int samp ); void MSR(float** luminance, float **originalLuminance, float **exLuminance, LUTf & mapcurve, bool &mapcontlutili, int width, int height, RetinexParams deh, const RetinextransmissionCurve & dehatransmissionCurve, float &minCD, float &maxCD, float &mini, float &maxi, float &Tmean, float &Tsigma, float &Tmin, float &Tmax); -// void MSR(LabImage* lab, int width, int height, int skip, RetinexParams deh, const RetinextransmissionCurve & dehatransmissionCurve); - - //void boxblur_resamp(float **red, float **green, float **blue, int H, int W, float thresh[3], float max[3], - // multi_array2D & hfsize, multi_array2D & hilite, int box ); void HLRecovery_inpaint (float** red, float** green, float** blue); - //void HLRecovery_inpaint (); - static void HLRecovery_Luminance (float* rin, float* gin, float* bin, float* rout, float* gout, float* bout, int width, float maxval); static void HLRecovery_CIELab (float* rin, float* gin, float* bin, float* rout, float* gout, float* bout, int width, float maxval, double cam[3][3], double icam[3][3]); static void HLRecovery_blend (float* rin, float* gin, float* bin, int width, float maxval, float* hlmax); @@ -277,8 +237,6 @@ protected: void jdl_interpolate_omp(); void igv_interpolate(int winw, int winh); void lmmse_interpolate_omp(int winw, int winh, int iterations); -// void MSR(LabImage* lab, int width, int height, int skip, const LCurveParams &lcur); - void amaze_demosaic_RT(int winx, int winy, int winw, int winh);//Emil's code for AMaZE void fast_demosaic(int winx, int winy, int winw, int winh );//Emil's code for fast demosaicing void dcb_demosaic(int iterations, bool dcb_enhance); @@ -303,7 +261,6 @@ protected: void xtransborder_interpolate (int border); void xtrans_interpolate (int passes, bool useCieLab); void fast_xtrans_interpolate (); - void transLine (float* red, float* green, float* blue, int i, Imagefloat* image, int tran, int imw, int imh, int fw); void hflip (Imagefloat* im); void vflip (Imagefloat* im); diff --git a/rtengine/shmap.cc b/rtengine/shmap.cc index 4f93c80d3..1783a6bc8 100644 --- a/rtengine/shmap.cc +++ b/rtengine/shmap.cc @@ -21,6 +21,7 @@ #include "rtengine.h" #include "rt_math.h" #include "rawimagesource.h" +#include "jaggedarray.h" #undef THREAD_PRIORITY_NORMAL #include "opthelper.h" @@ -114,7 +115,7 @@ void SHMap::update (Imagefloat* img, double radius, double lumi[3], bool hq, int rangefn[lutSize - 1] = 1e-15f; // We need one temporary buffer - float ** buffer = allocArray (W, H); + const JaggedArray buffer (W, H); // the final result has to be in map // for an even number of levels that means: map => buffer, buffer => map @@ -157,23 +158,6 @@ void SHMap::update (Imagefloat* img, double radius, double lumi[3], bool hq, int } dirpyr_shmap(dirpyrlo[indx], dirpyrlo[1 - indx], W, H, rangefn, level, scale ); - - freeArray(buffer, H); - - //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - /* - // anti-alias filtering the result - #ifdef _OPENMP - #pragma omp for - #endif - for (int i=0; i0 && j>0 && i (W, H); + const JaggedArray buffer (W, H); // the final result has to be in map // for an even number of levels that means: map => buffer, buffer => map @@ -306,23 +290,6 @@ void SHMap::updateL (float** L, double radius, bool hq, int skip) } dirpyr_shmap(dirpyrlo[indx], dirpyrlo[1 - indx], W, H, rangefn, level, scale ); - - freeArray(buffer, H); - - //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - /* - // anti-alias filtering the result - #ifdef _OPENMP - #pragma omp for - #endif - for (int i=0; i0 && j>0 && itrimValue(pp->retinex.limd); highl->trimValue(pp->retinex.highl); baselog->trimValue(pp->retinex.baselog); -// grbl->trimValue(pp->retinex.grbl); gam->trimValue(pp->retinex.gam); slope->trimValue(pp->retinex.slope); highlights->trimValue(pp->retinex.highlights); shadows->trimValue(pp->retinex.shadows); + } void Retinex::updateCurveBackgroundHistogram (LUTu & histToneCurve, LUTu & histLCurve, LUTu & histCCurve,/* LUTu & histCLurve, LUTu & histLLCurve,*/ LUTu & histLCAM, LUTu & histCCAM, LUTu & histRed, LUTu & histGreen, LUTu & histBlue, LUTu & histLuma, LUTu & histLRETI) {