1 /*=========================================================================
2
3 Program: GDCM (Grassroots DICOM). A DICOM library
4
5 Copyright (c) 2006-2011 Mathieu Malaterre
6 All rights reserved.
7 See Copyright.txt or http://gdcm.sourceforge.net/Copyright.html for details.
8
9 This software is distributed WITHOUT ANY WARRANTY; without even
10 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
11 PURPOSE. See the above copyright notice for more information.
12
13 =========================================================================*/
14 #include "gdcmSurfaceHelper.h"
15
16 #include <cmath>
17
18 namespace gdcm
19 {
20
RGBToXYZ(const std::vector<float> & RGB)21 std::vector< float > SurfaceHelper::RGBToXYZ(const std::vector<float> & RGB)
22 {
23 std::vector< float > XYZ(3);
24 float tmp[3];
25
26 tmp[0] = RGB[0];
27 tmp[1] = RGB[1];
28 tmp[2] = RGB[2];
29
30 const float A = 1.0f / 12.92f;
31 const float B = 1.0f / 1.055f;
32
33 if ( tmp[0] > 0.04045f ) tmp[0] = powf( ( tmp[0] + 0.055f ) * B, 2.4f);
34 else tmp[0] *= A;
35 if ( tmp[1] > 0.04045f ) tmp[1] = powf( ( tmp[1] + 0.055f ) * B, 2.4f);
36 else tmp[1] *= A;
37 if ( tmp[2] > 0.04045f ) tmp[2] = powf( ( tmp[2] + 0.055f ) * B, 2.4f);
38 else tmp[2] *= A;
39
40 tmp[0] *= 100;
41 tmp[1] *= 100;
42 tmp[2] *= 100;
43
44 //Observer. = 2°, Illuminant = D65
45 XYZ[0] = (tmp[0] * 0.4124f + tmp[1] * 0.3576f + tmp[2] * 0.1805f);
46 XYZ[1] = (tmp[0] * 0.2126f + tmp[1] * 0.7152f + tmp[2] * 0.0722f);
47 XYZ[2] = (tmp[0] * 0.0193f + tmp[1] * 0.1192f + tmp[2] * 0.9505f);
48
49 return XYZ;
50 }
51
XYZToRGB(const std::vector<float> & XYZ)52 std::vector< float > SurfaceHelper::XYZToRGB(const std::vector<float> & XYZ)
53 {
54 std::vector< float > RGB(3);
55 float tmp[3];
56
57 tmp[0] = XYZ[0];
58 tmp[1] = XYZ[1];
59 tmp[2] = XYZ[2];
60
61 // Divide by 100
62 tmp[0] *= 0.01f; //X from 0 to 95.047 (Observer = 2°, Illuminant = D65)
63 tmp[1] *= 0.01f; //Y from 0 to 100.000
64 tmp[2] *= 0.01f; //Z from 0 to 108.883
65
66 RGB[0] = (tmp[0] * 3.2406f + tmp[1] * -1.5372f + tmp[2] * -0.4986f);
67 RGB[1] = (tmp[0] * -0.9689f + tmp[1] * 1.8758f + tmp[2] * 0.0415f);
68 RGB[2] = (tmp[0] * 0.0557f + tmp[1] * -0.2040f + tmp[2] * 1.0570f);
69
70 const float A = 1.0f / 2.4f;
71
72 if ( RGB[0] > 0.0031308f ) RGB[0] = 1.055f * powf( RGB[0], A) - 0.055f;
73 else RGB[0] *= 12.92f;
74 if ( RGB[1] > 0.0031308f ) RGB[1] = 1.055f * powf( RGB[1], A) - 0.055f;
75 else RGB[1] *= 12.92f;
76 if ( RGB[2] > 0.0031308f ) RGB[2] = 1.055f * powf( RGB[2], A) - 0.055f;
77 else RGB[2] *= 12.92f;
78
79 return RGB;
80 }
81
XYZToCIELab(const std::vector<float> & XYZ)82 std::vector< float > SurfaceHelper::XYZToCIELab(const std::vector<float> & XYZ)
83 {
84 std::vector< float > CIELab(3);
85 float tmp[3];
86
87 tmp[0] = (XYZ[0] / 95.047f); //ref_X = 95.047 Observer= 2°, Illuminant= D65
88 tmp[1] = (XYZ[1] * 0.01f); //ref_Y = 100.000
89 tmp[2] = (XYZ[2] / 108.883f); //ref_Z = 108.883
90
91 const float A = 1.0f / 3.0f;
92 const float B = 16.0f / 116.0f;
93
94 if ( tmp[0] > 0.008856f ) tmp[0] = powf(tmp[0], A);
95 else tmp[0] = (7.787f * tmp[0] + B);
96 if ( tmp[1] > 0.008856f ) tmp[1] = powf(tmp[1], A);
97 else tmp[1] = (7.787f * tmp[1] + B);
98 if ( tmp[2] > 0.008856f ) tmp[2] = powf(tmp[2], A);
99 else tmp[2] = (7.787f * tmp[2] + B);
100
101 CIELab[0] = ( 116 * tmp[1] ) - 16;
102 CIELab[1] = 500 * ( tmp[0] - tmp[1] );
103 CIELab[2] = 200 * ( tmp[1] - tmp[2] );
104
105 return CIELab;
106 }
107
CIELabToXYZ(const std::vector<float> & CIELab)108 std::vector< float > SurfaceHelper::CIELabToXYZ(const std::vector<float> & CIELab)
109 {
110 std::vector< float > XYZ(3);
111 float tmp[3];
112
113 tmp[1] = (( CIELab[0] + 16 ) / 116.0f);
114 tmp[0] = (CIELab[1] * 0.002f + tmp[1]);
115 tmp[2] = (tmp[1] - CIELab[2] * 0.005f);
116
117 // Compute t
118 const float A = tmp[0]*tmp[0]*tmp[0];
119 const float B = tmp[1]*tmp[1]*tmp[1];
120 const float C = tmp[2]*tmp[2]*tmp[2];
121
122 const float D = 16.0f / 116.0f;
123
124 // Compute f(t)
125 if ( B > 0.008856f) tmp[1] = B;
126 else tmp[1] = ( tmp[1] - D ) / 7.787f;
127 if ( A > 0.008856f) tmp[0] = A;
128 else tmp[0] = ( tmp[0] - D ) / 7.787f;
129 if ( C > 0.008856f) tmp[2] = C;
130 else tmp[2] = ( tmp[2] - D ) / 7.787f;
131
132 XYZ[0] = (tmp[0] * 95.047f); //ref_X = 95.047 Observer= 2°, Illuminant= D65
133 XYZ[1] = (tmp[1] * 100.0f); //ref_Y = 100.000
134 XYZ[2] = (tmp[2] * 108.883f); //ref_Z = 108.883
135
136 return XYZ;
137 }
138
139 }
140