comparison Framework/ColorSpaces.cpp @ 330:c42083d50ddf

Added support for DICOM tag "Recommended Absent Pixel CIELab" (0048,0015)
author Sebastien Jodogne <s.jodogne@gmail.com>
date Fri, 18 Oct 2024 13:08:55 +0200
parents
children
comparison
equal deleted inserted replaced
329:ae2d769215d2 330:c42083d50ddf
1 /**
2 * Orthanc - A Lightweight, RESTful DICOM Store
3 * Copyright (C) 2012-2016 Sebastien Jodogne, Medical Physics
4 * Department, University Hospital of Liege, Belgium
5 * Copyright (C) 2017-2023 Osimis S.A., Belgium
6 * Copyright (C) 2024-2024 Orthanc Team SRL, Belgium
7 * Copyright (C) 2021-2024 Sebastien Jodogne, ICTEAM UCLouvain, Belgium
8 *
9 * This program is free software: you can redistribute it and/or
10 * modify it under the terms of the GNU Affero General Public License
11 * as published by the Free Software Foundation, either version 3 of
12 * the License, or (at your option) any later version.
13 *
14 * This program is distributed in the hope that it will be useful, but
15 * WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 * Affero General Public License for more details.
18 *
19 * You should have received a copy of the GNU Affero General Public License
20 * along with this program. If not, see <http://www.gnu.org/licenses/>.
21 **/
22
23
24 #include "ColorSpaces.h"
25
26 #include <SerializationToolbox.h>
27 #include <Toolbox.h>
28
29 #include <boost/math/special_functions/round.hpp>
30
31
32 namespace OrthancWSI
33 {
34 RGBColor::RGBColor(const sRGBColor& srgb)
35 {
36 if (srgb.GetR() < 0)
37 {
38 r_ = 0;
39 }
40 else if (srgb.GetR() >= 1)
41 {
42 r_ = 255;
43 }
44 else
45 {
46 r_ = boost::math::iround(srgb.GetR() * 255.0f);
47 }
48
49 if (srgb.GetG() < 0)
50 {
51 g_ = 0;
52 }
53 else if (srgb.GetG() >= 1)
54 {
55 g_ = 255;
56 }
57 else
58 {
59 g_ = boost::math::iround(srgb.GetG() * 255.0f);
60 }
61
62 if (srgb.GetB() < 0)
63 {
64 b_ = 0;
65 }
66 else if (srgb.GetB() >= 1)
67 {
68 b_ = 255;
69 }
70 else
71 {
72 b_ = boost::math::iround(srgb.GetB() * 255.0f);
73 }
74 }
75
76
77 sRGBColor::sRGBColor(const RGBColor& rgb)
78 {
79 r_ = static_cast<float>(rgb.GetR()) / 255.0f;
80 g_ = static_cast<float>(rgb.GetG()) / 255.0f;
81 b_ = static_cast<float>(rgb.GetB()) / 255.0f;
82 }
83
84
85 static float ApplyGammaXYZ(float value)
86 {
87 // https://www.image-engineering.de/library/technotes/958-how-to-convert-between-srgb-and-ciexyz
88 if (value <= 0.0031308f)
89 {
90 return value * 12.92f;
91 }
92 else
93 {
94 return 1.055f * powf(value, 1.0f / 2.4f) - 0.055f;
95 }
96 }
97
98
99 sRGBColor::sRGBColor(const XYZColor& xyz)
100 {
101 // https://en.wikipedia.org/wiki/SRGB#From_CIE_XYZ_to_sRGB
102 // https://www.image-engineering.de/library/technotes/958-how-to-convert-between-srgb-and-ciexyz
103 const float sr = 3.2404542f * xyz.GetX() - 1.5371385f * xyz.GetY() - 0.4985314f * xyz.GetZ();
104 const float sg = -0.9692660f * xyz.GetX() + 1.8760108f * xyz.GetY() + 0.0415560f * xyz.GetZ();
105 const float sb = 0.0556434f * xyz.GetX() - 0.2040259f * xyz.GetY() + 1.0572252f * xyz.GetZ();
106
107 r_ = ApplyGammaXYZ(sr);
108 g_ = ApplyGammaXYZ(sg);
109 b_ = ApplyGammaXYZ(sb);
110 }
111
112
113 static float LinearizeXYZ(float value)
114 {
115 // https://www.image-engineering.de/library/technotes/958-how-to-convert-between-srgb-and-ciexyz
116 if (value <= 0.04045f)
117 {
118 return value / 12.92f;
119 }
120 else
121 {
122 return powf((value + 0.055f) / 1.055f, 2.4f);
123 }
124 }
125
126
127 XYZColor::XYZColor(const sRGBColor& srgb)
128 {
129 // https://en.wikipedia.org/wiki/SRGB#From_sRGB_to_CIE_XYZ
130 // https://www.image-engineering.de/library/technotes/958-how-to-convert-between-srgb-and-ciexyz
131 const float linearizedR = LinearizeXYZ(srgb.GetR());
132 const float linearizedG = LinearizeXYZ(srgb.GetG());
133 const float linearizedB = LinearizeXYZ(srgb.GetB());
134
135 x_ = 0.4124564f * linearizedR + 0.3575761f * linearizedG + 0.1804375f * linearizedB;
136 y_ = 0.2126729f * linearizedR + 0.7151522f * linearizedG + 0.0721750f * linearizedB;
137 z_ = 0.0193339f * linearizedR + 0.1191920f * linearizedG + 0.9503041f * linearizedB;
138 }
139
140
141 static const float LAB_DELTA = 6.0f / 29.0f;
142
143 static float LABf_inv(float t)
144 {
145 if (t > LAB_DELTA)
146 {
147 return powf(t, 3.0f);
148 }
149 else
150 {
151 return 3.0f * LAB_DELTA * LAB_DELTA * (t - 4.0f / 29.0f);
152 }
153 }
154
155
156 // Those correspond to Standard Illuminant D65
157 // https://en.wikipedia.org/wiki/CIELAB_color_space#From_CIEXYZ_to_CIELAB
158 static const float X_N = 95.0489f;
159 static const float Y_N = 100.0f;
160 static const float Z_N = 108.8840f;
161
162 XYZColor::XYZColor(const LABColor& lab)
163 {
164 // https://en.wikipedia.org/wiki/CIELAB_color_space#From_CIELAB_to_CIEXYZ
165 const float shared = (lab.GetL() + 16.0f) / 116.0f;
166
167 x_ = X_N * LABf_inv(shared + lab.GetA() / 500.0f) / 100.0f;
168 y_ = Y_N * LABf_inv(shared) / 100.0f;
169 z_ = Z_N * LABf_inv(shared - lab.GetB() / 200.0f) / 100.0f;
170 }
171
172
173 static float LABf(float t)
174 {
175 if (t > powf(LAB_DELTA, 3.0f))
176 {
177 return powf(t, 1.0f / 3.0f);
178 }
179 else
180 {
181 return t / 3.0f * powf(LAB_DELTA, -2.0f) + 4.0f / 29.0f;
182 }
183 }
184
185
186 LABColor::LABColor(const XYZColor& xyz)
187 {
188 // https://en.wikipedia.org/wiki/CIELAB_color_space#From_CIEXYZ_to_CIELAB
189
190 const float fx = LABf(xyz.GetX() * 100.0f / X_N);
191 const float fy = LABf(xyz.GetY() * 100.0f / Y_N);
192 const float fz = LABf(xyz.GetZ() * 100.0f / Z_N);
193
194 l_ = 116.0f * fy - 16.0f;
195 a_ = 500.0f * (fx - fy);
196 b_ = 200.0f * (fy - fz);
197 }
198
199
200 static uint16_t EncodeUint16(float value,
201 float minValue,
202 float maxValue)
203 {
204 if (value <= minValue)
205 {
206 return 0;
207 }
208 else if (value >= maxValue)
209 {
210 return 0xffff;
211 }
212 else
213 {
214 float lambda = (value - minValue) / (maxValue - minValue);
215 assert(lambda >= 0 && lambda <= 1);
216 return static_cast<uint16_t>(boost::math::iround(lambda * static_cast<float>(0xffff)));
217 }
218 }
219
220
221 void LABColor::EncodeDicomRecommendedAbsentPixelCIELab(uint16_t target[3]) const
222 {
223 /**
224 * "An L value linearly scaled to 16 bits, such that 0x0000
225 * corresponds to an L of 0.0, and 0xFFFF corresponds to an L of
226 * 100.0."
227 **/
228 target[0] = EncodeUint16(GetL(), 0.0f, 100.0f);
229
230 /**
231 * "An a* then a b* value, each linearly scaled to 16 bits and
232 * offset to an unsigned range, such that 0x0000 corresponds to an
233 * a* or b* of -128.0, 0x8080 corresponds to an a* or b* of 0.0
234 * and 0xFFFF corresponds to an a* or b* of 127.0"
235 **/
236 target[1] = EncodeUint16(GetA(), -128.0f, 127.0f);
237 target[2] = EncodeUint16(GetB(), -128.0f, 127.0f);
238 }
239
240
241 LABColor LABColor::DecodeDicomRecommendedAbsentPixelCIELab(uint16_t l,
242 uint16_t a,
243 uint16_t b)
244 {
245 return LABColor(static_cast<float>(l) / static_cast<float>(0xffff) * 100.0f,
246 -128.0f + static_cast<float>(a) / static_cast<float>(0xffff) * 255.0f,
247 -128.0f + static_cast<float>(b) / static_cast<float>(0xffff) * 255.0f);
248 }
249
250
251 bool LABColor::DecodeDicomRecommendedAbsentPixelCIELab(LABColor& target,
252 const std::string& tag)
253 {
254 std::vector<std::string> channels;
255 Orthanc::Toolbox::TokenizeString(channels, tag, '\\');
256
257 unsigned int l, a, b;
258 if (channels.size() == 3 &&
259 Orthanc::SerializationToolbox::ParseUnsignedInteger32(l, channels[0]) &&
260 Orthanc::SerializationToolbox::ParseUnsignedInteger32(a, channels[1]) &&
261 Orthanc::SerializationToolbox::ParseUnsignedInteger32(b, channels[2]) &&
262 l <= 0xffffu &&
263 a <= 0xffffu &&
264 b <= 0xffffu)
265 {
266 target = LABColor::DecodeDicomRecommendedAbsentPixelCIELab(l, a, b);
267 return true;
268 }
269 else
270 {
271 return false;
272 }
273 }
274 }