1 /*=========================================================================
2
3 Program: Visualization Toolkit
4 Module: vtkUnstructuredGridLinearRayIntegrator.cxx
5
6 Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
7 All rights reserved.
8 See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
9
10 This software is distributed WITHOUT ANY WARRANTY; without even
11 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
12 PURPOSE. See the above copyright notice for more information.
13
14 =========================================================================*/
15
16 /*
17 * Copyright 2004 Sandia Corporation.
18 * Under the terms of Contract DE-AC04-94AL85000, there is a non-exclusive
19 * license for use of this work by or on behalf of the
20 * U.S. Government. Redistribution and use in source and binary forms, with
21 * or without modification, are permitted provided that this Notice and any
22 * statement of authorship are reproduced on all copies.
23 */
24
25 #include "vtkUnstructuredGridLinearRayIntegrator.h"
26
27 #include "vtkObjectFactory.h"
28 #include "vtkVolumeProperty.h"
29 #include "vtkVolume.h"
30 #include "vtkDoubleArray.h"
31 #include "vtkPiecewiseFunction.h"
32 #include "vtkColorTransferFunction.h"
33 #include "vtkMath.h"
34
35 #include <vector>
36 #include <set>
37 #include <algorithm>
38
39 #include <math.h>
40
41 #ifndef M_SQRTPI
42 #define M_SQRTPI 1.77245385090551602792981
43 #endif
44 #ifndef M_2_SQRTPI
45 #define M_2_SQRTPI 1.12837916709551257390
46 #endif
47 #ifndef M_1_SQRTPI
48 #define M_1_SQRTPI (0.5*M_2_SQRTPI)
49 #endif
50
51 //-----------------------------------------------------------------------------
52
53 // VTK's native classes for defining transfer functions is actually slow to
54 // access, so we have to cache it somehow. This class is straightforward
55 // copy of the transfer function.
56 class vtkLinearRayIntegratorTransferFunction
57 {
58 public:
59 vtkLinearRayIntegratorTransferFunction();
60 ~vtkLinearRayIntegratorTransferFunction();
61
62 void GetTransferFunction(vtkColorTransferFunction *color,
63 vtkPiecewiseFunction *opacity,
64 double unit_distance,
65 double scalar_range[2]);
66 void GetTransferFunction(vtkPiecewiseFunction *intensity,
67 vtkPiecewiseFunction *opacity,
68 double unit_distance,
69 double scalar_range[2]);
70
71 inline void GetColor(double x, double c[4]);
72
73 struct acolor {
74 double c[4];
75 };
76 double *ControlPoints;
77 int NumControlPoints;
78 acolor *Colors;
79
80 private:
81 vtkLinearRayIntegratorTransferFunction(const vtkLinearRayIntegratorTransferFunction&); // Not implemented.
82 void operator=(const vtkLinearRayIntegratorTransferFunction &); // Not implemented.
83 };
84
vtkLinearRayIntegratorTransferFunction()85 vtkLinearRayIntegratorTransferFunction::vtkLinearRayIntegratorTransferFunction()
86 {
87 this->ControlPoints = NULL;
88 this->Colors = NULL;
89
90 this->NumControlPoints = 0;
91 }
92
~vtkLinearRayIntegratorTransferFunction()93 vtkLinearRayIntegratorTransferFunction::~vtkLinearRayIntegratorTransferFunction()
94 {
95 delete[] this->ControlPoints;
96 delete[] this->Colors;
97 }
98
99 static const double huebends[6] = {
100 1.0/6.0, 1.0/3.0, 0.5, 2.0/3.0, 5.0/6.0, 1.0
101 };
102
GetTransferFunction(vtkColorTransferFunction * color,vtkPiecewiseFunction * opacity,double unit_distance,double scalar_range[2])103 void vtkLinearRayIntegratorTransferFunction::GetTransferFunction(
104 vtkColorTransferFunction *color,
105 vtkPiecewiseFunction *opacity,
106 double unit_distance,
107 double scalar_range[2])
108 {
109 std::set<double> cpset;
110
111 double *function_range = color->GetRange();
112 double *function = color->GetDataPointer();
113 while (1)
114 {
115 cpset.insert(function[0]);
116 if (function[0] == function_range[1]) break;
117 function += 4;
118 }
119
120 if (color->GetColorSpace() != VTK_CTF_RGB)
121 {
122 // If we are in an HSV color space, we must insert control points
123 // in places where the RGB bends.
124 double rgb[3], hsv[3];
125 double hue1, hue2;
126 double x1, x2;
127 std::set<double>::iterator i = cpset.begin();
128 x1 = *i;
129 color->GetColor(x1, rgb);
130 vtkMath::RGBToHSV(rgb, hsv);
131 hue1 = hsv[0];
132 for (i++; i != cpset.end(); i++)
133 {
134 x2 = *i;
135 color->GetColor(x2, rgb);
136 vtkMath::RGBToHSV(rgb, hsv);
137 hue2 = hsv[0];
138
139 // Are we crossing the 0/1 boundary?
140 if ( (color->GetColorSpace() == VTK_CTF_HSV && color->GetHSVWrap() )
141 && ((hue1 - hue2 > 0.5) || (hue2 - hue1 > 0.5)) )
142 {
143 // Yes, we are crossing the boundary.
144 if (hue1 > hue2)
145 {
146 int j;
147 for (j = 0; huebends[j] <= hue2; j++)
148 {
149 double interp = (1-hue1+huebends[j])/(1-hue1+hue2);
150 cpset.insert((x2-x1)*interp + x1);
151 }
152 while (huebends[j] < hue1) j++;
153 for ( ; j < 6; j++)
154 {
155 double interp = (huebends[j]-hue1)/(1-hue1+hue2);
156 cpset.insert((x2-x1)*interp + x1);
157 }
158 }
159 else
160 {
161 int j;
162 for (j = 0; huebends[j] <= hue1; j++)
163 {
164 double interp = (hue1-huebends[j])/(1-hue2+hue1);
165 cpset.insert((x2-x1)*interp + x1);
166 }
167 while (huebends[j] < hue2) j++;
168 for ( ; j < 6; j++)
169 {
170 double interp = (1-huebends[j]+hue1)/(1-hue2+hue1);
171 cpset.insert((x2-x1)*interp + x1);
172 }
173 }
174 }
175 else
176 {
177 // No, we are not crossing the boundary.
178 int j = 0;
179 double minh, maxh;
180 if (hue1 < hue2)
181 {
182 minh = hue1; maxh = hue2;
183 }
184 else
185 {
186 minh = hue2; maxh = hue1;
187 }
188 while (huebends[j] < minh) j++;
189 for (j = 0; huebends[j] < maxh; j++)
190 {
191 double interp = (huebends[j]-hue1)/(hue2-hue1);
192 cpset.insert((x2-x1)*interp + x1);
193 }
194 }
195
196 x1 = x2;
197 hue1 = hue2;
198 }
199 }
200
201 function_range = opacity->GetRange();
202 function = opacity->GetDataPointer();
203 while (1)
204 {
205 cpset.insert(function[0]);
206 if (function[0] == function_range[0]) break;
207 function += 2;
208 }
209
210 // Add the scalar at the beginning and end of the range so the interpolation
211 // is correct there.
212 cpset.insert(scalar_range[0]);
213 cpset.insert(scalar_range[1]);
214 // Make extra sure there are at least two entries in cpset.
215 if (cpset.size() < 2)
216 {
217 cpset.insert(0.0);
218 cpset.insert(1.0);
219 }
220
221 // Now record control points and colors.
222 delete[] this->ControlPoints;
223 delete[] this->Colors;
224 this->NumControlPoints = static_cast<int>(cpset.size());
225 this->ControlPoints = new double[this->NumControlPoints];
226 this->Colors = new acolor[this->NumControlPoints];
227
228 std::copy(cpset.begin(), cpset.end(), this->ControlPoints);
229 for (int i = 0; i < this->NumControlPoints; i++)
230 {
231 color->GetColor(this->ControlPoints[i], this->Colors[i].c);
232 this->Colors[i].c[3] = ( opacity->GetValue(this->ControlPoints[i])
233 / unit_distance);
234 }
235 }
236
GetTransferFunction(vtkPiecewiseFunction * intensity,vtkPiecewiseFunction * opacity,double unit_distance,double scalar_range[2])237 void vtkLinearRayIntegratorTransferFunction::GetTransferFunction(
238 vtkPiecewiseFunction *intensity,
239 vtkPiecewiseFunction *opacity,
240 double unit_distance,
241 double scalar_range[2])
242 {
243 std::set<double> cpset;
244
245 double *function_range = intensity->GetRange();
246 double *function = intensity->GetDataPointer();
247 while (1)
248 {
249 cpset.insert(function[0]);
250 if (function[0] == function_range[1]) break;
251 function += 2;
252 }
253
254 function_range = opacity->GetRange();
255 function = opacity->GetDataPointer();
256 while (1)
257 {
258 cpset.insert(function[0]);
259 if (function[0] == function_range[0]) break;
260 function += 2;
261 }
262
263 // Add the scalar at the beginning and end of the range so the interpolation
264 // is correct there.
265 cpset.insert(scalar_range[0]);
266 cpset.insert(scalar_range[1]);
267 // Make extra sure there are at least two entries in cpset.
268 if (cpset.size() < 2)
269 {
270 cpset.insert(0.0);
271 cpset.insert(1.0);
272 }
273
274 // Now record control points and colors.
275 delete[] this->ControlPoints;
276 delete[] this->Colors;
277 this->NumControlPoints = static_cast<int>(cpset.size());
278 this->ControlPoints = new double[this->NumControlPoints];
279 this->Colors = new acolor[this->NumControlPoints];
280
281 std::copy(cpset.begin(), cpset.end(), this->ControlPoints);
282 for (int i = 0; i < this->NumControlPoints; i++)
283 {
284 // Is setting all the colors to the same value the right thing to do?
285 this->Colors[i].c[0] = this->Colors[i].c[1] = this->Colors[i].c[2]
286 = intensity->GetValue(this->ControlPoints[i]);
287 this->Colors[i].c[3] = ( opacity->GetValue(this->ControlPoints[i])
288 / unit_distance);
289 }
290 }
291
GetColor(double x,double c[4])292 inline void vtkLinearRayIntegratorTransferFunction::GetColor(double x,
293 double c[4])
294 {
295 int i = 1;
296 while ((this->ControlPoints[i] < x) && (i < this->NumControlPoints-1)) i++;
297
298 double before = this->ControlPoints[i-1];
299 double after = this->ControlPoints[i];
300
301 double interp = (x-before)/(after-before);
302
303 double *beforec = this->Colors[i-1].c;
304 double *afterc = this->Colors[i].c;
305 c[0] = (1-interp)*beforec[0] + interp*afterc[0];
306 c[1] = (1-interp)*beforec[1] + interp*afterc[1];
307 c[2] = (1-interp)*beforec[2] + interp*afterc[2];
308 c[3] = (1-interp)*beforec[3] + interp*afterc[3];
309 }
310
311 //-----------------------------------------------------------------------------
312
313 vtkStandardNewMacro(vtkUnstructuredGridLinearRayIntegrator);
314
vtkUnstructuredGridLinearRayIntegrator()315 vtkUnstructuredGridLinearRayIntegrator::vtkUnstructuredGridLinearRayIntegrator()
316 {
317 this->Property = NULL;
318 this->TransferFunctions = NULL;
319 this->NumIndependentComponents = 0;
320 }
321
322 //-----------------------------------------------------------------------------
323
~vtkUnstructuredGridLinearRayIntegrator()324 vtkUnstructuredGridLinearRayIntegrator::~vtkUnstructuredGridLinearRayIntegrator()
325 {
326 delete[] this->TransferFunctions;
327 }
328
329 //-----------------------------------------------------------------------------
330
PrintSelf(ostream & os,vtkIndent indent)331 void vtkUnstructuredGridLinearRayIntegrator::PrintSelf(ostream &os,
332 vtkIndent indent)
333 {
334 this->Superclass::PrintSelf(os, indent);
335 }
336
337 //-----------------------------------------------------------------------------
338
Initialize(vtkVolume * volume,vtkDataArray * scalars)339 void vtkUnstructuredGridLinearRayIntegrator::Initialize(
340 vtkVolume *volume,
341 vtkDataArray *scalars)
342 {
343 vtkVolumeProperty *property = volume->GetProperty();
344
345 if ( (property == this->Property)
346 && (this->TransferFunctionsModified > property->GetMTime()) )
347 {
348 // Nothing has changed from the last time Initialize was run.
349 return;
350 }
351
352 int numcomponents = scalars->GetNumberOfComponents();
353
354 this->Property = property;
355 this->TransferFunctionsModified.Modified();
356
357 if (!property->GetIndependentComponents())
358 {
359 // The scalars actually hold material properties.
360 if ((numcomponents != 4) && (numcomponents != 2) )
361 {
362 vtkErrorMacro("Only 2-tuples and 4-tuples allowed for dependent components.");
363 }
364 return;
365 }
366
367 delete[] this->TransferFunctions;
368
369 this->NumIndependentComponents = numcomponents;
370 this->TransferFunctions
371 = new vtkLinearRayIntegratorTransferFunction[numcomponents];
372
373 for (int component = 0; component < numcomponents; component++)
374 {
375 if (property->GetColorChannels(component) == 1)
376 {
377 this->TransferFunctions[component]
378 .GetTransferFunction(property->GetGrayTransferFunction(component),
379 property->GetScalarOpacity(component),
380 property->GetScalarOpacityUnitDistance(component),
381 scalars->GetRange(component));
382 }
383 else
384 {
385 this->TransferFunctions[component]
386 .GetTransferFunction(property->GetRGBTransferFunction(component),
387 property->GetScalarOpacity(component),
388 property->GetScalarOpacityUnitDistance(component),
389 scalars->GetRange(component));
390 }
391 }
392 }
393
394 //-----------------------------------------------------------------------------
395
Integrate(vtkDoubleArray * intersectionLengths,vtkDataArray * nearIntersections,vtkDataArray * farIntersections,float color[4])396 void vtkUnstructuredGridLinearRayIntegrator::Integrate(
397 vtkDoubleArray *intersectionLengths,
398 vtkDataArray *nearIntersections,
399 vtkDataArray *farIntersections,
400 float color[4])
401 {
402 int numintersections = intersectionLengths->GetNumberOfTuples();
403 if (this->Property->GetIndependentComponents())
404 {
405 int numscalars = nearIntersections->GetNumberOfComponents();
406 double *nearScalars = new double[numscalars];
407 double *farScalars = new double[numscalars];
408 std::set<double> segments;
409 for (vtkIdType i = 0; i < numintersections; i++)
410 {
411 double total_length = intersectionLengths->GetValue(i);
412 nearIntersections->GetTuple(i, nearScalars);
413 farIntersections->GetTuple(i, farScalars);
414
415 // Split up segment on control points, because it is nonlinear in
416 // these regions.
417 segments.erase(segments.begin(), segments.end());
418 segments.insert(0.0);
419 segments.insert(1.0);
420 for (int j = 0; j < numscalars; j++)
421 {
422 double *cp = this->TransferFunctions[j].ControlPoints;
423 vtkIdType numcp = this->TransferFunctions[j].NumControlPoints;
424 double minscalar, maxscalar;
425 if (nearScalars[j] < farScalars[j])
426 {
427 minscalar = nearScalars[j]; maxscalar = farScalars[j];
428 }
429 else
430 {
431 minscalar = farScalars[j]; maxscalar = nearScalars[j];
432 }
433 for (int k = 0; k < numcp; k++)
434 {
435 if (cp[k] <= minscalar) continue;
436 if (cp[k] >= maxscalar) break;
437 // If we are here, we need to break the segment at the given scalar.
438 // Find the fraction between the near and far segment points.
439 segments.insert( (cp[k]-nearScalars[j])
440 / (farScalars[j]-nearScalars[j]));
441 }
442 }
443
444 // Iterate over all the segment pieces (from front to back) and
445 // integrate each piece.
446 std::set<double>::iterator segi = segments.begin();
447 double nearInterpolant = *segi;
448 for (segi++; segi != segments.end(); segi++)
449 {
450 double farInterpolant = *segi;
451 double nearcolor[4] = {0.0, 0.0, 0.0, 0.0};
452 double farcolor[4] = {0.0, 0.0, 0.0, 0.0};
453 double length = total_length*(farInterpolant-nearInterpolant);
454 // Here we handle the mixing of material properties. This never
455 // seems to be defined very clearly. I handle this by assuming
456 // that each scalar represents a cloud of particles of a certain
457 // color and a certain density. We mix the scalars in the same way
458 // as mixing these particles together. By necessity, the density
459 // becomes greater. The "opacity" parameter is really interpreted
460 // as the attenuation coefficient (which is proportional to
461 // density) and can therefore easily be greater than one. The
462 // opacity of the resulting color will, however, always be scaled
463 // between 0 and 1.
464 for (int j = 0; j < numscalars; j++)
465 {
466 double scalar
467 = (farScalars[j]-nearScalars[j])*nearInterpolant + nearScalars[j];
468 if (j == 0)
469 {
470 this->TransferFunctions[j].GetColor(scalar, nearcolor);
471 }
472 else
473 {
474 double c[4];
475 this->TransferFunctions[j].GetColor(scalar, c);
476 if (c[3] + nearcolor[3] > 1.0e-8f)
477 {
478 nearcolor[0] *= nearcolor[3]/(c[3] + nearcolor[3]);
479 nearcolor[1] *= nearcolor[3]/(c[3] + nearcolor[3]);
480 nearcolor[2] *= nearcolor[3]/(c[3] + nearcolor[3]);
481 nearcolor[0] += c[0]*c[3]/(c[3] + nearcolor[3]);
482 nearcolor[1] += c[1]*c[3]/(c[3] + nearcolor[3]);
483 nearcolor[2] += c[2]*c[3]/(c[3] + nearcolor[3]);
484 nearcolor[3] += c[3];
485 }
486 }
487
488 scalar
489 = (farScalars[j]-nearScalars[j])*farInterpolant + nearScalars[j];
490 if (j == 0)
491 {
492 this->TransferFunctions[j].GetColor(scalar, farcolor);
493 }
494 else
495 {
496 double c[4];
497 this->TransferFunctions[j].GetColor(scalar, c);
498 if (c[3] + farcolor[3] > 1.0e-8f)
499 {
500 farcolor[0] *= farcolor[3]/(c[3] + farcolor[3]);
501 farcolor[1] *= farcolor[3]/(c[3] + farcolor[3]);
502 farcolor[2] *= farcolor[3]/(c[3] + farcolor[3]);
503 farcolor[0] += c[0]*c[3]/(c[3] + farcolor[3]);
504 farcolor[1] += c[1]*c[3]/(c[3] + farcolor[3]);
505 farcolor[2] += c[2]*c[3]/(c[3] + farcolor[3]);
506 farcolor[3] += c[3];
507 }
508 }
509 }
510 this->IntegrateRay(length, nearcolor, nearcolor[3],
511 farcolor, farcolor[3], color);
512
513 nearInterpolant = farInterpolant;
514 }
515 }
516 delete[] nearScalars;
517 delete[] farScalars;
518 }
519 else
520 {
521 double unitdistance = this->Property->GetScalarOpacityUnitDistance();
522 if (nearIntersections->GetNumberOfComponents() == 4)
523 {
524 for (vtkIdType i = 0; i < numintersections; i++)
525 {
526 double length = intersectionLengths->GetValue(i);
527 double *nearcolor = nearIntersections->GetTuple(i);
528 double *farcolor = farIntersections->GetTuple(i);
529 this->IntegrateRay(length, nearcolor, nearcolor[3]/unitdistance,
530 farcolor, farcolor[3]/unitdistance, color);
531 }
532 }
533 else // Two components.
534 {
535 for (vtkIdType i = 0; i < numintersections; i++)
536 {
537 double length = intersectionLengths->GetValue(i);
538 double *nearcolor = nearIntersections->GetTuple(i);
539 double *farcolor = farIntersections->GetTuple(i);
540 this->IntegrateRay(length, nearcolor[0], nearcolor[1]/unitdistance,
541 farcolor[0], farcolor[1]/unitdistance, color);
542 }
543 }
544 }
545 }
546
547 //-----------------------------------------------------------------------------
548
IntegrateRay(double length,double intensity_front,double attenuation_front,double intensity_back,double attenuation_back,float color[4])549 void vtkUnstructuredGridLinearRayIntegrator::IntegrateRay(
550 double length,
551 double intensity_front,
552 double attenuation_front,
553 double intensity_back,
554 double attenuation_back,
555 float color[4])
556 {
557 float Psi = vtkUnstructuredGridLinearRayIntegrator::Psi(length,
558 attenuation_front,
559 attenuation_back);
560 float zeta = (float)exp(-0.5*length*(attenuation_front+attenuation_back));
561 float alpha = 1-zeta;
562
563 float newintensity = (1-color[3])*( intensity_front*(1-Psi)
564 + intensity_back*(Psi-zeta) );
565 // Is setting the RGB values the same the right thing to do?
566 color[0] += newintensity;
567 color[1] += newintensity;
568 color[2] += newintensity;
569 color[3] += (1-color[3])*alpha;
570 }
571
IntegrateRay(double length,const double color_front[3],double attenuation_front,const double color_back[3],double attenuation_back,float color[4])572 void vtkUnstructuredGridLinearRayIntegrator::IntegrateRay(
573 double length,
574 const double color_front[3],
575 double attenuation_front,
576 const double color_back[3],
577 double attenuation_back,
578 float color[4])
579 {
580 float Psi = vtkUnstructuredGridLinearRayIntegrator::Psi(length,
581 attenuation_front,
582 attenuation_back);
583 float zeta = (float)exp(-0.5*length*(attenuation_front+attenuation_back));
584 float alpha = 1-zeta;
585
586 color[0] += (1-color[3])*(color_front[0]*(1-Psi) + color_back[0]*(Psi-zeta));
587 color[1] += (1-color[3])*(color_front[1]*(1-Psi) + color_back[1]*(Psi-zeta));
588 color[2] += (1-color[3])*(color_front[2]*(1-Psi) + color_back[2]*(Psi-zeta));
589 color[3] += (1-color[3])*alpha;
590 }
591
592 //-----------------------------------------------------------------------------
593
erf_fitting_function(float u)594 static inline float erf_fitting_function(float u)
595 {
596 return
597 - 1.26551223 + u*(1.00002368 + u*(0.37409196 + u*(0.09678418 +
598 u*(-0.18628806 + u*(0.27886807 + u*(-1.13520398 + u*(1.48851587 +
599 u*(-0.82215223 + u*0.17087277))))))));
600 }
601
602 #if 0
603 // This function is not used directly. It is here for reference.
604 static inline float erf(float x)
605 {
606 /* Compute as described in Numerical Recipes in C++ by Press, et al. */
607 /* x = abs(x); In this application, x should always be >= 0. */
608 float u = 1/(1 + 0.5*x);
609 float ans = u*exp(-x*x + erf_fitting_function(u));
610 /* return (x >= 0 ? 1 - ans : ans - 1); x should always be >= 0. */
611 return 1 - ans;
612 }
613 #endif
614
615 /* Compute Dawson's integral as described in Numerical Recipes in C++ by
616 Press, et al. */
617 #define H 0.4
618 static const float dawson_constant0 = 0.852144;
619 static const float dawson_constant1 = 0.236928;
620 static const float dawson_constant2 = 0.0183156;
621 static const float dawson_constant3 = 0.000393669;
622 static const float dawson_constant4 = 2.35258e-6;
623 static const float dawson_constant5 = 3.90894e-9;
dawson(float x)624 static inline float dawson(float x)
625 {
626 if (x > 0.2)
627 {
628 /* x = abs(x); In this application, x should always be >= 0. */
629 int n0 = 2*(int)((0.5/H)*x + 0.5);
630 float xp = x - (float)n0*H;
631 float e1 = exp((2*H)*xp);
632 float e2 = e1*e1;
633 float d1 = n0 + 1;
634 float d2 = d1 - 2;
635 float sum = 0;
636 sum = dawson_constant0*(e1/d1 + 1/(d2*e1));
637 d1 += 2; d2 -= 2; e1 *= e2;
638 sum += dawson_constant1*(e1/d1 + 1/(d2*e1));
639 d1 += 2; d2 -= 2; e1 *= e2;
640 sum += dawson_constant2*(e1/d1 + 1/(d2*e1));
641 d1 += 2; d2 -= 2; e1 *= e2;
642 sum += dawson_constant3*(e1/d1 + 1/(d2*e1));
643 d1 += 2; d2 -= 2; e1 *= e2;
644 sum += dawson_constant4*(e1/d1 + 1/(d2*e1));
645 d1 += 2; d2 -= 2; e1 *= e2;
646 sum += dawson_constant5*(e1/d1 + 1/(d2*e1));
647 return M_1_SQRTPI*exp(-xp*xp)*sum;
648 }
649 else
650 {
651 float x2 = x*x;
652 return x*(1 - (2.0/3.0)*x2*(1 - .4*x2*(1 - (2.0/7.0)*x2)));
653 }
654 }
655
656 #if 0
657 // This function is not used directly. It is here for reference.
658 inline float erfi(float x)
659 {
660 return M_2_SQRTPI*exp(x*x)*dawson(x);
661 }
662 #endif
663
Psi(float length,float attenuation_front,float attenuation_back)664 float vtkUnstructuredGridLinearRayIntegrator::Psi(float length,
665 float attenuation_front,
666 float attenuation_back)
667 {
668 float difftauD = length*fabs(attenuation_back - attenuation_front);
669
670 if (difftauD < 1.0e-8f)
671 {
672 // Volume is homogeneous (with respect to attenuation).
673 float tauD = length * attenuation_front;
674 if (tauD < 1.0e-8f)
675 {
676 return 1;
677 }
678 else
679 {
680 return (1 - (float)exp(-tauD))/tauD;
681 }
682 }
683 else
684 {
685 float invsqrt2difftauD = 1/(float)sqrt(2*difftauD);
686 float frontterm = length*invsqrt2difftauD*attenuation_front;
687 float backterm = length*invsqrt2difftauD*attenuation_back;
688 if (attenuation_back > attenuation_front)
689 {
690 float u, Y;
691 u = 1/(1+0.5f*frontterm);
692 Y = u*(float)exp(erf_fitting_function(u));
693 u = 1/(1+0.5f*backterm);
694 Y += -u*exp( frontterm*frontterm-backterm*backterm
695 + erf_fitting_function(u));
696 Y *= M_SQRTPI*invsqrt2difftauD;
697 return Y;
698 }
699 else
700 {
701 float expterm = (float)exp(backterm*backterm-frontterm*frontterm);
702 return 2*invsqrt2difftauD*(dawson(frontterm) - expterm*dawson(backterm));
703 }
704 }
705 }
706