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

Speedup and reduced memory usage for lmmse demosaic, Issue 2665

Browse Source
This commit is contained in:
Ingo
2015-02-16 14:08:55 +01:00
parent ec2f378bbf
commit 9df2a5d8ed
1 changed files with 258 additions and 204 deletions
Download Patch File Download Diff File Expand all files Collapse all files

462
rtengine/demosaic_algos.cc
View File

@@ -1190,11 +1190,8 @@ void RawImageSource::lmmse_interpolate_omp(int winw, int winh, int iterations)
{
int c, ii;
float h0, h1, h2, h3, h4, hs;
float (*rix)[6];
float (*qix)[6];
char *buffer;
const int width=winw, height=winh;
const int ba = 10;
const int rr1 = height + 2*ba;
@@ -1226,24 +1223,38 @@ void RawImageSource::lmmse_interpolate_omp(int winw, int winh, int iterations)
plistener->setProgressStr (Glib::ustring::compose(M("TP_RAW_DMETHOD_PROGRESSBAR"), RAWParams::BayerSensor::methodstring[RAWParams::BayerSensor::lmmse]));
plistener->setProgress (0.0);
}
float (*image)[3];
float maxdata=0.f;
image = (float (*)[3]) calloc (width*height, sizeof *image);
float *imageBuffer = (float*)malloc (((((width+1)/2)*height) + (((height+1)/2)*(width+1))) * sizeof *imageBuffer);
float *image[3];
image[0] = imageBuffer;
image[1] = imageBuffer + ((height+1)/2)*((width+1)/2);
image[2] = imageBuffer + ((height+1)/2)*((width+1)/2) + (((width+1)/2)*height);
unsigned int a=0;
float *rix[6];
float *qix[6];
float *buffer = (float *)calloc(rr1*cc1*6*sizeof(float),1);
if(buffer == NULL) { // allocation of big block of memory failed, try to get 6 smaller ones
printf("lmmse_interpolate_omp: allocation of big memory block failed, try to get 6 smaller ones now...\n");
for(int i=0;i<6;i++)
qix[i] = (float *)calloc(rr1*cc1*sizeof(float),1);
} else {
qix[0] = buffer;
for(int i=1;i<6;i++)
qix[i] = qix[i-1] + rr1*cc1;
}
#ifdef _OPENMP
#pragma omp parallel for
#endif
for (int ii=0; ii<height; ii++)
for (int jj=0; jj<width; jj++) {
image[ii*width+jj][fc(ii,jj)] = CLIP(rawData[ii][jj]);
for (int row=0; row<height; row++)
for (int col=0; col<width; col++) {
int c = fc(row,col);
image[c][(((row*width)>>(1-(c&1)))+col)>>1] = CLIP(rawData[row][col]);
}
maxdata=65535.f;
if (applyGamma)
buffer = (char *)calloc(rr1*cc1*6*sizeof(float)+(int)(maxdata)*sizeof(float),1);
else
buffer = (char *)calloc(rr1*cc1*6*sizeof(float),1);
qix = (float (*)[6])buffer;
{
if (plistener) plistener->setProgress (0.1);
}
@@ -1255,18 +1266,15 @@ void RawImageSource::lmmse_interpolate_omp(int winw, int winh, int iterations)
#ifdef _OPENMP
#pragma omp for
#endif
for (int rrr=0; rrr < rr1; rrr++)
for (int ccc=0, row=rrr-ba; ccc < cc1; ccc++) {
for (int rrr=ba; rrr < rr1-ba; rrr++)
for (int ccc=ba, row=rrr-ba; ccc < cc1-ba; ccc++) {
int col = ccc - ba;
rix = qix + rrr*cc1 + ccc;
if ((row >= 0) & (row < height) & (col >= 0) & (col < width)) {
rix[4] = qix[4] + rrr*cc1 + ccc;
int c = fc(row,col);
if (applyGamma)
rix[0][4] = Color::gammatab_24_17a[image[row*width+col][FC(row,col)]];
rix[4][0] = Color::gammatab_24_17a[image[c][(((row*width)>>(1-(c&1)))+col)>>1]];
else
rix[0][4] = (float)image[row*width+col][FC(row,col)]/65535.0f;
}
else
rix[0][4] = 0.f;
rix[4][0] = (float)image[c][(((row*width)>>(1-(c&1)))+col)>>1]/65535.0f;
}
#ifdef _OPENMP
#pragma omp single
@@ -1278,37 +1286,41 @@ void RawImageSource::lmmse_interpolate_omp(int winw, int winh, int iterations)
// G-R(B)
#ifdef _OPENMP
#pragma omp for
#pragma omp for schedule(dynamic,16)
#endif
for (int rr=2; rr < rr1-2; rr++) {
// G-R(B) at R(B) location
for (int cc=2+(FC(rr,2)&1); cc < cc1-2; cc+=2) {
rix = qix + rr*cc1 + cc;
float v0 = x00625(rix[-w1-1][4]+rix[-w1+1][4]+rix[w1-1][4]+rix[w1+1][4]) +x0250(rix[0][4]);
rix[4] = qix[4] + rr*cc1 + cc;
float v0 = x00625(rix[4][-w1-1]+rix[4][-w1+1]+rix[4][w1-1]+rix[4][w1+1]) +x0250(rix[4][0]);
// horizontal
rix[0][0] = - x0250(rix[ -2][4] + rix[ 2][4])+ xdiv2f(rix[ -1][4] + rix[0][4] + rix[ 1][4]);
rix[0] = qix[0] + rr*cc1 + cc;
rix[0][0] = - x0250(rix[4][ -2] + rix[4][ 2])+ xdiv2f(rix[4][ -1] + rix[4][0] + rix[4][ 1]);
float Y = v0 + xdiv2f(rix[0][0]);
if (rix[0][4] > 1.75f*Y)
rix[0][0] = ULIM(rix[0][0],rix[ -1][4],rix[ 1][4]);
if (rix[4][0] > 1.75f*Y)
rix[0][0] = ULIM(rix[0][0],rix[4][ -1],rix[4][ 1]);
else
rix[0][0] = LIM(rix[0][0],0.0f,1.0f);
rix[0][0] -= rix[0][4];
rix[0][0] -= rix[4][0];
// vertical
rix[0][1] = -x0250(rix[-w2][4] + rix[w2][4])+ xdiv2f(rix[-w1][4] + rix[0][4] + rix[w1][4]);
Y = v0 + xdiv2f(rix[0][1]);
if (rix[0][4] > 1.75f*Y)
rix[0][1] = ULIM(rix[0][1],rix[-w1][4],rix[w1][4]);
rix[1] = qix[1] + rr*cc1 + cc;
rix[1][0] = -x0250(rix[4][-w2] + rix[4][w2])+ xdiv2f(rix[4][-w1] + rix[4][0] + rix[4][w1]);
Y = v0 + xdiv2f(rix[1][0]);
if (rix[4][0] > 1.75f*Y)
rix[1][0] = ULIM(rix[1][0],rix[4][-w1],rix[4][w1]);
else
rix[0][1] = LIM(rix[0][1],0.0f,1.0f);
rix[0][1] -= rix[0][4];
rix[1][0] = LIM(rix[1][0],0.0f,1.0f);
rix[1][0] -= rix[4][0];
}
// G-R(B) at G location
for (int ccc=2+(FC(rr,3)&1); ccc < cc1-2; ccc+=2) {
rix = qix + rr*cc1 + ccc;
rix[0][0] = x0250(rix[ -2][4] + rix[ 2][4])- xdiv2f(rix[ -1][4] + rix[0][4] + rix[ 1][4]);
rix[0][1] = x0250(rix[-w2][4] + rix[w2][4])- xdiv2f(rix[-w1][4] + rix[0][4] + rix[w1][4]);
rix[0][0] = LIM(rix[0][0],-1.0f,0.0f) + rix[0][4];
rix[0][1] = LIM(rix[0][1],-1.0f,0.0f) + rix[0][4];
rix[0] = qix[0] + rr*cc1 + ccc;
rix[1] = qix[1] + rr*cc1 + ccc;
rix[4] = qix[4] + rr*cc1 + ccc;
rix[0][0] = x0250(rix[4][ -2] + rix[4][ 2])- xdiv2f(rix[4][ -1] + rix[4][0] + rix[4][ 1]);
rix[1][0] = x0250(rix[4][-w2] + rix[4][w2])- xdiv2f(rix[4][-w1] + rix[4][0] + rix[4][w1]);
rix[0][0] = LIM(rix[0][0],-1.0f,0.0f) + rix[4][0];
rix[1][0] = LIM(rix[1][0],-1.0f,0.0f) + rix[4][0];
}
}
#ifdef _OPENMP
@@ -1325,9 +1337,12 @@ void RawImageSource::lmmse_interpolate_omp(int winw, int winh, int iterations)
#endif
for (int rr=4; rr < rr1-4; rr++)
for (int cc=4; cc < cc1-4; cc++) {
rix = qix + rr*cc1 + cc;
rix[0][2] = h0*rix[0][0] + h1*(rix[ -1][0] + rix[ 1][0]) + h2*(rix[ -2][0] + rix[ 2][0]) + h3*(rix[ -3][0] + rix[ 3][0]) + h4*(rix[ -4][0] + rix[ 4][0]);
rix[0][3] = h0*rix[0][1] + h1*(rix[-w1][1] + rix[w1][1]) + h2*(rix[-w2][1] + rix[w2][1]) + h3*(rix[-w3][1] + rix[w3][1]) + h4*(rix[-w4][1] + rix[w4][1]);
rix[0] = qix[0] + rr*cc1 + cc;
rix[2] = qix[2] + rr*cc1 + cc;
rix[2][0] = h0*rix[0][0] + h1*(rix[0][ -1] + rix[0][ 1]) + h2*(rix[0][ -2] + rix[0][ 2]) + h3*(rix[0][ -3] + rix[0][ 3]) + h4*(rix[0][ -4] + rix[0][ 4]);
rix[1] = qix[1] + rr*cc1 + cc;
rix[3] = qix[3] + rr*cc1 + cc;
rix[3][0] = h0*rix[1][0] + h1*(rix[1][-w1] + rix[1][w1]) + h2*(rix[1][-w2] + rix[1][w2]) + h3*(rix[1][-w3] + rix[1][w3]) + h4*(rix[1][-w4] + rix[1][w4]);
}
#ifdef _OPENMP
#pragma omp single
@@ -1343,57 +1358,61 @@ void RawImageSource::lmmse_interpolate_omp(int winw, int winh, int iterations)
#endif
for (int rr=4; rr < rr1-4; rr++)
for (int cc=4+(FC(rr,4)&1); cc < cc1-4; cc+=2) {
rix = qix + rr*cc1 + cc;
rix[0] = qix[0] + rr*cc1 + cc;
rix[1] = qix[1] + rr*cc1 + cc;
rix[2] = qix[2] + rr*cc1 + cc;
rix[3] = qix[3] + rr*cc1 + cc;
rix[4] = qix[4] + rr*cc1 + cc;
// horizontal
float mu = (rix[-4][2] + rix[-3][2] + rix[-2][2] + rix[-1][2] + rix[0][2]+ rix[ 1][2] + rix[ 2][2] + rix[ 3][2] + rix[ 4][2]) / 9.0f;
float p1 = rix[-4][2] - mu;
float p2 = rix[-3][2] - mu;
float p3 = rix[-2][2] - mu;
float p4 = rix[-1][2] - mu;
float p5 = rix[ 0][2] - mu;
float p6 = rix[ 1][2] - mu;
float p7 = rix[ 2][2] - mu;
float p8 = rix[ 3][2] - mu;
float p9 = rix[ 4][2] - mu;
float mu = (rix[2][-4] + rix[2][-3] + rix[2][-2] + rix[2][-1] + rix[2][0]+ rix[2][ 1] + rix[2][ 2] + rix[2][ 3] + rix[2][ 4]) / 9.0f;
float p1 = rix[2][-4] - mu;
float p2 = rix[2][-3] - mu;
float p3 = rix[2][-2] - mu;
float p4 = rix[2][-1] - mu;
float p5 = rix[2][ 0] - mu;
float p6 = rix[2][ 1] - mu;
float p7 = rix[2][ 2] - mu;
float p8 = rix[2][ 3] - mu;
float p9 = rix[2][ 4] - mu;
float vx = 1e-7+p1*p1+p2*p2+p3*p3+p4*p4+p5*p5+p6*p6+p7*p7+p8*p8+p9*p9;
p1 = rix[-4][0] - rix[-4][2];
p2 = rix[-3][0] - rix[-3][2];
p3 = rix[-2][0] - rix[-2][2];
p4 = rix[-1][0] - rix[-1][2];
p5 = rix[ 0][0] - rix[ 0][2];
p6 = rix[ 1][0] - rix[ 1][2];
p7 = rix[ 2][0] - rix[ 2][2];
p8 = rix[ 3][0] - rix[ 3][2];
p9 = rix[ 4][0] - rix[ 4][2];
p1 = rix[0][-4] - rix[2][-4];
p2 = rix[0][-3] - rix[2][-3];
p3 = rix[0][-2] - rix[2][-2];
p4 = rix[0][-1] - rix[2][-1];
p5 = rix[0][ 0] - rix[2][ 0];
p6 = rix[0][ 1] - rix[2][ 1];
p7 = rix[0][ 2] - rix[2][ 2];
p8 = rix[0][ 3] - rix[2][ 3];
p9 = rix[0][ 4] - rix[2][ 4];
float vn = 1e-7+p1*p1+p2*p2+p3*p3+p4*p4+p5*p5+p6*p6+p7*p7+p8*p8+p9*p9;
float xh = (rix[0][0]*vx + rix[0][2]*vn)/(vx + vn);
float xh = (rix[0][0]*vx + rix[2][0]*vn)/(vx + vn);
float vh = vx*vn/(vx + vn);
// vertical
mu = (rix[-w4][3] + rix[-w3][3] + rix[-w2][3] + rix[-w1][3] + rix[0][3]+rix[ w1][3] + rix[ w2][3] + rix[ w3][3] + rix[ w4][3]) / 9.0f;
p1 = rix[-w4][3] - mu;
p2 = rix[-w3][3] - mu;
p3 = rix[-w2][3] - mu;
p4 = rix[-w1][3] - mu;
p5 = rix[ 0][3] - mu;
p6 = rix[ w1][3] - mu;
p7 = rix[ w2][3] - mu;
p8 = rix[ w3][3] - mu;
p9 = rix[ w4][3] - mu;
mu = (rix[3][-w4] + rix[3][-w3] + rix[3][-w2] + rix[3][-w1] + rix[3][0]+rix[3][ w1] + rix[3][ w2] + rix[3][ w3] + rix[3][ w4]) / 9.0f;
p1 = rix[3][-w4] - mu;
p2 = rix[3][-w3] - mu;
p3 = rix[3][-w2] - mu;
p4 = rix[3][-w1] - mu;
p5 = rix[3][ 0] - mu;
p6 = rix[3][ w1] - mu;
p7 = rix[3][ w2] - mu;
p8 = rix[3][ w3] - mu;
p9 = rix[3][ w4] - mu;
vx = 1e-7+p1*p1+p2*p2+p3*p3+p4*p4+p5*p5+p6*p6+p7*p7+p8*p8+p9*p9;
p1 = rix[-w4][1] - rix[-w4][3];
p2 = rix[-w3][1] - rix[-w3][3];
p3 = rix[-w2][1] - rix[-w2][3];
p4 = rix[-w1][1] - rix[-w1][3];
p5 = rix[ 0][1] - rix[ 0][3];
p6 = rix[ w1][1] - rix[ w1][3];
p7 = rix[ w2][1] - rix[ w2][3];
p8 = rix[ w3][1] - rix[ w3][3];
p9 = rix[ w4][1] - rix[ w4][3];
p1 = rix[1][-w4] - rix[3][-w4];
p2 = rix[1][-w3] - rix[3][-w3];
p3 = rix[1][-w2] - rix[3][-w2];
p4 = rix[1][-w1] - rix[3][-w1];
p5 = rix[1][ 0] - rix[3][ 0];
p6 = rix[1][ w1] - rix[3][ w1];
p7 = rix[1][ w2] - rix[3][ w2];
p8 = rix[1][ w3] - rix[3][ w3];
p9 = rix[1][ w4] - rix[3][ w4];
vn = 1e-7+p1*p1+p2*p2+p3*p3+p4*p4+p5*p5+p6*p6+p7*p7+p8*p8+p9*p9;
float xv = (rix[0][1]*vx + rix[0][3]*vn)/(vx + vn);
float xv = (rix[1][0]*vx + rix[3][0]*vn)/(vx + vn);
float vv = vx*vn/(vx + vn);
// interpolated G-R(B)
rix[0][4] = (xh*vv + xv*vh)/(vh + vv);
rix[4][0] = (xh*vv + xv*vh)/(vh + vv);
}
#ifdef _OPENMP
#pragma omp single
@@ -1410,17 +1429,21 @@ void RawImageSource::lmmse_interpolate_omp(int winw, int winh, int iterations)
for (int rr=0; rr < rr1; rr++)
for (int cc=0, row=rr-ba; cc < cc1; cc++) {
int col=cc-ba;
rix = qix + rr*cc1 + cc;
int c = FC(rr,cc);
rix[c] = qix[c] + rr*cc1 + cc;
if ((row >= 0) & (row < height) & (col >= 0) & (col < width)) {
if (applyGamma)
rix[0][c] = Color::gammatab_24_17a[image[row*width+col][c]];
rix[c][0] = Color::gammatab_24_17a[image[c][(((row*width)>>(1-(c&1)))+col)>>1]];
else
rix[0][c] = (float)image[row*width+col][c]/65535.0f;
rix[c][0] = (float)image[c][(((row*width)>>(1-(c&1)))+col)>>1]/65535.0f;
}
else
rix[0][c] = 0.f;
if (c != 1) rix[0][1] = rix[0][c] + rix[0][4];
rix[c][0] = 0.f;
if (c != 1) {
rix[1] = qix[1] + rr*cc1 + cc;
rix[4] = qix[4] + rr*cc1 + cc;
rix[1][0] = rix[c][0] + rix[4][0];
}
}
#ifdef _OPENMP
#pragma omp single
@@ -1437,10 +1460,12 @@ void RawImageSource::lmmse_interpolate_omp(int winw, int winh, int iterations)
#endif
for (int rr=1; rr < rr1-1; rr++)
for (int cc=1+(FC(rr,2)&1), c=FC(rr,cc+1); cc < cc1-1; cc+=2) {
rix = qix + rr*cc1 + cc;
rix[0][c] = rix[0][1] + xdiv2f(rix[ -1][c] - rix[ -1][1] + rix[ 1][c] - rix[ 1][1]);
rix[c] = qix[c] + rr*cc1 + cc;
rix[1] = qix[1] + rr*cc1 + cc;
rix[c][0] = rix[1][0] + xdiv2f(rix[c][ -1] - rix[1][ -1] + rix[c][ 1] - rix[1][ 1]);
c = 2 - c;
rix[0][c] = rix[0][1]+ xdiv2f(rix[-w1][c] - rix[-w1][1] + rix[w1][c] - rix[w1][1]);
rix[c] = qix[c] + rr*cc1 + cc;
rix[c][0] = rix[1][0]+ xdiv2f(rix[c][-w1] - rix[1][-w1] + rix[c][w1] - rix[1][w1]);
c = 2 - c;
}
#ifdef _OPENMP
@@ -1457,8 +1482,9 @@ void RawImageSource::lmmse_interpolate_omp(int winw, int winh, int iterations)
#endif
for (int rr=1; rr < rr1-1; rr++)
for (int cc=1+(FC(rr,1)&1), c=2-FC(rr,cc); cc < cc1-1; cc+=2) {
rix = qix + rr*cc1 + cc;
rix[0][c] = rix[0][1]+ x0250(rix[-w1][c] - rix[-w1][1] + rix[ -1][c] - rix[ -1][1]+ rix[ 1][c] - rix[ 1][1] + rix[ w1][c] - rix[ w1][1]);
rix[c] = qix[c] + rr*cc1 + cc;
rix[1] = qix[1] + rr*cc1 + cc;
rix[c][0] = rix[1][0]+ x0250(rix[c][-w1] - rix[1][-w1] + rix[c][ -1] - rix[1][ -1]+ rix[c][ 1] - rix[1][ 1] + rix[c][ w1] - rix[1][ w1]);
}
#ifdef _OPENMP
#pragma omp single
@@ -1474,7 +1500,12 @@ void RawImageSource::lmmse_interpolate_omp(int winw, int winh, int iterations)
for (int c=0; c < 3; c+=2) {
// Compute median(R-G) and median(B-G)
int d = c + 3;
for (int ii=0; ii < rr1*cc1; ii++) qix[ii][d] = qix[ii][c] - qix[ii][1];
#ifdef _OPENMP
#pragma omp parallel for firstprivate (qix)
#endif
for (int ii=0; ii < rr1*cc1; ii++)
qix[4][ii] = qix[c][ii] - qix[1][ii];
// Apply 3x3 median filter
#ifdef _OPENMP
#pragma omp parallel for firstprivate (rix,qix)
@@ -1482,11 +1513,12 @@ void RawImageSource::lmmse_interpolate_omp(int winw, int winh, int iterations)
for (int rr=1; rr < rr1-1; rr++)
for (int cc=1; cc < cc1-1; cc++) {
float temp;
rix = qix + rr*cc1 + cc;
rix[d] = qix[d] + rr*cc1 + cc;
rix[4] = qix[4] + rr*cc1 + cc;
// Assign 3x3 differential color values
float p1 = rix[-w1-1][d]; float p2 = rix[-w1][d]; float p3 = rix[-w1+1][d];
float p4 = rix[ -1][d]; float p5 = rix[ 0][d]; float p6 = rix[ 1][d];
float p7 = rix[ w1-1][d]; float p8 = rix[ w1][d]; float p9 = rix[ w1+1][d];
float p1 = rix[4][-w1-1]; float p2 = rix[4][-w1]; float p3 = rix[4][-w1+1];
float p4 = rix[4][ -1]; float p5 = rix[4][ 0]; float p6 = rix[4][ 1];
float p7 = rix[4][ w1-1]; float p8 = rix[4][ w1]; float p9 = rix[4][ w1+1];
// Sort for median of 9 values
PIX_SORT(p2,p3); PIX_SORT(p5,p6); PIX_SORT(p8,p9);
PIX_SORT(p1,p2); PIX_SORT(p4,p5); PIX_SORT(p7,p8);
@@ -1495,23 +1527,36 @@ void RawImageSource::lmmse_interpolate_omp(int winw, int winh, int iterations)
PIX_SORT(p4,p7); PIX_SORT(p2,p5); PIX_SORT(p3,p6);
PIX_SORT(p5,p8); PIX_SORT(p5,p3); PIX_SORT(p7,p5);
PIX_SORT(p5,p3);
rix[0][4] = p5;
rix[d][0] = p5;
}
for (int ii=0; ii < rr1*cc1; ii++) qix[ii][d] = qix[ii][4];
}
// red/blue at GREEN pixel locations
#ifdef _OPENMP
#pragma omp parallel for firstprivate (rix,qix)
#endif
for (int rr=0; rr < rr1; rr++)
for (int cc=(FC(rr,1)&1) /*, c=FC(rr,cc+1)*/; cc < cc1; cc+=2) {
rix = qix + rr*cc1 + cc;
rix[0][0] = rix[0][1] + rix[0][3];
rix[0][2] = rix[0][1] + rix[0][5];
rix[0] = qix[0] + rr*cc1 + cc;
rix[1] = qix[1] + rr*cc1 + cc;
rix[2] = qix[2] + rr*cc1 + cc;
rix[3] = qix[3] + rr*cc1 + cc;
rix[5] = qix[5] + rr*cc1 + cc;
rix[0][0] = rix[1][0] + rix[3][0];
rix[2][0] = rix[1][0] + rix[5][0];
}
// red/blue and green at BLUE/RED pixel locations
#ifdef _OPENMP
#pragma omp parallel for firstprivate (rix,qix)
#endif
for (int rr=0; rr < rr1; rr++)
for (int cc=(FC(rr,0)&1), c=2-FC(rr,cc),d=c+3; cc < cc1; cc+=2) {
rix = qix + rr*cc1 + cc;
rix[0][c] = rix[0][1] + rix[0][d];
rix[0][1] = xdiv2f(rix[0][0] - rix[0][3] + rix[0][2] - rix[0][5]);
rix[0] = qix[0] + rr*cc1 + cc;
rix[1] = qix[1] + rr*cc1 + cc;
rix[2] = qix[2] + rr*cc1 + cc;
rix[3] = qix[3] + rr*cc1 + cc;
rix[5] = qix[5] + rr*cc1 + cc;
rix[c][0] = rix[1][0] + rix[d][0];
rix[1][0] = 0.5f * (rix[0][0] - rix[3][0] + rix[2][0] - rix[5][0]);
}
}
@@ -1524,41 +1569,63 @@ void RawImageSource::lmmse_interpolate_omp(int winw, int winh, int iterations)
#ifdef _OPENMP
#pragma omp for
#endif
for (int row=0; row < height; row++)
for (int row=0; row < height; row++) {
for (int col=0, rr=row+ba; col < width; col++) {
int cc = col+ba;
rix = qix + rr*cc1 + cc;
int c = FC(row,col);
float v0;
array2D<float> *nonInterpolated;
if(c==0)
nonInterpolated = &red;
else if(c==1)
nonInterpolated = &green;
else
nonInterpolated = &blue;
if (applyGamma) {
for (int ii=0; ii < 3; ii++)
if (ii != c) {
v0 = 65535.f*rix[0][ii];
image[row*width+col][ii]=Color::igammatab_24_17[v0];
array2D<float> *interpolated;
if(ii==0)
interpolated = &red;
else if(ii==1)
interpolated = &green;
else
interpolated = &blue;
rix[ii] = qix[ii] + rr*cc1 + cc;
float v0 = 65535.f*rix[ii][0];
(*interpolated)[row][col] = Color::igammatab_24_17[v0];
} else {
(*nonInterpolated)[row][col] = image[ii][(((row*width)>>(1-(ii&1)))+col)>>1];
}
}
else
for (int ii=0; ii < 3; ii++)
if (ii != c)
image[row*width+col][ii] = ((65535.0f*rix[0][ii] + 0.5f));
}
#ifdef _OPENMP
#pragma omp for
#endif
for (int ii=0; ii<height; ii++) {
for (int jj=0; jj<width; jj++){
red[ii][jj] = (image[ii*width+jj][0]);
green[ii][jj] = (image[ii*width+jj][1]);
blue[ii][jj] = (image[ii*width+jj][2]);
if (ii != c) {
array2D<float> *interpolated;
if(ii==0)
interpolated = &red;
else if(ii==1)
interpolated = &green;
else
interpolated = &blue;
rix[ii] = qix[ii] + rr*cc1 + cc;
(*interpolated)[row][col] = ((65535.0f*rix[ii][0] + 0.5f));
} else {
(*nonInterpolated)[row][col] = image[ii][(((row*width)>>(1-(ii&1)))+col)>>1];
}
}
}
}
// End of parallelization 2
if (plistener) plistener->setProgress (1.0);
free(buffer);
free(image);
if(buffer)
free(buffer);
else
for(int i=0;i<6;i++)
free(qix[i]);
free(imageBuffer);
//if(iterations > 4) refinement_lassus(passref);
if(iterations > 4 && iterations <=6) refinement(passref);
else if(iterations > 6) refinement_lassus(passref);
@@ -2388,101 +2455,88 @@ void RawImageSource::refinement(int PassCount)
{
MyTime t1e,t2e;
t1e.set();
float (*image)[3];
int width=W;
int height=H;
int w1 = width;
int w2 = 2*w1;
if (plistener) {
plistener->setProgressStr (M("TP_RAW_DMETHOD_PROGRESSBAR_REFINE"));
}
array2D<float> *rgb[3];
rgb[0] = &red;
rgb[1] = &green;
rgb[2] = &blue;
for (int b=0; b<PassCount; b++) {
if (plistener) {
plistener->setProgress ((float)b/PassCount);
}
image = (float(*)[3]) calloc(W*H, sizeof *image);
#ifdef _OPENMP
#pragma omp parallel shared(image)
#pragma omp parallel
#endif
{
{
float *pix[3];
/* Reinforce interpolated green pixels on RED/BLUE pixel locations */
#ifdef _OPENMP
#pragma omp for
#endif
for (int i=0;i<H;i++) {
for (int j=0;j<W;j++) {
image[i*W+j][0] = red [i][j];
image[i*W+j][1] = green[i][j];
image[i*W+j][2] = blue [i][j];
}
}
for (int b=0; b<PassCount; b++) {
if (plistener) {
plistener->setProgressStr (M("TP_RAW_DMETHOD_PROGRESSBAR_REFINE"));
plistener->setProgress ((float)b/PassCount);
}
/* Reinforce interpolated green pixels on RED/BLUE pixel locations */
#ifdef _OPENMP
#pragma omp for
#endif
for (int row=2; row < height-2; row++)
for (int col=2+(FC(row,2) & 1), c=FC(row,col); col < width-2; col+=2) {
int indx = row*width+col;
float (*pix)[3]= image + indx;
double dL = 1.0/(1.0+fabs(pix[ -2][c]-pix[0][c])+fabs(pix[ 1][1]-pix[ -1][1]));
double dR = 1.0/(1.0+fabs(pix[ 2][c]-pix[0][c])+fabs(pix[ 1][1]-pix[ -1][1]));
double dU = 1.0/(1.0+fabs(pix[-w2][c]-pix[0][c])+fabs(pix[w1][1]-pix[-w1][1]));
double dD = 1.0/(1.0+fabs(pix[ w2][c]-pix[0][c])+fabs(pix[w1][1]-pix[-w1][1]));
float v0 = (float)((double)pix[0][c] + 0.5 +((double)(pix[ -1][1]-pix[ -1][c])*dL +(double)(pix[ 1][1]-pix[ 1][c])*dR +(double)(pix[-w1][1]-pix[-w1][c])*dU +
(double)(pix[ w1][1]-pix[ w1][c])*dD ) / (dL+dR+dU+dD));
pix[0][1] = CLIP(v0);
}
for (int row=2; row < height-2; row++)
for (int col=2+(FC(row,2) & 1), c=FC(row,col); col < width-2; col+=2) {
int indx = row*width+col;
pix[c] = (float*)(*rgb[c]) + indx;
pix[1] = (float*)(*rgb[1]) + indx;
float dL = 1.f/(1.f+fabsf(pix[c][ -2]-pix[c][0])+fabsf(pix[1][ 1]-pix[1][ -1]));
float dR = 1.f/(1.f+fabsf(pix[c][ 2]-pix[c][0])+fabsf(pix[1][ 1]-pix[1][ -1]));
float dU = 1.f/(1.f+fabsf(pix[c][-w2]-pix[c][0])+fabsf(pix[1][w1]-pix[1][-w1]));
float dD = 1.f/(1.f+fabsf(pix[c][ w2]-pix[c][0])+fabsf(pix[1][w1]-pix[1][-w1]));
float v0 = (pix[c][0] + 0.5f +((pix[1][ -1]-pix[c][ -1])*dL +(pix[1][ 1]-pix[c][ 1])*dR +(pix[1][-w1]-pix[c][-w1])*dU + (pix[1][ w1]-pix[c][ w1])*dD ) / (dL+dR+dU+dD));
pix[1][0] = CLIP(v0);
}
/* Reinforce interpolated red/blue pixels on GREEN pixel locations */
#ifdef _OPENMP
#pragma omp for
#endif
for (int row=2; row < height-2; row++)
for (int col=2+(FC(row,3) & 1), c=FC(row,col+1); col < width-2; col+=2) {
int indx = row*width+col;
float (*pix)[3] = image + indx;
for (int i=0; i < 2; c=2-c, i++) {
double dL = 1.0/(1.0+fabs(pix[ -2][1]-pix[0][1])+fabs(pix[ 1][c]-pix[ -1][c]));
double dR = 1.0/(1.0+fabs(pix[ 2][1]-pix[0][1])+fabs(pix[ 1][c]-pix[ -1][c]));
double dU = 1.0/(1.0+fabs(pix[-w2][1]-pix[0][1])+fabs(pix[w1][c]-pix[-w1][c]));
double dD = 1.0/(1.0+fabs(pix[ w2][1]-pix[0][1])+fabs(pix[w1][c]-pix[-w1][c]));
float v0 = (float)((double)pix[0][1] + 0.5 -((double)(pix[ -1][1]-pix[ -1][c])*dL + (double)(pix[ 1][1]-pix[ 1][c])*dR +(double)(pix[-w1][1]-pix[-w1][c])*dU +
(double)(pix[ w1][1]-pix[ w1][c])*dD ) / (dL+dR+dU+dD));
pix[0][c] = CLIP(v0);
}
}
for (int row=2; row < height-2; row++)
for (int col=2+(FC(row,3) & 1), c=FC(row,col+1); col < width-2; col+=2) {
int indx = row*width+col;
pix[1] = (float*)(*rgb[1]) + indx;
for (int i=0; i < 2; c=2-c, i++) {
pix[c] = (float*)(*rgb[c]) + indx;
float dL = 1.f/(1.f+fabsf(pix[1][ -2]-pix[1][0])+fabsf(pix[c][ 1]-pix[c][ -1]));
float dR = 1.f/(1.f+fabsf(pix[1][ 2]-pix[1][0])+fabsf(pix[c][ 1]-pix[c][ -1]));
float dU = 1.f/(1.f+fabsf(pix[1][-w2]-pix[1][0])+fabsf(pix[c][w1]-pix[c][-w1]));
float dD = 1.f/(1.f+fabsf(pix[1][ w2]-pix[1][0])+fabsf(pix[c][w1]-pix[c][-w1]));
float v0 = (pix[1][0] + 0.5f -((pix[1][ -1]-pix[c][ -1])*dL + (pix[1][ 1]-pix[c][ 1])*dR +(pix[1][-w1]-pix[c][-w1])*dU + (pix[1][ w1]-pix[c][ w1])*dD ) / (dL+dR+dU+dD));
pix[c][0] = CLIP(v0);
}
}
/* Reinforce integrated red/blue pixels on BLUE/RED pixel locations */
/* Reinforce integrated red/blue pixels on BLUE/RED pixel locations */
#ifdef _OPENMP
#pragma omp for
#endif
for (int row=2; row < height-2; row++)
for (int col=2+(FC(row,2) & 1), c=2-FC(row,col); col < width-2; col+=2) {
int indx = row*width+col;
float (*pix)[3] = image + indx;
int d = 2 - c;
double dL = 1.0/(1.0+ABS(pix[ -2][d]-pix[0][d])+ABS(pix[ 1][1]-pix[ -1][1]));
double dR = 1.0/(1.0+ABS(pix[ 2][d]-pix[0][d])+ABS(pix[ 1][1]-pix[ -1][1]));
double dU = 1.0/(1.0+ABS(pix[-w2][d]-pix[0][d])+ABS(pix[w1][1]-pix[-w1][1]));
double dD = 1.0/(1.0+ABS(pix[ w2][d]-pix[0][d])+ABS(pix[w1][1]-pix[-w1][1]));
float v0 = (float)((double)pix[0][1] + 0.5 -((double)(pix[ -1][1]-pix[ -1][c])*dL +(double)(pix[ 1][1]-pix[ 1][c])*dR +(double)(pix[-w1][1]-pix[-w1][c])*dU +
(double)(pix[ w1][1]-pix[ w1][c])*dD ) / (dL+dR+dU+dD));
pix[0][c] = CLIP(v0);
}
}
#ifdef _OPENMP
#pragma omp for
#endif
for (int i=0;i<H;i++) {
for (int j=0; j<W; j++) {
red [i][j] =image[i*W+j][0];
green[i][j] =image[i*W+j][1];
blue [i][j] =image[i*W+j][2];
}
}
for (int row=2; row < height-2; row++)
for (int col=2+(FC(row,2) & 1), c=2-FC(row,col); col < width-2; col+=2) {
int indx = row*width+col;
pix[0] = (float*)(*rgb[0]) + indx;
pix[1] = (float*)(*rgb[1]) + indx;
pix[2] = (float*)(*rgb[2]) + indx;
int d = 2 - c;
float dL = 1.f/(1.f+fabsf(pix[d][ -2]-pix[d][0])+fabsf(pix[1][ 1]-pix[1][ -1]));
float dR = 1.f/(1.f+fabsf(pix[d][ 2]-pix[d][0])+fabsf(pix[1][ 1]-pix[1][ -1]));
float dU = 1.f/(1.f+fabsf(pix[d][-w2]-pix[d][0])+fabsf(pix[1][w1]-pix[1][-w1]));
float dD = 1.f/(1.f+fabsf(pix[d][ w2]-pix[d][0])+fabsf(pix[1][w1]-pix[1][-w1]));
float v0 = (pix[1][0] + 0.5f -((pix[1][ -1]-pix[c][ -1])*dL +(pix[1][ 1]-pix[c][ 1])*dR +(pix[1][-w1]-pix[c][-w1])*dU + (pix[1][ w1]-pix[c][ w1])*dD ) / (dL+dR+dU+dD));
pix[c][0] = CLIP(v0);
}
} // end parallel
}
free(image);
t2e.set();
if (settings->verbose) printf("Refinement Lee %d usec\n", t2e.etime(t1e));
}