1 /*****************************************************************************/
2 // Copyright 2006-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_color_spec.cpp#1 $ */
10 /* $DateTime: 2012/05/30 13:28:51 $ */
11 /* $Change: 832332 $ */
12 /* $Author: tknoll $ */
13 
14 #include "dng_color_spec.h"
15 
16 #include "dng_assertions.h"
17 #include "dng_camera_profile.h"
18 #include "dng_exceptions.h"
19 #include "dng_matrix.h"
20 #include "dng_negative.h"
21 #include "dng_temperature.h"
22 #include "dng_utils.h"
23 #include "dng_xy_coord.h"
24 
25 /*****************************************************************************/
26 
MapWhiteMatrix(const dng_xy_coord & white1,const dng_xy_coord & white2)27 dng_matrix_3by3 MapWhiteMatrix (const dng_xy_coord &white1,
28 						        const dng_xy_coord &white2)
29 	{
30 
31 	// Use the linearized Bradford adaptation matrix.
32 
33 	dng_matrix_3by3 Mb ( 0.8951,  0.2664, -0.1614,
34 		 		        -0.7502,  1.7135,  0.0367,
35 		  			     0.0389, -0.0685,  1.0296);
36 
37 	dng_vector_3 w1 = Mb * XYtoXYZ (white1);
38 	dng_vector_3 w2 = Mb * XYtoXYZ (white2);
39 
40 	// Negative white coordinates are kind of meaningless.
41 
42 	w1 [0] = Max_real64 (w1 [0], 0.0);
43 	w1 [1] = Max_real64 (w1 [1], 0.0);
44 	w1 [2] = Max_real64 (w1 [2], 0.0);
45 
46 	w2 [0] = Max_real64 (w2 [0], 0.0);
47 	w2 [1] = Max_real64 (w2 [1], 0.0);
48 	w2 [2] = Max_real64 (w2 [2], 0.0);
49 
50 	// Limit scaling to something reasonable.
51 
52 	dng_matrix_3by3 A;
53 
54 	A [0] [0] = Pin_real64 (0.1, w1 [0] > 0.0 ? w2 [0] / w1 [0] : 10.0, 10.0);
55 	A [1] [1] = Pin_real64 (0.1, w1 [1] > 0.0 ? w2 [1] / w1 [1] : 10.0, 10.0);
56 	A [2] [2] = Pin_real64 (0.1, w1 [2] > 0.0 ? w2 [2] / w1 [2] : 10.0, 10.0);
57 
58 	dng_matrix_3by3 B = Invert (Mb) * A * Mb;
59 
60 	return B;
61 
62 	}
63 
64 /******************************************************************************/
65 
dng_color_spec(const dng_negative & negative,const dng_camera_profile * profile)66 dng_color_spec::dng_color_spec (const dng_negative &negative,
67 							    const dng_camera_profile *profile)
68 
69 	:	fChannels (negative.ColorChannels ())
70 
71 	,	fTemperature1 (0.0)
72 	,	fTemperature2 (0.0)
73 
74 	,	fColorMatrix1 ()
75 	,	fColorMatrix2 ()
76 
77 	,	fForwardMatrix1 ()
78 	,	fForwardMatrix2 ()
79 
80 	,	fReductionMatrix1 ()
81 	,	fReductionMatrix2 ()
82 
83 	,	fCameraCalibration1 ()
84 	,	fCameraCalibration2 ()
85 
86 	,	fAnalogBalance ()
87 
88 	,	fWhiteXY ()
89 
90 	,	fCameraWhite ()
91 	,	fCameraToPCS ()
92 
93 	,	fPCStoCamera ()
94 
95 	{
96 
97 	if (fChannels > 1)
98 		{
99 
100 		if (!profile || !profile->IsValid (fChannels))
101 			{
102 			ThrowBadFormat ();
103 			}
104 
105 		if (profile->WasStubbed ())
106 			{
107 			ThrowProgramError ("Using stubbed profile");
108 			}
109 
110 		fTemperature1 = profile->CalibrationTemperature1 ();
111 		fTemperature2 = profile->CalibrationTemperature2 ();
112 
113 		fColorMatrix1 = profile->ColorMatrix1 ();
114 		fColorMatrix2 = profile->ColorMatrix2 ();
115 
116 		fForwardMatrix1 = profile->ForwardMatrix1 ();
117 		fForwardMatrix2 = profile->ForwardMatrix2 ();
118 
119 		fReductionMatrix1 = profile->ReductionMatrix1 ();
120 		fReductionMatrix2 = profile->ReductionMatrix2 ();
121 
122 		fCameraCalibration1.SetIdentity (fChannels);
123 		fCameraCalibration2.SetIdentity (fChannels);
124 
125 		if (negative. CameraCalibrationSignature () ==
126 			profile->ProfileCalibrationSignature ())
127 			{
128 
129 			if (negative.CameraCalibration1 ().Rows () == fChannels &&
130 				negative.CameraCalibration1 ().Cols () == fChannels)
131 				{
132 
133 				fCameraCalibration1 = negative.CameraCalibration1 ();
134 
135 				}
136 
137 			if (negative.CameraCalibration2 ().Rows () == fChannels &&
138 				negative.CameraCalibration2 ().Cols () == fChannels)
139 				{
140 
141 				fCameraCalibration2 = negative.CameraCalibration2 ();
142 
143 				}
144 
145 			}
146 
147 		fAnalogBalance = dng_matrix (fChannels, fChannels);
148 
149 		for (uint32 j = 0; j < fChannels; j++)
150 			{
151 
152 			fAnalogBalance [j] [j] = negative.AnalogBalance (j);
153 
154 			}
155 
156 		dng_camera_profile::NormalizeForwardMatrix (fForwardMatrix1);
157 
158 		fColorMatrix1 = fAnalogBalance * fCameraCalibration1 * fColorMatrix1;
159 
160 		if (!profile->HasColorMatrix2 () ||
161 				fTemperature1 <= 0.0 ||
162 				fTemperature2 <= 0.0 ||
163 				fTemperature1 == fTemperature2)
164 			{
165 
166 			fTemperature1 = 5000.0;
167 			fTemperature2 = 5000.0;
168 
169 			fColorMatrix2       = fColorMatrix1;
170 			fForwardMatrix2     = fForwardMatrix1;
171 			fReductionMatrix2   = fReductionMatrix1;
172 			fCameraCalibration2 = fCameraCalibration1;
173 
174 			}
175 
176 		else
177 			{
178 
179 			dng_camera_profile::NormalizeForwardMatrix (fForwardMatrix2);
180 
181 			fColorMatrix2 = fAnalogBalance * fCameraCalibration2 * fColorMatrix2;
182 
183 			// Swap values if temperatures are out of order.
184 
185 			if (fTemperature1 > fTemperature2)
186 				{
187 
188 				real64 temp   = fTemperature1;
189 				fTemperature1 = fTemperature2;
190 				fTemperature2 = temp;
191 
192 				dng_matrix T  = fColorMatrix1;
193 				fColorMatrix1 = fColorMatrix2;
194 				fColorMatrix2 = T;
195 
196 				T               = fForwardMatrix1;
197 				fForwardMatrix1 = fForwardMatrix2;
198 				fForwardMatrix2 = T;
199 
200 				T                 = fReductionMatrix1;
201 				fReductionMatrix1 = fReductionMatrix2;
202 				fReductionMatrix2 = T;
203 
204 				T                   = fCameraCalibration1;
205 				fCameraCalibration1 = fCameraCalibration2;
206 				fCameraCalibration2 = T;
207 
208 				}
209 
210 			}
211 
212 		}
213 
214 	}
215 
216 /*****************************************************************************/
217 
FindXYZtoCamera(const dng_xy_coord & white,dng_matrix * forwardMatrix,dng_matrix * reductionMatrix,dng_matrix * cameraCalibration)218 dng_matrix dng_color_spec::FindXYZtoCamera (const dng_xy_coord &white,
219 											dng_matrix *forwardMatrix,
220 									        dng_matrix *reductionMatrix,
221 											dng_matrix *cameraCalibration)
222 	{
223 
224 	// Convert to temperature/offset space.
225 
226 	dng_temperature td (white);
227 
228 	// Find fraction to weight the first calibration.
229 
230 	real64 g;
231 
232 	if (td.Temperature () <= fTemperature1)
233 		g = 1.0;
234 
235 	else if (td.Temperature () >= fTemperature2)
236 		g = 0.0;
237 
238 	else
239 		{
240 
241 		real64 invT = 1.0 / td.Temperature ();
242 
243 		g = (invT                  - (1.0 / fTemperature2)) /
244 		    ((1.0 / fTemperature1) - (1.0 / fTemperature2));
245 
246 		}
247 
248 	// Interpolate the color matrix.
249 
250 	dng_matrix colorMatrix;
251 
252 	if (g >= 1.0)
253 		colorMatrix = fColorMatrix1;
254 
255 	else if (g <= 0.0)
256 		colorMatrix = fColorMatrix2;
257 
258 	else
259 		colorMatrix = (g      ) * fColorMatrix1 +
260 					  (1.0 - g) * fColorMatrix2;
261 
262 	// Interpolate forward matrix, if any.
263 
264 	if (forwardMatrix)
265 		{
266 
267 		bool has1 = fForwardMatrix1.NotEmpty ();
268 		bool has2 = fForwardMatrix2.NotEmpty ();
269 
270 		if (has1 && has2)
271 			{
272 
273 			if (g >= 1.0)
274 				*forwardMatrix = fForwardMatrix1;
275 
276 			else if (g <= 0.0)
277 				*forwardMatrix = fForwardMatrix2;
278 
279 			else
280 				*forwardMatrix = (g      ) * fForwardMatrix1 +
281 								 (1.0 - g) * fForwardMatrix2;
282 
283 			}
284 
285 		else if (has1)
286 			{
287 
288 			*forwardMatrix = fForwardMatrix1;
289 
290 			}
291 
292 		else if (has2)
293 			{
294 
295 			*forwardMatrix = fForwardMatrix2;
296 
297 			}
298 
299 		else
300 			{
301 
302 			forwardMatrix->Clear ();
303 
304 			}
305 
306 		}
307 
308 	// Interpolate reduction matrix, if any.
309 
310 	if (reductionMatrix)
311 		{
312 
313 		bool has1 = fReductionMatrix1.NotEmpty ();
314 		bool has2 = fReductionMatrix2.NotEmpty ();
315 
316 		if (has1 && has2)
317 			{
318 
319 			if (g >= 1.0)
320 				*reductionMatrix = fReductionMatrix1;
321 
322 			else if (g <= 0.0)
323 				*reductionMatrix = fReductionMatrix2;
324 
325 			else
326 				*reductionMatrix = (g      ) * fReductionMatrix1 +
327 								   (1.0 - g) * fReductionMatrix2;
328 
329 			}
330 
331 		else if (has1)
332 			{
333 
334 			*reductionMatrix = fReductionMatrix1;
335 
336 			}
337 
338 		else if (has2)
339 			{
340 
341 			*reductionMatrix = fReductionMatrix2;
342 
343 			}
344 
345 		else
346 			{
347 
348 			reductionMatrix->Clear ();
349 
350 			}
351 
352 		}
353 
354 	// Interpolate camera calibration matrix.
355 
356 	if (cameraCalibration)
357 		{
358 
359 		if (g >= 1.0)
360 			*cameraCalibration = fCameraCalibration1;
361 
362 		else if (g <= 0.0)
363 			*cameraCalibration = fCameraCalibration2;
364 
365 		else
366 			*cameraCalibration = (g      ) * fCameraCalibration1 +
367 							     (1.0 - g) * fCameraCalibration2;
368 
369 		}
370 
371 	// Return the interpolated color matrix.
372 
373 	return colorMatrix;
374 
375 	}
376 
377 /*****************************************************************************/
378 
SetWhiteXY(const dng_xy_coord & white)379 void dng_color_spec::SetWhiteXY (const dng_xy_coord &white)
380 	{
381 
382 	fWhiteXY = white;
383 
384 	// Deal with monochrome cameras.
385 
386 	if (fChannels == 1)
387 		{
388 
389 		fCameraWhite.SetIdentity (1);
390 
391 		fCameraToPCS = PCStoXYZ ().AsColumn ();
392 
393 		return;
394 
395 		}
396 
397 	// Interpolate an matric values for this white point.
398 
399 	dng_matrix colorMatrix;
400 	dng_matrix forwardMatrix;
401 	dng_matrix reductionMatrix;
402 	dng_matrix cameraCalibration;
403 
404 	colorMatrix = FindXYZtoCamera (fWhiteXY,
405 								   &forwardMatrix,
406 								   &reductionMatrix,
407 								   &cameraCalibration);
408 
409 	// Find the camera white values.
410 
411 	fCameraWhite = colorMatrix * XYtoXYZ (fWhiteXY);
412 
413 	real64 cameraWhiteMaxEntry = MaxEntry (fCameraWhite);
414 	if (cameraWhiteMaxEntry == 0)
415 		{
416 		 ThrowBadFormat ();
417 		}
418 	real64 whiteScale = 1.0 / cameraWhiteMaxEntry;
419 
420 	for (uint32 j = 0; j < fChannels; j++)
421 		{
422 
423 		// We don't support non-positive values for camera neutral values.
424 
425 		fCameraWhite [j] = Pin_real64 (0.001,
426 									   whiteScale * fCameraWhite [j],
427 									   1.0);
428 
429 		}
430 
431 	// Find PCS to Camera transform. Scale matrix so PCS white can just be
432 	// reached when the first camera channel saturates
433 
434 	fPCStoCamera = colorMatrix * MapWhiteMatrix (PCStoXY (), fWhiteXY);
435 
436 	real64 scale = MaxEntry (fPCStoCamera * PCStoXYZ ());
437 
438 	if (scale == 0)
439 		{
440 		 ThrowBadFormat ();
441 		}
442 	fPCStoCamera = (1.0 / scale) * fPCStoCamera;
443 
444 	// If we have a forward matrix, then just use that.
445 
446 	if (forwardMatrix.NotEmpty ())
447 		{
448 
449 		dng_matrix individualToReference = Invert (fAnalogBalance * cameraCalibration);
450 
451 		dng_vector refCameraWhite = individualToReference * fCameraWhite;
452 
453 		fCameraToPCS = forwardMatrix *
454 					   Invert (refCameraWhite.AsDiagonal ()) *
455 					   individualToReference;
456 
457 		}
458 
459 	// Else we need to use the adapt in XYZ method.
460 
461 	else
462 		{
463 
464 		// Invert this PCS to camera matrix.  Note that if there are more than three
465 		// camera channels, this inversion is non-unique.
466 
467 		fCameraToPCS = Invert (fPCStoCamera, reductionMatrix);
468 
469 		}
470 
471 	}
472 
473 /*****************************************************************************/
474 
WhiteXY() const475 const dng_xy_coord & dng_color_spec::WhiteXY () const
476 	{
477 
478 	DNG_ASSERT (fWhiteXY.IsValid (), "Using invalid WhiteXY");
479 
480 	return fWhiteXY;
481 
482 	}
483 
484 /*****************************************************************************/
485 
CameraWhite() const486 const dng_vector & dng_color_spec::CameraWhite () const
487 	{
488 
489 	DNG_ASSERT (fCameraWhite.NotEmpty (), "Using invalid CameraWhite");
490 
491 	return fCameraWhite;
492 
493 	}
494 
495 /*****************************************************************************/
496 
CameraToPCS() const497 const dng_matrix & dng_color_spec::CameraToPCS () const
498 	{
499 
500 	DNG_ASSERT (fCameraToPCS.NotEmpty (), "Using invalid CameraToPCS");
501 
502 	return fCameraToPCS;
503 
504 	}
505 
506 /*****************************************************************************/
507 
PCStoCamera() const508 const dng_matrix & dng_color_spec::PCStoCamera () const
509 	{
510 
511 	DNG_ASSERT (fPCStoCamera.NotEmpty (), "Using invalid PCStoCamera");
512 
513 	return fPCStoCamera;
514 
515 	}
516 
517 /*****************************************************************************/
518 
NeutralToXY(const dng_vector & neutral)519 dng_xy_coord dng_color_spec::NeutralToXY (const dng_vector &neutral)
520 	{
521 
522 	const uint32 kMaxPasses = 30;
523 
524 	if (fChannels == 1)
525 		{
526 
527 		return PCStoXY ();
528 
529 		}
530 
531 	dng_xy_coord last = D50_xy_coord ();
532 
533 	for (uint32 pass = 0; pass < kMaxPasses; pass++)
534 		{
535 
536 		dng_matrix xyzToCamera = FindXYZtoCamera (last);
537 
538 		dng_xy_coord next = XYZtoXY (Invert (xyzToCamera) * neutral);
539 
540 		if (Abs_real64 (next.x - last.x) +
541 			Abs_real64 (next.y - last.y) < 0.0000001)
542 			{
543 
544 			return next;
545 
546 			}
547 
548 		// If we reach the limit without converging, we are most likely
549 		// in a two value oscillation.  So take the average of the last
550 		// two estimates and give up.
551 
552 		if (pass == kMaxPasses - 1)
553 			{
554 
555 			next.x = (last.x + next.x) * 0.5;
556 			next.y = (last.y + next.y) * 0.5;
557 
558 			}
559 
560 		last = next;
561 
562 		}
563 
564 	return last;
565 
566 	}
567 
568 /*****************************************************************************/
569