From 87b107a48ca4d4c2954a73a79c98233c5e72d65a Mon Sep 17 00:00:00 2001 From: torger Date: Mon, 29 Dec 2014 20:30:56 +0100 Subject: [PATCH] Added Hasselblad 3FR calibration data parser. --- rtengine/dcraw.cc | 325 +++++++++++++++++++++++++++++++ rtengine/dcraw.h | 6 + rtengine/dcraw.patch | 451 ++++++++++++++++++++++++++++++++++++++----- 3 files changed, 729 insertions(+), 53 deletions(-) diff --git a/rtengine/dcraw.cc b/rtengine/dcraw.cc index fb3f555f6..30750326a 100644 --- a/rtengine/dcraw.cc +++ b/rtengine/dcraw.cc @@ -1676,6 +1676,326 @@ void CLASS phase_one_load_raw_c() maximum = 0xfffc - ph1.black; } +void CLASS parse_hasselblad_gain() +{ + /* + Reverse-engineer's notes: + + The Hasselblad gain tag (0x19 in makernotes) is only available in the 3FR format and + is applied and removed when Hasselblad Phocus converts it to the FFF format. It's + always 0x300000 bytes large regardless of (tested) model, not all space in it is + used though. + + It contains individual calibration information from the factory to tune the sensor + performance. + + There is more calibration data in the tag than what is applied in conversion to FFF, + I've only cared to figure out the data which is actually used, but have some leads on + remaining info. + + The format is not equal between all models. Due to lack of 3FR files (harder to get + than FFF) only a subset of the models have been reverse-engineered. + + The header space is 512 bytes and is a mix of 16 and 32 bit values. Offset to blocks + are 32 bit values, but all information seems to be encoded in 16 bit values. Many + values in the header can be zeroed with no effect on FFF conversion, and their + meaning, if any, have not been further investigated. + + Formats: + hdr16[22] = raw width + hdr16[23] = raw height + hdr32[24] = offset to level corr block + Data block format. Seems to differ depending on sensor type. For tested H4D-50 + and H3D-31: 10 16 bit signed values per row + value 0: a correction factor (k) used on even columns, where the new pixel value is + calulated as follows: + new_value = old_value + (2 * ((k * (old_value_on_row_above-256)) / 32767) - 2) + note the connection to the value on the row above, seems to be some sort of signal + leakage correction. + value 1: same as value 0 but using old value on row below instead of above + value 2-3: same as value 0-1 but for odd columns + value 4-9: has some effect if non-zero (probably similar to the others) but not + investigated which, as it's seems to be always zero for the tested cameras. + + hdr32[25] = probably offset to "bad/unreliable pixels" info, always 512 as it starts + directly after the header. Not applied in FFF conversion (at least + typically). + Data block format guess: raw_height packets of + + + hdr32[27] = offset to unknown data (bad colulmns?), of the form: + <0> + packet: . + + hdr32[34] = curves offset, seems to be A/D curves (one per channel) on newer models + and some sort of a film curve on older. Not applied in FFF conversion. + + hdr32[36] = flatfield correction, not available in older models. Data format: + <1><11 * 2 pad> + packet: + + The header pad is not zeroed and might seem to contain some sort of + information, but it makes no difference if set to zero. See + hasselblad_correct() how the flatfield is applied. + + Applied in FFF conversion is flatfield correction, levels, and the bad columns(?) + data. A/D curves are surprisingly not applied, maybe pre-applied in hardware and + only available as information? + + Not all combinations/models have been tested so there may be gaps. + */ + + ushort raw_h, count, ch_count, u16; + int i, offset; + off_t base; + + base = ftell(ifp); + fseek(ifp, 2 * 23, SEEK_CUR); + raw_h = get2(); + fseek(ifp, 48, SEEK_CUR); + offset = get4(); + hbd.levels = offset ? base + offset : 0; + fseek(ifp, 8, SEEK_CUR); + offset = get4(); + hbd.unknown1 = offset ? base + offset : 0; + fseek(ifp, 32, SEEK_CUR); + offset = get4(); + hbd.flatfield = offset ? base + offset : 0; +} + +void CLASS hasselblad_correct() +{ + unsigned col, row; + + /* + + This function applies 3FR calibration data. At the time of writing it supports a + subset, so here's a todo list: + + TODO: + - Support all gain tag formats + - The 0x19 tag varies a bit between models. We don't have parsers for all models. + - The reference model used was during inital reverse-engineering was a H4D-50, + we probably support the Hasselblads with Kodak 31, 39, 40 and 50 megapixels + well, but more work is needed for Kodak 16 and 22, Dalsa 60 and Sony 50. + - Apply bad column(?) data (hbd.unknown1) + - It was left out in this initial release since the effect is very small. + - Apply black(?) data, tag 0x1a and 0x1b is not parsed, it has however marginal + effect (at least for shorter exposures) so it's not too important. + + While there are things to do, the current implementation supports the most + important aspects, the faltfield and levels calibrations applied can have strong + visible effects. + + */ + + if (hbd.flatfield) { + int bw, bh, ffrows, ffcols, i, c; + ushort ref[4], ref_max; + fseek(ifp, hbd.flatfield, SEEK_SET); + get2(); + bw = get2(); + bh = get2(); + ffcols = get2(); + ffrows = get2(); + fseek(ifp, hbd.flatfield + 16 * 2, SEEK_SET); + + ushort *ffmap = (ushort *)malloc(sizeof(*ffmap) * 4 * ffcols * ffrows); + for (i = 0; i < 4 * ffcols * ffrows; i++) ffmap[i] = get2(); + + /* Get reference values from center of field. This seems to be what Phocus does too, + but haven't cared to figure out exactly at which coordinate */ + i = 4 * (ffcols * ffrows / 2 + ffcols / 2); + ref[0] = ffmap[i+0]; + ref[1] = ffmap[i+1]; + ref[3] = ffmap[i+2]; // G2 = index 3 in dcraw, 2 in 3FR + ref[2] = ffmap[i+3]; + ref_max = 0; + FORC4 if (ref[c] > ref_max) ref_max = ref[c]; + if (ref_max == 0) ref[0] = ref[1] = ref[2] = ref[3] = ref_max = 10000; + + /* Convert measured flatfield values to actual multipliers. The measured flatfield + can have vignetting which should be normalized, only color cast should be corrected. */ + for (i = 0; i < 4 * ffcols * ffrows; i += 4) { + double base, min = 65535.0, max = 0; + double cur[4]; + cur[0] = (double)ffmap[i+0] / ref[0]; + cur[1] = (double)ffmap[i+1] / ref[1]; + cur[3] = (double)ffmap[i+2] / ref[3]; // G2 index differs in dcraw and 3FR + cur[2] = (double)ffmap[i+3] / ref[2]; + FORC4 { + if (cur[c] < min) min = cur[c]; + if (cur[c] > max) max = cur[c]; + } + if (max == 0) max = 1.0; + base = (cur[0]+cur[1]+cur[2]+cur[3])/(max*4); + FORC4 cur[c] = cur[c] == 0 ? 1.0 : (base * max) / cur[c]; + /* convert to integer multiplier and store back to ffmap, we limit + range to 4 (16384*4) which should be fine for flatfield */ + FORC4 { + cur[c] *= 16384.0; + if (cur[c] > 65535.0) cur[c] = 65535.0; + ffmap[i+c] = (ushort)cur[c]; + } + } + + // of the cameras we've tested we know the exact placement of the flatfield map + int row_offset, col_offset; + switch (raw_width) { + case 8282: // 50 megapixel Kodak + row_offset = 21; + col_offset = 71; + break; + default: + /* Default case for camera models we've not tested. We center the map, which may + not be exactly where it should be but close enough for the smooth flatfield */ + row_offset = (raw_height - bh * ffrows) / 2; + col_offset = (raw_width - bw * ffcols) / 2; + break; + } + + /* + Concerning smoothing between blocks in the map Phocus makes it only vertically, + probably because it's simpler and faster. Looking at actual flatfield data it seems + like it's better to smooth in both directions. Possibly flatfield could be used for + correcting tiling on Dalsa sensors (H4D-60) like partly done in Phase One IIQ format, + and then sharp edges may be beneficial at least at the tiling seams, but at the time + of writing I've had no H4D-60 3FR files to test with to verify that. + + Meanwhile we do both vertical and horizontal smoothing/blurring. + */ + + /* pre-calculate constants for blurring. We probably should make a more efficient + blur than this, but this does not need any buffer and makes nice-looking + radial gradients */ + ushort *corners_weight = (ushort *)malloc(bw*bh*9*sizeof(*corners_weight)); + const char corners_mix[9][4][2] = { { {0,0}, {0,1}, {1,0}, {1,1} }, + { {0,1}, {1,1}, {-1,-1}, {-1,-1} }, + { {0,1}, {0,2}, {1,1}, {1,2} }, + { {1,0}, {1,1}, {-1,-1}, {-1,-1} }, + { {1,1}, {-1,-1}, {-1,-1}, {-1,-1} }, + { {1,1}, {1,2}, {-1,-1}, {-1,-1} }, + { {1,0}, {1,1}, {2,0}, {2,1} }, + { {1,1}, {2,1}, {-1,-1}, {-1,-1} }, + { {1,1}, {1,2}, {2,1}, {2,2} } }; + const ushort corners_shift[9] = { 2, 1, 2, 1, 0, 1, 2, 1, 2 }; + for (row = 0; row < bh; row++) { + const ushort maxdist = bw < bh ? bw/2-1 : bh/2-1; + const unsigned corners[9][2] = {{0,0}, {0,bw/2}, {0,bw-1}, + {bh/2,0},{bh/2,bw/2},{bh/2,bw-1}, + {bh-1,0},{bh-1,bw/2},{bh-1,bw-1}}; + for (col = 0; col < bw; col++) { + for (i = 0; i < 9; i++) { + ushort dist = (ushort)sqrt(abs(corners[i][0] - row) * abs(corners[i][0] - row) + abs(corners[i][1] - col) * abs(corners[i][1] - col)); + ushort weight = dist > maxdist ? 0 : maxdist - dist; + corners_weight[9*(row*bw+col)+i] = weight; + } + } + } + + // apply flatfield +#pragma omp parallel for + for (int row = 0; row < raw_height; row++) { + int ffs, cur_ffr, i, c; + if (row < row_offset) { + cur_ffr = row_offset; + ffs = 0; + } else if (row >= row_offset + ffrows * bh) { + cur_ffr = row_offset + (ffrows-1) * bh; + ffs = 4 * ffcols * (ffrows-1); + } else { + cur_ffr = row_offset + bh * ((row - row_offset) / bh); + ffs = 4 * ffcols * ((row - row_offset) / bh); + } + int next_ffc = 0, cur_ffc = col_offset; + int ffc = ffs; + ushort *cur[3][3]; // points to local ffmap entries with center at cur[1][1] + for (int col = 0; col < raw_width; col++) { + if (col == next_ffc) { + int rowsub = ffs == 0 ? 0 : ffcols*4; + int rowadd = ffs == 4 * ffcols * (ffrows-1) ? 0 : ffcols * 4; + int colsub = ffc == ffs ? 0 : 4; + int coladd = ffc == ffs + 4 * (ffcols-1) ? 0 : 4; + if (col != 0) cur_ffc = next_ffc; + else next_ffc += col_offset; + next_ffc += bw; + + cur[0][0] = &ffmap[ffc-rowsub-colsub]; + cur[0][1] = &ffmap[ffc-rowsub]; + cur[0][2] = &ffmap[ffc-rowsub+coladd]; + + cur[1][0] = &ffmap[ffc-colsub]; + cur[1][1] = &ffmap[ffc]; + cur[1][2] = &ffmap[ffc+coladd]; + + cur[2][0] = &ffmap[ffc+rowadd-colsub]; + cur[2][1] = &ffmap[ffc+rowadd]; + cur[2][2] = &ffmap[ffc+rowadd+coladd]; + + ffc += 4; + if (ffc == ffs + 4 * ffcols) next_ffc += raw_width; // last col in map, avoid stepping further + } + unsigned v = RAW(row,col); + if (v > black) { + c = FC(row,col); + unsigned x = col < cur_ffc ? 0 : col - cur_ffc; + unsigned y = row < cur_ffr ? 0 : row - cur_ffr; + if (x >= bw) x = bw-1; + if (y >= bh) y = bh-1; + unsigned wsum = 0; + unsigned mul = 0; + for (i = 0; i < 9; i++) { + ushort cw = corners_weight[9*(y*bw+x)+i]; + if (cw) { + unsigned m = 0; + int j; + for (j = 0; j < 1 << corners_shift[i]; j++) { + int cr = corners_mix[i][j][0], cc = corners_mix[i][j][1]; + m += cur[cr][cc][c]; + } + m >>= corners_shift[i]; + mul += m * cw; + wsum += cw; + } + } + mul /= wsum; + v = black + ((v-black) * mul) / 16384; + RAW(row,col) = v > 65535 ? 65535 : v; + } + } + } + free(ffmap); + free(corners_weight); + } + + if (hbd.levels) { + int i; + fseek(ifp, hbd.levels, SEEK_SET); + /* skip the first set (not used as we don't apply on first/last row), we look at it though to see if + the levels format is one that we support (there are other formats on some models which is not + supported here) */ + short test[10]; + for (i = 0; i < 10; i++) test[i] = (short)get2(); + if (test[5] == 0 && test[6] == 0 && test[7] == 0 && test[8] == 0 && test[9] == 0) { + int corr[4]; + ushort *row_above = (ushort *)malloc(sizeof(ushort) * raw_width); // we need to cache row above as we write new values as we go + for (col = 0; col < raw_width; col++) row_above[col] = RAW(0,col); + for (row = 1; row < raw_height-1; row++) { + for (i = 0; i < 4; i++) corr[i] = (short)get2(); + fseek(ifp, 6 * 2, SEEK_CUR); + for (col = 0; col < raw_width; col++) { + unsigned v = RAW(row,col); + if (corr[((col & 1)<<1)+0] && row_above[col] > black) v += 2 * ((corr[((col & 1)<<1)+0] * (row_above[col]-(int)black)) / 32767) - 2; + if (corr[((col & 1)<<1)+1] && RAW(row+1,col) > black) v += 2 * ((corr[((col & 1)<<1)+1] * (RAW(row+1,col)-(int)black)) / 32767) - 2; + row_above[col] = RAW(row,col); + RAW(row,col) = CLIP(v); + } + } + free(row_above); + } + } +} + void CLASS hasselblad_load_raw() { struct jhead jh; @@ -3536,6 +3856,8 @@ void CLASS crop_masked_pixels() if (load_raw == &CLASS phase_one_load_raw || load_raw == &CLASS phase_one_load_raw_c) phase_one_correct(); + if (load_raw == &CLASS hasselblad_load_raw) // RT + hasselblad_correct(); // RT if (fuji_width) { for (row=0; row < raw_height-top_margin*2; row++) { for (col=0; col < fuji_width << !fuji_layout; col++) { @@ -4755,6 +5077,9 @@ nf: order = 0x4949; if (tag == 0x1b) tag = 0x1018; if (tag == 0x1c) tag = 0x1017; } + if (tag == 0x19 && !strcmp(make,"Hasselblad") && len == 0x300000) { // RT + parse_hasselblad_gain(); // RT + } // RT if (tag == 0x1d) while ((c = fgetc(ifp)) && c != EOF) serial = serial*10 + (isdigit(c) ? c - '0' : c % 10); diff --git a/rtengine/dcraw.h b/rtengine/dcraw.h index f9e6389ce..3bfda4c1f 100644 --- a/rtengine/dcraw.h +++ b/rtengine/dcraw.h @@ -62,6 +62,7 @@ public: ,pana_bits(ifp,load_flags) ,float_raw_image(NULL) { + memset(&hbd, 0, sizeof(hbd)); aber[0]=aber[1]=aber[2]=aber[3]=1; gamm[0]=0.45;gamm[1]=4.5;gamm[2]=gamm[3]=gamm[4]=gamm[5]=0; user_mul[0]=user_mul[1]=user_mul[2]=user_mul[3]=0; @@ -133,6 +134,9 @@ protected: int black, split_col, black_col, split_row, black_row; float tag_210; } ph1; + struct hbd { + off_t levels, unknown1, flatfield; + } hbd; struct jhead { int bits, high, wide, clrs, sraw, psv, restart, vpred[6]; @@ -254,6 +258,8 @@ private: ph1_bithuff_t ph1_bithuff; void phase_one_load_raw_c(); +void hasselblad_correct(); +void parse_hasselblad_gain(); void hasselblad_load_raw(); void leaf_hdr_load_raw(); void unpacked_load_raw(); diff --git a/rtengine/dcraw.patch b/rtengine/dcraw.patch index e4f0ec391..67e679847 100755 --- a/rtengine/dcraw.patch +++ b/rtengine/dcraw.patch @@ -1,5 +1,5 @@ --- dcraw.c 2014-07-24 16:15:36.700261700 +0200 -+++ dcraw.cc 2014-12-17 20:56:37.864321656 +0100 ++++ dcraw.cc 2014-12-29 20:24:33.459795490 +0100 @@ -1,3 +1,15 @@ +/*RT*/#include +/*RT*/#include @@ -245,7 +245,334 @@ unsigned c; if (nbits == -1) -@@ -1903,10 +1867,10 @@ +@@ -1712,6 +1676,326 @@ + maximum = 0xfffc - ph1.black; + } + ++void CLASS parse_hasselblad_gain() ++{ ++ /* ++ Reverse-engineer's notes: ++ ++ The Hasselblad gain tag (0x19 in makernotes) is only available in the 3FR format and ++ is applied and removed when Hasselblad Phocus converts it to the FFF format. It's ++ always 0x300000 bytes large regardless of (tested) model, not all space in it is ++ used though. ++ ++ It contains individual calibration information from the factory to tune the sensor ++ performance. ++ ++ There is more calibration data in the tag than what is applied in conversion to FFF, ++ I've only cared to figure out the data which is actually used, but have some leads on ++ remaining info. ++ ++ The format is not equal between all models. Due to lack of 3FR files (harder to get ++ than FFF) only a subset of the models have been reverse-engineered. ++ ++ The header space is 512 bytes and is a mix of 16 and 32 bit values. Offset to blocks ++ are 32 bit values, but all information seems to be encoded in 16 bit values. Many ++ values in the header can be zeroed with no effect on FFF conversion, and their ++ meaning, if any, have not been further investigated. ++ ++ Formats: ++ hdr16[22] = raw width ++ hdr16[23] = raw height ++ hdr32[24] = offset to level corr block ++ Data block format. Seems to differ depending on sensor type. For tested H4D-50 ++ and H3D-31: 10 16 bit signed values per row ++ value 0: a correction factor (k) used on even columns, where the new pixel value is ++ calulated as follows: ++ new_value = old_value + (2 * ((k * (old_value_on_row_above-256)) / 32767) - 2) ++ note the connection to the value on the row above, seems to be some sort of signal ++ leakage correction. ++ value 1: same as value 0 but using old value on row below instead of above ++ value 2-3: same as value 0-1 but for odd columns ++ value 4-9: has some effect if non-zero (probably similar to the others) but not ++ investigated which, as it's seems to be always zero for the tested cameras. ++ ++ hdr32[25] = probably offset to "bad/unreliable pixels" info, always 512 as it starts ++ directly after the header. Not applied in FFF conversion (at least ++ typically). ++ Data block format guess: raw_height packets of ++ ++ ++ hdr32[27] = offset to unknown data (bad colulmns?), of the form: ++ <0> ++ packet: . ++ ++ hdr32[34] = curves offset, seems to be A/D curves (one per channel) on newer models ++ and some sort of a film curve on older. Not applied in FFF conversion. ++ ++ hdr32[36] = flatfield correction, not available in older models. Data format: ++ <1><11 * 2 pad> ++ packet: ++ ++ The header pad is not zeroed and might seem to contain some sort of ++ information, but it makes no difference if set to zero. See ++ hasselblad_correct() how the flatfield is applied. ++ ++ Applied in FFF conversion is flatfield correction, levels, and the bad columns(?) ++ data. A/D curves are surprisingly not applied, maybe pre-applied in hardware and ++ only available as information? ++ ++ Not all combinations/models have been tested so there may be gaps. ++ */ ++ ++ ushort raw_h, count, ch_count, u16; ++ int i, offset; ++ off_t base; ++ ++ base = ftell(ifp); ++ fseek(ifp, 2 * 23, SEEK_CUR); ++ raw_h = get2(); ++ fseek(ifp, 48, SEEK_CUR); ++ offset = get4(); ++ hbd.levels = offset ? base + offset : 0; ++ fseek(ifp, 8, SEEK_CUR); ++ offset = get4(); ++ hbd.unknown1 = offset ? base + offset : 0; ++ fseek(ifp, 32, SEEK_CUR); ++ offset = get4(); ++ hbd.flatfield = offset ? base + offset : 0; ++} ++ ++void CLASS hasselblad_correct() ++{ ++ unsigned col, row; ++ ++ /* ++ ++ This function applies 3FR calibration data. At the time of writing it supports a ++ subset, so here's a todo list: ++ ++ TODO: ++ - Support all gain tag formats ++ - The 0x19 tag varies a bit between models. We don't have parsers for all models. ++ - The reference model used was during inital reverse-engineering was a H4D-50, ++ we probably support the Hasselblads with Kodak 31, 39, 40 and 50 megapixels ++ well, but more work is needed for Kodak 16 and 22, Dalsa 60 and Sony 50. ++ - Apply bad column(?) data (hbd.unknown1) ++ - It was left out in this initial release since the effect is very small. ++ - Apply black(?) data, tag 0x1a and 0x1b is not parsed, it has however marginal ++ effect (at least for shorter exposures) so it's not too important. ++ ++ While there are things to do, the current implementation supports the most ++ important aspects, the faltfield and levels calibrations applied can have strong ++ visible effects. ++ ++ */ ++ ++ if (hbd.flatfield) { ++ int bw, bh, ffrows, ffcols, i, c; ++ ushort ref[4], ref_max; ++ fseek(ifp, hbd.flatfield, SEEK_SET); ++ get2(); ++ bw = get2(); ++ bh = get2(); ++ ffcols = get2(); ++ ffrows = get2(); ++ fseek(ifp, hbd.flatfield + 16 * 2, SEEK_SET); ++ ++ ushort *ffmap = (ushort *)malloc(sizeof(*ffmap) * 4 * ffcols * ffrows); ++ for (i = 0; i < 4 * ffcols * ffrows; i++) ffmap[i] = get2(); ++ ++ /* Get reference values from center of field. This seems to be what Phocus does too, ++ but haven't cared to figure out exactly at which coordinate */ ++ i = 4 * (ffcols * ffrows / 2 + ffcols / 2); ++ ref[0] = ffmap[i+0]; ++ ref[1] = ffmap[i+1]; ++ ref[3] = ffmap[i+2]; // G2 = index 3 in dcraw, 2 in 3FR ++ ref[2] = ffmap[i+3]; ++ ref_max = 0; ++ FORC4 if (ref[c] > ref_max) ref_max = ref[c]; ++ if (ref_max == 0) ref[0] = ref[1] = ref[2] = ref[3] = ref_max = 10000; ++ ++ /* Convert measured flatfield values to actual multipliers. The measured flatfield ++ can have vignetting which should be normalized, only color cast should be corrected. */ ++ for (i = 0; i < 4 * ffcols * ffrows; i += 4) { ++ double base, min = 65535.0, max = 0; ++ double cur[4]; ++ cur[0] = (double)ffmap[i+0] / ref[0]; ++ cur[1] = (double)ffmap[i+1] / ref[1]; ++ cur[3] = (double)ffmap[i+2] / ref[3]; // G2 index differs in dcraw and 3FR ++ cur[2] = (double)ffmap[i+3] / ref[2]; ++ FORC4 { ++ if (cur[c] < min) min = cur[c]; ++ if (cur[c] > max) max = cur[c]; ++ } ++ if (max == 0) max = 1.0; ++ base = (cur[0]+cur[1]+cur[2]+cur[3])/(max*4); ++ FORC4 cur[c] = cur[c] == 0 ? 1.0 : (base * max) / cur[c]; ++ /* convert to integer multiplier and store back to ffmap, we limit ++ range to 4 (16384*4) which should be fine for flatfield */ ++ FORC4 { ++ cur[c] *= 16384.0; ++ if (cur[c] > 65535.0) cur[c] = 65535.0; ++ ffmap[i+c] = (ushort)cur[c]; ++ } ++ } ++ ++ // of the cameras we've tested we know the exact placement of the flatfield map ++ int row_offset, col_offset; ++ switch (raw_width) { ++ case 8282: // 50 megapixel Kodak ++ row_offset = 21; ++ col_offset = 71; ++ break; ++ default: ++ /* Default case for camera models we've not tested. We center the map, which may ++ not be exactly where it should be but close enough for the smooth flatfield */ ++ row_offset = (raw_height - bh * ffrows) / 2; ++ col_offset = (raw_width - bw * ffcols) / 2; ++ break; ++ } ++ ++ /* ++ Concerning smoothing between blocks in the map Phocus makes it only vertically, ++ probably because it's simpler and faster. Looking at actual flatfield data it seems ++ like it's better to smooth in both directions. Possibly flatfield could be used for ++ correcting tiling on Dalsa sensors (H4D-60) like partly done in Phase One IIQ format, ++ and then sharp edges may be beneficial at least at the tiling seams, but at the time ++ of writing I've had no H4D-60 3FR files to test with to verify that. ++ ++ Meanwhile we do both vertical and horizontal smoothing/blurring. ++ */ ++ ++ /* pre-calculate constants for blurring. We probably should make a more efficient ++ blur than this, but this does not need any buffer and makes nice-looking ++ radial gradients */ ++ ushort *corners_weight = (ushort *)malloc(bw*bh*9*sizeof(*corners_weight)); ++ const char corners_mix[9][4][2] = { { {0,0}, {0,1}, {1,0}, {1,1} }, ++ { {0,1}, {1,1}, {-1,-1}, {-1,-1} }, ++ { {0,1}, {0,2}, {1,1}, {1,2} }, ++ { {1,0}, {1,1}, {-1,-1}, {-1,-1} }, ++ { {1,1}, {-1,-1}, {-1,-1}, {-1,-1} }, ++ { {1,1}, {1,2}, {-1,-1}, {-1,-1} }, ++ { {1,0}, {1,1}, {2,0}, {2,1} }, ++ { {1,1}, {2,1}, {-1,-1}, {-1,-1} }, ++ { {1,1}, {1,2}, {2,1}, {2,2} } }; ++ const ushort corners_shift[9] = { 2, 1, 2, 1, 0, 1, 2, 1, 2 }; ++ for (row = 0; row < bh; row++) { ++ const ushort maxdist = bw < bh ? bw/2-1 : bh/2-1; ++ const unsigned corners[9][2] = {{0,0}, {0,bw/2}, {0,bw-1}, ++ {bh/2,0},{bh/2,bw/2},{bh/2,bw-1}, ++ {bh-1,0},{bh-1,bw/2},{bh-1,bw-1}}; ++ for (col = 0; col < bw; col++) { ++ for (i = 0; i < 9; i++) { ++ ushort dist = (ushort)sqrt(abs(corners[i][0] - row) * abs(corners[i][0] - row) + abs(corners[i][1] - col) * abs(corners[i][1] - col)); ++ ushort weight = dist > maxdist ? 0 : maxdist - dist; ++ corners_weight[9*(row*bw+col)+i] = weight; ++ } ++ } ++ } ++ ++ // apply flatfield ++#pragma omp parallel for ++ for (int row = 0; row < raw_height; row++) { ++ int ffs, cur_ffr, i, c; ++ if (row < row_offset) { ++ cur_ffr = row_offset; ++ ffs = 0; ++ } else if (row >= row_offset + ffrows * bh) { ++ cur_ffr = row_offset + (ffrows-1) * bh; ++ ffs = 4 * ffcols * (ffrows-1); ++ } else { ++ cur_ffr = row_offset + bh * ((row - row_offset) / bh); ++ ffs = 4 * ffcols * ((row - row_offset) / bh); ++ } ++ int next_ffc = 0, cur_ffc = col_offset; ++ int ffc = ffs; ++ ushort *cur[3][3]; // points to local ffmap entries with center at cur[1][1] ++ for (int col = 0; col < raw_width; col++) { ++ if (col == next_ffc) { ++ int rowsub = ffs == 0 ? 0 : ffcols*4; ++ int rowadd = ffs == 4 * ffcols * (ffrows-1) ? 0 : ffcols * 4; ++ int colsub = ffc == ffs ? 0 : 4; ++ int coladd = ffc == ffs + 4 * (ffcols-1) ? 0 : 4; ++ if (col != 0) cur_ffc = next_ffc; ++ else next_ffc += col_offset; ++ next_ffc += bw; ++ ++ cur[0][0] = &ffmap[ffc-rowsub-colsub]; ++ cur[0][1] = &ffmap[ffc-rowsub]; ++ cur[0][2] = &ffmap[ffc-rowsub+coladd]; ++ ++ cur[1][0] = &ffmap[ffc-colsub]; ++ cur[1][1] = &ffmap[ffc]; ++ cur[1][2] = &ffmap[ffc+coladd]; ++ ++ cur[2][0] = &ffmap[ffc+rowadd-colsub]; ++ cur[2][1] = &ffmap[ffc+rowadd]; ++ cur[2][2] = &ffmap[ffc+rowadd+coladd]; ++ ++ ffc += 4; ++ if (ffc == ffs + 4 * ffcols) next_ffc += raw_width; // last col in map, avoid stepping further ++ } ++ unsigned v = RAW(row,col); ++ if (v > black) { ++ c = FC(row,col); ++ unsigned x = col < cur_ffc ? 0 : col - cur_ffc; ++ unsigned y = row < cur_ffr ? 0 : row - cur_ffr; ++ if (x >= bw) x = bw-1; ++ if (y >= bh) y = bh-1; ++ unsigned wsum = 0; ++ unsigned mul = 0; ++ for (i = 0; i < 9; i++) { ++ ushort cw = corners_weight[9*(y*bw+x)+i]; ++ if (cw) { ++ unsigned m = 0; ++ int j; ++ for (j = 0; j < 1 << corners_shift[i]; j++) { ++ int cr = corners_mix[i][j][0], cc = corners_mix[i][j][1]; ++ m += cur[cr][cc][c]; ++ } ++ m >>= corners_shift[i]; ++ mul += m * cw; ++ wsum += cw; ++ } ++ } ++ mul /= wsum; ++ v = black + ((v-black) * mul) / 16384; ++ RAW(row,col) = v > 65535 ? 65535 : v; ++ } ++ } ++ } ++ free(ffmap); ++ free(corners_weight); ++ } ++ ++ if (hbd.levels) { ++ int i; ++ fseek(ifp, hbd.levels, SEEK_SET); ++ /* skip the first set (not used as we don't apply on first/last row), we look at it though to see if ++ the levels format is one that we support (there are other formats on some models which is not ++ supported here) */ ++ short test[10]; ++ for (i = 0; i < 10; i++) test[i] = (short)get2(); ++ if (test[5] == 0 && test[6] == 0 && test[7] == 0 && test[8] == 0 && test[9] == 0) { ++ int corr[4]; ++ ushort *row_above = (ushort *)malloc(sizeof(ushort) * raw_width); // we need to cache row above as we write new values as we go ++ for (col = 0; col < raw_width; col++) row_above[col] = RAW(0,col); ++ for (row = 1; row < raw_height-1; row++) { ++ for (i = 0; i < 4; i++) corr[i] = (short)get2(); ++ fseek(ifp, 6 * 2, SEEK_CUR); ++ for (col = 0; col < raw_width; col++) { ++ unsigned v = RAW(row,col); ++ if (corr[((col & 1)<<1)+0] && row_above[col] > black) v += 2 * ((corr[((col & 1)<<1)+0] * (row_above[col]-(int)black)) / 32767) - 2; ++ if (corr[((col & 1)<<1)+1] && RAW(row+1,col) > black) v += 2 * ((corr[((col & 1)<<1)+1] * (RAW(row+1,col)-(int)black)) / 32767) - 2; ++ row_above[col] = RAW(row,col); ++ RAW(row,col) = CLIP(v); ++ } ++ } ++ free(row_above); ++ } ++ } ++} ++ + void CLASS hasselblad_load_raw() + { + struct jhead jh; +@@ -1903,10 +2187,10 @@ maximum = curve[0x3ff]; } @@ -259,7 +586,7 @@ int byte; if (!nbits) return vbits=0; -@@ -2195,11 +2159,11 @@ +@@ -2195,11 +2479,11 @@ METHODDEF(boolean) fill_input_buffer (j_decompress_ptr cinfo) { @@ -273,7 +600,7 @@ cinfo->src->next_input_byte = jpeg_buffer; cinfo->src->bytes_in_buffer = nbytes; return TRUE; -@@ -2524,10 +2488,9 @@ +@@ -2524,10 +2808,9 @@ maximum = (1 << (thumb_misc & 31)) - 1; } @@ -286,7 +613,7 @@ if (start) { for (p=0; p < 4; p++) pad[p] = key = key * 48828125 + 1; -@@ -2612,11 +2575,13 @@ +@@ -2612,11 +2895,13 @@ bit += 7; } for (i=0; i < 16; i++, col+=2) @@ -301,7 +628,7 @@ } void CLASS samsung_load_raw() -@@ -2863,7 +2828,7 @@ +@@ -2863,7 +3148,7 @@ void CLASS foveon_decoder (unsigned size, unsigned code) { @@ -310,7 +637,16 @@ struct decode *cur; int i, len; -@@ -3586,10 +3551,13 @@ +@@ -3571,6 +3856,8 @@ + if (load_raw == &CLASS phase_one_load_raw || + load_raw == &CLASS phase_one_load_raw_c) + phase_one_correct(); ++ if (load_raw == &CLASS hasselblad_load_raw) // RT ++ hasselblad_correct(); // RT + if (fuji_width) { + for (row=0; row < raw_height-top_margin*2; row++) { + for (col=0; col < fuji_width << !fuji_layout; col++) { +@@ -3586,10 +3873,13 @@ } } } else { @@ -326,7 +662,7 @@ if (mask[0][3] > 0) goto mask_set; if (load_raw == &CLASS canon_load_raw || load_raw == &CLASS lossless_jpeg_load_raw) { -@@ -4191,239 +4159,8 @@ +@@ -4191,239 +4481,8 @@ } } @@ -375,7 +711,8 @@ - -/* - This algorithm is officially called: -- ++/* RT: delete interpolation functions */ + - "Interpolation using a Threshold-based variable number of gradients" - - described in http://scien.stanford.edu/pages/labsite/1999/psych221/projects/99/tingchen/algodep/vargra.html @@ -521,8 +858,7 @@ - - border_interpolate(3); - if (verbose) fprintf (stderr,_("PPG interpolation...\n")); -+/* RT: delete interpolation functions */ - +- -/* Fill in the green layer with gradients and pattern recognition: */ - for (row=3; row < height-3; row++) - for (col=3+(FC(row,3) & 1), c=FC(row,col); col < width-3; col+=2) { @@ -567,7 +903,7 @@ void CLASS cielab (ushort rgb[3], short lab[3]) { -@@ -4684,112 +4421,7 @@ +@@ -4684,112 +4743,7 @@ } #undef fcol @@ -585,7 +921,7 @@ - char (*homo)[TS][TS], *buffer; - - if (verbose) fprintf (stderr,_("AHD interpolation...\n")); - +- - cielab (0,0); - border_interpolate(5); - buffer = (char *) malloc (26*TS*TS); @@ -593,7 +929,7 @@ - rgb = (ushort(*)[TS][TS][3]) buffer; - lab = (short (*)[TS][TS][3])(buffer + 12*TS*TS); - homo = (char (*)[TS][TS]) (buffer + 24*TS*TS); -- + - for (top=2; top < height-5; top += TS-6) - for (left=2; left < width-5; left += TS-6) { - @@ -680,7 +1016,7 @@ #undef TS void CLASS median_filter() -@@ -4959,7 +4591,7 @@ +@@ -4959,7 +4913,7 @@ } } @@ -689,7 +1025,7 @@ void CLASS parse_makernote (int base, int uptag) { -@@ -5116,7 +4748,8 @@ +@@ -5116,12 +5070,16 @@ cam_mul[2] = get4() << 2; } } @@ -699,7 +1035,15 @@ fread (model, 64, 1, ifp); if (strstr(make,"PENTAX")) { if (tag == 0x1b) tag = 0x1018; -@@ -5367,7 +5000,7 @@ + if (tag == 0x1c) tag = 0x1017; + } ++ if (tag == 0x19 && !strcmp(make,"Hasselblad") && len == 0x300000) { // RT ++ parse_hasselblad_gain(); // RT ++ } // RT + if (tag == 0x1d) + while ((c = fgetc(ifp)) && c != EOF) + serial = serial*10 + (isdigit(c) ? c - '0' : c % 10); +@@ -5367,7 +5325,7 @@ { "","DCB2","Volare","Cantare","CMost","Valeo 6","Valeo 11","Valeo 22", "Valeo 11p","Valeo 17","","Aptus 17","Aptus 22","Aptus 75","Aptus 65", "Aptus 54S","Aptus 65S","Aptus 75S","AFi 5","AFi 6","AFi 7", @@ -708,7 +1052,7 @@ "","","","","Aptus-II 10R","Aptus-II 8","","Aptus-II 12","","AFi-II 12" }; float romm_cam[3][3]; -@@ -5456,6 +5089,8 @@ +@@ -5456,6 +5414,8 @@ wbi = -2; } if (tag == 2118) wbtemp = getint(type); @@ -717,7 +1061,7 @@ if (tag == 2130 + wbi) FORC3 mul[c] = getreal(type); if (tag == 2140 + wbi && wbi >= 0) -@@ -5475,8 +5110,8 @@ +@@ -5475,8 +5435,8 @@ } } @@ -728,7 +1072,7 @@ int CLASS parse_tiff_ifd (int base) { -@@ -5489,7 +5124,7 @@ +@@ -5489,7 +5449,7 @@ unsigned sony_curve[] = { 0,0,0,0,0,4095 }; unsigned *buf, sony_offset=0, sony_length=0, sony_key=0; struct jhead jh; @@ -737,7 +1081,7 @@ if (tiff_nifds >= sizeof tiff_ifd / sizeof tiff_ifd[0]) return 1; -@@ -5564,7 +5199,8 @@ +@@ -5564,7 +5524,8 @@ fgets (make, 64, ifp); break; case 272: /* Model */ @@ -747,7 +1091,7 @@ break; case 280: /* Panasonic RW2 offset */ if (type != 4) break; -@@ -5620,6 +5256,9 @@ +@@ -5620,6 +5581,9 @@ case 315: /* Artist */ fread (artist, 64, 1, ifp); break; @@ -757,7 +1101,7 @@ case 322: /* TileWidth */ tiff_ifd[ifd].tile_width = getint(type); break; -@@ -5635,6 +5274,9 @@ +@@ -5635,6 +5599,9 @@ is_raw = 5; } break; @@ -767,7 +1111,7 @@ case 330: /* SubIFDs */ if (!strcmp(model,"DSLR-A100") && tiff_ifd[ifd].width == 3872) { load_raw = &CLASS sony_arw_load_raw; -@@ -5648,6 +5290,9 @@ +@@ -5648,6 +5615,9 @@ fseek (ifp, i+4, SEEK_SET); } break; @@ -777,7 +1121,7 @@ case 400: strcpy (make, "Sarnoff"); maximum = 0xfff; -@@ -5828,6 +5473,9 @@ +@@ -5828,6 +5798,9 @@ if (!make[0]) strcpy (make, "DNG"); is_raw = 1; break; @@ -787,7 +1131,7 @@ case 50710: /* CFAPlaneColor */ if (filters == 9) break; if (len > 4) len = 4; -@@ -5859,10 +5507,21 @@ +@@ -5859,10 +5832,21 @@ case 61450: cblack[4] = cblack[5] = MIN(sqrt(len),64); case 50714: /* BlackLevel */ @@ -813,7 +1157,7 @@ case 50715: /* BlackLevelDeltaH */ case 50716: /* BlackLevelDeltaV */ for (num=i=0; i < len; i++) -@@ -5940,12 +5599,15 @@ +@@ -5940,12 +5924,15 @@ fread (buf, sony_length, 1, ifp); sony_decrypt (buf, sony_length/4, 1, sony_key); sfp = ifp; @@ -835,7 +1179,7 @@ ifp = sfp; free (buf); } -@@ -5969,6 +5631,7 @@ +@@ -5969,6 +5956,7 @@ int CLASS parse_tiff (int base) { int doff; @@ -843,7 +1187,7 @@ fseek (ifp, base, SEEK_SET); order = get2(); -@@ -6046,7 +5709,7 @@ +@@ -6046,7 +6034,7 @@ case 8: load_raw = &CLASS eight_bit_load_raw; break; case 12: if (tiff_ifd[raw].phint == 2) load_flags = 6; @@ -852,7 +1196,7 @@ case 14: load_flags = 0; case 16: load_raw = &CLASS unpacked_load_raw; if (!strncmp(make,"OLYMPUS",7) && -@@ -6079,6 +5742,7 @@ +@@ -6079,6 +6067,7 @@ case 32803: load_raw = &CLASS kodak_65000_load_raw; } case 32867: case 34892: break; @@ -860,7 +1204,7 @@ default: is_raw = 0; } if (!dng_version) -@@ -6164,7 +5828,7 @@ +@@ -6164,7 +6153,7 @@ { const char *file, *ext; char *jname, *jfile, *jext; @@ -869,7 +1213,7 @@ ext = strrchr (ifname, '.'); file = strrchr (ifname, '/'); -@@ -6186,13 +5850,14 @@ +@@ -6186,13 +6175,14 @@ } else while (isdigit(*--jext)) { if (*jext != '9') { @@ -886,7 +1230,7 @@ if (verbose) fprintf (stderr,_("Reading metadata from %s ...\n"), jname); parse_tiff (12); -@@ -6537,7 +6202,11 @@ +@@ -6537,7 +6527,11 @@ order = get2(); hlen = get4(); if (get4() == 0x48454150) /* "HEAP" */ @@ -898,7 +1242,7 @@ if (parse_tiff (save+6)) apply_tiff(); fseek (ifp, save+len, SEEK_SET); } -@@ -6789,7 +6458,8 @@ +@@ -6789,7 +6783,8 @@ { static const struct { const char *prefix; @@ -908,7 +1252,7 @@ } table[] = { { "AgfaPhoto DC-833m", 0, 0, /* DJC */ { 11438,-3762,-1115,-2409,9914,2497,-1227,2295,5300 } }, -@@ -7690,6 +7360,27 @@ +@@ -7690,6 +7685,27 @@ } break; } @@ -936,7 +1280,7 @@ } void CLASS simple_coeff (int index) -@@ -7967,7 +7658,7 @@ +@@ -7967,7 +7983,7 @@ tiff_flip = flip = filters = UINT_MAX; /* unknown */ raw_height = raw_width = fuji_width = fuji_layout = cr2_slice[0] = 0; maximum = height = width = top_margin = left_margin = 0; @@ -945,7 +1289,7 @@ iso_speed = shutter = aperture = focal_len = unique_id = 0; tiff_nifds = 0; memset (tiff_ifd, 0, sizeof tiff_ifd); -@@ -7999,13 +7690,20 @@ +@@ -7999,13 +8015,20 @@ fread (head, 1, 32, ifp); fseek (ifp, 0, SEEK_END); flen = fsize = ftell(ifp); @@ -968,7 +1312,7 @@ parse_ciff (hlen, flen-hlen, 0); load_raw = &CLASS canon_load_raw; } else if (parse_tiff(0)) apply_tiff(); -@@ -8051,6 +7749,7 @@ +@@ -8051,6 +8074,7 @@ fseek (ifp, 100+28*(shot_select > 0), SEEK_SET); parse_tiff (data_offset = get4()); parse_tiff (thumb_offset+12); @@ -976,7 +1320,7 @@ apply_tiff(); } else if (!memcmp (head,"RIFF",4)) { fseek (ifp, 0, SEEK_SET); -@@ -8160,15 +7859,18 @@ +@@ -8160,15 +8184,18 @@ if (make[0] == 0) parse_smal (0, flen); if (make[0] == 0) { parse_jpeg(0); @@ -1004,7 +1348,7 @@ } for (i=0; i < sizeof corp / sizeof *corp; i++) -@@ -8201,7 +7903,7 @@ +@@ -8201,7 +8228,7 @@ if (height == 3136 && width == 4864) /* Pentax K20D and Samsung GX20 */ { height = 3124; width = 4688; filters = 0x16161616; } if (width == 4352 && (!strcmp(model,"K-r") || !strcmp(model,"K-x"))) @@ -1013,7 +1357,7 @@ if (width >= 4960 && !strncmp(model,"K-5",3)) { left_margin = 10; width = 4950; filters = 0x16161616; } if (width == 4736 && !strcmp(model,"K-7")) -@@ -8220,6 +7922,7 @@ +@@ -8220,6 +8247,7 @@ switch (tiff_compress) { case 1: load_raw = &CLASS packed_dng_load_raw; break; case 7: load_raw = &CLASS lossless_dng_load_raw; break; @@ -1021,7 +1365,7 @@ case 34892: load_raw = &CLASS lossy_dng_load_raw; break; default: load_raw = 0; } -@@ -8347,7 +8050,7 @@ +@@ -8347,7 +8375,7 @@ width -= 44; } else if (!strcmp(model,"D3200") || !strcmp(model,"D600") || @@ -1030,7 +1374,7 @@ width -= 46; } else if (!strcmp(model,"D4") || !strcmp(model,"Df")) { -@@ -8567,24 +8270,46 @@ +@@ -8567,24 +8595,46 @@ if (load_raw == &CLASS lossless_jpeg_load_raw) load_raw = &CLASS hasselblad_load_raw; if (raw_width == 7262) { @@ -1079,7 +1423,7 @@ } else if (raw_width == 4090) { strcpy (model, "V96C"); height -= (top_margin = 6); -@@ -8637,6 +8362,7 @@ +@@ -8637,6 +8687,7 @@ filters = 0x16161616; } } else if (!strcmp(make,"Leica") || !strcmp(make,"Panasonic")) { @@ -1087,7 +1431,7 @@ if ((flen - data_offset) / (raw_width*8/7) == raw_height) load_raw = &CLASS panasonic_load_raw; if (!load_raw) { -@@ -8654,6 +8380,7 @@ +@@ -8654,6 +8705,7 @@ } filters = 0x01010101 * (uchar) "\x94\x61\x49\x16" [((filters-1) ^ (left_margin & 1) ^ (top_margin << 1)) & 3]; @@ -1095,7 +1439,7 @@ } else if (!strcmp(model,"C770UZ")) { height = 1718; width = 2304; -@@ -8883,6 +8610,14 @@ +@@ -8883,6 +8935,14 @@ memcpy (rgb_cam, cmatrix, sizeof cmatrix); raw_color = 0; } @@ -1110,7 +1454,7 @@ if (raw_color) adobe_coeff (make, model); if (load_raw == &CLASS kodak_radc_load_raw) if (raw_color) adobe_coeff ("Apple","Quicktake"); -@@ -8899,7 +8634,7 @@ +@@ -8899,7 +8959,7 @@ if (!tiff_bps) tiff_bps = 12; if (!maximum) maximum = (1 << tiff_bps) - 1; if (!load_raw || height < 22 || width < 22 || @@ -1119,7 +1463,7 @@ is_raw = 0; #ifdef NO_JASPER if (load_raw == &CLASS redcine_load_raw) { -@@ -8978,194 +8713,249 @@ +@@ -8978,195 +9038,250 @@ } #endif @@ -1532,19 +1876,20 @@ + + delete [] cBuffer; + delete [] uBuffer; - } ++} + } + + if (ifd->sample_format == 3) { // Floating point data + copyFloatDataToInt(float_raw_image, raw_image, raw_width*raw_height, maximum); + } -+} -+ -+/* RT: removed unused functions */ + } ++/* RT: removed unused functions */ ++ struct tiff_tag { ushort tag, type; -@@ -9188,585 +8978,12 @@ + int count; +@@ -9188,585 +9303,12 @@ unsigned gps[26]; char desc[512], make[64], model[64], soft[32], date[20], artist[64]; };