www.mooseframework.org
Euler2RGB.C
Go to the documentation of this file.
1 /****************************************************************/
2 /* MOOSE - Multiphysics Object Oriented Simulation Environment */
3 /* */
4 /* All contents are licensed under LGPL V2.1 */
5 /* See LICENSE for full restrictions */
6 /****************************************************************/
7 #include "Euler2RGB.h"
8 #include "MathUtils.h"
9 #include "libmesh/utility.h"
10 
42 Point
43 euler2RGB(unsigned int sd, Real phi1, Real PHI, Real phi2, unsigned int phase, unsigned int sym)
44 {
45  // Define Constants
46  const Real pi = libMesh::pi;
47  const Real pi_x2 = 2.0 * pi;
48  const Real a = std::sqrt(3.0) / 2.0;
49 
50  // Preallocate and zero variables
51  unsigned int index = 0;
52  unsigned int nsym = 1;
53 
54  Real blue = 0.0;
55  Real chi = 0.0;
56  Real chi_min = 0.0;
57  Real chi_max = 0.0;
58  Real chi_max2 = 0.0;
59  Real eta = 0.0;
60  Real eta_min = 0.0;
61  Real eta_max = 0.0;
62  Real green = 0.0;
63  Real maxRGB = 0.0;
64  Real red = 0.0;
65 
66  unsigned int ref_dir[3] = {0, 0, 0};
67  Real g[3][3] = {{0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
68  Real hkl[3] = {0.0, 0.0, 0.0};
69  Real S[3][3] = {{0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
70  const Real(*SymOps)[3][3];
71 
72  Point RGB;
73 
74  // Assign reference sample direction
75  switch (sd)
76  {
77  case 1: // 100
78  ref_dir[0] = 1;
79  ref_dir[1] = 0;
80  ref_dir[2] = 0;
81  break;
82 
83  case 2: // 010
84  ref_dir[0] = 0;
85  ref_dir[1] = 1;
86  ref_dir[2] = 0;
87  break;
88 
89  case 3: // 001
90  ref_dir[0] = 0;
91  ref_dir[1] = 0;
92  ref_dir[2] = 1;
93  break;
94  };
95 
96  // Define symmetry operators for each of the seven crystal systems
97  const Real SymOpsCubic[24][3][3] = {
98  {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}, {{0, 0, 1}, {1, 0, 0}, {0, 1, 0}},
99  {{0, 1, 0}, {0, 0, 1}, {1, 0, 0}}, {{0, -1, 0}, {0, 0, 1}, {-1, 0, 0}},
100  {{0, -1, 0}, {0, 0, -1}, {1, 0, 0}}, {{0, 1, 0}, {0, 0, -1}, {-1, 0, 0}},
101  {{0, 0, -1}, {1, 0, 0}, {0, -1, 0}}, {{0, 0, -1}, {-1, 0, 0}, {0, 1, 0}},
102  {{0, 0, 1}, {-1, 0, 0}, {0, -1, 0}}, {{-1, 0, 0}, {0, 1, 0}, {0, 0, -1}},
103  {{-1, 0, 0}, {0, -1, 0}, {0, 0, 1}}, {{1, 0, 0}, {0, -1, 0}, {0, 0, -1}},
104  {{0, 0, -1}, {0, -1, 0}, {-1, 0, 0}}, {{0, 0, 1}, {0, -1, 0}, {1, 0, 0}},
105  {{0, 0, 1}, {0, 1, 0}, {-1, 0, 0}}, {{0, 0, -1}, {0, 1, 0}, {1, 0, 0}},
106  {{-1, 0, 0}, {0, 0, -1}, {0, -1, 0}}, {{1, 0, 0}, {0, 0, -1}, {0, 1, 0}},
107  {{1, 0, 0}, {0, 0, 1}, {0, -1, 0}}, {{-1, 0, 0}, {0, 0, 1}, {0, 1, 0}},
108  {{0, -1, 0}, {-1, 0, 0}, {0, 0, -1}}, {{0, 1, 0}, {-1, 0, 0}, {0, 0, -1}},
109  {{0, 1, 0}, {1, 0, 0}, {0, 0, -1}}, {{0, -1, 0}, {1, 0, 0}, {0, 0, 1}}};
110 
111  const Real SymOpsHexagonal[12][3][3] = {{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}},
112  {{-0.5, a, 0}, {-a, -0.5, 0}, {0, 0, 1}},
113  {{-0.5, -a, 0}, {a, -0.5, 0}, {0, 0, 1}},
114  {{0.5, a, 0}, {-a, 0.5, 0}, {0, 0, 1}},
115  {{-1, 0, 0}, {0, -1, 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 
124  const Real SymOpsTetragonal[8][3][3] = {{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}},
125  {{-1, 0, 0}, {0, 1, 0}, {0, 0, -1}},
126  {{1, 0, 0}, {0, -1, 0}, {0, 0, -1}},
127  {{-1, 0, 0}, {0, -1, 0}, {0, 0, 1}},
128  {{0, 1, 0}, {-1, 0, 0}, {0, 0, 1}},
129  {{0, -1, 0}, {1, 0, 0}, {0, 0, 1}},
130  {{0, 1, 0}, {1, 0, 0}, {0, 0, -1}},
131  {{0, -1, 0}, {-1, 0, 0}, {0, 0, -1}}};
132 
133  const Real SymOpsTrigonal[6][3][3] = {{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}},
134  {{-0.5, a, 0}, {-a, -0.5, 0}, {0, 0, 1}},
135  {{-0.5, -a, 0}, {a, -0.5, 0}, {0, 0, 1}},
136  {{0.5, a, 0}, {a, -0.5, 0}, {0, 0, -1}},
137  {{-1, 0, 0}, {0, 1, 0}, {0, 0, 1}},
138  {{0.5, -a, 0}, {-a, -0.5, 0}, {0, 0, -1}}};
139 
140  const Real SymOpsOrthorhombic[4][3][3] = {{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}},
141  {{-1, 0, 0}, {0, 1, 0}, {0, 0, -1}},
142  {{1, 0, 0}, {0, -1, 0}, {0, 0, 1}},
143  {{-1, 0, 0}, {0, -1, 0}, {0, 0, 1}}};
144 
145  const Real SymOpsMonoclinic[2][3][3] = {{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}},
146  {{-1, 0, 0}, {0, 1, 0}, {0, 0, -1}}};
147 
148  const Real SymOpsTriclinic[1][3][3] = {{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}};
149 
150  // Assign parameters based on crystal class (sym)
151  // Load cubic parameters (class 432)
152  if (sym == 43)
153  {
154  nsym = 24;
155  SymOps = SymOpsCubic;
156  eta_min = 0 * (pi / 180);
157  eta_max = 45 * (pi / 180);
158  chi_min = 0 * (pi / 180);
159  chi_max = std::acos(std::sqrt(1.0 / (2.0 + (std::tan(Utility::pow<2>(eta_max))))));
160  }
161 
162  // Load hexagonal parameters (class 622)
163  else if (sym == 62)
164  {
165  nsym = 12;
166  SymOps = SymOpsHexagonal;
167  eta_min = 0 * (pi / 180);
168  eta_max = 30 * (pi / 180);
169  chi_min = 0 * (pi / 180);
170  chi_max = pi / 2;
171  }
172 
173  // Load tetragonal parameters (class 422)
174  else if (sym == 42)
175  {
176  nsym = 8;
177  SymOps = SymOpsTetragonal;
178  eta_min = 0 * (pi / 180);
179  eta_max = 45 * (pi / 180);
180  chi_min = 0 * (pi / 180);
181  chi_max = pi / 2;
182  }
183 
184  // Load trigonal parameters (class 32)
185  else if (sym == 32)
186  {
187  nsym = 6;
188  SymOps = SymOpsTrigonal;
189  eta_min = 0 * (pi / 180);
190  eta_max = 60 * (pi / 180);
191  chi_min = 0 * (pi / 180);
192  chi_max = pi / 2;
193  }
194 
195  // Load orthorhombic parameters (class 22)
196  else if (sym == 22)
197  {
198  nsym = 4;
199  SymOps = SymOpsOrthorhombic;
200  eta_min = 0 * (pi / 180);
201  eta_max = 90 * (pi / 180);
202  chi_min = 0 * (pi / 180);
203  chi_max = pi / 2;
204  }
205 
206  // Load monoclinic parameters (class 2)
207  else if (sym == 2)
208  {
209  nsym = 2;
210  SymOps = SymOpsMonoclinic;
211  eta_min = 0 * (pi / 180);
212  eta_max = 180 * (pi / 180);
213  chi_min = 0 * (pi / 180);
214  chi_max = pi / 2;
215  }
216 
217  // Load triclinic parameters (class 1)
218  else if (sym == 1)
219  {
220  nsym = 1;
221  SymOps = SymOpsTriclinic;
222  eta_min = 0 * (pi / 180);
223  eta_max = 360 * (pi / 180);
224  chi_min = 0 * (pi / 180);
225  chi_max = pi / 2;
226  }
227 
228  // Accomodate non-conforming (bad) data points
229  else
230  {
231  nsym = 0;
232  }
233 
234  // Start of main routine //
235  // Assign black RGB values for bad data points (nsym = 0) or voids (phase = 0)
236  if (nsym == 0 || phase == 0)
237  RGB = 0;
238 
239  // Assign black RGB value for Euler angles outside of allowable range
240  else if (phi1 > pi_x2 || PHI > pi || phi2 > pi_x2)
241  RGB = 0;
242 
243  // Routine for valid set of Euler angles
244  else
245  {
246  // Construct 3X3 orientation matrix
247  g[0][0] = std::cos(phi1) * std::cos(phi2) - std::sin(phi1) * std::cos(PHI) * std::sin(phi2);
248  g[0][1] = std::sin(phi1) * std::cos(phi2) + std::cos(phi1) * std::cos(PHI) * std::sin(phi2);
249  g[0][2] = std::sin(phi2) * std::sin(PHI);
250  g[1][0] = -std::cos(phi1) * std::sin(phi2) - std::sin(phi1) * std::cos(PHI) * std::cos(phi2);
251  g[1][1] = -std::sin(phi1) * std::sin(phi2) + std::cos(phi1) * std::cos(PHI) * std::cos(phi2);
252  g[1][2] = std::cos(phi2) * std::sin(PHI);
253  g[2][0] = std::sin(phi1) * std::sin(PHI);
254  g[2][1] = -std::cos(phi1) * std::sin(PHI);
255  g[2][2] = std::cos(PHI);
256 
257  // Loop to sort Euler angles into standard stereographic triangle (SST)
258  index = 0;
259  while (index < nsym)
260  {
261  // Form orientation matrix
262  for (unsigned int i = 0; i < 3; ++i)
263  for (unsigned int j = 0; j < 3; ++j)
264  {
265  S[i][j] = 0.0;
266  for (unsigned int k = 0; k < 3; ++k)
267  S[i][j] += SymOps[index][i][k] * g[k][j];
268  }
269 
270  // Multiple orientation matrix by reference sample direction
271  for (unsigned int i = 0; i < 3; ++i)
272  {
273  hkl[i] = 0;
274  for (unsigned int j = 0; j < 3; ++j)
275  hkl[i] += S[i][j] * ref_dir[j];
276  }
277 
278  // Convert to spherical coordinates (ignore "r" variable since r=1)
279  eta = std::abs(std::atan2(hkl[1], hkl[0]));
280  chi = std::acos(std::abs(hkl[2]));
281 
282  // Continue if eta and chi values are within the SST
283  if (eta >= eta_min && eta < eta_max && chi >= chi_min && chi < chi_max)
284  break;
285 
286  // Increment to next symmetry operator if not in SST
287  else
288  index++;
289  }
290 
291  // Adjust maximum chi value to ensure it falls within the SST (cubic materials only)
292  if (sym == 43)
293  chi_max2 = std::acos(std::sqrt(1.0 / (2.0 + (std::tan(Utility::pow<2>(eta))))));
294  else
295  chi_max2 = pi / 2;
296 
297  // Calculate the RGB color values and make adjustments to maximize colorspace
298  red = std::abs(1.0 - (chi / chi_max2));
299  blue = std::abs((eta - eta_min) / (eta_max - eta_min));
300  green = 1.0 - blue;
301 
302  blue = blue * (chi / chi_max2);
303  green = green * (chi / chi_max2);
304 
305  RGB(0) = std::sqrt(red);
306  RGB(1) = std::sqrt(green);
307  RGB(2) = std::sqrt(blue);
308 
309  // Find maximum value of red, green, or blue
310  maxRGB = std::max(RGB(0), std::max(RGB(1), RGB(2)));
311 
312  // Adjust RGB values to enforce white center point instead of black
313  RGB /= maxRGB;
314  }
315 
316  return RGB;
317 }
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:43