← Back to branch summary

~ubuntu-branches/ubuntu/precise/openwalnut/precise

« back to all changes in this revision

Viewing changes to src/core/graphicsEngine/shaders/shaders/WGETensorTools.glsl

  • Committer: Bazaar Package Importer
  • Author(s): Sebastian Eichelbaum
  • Date: 2011-06-21 10:26:54 UTC
  • Revision ID: james.westby@ubuntu.com-20110621102654-rq0zf436q949biih
Tags: upstream-1.2.5
Import upstream version 1.2.5

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/core/graphicsEngine/shaders/shaders/WGETensorTools.glsl
 
1
//---------------------------------------------------------------------------
 
2
//
 
3
// Project: OpenWalnut ( http://www.openwalnut.org )
 
4
//
 
5
// Copyright 2009 OpenWalnut Community, BSV-Leipzig and CNCF-CBS
 
6
// For more information see http://www.openwalnut.org/copying
 
7
//
 
8
// This file is part of OpenWalnut.
 
9
//
 
10
// OpenWalnut is free software: you can redistribute it and/or modify
 
11
// it under the terms of the GNU Lesser General Public License as published by
 
12
// the Free Software Foundation, either version 3 of the License, or
 
13
// (at your option) any later version.
 
14
//
 
15
// OpenWalnut is distributed in the hope that it will be useful,
 
16
// but WITHOUT ANY WARRANTY; without even the implied warranty of
 
17
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
18
// GNU Lesser General Public License for more details.
 
19
//
 
20
// You should have received a copy of the GNU Lesser General Public License
 
21
// along with OpenWalnut. If not, see <http://www.gnu.org/licenses/>.
 
22
//
 
23
//---------------------------------------------------------------------------
 
24
 
 
25
#ifndef WGETENSORTOOLS_GLSL
 
26
#define WGETENSORTOOLS_GLSL
 
27
 
 
28
#version 120
 
29
 
 
30
// (c) 2007 by Mario Hlawitschka
 
31
 
 
32
// tensors have to be stored in the following way:
 
33
// diag.x = Dxx     Diagonal elements are stored in one vector
 
34
// diag.y = Dyy
 
35
// diag.z = Dzz
 
36
// offdiag.x = Dyz  Off diagonal elements are stored in another vector
 
37
// offdiag.y = Dxz      where all components are stored in that location
 
38
// offdiag.z = Dxy      that is indexed by the coordinate that is not in the tensor index
 
39
 
 
40
 
 
41
// compute the eigenvalues of a tensor
 
42
// (Hasan et al. 2001)
 
43
vec3 getEigenvaluesHasan( vec3 diag, vec3 offdiag )
 
44
{
 
45
  const float PiThird = 3.14159265358979323846264 / 3.;
 
46
  float I1 = diag.x + diag.y + diag.z;
 
47
  float I2 = diag.x * diag.y + diag.x * diag.z + diag.y * diag.z - dot( offdiag, offdiag );
 
48
  float I3 = diag.x * diag.y * diag.z + 2. * offdiag.x * offdiag.y * offdiag.z
 
49
      - ( diag.x * offdiag.x * offdiag.x + diag.y * offdiag.y * offdiag.y + diag.z * offdiag.z * offdiag.z );
 
50
 
 
51
  const float third = 1. / 3.;
 
52
  float I1third = I1 * third;
 
53
  float I1thirdsqr = I1third * I1third;
 
54
  float I2third = I2 * third;
 
55
  float v  = I1thirdsqr - I2third;
 
56
 
 
57
  vec3 lambda;
 
58
 
 
59
 
 
60
  // debug
 
61
  //if( v < 0.000001 ) lambda = vec3(1.,0.,0.);
 
62
  //else
 
63
  {
 
64
      float s  = I1thirdsqr * I1third - I1 * I2/6. + 0.5 * I3;
 
65
 
 
66
      // for real eigenvalues: v>0, s^2 < v^3
 
67
      float sqrtv = sqrt( v );
 
68
      float phi= acos( s / ( v * sqrtv ) ) * third;
 
69
 
 
70
      float sqrtv2 = 2. * sqrtv;
 
71
 
 
72
      // due to the cosine function and the fact, that 0 <= phi <= pi/3
 
73
      // it is obvious that the eigenvalues need no further sorting
 
74
      lambda.x = I1third + sqrtv2 * cos( phi );
 
75
      lambda.y = I1third - sqrtv2 * cos( PiThird + phi );
 
76
      lambda.z = I1third - sqrtv2 * cos( PiThird - phi );
 
77
  }
 
78
  return lambda;
 
79
}
 
