/src/dcmtk/dcmimgle/libsrc/digsdlut.cc
Line | Count | Source |
1 | | /* |
2 | | * |
3 | | * Copyright (C) 1999-2025, OFFIS e.V. |
4 | | * All rights reserved. See COPYRIGHT file for details. |
5 | | * |
6 | | * This software and supporting documentation were developed by |
7 | | * |
8 | | * OFFIS e.V. |
9 | | * R&D Division Health |
10 | | * Escherweg 2 |
11 | | * D-26121 Oldenburg, Germany |
12 | | * |
13 | | * |
14 | | * Module: dcmimgle |
15 | | * |
16 | | * Author: Joerg Riesmeier |
17 | | * |
18 | | * Purpose: DicomGSDFLUT (Source) |
19 | | * |
20 | | */ |
21 | | |
22 | | |
23 | | #include "dcmtk/config/osconfig.h" |
24 | | |
25 | | #include "dcmtk/dcmimgle/digsdlut.h" |
26 | | #include "dcmtk/dcmimgle/displint.h" |
27 | | |
28 | | #include <cmath> |
29 | | |
30 | | /*----------------* |
31 | | * constructors * |
32 | | *----------------*/ |
33 | | |
34 | | DiGSDFLUT::DiGSDFLUT(const unsigned long count, |
35 | | const Uint16 max, |
36 | | const Uint16 *ddl_tab, |
37 | | const double *val_tab, |
38 | | const unsigned long ddl_cnt, |
39 | | const double *gsdf_tab, |
40 | | const double *gsdf_spl, |
41 | | const unsigned int gsdf_cnt, |
42 | | const double jnd_min, |
43 | | const double jnd_max, |
44 | | const double lum_min, |
45 | | const double lum_max, |
46 | | const double amb, |
47 | | const double illum, |
48 | | const OFBool inverse, |
49 | | STD_NAMESPACE ostream *stream, |
50 | | const OFBool printMode) |
51 | 0 | : DiDisplayLUT(count, max, amb, illum) |
52 | 0 | { |
53 | 0 | if ((Count > 0) && (Bits > 0)) |
54 | 0 | { |
55 | 0 | DCMIMGLE_DEBUG("new GSDF LUT with " << Bits << " bits output and " << Count << " entries created"); |
56 | 0 | if (jnd_min >= jnd_max) |
57 | 0 | { |
58 | 0 | DCMIMGLE_ERROR("invalid JND range for GSDF LUT creation (" << jnd_min << " - " << jnd_max << ")"); |
59 | 0 | } |
60 | | /* create the lookup table */ |
61 | 0 | Valid = createLUT(ddl_tab, val_tab, ddl_cnt, gsdf_tab, gsdf_spl, gsdf_cnt, |
62 | 0 | jnd_min, jnd_max, lum_min, lum_max, inverse, stream, printMode); |
63 | 0 | } |
64 | 0 | } |
65 | | |
66 | | |
67 | | /*--------------* |
68 | | * destructor * |
69 | | *--------------*/ |
70 | | |
71 | | DiGSDFLUT::~DiGSDFLUT() |
72 | 0 | { |
73 | 0 | } |
74 | | |
75 | | |
76 | | /********************************************************************/ |
77 | | |
78 | | |
79 | | int DiGSDFLUT::createLUT(const Uint16 *ddl_tab, |
80 | | const double *val_tab, |
81 | | const unsigned long ddl_cnt, |
82 | | const double *gsdf_tab, |
83 | | const double *gsdf_spl, |
84 | | const unsigned int gsdf_cnt, |
85 | | const double jnd_min, |
86 | | const double jnd_max, |
87 | | const double lum_min, |
88 | | const double lum_max, |
89 | | const OFBool inverse, |
90 | | STD_NAMESPACE ostream *stream, |
91 | | const OFBool printMode) |
92 | 0 | { |
93 | | /* check for valid parameters */ |
94 | 0 | if ((ddl_tab != NULL) && (val_tab != NULL) && (ddl_cnt > 0) && (gsdf_tab != NULL) && (gsdf_spl != NULL) && |
95 | 0 | (gsdf_cnt > 0) && (jnd_min < jnd_max)) |
96 | 0 | { |
97 | 0 | int status = 0; |
98 | 0 | const unsigned long gin_ctn = (inverse) ? ddl_cnt : Count; // number of points to be interpolated |
99 | 0 | double *jidx = new double[gin_ctn]; |
100 | 0 | if (jidx != NULL) |
101 | 0 | { |
102 | 0 | const double dist = (jnd_max - jnd_min) / (gin_ctn - 1); // distance between two entries |
103 | 0 | unsigned long i; |
104 | 0 | double *s = jidx; |
105 | 0 | double value = jnd_min; // first value is fixed ! |
106 | 0 | for (i = gin_ctn; i > 1; --i) // initialize scaled JND index array |
107 | 0 | { |
108 | 0 | *(s++) = value; |
109 | 0 | value += dist; // add step by step ... |
110 | 0 | } |
111 | 0 | *s = jnd_max; // last value is fixed ! |
112 | 0 | double *jnd_idx = new double[gsdf_cnt]; |
113 | 0 | if (jnd_idx != NULL) |
114 | 0 | { |
115 | 0 | s = jnd_idx; |
116 | 0 | for (i = 0; i < gsdf_cnt; ++i) // initialize JND index array |
117 | 0 | *(s++) = i + 1; |
118 | 0 | double *gsdf = new double[gin_ctn]; // interpolated GSDF |
119 | 0 | if (gsdf != NULL) |
120 | 0 | { |
121 | 0 | if (DiCubicSpline<double, double>::Interpolation(jnd_idx, gsdf_tab, gsdf_spl, gsdf_cnt, jidx, gsdf, |
122 | 0 | OFstatic_cast(unsigned int, gin_ctn))) |
123 | 0 | { |
124 | 0 | DataBuffer = new Uint16[Count]; |
125 | 0 | if (DataBuffer != NULL) |
126 | 0 | { |
127 | 0 | const double amb = getAmbientLightValue(); |
128 | 0 | Uint16 *q = DataBuffer; |
129 | 0 | unsigned long j = 0; |
130 | | /* check whether to apply the inverse transformation */ |
131 | 0 | if (inverse) |
132 | 0 | { |
133 | 0 | double v; |
134 | 0 | const double factor = OFstatic_cast(double, ddl_cnt - 1) / OFstatic_cast(double, Count - 1); |
135 | | /* convert DDL to P-Value */ |
136 | 0 | for (i = 0; i < Count; ++i) |
137 | 0 | { |
138 | 0 | v = val_tab[OFstatic_cast(int, i * factor)] + amb; // need to scale index to range of value table |
139 | 0 | while ((j + 1 < ddl_cnt) && (gsdf[j] < v)) // search for closest index, assuming monotony |
140 | 0 | ++j; |
141 | 0 | if ((j > 0) && (fabs(gsdf[j - 1] - v) < fabs(gsdf[j] - v))) |
142 | 0 | --j; |
143 | 0 | *(q++) = ddl_tab[j]; |
144 | 0 | } |
145 | 0 | } else { |
146 | | /* initial DDL boundaries */ |
147 | 0 | unsigned long ddl_min = 0; |
148 | 0 | unsigned long ddl_max= ddl_cnt - 1; |
149 | | /* check whether minimum luminance is specified */ |
150 | 0 | if (lum_min >= 0) |
151 | 0 | { |
152 | 0 | j = ddl_min; |
153 | | /* determine corresponding minimum DDL value */ |
154 | 0 | while ((j < ddl_max) && (val_tab[j] + amb < lum_min)) |
155 | 0 | ++j; |
156 | 0 | ddl_min = j; |
157 | 0 | } |
158 | | /* check whether maximum luminance is specified */ |
159 | 0 | if (lum_max >= 0) |
160 | 0 | { |
161 | 0 | j = ddl_max; |
162 | | /* determine corresponding maximum DDL value */ |
163 | 0 | while ((j > ddl_min) && (val_tab[j] + amb > lum_max)) |
164 | 0 | --j; |
165 | 0 | ddl_max = j; |
166 | 0 | } |
167 | 0 | j = ddl_min; |
168 | 0 | const double *r = gsdf; |
169 | | /* convert P-Value to DDL */ |
170 | 0 | for (i = Count; i != 0; --i, ++r) |
171 | 0 | { |
172 | 0 | while ((j < ddl_max) && (val_tab[j] + amb < *r)) // search for closest index, assuming monotony |
173 | 0 | ++j; |
174 | 0 | if ((j > 0) && (fabs(val_tab[j - 1] + amb - *r) < fabs(val_tab[j] + amb - *r))) |
175 | 0 | --j; |
176 | 0 | *(q++) = ddl_tab[j]; |
177 | 0 | } |
178 | 0 | } |
179 | 0 | Data = DataBuffer; |
180 | 0 | if (stream != NULL) // write curve data to file |
181 | 0 | { |
182 | 0 | if (Count == ddl_cnt) // check whether GSDF LUT fits exactly to DISPLAY file |
183 | 0 | { |
184 | | /* use the standard "C" locale for proper decimal point */ |
185 | 0 | const STD_NAMESPACE locale oldLocale = stream->imbue(STD_NAMESPACE locale("C")); |
186 | 0 | for (i = 0; i < ddl_cnt; ++i) |
187 | 0 | { |
188 | 0 | (*stream) << ddl_tab[i]; // DDL |
189 | 0 | stream->setf(STD_NAMESPACE ios::fixed, STD_NAMESPACE ios::floatfield); |
190 | 0 | if (printMode) |
191 | 0 | (*stream) << "\t" << val_tab[i] + amb; // CC |
192 | 0 | (*stream) << "\t" << gsdf[i]; // GSDF |
193 | 0 | if (printMode) |
194 | 0 | { |
195 | 0 | if (inverse) |
196 | 0 | (*stream) << "\t" << gsdf[Data[i]]; // PSC' |
197 | 0 | else |
198 | 0 | (*stream) << "\t" << val_tab[Data[i]] + amb; // PSC |
199 | 0 | } |
200 | 0 | (*stream) << OFendl; |
201 | 0 | } |
202 | | /* reset locale */ |
203 | 0 | stream->imbue(oldLocale); |
204 | 0 | } else { |
205 | 0 | DCMIMGLE_WARN("can't write curve data, wrong DISPLAY file or GSDF LUT"); |
206 | 0 | } |
207 | 0 | } |
208 | 0 | status = 1; |
209 | 0 | } |
210 | 0 | } |
211 | 0 | } |
212 | 0 | delete[] gsdf; |
213 | 0 | } |
214 | 0 | delete[] jnd_idx; |
215 | 0 | } |
216 | 0 | delete[] jidx; |
217 | 0 | return status; |
218 | 0 | } |
219 | 0 | return 0; |
220 | 0 | } |