Tuesday, June 30, 2009

pseudoHDR - take 0.1.

I've been experimenting with reprocessing RAW files from my camera using dcraw and then taking the RGB data and reprocessing it.

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: