www.mooseframework.org
Euler2RGB.C
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://www.mooseframework.org
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 
45 Point
46 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  unsigned int ref_dir[3] = {0, 0, 0};
70  Real g[3][3] = {{0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
71  Real hkl[3] = {0.0, 0.0, 0.0};
72  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  switch (sd)
79  {
80  case 1: // 100
81  ref_dir[0] = 1;
82  ref_dir[1] = 0;
83  ref_dir[2] = 0;
84  break;
85 
86  case 2: // 010
87  ref_dir[0] = 0;
88  ref_dir[1] = 1;
89  ref_dir[2] = 0;
90  break;
91 
92  case 3: // 001
93  ref_dir[0] = 0;
94  ref_dir[1] = 0;
95  ref_dir[2] = 1;
96  break;
97  };
98 
99  // Define symmetry operators for each of the seven crystal systems
100  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  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  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  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  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  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  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  if (nsym == 0 || phase == 0)
240  RGB = 0;
241 
242  // Assign black RGB value for Euler angles outside of allowable range
243  else if (phi1 > pi_x2 || PHI > pi || phi2 > pi_x2)
244  RGB = 0;
245 
246  // Routine for valid set of Euler angles
247  else
248  {
249  // Construct 3X3 orientation matrix
250  g[0][0] = std::cos(phi1) * std::cos(phi2) - std::sin(phi1) * std::cos(PHI) * std::sin(phi2);
251  g[0][1] = std::sin(phi1) * std::cos(phi2) + std::cos(phi1) * std::cos(PHI) * std::sin(phi2);
252  g[0][2] = std::sin(phi2) * std::sin(PHI);
253  g[1][0] = -std::cos(phi1) * std::sin(phi2) - std::sin(phi1) * std::cos(PHI) * std::cos(phi2);
254  g[1][1] = -std::sin(phi1) * std::sin(phi2) + std::cos(phi1) * std::cos(PHI) * std::cos(phi2);
255  g[1][2] = std::cos(phi2) * std::sin(PHI);
256  g[2][0] = std::sin(phi1) * std::sin(PHI);
257  g[2][1] = -std::cos(phi1) * std::sin(PHI);
258  g[2][2] = std::cos(PHI);
259 
260  // Loop to sort Euler angles into standard stereographic triangle (SST)
261  index = 0;
262  while (index < nsym)
263  {
264  // Form orientation matrix
265  for (unsigned int i = 0; i < 3; ++i)
266  for (unsigned int j = 0; j < 3; ++j)
267  {
268  S[i][j] = 0.0;
269  for (unsigned int k = 0; k < 3; ++k)
270  S[i][j] += SymOps[index][i][k] * g[k][j];
271  }
272 
273  // Multiple orientation matrix by reference sample direction
274  for (unsigned int i = 0; i < 3; ++i)
275  {
276  hkl[i] = 0;
277  for (unsigned int j = 0; j < 3; ++j)
278  hkl[i] += S[i][j] * ref_dir[j];
279  }
280 
281  // Convert to spherical coordinates (ignore "r" variable since r=1)
282  eta = std::abs(std::atan2(hkl[1], hkl[0]));
283  chi = std::acos(std::abs(hkl[2]));
284 
285  // Continue if eta and chi values are within the SST
286  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  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  if (sym == 43)
299  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  red = std::abs(1.0 - (chi / chi_max2));
305  blue = std::abs((eta - eta_min) / (eta_max - eta_min));
306  green = 1.0 - blue;
307 
308  blue = blue * (chi / chi_max2);
309  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  RGB(0) = std::sqrt(red);
315  RGB(1) = std::sqrt(green);
316  RGB(2) = std::sqrt(blue);
317 
318  // Find maximum value of red, green, or blue
319  maxRGB = std::max({RGB(0), RGB(1), RGB(2)});
320 
321  // Normalize position of SST center point
322  RGB /= maxRGB;
323  }
324 
325  return RGB;
326 }
ADRealEigenVector< T, D, asd > sqrt(const ADRealEigenVector< T, D, asd > &)
ADRealEigenVector< T, D, asd > abs(const ADRealEigenVector< T, D, asd > &)
static const std::string S
Definition: NS.h:148
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
static const std::complex< double > j(0, 1)
Complex number "j" (also known as "i")
static const std::string k
Definition: NS.h:124
Point euler2RGB(unsigned int sd, Real phi1, Real PHI, Real phi2, unsigned int phase, unsigned int sym)
This function rotates a set of three Bunge Euler angles into the standard Stereographic triangle...
Definition: Euler2RGB.C:46
const Real pi