1 /*=========================================================================
2
3 Program: Visualization Toolkit
4 Module: vtkUnstructuredGridHomogeneousRayIntegrator.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 "vtkUnstructuredGridHomogeneousRayIntegrator.h"
26
27 #include "vtkAbstractVolumeMapper.h"
28 #include "vtkColorTransferFunction.h"
29 #include "vtkDoubleArray.h"
30 #include "vtkObjectFactory.h"
31 #include "vtkPiecewiseFunction.h"
32 #include "vtkUnstructuredGrid.h"
33 #include "vtkVolume.h"
34 #include "vtkVolumeProperty.h"
35
36 #include <cmath>
37
38 //------------------------------------------------------------------------------
39
40 vtkStandardNewMacro(vtkUnstructuredGridHomogeneousRayIntegrator);
41
42 //------------------------------------------------------------------------------
43
vtkUnstructuredGridHomogeneousRayIntegrator()44 vtkUnstructuredGridHomogeneousRayIntegrator::vtkUnstructuredGridHomogeneousRayIntegrator()
45 {
46 this->Property = nullptr;
47
48 this->NumComponents = 0;
49 this->ColorTable = nullptr;
50 this->AttenuationTable = nullptr;
51 this->TableShift = nullptr;
52 this->TableScale = nullptr;
53
54 this->UseAverageColor = 0;
55 this->TransferFunctionTableSize = 1024;
56 }
57
~vtkUnstructuredGridHomogeneousRayIntegrator()58 vtkUnstructuredGridHomogeneousRayIntegrator::~vtkUnstructuredGridHomogeneousRayIntegrator()
59 {
60 for (int i = 0; i < this->NumComponents; i++)
61 {
62 delete[] this->ColorTable[i];
63 delete[] this->AttenuationTable[i];
64 }
65 delete[] this->ColorTable;
66 delete[] this->AttenuationTable;
67 delete[] this->TableShift;
68 delete[] this->TableScale;
69 }
70
PrintSelf(ostream & os,vtkIndent indent)71 void vtkUnstructuredGridHomogeneousRayIntegrator::PrintSelf(ostream& os, vtkIndent indent)
72 {
73 this->Superclass::PrintSelf(os, indent);
74
75 os << indent << "UseAverageColor: " << this->UseAverageColor << endl;
76 os << indent << "TransferFunctionTableSize: " << this->TransferFunctionTableSize << endl;
77 }
78
79 //------------------------------------------------------------------------------
80
GetTransferFunctionTables(vtkDataArray * scalars)81 void vtkUnstructuredGridHomogeneousRayIntegrator::GetTransferFunctionTables(vtkDataArray* scalars)
82 {
83 for (int i = 0; i < this->NumComponents; i++)
84 {
85 delete[] this->ColorTable[i];
86 delete[] this->AttenuationTable[i];
87 }
88 delete[] this->ColorTable;
89 delete[] this->AttenuationTable;
90 delete[] this->TableShift;
91 delete[] this->TableScale;
92
93 this->NumComponents = scalars->GetNumberOfComponents();
94
95 this->ColorTable = new float*[this->NumComponents];
96 this->AttenuationTable = new float*[this->NumComponents];
97 this->TableShift = new double[this->NumComponents];
98 this->TableScale = new double[this->NumComponents];
99
100 for (int c = 0; c < this->NumComponents; c++)
101 {
102 double range[2];
103 scalars->GetRange(range, c);
104 if (range[0] >= range[1])
105 {
106 range[1] = range[0] + 1;
107 }
108 this->TableScale[c] = this->TransferFunctionTableSize / (range[1] - range[0]);
109 this->TableShift[c] = -range[0] * this->TransferFunctionTableSize / (range[1] - range[0]);
110
111 this->ColorTable[c] = new float[3 * this->TransferFunctionTableSize];
112 if (this->Property->GetColorChannels(c) == 1)
113 {
114 // Get gray values. Store temporarily in allocated RGB array.
115 this->Property->GetGrayTransferFunction(c)->GetTable(
116 range[0], range[1], this->TransferFunctionTableSize, this->ColorTable[c]);
117 // Convert gray into RGB. Copy backward so that we can use the same
118 // array.
119 for (int i = this->TransferFunctionTableSize - 1; i >= 0; i--)
120 {
121 this->ColorTable[c][3 * i + 0] = this->ColorTable[c][3 * i + 1] =
122 this->ColorTable[c][3 * i + 2] = this->ColorTable[c][i];
123 }
124 }
125 else
126 {
127 this->Property->GetRGBTransferFunction(c)->GetTable(
128 range[0], range[1], this->TransferFunctionTableSize, this->ColorTable[c]);
129 }
130
131 this->AttenuationTable[c] = new float[this->TransferFunctionTableSize];
132 this->Property->GetScalarOpacity(c)->GetTable(
133 range[0], range[1], this->TransferFunctionTableSize, this->AttenuationTable[c]);
134
135 // Adjust attenuation by scalar unit length. This will make the unit
136 // length the same as the model.
137 float unitlength = this->Property->GetScalarOpacityUnitDistance(c);
138 for (int i = 0; i < this->TransferFunctionTableSize; i++)
139 {
140 this->AttenuationTable[c][i] /= unitlength;
141 }
142 }
143
144 this->TablesBuilt.Modified();
145 }
146
147 //------------------------------------------------------------------------------
148
Initialize(vtkVolume * volume,vtkDataArray * scalars)149 void vtkUnstructuredGridHomogeneousRayIntegrator::Initialize(
150 vtkVolume* volume, vtkDataArray* scalars)
151 {
152 vtkVolumeProperty* property = volume->GetProperty();
153
154 if ((property == this->Property) && (this->TablesBuilt > property->GetMTime()) &&
155 (this->TablesBuilt > this->MTime))
156 {
157 // Nothing changed from the last time Initialize was run.
158 return;
159 }
160
161 this->Property = property;
162 this->Volume = volume;
163
164 if (property->GetIndependentComponents())
165 {
166 this->GetTransferFunctionTables(scalars);
167 }
168 }
169
170 //------------------------------------------------------------------------------
171
Integrate(vtkDoubleArray * intersectionLengths,vtkDataArray * nearIntersections,vtkDataArray * vtkNotUsed (farIntersections),float color[4])172 void vtkUnstructuredGridHomogeneousRayIntegrator::Integrate(vtkDoubleArray* intersectionLengths,
173 vtkDataArray* nearIntersections, vtkDataArray* vtkNotUsed(farIntersections), float color[4])
174 {
175 vtkIdType numIntersections = intersectionLengths->GetNumberOfTuples();
176
177 if (this->Property->GetIndependentComponents())
178 {
179 if (this->NumComponents == 1)
180 {
181 // Optimize for what I think is one of the most common uses.
182 for (vtkIdType i = 0; i < numIntersections; i++)
183 {
184 int table_index =
185 (int)(this->TableScale[0] * nearIntersections->GetComponent(i, 0) + this->TableShift[0]);
186 if (table_index < 0)
187 {
188 table_index = 0;
189 }
190 if (table_index >= this->TransferFunctionTableSize)
191 {
192 table_index = this->TransferFunctionTableSize - 1;
193 }
194 float* c = this->ColorTable[0] + 3 * table_index;
195 float tau = this->AttenuationTable[0][table_index];
196 float alpha = 1 - (float)exp(-intersectionLengths->GetComponent(i, 0) * tau);
197 color[0] += c[0] * alpha * (1 - color[3]);
198 color[1] += c[1] * alpha * (1 - color[3]);
199 color[2] += c[2] * alpha * (1 - color[3]);
200 color[3] += alpha * (1 - color[3]);
201 }
202 }
203 else
204 {
205 // Generic case.
206 for (vtkIdType i = 0; i < numIntersections; i++)
207 {
208 float newcolor[4];
209 int table_index =
210 (int)(this->TableScale[0] * nearIntersections->GetComponent(i, 0) + this->TableShift[0]);
211 if (table_index < 0)
212 {
213 table_index = 0;
214 }
215 if (table_index >= this->TransferFunctionTableSize)
216 {
217 table_index = this->TransferFunctionTableSize - 1;
218 }
219 float* c = this->ColorTable[0] + 3 * table_index;
220 float tau = this->AttenuationTable[0][table_index];
221 newcolor[0] = c[0];
222 newcolor[1] = c[1];
223 newcolor[2] = c[2];
224 newcolor[3] = tau;
225 for (int component = 1; component < this->NumComponents; component++)
226 {
227 table_index =
228 (int)(this->TableScale[component] * nearIntersections->GetComponent(i, component) +
229 this->TableShift[component]);
230 if (table_index < 0)
231 {
232 table_index = 0;
233 }
234 if (table_index >= this->TransferFunctionTableSize)
235 {
236 table_index = this->TransferFunctionTableSize - 1;
237 }
238 c = this->ColorTable[component] + 3 * table_index;
239 tau = this->AttenuationTable[component][table_index];
240 // Here we handle the mixing of material properties. This never
241 // seems to be defined very clearly. I handle this by assuming
242 // that each scalar represents a cloud of particles of a certain
243 // color and a certain density. We mix the scalars in the same
244 // way as mixing these particles together. By necessity, the
245 // density becomes greater. The "opacity" parameter is really
246 // interpreted as the attenuation coefficient (which is
247 // proportional to density) and can therefore easily be greater
248 // than one. The opacity of the resulting color will, however,
249 // always be scaled between 0 and 1.
250 if (tau + newcolor[3] > 1.0e-8f)
251 {
252 newcolor[0] *= newcolor[3] / (tau + newcolor[3]);
253 newcolor[1] *= newcolor[3] / (tau + newcolor[3]);
254 newcolor[2] *= newcolor[3] / (tau + newcolor[3]);
255 newcolor[0] += c[0] * tau / (tau + newcolor[3]);
256 newcolor[1] += c[1] * tau / (tau + newcolor[3]);
257 newcolor[2] += c[2] * tau / (tau + newcolor[3]);
258 newcolor[3] += tau;
259 }
260 }
261 float alpha = 1 - (float)exp(-intersectionLengths->GetComponent(i, 0) * newcolor[3]);
262 color[0] += newcolor[0] * alpha * (1 - color[3]);
263 color[1] += newcolor[1] * alpha * (1 - color[3]);
264 color[2] += newcolor[2] * alpha * (1 - color[3]);
265 color[3] += alpha * (1 - color[3]);
266 }
267 }
268 }
269 else
270 {
271 int numComponents = nearIntersections->GetNumberOfComponents();
272 for (vtkIdType i = 0; i < numIntersections; i++)
273 {
274 double c[4];
275 if (numComponents == 4)
276 {
277 nearIntersections->GetTuple(i, c);
278 }
279 else
280 {
281 double* lt = nearIntersections->GetTuple(i);
282 c[0] = c[1] = c[2] = lt[0];
283 c[3] = lt[1];
284 }
285 float alpha = 1 - (float)exp(-intersectionLengths->GetComponent(i, 0) * c[3]);
286 color[0] += (float)c[0] * alpha * (1 - color[3]);
287 color[1] += (float)c[1] * alpha * (1 - color[3]);
288 color[2] += (float)c[2] * alpha * (1 - color[3]);
289 color[3] += alpha * (1 - color[3]);
290 }
291 }
292 }
293