1/*****************************************************************************/
2// Copyright 2006 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_resample.cpp#1 $ */
10/* $DateTime: 2012/05/30 13:28:51 $ */
11/* $Change: 832332 $ */
12/* $Author: tknoll $ */
13
14/*****************************************************************************/
15
16#include "dng_resample.h"
17
18#include "dng_assertions.h"
19#include "dng_bottlenecks.h"
20#include "dng_filter_task.h"
21#include "dng_host.h"
22#include "dng_image.h"
23#include "dng_memory.h"
24#include "dng_pixel_buffer.h"
25#include "dng_safe_arithmetic.h"
26#include "dng_tag_types.h"
27#include "dng_utils.h"
28
29/******************************************************************************/
30
31real64 dng_resample_bicubic::Extent () const
32	{
33
34	return 2.0;
35
36	}
37
38/******************************************************************************/
39
40real64 dng_resample_bicubic::Evaluate (real64 x) const
41	{
42
43	const real64 A = -0.75;
44
45	x = Abs_real64 (x);
46
47    if (x >= 2.0)
48        return 0.0;
49
50    else if (x >= 1.0)
51        return (((A * x - 5.0 * A) * x + 8.0 * A) * x - 4.0 * A);
52
53    else
54        return (((A + 2.0) * x - (A + 3.0)) * x * x + 1.0);
55
56	}
57
58/******************************************************************************/
59
60const dng_resample_function & dng_resample_bicubic::Get ()
61	{
62
63	static dng_resample_bicubic static_dng_resample_bicubic;
64
65	return static_dng_resample_bicubic;
66
67	}
68
69/*****************************************************************************/
70
71dng_resample_coords::dng_resample_coords ()
72
73	:	fOrigin (0)
74	,	fCoords ()
75
76	{
77
78	}
79
80/*****************************************************************************/
81
82dng_resample_coords::~dng_resample_coords ()
83	{
84
85	}
86
87/*****************************************************************************/
88
89void dng_resample_coords::Initialize (int32 srcOrigin,
90									  int32 dstOrigin,
91									  uint32 srcCount,
92									  uint32 dstCount,
93									  dng_memory_allocator &allocator)
94	{
95
96	fOrigin = dstOrigin;
97
98	uint32 dstEntries = 0;
99	uint32 bufferSize = 0;
100	if (!RoundUpUint32ToMultiple(dstCount, 8, &dstEntries) ||
101	    !SafeUint32Mult(dstEntries, sizeof(int32), &bufferSize)) {
102		ThrowMemoryFull("Arithmetic overflow computing size for coordinate "
103						"buffer");
104	}
105	fCoords.Reset (allocator.Allocate (bufferSize));
106
107	int32 *coords = fCoords->Buffer_int32 ();
108
109	real64 invScale = (real64) srcCount /
110					  (real64) dstCount;
111
112	for (uint32 j = 0; j < dstCount; j++)
113		{
114
115		real64 x = (real64) j + 0.5;
116
117		real64 y = x * invScale - 0.5 + (real64) srcOrigin;
118
119		coords [j] = Round_int32 (y * (real64) kResampleSubsampleCount);
120
121		}
122
123	// Pad out table by replicating last entry.
124
125	for (uint32 k = dstCount; k < dstEntries; k++)
126		{
127
128		coords [k] = coords [dstCount - 1];
129
130		}
131
132	}
133
134/*****************************************************************************/
135
136dng_resample_weights::dng_resample_weights ()
137
138	:	fRadius (0)
139
140	,	fWeightStep (0)
141
142	,	fWeights32 ()
143	,	fWeights16 ()
144
145	{
146
147	}
148
149/*****************************************************************************/
150
151dng_resample_weights::~dng_resample_weights ()
152	{
153
154	}
155
156/*****************************************************************************/
157
158void dng_resample_weights::Initialize (real64 scale,
159									   const dng_resample_function &kernel,
160									   dng_memory_allocator &allocator)
161	{
162
163	uint32 j;
164
165	// We only adjust the kernel size for scale factors less than 1.0.
166
167	scale = Min_real64 (scale, 1.0);
168
169	// Find radius of this kernel.
170
171	fRadius = (uint32) (kernel.Extent () / scale + 0.9999);
172
173	// Width is twice the radius.
174
175	uint32 width = fRadius * 2;
176
177	// Round to each set to weights to a multiple of 8 entries.
178
179	if (!RoundUpUint32ToMultiple (width, 8, &fWeightStep))
180		{
181
182		ThrowMemoryFull ("Arithmetic overflow computing fWeightStep");
183
184		}
185
186	// Allocate and zero weight tables.
187
188	uint32 bufferSize = 0;
189
190	if (!SafeUint32Mult (fWeightStep, kResampleSubsampleCount, &bufferSize) ||
191		 !SafeUint32Mult (bufferSize, (uint32) sizeof (real32), &bufferSize))
192		{
193
194		ThrowMemoryFull("Arithmetic overflow computing buffer size.");
195
196		}
197
198	fWeights32.Reset (allocator.Allocate (bufferSize));
199
200	DoZeroBytes (fWeights32->Buffer		 (),
201				 fWeights32->LogicalSize ());
202
203	if (!SafeUint32Mult (fWeightStep, kResampleSubsampleCount, &bufferSize) ||
204		 !SafeUint32Mult (bufferSize, (uint32) sizeof (int16), &bufferSize))
205		{
206
207		ThrowMemoryFull("Arithmetic overflow computing buffer size.");
208
209		}
210
211	fWeights16.Reset (allocator.Allocate (bufferSize));
212
213	DoZeroBytes (fWeights16->Buffer		 (),
214				 fWeights16->LogicalSize ());
215
216	// Compute kernel for each subsample values.
217
218	for (uint32 sample = 0; sample < kResampleSubsampleCount; sample++)
219		{
220
221		real64 fract = sample * (1.0 / (real64) kResampleSubsampleCount);
222
223		real32 *w32 = fWeights32->Buffer_real32 () + fWeightStep * sample;
224
225		// Evaluate kernel function for 32 bit weights.
226
227			{
228
229			real64 t32 = 0.0;
230
231			for (j = 0; j < width; j++)
232				{
233
234				int32 k = (int32) j - (int32) fRadius + 1;
235
236				real64 x = (k - fract) * scale;
237
238				w32 [j] = (real32) kernel.Evaluate (x);
239
240				t32 += w32 [j];
241
242				}
243
244			// Scale 32 bit weights so total of weights is 1.0.
245
246			real32 s32 = (real32) (1.0 / t32);
247
248			for (j = 0; j < width; j++)
249				{
250
251				w32 [j] *= s32;
252
253				}
254
255			}
256
257		// Round off 32 bit weights to 16 bit weights.
258
259			{
260
261			int16 *w16 = fWeights16->Buffer_int16 () + fWeightStep * sample;
262
263			int32 t16 = 0;
264
265			for (j = 0; j < width; j++)
266				{
267
268				w16 [j] = (int16) Round_int32 (w32 [j] * 16384.0);
269
270				t16 += w16 [j];
271
272				}
273
274			// Adjust center entry for any round off error so total is
275			// exactly 16384.
276
277			w16 [fRadius - (fract >= 0.5 ? 0 : 1)] += (int16) (16384 - t16);
278
279			}
280
281		}
282
283	}
284
285/*****************************************************************************/
286
287dng_resample_weights_2d::dng_resample_weights_2d ()
288
289	:	fRadius (0)
290
291	,	fRowStep (0)
292	,	fColStep (0)
293
294	,	fWeights32 ()
295	,	fWeights16 ()
296
297	{
298
299	}
300
301/*****************************************************************************/
302
303dng_resample_weights_2d::~dng_resample_weights_2d ()
304	{
305
306	}
307
308/*****************************************************************************/
309
310void dng_resample_weights_2d::Initialize (const dng_resample_function &kernel,
311										  dng_memory_allocator &allocator)
312	{
313
314	// Find radius of this kernel. Unlike with 1d resample weights (see
315	// dng_resample_weights), we never scale up the kernel size.
316
317	fRadius = (uint32) (kernel.Extent () + 0.9999);
318
319	// Width is twice the radius.
320
321	uint32 width    = 0;
322	uint32 widthSqr = 0;
323	uint32 step = 0;
324
325	if (!SafeUint32Mult (fRadius, 2, &width) ||
326		 !SafeUint32Mult (width, width, &widthSqr) ||
327		 !RoundUpUint32ToMultiple (widthSqr, 8, &step) ||
328		 !SafeUint32Mult (step, kResampleSubsampleCount2D, &fRowStep))
329			{
330
331			ThrowMemoryFull ("Arithmetic overflow computing row step.");
332
333			}
334
335	fColStep = step;
336
337	// Allocate and zero weight tables.
338
339	uint32 bufferSize = 0;
340
341	if (!SafeUint32Mult (step, kResampleSubsampleCount2D, &bufferSize) ||
342		 !SafeUint32Mult (bufferSize, kResampleSubsampleCount2D, &bufferSize) ||
343		 !SafeUint32Mult (bufferSize, (uint32) sizeof (real32), &bufferSize))
344		{
345
346		ThrowMemoryFull ("Arithmetic overflow computing buffer size.");
347
348		}
349
350	fWeights32.Reset (allocator.Allocate (bufferSize));
351
352	DoZeroBytes (fWeights32->Buffer		 (),
353				 fWeights32->LogicalSize ());
354
355
356	if (!SafeUint32Mult (step, kResampleSubsampleCount2D, &bufferSize) ||
357		 !SafeUint32Mult (bufferSize, kResampleSubsampleCount2D, &bufferSize) ||
358		 !SafeUint32Mult (bufferSize, (uint32) sizeof (int16), &bufferSize))
359		{
360
361		ThrowMemoryFull ("Arithmetic overflow computing buffer size.");
362
363		}
364
365	fWeights16.Reset (allocator.Allocate (bufferSize));
366
367	DoZeroBytes (fWeights16->Buffer		 (),
368				 fWeights16->LogicalSize ());
369
370	// Compute kernel for each subsample values.
371
372	for (uint32 y = 0; y < kResampleSubsampleCount2D; y++)
373		{
374
375		real64 yFract = y * (1.0 / (real64) kResampleSubsampleCount2D);
376
377		for (uint32 x = 0; x < kResampleSubsampleCount2D; x++)
378			{
379
380			real64 xFract = x * (1.0 / (real64) kResampleSubsampleCount2D);
381
382			real32 *w32 = (real32 *) Weights32 (dng_point ((int32) y,
383														   (int32) x));
384
385			// Evaluate kernel function for 32 bit weights.
386
387				{
388
389				real64 t32 = 0.0;
390
391				uint32 index = 0;
392
393				for (uint32 i = 0; i < width; i++)
394					{
395
396					int32 yInt = ((int32) i) - (int32) fRadius + 1;
397					real64 yPos = yInt - yFract;
398
399					for (uint32 j = 0; j < width; j++)
400						{
401
402						int32 xInt = ((int32) j) - (int32) fRadius + 1;
403						real64 xPos = xInt - xFract;
404
405						#if 0
406
407						// Radial.
408
409						real64 dy2 = yPos * yPos;
410						real64 dx2 = xPos * xPos;
411
412						real64 r = sqrt (dx2 + dy2);
413
414						w32 [index] = (real32) kernel.Evaluate (r);
415
416						#else
417
418						// Separable.
419
420						w32 [index] = (real32) kernel.Evaluate (xPos) *
421							          (real32) kernel.Evaluate (yPos);
422
423						#endif
424
425						t32 += w32 [index];
426
427						index++;
428
429						}
430
431					}
432
433				// Scale 32 bit weights so total of weights is 1.0.
434
435				const real32 s32 = (real32) (1.0 / t32);
436
437				for (uint32 i = 0; i < widthSqr; i++)
438					{
439
440					w32 [i] *= s32;
441
442					}
443
444				}
445
446			// Round off 32 bit weights to 16 bit weights.
447
448				{
449
450				int16 *w16 = (int16 *) Weights16 (dng_point ((int32) y,
451															 (int32) x));
452
453				int32 t16 = 0;
454
455				for (uint32 j = 0; j < widthSqr; j++)
456					{
457
458					w16 [j] = (int16) Round_int32 (w32 [j] * 16384.0);
459
460					t16 += w16 [j];
461
462					}
463
464				// Adjust one of the center entries for any round off error so total
465				// is exactly 16384.
466
467				const uint32 xOffset      = fRadius - ((xFract >= 0.5) ? 0 : 1);
468				const uint32 yOffset      = fRadius - ((yFract >= 0.5) ? 0 : 1);
469				const uint32 centerOffset = width * yOffset + xOffset;
470
471				w16 [centerOffset] += (int16) (16384 - t16);
472
473				}
474
475			}
476
477		}
478
479	}
480
481/*****************************************************************************/
482
483class dng_resample_task: public dng_filter_task
484	{
485
486	protected:
487
488		dng_rect fSrcBounds;
489		dng_rect fDstBounds;
490
491		const dng_resample_function &fKernel;
492
493		real64 fRowScale;
494		real64 fColScale;
495
496		dng_resample_coords fRowCoords;
497		dng_resample_coords fColCoords;
498
499		dng_resample_weights fWeightsV;
500		dng_resample_weights fWeightsH;
501
502		dng_point fSrcTileSize;
503
504		AutoPtr<dng_memory_block> fTempBuffer [kMaxMPThreads];
505
506	public:
507
508		dng_resample_task (const dng_image &srcImage,
509						   dng_image &dstImage,
510						   const dng_rect &srcBounds,
511						   const dng_rect &dstBounds,
512						   const dng_resample_function &kernel);
513
514		virtual dng_rect SrcArea (const dng_rect &dstArea);
515
516		virtual dng_point SrcTileSize (const dng_point &dstTileSize);
517
518		virtual void Start (uint32 threadCount,
519							const dng_point &tileSize,
520							dng_memory_allocator *allocator,
521							dng_abort_sniffer *sniffer);
522
523		virtual void ProcessArea (uint32 threadIndex,
524								  dng_pixel_buffer &srcBuffer,
525								  dng_pixel_buffer &dstBuffer);
526
527	};
528
529/*****************************************************************************/
530
531dng_resample_task::dng_resample_task (const dng_image &srcImage,
532						   			  dng_image &dstImage,
533						   			  const dng_rect &srcBounds,
534						   			  const dng_rect &dstBounds,
535									  const dng_resample_function &kernel)
536
537	:	dng_filter_task (srcImage,
538						 dstImage)
539
540	,	fSrcBounds (srcBounds)
541	,	fDstBounds (dstBounds)
542
543	,	fKernel (kernel)
544
545	,	fRowScale ((srcBounds.H () != 0) ? dstBounds.H () / (real64) srcBounds.H () : 0)
546	,	fColScale ((srcBounds.W () != 0) ? dstBounds.W () / (real64) srcBounds.W () : 0)
547
548	,	fRowCoords ()
549	,	fColCoords ()
550
551	,	fWeightsV ()
552	,	fWeightsH ()
553
554	,	fSrcTileSize ()
555
556	{
557	if (fRowScale == 0 || fColScale == 0)
558		{
559		 ThrowBadFormat ();
560	}
561
562	if (srcImage.PixelSize  () <= 2 &&
563		dstImage.PixelSize  () <= 2 &&
564		srcImage.PixelRange () == dstImage.PixelRange ())
565		{
566		fSrcPixelType = ttShort;
567		fDstPixelType = ttShort;
568		}
569
570	else
571		{
572		fSrcPixelType = ttFloat;
573		fDstPixelType = ttFloat;
574		}
575
576	fUnitCell = dng_point (8, 8);
577
578	fMaxTileSize.v = Pin_int32 (fUnitCell.v,
579								Round_int32 (fMaxTileSize.v * fRowScale),
580								fMaxTileSize.v);
581
582	fMaxTileSize.h = Pin_int32 (fUnitCell.h,
583								Round_int32 (fMaxTileSize.h * fColScale),
584								fMaxTileSize.h);
585
586	}
587
588/*****************************************************************************/
589
590dng_rect dng_resample_task::SrcArea (const dng_rect &dstArea)
591	{
592
593	int32 offsetV = fWeightsV.Offset ();
594	int32 offsetH = fWeightsH.Offset ();
595
596	uint32 widthV = fWeightsV.Width ();
597	uint32 widthH = fWeightsH.Width ();
598
599	dng_rect srcArea;
600
601	srcArea.t = SafeInt32Add (fRowCoords.Pixel (dstArea.t), offsetV);
602	srcArea.l = SafeInt32Add (fColCoords.Pixel (dstArea.l), offsetH);
603
604	srcArea.b = SafeInt32Add (SafeInt32Add (
605						fRowCoords.Pixel (SafeInt32Sub (dstArea.b, 1)),
606						offsetV),
607					ConvertUint32ToInt32 (widthV));;
608	srcArea.r = SafeInt32Add(SafeInt32Add(
609						fColCoords.Pixel (SafeInt32Sub (dstArea.r, 1)),
610						offsetH),
611					ConvertUint32ToInt32(widthH));;
612
613	return srcArea;
614
615	}
616
617/*****************************************************************************/
618
619dng_point dng_resample_task::SrcTileSize (const dng_point & /* dstTileSize */)
620	{
621
622	return fSrcTileSize;
623
624	}
625
626/*****************************************************************************/
627
628void dng_resample_task::Start (uint32 threadCount,
629							   const dng_point &tileSize,
630							   dng_memory_allocator *allocator,
631							   dng_abort_sniffer *sniffer)
632	{
633
634	// Compute sub-pixel resolution coordinates in the source image for
635	// each row and column of the destination area.
636
637	fRowCoords.Initialize (fSrcBounds.t,
638						   fDstBounds.t,
639						   fSrcBounds.H (),
640						   fDstBounds.H (),
641						   *allocator);
642
643	fColCoords.Initialize (fSrcBounds.l,
644						   fDstBounds.l,
645						   fSrcBounds.W (),
646						   fDstBounds.W (),
647						   *allocator);
648
649	// Compute resampling kernels.
650
651	fWeightsV.Initialize (fRowScale,
652						  fKernel,
653						  *allocator);
654
655	fWeightsH.Initialize (fColScale,
656						  fKernel,
657						  *allocator);
658
659	// Find upper bound on source source tile.
660
661	fSrcTileSize.v = Round_int32 (tileSize.v / fRowScale) + fWeightsV.Width () + 2;
662	fSrcTileSize.h = Round_int32 (tileSize.h / fColScale) + fWeightsH.Width () + 2;
663
664	// Allocate temp buffers.
665
666	uint32 tempBufferSize = 0;
667	if (!RoundUpUint32ToMultiple (fSrcTileSize.h, 8, &tempBufferSize) ||
668		 !SafeUint32Mult (tempBufferSize,
669			 static_cast<uint32> (sizeof (real32)),
670			 &tempBufferSize))
671		{
672
673		ThrowMemoryFull("Arithmetic overflow computing buffer size.");
674
675		}
676
677	for (uint32 threadIndex = 0; threadIndex < threadCount; threadIndex++)
678		{
679
680		fTempBuffer [threadIndex] . Reset (allocator->Allocate (tempBufferSize));
681
682		}
683
684	// Allocate the pixel buffers.
685
686	dng_filter_task::Start (threadCount,
687							tileSize,
688							allocator,
689							sniffer);
690
691	}
692
693/*****************************************************************************/
694
695void dng_resample_task::ProcessArea (uint32 threadIndex,
696								     dng_pixel_buffer &srcBuffer,
697								     dng_pixel_buffer &dstBuffer)
698	{
699
700	dng_rect srcArea = srcBuffer.fArea;
701	dng_rect dstArea = dstBuffer.fArea;
702
703	uint32 srcCols = srcArea.W ();
704	uint32 dstCols = dstArea.W ();
705
706	uint32 widthV = fWeightsV.Width ();
707	uint32 widthH = fWeightsH.Width ();
708
709	int32 offsetV = fWeightsV.Offset ();
710	int32 offsetH = fWeightsH.Offset ();
711
712	uint32 stepH = fWeightsH.Step ();
713
714	const int32 *rowCoords = fRowCoords.Coords (0        );
715	const int32 *colCoords = fColCoords.Coords (dstArea.l);
716
717	if (fSrcPixelType == ttFloat)
718		{
719
720		const real32 *weightsH = fWeightsH.Weights32 (0);
721
722		real32 *tPtr = fTempBuffer [threadIndex]->Buffer_real32 ();
723
724		real32 *ttPtr = tPtr + offsetH - srcArea.l;
725
726		for (int32 dstRow = dstArea.t; dstRow < dstArea.b; dstRow++)
727			{
728
729			int32 rowCoord = rowCoords [dstRow];
730
731			int32 rowFract = rowCoord & kResampleSubsampleMask;
732
733			const real32 *weightsV = fWeightsV.Weights32 (rowFract);
734
735			int32 srcRow = (rowCoord >> kResampleSubsampleBits) + offsetV;
736
737			for (uint32 plane = 0; plane < dstBuffer.fPlanes; plane++)
738				{
739
740				const real32 *sPtr = srcBuffer.ConstPixel_real32 (srcRow,
741																  srcArea.l,
742																  plane);
743
744				DoResampleDown32 (sPtr,
745								  tPtr,
746								  srcCols,
747								  srcBuffer.fRowStep,
748								  weightsV,
749								  widthV);
750
751				real32 *dPtr = dstBuffer.DirtyPixel_real32 (dstRow,
752															dstArea.l,
753															plane);
754
755				DoResampleAcross32 (ttPtr,
756									dPtr,
757									dstCols,
758									colCoords,
759									weightsH,
760									widthH,
761									stepH);
762
763				}
764
765			}
766
767		}
768
769	else
770		{
771
772		const int16 *weightsH = fWeightsH.Weights16 (0);
773
774		uint16 *tPtr = fTempBuffer [threadIndex]->Buffer_uint16 ();
775
776		uint16 *ttPtr = tPtr + offsetH - srcArea.l;
777
778		uint32 pixelRange = fDstImage.PixelRange ();
779
780		for (int32 dstRow = dstArea.t; dstRow < dstArea.b; dstRow++)
781			{
782
783			int32 rowCoord = rowCoords [dstRow];
784
785			int32 rowFract = rowCoord & kResampleSubsampleMask;
786
787			const int16 *weightsV = fWeightsV.Weights16 (rowFract);
788
789			int32 srcRow = (rowCoord >> kResampleSubsampleBits) + offsetV;
790
791			for (uint32 plane = 0; plane < dstBuffer.fPlanes; plane++)
792				{
793
794				const uint16 *sPtr = srcBuffer.ConstPixel_uint16 (srcRow,
795																  srcArea.l,
796																  plane);
797
798				DoResampleDown16 (sPtr,
799								  tPtr,
800								  srcCols,
801								  srcBuffer.fRowStep,
802								  weightsV,
803								  widthV,
804								  pixelRange);
805
806				uint16 *dPtr = dstBuffer.DirtyPixel_uint16 (dstRow,
807															dstArea.l,
808															plane);
809
810				DoResampleAcross16 (ttPtr,
811									dPtr,
812									dstCols,
813									colCoords,
814									weightsH,
815									widthH,
816									stepH,
817									pixelRange);
818
819				}
820
821			}
822
823		}
824
825	}
826
827/*****************************************************************************/
828
829void ResampleImage (dng_host &host,
830					const dng_image &srcImage,
831					dng_image &dstImage,
832					const dng_rect &srcBounds,
833					const dng_rect &dstBounds,
834					const dng_resample_function &kernel)
835	{
836
837	dng_resample_task task (srcImage,
838							dstImage,
839							srcBounds,
840							dstBounds,
841							kernel);
842
843	host.PerformAreaTask (task,
844						  dstBounds);
845
846	}
847
848/*****************************************************************************/
849