Additional optimizations for Phase One P45+ and P65+

This commit is contained in:
heckflosse 2017-08-30 18:55:30 +02:00
parent 61a69ba193
commit fd1909e4da

View File

@ -1501,204 +1501,215 @@ void CLASS phase_one_flat_field (int is_float, int nc)
void CLASS phase_one_correct() void CLASS phase_one_correct()
{ {
BENCHFUN BENCHFUN
unsigned entries, tag, data, save, col, row, type; unsigned entries, tag, data, save, col, row, type;
int len, i, j, k, cip, val[4], dev[4], sum, max; int len, i, j, k, cip, val[4], dev[4], sum, max;
int head[9], diff, mindiff=INT_MAX, off_412=0; int head[9], diff, mindiff=INT_MAX, off_412=0;
static const signed char dir[12][2] = static const signed char dir[12][2] = { {-1,-1}, {-1,1}, {1,-1}, {1,1}, {-2,0}, {0,-2}, {0,2}, {2,0}, {-2,-2}, {-2,2}, {2,-2}, {2,2} };
{ {-1,-1}, {-1,1}, {1,-1}, {1,1}, {-2,0}, {0,-2}, {0,2}, {2,0}, float poly[8], num, cfrac, frac, mult[2], *yval[2];
{-2,-2}, {-2,2}, {2,-2}, {2,2} }; ushort *xval[2];
float poly[8], num, cfrac, frac, mult[2], *yval[2]; int qmult_applied = 0, qlin_applied = 0;
ushort *xval[2];
int qmult_applied = 0, qlin_applied = 0;
if (half_size || !meta_length) return; if (half_size || !meta_length) {
if (verbose) fprintf (stderr,_("Phase One correction...\n")); return;
fseek (ifp, meta_offset, SEEK_SET); }
order = get2(); if (verbose) {
fseek (ifp, 6, SEEK_CUR); fprintf (stderr,_("Phase One correction...\n"));
fseek (ifp, meta_offset+get4(), SEEK_SET); }
entries = get4(); get4(); fseek (ifp, meta_offset, SEEK_SET);
while (entries--) { order = get2();
tag = get4(); fseek (ifp, 6, SEEK_CUR);
len = get4(); fseek (ifp, meta_offset+get4(), SEEK_SET);
data = get4(); entries = get4(); get4();
save = ftell(ifp); while (entries--) {
fseek (ifp, meta_offset+data, SEEK_SET); tag = get4();
if (tag == 0x419) { /* Polynomial curve */ len = get4();
for (get4(), i=0; i < 8; i++) data = get4();
poly[i] = getreal(11); save = ftell(ifp);
poly[3] += (ph1.tag_210 - poly[7]) * poly[6] + 1; fseek (ifp, meta_offset+data, SEEK_SET);
for (i=0; i < 0x10000; i++) { if (tag == 0x419) { /* Polynomial curve */
num = (poly[5]*i + poly[3])*i + poly[1]; for (get4(), i=0; i < 8; i++) {
curve[i] = LIM(num,0,65535); poly[i] = getreal(11);
} goto apply; /* apply to right half */ }
} else if (tag == 0x41a) { /* Polynomial curve */ poly[3] += (ph1.tag_210 - poly[7]) * poly[6] + 1;
for (i=0; i < 4; i++) for (i=0; i < 0x10000; i++) {
poly[i] = getreal(11); num = (poly[5]*i + poly[3])*i + poly[1];
for (i=0; i < 0x10000; i++) { curve[i] = LIM(num,0,65535);
for (num=0, j=4; j--; ) }
num = num * i + poly[j]; goto apply; /* apply to right half */
curve[i] = LIM(num+i,0,65535); } else if (tag == 0x41a) { /* Polynomial curve */
} apply: /* apply to whole image */ for (i=0; i < 4; i++) {
for (row=0; row < raw_height; row++) poly[i] = getreal(11);
for (col = (tag & 1)*ph1.split_col; col < raw_width; col++) }
RAW(row,col) = curve[RAW(row,col)]; for (i=0; i < 0x10000; i++) {
} else if (tag == 0x400) { /* Sensor defects */ for (num=0, j=4; j--;) {
while ((len -= 8) >= 0) { num = num * i + poly[j];
col = get2(); }
row = get2(); curve[i] = LIM(num+i,0,65535);
type = get2(); get2(); }
if (col >= raw_width) continue; apply: /* apply to whole image */
if (type == 131 || type == 137) /* Bad column */ #pragma omp parallel for schedule(dynamic,16)
for (row=0; row < raw_height; row++) for (int row=0; row < raw_height; row++) {
if (FC(row-top_margin,col-left_margin) == 1) { for (int col = (tag & 1)*ph1.split_col; col < raw_width; col++) {
for (sum=i=0; i < 4; i++) RAW(row,col) = curve[RAW(row,col)];
sum += val[i] = raw (row+dir[i][0], col+dir[i][1]); }
for (max=i=0; i < 4; i++) { }
dev[i] = abs((val[i] << 2) - sum); } else if (tag == 0x400) { /* Sensor defects */
if (dev[max] < dev[i]) max = i; while ((len -= 8) >= 0) {
} col = get2();
RAW(row,col) = (sum - val[max])/3.0 + 0.5; row = get2();
} else { type = get2();
for (sum=0, i=8; i < 12; i++) get2();
sum += raw (row+dir[i][0], col+dir[i][1]); if (col >= raw_width) continue;
RAW(row,col) = 0.5 + sum * 0.0732233 + if (type == 131 || type == 137) { /* Bad column */
(raw(row,col-2) + raw(row,col+2)) * 0.3535534; for (row=0; row < raw_height; row++) {
} if (FC(row-top_margin,col-left_margin) == 1) {
else if (type == 129) { /* Bad pixel */ for (sum=i=0; i < 4; i++)
if (row >= raw_height) continue; sum += val[i] = raw (row+dir[i][0], col+dir[i][1]);
j = (FC(row-top_margin,col-left_margin) != 1) * 4; for (max=i=0; i < 4; i++) {
for (sum=0, i=j; i < j+8; i++) dev[i] = abs((val[i] << 2) - sum);
sum += raw (row+dir[i][0], col+dir[i][1]); if (dev[max] < dev[i]) max = i;
RAW(row,col) = (sum + 4) >> 3; }
} RAW(row,col) = (sum - val[max])/3.0 + 0.5;
} } else {
} else if (tag == 0x401) { /* All-color flat fields */ for (sum=0, i=8; i < 12; i++)
phase_one_flat_field (1, 2); sum += raw (row+dir[i][0], col+dir[i][1]);
} else if (tag == 0x416 || tag == 0x410) { RAW(row,col) = 0.5 + sum * 0.0732233 + (raw(row,col-2) + raw(row,col+2)) * 0.3535534;
phase_one_flat_field (0, 2); }
} else if (tag == 0x40b) { /* Red+blue flat field */ }
phase_one_flat_field (0, 4); } else if (type == 129) { /* Bad pixel */
} else if (tag == 0x412) { if (row >= raw_height) continue;
fseek (ifp, 36, SEEK_CUR); j = (FC(row-top_margin,col-left_margin) != 1) * 4;
diff = abs (get2() - ph1.tag_21a); for (sum=0, i=j; i < j+8; i++)
if (mindiff > diff) { sum += raw (row+dir[i][0], col+dir[i][1]);
mindiff = diff; RAW(row,col) = (sum + 4) >> 3;
off_412 = ftell(ifp) - 38; }
} }
} else if (tag == 0x41f && !qlin_applied) { /* Quadrant linearization */ } else if (tag == 0x401) { /* All-color flat fields */
ushort lc[2][2][16], ref[16]; phase_one_flat_field (1, 2);
int qr, qc; } else if (tag == 0x416 || tag == 0x410) {
for (qr = 0; qr < 2; qr++) phase_one_flat_field (0, 2);
for (qc = 0; qc < 2; qc++) } else if (tag == 0x40b) { /* Red+blue flat field */
for (i = 0; i < 16; i++) phase_one_flat_field (0, 4);
lc[qr][qc][i] = get4(); } else if (tag == 0x412) {
for (i = 0; i < 16; i++) { fseek (ifp, 36, SEEK_CUR);
int v = 0; diff = abs (get2() - ph1.tag_21a);
for (qr = 0; qr < 2; qr++) if (mindiff > diff) {
for (qc = 0; qc < 2; qc++) mindiff = diff;
v += lc[qr][qc][i]; off_412 = ftell(ifp) - 38;
ref[i] = (v + 2) >> 2; }
} } else if (tag == 0x41f && !qlin_applied) { /* Quadrant linearization */
for (qr = 0; qr < 2; qr++) { ushort lc[2][2][16], ref[16];
for (qc = 0; qc < 2; qc++) { int qr, qc;
int cx[19], cf[19]; for (qr = 0; qr < 2; qr++)
for (i = 0; i < 16; i++) { for (qc = 0; qc < 2; qc++)
cx[1+i] = lc[qr][qc][i]; for (i = 0; i < 16; i++)
cf[1+i] = ref[i]; lc[qr][qc][i] = get4();
} for (i = 0; i < 16; i++) {
cx[0] = cf[0] = 0; int v = 0;
cx[17] = cf[17] = ((unsigned) ref[15] * 65535) / lc[qr][qc][15]; for (qr = 0; qr < 2; qr++)
cx[18] = cf[18] = 65535; for (qc = 0; qc < 2; qc++)
cubic_spline(cx, cf, 19); v += lc[qr][qc][i];
for (row = (qr ? ph1.split_row : 0); ref[i] = (v + 2) >> 2;
row < (qr ? raw_height : ph1.split_row); row++) }
for (col = (qc ? ph1.split_col : 0); for (qr = 0; qr < 2; qr++) {
col < (qc ? raw_width : ph1.split_col); col++) for (qc = 0; qc < 2; qc++) {
RAW(row,col) = curve[RAW(row,col)]; int cx[19], cf[19];
} for (i = 0; i < 16; i++) {
} cx[1+i] = lc[qr][qc][i];
qlin_applied = 1; cf[1+i] = ref[i];
} else if (tag == 0x41e && !qmult_applied) { /* Quadrant multipliers */ }
float qmult[2][2] = { { 1, 1 }, { 1, 1 } }; cx[0] = cf[0] = 0;
get4(); get4(); get4(); get4(); cx[17] = cf[17] = ((unsigned) ref[15] * 65535) / lc[qr][qc][15];
qmult[0][0] = 1.0 + getreal(11); cx[18] = cf[18] = 65535;
get4(); get4(); get4(); get4(); get4(); cubic_spline(cx, cf, 19);
qmult[0][1] = 1.0 + getreal(11); #pragma omp parallel for schedule(dynamic,16)
get4(); get4(); get4(); for (int row = (qr ? ph1.split_row : 0); row < (qr ? raw_height : ph1.split_row); row++)
qmult[1][0] = 1.0 + getreal(11); for (int col = (qc ? ph1.split_col : 0); col < (qc ? raw_width : ph1.split_col); col++)
get4(); get4(); get4(); RAW(row,col) = curve[RAW(row,col)];
qmult[1][1] = 1.0 + getreal(11); }
for (row=0; row < raw_height; row++) }
for (col=0; col < raw_width; col++) { qlin_applied = 1;
i = qmult[row >= ph1.split_row][col >= ph1.split_col] * RAW(row,col); } else if (tag == 0x41e && !qmult_applied) { /* Quadrant multipliers */
RAW(row,col) = LIM(i,0,65535); float qmult[2][2] = { { 1, 1 }, { 1, 1 } };
} get4(); get4(); get4(); get4();
qmult_applied = 1; qmult[0][0] = 1.0 + getreal(11);
} else if (tag == 0x431 && !qmult_applied) { /* Quadrant combined */ get4(); get4(); get4(); get4(); get4();
ushort lc[2][2][7], ref[7]; qmult[0][1] = 1.0 + getreal(11);
int qr, qc; get4(); get4(); get4();
for (i = 0; i < 7; i++) qmult[1][0] = 1.0 + getreal(11);
ref[i] = get4(); get4(); get4(); get4();
for (qr = 0; qr < 2; qr++) qmult[1][1] = 1.0 + getreal(11);
for (qc = 0; qc < 2; qc++) #pragma omp parallel for schedule(dynamic,16)
for (i = 0; i < 7; i++) for (int row=0; row < raw_height; row++) {
lc[qr][qc][i] = get4(); for (int col=0; col < raw_width; col++) {
for (qr = 0; qr < 2; qr++) { int i = qmult[row >= ph1.split_row][col >= ph1.split_col] * RAW(row,col);
for (qc = 0; qc < 2; qc++) { RAW(row,col) = LIM(i,0,65535);
int cx[9], cf[9]; }
for (i = 0; i < 7; i++) { }
cx[1+i] = ref[i]; qmult_applied = 1;
cf[1+i] = ((unsigned) ref[i] * lc[qr][qc][i]) / 10000; } else if (tag == 0x431 && !qmult_applied) { /* Quadrant combined */
} ushort lc[2][2][7], ref[7];
cx[0] = cf[0] = 0; int qr, qc;
cx[8] = cf[8] = 65535; for (i = 0; i < 7; i++)
cubic_spline(cx, cf, 9); ref[i] = get4();
for (row = (qr ? ph1.split_row : 0); for (qr = 0; qr < 2; qr++)
row < (qr ? raw_height : ph1.split_row); row++) for (qc = 0; qc < 2; qc++)
for (col = (qc ? ph1.split_col : 0); for (i = 0; i < 7; i++)
col < (qc ? raw_width : ph1.split_col); col++) lc[qr][qc][i] = get4();
RAW(row,col) = curve[RAW(row,col)]; for (qr = 0; qr < 2; qr++) {
} for (qc = 0; qc < 2; qc++) {
} int cx[9], cf[9];
qmult_applied = 1; for (i = 0; i < 7; i++) {
qlin_applied = 1; cx[1+i] = ref[i];
cf[1+i] = ((unsigned) ref[i] * lc[qr][qc][i]) / 10000;
}
cx[0] = cf[0] = 0;
cx[8] = cf[8] = 65535;
cubic_spline(cx, cf, 9);
for (row = (qr ? ph1.split_row : 0); row < (qr ? raw_height : ph1.split_row); row++)
for (col = (qc ? ph1.split_col : 0); col < (qc ? raw_width : ph1.split_col); col++)
RAW(row,col) = curve[RAW(row,col)];
}
}
qmult_applied = 1;
qlin_applied = 1;
}
fseek (ifp, save, SEEK_SET);
}
if (off_412) {
fseek (ifp, off_412, SEEK_SET);
for (i=0; i < 9; i++)
head[i] = get4() & 0x7fff;
yval[0] = (float *) calloc (head[1]*head[3] + head[2]*head[4], 6);
merror (yval[0], "phase_one_correct()");
yval[1] = (float *) (yval[0] + head[1]*head[3]);
xval[0] = (ushort *) (yval[1] + head[2]*head[4]);
xval[1] = (ushort *) (xval[0] + head[1]*head[3]);
get2();
for (i=0; i < 2; i++)
for (j=0; j < head[i+1]*head[i+3]; j++)
yval[i][j] = getreal(11);
for (i=0; i < 2; i++)
for (j=0; j < head[i+1]*head[i+3]; j++)
xval[i][j] = get2();
for (row=0; row < raw_height; row++)
for (col=0; col < raw_width; col++) {
cfrac = (float) col * head[3] / raw_width;
cfrac -= cip = cfrac;
num = RAW(row,col) * 0.5;
for (i=cip; i < cip+2; i++) {
for (k=j=0; j < head[1]; j++)
if (num < xval[0][k = head[1]*i+j])
break;
frac = (j == 0 || j == head[1]) ? 0 : (xval[0][k] - num) / (xval[0][k] - xval[0][k-1]);
mult[i-cip] = yval[0][k-1] * frac + yval[0][k] * (1-frac);
}
i = ((mult[0] * (1-cfrac) + mult[1] * cfrac) * row + num) * 2;
RAW(row,col) = LIM(i,0,65535);
}
free (yval[0]);
} }
fseek (ifp, save, SEEK_SET);
}
if (off_412) {
fseek (ifp, off_412, SEEK_SET);
for (i=0; i < 9; i++) head[i] = get4() & 0x7fff;
yval[0] = (float *) calloc (head[1]*head[3] + head[2]*head[4], 6);
merror (yval[0], "phase_one_correct()");
yval[1] = (float *) (yval[0] + head[1]*head[3]);
xval[0] = (ushort *) (yval[1] + head[2]*head[4]);
xval[1] = (ushort *) (xval[0] + head[1]*head[3]);
get2();
for (i=0; i < 2; i++)
for (j=0; j < head[i+1]*head[i+3]; j++)
yval[i][j] = getreal(11);
for (i=0; i < 2; i++)
for (j=0; j < head[i+1]*head[i+3]; j++)
xval[i][j] = get2();
for (row=0; row < raw_height; row++)
for (col=0; col < raw_width; col++) {
cfrac = (float) col * head[3] / raw_width;
cfrac -= cip = cfrac;
num = RAW(row,col) * 0.5;
for (i=cip; i < cip+2; i++) {
for (k=j=0; j < head[1]; j++)
if (num < xval[0][k = head[1]*i+j]) break;
frac = (j == 0 || j == head[1]) ? 0 :
(xval[0][k] - num) / (xval[0][k] - xval[0][k-1]);
mult[i-cip] = yval[0][k-1] * frac + yval[0][k] * (1-frac);
}
i = ((mult[0] * (1-cfrac) + mult[1] * cfrac) * row + num) * 2;
RAW(row,col) = LIM(i,0,65535);
}
free (yval[0]);
}
} }
void CLASS phase_one_load_raw() void CLASS phase_one_load_raw()