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
36dng_warp_params::dng_warp_params ()
37
38 : fPlanes (1)
39 , fCenter (0.5, 0.5)
40
41 {
42
43 }
44
45/*****************************************************************************/
46
47dng_warp_params::dng_warp_params (uint32 planes,
48 const dng_point_real64 &center)
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
68dng_warp_params::~dng_warp_params ()
69 {
70
71 }
72
73/*****************************************************************************/
74
75bool dng_warp_params::IsNOPAll () const
76 {
77
78 return IsRadNOPAll () &&
79 IsTanNOPAll ();
80
81 }
82
83/*****************************************************************************/
84
85bool dng_warp_params::IsNOP (uint32 plane) const
86 {
87
88 return IsRadNOP (plane) &&
89 IsTanNOP (plane);
90
91 }
92
93/*****************************************************************************/
94
95bool 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
114bool dng_warp_params::IsRadNOP (uint32 /* plane */) const
115 {
116
117 return false;
118
119 }
120
121/*****************************************************************************/
122
123bool 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
142bool dng_warp_params::IsTanNOP (uint32 /* plane */) const
143 {
144
145 return false;
146
147 }
148
149/*****************************************************************************/
150
151bool 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
177bool 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
197real64 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
241dng_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
262dng_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
279void 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
298dng_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
318dng_warp_params_rectilinear::dng_warp_params_rectilinear (uint32 planes,
319 const dng_vector radParams [],
320 const dng_vector tanParams [],
321 const dng_point_real64 &center)
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
338dng_warp_params_rectilinear::~dng_warp_params_rectilinear ()
339 {
340
341 }
342
343/*****************************************************************************/
344
345bool 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
361bool 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
375bool 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
399void 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
414real64 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
428real64 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
440dng_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
462real64 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
634dng_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
688void 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
718dng_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
735dng_warp_params_fisheye::dng_warp_params_fisheye (uint32 planes,
736 const dng_vector radParams [],
737 const dng_point_real64 &center)
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
754dng_warp_params_fisheye::~dng_warp_params_fisheye ()
755 {
756
757 }
758
759/*****************************************************************************/
760
761bool dng_warp_params_fisheye::IsRadNOP (uint32 /* plane */) const
762 {
763
764 return false;
765
766 }
767
768/*****************************************************************************/
769
770bool dng_warp_params_fisheye::IsTanNOP (uint32 /* plane */) const
771 {
772
773 return true;
774
775 }
776
777/*****************************************************************************/
778
779bool 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
798void 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
812real64 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
828real64 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
851dng_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
867real64 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
919dng_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
931void 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
957class 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> &params);
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
1001dng_filter_warp::dng_filter_warp (const dng_image &srcImage,
1002 dng_image &dstImage,
1003 const dng_negative &negative,
1004 AutoPtr<dng_warp_params> &params)
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
1101void 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
1115dng_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
1224dng_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
1293void 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
1435dng_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
1507dng_opcode_WarpRectilinear::dng_opcode_WarpRectilinear (const dng_warp_params_rectilinear &params,
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
1527dng_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
1600bool dng_opcode_WarpRectilinear::IsNOP () const
1601 {
1602
1603 return fWarpParams.IsNOPAll ();
1604
1605 }
1606
1607/*****************************************************************************/
1608
1609bool dng_opcode_WarpRectilinear::IsValidForNegative (const dng_negative &negative) const
1610 {
1611
1612 return fWarpParams.IsValidForNegative (negative);
1613
1614 }
1615
1616/*****************************************************************************/
1617
1618void 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
1647void 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
1686uint32 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
1697dng_opcode_WarpFisheye::dng_opcode_WarpFisheye (const dng_warp_params_fisheye &params,
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
1717dng_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
1787bool dng_opcode_WarpFisheye::IsNOP () const
1788 {
1789
1790 return fWarpParams.IsNOPAll ();
1791
1792 }
1793
1794/*****************************************************************************/
1795
1796bool dng_opcode_WarpFisheye::IsValidForNegative (const dng_negative &negative) const
1797 {
1798
1799 return fWarpParams.IsValidForNegative (negative);
1800
1801 }
1802
1803/*****************************************************************************/
1804
1805void 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
1837void 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
1876uint32 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
1887dng_vignette_radial_params::dng_vignette_radial_params ()
1888
1889 : fParams (kNumTerms)
1890 , fCenter (0.5, 0.5)
1891
1892 {
1893
1894 }
1895
1896/*****************************************************************************/
1897
1898dng_vignette_radial_params::dng_vignette_radial_params (const dng_std_vector<real64> &params,
1899 const dng_point_real64 &center)
1900
1901 : fParams (params)
1902 , fCenter (center)
1903
1904 {
1905
1906 }
1907
1908/*****************************************************************************/
1909
1910bool 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
1929bool 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
1951void 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
1981class 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 &params)
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
2027dng_opcode_FixVignetteRadial::dng_opcode_FixVignetteRadial (const dng_vignette_radial_params &params,
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
2060dng_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
2130bool dng_opcode_FixVignetteRadial::IsNOP () const
2131 {
2132
2133 return fParams.IsNOP ();
2134
2135 }
2136
2137/*****************************************************************************/
2138
2139bool dng_opcode_FixVignetteRadial::IsValidForNegative (const dng_negative & /* negative */) const
2140 {
2141
2142 return fParams.IsValid ();
2143
2144 }
2145
2146/*****************************************************************************/
2147
2148void 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
2170void 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
2313void 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
2355uint32 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