/* linear140.cpp - tests linearity of raw images on 140 checker */ static char *usage ="\n usage: %s in.png [plot.png]\n"; #include "imageLib.h" #include "math.h" #define VERBOSE 0 #define ROUND(x) ((x) >= 0? ((int)((x) + .5)) : ((int)((x) - .5))) #define SATURATION 251 /* ===================================== xyz-XRite2-d65.csv E6 xyz=(0.01 0.00740 0.01) L = 6.68 sRGB=( 20 21 21) J5 xyz=(0.02 0.02090 0.02) L = 15.95 sRGB=( 39 40 41) F6 xyz=(0.06 0.06510 0.07) L = 30.66 sRGB=( 72 72 73) I5 xyz=(0.08 0.08630 0.09) L = 35.26 sRGB=( 83 83 83) G6 xyz=(0.11 0.11340 0.12) L = 40.15 sRGB=( 94 95 95) K7 xyz=(0.14 0.14970 0.16) L = 45.59 sRGB=(108 108 107) H5 xyz=(0.17 0.18180 0.20) L = 49.71 sRGB=(118 118 118) H6 xyz=(0.27 0.28430 0.31) L = 60.28 sRGB=(145 145 146) G5 xyz=(0.32 0.34120 0.37) L = 65.06 sRGB=(158 158 158) K6 xyz=(0.40 0.41840 0.46) L = 70.76 sRGB=(173 173 174) I6 xyz=(0.46 0.48540 0.53) L = 75.16 sRGB=(185 185 185) F5 xyz=(0.53 0.55670 0.61) L = 79.43 sRGB=(197 197 197) J6 xyz=(0.70 0.74230 0.81) L = 89.03 sRGB=(222 224 224) E5 xyz=(0.86 0.91300 0.98) L = 96.53 sRGB=(243 246 243) ===================================== xyz-Nikolay1-d65.csv E6, 0.0066, 0.0070, 0.0074 J5, 0.0216, 0.0229, 0.0264 F6, 0.0628, 0.0662, 0.0730 I5, 0.0850, 0.0897, 0.0978 G6, 0.1110, 0.1172, 0.1286 K7, 0.1471, 0.1551, 0.1665 H5, 0.1742, 0.1836, 0.1970 H6, 0.2761, 0.2908, 0.3182 G5, 0.3281, 0.3460, 0.3766 K6, 0.4032, 0.4250, 0.4624 I6, 0.4680, 0.4926, 0.5347 F5, 0.5305, 0.5583, 0.6059 J6, 0.7053, 0.7451, 0.8101 E5, 0.8587, 0.9117, 0.9726 */ #define N 14 char *coords[] = {"E6", "J5", "F6", "I5", "G6", "K7", "H5", "H6", "G5", "K6", "I6", "F5", "J6", "E5"}; // XRite2 float lum1[] = {0.0074, 0.0209, 0.0651, 0.0863, 0.1134, 0.1497, 0.1818, 0.2843, 0.3412, 0.4184, 0.4854, 0.5567, 0.7423, 0.91300}; // Nikolay float lum2[] = {0.0070, 0.0229, 0.0662, 0.0897, 0.1172, 0.1551, 0.1836, 0.2908, 0.3460, 0.4250, 0.4926, 0.5583, 0.7451, 0.9117}; int main(int argc, char *argv[]) { float lum[N]; for (int i=0; i 2) plotfile = argv[2]; CByteImage im; ReadImageVerb(im, infile, VERBOSE); CShape sh = im.Shape(); int w = sh.width, h = sh.height; int b; int s = w/14; if (14*s != w || 10*s != h) throw CError("14x10 aspect ratio expected"); CShape sh2(N, 256, 3); CIntImage count(sh2); count.ClearPixels(); int maxcount = 0; // noise counters double totstd[] = {0, 0, 0}; int totn[] = {0, 0, 0}; // slope counters double sxx[] = {0, 0, 0}; double sxy[] = {0, 0, 0}; for (int k = 0; k < N; k++) { int cx = coords[k][0] - 'A'; int cy = coords[k][1] - '1'; // estimate noise (std) of each square: int nv[] = {0, 0, 0}; double nsv[] = {0, 0, 0}; double nsvv[] = {0, 0, 0}; for(int y = 0; y < s; y++){ for(int x = 0; x < s; x++){ int xi = x + s * cx; int yi = y + s * cy; for(b = 0; b < 3; b++){ int v = im.Pixel(xi, yi, b); if (v >= SATURATION) // skip saturated points continue; nv[b]++; nsv[b] += v; nsvv[b] += v*v; int c = count.Pixel(k, v, b); c++; if (c > maxcount) maxcount = c; count.Pixel(k, v, b) = c; double lu = lum[k] * 255; sxx[b] += lu * lu; sxy[b] += lu * v; } } } // accumulate total noise stddev (RMS from avg) for(b = 0; b < 3; b++) { int n = nv[b]; if (n > 0) { double avg = nsv[b] / n; double std = sqrt(nsvv[b]/n - avg*avg); totstd[b] += n * std; totn[b] += n; } } } // compute slopes and noise double sl[3], noise[3]; for(b = 0; b < 3; b++) { sl[b] = sxy[b] / sxx[b]; noise[b] = totstd[b] / totn[b]; } // compute RMS error of fit double serr[] = {0, 0, 0}; int sn[] = {0, 0, 0}; for (int k = 0; k < N; k++) { double x = (lum[k] * 255); for(b = 0; b < 3; b++){ double yfit = (sl[b]*x); for (int y = 0; y < 256; y++) { double err = y - yfit; int c = count.Pixel(k, y, b); serr[b] += c * err * err; sn[b] += c; } } } /* // for verification, different (long) way of computing RMSE double serr2[] = {0, 0, 0}; int sn2[] = {0, 0, 0}; for (int k = 0; k < N; k++) { double kx = (lum[k] * 255); int cx = coords[k][0] - 'A'; int cy = coords[k][1] - '1'; for(int y = 0; y < s; y++){ for(int x = 0; x < s; x++){ int xi = x + s * cx; int yi = y + s * cy; for(b = 0; b < 3; b++){ int v = im.Pixel(xi, yi, b); if (v >= SATURATION) // skip saturated points continue; double vfit = (sl[b]*kx); double err = v - vfit; sn2[b]++; serr2[b] += err * err; } } } } */ double rms[3]; for(b = 2; b >= 0; b--){ // 2==R, 1==G, 0==B rms[b] = sqrt(serr[b] / sn[b]); //double rms2 = sqrt(serr2[b] / sn2[b]); //printf("RMS=%3.1f/n=%3.1f: %3.1f ", // rms[b], noise[b], rms[b]/noise[b]); printf("%3.3f %3.3f %3.3f ", // RMS, noise, ratio rms[b], noise[b], rms[b]/noise[b]); } printf("\n"); if (plotfile) { CShape sh3(256, 256, 3); CByteImage plot(sh3); plot.ClearPixels(); // plot slopes: for(b = 0; b < 3; b++){ for (int x = 0; x < 256; x++) { int y = ROUND(sl[b]*x); if (y >= 0 && y < 256) // draw lines, brighter for blue plot.Pixel(x, 255-y, b) = (b==0) ? 100 : 70; // for paper make brighter: //plot.Pixel(x, 255-y, b) = 2.2*((b==0) ? 100 : 70); } } // plot datapoints for (int k = 0; k < N; k++) { int x = ROUND(lum[k] * 255); for (int y = 0; y < 256; y++) { for(b = 0; b < 3; b++){ double c = (double)count.Pixel(k, y, b) / maxcount; c = pow(c, 1.0/3.0); //c = c*2; // for paper make brighter int v = __min(255, ROUND(255 * c)); int x1 = __max(0, x-1); int x2 = __min(255, x+1); plot.Pixel(x, 255-y, b) = v; plot.Pixel(x1, 255-y, b) = v; plot.Pixel(x2, 255-y, b) = v; } } } WriteImageVerb(plot, plotfile, VERBOSE); } } catch (CError &err) { fprintf(stderr, err.message); fprintf(stderr, "\n"); return -1; } return 0; }