The core idea is that my dSLR has 14-bits of range (and your typical .jpg? 8-bit.) So that gives us extra dynamic range to reprocess the picture and do an HDR-type effect.
Here's my first pass at it. Yes, this code is structured horribly (if you're a future prospective employer, I can do better than this. Really!) I might or might not rewrite it properly once I get back from vacation. I've got other things I wanna throw together, too.
To use it, run dcraw with something like these settings: dcraw -h -a -4 -n 100 .CR2 and then pipe it through this program.
---
#include
#include
#include
double Kb = 0.114, Kr = 0.299;
char line1[512], line2[512], line3[512];
int i = 0, x, y, px, py;
unsigned short int *pic, *out;
double *fpic;
double *Y, *Yorig, *Y2, *Y3, *Cr, *Cg, *Cb;
double minY = 65536, maxY = 0.0;
void process(int inverse, double f1, double f2)
{
#define R 4
for (py = 0; py < y; py++) {
// fprintf(stderr, "%d\n", py);
for (px = 0; px < x; px++) {
double total = 0, peak = 0, mult, factor = 0, tfactor = 0;
int ty;
for (ty = ((py - R) > 0) ? (py - R) : 0; (ty < (y - 1)) && ((ty - py) < R); ty++) {
if (Y[(ty * x) + px] > peak) peak = Y[(ty * x) + px];
factor = 1.0 / ((abs(ty - py)) + 1);
total += (Y[(ty * x) + px] * factor);
tfactor += factor;
}
mult = (65536.0 - (total / tfactor)) / 65536.0;
mult = 1 - ((mult * mult) * f1);
if (mult > 4) mult = 4;
if (mult < 0.001) mult = 0.001;
Y2[(py * x) + px] = Y[(py * x) + px] / mult;
}
}
for (py = 0; py < y; py++) {
// fprintf(stderr, "%d\n", py);
for (px = 0; px < x; px++) {
double total = 0, peak = 0, mult, factor = 0, tfactor = 0;
int tx;
for (tx = ((px - R) > 0) ? (px - R) : 0; (tx < (x - 1)) && ((tx - px) < R); tx++) {
if (Y2[(py * x) + tx] > peak) peak = Y2[(py * x) + tx];
factor = 1.0 / ((abs(tx - px)) + 1);
total += (Y2[(py * x) + tx] * factor);
tfactor += factor;
}
mult = (65536.0 - (total / tfactor)) / 65536.0;
mult = 1 - ((mult * mult) * f2);
if (mult > 4) mult = 4;
if (mult < 0.001) mult = 0.001;
Y3[(py * x) + px] = Y2[(py * x) + px] / mult;
// fprintf(stderr, "%lf %lf %lf %lf\n", Y[(py * x) + px], Y2[(py * x)+ px], (total / tfactor), mult);
}
}
}void processdark(int inverse, double f1, double f2)
{
#if 0
for (py = 0; py < y; py++) {
// fprintf(stderr, "%d\n", py);
for (px = 0; px < x; px++) {
double total = 0, peak = 0, mult, factor = 0, tfactor = 0;
int ty;
for (ty = ((py - R) > 0) ? (py - R) : 0; (ty < (y - 1)) && ((ty - py) < R); ty++) {
if (Y[(ty * x) + px] > peak) peak = Y[(ty * x) + px];
factor = 1.0 / ((abs(ty - py)) + 1);
total += ((65536 - Y[(ty * x) + px]) * factor);
tfactor += factor;
}
mult = (65536.0 - (total / tfactor)) / 65536.0;
mult = 1 - ((mult * mult) * f1);
if (mult > 4) mult = 4;
if (mult <= 0) mult = 0;
Y2[(py * x) + px] = Y[(py * x) + px] * mult;
}
}
#else
memcpy(Y2, Y, sizeof(double) * x * y);
#endif
for (py = 0; py < y; py++) {
// fprintf(stderr, "%d\n", py);
for (px = 0; px < x; px++) {
double total = 0, peak = 0, mult, factor = 0, tfactor = 0;
int tx;
for (tx = ((px - R) > 0) ? (px - R) : 0; (tx < (x - 1)) && ((tx - px) < R); tx++) {
if (Y2[(py * x) + tx] > peak) peak = Y2[(py * x) + tx];
factor = 1.0 / ((abs(tx - px)) + 1);
total += ((Y[(py * x) + tx]) * factor);
tfactor += factor;
}
mult = (65536.0 - (total / tfactor)) / 65536.0;
mult = 1 - ((mult * mult) * f2);
if (mult > 4) mult = 4;
if (mult <= 0) mult = 0;
Y3[(py * x) + px] = Y2[(py * x) + px] * mult;
// fprintf(stderr, "%lf %lf %lf %lf\n", Y2[(py * x) + px], Y3[(py * x) + px], (total / tfactor), mult);
}
}
}
void processdark2(int inverse, double f1, double f2)
{
double total = 0, avg;
maxY = 0.0; minY = 65536.0;
for (i = 0; i < x * y; i++) {
if (Y[i] > maxY) maxY = Y[i];
if (Y[i] < minY) minY = Y[i];
if (Y[i] > 65535.0)
Y2[i] = 0.0;
else
Y2[i] = 65535.0 - Y[i];
}
maxY = 0.0; minY = 65536.0;
for (i = 0; i < x * y; i++) {
if (Y2[i] > maxY) maxY = Y2[i];
if (Y2[i] < minY) minY = Y2[i];
}
for (i = 0; i < x * y; i++) {
total += Y3[i] = 65535.0 - (Y2[i] * (65535.0 / maxY));
}
avg = total / (x * y);
for (i = 0; i < x * y; i++) {
if (Y3[i] > (avg * 4)) {
Y[i] = avg * 4;
} else {
Y[i] = Y3[i];
}
Y[i] = Y3[i];
}
}
int main(int argc, char *argv[])
{
double p1 = 0.8, p2 = 0.8;
if (argc >= 2) sscanf(argv[1], "%lf", &p1);
if (argc >= 3) sscanf(argv[2], "%lf", &p2);
memset(line1, 512, 0);
/* read the first line */
while (read(0, &line1[i], 1)) {
if (line1[i] == '\n') break;
i++;
}
line1[i + 1] = 0;
memset(line2, 512, 0);
i = 0;
/* read the second line */
while (read(0, &line2[i], 1)) {
if (line2[i] == '\n') break;
i++;
}
line2[i + 1] = 0;
sscanf(line2, "%d %d", &x, &y);
memset(line3, 512, 0);
i = 0;
/* read the third line */
while (read(0, &line3[i], 1)) {
if (line3[i] == '\n') break;
i++;
}
line3[i + 1] = 0;
pic = (unsigned short *)malloc(x * y * 6);
out = (unsigned short *)malloc(x * y * 6);
fpic = (double *)malloc(x * y * 3 * sizeof(double));
read(0, pic, (x * y * 6));
memcpy(out, pic, (x * y * 6));
Y = (double *)malloc(x * y * sizeof(double));
Yorig = (double *)malloc(x * y * sizeof(double));
Y2 = (double *)malloc(x * y * sizeof(double));
Y3 = (double *)malloc(x * y * sizeof(double));
Cr = (double *)malloc(x * y * sizeof(double));
Cg = (double *)malloc(x * y * sizeof(double));
Cb = (double *)malloc(x * y * sizeof(double));
for (i = 0; i < (x * y * 3); i++) {
fpic[i] = ntohs(pic[i]);
}
for (i = 0; i < x * y; i++) {
double r = fpic[(i * 3)], g = fpic[(i * 3) + 1], b = fpic[(i * 3) + 2];
Y[i] = (0.299 * r) + (0.587 * g) + (0.114 * b);
Yorig[i] = (0.299 * r) + (0.587 * g) + (0.114 * b);
Cr[i] = r / Y[i];
Cg[i] = g / Y[i];
Cb[i] = b / Y[i];
// Cr[i] = -(0.168736 * r) - (0.331264 * g) + (0.5 * b);
// Cb[i] = +(0.5 * r) - (0.418688 * g) - (0.081312 * b);
if (Y[i] > maxY) maxY = Y[i];
if (Y[i] < minY) minY = Y[i];
}
fprintf(stderr, "%lf %lf\n", maxY, minY);
processdark2(1, 0.8, 0.8);
process(0, p1, p1);
memcpy(Y, Y3, (x * y * sizeof(double)));
maxY = 0.0; minY = 65536.0;
for (i = 0; i < x * y; i++) {
if (Y[i] > maxY) maxY = Y[i];
if (Y[i] < minY) minY = Y[i];
}
fprintf(stderr, "%lf %lf\n", maxY, minY);
processdark(1, p2, p2);
for (i = 0; i < x * y; i++) {
double r1 = fpic[i * 3], g1 = fpic[(i * 3) + 1], b1 = fpic[(i * 3) + 2];
// fprintf(stderr, "%lf %lf %lf ", fpic[i * 3], fpic[(i * 3) + 1], fpic[(i * 3) + 2]);
fpic[(i * 3)] = r1 * (Y3[i] / Yorig[i]);
fpic[(i * 3) + 1] = g1 * (Y3[i] / Yorig[i]);
fpic[(i * 3) + 2] = b1 * (Y3[i] / Yorig[i]);
double r2 = fpic[i * 3], g2 = fpic[(i * 3) + 1], b2 = fpic[(i * 3) + 2];
double max = r2;
if (g2 > max) max = g2;
if (b2 > max) max = b2;
if (max > 65535.0) {
Y3[i] *= (65535.0 / max);
fpic[(i * 3)] = r1 * (Y3[i] / Yorig[i]);
fpic[(i * 3) + 1] = g1 * (Y3[i] / Yorig[i]);
fpic[(i * 3) + 2] = b1 * (Y3[i] / Yorig[i]);
// fprintf(stderr, "%lf %lf %lf %lf, ", Y3[i], fpic[(i * 3)], fpic[(i* 3) + 1], fpic[(i * 3) + 2]);
// fprintf(stderr, "%lf %lf %lf, %d %d %d\n", fpic[(i * 3)], fpic[(i * 3) + 1], fpic[(i * 3) + 2], pic[i * 3], pic[(i * 3) + 1], pic[(i * 3) + 2]);
}
// if (fabs(g2 - g1) > 33)
// fprintf(stderr, "%lf %lf %lf\n", r2 - r1, g2 - g1, b2 - b1);
}
maxY = 0.0; minY = 65536.0;
for (i = 0; i < x * y; i++) {
if (Y3[i] > maxY) maxY = Y3[i];
if (Y3[i] < minY) minY = Y3[i];
}
fprintf(stderr, "%lf %lf\n", maxY, minY);
for (i = 0; i < (x * y * 3); i++) {
if (fpic[i] < 0) fpic[i] = 0;
if (fpic[i] > 65535) fpic[i] = 65535;
out[i] = htons((unsigned short)fpic[i]);
// fprintf(stderr, "%d %d %lf\n", pic[i], out[i], fpic[i]);
}
write(1, line1, strlen(line1));
write(1, line2, strlen(line2));
write(1, line3, strlen(line3));
write(1, out, x * y * 6);
return 0;
}
No comments:
Post a Comment