Bugfix for AMaZE.

This commit is contained in:
Emil Martinec
2010-07-17 16:28:50 -05:00
parent eced905191
commit a9f2594602

View File

@@ -28,18 +28,14 @@
void RawImageSource::amaze_demosaic_RT() { void RawImageSource::amaze_demosaic_RT() {
#undef MAX
#undef MIN
#undef CLIP
#define SQR(x) ((x)*(x)) #define SQR(x) ((x)*(x))
#define MIN(a,b) ((a) < (b) ? (a) : (b)) //#define MIN(a,b) ((a) < (b) ? (a) : (b))
#define MAX(a,b) ((a) > (b) ? (a) : (b)) //#define MAX(a,b) ((a) > (b) ? (a) : (b))
#define LIM(x,min,max) MAX(min,MIN(x,max)) #define LIM(x,min,max) MAX(min,MIN(x,max))
#define ULIM(x,y,z) ((y) < (z) ? LIM(x,y,z) : LIM(x,z,y)) #define ULIM(x,y,z) ((y) < (z) ? LIM(x,y,z) : LIM(x,z,y))
#define CLIP(x) LIM(x,0,65535) //#define CLIP(x) LIM(x,0,65535)
//allocate outpute arrays
int width=W, height=H; int width=W, height=H;
red = new unsigned short*[H]; red = new unsigned short*[H];
for (int i=0; i<H; i++) { for (int i=0; i<H; i++) {
@@ -58,62 +54,109 @@ void RawImageSource::amaze_demosaic_RT() {
// local variables // local variables
//position of top/left corner of the tile
int top, left; int top, left;
//offset of R pixel within a Bayer quartet
int ex, ey; int ex, ey;
//shifts of pointer value to access pixels in vertical and diagonal directions
static const int v1=TS, v2=2*TS, v3=3*TS, p1=-TS+1, p2=-2*TS+2, p3=-3*TS+3, m1=TS+1, m2=2*TS+2, m3=3*TS+3; static const int v1=TS, v2=2*TS, v3=3*TS, p1=-TS+1, p2=-2*TS+2, p3=-3*TS+3, m1=TS+1, m2=2*TS+2, m3=3*TS+3;
//neighborhood of a pixel
static const int nbr[5] = {-v2,-2,2,v2,0}; static const int nbr[5] = {-v2,-2,2,v2,0};
static const float eps=1e-5; //tolerance to avoid dividing by zero //tolerance to avoid dividing by zero
static const float epssq=1e-10; //tolerance to avoid dividing by zero static const float eps=1e-5, epssq=1e-10; //tolerance to avoid dividing by zero
//adaptive ratios threshold
static const float arthresh=0.75; static const float arthresh=0.75;
static const float nyqthresh=0.5;//0.5 //nyquist texture test threshold
static const float pmthresh=0.25;//0.25 static const float nyqthresh=0.5;
static const float lbd=1.0, ubd=1.0;//lbd=0.66, ubd=1.5; //diagonal interpolation test threshold
static const float pmthresh=0.25;
//factors for bounding interpolation in saturated regions
static const float lbd=1.0, ubd=1.0; //lbd=0.66, ubd=1.5 alternative values;
static const float gaussodd[4] = {0.14659727707323927f, 0.103592713382435f, 0.0732036125103057f, 0.0365543548389495f};//gaussian on 5x5 quincunx, sigma=1.2 //gaussian on 5x5 quincunx, sigma=1.2
static const float gaussodd[4] = {0.14659727707323927f, 0.103592713382435f, 0.0732036125103057f, 0.0365543548389495f};
//gaussian on 5x5, sigma=1.2
static const float gaussgrad[6] = {0.07384411893421103f, 0.06207511968171489f, 0.0521818194747806f, \ static const float gaussgrad[6] = {0.07384411893421103f, 0.06207511968171489f, 0.0521818194747806f, \
0.03687419286733595f, 0.03099732204057846f, 0.018413194161458882f};//gaussian on 5x5, sigma=1.2 0.03687419286733595f, 0.03099732204057846f, 0.018413194161458882f};
static const float gauss1[3] = {0.3376688223162362f, 0.12171198028231786f, 0.04387081413862306f};//gaussian on 3x3, sigma =0.7 //gaussian on 3x3, sigma =0.7
static const float gausseven[2] = {0.13719494435797422f, 0.05640252782101291f};//gaussian on 5x5 alt quincunx, sigma=1.5 static const float gauss1[3] = {0.3376688223162362f, 0.12171198028231786f, 0.04387081413862306f};
//gaussian on 5x5 alt quincunx, sigma=1.5
static const float gausseven[2] = {0.13719494435797422f, 0.05640252782101291f};
//guassian on quincunx grid
static const float gquinc[4] = {0.169917f, 0.108947f, 0.069855f, 0.0287182f}; static const float gquinc[4] = {0.169917f, 0.108947f, 0.069855f, 0.0287182f};
// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
char *buffer; // TS*TS*168 // beginning of storage block for tile
float (*rgb)[3]; // TS*TS*12 char *buffer;
float (*delh); // TS*TS*4 // rgb values
float (*delv); // TS*TS*4 float (*rgb)[3];
float (*delhsq); // TS*TS*4 // horizontal gradient
float (*delvsq); // TS*TS*4 float (*delh);
float (*dirwts)[2]; // TS*TS*8 // vertical gradient
float (*vcd); // TS*TS*4 float (*delv);
float (*hcd); // TS*TS*4 // square of delh
float (*vcdalt); // TS*TS*4 float (*delhsq);
float (*hcdalt); // TS*TS*4 // square of delv
float (*vcdsq); // TS*TS*4 float (*delvsq);
float (*hcdsq); // TS*TS*4 // gradient based directional weights for interpolation
float (*cddiffsq); // TS*TS*4 float (*dirwts)[2];
float (*hvwt); // TS*TS*4 // vertically interpolated color differences G-R, G-B
float (*Dgrb)[2]; // TS*TS*8 float (*vcd);
float (*delp); // TS*TS*4 // horizontally interpolated color differences
float (*delm); // TS*TS*4 float (*hcd);
float (*rbint); // TS*TS*4 // alternative vertical interpolation
float (*Dgrbh2); // TS*TS*4 float (*vcdalt);
float (*Dgrbv2); // TS*TS*4 // alternative horizontal interpolation
float (*dgintv); // TS*TS*4 float (*hcdalt);
float (*dginth); // TS*TS*4 // square of vcd
float (*Dgrbp1); // TS*TS*4 float (*vcdsq);
float (*Dgrbm1); // TS*TS*4 // square of hcd
float (*Dgrbpsq1); // TS*TS*4 float (*hcdsq);
float (*Dgrbmsq1); // TS*TS*4 // square of average color difference
float (*cfa); // TS*TS*4 float (*cddiffsq);
float (*pmwt); // TS*TS*4 // weight to give horizontal vs vertical interpolation
float (*rbp); // TS*TS*4 float (*hvwt);
float (*rbm); // TS*TS*4 // final interpolated color difference
float (*Dgrb)[2];
// gradient in plus (NE/SW) direction
float (*delp);
// gradient in minus (NW/SE) direction
float (*delm);
// diagonal interpolation of R+B
float (*rbint);
// horizontal curvature of interpolated G (used to refine interpolation in Nyquist texture regions)
float (*Dgrbh2);
// vertical curvature of interpolated G
float (*Dgrbv2);
// difference between up/down interpolations of G
float (*dgintv);
// difference between left/right interpolations of G
float (*dginth);
// diagonal (plus) color difference R-B or G1-G2
float (*Dgrbp1);
// diagonal (minus) color difference R-B or G1-G2
float (*Dgrbm1);
// square of diagonal color difference
float (*Dgrbpsq1);
// square of diagonal color difference
float (*Dgrbmsq1);
// tile raw data
float (*cfa);
// relative weight for combining plus and minus diagonal interpolations
float (*pmwt);
// interpolated color difference R-B in plus direction
float (*rbp);
// interpolated color difference R-B in minus direction
float (*rbm);
int (*nyquist); // TS*TS*4 // nyquist texture flag 1=nyquist, 0=not nyquist
int (*nyquist);
// assign working space // assign working space
@@ -121,7 +164,7 @@ void RawImageSource::amaze_demosaic_RT() {
//merror(buffer,"amaze_interpolate()"); //merror(buffer,"amaze_interpolate()");
memset(buffer,0,(34*sizeof(float)+sizeof(int))*TS*TS); memset(buffer,0,(34*sizeof(float)+sizeof(int))*TS*TS);
// rgb array // rgb array
rgb = (float (*)[3]) buffer; //pointers to array rgb = (float (*)[3]) buffer; //pointers to array
delh = (float (*)) (buffer + 3*sizeof(float)*TS*TS); delh = (float (*)) (buffer + 3*sizeof(float)*TS*TS);
delv = (float (*)) (buffer + 4*sizeof(float)*TS*TS); delv = (float (*)) (buffer + 4*sizeof(float)*TS*TS);
delhsq = (float (*)) (buffer + 5*sizeof(float)*TS*TS); delhsq = (float (*)) (buffer + 5*sizeof(float)*TS*TS);
@@ -195,30 +238,81 @@ void RawImageSource::amaze_demosaic_RT() {
//code is openmp ready; just have to pull local tile variable declarations inside the tile loop //code is openmp ready; just have to pull local tile variable declarations inside the tile loop
for (top=-16; top < height; top += TS-32) for (top=-16; top < height; top += TS-32)
for (left=-16; left < width; left += TS-32) { for (left=-16; left < width; left += TS-32) {
//location of tile bottom edge
int bottom = MIN( top+TS,height+16); int bottom = MIN( top+TS,height+16);
//location of tile right edge
int right = MIN(left+TS, width+16); int right = MIN(left+TS, width+16);
//tile width (=TS except for right edge of image)
int rr1 = bottom - top; int rr1 = bottom - top;
//tile height (=TS except for bottom edge of image)
int cc1 = right - left; int cc1 = right - left;
//tile vars //tile vars
//counters for pixel location in the image
int row, col; int row, col;
//min and max row/column in the tile
int rrmin, rrmax, ccmin, ccmax; int rrmin, rrmax, ccmin, ccmax;
int rr, cc, c, indx, indx1, dir, i, j, sgn; //counters for pixel location within the tile
int rr, cc;
float crse, crnw, crne, crsw, rbse, rbnw, rbne, rbsw, wtse, wtnw, wtsw, wtne; //color index 0=R, 1=G, 2=B
float pmwtalt; int c;
//pointer counters within the tile
int indx, indx1;
//direction counter for nbrs[]
int dir;
//dummy indices
int i, j;
// +1 or -1
int sgn;
//color ratios in up/down/left/right directions
float cru, crd, crl, crr; float cru, crd, crl, crr;
float vwt, hwt, pwt, mwt, Gintv, Ginth; //adaptive weights for vertical/horizontal/plus/minus directions
float guar, gdar, glar, grar, guha, gdha, glha, grha, Ginthar, Ginthha, Gintvar, Gintvha, hcdaltvar, vcdaltvar; float vwt, hwt, pwt, mwt;
//vertical and horizontal G interpolations
float Gintv, Ginth;
//G interpolated in vert/hor directions using adaptive ratios
float guar, gdar, glar, grar;
//G interpolated in vert/hor directions using Hamilton-Adams method
float guha, gdha, glha, grha;
//interpolated G from fusing left/right or up/down
float Ginthar, Ginthha, Gintvar, Gintvha;
//color difference (G-R or G-B) variance in up/down/left/right directions
float Dgrbvvaru, Dgrbvvard, Dgrbhvarl, Dgrbhvarr; float Dgrbvvaru, Dgrbvvard, Dgrbhvarl, Dgrbhvarr;
float sumh, sumv, sumsqh, sumsqv, areawt; //gradients in various directions
float nyqtest, vcdvar, hcdvar, hvwtalt, vo, ve, gradp, gradm, gradv, gradh, gradpm, gradhv; float gradp, gradm, gradv, gradh, gradpm, gradhv;
float vcdvar1, hcdvar1, varwt, diffwt; //color difference variances in vertical and horizontal directions
float rbvarp, rbvarm; float vcdvar, hcdvar, vcdvar1, hcdvar1, hcdaltvar, vcdaltvar;
//adaptive interpolation weight using variance of color differences
float varwt;
//adaptive interpolation weight using difference of left-right and up-down G interpolations
float diffwt;
//alternative adaptive weight for combining horizontal/vertical interpolations
float hvwtalt;
//temporary variables for combining interpolation weights at R and B sites
float vo, ve;
//interpolation of G in four directions
float gu, gd, gl, gr; float gu, gd, gl, gr;
//variance of G in vertical/horizontal directions
float gvarh, gvarv; float gvarh, gvarv;
//Nyquist texture test
float nyqtest;
//accumulators for Nyquist texture interpolation
float sumh, sumv, sumsqh, sumsqv, areawt;
//color ratios in diagonal directions
float crse, crnw, crne, crsw;
//color differences in diagonal directions
float rbse, rbnw, rbne, rbsw;
//adaptive weights for combining diagonal interpolations
float wtse, wtnw, wtsw, wtne;
//alternate weight for combining diagonal interpolations
float pmwtalt;
//variance of R-B in plus/minus directions
float rbvarp, rbvarm;
// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
/* /*
char *buffer; // TS*TS*168 char *buffer; // TS*TS*168
@@ -300,6 +394,9 @@ void RawImageSource::amaze_demosaic_RT() {
// rgb from input CFA data // rgb from input CFA data
// rgb values should be floating point number between 0 and 1 // rgb values should be floating point number between 0 and 1
// after white balance multipliers are applied // after white balance multipliers are applied
// a 16 pixel border is added to each side of the image
// bookkeeping for borders
if (top<0) {rrmin=16;} else {rrmin=0;} if (top<0) {rrmin=16;} else {rrmin=0;}
if (left<0) {ccmin=16;} else {ccmin=0;} if (left<0) {ccmin=16;} else {ccmin=0;}
if (bottom>height) {rrmax=height-top;} else {rrmax=rr1;} if (bottom>height) {rrmax=height-top;} else {rrmax=rr1;}
@@ -442,10 +539,10 @@ void RawImageSource::amaze_demosaic_RT() {
rbint[indx]=0; rbint[indx]=0;
//color ratios in each cardinal direction //color ratios in each cardinal direction
cru = cfa[indx-v1]*(dirwts[indx-v2][0]+dirwts[indx][0])/(dirwts[indx-v2][0]*cfa[indx]+dirwts[indx][0]*cfa[indx-v2]); cru = cfa[indx-v1]*(dirwts[indx-v2][0]+dirwts[indx][0])/(eps+dirwts[indx-v2][0]*cfa[indx]+dirwts[indx][0]*cfa[indx-v2]);
crd = cfa[indx+v1]*(dirwts[indx+v2][0]+dirwts[indx][0])/(dirwts[indx+v2][0]*cfa[indx]+dirwts[indx][0]*cfa[indx+v2]); crd = cfa[indx+v1]*(dirwts[indx+v2][0]+dirwts[indx][0])/(eps+dirwts[indx+v2][0]*cfa[indx]+dirwts[indx][0]*cfa[indx+v2]);
crl = cfa[indx-1]*(dirwts[indx-2][1]+dirwts[indx][1])/(dirwts[indx-2][1]*cfa[indx]+dirwts[indx][1]*cfa[indx-2]); crl = cfa[indx-1]*(dirwts[indx-2][1]+dirwts[indx][1])/(eps+dirwts[indx-2][1]*cfa[indx]+dirwts[indx][1]*cfa[indx-2]);
crr = cfa[indx+1]*(dirwts[indx+2][1]+dirwts[indx][1])/(dirwts[indx+2][1]*cfa[indx]+dirwts[indx][1]*cfa[indx+2]); crr = cfa[indx+1]*(dirwts[indx+2][1]+dirwts[indx][1])/(eps+dirwts[indx+2][1]*cfa[indx]+dirwts[indx][1]*cfa[indx+2]);
guha=cfa[indx-v1]+0.5*(cfa[indx]-cfa[indx-v2]); guha=cfa[indx-v1]+0.5*(cfa[indx]-cfa[indx-v2]);
gdha=cfa[indx+v1]+0.5*(cfa[indx]-cfa[indx+v2]); gdha=cfa[indx+v1]+0.5*(cfa[indx]-cfa[indx+v2]);