1 | /*****************************************************************************/ |
2 | // Copyright 2008 Adobe Systems Incorporated |
3 | // All Rights Reserved. |
4 | // |
5 | // NOTICE: Adobe permits you to use, modify, and distribute this file in |
6 | // accordance with the terms of the Adobe license agreement accompanying it. |
7 | /*****************************************************************************/ |
8 | |
9 | /* $Id: //mondo/dng_sdk_1_4/dng_sdk/source/dng_lens_correction.cpp#1 $ */ |
10 | /* $DateTime: 2012/05/30 13:28:51 $ */ |
11 | /* $Change: 832332 $ */ |
12 | /* $Author: tknoll $ */ |
13 | |
14 | /*****************************************************************************/ |
15 | |
16 | #include <cfloat> |
17 | #include <limits.h> |
18 | |
19 | #include "dng_1d_table.h" |
20 | #include "dng_assertions.h" |
21 | #include "dng_bottlenecks.h" |
22 | #include "dng_exceptions.h" |
23 | #include "dng_filter_task.h" |
24 | #include "dng_globals.h" |
25 | #include "dng_host.h" |
26 | #include "dng_image.h" |
27 | #include "dng_lens_correction.h" |
28 | #include "dng_negative.h" |
29 | #include "dng_safe_arithmetic.h" |
30 | #include "dng_sdk_limits.h" |
31 | #include "dng_tag_values.h" |
32 | #include "dng_utils.h" |
33 | |
34 | /*****************************************************************************/ |
35 | |
36 | dng_warp_params::dng_warp_params () |
37 | |
38 | : fPlanes (1) |
39 | , fCenter (0.5, 0.5) |
40 | |
41 | { |
42 | |
43 | } |
44 | |
45 | /*****************************************************************************/ |
46 | |
47 | dng_warp_params::dng_warp_params (uint32 planes, |
48 | const dng_point_real64 ¢er) |
49 | |
50 | : fPlanes (planes) |
51 | , fCenter (center) |
52 | |
53 | { |
54 | |
55 | DNG_ASSERT (planes >= 1, "Too few planes." ); |
56 | DNG_ASSERT (planes <= kMaxColorPlanes, "Too many planes." ); |
57 | |
58 | DNG_ASSERT (fCenter.h >= 0.0 && fCenter.h <= 1.0, |
59 | "Center (horizontal) out of range." ); |
60 | |
61 | DNG_ASSERT (fCenter.v >= 0.0 && fCenter.v <= 1.0, |
62 | "Center (vertical) out of range." ); |
63 | |
64 | } |
65 | |
66 | /*****************************************************************************/ |
67 | |
68 | dng_warp_params::~dng_warp_params () |
69 | { |
70 | |
71 | } |
72 | |
73 | /*****************************************************************************/ |
74 | |
75 | bool dng_warp_params::IsNOPAll () const |
76 | { |
77 | |
78 | return IsRadNOPAll () && |
79 | IsTanNOPAll (); |
80 | |
81 | } |
82 | |
83 | /*****************************************************************************/ |
84 | |
85 | bool dng_warp_params::IsNOP (uint32 plane) const |
86 | { |
87 | |
88 | return IsRadNOP (plane) && |
89 | IsTanNOP (plane); |
90 | |
91 | } |
92 | |
93 | /*****************************************************************************/ |
94 | |
95 | bool dng_warp_params::IsRadNOPAll () const |
96 | { |
97 | |
98 | for (uint32 plane = 0; plane < fPlanes; plane++) |
99 | { |
100 | |
101 | if (!IsRadNOP (plane)) |
102 | { |
103 | return false; |
104 | } |
105 | |
106 | } |
107 | |
108 | return true; |
109 | |
110 | } |
111 | |
112 | /*****************************************************************************/ |
113 | |
114 | bool dng_warp_params::IsRadNOP (uint32 /* plane */) const |
115 | { |
116 | |
117 | return false; |
118 | |
119 | } |
120 | |
121 | /*****************************************************************************/ |
122 | |
123 | bool dng_warp_params::IsTanNOPAll () const |
124 | { |
125 | |
126 | for (uint32 plane = 0; plane < fPlanes; plane++) |
127 | { |
128 | |
129 | if (!IsTanNOP (plane)) |
130 | { |
131 | return false; |
132 | } |
133 | |
134 | } |
135 | |
136 | return true; |
137 | |
138 | } |
139 | |
140 | /*****************************************************************************/ |
141 | |
142 | bool dng_warp_params::IsTanNOP (uint32 /* plane */) const |
143 | { |
144 | |
145 | return false; |
146 | |
147 | } |
148 | |
149 | /*****************************************************************************/ |
150 | |
151 | bool dng_warp_params::IsValid () const |
152 | { |
153 | |
154 | if (fPlanes < 1 || fPlanes > kMaxColorPlanes) |
155 | { |
156 | |
157 | return false; |
158 | |
159 | } |
160 | |
161 | if (fCenter.h < 0.0 || |
162 | fCenter.h > 1.0 || |
163 | fCenter.v < 0.0 || |
164 | fCenter.v > 1.0) |
165 | { |
166 | |
167 | return false; |
168 | |
169 | } |
170 | |
171 | return true; |
172 | |
173 | } |
174 | |
175 | /*****************************************************************************/ |
176 | |
177 | bool dng_warp_params::IsValidForNegative (const dng_negative &negative) const |
178 | { |
179 | |
180 | if (!IsValid ()) |
181 | { |
182 | return false; |
183 | } |
184 | |
185 | if ((fPlanes != 1) && |
186 | (fPlanes != negative.ColorChannels ())) |
187 | { |
188 | return false; |
189 | } |
190 | |
191 | return true; |
192 | |
193 | } |
194 | |
195 | /*****************************************************************************/ |
196 | |
197 | real64 dng_warp_params::EvaluateInverse (uint32 plane, |
198 | real64 y) const |
199 | { |
200 | |
201 | const uint32 kMaxIterations = 30; |
202 | const real64 kNearZero = 1.0e-10; |
203 | |
204 | real64 x0 = 0.0; |
205 | real64 y0 = Evaluate (plane, |
206 | x0); |
207 | |
208 | real64 x1 = 1.0; |
209 | real64 y1 = Evaluate (plane, |
210 | x1); |
211 | |
212 | for (uint32 iteration = 0; iteration < kMaxIterations; iteration++) |
213 | { |
214 | |
215 | if (Abs_real64 (y1 - y0) < kNearZero) |
216 | { |
217 | break; |
218 | } |
219 | |
220 | const real64 x2 = Pin_real64 (0.0, |
221 | x1 + (y - y1) * (x1 - x0) / (y1 - y0), |
222 | 1.0); |
223 | |
224 | const real64 y2 = Evaluate (plane, |
225 | x2); |
226 | |
227 | x0 = x1; |
228 | y0 = y1; |
229 | |
230 | x1 = x2; |
231 | y1 = y2; |
232 | |
233 | } |
234 | |
235 | return x1; |
236 | |
237 | } |
238 | |
239 | /*****************************************************************************/ |
240 | |
241 | dng_point_real64 dng_warp_params::EvaluateTangential2 (uint32 plane, |
242 | const dng_point_real64 &diff) const |
243 | { |
244 | |
245 | const real64 dvdv = diff.v * diff.v; |
246 | const real64 dhdh = diff.h * diff.h; |
247 | |
248 | const real64 rr = dvdv + dhdh; |
249 | |
250 | dng_point_real64 diffSqr (dvdv, |
251 | dhdh); |
252 | |
253 | return EvaluateTangential (plane, |
254 | rr, |
255 | diff, |
256 | diffSqr); |
257 | |
258 | } |
259 | |
260 | /*****************************************************************************/ |
261 | |
262 | dng_point_real64 dng_warp_params::EvaluateTangential3 (uint32 plane, |
263 | real64 r2, |
264 | const dng_point_real64 &diff) const |
265 | { |
266 | |
267 | dng_point_real64 diffSqr (diff.v * diff.v, |
268 | diff.h * diff.h); |
269 | |
270 | return EvaluateTangential (plane, |
271 | r2, |
272 | diff, |
273 | diffSqr); |
274 | |
275 | } |
276 | |
277 | /*****************************************************************************/ |
278 | |
279 | void dng_warp_params::Dump () const |
280 | { |
281 | |
282 | #if qDNGValidate |
283 | |
284 | printf ("Planes: %u\n" , (unsigned) fPlanes); |
285 | |
286 | printf (" Optical center:\n" |
287 | " h = %.6lf\n" |
288 | " v = %.6lf\n" , |
289 | fCenter.h, |
290 | fCenter.v); |
291 | |
292 | #endif |
293 | |
294 | } |
295 | |
296 | /*****************************************************************************/ |
297 | |
298 | dng_warp_params_rectilinear::dng_warp_params_rectilinear () |
299 | |
300 | : dng_warp_params () |
301 | |
302 | { |
303 | |
304 | for (uint32 plane = 0; plane < kMaxColorPlanes; plane++) |
305 | { |
306 | |
307 | fRadParams [plane] = dng_vector (4); |
308 | fTanParams [plane] = dng_vector (2); |
309 | |
310 | fRadParams [plane][0] = 1.0; |
311 | |
312 | } |
313 | |
314 | } |
315 | |
316 | /*****************************************************************************/ |
317 | |
318 | dng_warp_params_rectilinear::dng_warp_params_rectilinear (uint32 planes, |
319 | const dng_vector radParams [], |
320 | const dng_vector tanParams [], |
321 | const dng_point_real64 ¢er) |
322 | |
323 | : dng_warp_params (planes, |
324 | center) |
325 | |
326 | { |
327 | |
328 | for (uint32 i = 0; i < fPlanes; i++) |
329 | { |
330 | fRadParams [i] = radParams [i]; |
331 | fTanParams [i] = tanParams [i]; |
332 | } |
333 | |
334 | } |
335 | |
336 | /*****************************************************************************/ |
337 | |
338 | dng_warp_params_rectilinear::~dng_warp_params_rectilinear () |
339 | { |
340 | |
341 | } |
342 | |
343 | /*****************************************************************************/ |
344 | |
345 | bool dng_warp_params_rectilinear::IsRadNOP (uint32 plane) const |
346 | { |
347 | |
348 | DNG_ASSERT (plane < fPlanes, "plane out of range." ); |
349 | |
350 | const dng_vector &r = fRadParams [plane]; |
351 | |
352 | return (r [0] == 1.0 && |
353 | r [1] == 0.0 && |
354 | r [2] == 0.0 && |
355 | r [3] == 0.0); |
356 | |
357 | } |
358 | |
359 | /*****************************************************************************/ |
360 | |
361 | bool dng_warp_params_rectilinear::IsTanNOP (uint32 plane) const |
362 | { |
363 | |
364 | DNG_ASSERT (plane < fPlanes, "plane out of range." ); |
365 | |
366 | const dng_vector &t = fTanParams [plane]; |
367 | |
368 | return (t [0] == 0.0 && |
369 | t [1] == 0.0); |
370 | |
371 | } |
372 | |
373 | /*****************************************************************************/ |
374 | |
375 | bool dng_warp_params_rectilinear::IsValid () const |
376 | { |
377 | |
378 | for (uint32 plane = 0; plane < fPlanes; plane++) |
379 | { |
380 | |
381 | if (fRadParams [plane].Count () != 4) |
382 | { |
383 | return false; |
384 | } |
385 | |
386 | if (fTanParams [plane].Count () < 2) |
387 | { |
388 | return false; |
389 | } |
390 | |
391 | } |
392 | |
393 | return dng_warp_params::IsValid (); |
394 | |
395 | } |
396 | |
397 | /*****************************************************************************/ |
398 | |
399 | void dng_warp_params_rectilinear::PropagateToAllPlanes (uint32 totalPlanes) |
400 | { |
401 | |
402 | for (uint32 plane = fPlanes; plane < totalPlanes; plane++) |
403 | { |
404 | |
405 | fRadParams [plane] = fRadParams [0]; |
406 | fTanParams [plane] = fTanParams [0]; |
407 | |
408 | } |
409 | |
410 | } |
411 | |
412 | /*****************************************************************************/ |
413 | |
414 | real64 dng_warp_params_rectilinear::Evaluate (uint32 plane, |
415 | real64 x) const |
416 | { |
417 | |
418 | const dng_vector &K = fRadParams [plane]; // Coefficients. |
419 | |
420 | const real64 x2 = x * x; |
421 | |
422 | return x * (K [0] + x2 * (K [1] + x2 * (K [2] + x2 * K [3]))); |
423 | |
424 | } |
425 | |
426 | /*****************************************************************************/ |
427 | |
428 | real64 dng_warp_params_rectilinear::EvaluateRatio (uint32 plane, |
429 | real64 r2) const |
430 | { |
431 | |
432 | const dng_vector &K = fRadParams [plane]; // Coefficients. |
433 | |
434 | return K [0] + r2 * (K [1] + r2 * (K [2] + r2 * K [3])); |
435 | |
436 | } |
437 | |
438 | /*****************************************************************************/ |
439 | |
440 | dng_point_real64 dng_warp_params_rectilinear::EvaluateTangential (uint32 plane, |
441 | real64 r2, |
442 | const dng_point_real64 &diff, |
443 | const dng_point_real64 &diff2) const |
444 | { |
445 | |
446 | const real64 kt0 = fTanParams [plane][0]; |
447 | const real64 kt1 = fTanParams [plane][1]; |
448 | |
449 | const real64 dh = diff.h; |
450 | const real64 dv = diff.v; |
451 | |
452 | const real64 dhdh = diff2.h; |
453 | const real64 dvdv = diff2.v; |
454 | |
455 | return dng_point_real64 (kt0 * (r2 + 2.0 * dvdv) + (2.0 * kt1 * dh * dv), // v |
456 | kt1 * (r2 + 2.0 * dhdh) + (2.0 * kt0 * dh * dv)); // h |
457 | |
458 | } |
459 | |
460 | /*****************************************************************************/ |
461 | |
462 | real64 dng_warp_params_rectilinear::MaxSrcRadiusGap (real64 maxDstGap) const |
463 | { |
464 | |
465 | real64 maxSrcGap = 0.0; |
466 | |
467 | for (uint32 plane = 0; plane < fPlanes; plane++) |
468 | { |
469 | |
470 | const dng_vector &coefs = fRadParams [plane]; |
471 | |
472 | const real64 k3 = coefs [1]; |
473 | const real64 k5 = coefs [2]; |
474 | const real64 k7 = coefs [3]; |
475 | |
476 | // |
477 | // Let f (r) be the radius warp function. Consider the function |
478 | // |
479 | // gap (r) = f (r + maxDstGap) - f (r) |
480 | // |
481 | // We wish to maximize gap (r) over the domain [0, 1 - maxDstGap]. This will |
482 | // give us a reasonable upper bound on the src tile size, independent of |
483 | // dstBounds. |
484 | // |
485 | // As usual, we maximize gap (r) by examining its critical points, i.e., by |
486 | // considering the roots of its derivative which lie in the domain [0, 1 - |
487 | // maxDstGap]. gap' (r) is a 5th-order polynomial. One of its roots is |
488 | // -maxDstGap / 2, which is negative and hence lies outside the domain of |
489 | // interest. This leaves 4 other possible roots. We solve for these |
490 | // analytically below. |
491 | // |
492 | |
493 | real64 roots [4]; |
494 | uint32 numRoots = 0; |
495 | |
496 | if (k7 == 0.0) |
497 | { |
498 | |
499 | if (k5 == 0.0) |
500 | { |
501 | |
502 | // No roots in [0,1]. |
503 | |
504 | } |
505 | |
506 | else |
507 | { |
508 | |
509 | // k7 is zero, but k5 is non-zero. At most two real roots. |
510 | |
511 | const real64 discrim = 25.0 * (-6.0 * k3 * k5 - 5.0 * k5 * maxDstGap * maxDstGap); |
512 | |
513 | if (discrim >= 0.0) |
514 | { |
515 | |
516 | // Two real roots. |
517 | |
518 | const real64 scale = 0.1 * k5; |
519 | const real64 offset = -5.0 * maxDstGap * k5; |
520 | const real64 sDiscrim = sqrt (discrim); |
521 | |
522 | roots [numRoots++] = scale * (offset + sDiscrim); |
523 | roots [numRoots++] = scale * (offset - sDiscrim); |
524 | |
525 | } |
526 | |
527 | } |
528 | |
529 | } |
530 | |
531 | else |
532 | { |
533 | |
534 | // k7 is non-zero. Up to 4 real roots. |
535 | |
536 | const real64 d = maxDstGap; |
537 | const real64 d2 = d * d; |
538 | const real64 d4 = d2 * d2; |
539 | |
540 | const real64 discrim = 25.0 * k5 * k5 |
541 | - 63.0 * k3 * k7 |
542 | + 35.0 * d2 * k5 * k7 |
543 | + 49.0 * d4 * k7 * k7; |
544 | |
545 | if (discrim >= 0.0) |
546 | { |
547 | |
548 | const real64 sDiscrim = 4.0 * k7 * sqrt (discrim); |
549 | |
550 | const real64 offset = -20.0 * k5 * k7 - 35.0 * d2 * k7 * k7; |
551 | |
552 | const real64 discrim1 = offset - sDiscrim; |
553 | const real64 discrim2 = offset + sDiscrim; |
554 | |
555 | const real64 scale = sqrt (21.0) / (42.0 * k7); |
556 | |
557 | if (discrim1 >= 0.0) |
558 | { |
559 | |
560 | const real64 offset1 = -d * 0.5; |
561 | const real64 sDiscrim1 = scale * sqrt (discrim1); |
562 | |
563 | roots [numRoots++] = offset1 + sDiscrim1; |
564 | roots [numRoots++] = offset1 - sDiscrim1; |
565 | |
566 | } |
567 | |
568 | if (discrim2 >= 0.0) |
569 | { |
570 | |
571 | const real64 offset2 = -d * 0.5; |
572 | const real64 sDiscrim2 = scale * sqrt (discrim2); |
573 | |
574 | roots [numRoots++] = offset2 + sDiscrim2; |
575 | roots [numRoots++] = offset2 - sDiscrim2; |
576 | |
577 | } |
578 | |
579 | } |
580 | |
581 | } |
582 | |
583 | real64 planeMaxSrcGap = 0.0; |
584 | |
585 | // Examine the endpoints. |
586 | |
587 | { |
588 | |
589 | // Check left endpoint: f (maxDstGap) - f (0). Remember that f (0) == 0. |
590 | |
591 | const real64 gap1 = Evaluate (plane, maxDstGap); |
592 | |
593 | planeMaxSrcGap = Max_real64 (planeMaxSrcGap, gap1); |
594 | |
595 | // Check right endpoint: f (1) - f (1 - maxDstGap). |
596 | |
597 | const real64 gap2 = Evaluate (plane, 1.0) |
598 | - Evaluate (plane, 1.0 - maxDstGap); |
599 | |
600 | planeMaxSrcGap = Max_real64 (planeMaxSrcGap, gap2); |
601 | |
602 | } |
603 | |
604 | // Examine the roots we found, if any. |
605 | |
606 | for (uint32 i = 0; i < numRoots; i++) |
607 | { |
608 | |
609 | const real64 r = roots [i]; |
610 | |
611 | if (r > 0.0 && r < 1.0 - maxDstGap) |
612 | { |
613 | |
614 | const real64 gap = Evaluate (plane, r + maxDstGap) |
615 | - Evaluate (plane, r); |
616 | |
617 | planeMaxSrcGap = Max_real64 (planeMaxSrcGap, gap); |
618 | |
619 | } |
620 | |
621 | } |
622 | |
623 | maxSrcGap = Max_real64 (maxSrcGap, |
624 | planeMaxSrcGap); |
625 | |
626 | } |
627 | |
628 | return maxSrcGap; |
629 | |
630 | } |
631 | |
632 | /*****************************************************************************/ |
633 | |
634 | dng_point_real64 dng_warp_params_rectilinear::MaxSrcTanGap (dng_point_real64 minDst, |
635 | dng_point_real64 maxDst) const |
636 | { |
637 | |
638 | const real64 v [] = { minDst.v, maxDst.v, 0.0 }; |
639 | const real64 h [] = { minDst.h, maxDst.h, 0.0 }; |
640 | |
641 | dng_point_real64 maxGap; |
642 | |
643 | for (uint32 plane = 0; plane < fPlanes; plane++) |
644 | { |
645 | |
646 | real64 hMin = +FLT_MAX; |
647 | real64 hMax = -FLT_MAX; |
648 | |
649 | real64 vMin = +FLT_MAX; |
650 | real64 vMax = -FLT_MAX; |
651 | |
652 | for (uint32 i = 0; i < 3; i++) |
653 | { |
654 | |
655 | for (uint32 j = 0; j < 3; j++) |
656 | { |
657 | |
658 | dng_point_real64 dstDiff (v [i], |
659 | h [j]); |
660 | |
661 | dng_point_real64 srcDiff = EvaluateTangential2 (plane, |
662 | dstDiff); |
663 | |
664 | hMin = Min_real64 (hMin, srcDiff.h); |
665 | hMax = Max_real64 (hMax, srcDiff.h); |
666 | |
667 | vMin = Min_real64 (vMin, srcDiff.v); |
668 | vMax = Max_real64 (vMax, srcDiff.v); |
669 | |
670 | } |
671 | |
672 | } |
673 | |
674 | const real64 hGap = hMax - hMin; |
675 | const real64 vGap = vMax - vMin; |
676 | |
677 | maxGap.h = Max_real64 (maxGap.h, hGap); |
678 | maxGap.v = Max_real64 (maxGap.v, vGap); |
679 | |
680 | } |
681 | |
682 | return maxGap; |
683 | |
684 | } |
685 | |
686 | /*****************************************************************************/ |
687 | |
688 | void dng_warp_params_rectilinear::Dump () const |
689 | { |
690 | |
691 | #if qDNGValidate |
692 | |
693 | dng_warp_params::Dump (); |
694 | |
695 | for (uint32 plane = 0; plane < fPlanes; plane++) |
696 | { |
697 | |
698 | printf (" Plane %u:\n" , (unsigned) plane); |
699 | |
700 | printf (" Radial params: %.6lf, %.6lf, %.6lf, %.6lf\n" , |
701 | fRadParams [plane][0], |
702 | fRadParams [plane][1], |
703 | fRadParams [plane][2], |
704 | fRadParams [plane][3]); |
705 | |
706 | printf (" Tangential params: %.6lf, %.6lf\n" , |
707 | fTanParams [plane][0], |
708 | fTanParams [plane][1]); |
709 | |
710 | } |
711 | |
712 | #endif |
713 | |
714 | } |
715 | |
716 | /*****************************************************************************/ |
717 | |
718 | dng_warp_params_fisheye::dng_warp_params_fisheye () |
719 | |
720 | : dng_warp_params () |
721 | |
722 | { |
723 | |
724 | for (uint32 plane = 0; plane < kMaxColorPlanes; plane++) |
725 | { |
726 | |
727 | fRadParams [plane] = dng_vector (4); |
728 | |
729 | } |
730 | |
731 | } |
732 | |
733 | /*****************************************************************************/ |
734 | |
735 | dng_warp_params_fisheye::dng_warp_params_fisheye (uint32 planes, |
736 | const dng_vector radParams [], |
737 | const dng_point_real64 ¢er) |
738 | |
739 | : dng_warp_params (planes, center) |
740 | |
741 | { |
742 | |
743 | for (uint32 i = 0; i < fPlanes; i++) |
744 | { |
745 | |
746 | fRadParams [i] = radParams [i]; |
747 | |
748 | } |
749 | |
750 | } |
751 | |
752 | /*****************************************************************************/ |
753 | |
754 | dng_warp_params_fisheye::~dng_warp_params_fisheye () |
755 | { |
756 | |
757 | } |
758 | |
759 | /*****************************************************************************/ |
760 | |
761 | bool dng_warp_params_fisheye::IsRadNOP (uint32 /* plane */) const |
762 | { |
763 | |
764 | return false; |
765 | |
766 | } |
767 | |
768 | /*****************************************************************************/ |
769 | |
770 | bool dng_warp_params_fisheye::IsTanNOP (uint32 /* plane */) const |
771 | { |
772 | |
773 | return true; |
774 | |
775 | } |
776 | |
777 | /*****************************************************************************/ |
778 | |
779 | bool dng_warp_params_fisheye::IsValid () const |
780 | { |
781 | |
782 | for (uint32 plane = 0; plane < fPlanes; plane++) |
783 | { |
784 | |
785 | if (fRadParams [plane].Count () != 4) |
786 | { |
787 | return false; |
788 | } |
789 | |
790 | } |
791 | |
792 | return dng_warp_params::IsValid (); |
793 | |
794 | } |
795 | |
796 | /*****************************************************************************/ |
797 | |
798 | void dng_warp_params_fisheye::PropagateToAllPlanes (uint32 totalPlanes) |
799 | { |
800 | |
801 | for (uint32 plane = fPlanes; plane < totalPlanes; plane++) |
802 | { |
803 | |
804 | fRadParams [plane] = fRadParams [0]; |
805 | |
806 | } |
807 | |
808 | } |
809 | |
810 | /*****************************************************************************/ |
811 | |
812 | real64 dng_warp_params_fisheye::Evaluate (uint32 plane, |
813 | real64 r) const |
814 | { |
815 | |
816 | const real64 t = atan (r); |
817 | |
818 | const dng_vector &K = fRadParams [plane]; |
819 | |
820 | const real64 t2 = t * t; |
821 | |
822 | return t * (K [0] + t2 * (K [1] + t2 * (K [2] + t2 * K [3]))); |
823 | |
824 | } |
825 | |
826 | /*****************************************************************************/ |
827 | |
828 | real64 dng_warp_params_fisheye::EvaluateRatio (uint32 plane, |
829 | real64 rSqr) const |
830 | { |
831 | |
832 | const real64 eps = 1.0e-12; |
833 | |
834 | if (rSqr < eps) |
835 | { |
836 | |
837 | // r is very close to zero. |
838 | |
839 | return 1.0; |
840 | |
841 | } |
842 | |
843 | const real64 r = sqrt (rSqr); |
844 | |
845 | return Evaluate (plane, r) / r; |
846 | |
847 | } |
848 | |
849 | /*****************************************************************************/ |
850 | |
851 | dng_point_real64 dng_warp_params_fisheye::EvaluateTangential (uint32 /* plane */, |
852 | real64 /* r2 */, |
853 | const dng_point_real64 & /* diff */, |
854 | const dng_point_real64 & /* diff2 */) const |
855 | { |
856 | |
857 | // This fisheye model does not support tangential warping. |
858 | |
859 | ThrowProgramError (); |
860 | |
861 | return dng_point_real64 (0.0, 0.0); |
862 | |
863 | } |
864 | |
865 | /*****************************************************************************/ |
866 | |
867 | real64 dng_warp_params_fisheye::MaxSrcRadiusGap (real64 maxDstGap) const |
868 | { |
869 | |
870 | // |
871 | // Let f (r) be the radius warp function. Consider the function |
872 | // |
873 | // gap (r) = f (r + maxDstGap) - f (r) |
874 | // |
875 | // We wish to maximize gap (r) over the domain [0, 1 - maxDstGap]. This will |
876 | // give us a reasonable upper bound on the src tile size, independent of |
877 | // dstBounds. |
878 | // |
879 | // Ideally, we'd like to maximize gap (r) by examining its critical points, |
880 | // i.e., by considering the roots of its derivative which lie in the domain [0, |
881 | // 1 - maxDstGap]. However, gap' (r) is too complex to find its roots |
882 | // analytically. |
883 | // |
884 | |
885 | real64 maxSrcGap = 0.0; |
886 | |
887 | DNG_REQUIRE (maxDstGap > 0.0, "maxDstGap must be positive." ); |
888 | |
889 | const real64 kMaxValue = 1.0 - maxDstGap; |
890 | |
891 | const uint32 kSteps = 128; |
892 | |
893 | const real64 kNorm = kMaxValue / real64 (kSteps - 1); |
894 | |
895 | for (uint32 plane = 0; plane < fPlanes; plane++) |
896 | { |
897 | |
898 | for (uint32 i = 0; i < kSteps; i++) |
899 | { |
900 | |
901 | const real64 tt = i * kNorm; |
902 | |
903 | const real64 gap = Evaluate (plane, tt + maxDstGap) |
904 | - Evaluate (plane, tt); |
905 | |
906 | maxSrcGap = Max_real64 (maxSrcGap, |
907 | gap); |
908 | |
909 | } |
910 | |
911 | } |
912 | |
913 | return maxSrcGap; |
914 | |
915 | } |
916 | |
917 | /*****************************************************************************/ |
918 | |
919 | dng_point_real64 dng_warp_params_fisheye::MaxSrcTanGap (dng_point_real64 /* minDst */, |
920 | dng_point_real64 /* maxDst */) const |
921 | { |
922 | |
923 | // This fisheye model does not support tangential distortion. |
924 | |
925 | return dng_point_real64 (0.0, 0.0); |
926 | |
927 | } |
928 | |
929 | /*****************************************************************************/ |
930 | |
931 | void dng_warp_params_fisheye::Dump () const |
932 | { |
933 | |
934 | #if qDNGValidate |
935 | |
936 | dng_warp_params::Dump (); |
937 | |
938 | for (uint32 plane = 0; plane < fPlanes; plane++) |
939 | { |
940 | |
941 | printf (" Plane %u:\n" , (unsigned) plane); |
942 | |
943 | printf (" Radial params: %.6lf, %.6lf, %.6lf, %.6lf\n" , |
944 | fRadParams [plane][0], |
945 | fRadParams [plane][1], |
946 | fRadParams [plane][2], |
947 | fRadParams [plane][3]); |
948 | |
949 | } |
950 | |
951 | #endif |
952 | |
953 | } |
954 | |
955 | /*****************************************************************************/ |
956 | |
957 | class dng_filter_warp: public dng_filter_task |
958 | { |
959 | |
960 | protected: |
961 | |
962 | AutoPtr<dng_warp_params> fParams; |
963 | |
964 | dng_point_real64 fCenter; |
965 | |
966 | dng_resample_weights_2d fWeights; |
967 | |
968 | real64 fNormRadius; |
969 | real64 fInvNormRadius; |
970 | |
971 | bool fIsRadNOP; |
972 | bool fIsTanNOP; |
973 | |
974 | const real64 fPixelScaleV; |
975 | const real64 fPixelScaleVInv; |
976 | |
977 | public: |
978 | |
979 | dng_filter_warp (const dng_image &srcImage, |
980 | dng_image &dstImage, |
981 | const dng_negative &negative, |
982 | AutoPtr<dng_warp_params> ¶ms); |
983 | |
984 | virtual void Initialize (dng_host &host); |
985 | |
986 | virtual dng_rect SrcArea (const dng_rect &dstArea); |
987 | |
988 | virtual dng_point SrcTileSize (const dng_point &dstTileSize); |
989 | |
990 | virtual void ProcessArea (uint32 threadIndex, |
991 | dng_pixel_buffer &srcBuffer, |
992 | dng_pixel_buffer &dstBuffer); |
993 | |
994 | virtual dng_point_real64 GetSrcPixelPosition (const dng_point_real64 &dst, |
995 | uint32 plane); |
996 | |
997 | }; |
998 | |
999 | /*****************************************************************************/ |
1000 | |
1001 | dng_filter_warp::dng_filter_warp (const dng_image &srcImage, |
1002 | dng_image &dstImage, |
1003 | const dng_negative &negative, |
1004 | AutoPtr<dng_warp_params> ¶ms) |
1005 | |
1006 | : dng_filter_task (srcImage, |
1007 | dstImage) |
1008 | |
1009 | , fParams (params.Release ()) |
1010 | |
1011 | , fCenter () |
1012 | |
1013 | , fWeights () |
1014 | |
1015 | , fNormRadius (1.0) |
1016 | , fInvNormRadius (1.0) |
1017 | |
1018 | , fIsRadNOP (false) |
1019 | , fIsTanNOP (false) |
1020 | |
1021 | , fPixelScaleV (1.0 / negative.PixelAspectRatio ()) |
1022 | , fPixelScaleVInv (1.0 / fPixelScaleV) |
1023 | |
1024 | { |
1025 | |
1026 | // Force float processing. |
1027 | |
1028 | fSrcPixelType = ttFloat; |
1029 | fDstPixelType = ttFloat; |
1030 | |
1031 | fIsRadNOP = fParams->IsRadNOPAll (); |
1032 | fIsTanNOP = fParams->IsTanNOPAll (); |
1033 | |
1034 | const uint32 negPlanes = negative.ColorChannels (); |
1035 | |
1036 | DNG_ASSERT (negPlanes >= 1, "Too few planes." ); |
1037 | DNG_ASSERT (negPlanes <= kMaxColorPlanes, "Too many planes." ); |
1038 | |
1039 | (void) negPlanes; |
1040 | |
1041 | // At least one set of params must do something interesting. |
1042 | |
1043 | if (fIsRadNOP && fIsTanNOP) |
1044 | { |
1045 | ThrowProgramError (); |
1046 | } |
1047 | |
1048 | // Make sure the warp params are valid for this negative. |
1049 | |
1050 | if (!fParams->IsValidForNegative (negative)) |
1051 | { |
1052 | ThrowBadFormat (); |
1053 | } |
1054 | |
1055 | // Compute center. |
1056 | |
1057 | const dng_rect bounds = srcImage.Bounds (); |
1058 | |
1059 | fCenter.h = Lerp_real64 ((real64) bounds.l, |
1060 | (real64) bounds.r, |
1061 | fParams->fCenter.h); |
1062 | |
1063 | fCenter.v = Lerp_real64 ((real64) bounds.t, |
1064 | (real64) bounds.b, |
1065 | fParams->fCenter.v); |
1066 | |
1067 | // Compute max pixel radius and derive various normalized radius values. Note |
1068 | // that when computing the max pixel radius, we must take into account the pixel |
1069 | // aspect ratio. |
1070 | |
1071 | { |
1072 | |
1073 | dng_rect squareBounds (bounds); |
1074 | |
1075 | squareBounds.b = squareBounds.t + |
1076 | Round_int32 (fPixelScaleV * (real64) squareBounds.H ()); |
1077 | |
1078 | const dng_point_real64 squareCenter (Lerp_real64 ((real64) squareBounds.t, |
1079 | (real64) squareBounds.b, |
1080 | fParams->fCenter.v), |
1081 | |
1082 | Lerp_real64 ((real64) squareBounds.l, |
1083 | (real64) squareBounds.r, |
1084 | fParams->fCenter.h)); |
1085 | |
1086 | fNormRadius = MaxDistancePointToRect (squareCenter, |
1087 | squareBounds); |
1088 | |
1089 | fInvNormRadius = 1.0 / fNormRadius; |
1090 | |
1091 | } |
1092 | |
1093 | // Propagate warp params to other planes. |
1094 | |
1095 | fParams->PropagateToAllPlanes (fDstPlanes); |
1096 | |
1097 | } |
1098 | |
1099 | /*****************************************************************************/ |
1100 | |
1101 | void dng_filter_warp::Initialize (dng_host &host) |
1102 | { |
1103 | |
1104 | // Make resample weights. |
1105 | |
1106 | const dng_resample_function &kernel = dng_resample_bicubic::Get (); |
1107 | |
1108 | fWeights.Initialize (kernel, |
1109 | host.Allocator ()); |
1110 | |
1111 | } |
1112 | |
1113 | /*****************************************************************************/ |
1114 | |
1115 | dng_rect dng_filter_warp::SrcArea (const dng_rect &dstArea) |
1116 | { |
1117 | |
1118 | // Walk each pixel of the boundary of dstArea, map it to the uncorrected src |
1119 | // pixel position, and return the rectangle that contains all such src pixels. |
1120 | |
1121 | int32 xMin = INT_MAX; |
1122 | int32 xMax = INT_MIN; |
1123 | int32 yMin = INT_MAX; |
1124 | int32 yMax = INT_MIN; |
1125 | |
1126 | for (uint32 plane = 0; plane < fDstPlanes; plane++) |
1127 | { |
1128 | |
1129 | // Top and bottom edges. |
1130 | |
1131 | for (int32 c = dstArea.l; c < dstArea.r; c++) |
1132 | { |
1133 | |
1134 | // Top edge. |
1135 | |
1136 | { |
1137 | |
1138 | const dng_point_real64 dst (dstArea.t, c); |
1139 | |
1140 | const dng_point_real64 src = GetSrcPixelPosition (dst, plane); |
1141 | |
1142 | const int32 y = ConvertDoubleToInt32 (floor (src.v)); |
1143 | |
1144 | yMin = Min_int32 (yMin, y); |
1145 | |
1146 | } |
1147 | |
1148 | // Bottom edge. |
1149 | |
1150 | { |
1151 | |
1152 | const dng_point_real64 dst (dstArea.b - 1, c); |
1153 | |
1154 | const dng_point_real64 src = GetSrcPixelPosition (dst, plane); |
1155 | |
1156 | const int32 y = ConvertDoubleToInt32 (ceil (src.v)); |
1157 | |
1158 | yMax = Max_int32 (yMax, y); |
1159 | |
1160 | } |
1161 | |
1162 | } |
1163 | |
1164 | // Left and right edges. |
1165 | |
1166 | for (int32 r = dstArea.t; r < dstArea.b; r++) |
1167 | { |
1168 | |
1169 | // Left edge. |
1170 | |
1171 | { |
1172 | |
1173 | const dng_point_real64 dst (r, dstArea.l); |
1174 | |
1175 | const dng_point_real64 src = GetSrcPixelPosition (dst, plane); |
1176 | |
1177 | const int32 x = ConvertDoubleToInt32 (floor (src.h)); |
1178 | |
1179 | xMin = Min_int32 (xMin, x); |
1180 | |
1181 | } |
1182 | |
1183 | // Right edge. |
1184 | |
1185 | { |
1186 | |
1187 | const dng_point_real64 dst (r, dstArea.r - 1); |
1188 | |
1189 | const dng_point_real64 src = GetSrcPixelPosition (dst, plane); |
1190 | |
1191 | const int32 x = ConvertDoubleToInt32 (ceil (src.h)); |
1192 | |
1193 | xMax = Max_int32 (xMax, x); |
1194 | |
1195 | } |
1196 | |
1197 | } |
1198 | |
1199 | } |
1200 | |
1201 | // Pad each side by filter radius. |
1202 | |
1203 | const int32 pad = ConvertUint32ToInt32(fWeights.Radius()); |
1204 | |
1205 | xMin = SafeInt32Sub(xMin, pad); |
1206 | yMin = SafeInt32Sub(yMin, pad); |
1207 | xMax = SafeInt32Add(xMax, pad); |
1208 | yMax = SafeInt32Add(yMax, pad); |
1209 | |
1210 | xMax = SafeInt32Add(xMax, 1); |
1211 | yMax = SafeInt32Add(yMax, 1); |
1212 | |
1213 | const dng_rect srcArea (yMin, |
1214 | xMin, |
1215 | yMax, |
1216 | xMax); |
1217 | |
1218 | return srcArea; |
1219 | |
1220 | } |
1221 | |
1222 | /*****************************************************************************/ |
1223 | |
1224 | dng_point dng_filter_warp::SrcTileSize (const dng_point &dstTileSize) |
1225 | { |
1226 | |
1227 | // Obtain an upper bound on the source tile size. We'll do this by considering |
1228 | // upper bounds on the radial and tangential warp components separately, then add |
1229 | // them together. This is not a tight bound but is good enough because the |
1230 | // tangential terms are usually quite small. |
1231 | |
1232 | // Get upper bound on src tile size from radial warp. |
1233 | |
1234 | DNG_REQUIRE (dstTileSize.v > 0, "Invalid tile height." ); |
1235 | DNG_REQUIRE (dstTileSize.h > 0, "Invalid tile width." ); |
1236 | |
1237 | const real64 maxDstGap = fInvNormRadius * hypot ((real64) dstTileSize.h, |
1238 | (real64) dstTileSize.v); |
1239 | |
1240 | dng_point srcTileSize; |
1241 | |
1242 | if (maxDstGap >= 1.0) |
1243 | { |
1244 | |
1245 | // The proposed tile size is unusually large, i.e., its hypotenuse is larger |
1246 | // than the maximum radius. Bite the bullet and just return a tile size big |
1247 | // enough to process the whole image. |
1248 | |
1249 | srcTileSize = SrcArea (fDstImage.Bounds ()).Size (); |
1250 | |
1251 | } |
1252 | |
1253 | else |
1254 | { |
1255 | |
1256 | // maxDstGap < 1.0. |
1257 | |
1258 | const real64 maxSrcGap = fParams->MaxSrcRadiusGap (maxDstGap); |
1259 | |
1260 | const int32 dim = ConvertDoubleToInt32 (ceil (maxSrcGap * fNormRadius)); |
1261 | |
1262 | srcTileSize = dng_point (dim, dim); |
1263 | |
1264 | } |
1265 | |
1266 | srcTileSize.h += ConvertUint32ToInt32(fWeights.Width()); |
1267 | srcTileSize.v += ConvertUint32ToInt32(fWeights.Width()); |
1268 | |
1269 | // Get upper bound on src tile size from tangential warp. |
1270 | |
1271 | const dng_rect_real64 bounds (fSrcImage.Bounds ()); |
1272 | |
1273 | const dng_point_real64 minDst ((bounds.t - fCenter.v) * fInvNormRadius, |
1274 | (bounds.l - fCenter.h) * fInvNormRadius); |
1275 | |
1276 | const dng_point_real64 maxDst ((bounds.b - 1.0 - fCenter.v) * fInvNormRadius, |
1277 | (bounds.r - 1.0 - fCenter.h) * fInvNormRadius); |
1278 | |
1279 | const dng_point_real64 srcTanGap = fParams->MaxSrcTanGap (minDst, |
1280 | maxDst); |
1281 | |
1282 | // Add the two bounds together. |
1283 | |
1284 | srcTileSize.v += ConvertDoubleToInt32 (ceil (srcTanGap.v * fNormRadius)); |
1285 | srcTileSize.h += ConvertDoubleToInt32 (ceil (srcTanGap.h * fNormRadius)); |
1286 | |
1287 | return srcTileSize; |
1288 | |
1289 | } |
1290 | |
1291 | /*****************************************************************************/ |
1292 | |
1293 | void dng_filter_warp::ProcessArea (uint32 /* threadIndex */, |
1294 | dng_pixel_buffer &srcBuffer, |
1295 | dng_pixel_buffer &dstBuffer) |
1296 | { |
1297 | |
1298 | // Prepare resample constants. |
1299 | |
1300 | const int32 wCount = fWeights.Width (); |
1301 | |
1302 | const dng_point srcOffset (fWeights.Offset (), |
1303 | fWeights.Offset ()); |
1304 | |
1305 | const real64 numSubsamples = (real64) kResampleSubsampleCount2D; |
1306 | |
1307 | // Prepare area and step constants. |
1308 | |
1309 | const dng_rect srcArea = srcBuffer.fArea; |
1310 | const dng_rect dstArea = dstBuffer.fArea; |
1311 | |
1312 | const int32 srcRowStep = (int32) srcBuffer.RowStep (); |
1313 | |
1314 | const int32 hMin = srcArea.l; |
1315 | const int32 hMax = SafeInt32Sub (SafeInt32Sub (srcArea.r, wCount), 1); |
1316 | |
1317 | const int32 vMin = srcArea.t; |
1318 | const int32 vMax = SafeInt32Sub (SafeInt32Sub (srcArea.b, wCount), 1); |
1319 | |
1320 | if (hMax < hMin || vMax < vMin) |
1321 | { |
1322 | |
1323 | ThrowBadFormat ("Empty source area in dng_filter_warp." ); |
1324 | |
1325 | } |
1326 | |
1327 | // Warp each plane. |
1328 | |
1329 | for (uint32 plane = 0; plane < dstBuffer.fPlanes; plane++) |
1330 | { |
1331 | |
1332 | real32 *dPtr = dstBuffer.DirtyPixel_real32 (dstArea.t, |
1333 | dstArea.l, |
1334 | plane); |
1335 | |
1336 | for (int32 dstRow = dstArea.t; dstRow < dstArea.b; dstRow++) |
1337 | { |
1338 | |
1339 | uint32 dstIndex = 0; |
1340 | |
1341 | for (int32 dstCol = dstArea.l; dstCol < dstArea.r; dstCol++, dstIndex++) |
1342 | { |
1343 | |
1344 | // Get destination (corrected) pixel position. |
1345 | |
1346 | const dng_point_real64 dPos ((real64) dstRow, |
1347 | (real64) dstCol); |
1348 | |
1349 | // Warp to source (uncorrected) pixel position. |
1350 | |
1351 | const dng_point_real64 sPos = GetSrcPixelPosition (dPos, |
1352 | plane); |
1353 | |
1354 | // Decompose into integer and fractional parts. |
1355 | |
1356 | dng_point sInt (ConvertDoubleToInt32 (floor (sPos.v)), |
1357 | ConvertDoubleToInt32 (floor (sPos.h))); |
1358 | |
1359 | dng_point sFct (ConvertDoubleToInt32 ((sPos.v - (real64) sInt.v) * numSubsamples), |
1360 | ConvertDoubleToInt32 ((sPos.h - (real64) sInt.h) * numSubsamples)); |
1361 | |
1362 | // Add resample offset. |
1363 | |
1364 | sInt = sInt + srcOffset; |
1365 | |
1366 | // Clip. |
1367 | |
1368 | if (sInt.h < hMin) |
1369 | { |
1370 | sInt.h = hMin; |
1371 | sFct.h = 0; |
1372 | } |
1373 | |
1374 | else if (sInt.h > hMax) |
1375 | { |
1376 | sInt.h = hMax; |
1377 | sFct.h = 0; |
1378 | } |
1379 | |
1380 | if (sInt.v < vMin) |
1381 | { |
1382 | sInt.v = vMin; |
1383 | sFct.v = 0; |
1384 | } |
1385 | |
1386 | else if (sInt.v > vMax) |
1387 | { |
1388 | sInt.v = vMax; |
1389 | sFct.v = 0; |
1390 | } |
1391 | |
1392 | // Perform 2D resample. |
1393 | |
1394 | const real32 *w = fWeights.Weights32 (sFct); |
1395 | |
1396 | const real32 *s = srcBuffer.ConstPixel_real32 (sInt.v, |
1397 | sInt.h, |
1398 | plane); |
1399 | |
1400 | real32 total = 0.0f; |
1401 | |
1402 | for (int32 i = 0; i < wCount; i++) |
1403 | { |
1404 | |
1405 | for (int32 j = 0; j < wCount; j++) |
1406 | { |
1407 | |
1408 | total += w [j] * s [j]; |
1409 | |
1410 | } |
1411 | |
1412 | w += wCount; |
1413 | s += srcRowStep; |
1414 | |
1415 | } |
1416 | |
1417 | // Store final pixel value. |
1418 | |
1419 | dPtr [dstIndex] = Pin_real32 (total); |
1420 | |
1421 | } |
1422 | |
1423 | // Advance to next row. |
1424 | |
1425 | dPtr += dstBuffer.RowStep (); |
1426 | |
1427 | } |
1428 | |
1429 | } |
1430 | |
1431 | } |
1432 | |
1433 | /*****************************************************************************/ |
1434 | |
1435 | dng_point_real64 dng_filter_warp::GetSrcPixelPosition (const dng_point_real64 &dst, |
1436 | uint32 plane) |
1437 | { |
1438 | |
1439 | const dng_point_real64 diff = dst - fCenter; |
1440 | |
1441 | const dng_point_real64 diffNorm (diff.v * fInvNormRadius, |
1442 | diff.h * fInvNormRadius); |
1443 | |
1444 | const dng_point_real64 diffNormScaled (diffNorm.v * fPixelScaleV, |
1445 | diffNorm.h); |
1446 | |
1447 | const dng_point_real64 diffNormSqr (diffNormScaled.v * diffNormScaled.v, |
1448 | diffNormScaled.h * diffNormScaled.h); |
1449 | |
1450 | const real64 rr = Min_real64 (diffNormSqr.v + diffNormSqr.h, |
1451 | 1.0); |
1452 | |
1453 | dng_point_real64 dSrc; |
1454 | |
1455 | if (fIsTanNOP) |
1456 | { |
1457 | |
1458 | // Radial only. |
1459 | |
1460 | const real64 ratio = fParams->EvaluateRatio (plane, |
1461 | rr); |
1462 | |
1463 | dSrc.h = diff.h * ratio; |
1464 | dSrc.v = diff.v * ratio; |
1465 | |
1466 | } |
1467 | |
1468 | else if (fIsRadNOP) |
1469 | { |
1470 | |
1471 | // Tangential only. |
1472 | |
1473 | const dng_point_real64 tan = fParams->EvaluateTangential (plane, |
1474 | rr, |
1475 | diffNormScaled, |
1476 | diffNormSqr); |
1477 | |
1478 | dSrc.h = diff.h + (fNormRadius * tan.h); |
1479 | dSrc.v = diff.v + (fNormRadius * tan.v * fPixelScaleVInv); |
1480 | |
1481 | } |
1482 | |
1483 | else |
1484 | { |
1485 | |
1486 | // Radial and tangential. |
1487 | |
1488 | const real64 ratio = fParams->EvaluateRatio (plane, |
1489 | rr); |
1490 | |
1491 | const dng_point_real64 tan = fParams->EvaluateTangential (plane, |
1492 | rr, |
1493 | diffNormScaled, |
1494 | diffNormSqr); |
1495 | |
1496 | dSrc.h = fNormRadius * (diffNorm.h * ratio + tan.h); |
1497 | dSrc.v = fNormRadius * (diffNorm.v * ratio + tan.v * fPixelScaleVInv); |
1498 | |
1499 | } |
1500 | |
1501 | return fCenter + dSrc; |
1502 | |
1503 | } |
1504 | |
1505 | /*****************************************************************************/ |
1506 | |
1507 | dng_opcode_WarpRectilinear::dng_opcode_WarpRectilinear (const dng_warp_params_rectilinear ¶ms, |
1508 | uint32 flags) |
1509 | |
1510 | : dng_opcode (dngOpcode_WarpRectilinear, |
1511 | dngVersion_1_3_0_0, |
1512 | flags) |
1513 | |
1514 | , fWarpParams (params) |
1515 | |
1516 | { |
1517 | |
1518 | if (!params.IsValid ()) |
1519 | { |
1520 | ThrowBadFormat (); |
1521 | } |
1522 | |
1523 | } |
1524 | |
1525 | /*****************************************************************************/ |
1526 | |
1527 | dng_opcode_WarpRectilinear::dng_opcode_WarpRectilinear (dng_stream &stream) |
1528 | |
1529 | : dng_opcode (dngOpcode_WarpRectilinear, |
1530 | stream, |
1531 | "WarpRectilinear" ) |
1532 | |
1533 | , fWarpParams () |
1534 | |
1535 | { |
1536 | |
1537 | // Grab the size in bytes. |
1538 | |
1539 | const uint32 bytes = stream.Get_uint32 (); |
1540 | |
1541 | // Grab the number of planes to warp. |
1542 | |
1543 | fWarpParams.fPlanes = stream.Get_uint32 (); |
1544 | |
1545 | // Verify number of planes. |
1546 | |
1547 | if (fWarpParams.fPlanes == 0 || |
1548 | fWarpParams.fPlanes > kMaxColorPlanes) |
1549 | { |
1550 | ThrowBadFormat (); |
1551 | } |
1552 | |
1553 | // Verify the size. |
1554 | |
1555 | if (bytes != ParamBytes (fWarpParams.fPlanes)) |
1556 | { |
1557 | ThrowBadFormat (); |
1558 | } |
1559 | |
1560 | // Read warp parameters for each plane. |
1561 | |
1562 | for (uint32 plane = 0; plane < fWarpParams.fPlanes; plane++) |
1563 | { |
1564 | |
1565 | fWarpParams.fRadParams [plane][0] = stream.Get_real64 (); |
1566 | fWarpParams.fRadParams [plane][1] = stream.Get_real64 (); |
1567 | fWarpParams.fRadParams [plane][2] = stream.Get_real64 (); |
1568 | fWarpParams.fRadParams [plane][3] = stream.Get_real64 (); |
1569 | |
1570 | fWarpParams.fTanParams [plane][0] = stream.Get_real64 (); |
1571 | fWarpParams.fTanParams [plane][1] = stream.Get_real64 (); |
1572 | |
1573 | } |
1574 | |
1575 | // Read the image center. |
1576 | |
1577 | fWarpParams.fCenter.h = stream.Get_real64 (); |
1578 | fWarpParams.fCenter.v = stream.Get_real64 (); |
1579 | |
1580 | #if qDNGValidate |
1581 | |
1582 | if (gVerbose) |
1583 | { |
1584 | |
1585 | fWarpParams.Dump (); |
1586 | |
1587 | } |
1588 | |
1589 | #endif |
1590 | |
1591 | if (!fWarpParams.IsValid ()) |
1592 | { |
1593 | ThrowBadFormat (); |
1594 | } |
1595 | |
1596 | } |
1597 | |
1598 | /*****************************************************************************/ |
1599 | |
1600 | bool dng_opcode_WarpRectilinear::IsNOP () const |
1601 | { |
1602 | |
1603 | return fWarpParams.IsNOPAll (); |
1604 | |
1605 | } |
1606 | |
1607 | /*****************************************************************************/ |
1608 | |
1609 | bool dng_opcode_WarpRectilinear::IsValidForNegative (const dng_negative &negative) const |
1610 | { |
1611 | |
1612 | return fWarpParams.IsValidForNegative (negative); |
1613 | |
1614 | } |
1615 | |
1616 | /*****************************************************************************/ |
1617 | |
1618 | void dng_opcode_WarpRectilinear::PutData (dng_stream &stream) const |
1619 | { |
1620 | |
1621 | const uint32 bytes = ParamBytes (fWarpParams.fPlanes); |
1622 | |
1623 | stream.Put_uint32 (bytes); |
1624 | |
1625 | stream.Put_uint32 (fWarpParams.fPlanes); |
1626 | |
1627 | for (uint32 plane = 0; plane < fWarpParams.fPlanes; plane++) |
1628 | { |
1629 | |
1630 | stream.Put_real64 (fWarpParams.fRadParams [plane][0]); |
1631 | stream.Put_real64 (fWarpParams.fRadParams [plane][1]); |
1632 | stream.Put_real64 (fWarpParams.fRadParams [plane][2]); |
1633 | stream.Put_real64 (fWarpParams.fRadParams [plane][3]); |
1634 | |
1635 | stream.Put_real64 (fWarpParams.fTanParams [plane][0]); |
1636 | stream.Put_real64 (fWarpParams.fTanParams [plane][1]); |
1637 | |
1638 | } |
1639 | |
1640 | stream.Put_real64 (fWarpParams.fCenter.h); |
1641 | stream.Put_real64 (fWarpParams.fCenter.v); |
1642 | |
1643 | } |
1644 | |
1645 | /*****************************************************************************/ |
1646 | |
1647 | void dng_opcode_WarpRectilinear::Apply (dng_host &host, |
1648 | dng_negative &negative, |
1649 | AutoPtr<dng_image> &image) |
1650 | { |
1651 | |
1652 | #if qDNGValidate |
1653 | |
1654 | dng_timer timer ("WarpRectilinear time" ); |
1655 | |
1656 | #endif |
1657 | |
1658 | // Allocate destination image. |
1659 | |
1660 | AutoPtr<dng_image> dstImage (host.Make_dng_image (image->Bounds (), |
1661 | image->Planes (), |
1662 | image->PixelType ())); |
1663 | |
1664 | // Warp the image. |
1665 | |
1666 | AutoPtr<dng_warp_params> params (new dng_warp_params_rectilinear (fWarpParams)); |
1667 | |
1668 | dng_filter_warp filter (*image, |
1669 | *dstImage, |
1670 | negative, |
1671 | params); |
1672 | |
1673 | filter.Initialize (host); |
1674 | |
1675 | host.PerformAreaTask (filter, |
1676 | image->Bounds ()); |
1677 | |
1678 | // Return the new image. |
1679 | |
1680 | image.Reset (dstImage.Release ()); |
1681 | |
1682 | } |
1683 | |
1684 | /*****************************************************************************/ |
1685 | |
1686 | uint32 dng_opcode_WarpRectilinear::ParamBytes (uint32 planes) |
1687 | { |
1688 | |
1689 | return (1 * (uint32) sizeof (uint32) ) + // Number of planes. |
1690 | (6 * (uint32) sizeof (real64) * planes) + // Warp coefficients. |
1691 | (2 * (uint32) sizeof (real64) ); // Optical center. |
1692 | |
1693 | } |
1694 | |
1695 | /*****************************************************************************/ |
1696 | |
1697 | dng_opcode_WarpFisheye::dng_opcode_WarpFisheye (const dng_warp_params_fisheye ¶ms, |
1698 | uint32 flags) |
1699 | |
1700 | : dng_opcode (dngOpcode_WarpFisheye, |
1701 | dngVersion_1_3_0_0, |
1702 | flags) |
1703 | |
1704 | , fWarpParams (params) |
1705 | |
1706 | { |
1707 | |
1708 | if (!params.IsValid ()) |
1709 | { |
1710 | ThrowBadFormat (); |
1711 | } |
1712 | |
1713 | } |
1714 | |
1715 | /*****************************************************************************/ |
1716 | |
1717 | dng_opcode_WarpFisheye::dng_opcode_WarpFisheye (dng_stream &stream) |
1718 | |
1719 | : dng_opcode (dngOpcode_WarpFisheye, |
1720 | stream, |
1721 | "WarpFisheye" ) |
1722 | |
1723 | , fWarpParams () |
1724 | |
1725 | { |
1726 | |
1727 | // Grab the size in bytes. |
1728 | |
1729 | const uint32 bytes = stream.Get_uint32 (); |
1730 | |
1731 | // Grab the number of planes to warp. |
1732 | |
1733 | fWarpParams.fPlanes = stream.Get_uint32 (); |
1734 | |
1735 | // Verify number of planes. |
1736 | |
1737 | if (fWarpParams.fPlanes == 0 || |
1738 | fWarpParams.fPlanes > kMaxColorPlanes) |
1739 | { |
1740 | ThrowBadFormat (); |
1741 | } |
1742 | |
1743 | // Verify the size. |
1744 | |
1745 | if (bytes != ParamBytes (fWarpParams.fPlanes)) |
1746 | { |
1747 | ThrowBadFormat (); |
1748 | } |
1749 | |
1750 | // Read warp parameters for each plane. |
1751 | |
1752 | for (uint32 plane = 0; plane < fWarpParams.fPlanes; plane++) |
1753 | { |
1754 | |
1755 | fWarpParams.fRadParams [plane][0] = stream.Get_real64 (); |
1756 | fWarpParams.fRadParams [plane][1] = stream.Get_real64 (); |
1757 | fWarpParams.fRadParams [plane][2] = stream.Get_real64 (); |
1758 | fWarpParams.fRadParams [plane][3] = stream.Get_real64 (); |
1759 | |
1760 | } |
1761 | |
1762 | // Read the image center. |
1763 | |
1764 | fWarpParams.fCenter.h = stream.Get_real64 (); |
1765 | fWarpParams.fCenter.v = stream.Get_real64 (); |
1766 | |
1767 | #if qDNGValidate |
1768 | |
1769 | if (gVerbose) |
1770 | { |
1771 | |
1772 | fWarpParams.Dump (); |
1773 | |
1774 | } |
1775 | |
1776 | #endif |
1777 | |
1778 | if (!fWarpParams.IsValid ()) |
1779 | { |
1780 | ThrowBadFormat (); |
1781 | } |
1782 | |
1783 | } |
1784 | |
1785 | /*****************************************************************************/ |
1786 | |
1787 | bool dng_opcode_WarpFisheye::IsNOP () const |
1788 | { |
1789 | |
1790 | return fWarpParams.IsNOPAll (); |
1791 | |
1792 | } |
1793 | |
1794 | /*****************************************************************************/ |
1795 | |
1796 | bool dng_opcode_WarpFisheye::IsValidForNegative (const dng_negative &negative) const |
1797 | { |
1798 | |
1799 | return fWarpParams.IsValidForNegative (negative); |
1800 | |
1801 | } |
1802 | |
1803 | /*****************************************************************************/ |
1804 | |
1805 | void dng_opcode_WarpFisheye::PutData (dng_stream &stream) const |
1806 | { |
1807 | |
1808 | const uint32 bytes = ParamBytes (fWarpParams.fPlanes); |
1809 | |
1810 | stream.Put_uint32 (bytes); |
1811 | |
1812 | // Write the number of planes. |
1813 | |
1814 | stream.Put_uint32 (fWarpParams.fPlanes); |
1815 | |
1816 | // Write the warp coefficients. |
1817 | |
1818 | for (uint32 plane = 0; plane < fWarpParams.fPlanes; plane++) |
1819 | { |
1820 | |
1821 | stream.Put_real64 (fWarpParams.fRadParams [plane][0]); |
1822 | stream.Put_real64 (fWarpParams.fRadParams [plane][1]); |
1823 | stream.Put_real64 (fWarpParams.fRadParams [plane][2]); |
1824 | stream.Put_real64 (fWarpParams.fRadParams [plane][3]); |
1825 | |
1826 | } |
1827 | |
1828 | // Write the optical center. |
1829 | |
1830 | stream.Put_real64 (fWarpParams.fCenter.h); |
1831 | stream.Put_real64 (fWarpParams.fCenter.v); |
1832 | |
1833 | } |
1834 | |
1835 | /*****************************************************************************/ |
1836 | |
1837 | void dng_opcode_WarpFisheye::Apply (dng_host &host, |
1838 | dng_negative &negative, |
1839 | AutoPtr<dng_image> &image) |
1840 | { |
1841 | |
1842 | #if qDNGValidate |
1843 | |
1844 | dng_timer timer ("WarpFisheye time" ); |
1845 | |
1846 | #endif |
1847 | |
1848 | // Allocate destination image. |
1849 | |
1850 | AutoPtr<dng_image> dstImage (host.Make_dng_image (image->Bounds (), |
1851 | image->Planes (), |
1852 | image->PixelType ())); |
1853 | |
1854 | // Warp the image. |
1855 | |
1856 | AutoPtr<dng_warp_params> params (new dng_warp_params_fisheye (fWarpParams)); |
1857 | |
1858 | dng_filter_warp filter (*image, |
1859 | *dstImage, |
1860 | negative, |
1861 | params); |
1862 | |
1863 | filter.Initialize (host); |
1864 | |
1865 | host.PerformAreaTask (filter, |
1866 | image->Bounds ()); |
1867 | |
1868 | // Return the new image. |
1869 | |
1870 | image.Reset (dstImage.Release ()); |
1871 | |
1872 | } |
1873 | |
1874 | /*****************************************************************************/ |
1875 | |
1876 | uint32 dng_opcode_WarpFisheye::ParamBytes (uint32 planes) |
1877 | { |
1878 | |
1879 | return (1 * (uint32) sizeof (uint32) ) + // Number of planes. |
1880 | (4 * (uint32) sizeof (real64) * planes) + // Warp coefficients. |
1881 | (2 * (uint32) sizeof (real64) ); // Optical center. |
1882 | |
1883 | } |
1884 | |
1885 | /*****************************************************************************/ |
1886 | |
1887 | dng_vignette_radial_params::dng_vignette_radial_params () |
1888 | |
1889 | : fParams (kNumTerms) |
1890 | , fCenter (0.5, 0.5) |
1891 | |
1892 | { |
1893 | |
1894 | } |
1895 | |
1896 | /*****************************************************************************/ |
1897 | |
1898 | dng_vignette_radial_params::dng_vignette_radial_params (const dng_std_vector<real64> ¶ms, |
1899 | const dng_point_real64 ¢er) |
1900 | |
1901 | : fParams (params) |
1902 | , fCenter (center) |
1903 | |
1904 | { |
1905 | |
1906 | } |
1907 | |
1908 | /*****************************************************************************/ |
1909 | |
1910 | bool dng_vignette_radial_params::IsNOP () const |
1911 | { |
1912 | |
1913 | for (uint32 i = 0; i < fParams.size (); i++) |
1914 | { |
1915 | |
1916 | if (fParams [i] != 0.0) |
1917 | { |
1918 | return false; |
1919 | } |
1920 | |
1921 | } |
1922 | |
1923 | return true; |
1924 | |
1925 | } |
1926 | |
1927 | /*****************************************************************************/ |
1928 | |
1929 | bool dng_vignette_radial_params::IsValid () const |
1930 | { |
1931 | |
1932 | if (fParams.size () != kNumTerms) |
1933 | { |
1934 | return false; |
1935 | } |
1936 | |
1937 | if (fCenter.h < 0.0 || |
1938 | fCenter.h > 1.0 || |
1939 | fCenter.v < 0.0 || |
1940 | fCenter.v > 1.0) |
1941 | { |
1942 | return false; |
1943 | } |
1944 | |
1945 | return true; |
1946 | |
1947 | } |
1948 | |
1949 | /*****************************************************************************/ |
1950 | |
1951 | void dng_vignette_radial_params::Dump () const |
1952 | { |
1953 | |
1954 | #if qDNGValidate |
1955 | |
1956 | printf (" Radial vignette params: " ); |
1957 | |
1958 | for (uint32 i = 0; i < fParams.size (); i++) |
1959 | { |
1960 | |
1961 | printf ("%s%.6lf" , |
1962 | (i == 0) ? "" : ", " , |
1963 | fParams [i]); |
1964 | |
1965 | } |
1966 | |
1967 | printf ("\n" ); |
1968 | |
1969 | printf (" Optical center:\n" |
1970 | " h = %.6lf\n" |
1971 | " v = %.6lf\n" , |
1972 | fCenter.h, |
1973 | fCenter.v); |
1974 | |
1975 | #endif |
1976 | |
1977 | } |
1978 | |
1979 | /*****************************************************************************/ |
1980 | |
1981 | class dng_vignette_radial_function: public dng_1d_function |
1982 | { |
1983 | |
1984 | protected: |
1985 | |
1986 | const dng_vignette_radial_params fParams; |
1987 | |
1988 | public: |
1989 | |
1990 | explicit dng_vignette_radial_function (const dng_vignette_radial_params ¶ms) |
1991 | |
1992 | : fParams (params) |
1993 | |
1994 | { |
1995 | |
1996 | } |
1997 | |
1998 | // x represents r^2, where r is the normalized Euclidean distance from the |
1999 | // optical center to a pixel. r lies in [0,1], so r^2 (and hence x) also lies |
2000 | // in [0,1]. |
2001 | |
2002 | virtual real64 Evaluate (real64 x) const |
2003 | { |
2004 | |
2005 | DNG_REQUIRE (fParams.fParams.size () == dng_vignette_radial_params::kNumTerms, |
2006 | "Bad number of vignette opcode coefficients." ); |
2007 | |
2008 | real64 sum = 0.0; |
2009 | |
2010 | const dng_std_vector<real64> &v = fParams.fParams; |
2011 | |
2012 | for (dng_std_vector<real64>::const_reverse_iterator i = v.rbegin (); i != v.rend (); i++) |
2013 | { |
2014 | sum = x * ((*i) + sum); |
2015 | } |
2016 | |
2017 | sum += 1.0; |
2018 | |
2019 | return sum; |
2020 | |
2021 | } |
2022 | |
2023 | }; |
2024 | |
2025 | /*****************************************************************************/ |
2026 | |
2027 | dng_opcode_FixVignetteRadial::dng_opcode_FixVignetteRadial (const dng_vignette_radial_params ¶ms, |
2028 | uint32 flags) |
2029 | |
2030 | : dng_inplace_opcode (dngOpcode_FixVignetteRadial, |
2031 | dngVersion_1_3_0_0, |
2032 | flags) |
2033 | |
2034 | , fParams (params) |
2035 | |
2036 | , fImagePlanes (1) |
2037 | |
2038 | , fSrcOriginH (0) |
2039 | , fSrcOriginV (0) |
2040 | |
2041 | , fSrcStepH (0) |
2042 | , fSrcStepV (0) |
2043 | |
2044 | , fTableInputBits (0) |
2045 | , fTableOutputBits (0) |
2046 | |
2047 | , fGainTable () |
2048 | |
2049 | { |
2050 | |
2051 | if (!params.IsValid ()) |
2052 | { |
2053 | ThrowBadFormat (); |
2054 | } |
2055 | |
2056 | } |
2057 | |
2058 | /*****************************************************************************/ |
2059 | |
2060 | dng_opcode_FixVignetteRadial::dng_opcode_FixVignetteRadial (dng_stream &stream) |
2061 | |
2062 | : dng_inplace_opcode (dngOpcode_FixVignetteRadial, |
2063 | stream, |
2064 | "FixVignetteRadial" ) |
2065 | |
2066 | , fParams () |
2067 | |
2068 | , fImagePlanes (1) |
2069 | |
2070 | , fSrcOriginH (0) |
2071 | , fSrcOriginV (0) |
2072 | |
2073 | , fSrcStepH (0) |
2074 | , fSrcStepV (0) |
2075 | |
2076 | , fTableInputBits (0) |
2077 | , fTableOutputBits (0) |
2078 | |
2079 | , fGainTable () |
2080 | |
2081 | { |
2082 | |
2083 | // Grab the size in bytes. |
2084 | |
2085 | const uint32 bytes = stream.Get_uint32 (); |
2086 | |
2087 | // Verify the size. |
2088 | |
2089 | if (bytes != ParamBytes ()) |
2090 | { |
2091 | ThrowBadFormat (); |
2092 | } |
2093 | |
2094 | // Read vignette coefficients. |
2095 | |
2096 | fParams.fParams = dng_std_vector<real64> (dng_vignette_radial_params::kNumTerms); |
2097 | |
2098 | for (uint32 i = 0; i < dng_vignette_radial_params::kNumTerms; i++) |
2099 | { |
2100 | fParams.fParams [i] = stream.Get_real64 (); |
2101 | } |
2102 | |
2103 | // Read the image center. |
2104 | |
2105 | fParams.fCenter.h = stream.Get_real64 (); |
2106 | fParams.fCenter.v = stream.Get_real64 (); |
2107 | |
2108 | // Debug. |
2109 | |
2110 | #if qDNGValidate |
2111 | |
2112 | if (gVerbose) |
2113 | { |
2114 | |
2115 | fParams.Dump (); |
2116 | |
2117 | } |
2118 | |
2119 | #endif |
2120 | |
2121 | if (!fParams.IsValid ()) |
2122 | { |
2123 | ThrowBadFormat (); |
2124 | } |
2125 | |
2126 | } |
2127 | |
2128 | /*****************************************************************************/ |
2129 | |
2130 | bool dng_opcode_FixVignetteRadial::IsNOP () const |
2131 | { |
2132 | |
2133 | return fParams.IsNOP (); |
2134 | |
2135 | } |
2136 | |
2137 | /*****************************************************************************/ |
2138 | |
2139 | bool dng_opcode_FixVignetteRadial::IsValidForNegative (const dng_negative & /* negative */) const |
2140 | { |
2141 | |
2142 | return fParams.IsValid (); |
2143 | |
2144 | } |
2145 | |
2146 | /*****************************************************************************/ |
2147 | |
2148 | void dng_opcode_FixVignetteRadial::PutData (dng_stream &stream) const |
2149 | { |
2150 | |
2151 | const uint32 bytes = ParamBytes (); |
2152 | |
2153 | stream.Put_uint32 (bytes); |
2154 | |
2155 | DNG_REQUIRE (fParams.fParams.size () == dng_vignette_radial_params::kNumTerms, |
2156 | "Bad number of vignette opcode coefficients." ); |
2157 | |
2158 | for (uint32 i = 0; i < dng_vignette_radial_params::kNumTerms; i++) |
2159 | { |
2160 | stream.Put_real64 (fParams.fParams [i]); |
2161 | } |
2162 | |
2163 | stream.Put_real64 (fParams.fCenter.h); |
2164 | stream.Put_real64 (fParams.fCenter.v); |
2165 | |
2166 | } |
2167 | |
2168 | /*****************************************************************************/ |
2169 | |
2170 | void dng_opcode_FixVignetteRadial::Prepare (dng_negative &negative, |
2171 | uint32 threadCount, |
2172 | const dng_point &tileSize, |
2173 | const dng_rect &imageBounds, |
2174 | uint32 imagePlanes, |
2175 | uint32 bufferPixelType, |
2176 | dng_memory_allocator &allocator) |
2177 | { |
2178 | |
2179 | // This opcode is restricted to 32-bit images. |
2180 | |
2181 | if (bufferPixelType != ttFloat) |
2182 | { |
2183 | ThrowBadFormat (); |
2184 | } |
2185 | |
2186 | // Sanity check number of planes. |
2187 | |
2188 | DNG_ASSERT (imagePlanes >= 1 && imagePlanes <= kMaxColorPlanes, |
2189 | "Bad number of planes." ); |
2190 | |
2191 | if (imagePlanes < 1 || imagePlanes > kMaxColorPlanes) |
2192 | { |
2193 | ThrowProgramError (); |
2194 | } |
2195 | |
2196 | fImagePlanes = imagePlanes; |
2197 | |
2198 | // Set the vignette correction curve. |
2199 | |
2200 | const dng_vignette_radial_function curve (fParams); |
2201 | |
2202 | // Grab the destination image area. |
2203 | |
2204 | const dng_rect_real64 bounds (imageBounds); |
2205 | |
2206 | // Determine the optical center and maximum radius in pixel coordinates. |
2207 | |
2208 | const dng_point_real64 centerPixel (Lerp_real64 (bounds.t, |
2209 | bounds.b, |
2210 | fParams.fCenter.v), |
2211 | |
2212 | Lerp_real64 (bounds.l, |
2213 | bounds.r, |
2214 | fParams.fCenter.h)); |
2215 | |
2216 | const real64 pixelScaleV = 1.0 / negative.PixelAspectRatio (); |
2217 | |
2218 | const real64 maxRadius = hypot (Max_real64 (Abs_real64 (centerPixel.v - bounds.t), |
2219 | Abs_real64 (centerPixel.v - bounds.b)) * pixelScaleV, |
2220 | |
2221 | Max_real64 (Abs_real64 (centerPixel.h - bounds.l), |
2222 | Abs_real64 (centerPixel.h - bounds.r))); |
2223 | |
2224 | const dng_point_real64 radius (maxRadius, |
2225 | maxRadius); |
2226 | |
2227 | // Find origin and scale. |
2228 | |
2229 | const real64 pixelScaleH = 1.0; |
2230 | |
2231 | fSrcOriginH = Real64ToFixed64 (-centerPixel.h * pixelScaleH / radius.h); |
2232 | fSrcOriginV = Real64ToFixed64 (-centerPixel.v * pixelScaleV / radius.v); |
2233 | |
2234 | fSrcStepH = Real64ToFixed64 (pixelScaleH / radius.h); |
2235 | fSrcStepV = Real64ToFixed64 (pixelScaleV / radius.v); |
2236 | |
2237 | // Adjust for pixel centers. |
2238 | |
2239 | fSrcOriginH += fSrcStepH >> 1; |
2240 | fSrcOriginV += fSrcStepV >> 1; |
2241 | |
2242 | // Evaluate 32-bit vignette correction table. |
2243 | |
2244 | dng_1d_table table32; |
2245 | |
2246 | table32.Initialize (allocator, |
2247 | curve, |
2248 | false); |
2249 | |
2250 | // Find maximum scale factor. |
2251 | |
2252 | const real64 maxScale = Max_real32 (table32.Interpolate (0.0f), |
2253 | table32.Interpolate (1.0f)); |
2254 | |
2255 | // Find table input bits. |
2256 | |
2257 | fTableInputBits = 16; |
2258 | |
2259 | // Find table output bits. |
2260 | |
2261 | fTableOutputBits = 15; |
2262 | |
2263 | while ((1 << fTableOutputBits) * maxScale > 65535.0) |
2264 | { |
2265 | fTableOutputBits--; |
2266 | } |
2267 | |
2268 | // Allocate 16-bit scale table. |
2269 | |
2270 | const uint32 tableEntries = (1 << fTableInputBits) + 1; |
2271 | |
2272 | fGainTable.Reset (allocator.Allocate (tableEntries * (uint32) sizeof (uint16))); |
2273 | |
2274 | uint16 *table16 = fGainTable->Buffer_uint16 (); |
2275 | |
2276 | // Interpolate 32-bit table into 16-bit table. |
2277 | |
2278 | const real32 scale0 = 1.0f / (1 << fTableInputBits ); |
2279 | const real32 scale1 = 1.0f * (1 << fTableOutputBits); |
2280 | |
2281 | for (uint32 index = 0; index < tableEntries; index++) |
2282 | { |
2283 | |
2284 | real32 x = index * scale0; |
2285 | |
2286 | real32 y = table32.Interpolate (x) * scale1; |
2287 | |
2288 | table16 [index] = (uint16) Round_uint32 (y); |
2289 | |
2290 | } |
2291 | |
2292 | // Prepare vignette mask buffers. |
2293 | |
2294 | { |
2295 | |
2296 | const uint32 pixelType = ttShort; |
2297 | const uint32 bufferSize = ComputeBufferSize(pixelType, tileSize, |
2298 | imagePlanes, pad16Bytes); |
2299 | |
2300 | for (uint32 threadIndex = 0; threadIndex < threadCount; threadIndex++) |
2301 | { |
2302 | |
2303 | fMaskBuffers [threadIndex] . Reset (allocator.Allocate (bufferSize)); |
2304 | |
2305 | } |
2306 | |
2307 | } |
2308 | |
2309 | } |
2310 | |
2311 | /*****************************************************************************/ |
2312 | |
2313 | void dng_opcode_FixVignetteRadial::ProcessArea (dng_negative & /* negative */, |
2314 | uint32 threadIndex, |
2315 | dng_pixel_buffer &buffer, |
2316 | const dng_rect &dstArea, |
2317 | const dng_rect & /* imageBounds */) |
2318 | { |
2319 | |
2320 | // Setup mask pixel buffer. |
2321 | |
2322 | dng_pixel_buffer maskPixelBuffer(dstArea, 0, fImagePlanes, ttShort, |
2323 | pcRowInterleavedAlign16, |
2324 | fMaskBuffers [threadIndex]->Buffer ()); |
2325 | |
2326 | // Compute mask. |
2327 | |
2328 | DoVignetteMask16 (maskPixelBuffer.DirtyPixel_uint16 (dstArea.t, dstArea.l), |
2329 | dstArea.H (), |
2330 | dstArea.W (), |
2331 | maskPixelBuffer.RowStep (), |
2332 | fSrcOriginH + fSrcStepH * dstArea.l, |
2333 | fSrcOriginV + fSrcStepV * dstArea.t, |
2334 | fSrcStepH, |
2335 | fSrcStepV, |
2336 | fTableInputBits, |
2337 | fGainTable->Buffer_uint16 ()); |
2338 | |
2339 | // Apply mask. |
2340 | |
2341 | DoVignette32 (buffer.DirtyPixel_real32 (dstArea.t, dstArea.l), |
2342 | maskPixelBuffer.ConstPixel_uint16 (dstArea.t, dstArea.l), |
2343 | dstArea.H (), |
2344 | dstArea.W (), |
2345 | fImagePlanes, |
2346 | buffer.RowStep (), |
2347 | buffer.PlaneStep (), |
2348 | maskPixelBuffer.RowStep (), |
2349 | fTableOutputBits); |
2350 | |
2351 | } |
2352 | |
2353 | /*****************************************************************************/ |
2354 | |
2355 | uint32 dng_opcode_FixVignetteRadial::ParamBytes () |
2356 | { |
2357 | |
2358 | const uint32 N = dng_vignette_radial_params::kNumTerms; |
2359 | |
2360 | return ((N * sizeof (real64)) + // Vignette coefficients. |
2361 | (2 * sizeof (real64))); // Optical center. |
2362 | |
2363 | } |
2364 | |
2365 | /*****************************************************************************/ |
2366 | |