80
 
 
81
 
 
82
// Method similar to the above by Hasan proposed by Cardano et al.
 
83
// implementation is more stable than the above one, so use this as
 
84
// default in all applications
 
85
float sqr( float a )
 
86
{
 
87
    return a * a;
 
88
}
 
89
 
 
90
vec3 getEigenvaluesCardano( in vec3 diag, in vec3 offdiag )
 
91
{
 
92
  const float M_SQRT3 = 1.73205080756887729352744634151;
 
93
  float de = offdiag.z * offdiag.x;
 
94
  float dd = sqr( offdiag.z );
 
95
  float ee = sqr( offdiag.x );
 
96
  float ff = sqr( offdiag.y );
 
97
  float m  = diag.x + diag.y + diag.z;
 
98
  float c1 = diag.x * diag.y + diag.x * diag.z + diag.y * diag.z
 
99
             - ( dd + ee + ff );
 
100
  float c0 = diag.z * dd + diag.x * ee + diag.y * ff - diag.x * diag.y * diag.z - 2. * offdiag.y * de;
 
101
 
 
102
  float p, sqrt_p, q, c, s, phi;
 
103
  p = sqr( m ) - 3. * c1;
 
104
  q = m * ( p - ( 3. / 2. ) * c1 ) - ( 27. / 2. ) * c0;
 
105
  sqrt_p = sqrt( abs( p ) );
 
106
 
 
107
  phi = 27. * ( 0.25 * sqr( c1 ) * ( p - c1 ) + c0 * ( q + 27. / 4. * c0 ) );
 
108
  phi = ( 1. / 3. ) * atan( sqrt( abs( phi ) ), q );
 
109
 
 
110
  c = sqrt_p * cos( phi );
 
111
  s = ( 1. / M_SQRT3 ) * sqrt_p * sin( phi );
 
112
 
 
113
  vec3 w;
 
114
  w[2] = ( 1. / 3. ) * ( m - c );
 
115
  w[1] = w[2] + s;
 
116
  w[0] = w[2] + c;
 
117
  w[2] -= s;
 
118
  return w;
 
119
}
 
120
 
 
121
vec3 getEigenvalues( in vec3 diag, in vec3 offdiag )
 
122
{
 
123
  return getEigenvaluesCardano( diag, offdiag );
 
124
}
 
125
 
 
126
float getMajorEigenvalue( vec3 diag, vec3 offdiag )
 
127
{
 
128
  return getEigenvalues( diag, offdiag )[0];
 
129
}
 
130
 
 
131
/**
 
132
 * use the above if more than one eigenvector is required and only use this if you need exactly one value
 
133
 */
 
134
float getMajorEigenvalueOld( vec3 diag, vec3 offdiag )
 
