Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #include "Euler2RGB.h"
11 : #include "MathUtils.h"
12 : #include "libmesh/utility.h"
13 :
14 : /**
15 : * This function rotates a set of three Bunge Euler angles into
16 : * the standard Stereographic triangle, interpolates the RGB color
17 : * value based on a user selected reference sample direction, and
18 : * outputs an integer value representing the RGB tuplet usable for
19 : * plotting inverse pole figure colored grain maps. The program
20 : * can accommodate any of the seven crystal systems.
21 : *
22 : * Inputs: 1) Reference sample direction (sd) input as an integer between
23 : * 1 and 3. Options are [100], [010], and [001] with [001]
24 : * being the most common choice found in the literature. Input
25 : * "1" to select [100], "2" to select [010], and "3" to select
26 : * the [001] sample direction.
27 : * 2) Set of three Euler angles (phi1, PHI, phi2) in the Bunge
28 : * notation (must be in radians)
29 : * 3) Phase number "phase" used to assign black color values to
30 : * voids (phase = 0)
31 : * 4) Integer value of crystal class as defined by OIM Software:
32 : * "43" for cubic
33 : * "62" for hexagonal
34 : * "42" for tetragonal
35 : * "32" for trigonal
36 : * "22" for orthorhombic
37 : * "2" for monoclinic
38 : * "1" for triclinic
39 : * "0" for unindexed points (bad data point)
40 : *
41 : * Output: Integer value of RGB calculated from:
42 : * RGB = red*256^2 + green*256 + blue (where red, green, and blue
43 : * are integer values between 0 and 255)
44 : */
45 : Point
46 150034 : euler2RGB(unsigned int sd, Real phi1, Real PHI, Real phi2, unsigned int phase, unsigned int sym)
47 : {
48 : // Define Constants
49 : const Real pi = libMesh::pi;
50 : const Real pi_x2 = 2.0 * pi;
51 : const Real a = std::sqrt(3.0) / 2.0;
52 :
53 : // Preallocate and zero variables
54 : unsigned int index = 0;
55 : unsigned int nsym = 1;
56 :
57 : Real blue = 0.0;
58 : Real chi = 0.0;
59 : Real chi_min = 0.0;
60 : Real chi_max = 0.0;
61 : Real chi_max2 = 0.0;
62 : Real eta = 0.0;
63 : Real eta_min = 0.0;
64 : Real eta_max = 0.0;
65 : Real green = 0.0;
66 : Real maxRGB = 0.0;
67 : Real red = 0.0;
68 :
69 150034 : unsigned int ref_dir[3] = {0, 0, 0};
70 150034 : Real g[3][3] = {{0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
71 150034 : Real hkl[3] = {0.0, 0.0, 0.0};
72 150034 : Real S[3][3] = {{0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
73 : const Real(*SymOps)[3][3];
74 :
75 : Point RGB;
76 :
77 : // Assign reference sample direction
78 150034 : switch (sd)
79 : {
80 4 : case 1: // 100
81 4 : ref_dir[0] = 1;
82 : ref_dir[1] = 0;
83 : ref_dir[2] = 0;
84 4 : break;
85 :
86 0 : case 2: // 010
87 : ref_dir[0] = 0;
88 0 : ref_dir[1] = 1;
89 : ref_dir[2] = 0;
90 0 : break;
91 :
92 150030 : case 3: // 001
93 : ref_dir[0] = 0;
94 : ref_dir[1] = 0;
95 150030 : ref_dir[2] = 1;
96 150030 : break;
97 : };
98 :
99 : // Define symmetry operators for each of the seven crystal systems
100 150034 : const Real SymOpsCubic[24][3][3] = {
101 : {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}, {{0, 0, 1}, {1, 0, 0}, {0, 1, 0}},
102 : {{0, 1, 0}, {0, 0, 1}, {1, 0, 0}}, {{0, -1, 0}, {0, 0, 1}, {-1, 0, 0}},
103 : {{0, -1, 0}, {0, 0, -1}, {1, 0, 0}}, {{0, 1, 0}, {0, 0, -1}, {-1, 0, 0}},
104 : {{0, 0, -1}, {1, 0, 0}, {0, -1, 0}}, {{0, 0, -1}, {-1, 0, 0}, {0, 1, 0}},
105 : {{0, 0, 1}, {-1, 0, 0}, {0, -1, 0}}, {{-1, 0, 0}, {0, 1, 0}, {0, 0, -1}},
106 : {{-1, 0, 0}, {0, -1, 0}, {0, 0, 1}}, {{1, 0, 0}, {0, -1, 0}, {0, 0, -1}},
107 : {{0, 0, -1}, {0, -1, 0}, {-1, 0, 0}}, {{0, 0, 1}, {0, -1, 0}, {1, 0, 0}},
108 : {{0, 0, 1}, {0, 1, 0}, {-1, 0, 0}}, {{0, 0, -1}, {0, 1, 0}, {1, 0, 0}},
109 : {{-1, 0, 0}, {0, 0, -1}, {0, -1, 0}}, {{1, 0, 0}, {0, 0, -1}, {0, 1, 0}},
110 : {{1, 0, 0}, {0, 0, 1}, {0, -1, 0}}, {{-1, 0, 0}, {0, 0, 1}, {0, 1, 0}},
111 : {{0, -1, 0}, {-1, 0, 0}, {0, 0, -1}}, {{0, 1, 0}, {-1, 0, 0}, {0, 0, -1}},
112 : {{0, 1, 0}, {1, 0, 0}, {0, 0, -1}}, {{0, -1, 0}, {1, 0, 0}, {0, 0, 1}}};
113 :
114 150034 : const Real SymOpsHexagonal[12][3][3] = {{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}},
115 : {{-0.5, a, 0}, {-a, -0.5, 0}, {0, 0, 1}},
116 : {{-0.5, -a, 0}, {a, -0.5, 0}, {0, 0, 1}},
117 : {{0.5, a, 0}, {-a, 0.5, 0}, {0, 0, 1}},
118 : {{-1, 0, 0}, {0, -1, 0}, {0, 0, 1}},
119 : {{0.5, -a, 0}, {a, 0.5, 0}, {0, 0, 1}},
120 : {{-0.5, -a, 0}, {-a, 0.5, 0}, {0, 0, -1}},
121 : {{1, 0, 0}, {0, -1, 0}, {0, 0, -1}},
122 : {{-0.5, a, 0}, {a, 0.5, 0}, {0, 0, -1}},
123 : {{0.5, a, 0}, {a, -0.5, 0}, {0, 0, -1}},
124 : {{-1, 0, 0}, {0, 1, 0}, {0, 0, -1}},
125 : {{0.5, -a, 0}, {-a, -0.5, 0}, {0, 0, -1}}};
126 :
127 150034 : const Real SymOpsTetragonal[8][3][3] = {{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}},
128 : {{-1, 0, 0}, {0, 1, 0}, {0, 0, -1}},
129 : {{1, 0, 0}, {0, -1, 0}, {0, 0, -1}},
130 : {{-1, 0, 0}, {0, -1, 0}, {0, 0, 1}},
131 : {{0, 1, 0}, {-1, 0, 0}, {0, 0, 1}},
132 : {{0, -1, 0}, {1, 0, 0}, {0, 0, 1}},
133 : {{0, 1, 0}, {1, 0, 0}, {0, 0, -1}},
134 : {{0, -1, 0}, {-1, 0, 0}, {0, 0, -1}}};
135 :
136 150034 : const Real SymOpsTrigonal[6][3][3] = {{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}},
137 : {{-0.5, a, 0}, {-a, -0.5, 0}, {0, 0, 1}},
138 : {{-0.5, -a, 0}, {a, -0.5, 0}, {0, 0, 1}},
139 : {{0.5, a, 0}, {a, -0.5, 0}, {0, 0, -1}},
140 : {{-1, 0, 0}, {0, 1, 0}, {0, 0, 1}},
141 : {{0.5, -a, 0}, {-a, -0.5, 0}, {0, 0, -1}}};
142 :
143 150034 : const Real SymOpsOrthorhombic[4][3][3] = {{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}},
144 : {{-1, 0, 0}, {0, 1, 0}, {0, 0, -1}},
145 : {{1, 0, 0}, {0, -1, 0}, {0, 0, 1}},
146 : {{-1, 0, 0}, {0, -1, 0}, {0, 0, 1}}};
147 :
148 150034 : const Real SymOpsMonoclinic[2][3][3] = {{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}},
149 : {{-1, 0, 0}, {0, 1, 0}, {0, 0, -1}}};
150 :
151 150034 : const Real SymOpsTriclinic[1][3][3] = {{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}};
152 :
153 : // Assign parameters based on crystal class (sym)
154 : // Load cubic parameters (class 432)
155 : if (sym == 43)
156 : {
157 : nsym = 24;
158 : SymOps = SymOpsCubic;
159 : eta_min = 0 * (pi / 180);
160 : eta_max = 45 * (pi / 180);
161 : chi_min = 0 * (pi / 180);
162 : chi_max = std::acos(std::sqrt(1.0 / (2.0 + (Utility::pow<2>(std::tan(eta_max))))));
163 : }
164 :
165 : // Load hexagonal parameters (class 622)
166 : else if (sym == 62)
167 : {
168 : nsym = 12;
169 : SymOps = SymOpsHexagonal;
170 : eta_min = 0 * (pi / 180);
171 : eta_max = 30 * (pi / 180);
172 : chi_min = 0 * (pi / 180);
173 : chi_max = pi / 2;
174 : }
175 :
176 : // Load tetragonal parameters (class 422)
177 : else if (sym == 42)
178 : {
179 : nsym = 8;
180 : SymOps = SymOpsTetragonal;
181 : eta_min = 0 * (pi / 180);
182 : eta_max = 45 * (pi / 180);
183 : chi_min = 0 * (pi / 180);
184 : chi_max = pi / 2;
185 : }
186 :
187 : // Load trigonal parameters (class 32)
188 : else if (sym == 32)
189 : {
190 : nsym = 6;
191 : SymOps = SymOpsTrigonal;
192 : eta_min = 0 * (pi / 180);
193 : eta_max = 60 * (pi / 180);
194 : chi_min = 0 * (pi / 180);
195 : chi_max = pi / 2;
196 : }
197 :
198 : // Load orthorhombic parameters (class 22)
199 : else if (sym == 22)
200 : {
201 : nsym = 4;
202 : SymOps = SymOpsOrthorhombic;
203 : eta_min = 0 * (pi / 180);
204 : eta_max = 90 * (pi / 180);
205 : chi_min = 0 * (pi / 180);
206 : chi_max = pi / 2;
207 : }
208 :
209 : // Load monoclinic parameters (class 2)
210 : else if (sym == 2)
211 : {
212 : nsym = 2;
213 : SymOps = SymOpsMonoclinic;
214 : eta_min = 0 * (pi / 180);
215 : eta_max = 180 * (pi / 180);
216 : chi_min = 0 * (pi / 180);
217 : chi_max = pi / 2;
218 : }
219 :
220 : // Load triclinic parameters (class 1)
221 : else if (sym == 1)
222 : {
223 : nsym = 1;
224 : SymOps = SymOpsTriclinic;
225 : eta_min = 0 * (pi / 180);
226 : eta_max = 360 * (pi / 180);
227 : chi_min = 0 * (pi / 180);
228 : chi_max = pi / 2;
229 : }
230 :
231 : // Accomodate non-conforming (bad) data points
232 : else
233 : {
234 : nsym = 0;
235 : }
236 :
237 : // Start of main routine //
238 : // Assign black RGB values for bad data points (nsym = 0) or voids (phase = 0)
239 150034 : if (nsym == 0 || phase == 0)
240 258 : RGB = 0;
241 :
242 : // Assign black RGB value for Euler angles outside of allowable range
243 149776 : else if (phi1 > pi_x2 || PHI > pi || phi2 > pi_x2)
244 0 : RGB = 0;
245 :
246 : // Routine for valid set of Euler angles
247 : else
248 : {
249 : // Construct 3X3 orientation matrix
250 149776 : g[0][0] = std::cos(phi1) * std::cos(phi2) - std::sin(phi1) * std::cos(PHI) * std::sin(phi2);
251 149776 : g[0][1] = std::sin(phi1) * std::cos(phi2) + std::cos(phi1) * std::cos(PHI) * std::sin(phi2);
252 149776 : g[0][2] = std::sin(phi2) * std::sin(PHI);
253 149776 : g[1][0] = -std::cos(phi1) * std::sin(phi2) - std::sin(phi1) * std::cos(PHI) * std::cos(phi2);
254 149776 : g[1][1] = -std::sin(phi1) * std::sin(phi2) + std::cos(phi1) * std::cos(PHI) * std::cos(phi2);
255 149776 : g[1][2] = std::cos(phi2) * std::sin(PHI);
256 149776 : g[2][0] = std::sin(phi1) * std::sin(PHI);
257 149776 : g[2][1] = -std::cos(phi1) * std::sin(PHI);
258 149776 : g[2][2] = std::cos(PHI);
259 :
260 : // Loop to sort Euler angles into standard stereographic triangle (SST)
261 : index = 0;
262 1480525 : while (index < nsym)
263 : {
264 : // Form orientation matrix
265 5922100 : for (unsigned int i = 0; i < 3; ++i)
266 17766300 : for (unsigned int j = 0; j < 3; ++j)
267 : {
268 13324725 : S[i][j] = 0.0;
269 53298900 : for (unsigned int k = 0; k < 3; ++k)
270 39974175 : S[i][j] += SymOps[index][i][k] * g[k][j];
271 : }
272 :
273 : // Multiple orientation matrix by reference sample direction
274 5922100 : for (unsigned int i = 0; i < 3; ++i)
275 : {
276 4441575 : hkl[i] = 0;
277 17766300 : for (unsigned int j = 0; j < 3; ++j)
278 13324725 : hkl[i] += S[i][j] * ref_dir[j];
279 : }
280 :
281 : // Convert to spherical coordinates (ignore "r" variable since r=1)
282 1480525 : eta = std::abs(std::atan2(hkl[1], hkl[0]));
283 1480525 : chi = std::acos(std::abs(hkl[2]));
284 :
285 : // Continue if eta and chi values are within the SST
286 1480525 : if (eta >= eta_min && eta < eta_max && chi >= chi_min && chi < chi_max)
287 : break;
288 :
289 : // Increment to next symmetry operator if not in SST
290 : else
291 1330749 : index++;
292 :
293 : // Check for solution
294 : mooseAssert(index != nsym, "Euler2RGB failed to map the supplied Euler angle into the SST!");
295 : }
296 :
297 : // Adjust maximum chi value to ensure it falls within the SST (cubic materials only)
298 149776 : if (sym == 43)
299 149775 : chi_max2 = std::acos(std::sqrt(1.0 / (2.0 + (Utility::pow<2>(std::tan(eta))))));
300 : else
301 : chi_max2 = pi / 2;
302 :
303 : // Calculate the RGB color values and make adjustments to maximize colorspace
304 149776 : red = std::abs(1.0 - (chi / chi_max2));
305 149776 : blue = std::abs((eta - eta_min) / (eta_max - eta_min));
306 149776 : green = 1.0 - blue;
307 :
308 149776 : blue = blue * (chi / chi_max2);
309 149776 : green = green * (chi / chi_max2);
310 :
311 : // Check for negative RGB values before taking square root
312 : mooseAssert(red >= 0 || green >= 0 || blue >= 0, "RGB component values must be positive!");
313 :
314 149776 : RGB(0) = std::sqrt(red);
315 149776 : RGB(1) = std::sqrt(green);
316 149776 : RGB(2) = std::sqrt(blue);
317 :
318 : // Find maximum value of red, green, or blue
319 149776 : maxRGB = std::max({RGB(0), RGB(1), RGB(2)});
320 :
321 : // Normalize position of SST center point
322 : RGB /= maxRGB;
323 : }
324 :
325 150034 : return RGB;
326 : }
|
