From 77914ebbf8f9d362c8eed49497f81ed3908bcc2a Mon Sep 17 00:00:00 2001 From: Emil Martinec Date: Thu, 3 Jun 2010 22:41:24 -0500 Subject: [PATCH] CA auto correction routine --- CMakeLists.txt | 12 +- rtdata/languages/English (UK) | 2 + rtdata/languages/English (US) | 2 + rtengine/CA_correct_RT.cc | 678 ++++++++++++++++++++++++++++++++++ rtengine/rawimagesource.cc | 21 +- rtengine/rawimagesource.h | 3 + rtengine/settings.h | 1 + rtgui/options.cc | 2 + rtgui/preferences.cc | 5 + rtgui/preferences.h | 1 + 10 files changed, 718 insertions(+), 9 deletions(-) create mode 100644 rtengine/CA_correct_RT.cc diff --git a/CMakeLists.txt b/CMakeLists.txt index 8cd43729f..7152ca841 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -117,12 +117,12 @@ if (WITH_RAWZOR) endif (WIN32) endif (WITH_RAWZOR) -if (OPTION_OMP) - find_package(OpenMP) - if (OPENMP_FOUND) - set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}") - endif (OPENMP_FOUND) -endif (OPTION_OMP) +#if (OPTION_OMP) +# find_package(OpenMP) +# if (OPENMP_FOUND) +# set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}") +# endif (OPENMP_FOUND) +#endif (OPTION_OMP) add_subdirectory (rtexif) add_subdirectory (rtengine) diff --git a/rtdata/languages/English (UK) b/rtdata/languages/English (UK) index d992de66b..37b52708d 100644 --- a/rtdata/languages/English (UK) +++ b/rtdata/languages/English (UK) @@ -363,6 +363,8 @@ PREFERENCES_DATEFORMAT;Date Format PREFERENCES_DATEFORMATHINT;You can use the following formatting strings:\n%y : year\n%m : month\n%d : day\n\nFor example, the hungarian date format is:\n%y/%m/%d PREFERENCES_DCBENHANCE;Apply DCB enhancment step PREFERENCES_DCBITERATIONS;Number of DCB iterations +#Emil's CA autocorrection +PREFERENCES_CACORRECTION;Apply CA auto correction PREFERENCES_DEFAULTLANG;Default language PREFERENCES_DEFAULTTHEME;Default theme PREFERENCES_DEMOSAICINGALGO;Demosaicing Algorithm diff --git a/rtdata/languages/English (US) b/rtdata/languages/English (US) index b87063e3e..492a1b35c 100644 --- a/rtdata/languages/English (US) +++ b/rtdata/languages/English (US) @@ -365,6 +365,8 @@ PREFERENCES_DATEFORMAT;Date Format PREFERENCES_DATEFORMATHINT;You can use the following formatting strings:\n%y : year\n%m : month\n%d : day\n\nFor example, the hungarian date format is:\n%y/%m/%d PREFERENCES_DCBENHANCE;Apply DCB enhancement step PREFERENCES_DCBITERATIONS;Number of DCB iterations +#Emil's CA autocorrection +PREFERENCES_CACORRECTION;Apply CA auto correction PREFERENCES_DEFAULTLANG;Default language PREFERENCES_DEFAULTTHEME;Default theme PREFERENCES_DEMOSAICINGALGO;Demosaicing Algorithm diff --git a/rtengine/CA_correct_RT.cc b/rtengine/CA_correct_RT.cc new file mode 100644 index 000000000..820edab40 --- /dev/null +++ b/rtengine/CA_correct_RT.cc @@ -0,0 +1,678 @@ +//////////////////////////////////////////////////////////////// +// +// Chromatic Aberration Auto-correction +// +// copyright (c) 2008-2010 Emil Martinec +// +// +// code dated: June 2, 2010 +// +// CA_correct_RT.cc 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. +// +// This program 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 this program. If not, see . +// +//////////////////////////////////////////////////////////////// + + + +//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +int RawImageSource::LinEqSolve(int c, int dir, int nDim, float* pfMatr, float* pfVect, float* pfSolution) +{ +//============================================================================== +// return 1 if system not solving, 0 if system solved +// nDim - system dimension +// pfMatr - matrix with coefficients +// pfVect - vector with free members +// pfSolution - vector with system solution +// pfMatr becames trianglular after function call +// pfVect changes after function call +// +// Developer: Henry Guennadi Levkin +// +//============================================================================== + + float fMaxElem; + float fAcc; + + int i, j, k, m; + + for(k=0; k<(nDim-1); k++) {// base row of matrix + // search of line with max element + fMaxElem = fabs( pfMatr[k*nDim + k] ); + m = k; + for (i=k+1; i=0; k--) { + pfSolution[k] = pfVect[k]; + for(i=(k+1); i(b)) {temp=(a);(a)=(b);(b)=temp;} } + #define SQR(x) ((x)*(x)) + + // local variables + int width=W, height=H; + float (*Gtmp); + Gtmp = (float (*)) calloc ((height)*(width), sizeof *Gtmp); + + int polyord=4, numpar=16, numblox[3]={0,0,0}; + + //static const int border=8; + int rrmin, rrmax, ccmin, ccmax; + int top, bottom, left, right, row, col; + int rr, cc, rr1, cc1, c, indx, indx1, i, j, k, m, n, dir; + int GRBdir[2][3], offset[2][3]; + int shifthfloor[3], shiftvfloor[3], shifthceil[3], shiftvceil[3]; + int vblsz, hblsz, vblock, hblock; + //int verbose=1; + int res; + int v1=TS, v2=2*TS, v3=3*TS, v4=4*TS;//, p1=-TS+1, p2=-2*TS+2, p3=-3*TS+3, m1=TS+1, m2=2*TS+2, m3=3*TS+3; + + float eps=1e-10; //tolerance to avoid dividing by zero + + float wtu, wtd, wtl, wtr; + float coeff[2][3][3], CAshift[2][3]; + float polymat[3][2][1296], shiftmat[3][2][36], fitparams[3][2][36]; + float shifthfrac[3], shiftvfrac[3], temp, p[9]; + float gdiff, deltgrb, denom, Ginthfloor, Ginthceil, Gint, gradwt; + float grbdiffinthfloor, grbdiffinthceil, grbdiffint; + float blockgave, blockgsqave; + + static const float gaussg[5] = {0.171582, 0.15839, 0.124594, 0.083518, 0.0477063}; + static const float gaussrb[3] = {0.332406, 0.241376, 0.0924212}; + + //char *buffer1; // vblsz*hblsz*3*2 + //float (*blockshifts)[3][2]; // vblsz*hblsz*3*2 + float blockshifts[10000][3][2]; //fixed memory allocation + float blockwt[10000]; //fixed memory allocation + + + char *buffer; // TS*TS*16 + float (*rgb)[3]; // TS*TS*12 + float (*grbdiff); // TS*TS*4 + float (*gshift); // TS*TS*4 + float (*ghpfh); // TS*TS*4 + float (*ghpfv); // TS*TS*4 + float (*rbhpfh); // TS*TS*4 + float (*rbhpfv); // TS*TS*4 + + + /* assign working space; this would not be necessary + if the algorithm is part of the larger pre-interpolation processing */ + buffer = (char *) malloc(36*TS*TS); + //merror(buffer,"CA_correct()"); + memset(buffer,0,36*TS*TS); + + // rgb array + rgb = (float (*)[3]) buffer; + grbdiff = (float (*)) (buffer + 12*TS*TS); + gshift = (float (*)) (buffer + 16*TS*TS); + ghpfh = (float (*)) (buffer + 20*TS*TS); + ghpfv = (float (*)) (buffer + 24*TS*TS); + rbhpfh = (float (*)) (buffer + 28*TS*TS); + rbhpfv = (float (*)) (buffer + 32*TS*TS); + + + vblsz=ceil((float)(height+border2)/(TS-border2))+2; + hblsz=ceil((float)(width+border2)/(TS-border2))+2; + + /*buffer1 = (char *) malloc(4*vblsz*hblsz*3*2); + merror(buffer1,"CA_correct()"); + memset(buffer1,0,4*vblsz*hblsz*3*2); + // block CA shifts + blockshifts = (float (*)[3][2]) buffer1;*/ + + + + + // Main algorithm: Tile loop + //#pragma omp parallel for shared(image,height,width) private(top,left,indx,indx1) schedule(dynamic) + for (top=-border, vblock=1; top < height; top += TS-border2, vblock++) + for (left=-border, hblock=1; left < width; left += TS-border2, hblock++) { + bottom = MIN( top+TS,height+border); + right = MIN(left+TS, width+border); + rr1 = bottom - top; + cc1 = right - left; + //t1_init = clock(); + // rgb from input CFA data + // rgb values should be floating point number between 0 and 1 + // after white balance multipliers are applied + if (top<0) {rrmin=border;} else {rrmin=0;} + if (left<0) {ccmin=border;} else {ccmin=0;} + if (bottom>height) {rrmax=height-top;} else {rrmax=rr1;} + if (right>width) {ccmax=width-left;} else {ccmax=cc1;} + + blockgave=blockgsqave=denom=0; + for (rr=rrmin; rr < rrmax; rr++) + for (row=rr+top, cc=ccmin; cc < ccmax; cc++) { + col = cc+left; + c = FC(rr,cc); + indx=row*width+col; + indx1=rr*TS+cc; + rgb[indx1][c] = (ri->data[row][col])/65535.0f; + } + blockwt[vblock*hblsz+hblock]=0; + //blockwt[vblock*hblsz+hblock]=1; + + // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + //fill borders + if (rrmin>0) { + for (rr=0; rrdata[(height-rr-2)][left+cc])/65535.0f; + } + } + if (ccmin>0) { + for (rr=rrmin; rrdata[(top+rr)][(width-cc-2)])/65535.0f; + } + } + + //also, fill the image corners + if (rrmin>0 && ccmin>0) { + for (rr=0; rrdata[border2-rr][border2-cc])/65535.0f; + } + } + if (rrmaxdata[(height-rr-2)][(width-cc-2)])/65535.0f; + } + } + if (rrmin>0 && ccmaxdata[(border2-rr)][(width-cc-2)])/65535.0f; + } + } + if (rrmax0) { + for (rr=0; rrdata[(height-rr-2)][(border2-cc)])/65535.0f; + } + } + + //end of border fill + // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + for (j=0; j<2; j++) + for (k=0; k<3; k++) + for (c=0; c<3; c+=2) coeff[j][k][c]=0; + + //end of initialization + + + for (rr=3; rr < rr1-3; rr++) + for (row=rr+top, cc=3, indx=rr*TS+cc; cc < cc1-3; cc++, indx++) { + col = cc+left; + c = FC(rr,cc); + + if (c!=1) { + //compute directional weights using image gradients + wtu=1/SQR(eps+fabs(rgb[(rr+1)*TS+cc][1]-rgb[(rr-1)*TS+cc][1])+fabs(rgb[(rr)*TS+cc][c]-rgb[(rr-2)*TS+cc][c])+fabs(rgb[(rr-1)*TS+cc][1]-rgb[(rr-3)*TS+cc][1])); + wtd=1/SQR(eps+fabs(rgb[(rr-1)*TS+cc][1]-rgb[(rr+1)*TS+cc][1])+fabs(rgb[(rr)*TS+cc][c]-rgb[(rr+2)*TS+cc][c])+fabs(rgb[(rr+1)*TS+cc][1]-rgb[(rr+3)*TS+cc][1])); + wtl=1/SQR(eps+fabs(rgb[(rr)*TS+cc+1][1]-rgb[(rr)*TS+cc-1][1])+fabs(rgb[(rr)*TS+cc][c]-rgb[(rr)*TS+cc-2][c])+fabs(rgb[(rr)*TS+cc-1][1]-rgb[(rr)*TS+cc-3][1])); + wtr=1/SQR(eps+fabs(rgb[(rr)*TS+cc-1][1]-rgb[(rr)*TS+cc+1][1])+fabs(rgb[(rr)*TS+cc][c]-rgb[(rr)*TS+cc+2][c])+fabs(rgb[(rr)*TS+cc+1][1]-rgb[(rr)*TS+cc+3][1])); + + + //store in rgb array the interpolated G value at R/B grid points using directional weighted average + rgb[indx][1]=(wtu*rgb[indx-v1][1]+wtd*rgb[indx+v1][1]+wtl*rgb[indx-1][1]+wtr*rgb[indx+1][1])/(wtu+wtd+wtl+wtr); + } + if (row>-1 && row-1 && col1.0 || fabs(blockshifts[(vblock)*hblsz+hblock][c][1])>1.0) continue; + numblox[c] += 1; + for (dir=0; dir<2; dir++) { + for (i=0; iheight) {rrmax=height-top;} else {rrmax=rr1;} + if (right>width) {ccmax=width-left;} else {ccmax=cc1;} + + for (rr=rrmin; rr < rrmax; rr++) + for (row=rr+top, cc=ccmin; cc < ccmax; cc++) { + col = cc+left; + c = FC(rr,cc); + indx=row*width+col; + indx1=rr*TS+cc; + //rgb[indx1][c] = image[indx][c]/65535.0f; + rgb[indx1][c] = (ri->data[row][col])/65535.0f; + if ((c&1)==0) rgb[indx1][1] = Gtmp[indx]; + } + // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + //fill borders + if (rrmin>0) { + for (rr=0; rrdata[(height-rr-2)][left+cc])/65535.0f; + rgb[(rrmax+rr)*TS+cc][1] = Gtmp[(height-rr-2)*width+left+cc]; + } + } + if (ccmin>0) { + for (rr=rrmin; rrdata[(top+rr)][(width-cc-2)])/65535.0f; + rgb[rr*TS+ccmax+cc][1] = Gtmp[(top+rr)*width+(width-cc-2)]; + } + } + + //also, fill the image corners + if (rrmin>0 && ccmin>0) { + for (rr=0; rrdata[border2-rr][border2-cc])/65535.0f; + rgb[(rr)*TS+cc][1] = Gtmp[(border2-rr)*width+border2-cc]; + } + } + if (rrmaxdata[(height-rr-2)][(width-cc-2)])/65535.0f; + rgb[(rrmax+rr)*TS+ccmax+cc][1] = Gtmp[(height-rr-2)*width+(width-cc-2)]; + } + } + if (rrmin>0 && ccmaxdata[(border2-rr)][(width-cc-2)])/65535.0f; + rgb[(rr)*TS+ccmax+cc][1] = Gtmp[(border2-rr)*width+(width-cc-2)]; + } + } + if (rrmax0) { + for (rr=0; rrdata[(height-rr-2)][(border2-cc)])/65535.0f; + rgb[(rrmax+rr)*TS+cc][1] = Gtmp[(height-rr-2)*width+(border2-cc)]; + } + } + + //end of border fill + // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + + blockshifts[(vblock)*hblsz+hblock][0][0] = blockshifts[(vblock)*hblsz+hblock][0][1] = 0; + blockshifts[(vblock)*hblsz+hblock][2][0] = blockshifts[(vblock)*hblsz+hblock][2][1] = 0; + for (i=0; idata[row][col] = CLIP((int)(65535.0f*rgb[(rr)*TS+cc][c] + 0.5f)); + + } + } + + // clean up + free(buffer); + free(Gtmp); + //free(buffer1); + + + +#undef TS +#undef polyord +#undef numpar +#undef border +#undef PIX_SORT +#undef SQR + +} + + + diff --git a/rtengine/rawimagesource.cc b/rtengine/rawimagesource.cc index 4f6c23d36..8ef9b8a70 100644 --- a/rtengine/rawimagesource.cc +++ b/rtengine/rawimagesource.cc @@ -772,6 +772,20 @@ int RawImageSource::load (Glib::ustring fname) { delete corr_data; */ } + + +//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +//Emil's CA auto correction + + if (settings->ca_autocorrect) { + if (plistener) { + plistener->setProgressStr ("CA Auto Correction..."); + plistener->setProgress (0.0); + } + + CA_correct_RT(); + } +//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% if (ri->filters) { // demosaic @@ -3127,10 +3141,11 @@ void RawImageSource::dcb_demosaic(int iterations, int dcb_enhance) free(image); } - //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - //Emil's code for AMaZE +//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +//Emil's code for AMaZE #include "amaze_interpolate_RT.cc"//AMaZE demosaic - //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +#include "CA_correct_RT.cc"//Emil's CA auto correction +//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% } /* namespace */ diff --git a/rtengine/rawimagesource.h b/rtengine/rawimagesource.h index 0a9c7d31e..bab0bafd4 100644 --- a/rtengine/rawimagesource.h +++ b/rtengine/rawimagesource.h @@ -135,6 +135,9 @@ class RawImageSource : public ImageSource { inline void interpolate_row_rb (unsigned short* ar, unsigned short* ab, unsigned short* pg, unsigned short* cg, unsigned short* ng, int i); inline void interpolate_row_rb_mul_pp (unsigned short* ar, unsigned short* ab, unsigned short* pg, unsigned short* cg, unsigned short* ng, int i, double r_mul, double g_mul, double b_mul, int x1, int width, int skip); + int LinEqSolve (int c, int dir, int nDim, float* pfMatr, float* pfVect, float* pfSolution);//Emil's CA auto correction + void CA_correct_RT (); + void eahd_demosaic (); void hphd_demosaic (); void vng4_demosaic (); diff --git a/rtengine/settings.h b/rtengine/settings.h index bfe53173e..cd2dffa07 100644 --- a/rtengine/settings.h +++ b/rtengine/settings.h @@ -33,6 +33,7 @@ namespace rtengine { bool verbose; int dcb_iterations; // number of dcb iterations bool dcb_enhance; // whether to do image refinment + bool ca_autocorrect; // Emil's auto CA correction /** Creates a new instance of Settings. * @return a pointer to the new Settings instance. */ diff --git a/rtgui/options.cc b/rtgui/options.cc index 01d7e36b8..c28542616 100644 --- a/rtgui/options.cc +++ b/rtgui/options.cc @@ -120,6 +120,7 @@ void Options::setDefaults () { rtSettings.dualThreadEnabled = true; rtSettings.demosaicMethod = "amaze";//Emil's code for AMaZE + rtSettings.ca_autocorrect = false;//Emil's CA correction rtSettings.colorCorrectionSteps = 0; rtSettings.dcb_iterations = 2; rtSettings.dcb_enhance = true; @@ -384,6 +385,7 @@ int Options::saveToFile (Glib::ustring fname) { keyFile.set_integer ("Algorithms", "ColorCorrection", rtSettings.colorCorrectionSteps); keyFile.set_integer ("Algorithms", "DCBIterations", rtSettings.dcb_iterations); keyFile.set_boolean ("Algorithms", "DCBEnhance", rtSettings.dcb_enhance); + keyFile.set_boolean ("Algorithms", "CACorrect", rtSettings.ca_autocorrect);//Emil's CA correction keyFile.set_integer ("Crop Settings", "DPI", cropDPI); diff --git a/rtgui/preferences.cc b/rtgui/preferences.cc index acf7ffd3b..1a3ed53b0 100644 --- a/rtgui/preferences.cc +++ b/rtgui/preferences.cc @@ -286,11 +286,14 @@ Gtk::Widget* Preferences::getProcParamsPanel () { hb13->pack_start (*dcbIterations); dcbEnhance = Gtk::manage(new Gtk::CheckButton((M("PREFERENCES_DCBENHANCE")))); + + caAutoCorrect = Gtk::manage(new Gtk::CheckButton((M("PREFERENCES_CACORRECTION"))));//Emil's CA correction fdb->pack_start (*hb11, Gtk::PACK_SHRINK, 4); fdb->pack_start (*hb12, Gtk::PACK_SHRINK, 4); fdb->pack_start (*hb13, Gtk::PACK_SHRINK, 4); fdb->pack_start (*dcbEnhance, Gtk::PACK_SHRINK, 4); + fdb->pack_start (*caAutoCorrect, Gtk::PACK_SHRINK, 4);//Emil's CA correction mvbpp->pack_start (*fdem, Gtk::PACK_SHRINK, 4); mvbpp->set_border_width (4); // drlab->set_size_request (drimg->get_width(), -1); @@ -718,6 +721,7 @@ void Preferences::storePreferences () { moptions.rtSettings.demosaicMethod = "ahd"; moptions.rtSettings.dcb_iterations=(int)dcbIterations->get_value(); moptions.rtSettings.dcb_enhance=dcbEnhance->get_active(); + moptions.rtSettings.ca_autocorrect=caAutoCorrect->get_active();//Emil's CA correction if (sdcurrent->get_active ()) moptions.startupDir = STARTUPDIR_CURRENT; @@ -816,6 +820,7 @@ void Preferences::fillPreferences () { dcbEnhance->set_sensitive(moptions.rtSettings.demosaicMethod=="dcb"); dcbIterations->set_sensitive(moptions.rtSettings.demosaicMethod=="dcb"); dcbIterationsLabel->set_sensitive(moptions.rtSettings.demosaicMethod=="dcb"); + caAutoCorrect->set_active(moptions.rtSettings.ca_autocorrect);//Emil's CA Auto Correction if (moptions.startupDir==STARTUPDIR_CURRENT) sdcurrent->set_active (); diff --git a/rtgui/preferences.h b/rtgui/preferences.h index ebf2d9956..b721cb01e 100644 --- a/rtgui/preferences.h +++ b/rtgui/preferences.h @@ -75,6 +75,7 @@ class Preferences : public Gtk::Dialog { Gtk::Label* dcbIterationsLabel; Gtk::SpinButton* dcbIterations; Gtk::CheckButton* dcbEnhance; + Gtk::CheckButton* caAutoCorrect;//Emil's CA correction Gtk::FileChooserButton* iccDir; Gtk::FileChooserButton* monProfile;