135
{
 
136
  const float PiThird = 3.14159265358979323846264/3.;
 
137
  float I1 = diag.x + diag.y + diag.z;
 
138
  float I2 = diag.x * diag.y + diag.x * diag.z + diag.y * diag.z - dot( offdiag, offdiag );
 
139
  float I3 = diag.x * diag.y * diag.z + 2. * offdiag.x * offdiag.y * offdiag.z
 
140
      - ( diag.x * offdiag.x * offdiag.x + diag.y * offdiag.y * offdiag.y + diag.z * offdiag.z * offdiag.z );
 
141
 
 
142
  const float third = 1. / 3.;
 
143
  float I1third = I1 * third;
 
144
  float I1thirdsqr = I1third * I1third;
 
145
  float I2third = I2 * third;
 
146
  float v  = I1thirdsqr - I2third;
 
147
 
 
148
  vec3 lambda;
 
149
 
 
150
  float s  = I1thirdsqr * I1third - I1 * I2 / 6. + 0.5 * I3;
 
151
 
 
152
  // for real eigenvalues: v>0, s^2 < v^3
 
153
  float sqrtv = sqrt( v );
 
154
  float phi= acos( s / ( v * sqrtv ) ) * third;
 
155
 
 
156
  float sqrtv2 = 2. * sqrtv;
 
157
 
 
158
  // due to the cosine function and the fact, that 0 <= phi <= pi/3
 
159
  // it is obvious that the eigenvalues need no further sorting
 
160
  return I1third + sqrtv2 * cos( phi );
 
161
}
 
162
 
 
163
// compute vector direction depending on information computed by getEigenvalues
 
164
// before (Hasan et al. 2001)
 
165
vec3 getEigenvector( vec3 ABC /*diag without eigenalue i*/, vec3 offdiag )
 
166
{
 
167
  vec3 vec;
 
168
  vec.x = ( offdiag.z * offdiag.x - ABC.y * offdiag.y ) * ( offdiag.y * offdiag.x - ABC.z * offdiag.z ); // FIXME
 
169
  //< last component is missing in the paper! there is only a Dx?
 
170
  vec.y = ( offdiag.y * offdiag.x - ABC.z * offdiag.z ) * ( offdiag.y * offdiag.z - ABC.x * offdiag.x );
 
171
  vec.z = ( offdiag.z * offdiag.x - ABC.y * offdiag.y ) * ( offdiag.y * offdiag.z - ABC.x * offdiag.x );
 
172
 
 
173
  return normalize( vec );
 
174
}
 
175
 
 
176
// FA as used everywhere (Kondratieva et al 2005)
 
177
// computes FA from eigenvalues
 
178
float getFA( vec3 evalues )
 
179
{
 
180
  float sqrt32 = sqrt( 3. / 2. ); // would make this const, if compiler let me
 
181
  float I1 = evalues.x + evalues.y + evalues.z;
 
182
  float mu = I1 / 3.;
 
183
  vec3 deriv;
 
184
  deriv.x = ( evalues.x - mu );
 
185
  deriv.y = ( evalues.y - mu );
 
186
  deriv.z = ( evalues.z - mu );
 
187
  float FA = sqrt32 * length( deriv ) / length( evalues );
 
188
  return FA;
 
189
}
 
190
 
 
191
// Compute FA without doing eigenvalue decomposition
 
192
float getFA( vec3 diag, vec3 offdiag )
 
193
{
 
194
  float cross = dot( offdiag, offdiag );
 
195
  float j2 = diag.x * diag.y + diag.y * diag.z + diag.z * diag.x - cross;
 
196
  float j4 = dot( diag, diag ) + 2. * cross;
 
197
  return sqrt( ( j4 - j2 ) / j4 );
 
198
}
 
199
 
 
200
 
 
201
// simple default color map, but with precomputed FA values
 
202
vec4 getColor( vec3 dir, float FA )
 
203
{
 
204
  vec4 color;
 
205
  //dir = normalize( dir );
 
206
  //dir = clamp( dir, 1., 1.);
 
207
  color = vec4( abs( dir ) * FA, 1. );
 
208
  return color;
 
209
}
 
210
 
 
211
// compute the default RGB color scheme
 
212
vec4 getColor( vec3 evec, vec3 evalues )
 
213
{
 
214
  float fa = getFA( evalues );
 
215
  return getColor( evec, fa );
 
216
}
 
217
 
 
218
// make things transparent, proportional to FA
 
219
vec4 getColorAlpha( vec3 evec, vec3 evalues )
 
220
{
 
221
  float fa = getFA( evalues );
 
222
  //evec = clamp( evec, 1., 1.);
 
223
  return vec4( abs( evec ), 1. ) * fa;
 
224
}
 
225
 
 
226
#endif // WGETENSORTOOLS_GLSL
 
227