-
sparker
thiago
Manta Interactive Ray Tracer Development Mailing List
Text archives Help
- From: brownlee@sci.utah.edu
- To: manta@sci.utah.edu
- Subject: [Manta] r1792 - in trunk: Model/Materials scenes
- Date: Tue, 30 Oct 2007 02:53:31 -0600 (MDT)
Author: brownlee
Date: Tue Oct 30 02:53:28 2007
New Revision: 1792
Added:
trunk/Model/Materials/Volume.cc
trunk/Model/Materials/Volume.h
trunk/scenes/volumeTest.cc
Modified:
trunk/Model/Materials/CMakeLists.txt
trunk/scenes/CMakeLists.txt
Log:
added basic single ray volume renderer, need TEEM to build it
Modified: trunk/Model/Materials/CMakeLists.txt
==============================================================================
--- trunk/Model/Materials/CMakeLists.txt (original)
+++ trunk/Model/Materials/CMakeLists.txt Tue Oct 30 02:53:28 2007
@@ -32,3 +32,42 @@
Materials/Transparent.cc
Materials/Transparent.h
)
+
+IF(FOUND_TEEM)
+SET (Manta_Materials_SRCS
+ Materials/Volume.h
+ Materials/Volume.cc
+ Materials/AmbientOcclusion.h
+ Materials/AmbientOcclusion.cc
+ Materials/Checker.h
+ Materials/Checker.cc
+ Materials/CopyTextureMaterial.cc
+ Materials/CopyTextureMaterial.h
+ Materials/Dielectric.h
+ Materials/Dielectric.cc
+ Materials/Flat.h
+ Materials/Flat.cc
+ Materials/Lambertian.h
+ Materials/Lambertian.cc
+ Materials/OrenNayar.h
+ Materials/OrenNayar.cc
+ Materials/LitMaterial.h
+ Materials/LitMaterial.cc
+ Materials/MaterialTable.h
+ Materials/MaterialTable.cc
+ Materials/MetalMaterial.h
+ Materials/MetalMaterial.cc
+ Materials/NDotL.h
+ Materials/NDotL.cc
+ Materials/NullMaterial.h
+ Materials/OpaqueShadower.h
+ Materials/OpaqueShadower.cc
+ Materials/Phong.h
+ Materials/Phong.cc
+ Materials/ThinDielectric.h
+ Materials/ThinDielectric.cc
+ Materials/Transparent.cc
+ Materials/Transparent.h
+ )
+ENDIF(FOUND_TEEM)
+
Added: trunk/Model/Materials/Volume.cc
==============================================================================
--- (empty file)
+++ trunk/Model/Materials/Volume.cc Tue Oct 30 02:53:28 2007
@@ -0,0 +1,1706 @@
+#include <Interface/Scene.h>
+#include <Interface/Renderer.h>
+#include <Interface/RayPacket.h>
+#include <teem/nrrd.h>
+#include <cassert>
+
+#include <SCIRun/Core/Thread/Time.h>
+#include <SCIRun/Core/Thread/Thread.h>
+#include <Core/Geometry/BBox.h>
+#include <Model/Materials/Volume.h>
+using namespace std;
+
+namespace Manta
+{
+
+
+
+ bool lessThan(const ColorSlice& left, const ColorSlice& right)
+ {
+ return left.value <= right.value;
+ }
+
+
+
+
+ RGBAColorMap::RGBAColorMap(unsigned int type, int numSlices)
+ :oldTInc(0.01), _colorScaled(false)
+ {
+ fillColor(type);
+ computeHash();
+ }
+
+ RGBAColorMap::~RGBAColorMap()
+ {}
+
+ void RGBAColorMap::fillColor(unsigned int type)
+ {
+// Add colors base on color map type
+ switch (type) {
+ case InvRainbow:
+ {
+ float p = 1.0f/11.0f;
+ _slices.push_back(ColorSlice(p*0,
RGBAColor(Color(RGB(0, 0, 1)), 1)));
+ _slices.push_back(ColorSlice(p*1,
RGBAColor(Color(RGB(0, 0.40000001, 1)), 1)));
+ _slices.push_back(ColorSlice(p*2,
RGBAColor(Color(RGB(0, 0.80000001, 1)), 1)));
+ _slices.push_back(ColorSlice(p*3,
RGBAColor(Color(RGB(0, 1, 0.80000001)), 1)));
+ _slices.push_back(ColorSlice(p*4,
RGBAColor(Color(RGB(0, 1, 0.40000001)), 1)));
+ _slices.push_back(ColorSlice(p*5,
RGBAColor(Color(RGB(0, 1, 0)), 1)));
+ _slices.push_back(ColorSlice(p*6,
RGBAColor(Color(RGB(0.40000001, 1, 0)), 1)));
+ _slices.push_back(ColorSlice(p*7,
RGBAColor(Color(RGB(0.80000001, 1, 0)), 1)));
+ _slices.push_back(ColorSlice(p*8,
RGBAColor(Color(RGB(1, 0.91764706, 0)), 1)));
+ _slices.push_back(ColorSlice(p*9,
RGBAColor(Color(RGB(1, 0.80000001, 0)), 1)));
+ _slices.push_back(ColorSlice(p*10,
RGBAColor(Color(RGB(1, 0.40000001, 0)), 1)));
+ _slices.push_back(ColorSlice(p*11,
RGBAColor(Color(RGB(1, 0, 0)), 1)));
+ break;
+ }
+ case Rainbow:
+ {
+ float p = 1.0f/11.0f;
+ _slices.push_back(ColorSlice(p*0,
RGBAColor(Color(RGB(1, 0, 0)), 1)));
+ _slices.push_back(ColorSlice(p*1,
RGBAColor(Color(RGB(1, 0.40000001, 0)), 1)));
+ _slices.push_back(ColorSlice(p*2,
RGBAColor(Color(RGB(0, 0.80000001, 0)), 1)));
+ _slices.push_back(ColorSlice(p*3,
RGBAColor(Color(RGB(1, 0.91764706, 0)), 1)));
+ _slices.push_back(ColorSlice(p*4,
RGBAColor(Color(RGB(0.80000001, 1, 0)), 1)));
+ _slices.push_back(ColorSlice(p*5,
RGBAColor(Color(RGB(0.40000001, 1, 0)), 1)));
+ _slices.push_back(ColorSlice(p*6,
RGBAColor(Color(RGB(0, 1, 0)), 1)));
+ _slices.push_back(ColorSlice(p*7,
RGBAColor(Color(RGB(0, 1, 0.40000001)), 1)));
+ _slices.push_back(ColorSlice(p*8,
RGBAColor(Color(RGB(0, 1, 0.80000001)), 1)));
+ _slices.push_back(ColorSlice(p*9,
RGBAColor(Color(RGB(0, 0.80000001, 1)), 1)));
+ _slices.push_back(ColorSlice(p*10,
RGBAColor(Color(RGB(0, 0.40000001, 1)), 1)));
+ _slices.push_back(ColorSlice(p*11,
RGBAColor(Color(RGB(1, 0, 0)), 1)));
+ break;
+ }
+ case InvBlackBody:
+ {
+ float p = 1.0f/12.0f;
+ _slices.push_back(ColorSlice(p*0,
RGBAColor(Color(RGB(1, 1, 1)), 1)));
+ _slices.push_back(ColorSlice(p*1,
RGBAColor(Color(RGB(1, 1, 0.70588237)), 1)));
+ _slices.push_back(ColorSlice(p*2,
RGBAColor(Color(RGB(1, 0.96862745, 0.47058824)), 1)));
+ _slices.push_back(ColorSlice(p*3,
RGBAColor(Color(RGB(1, 0.89411765, 0.3137255)), 1)));
+ _slices.push_back(ColorSlice(p*4,
RGBAColor(Color(RGB(1, 0.80000001, 0.21568628)), 1)));
+ _slices.push_back(ColorSlice(p*5,
RGBAColor(Color(RGB(1, 0.63921571, 0.078431375)), 1)));
+ _slices.push_back(ColorSlice(p*6,
RGBAColor(Color(RGB(1, 0.47058824, 0)), 1)));
+ _slices.push_back(ColorSlice(p*7,
RGBAColor(Color(RGB(0.90196079, 0.27843139, 0)), 1)));
+ _slices.push_back(ColorSlice(p*8,
RGBAColor(Color(RGB(0.78431374, 0.16078432, 0)), 1)));
+ _slices.push_back(ColorSlice(p*9,
RGBAColor(Color(RGB(0.60000002, 0.070588239, 0)), 1)));
+ _slices.push_back(ColorSlice(p*10,
RGBAColor(Color(RGB(0.40000001, 0.0078431377, 0)), 1)));
+ _slices.push_back(ColorSlice(p*11,
RGBAColor(Color(RGB(0.20392157, 0, 0)), 1)));
+ _slices.push_back(ColorSlice(p*12,
RGBAColor(Color(RGB(0, 0, 0)), 1)));
+ break;
+ }
+
+ case BlackBody:
+ {
+ float p = 1.0f/12.0f;
+ _slices.push_back(ColorSlice(p*0,
RGBAColor(Color(RGB(0, 0, 0)), 1)));
+ _slices.push_back(ColorSlice(p*1,
RGBAColor(Color(RGB(0.20392157, 0, 0)), 1)));
+ _slices.push_back(ColorSlice(p*2,
RGBAColor(Color(RGB(0.40000001, 0.0078431377, 0)), 1)));
+ _slices.push_back(ColorSlice(p*3,
RGBAColor(Color(RGB(0.60000002, 0.070588239, 0)), 1)));
+ _slices.push_back(ColorSlice(p*4,
RGBAColor(Color(RGB(0.78431374, 0.16078432, 0)), 1)));
+ _slices.push_back(ColorSlice(p*5,
RGBAColor(Color(RGB(0.90196079, 0.27843139, 0)), 1)));
+ _slices.push_back(ColorSlice(p*6,
RGBAColor(Color(RGB(1, 0.47058824, 0)), 1)));
+ _slices.push_back(ColorSlice(p*7,
RGBAColor(Color(RGB(1, 0.63921571, 0.078431375)), 1)));
+ _slices.push_back(ColorSlice(p*8,
RGBAColor(Color(RGB(1, 0.80000001, 0.21568628)), 1)));
+ _slices.push_back(ColorSlice(p*9,
RGBAColor(Color(RGB(1, 0.89411765, 0.3137255)), 1)));
+ _slices.push_back(ColorSlice(p*10,
RGBAColor(Color(RGB(1, 0.96862745, 0.47058824)), 1)));
+ _slices.push_back(ColorSlice(p*11,
RGBAColor(Color(RGB(1, 1, 0.70588237)), 1)));
+ _slices.push_back(ColorSlice(p*11,
RGBAColor(Color(RGB(1, 1, 1)), 1)));
+ break;
+ }
+ case GreyScale:
+ _slices.push_back(ColorSlice(0,
RGBAColor(Color(RGB(0,0,0)), 1)));
+
_slices.push_back(ColorSlice(1,RGBAColor(Color(RGB(1,1,1)), 1)));
+ break;
+ case InvGreyScale:
+ _slices.push_back(ColorSlice(0,
RGBAColor(Color(RGB(1,1,1)), 1)));
+
_slices.push_back(ColorSlice(1,RGBAColor(Color(RGB(0,0,0)), 1)));
+ break;
+ case InvRainbowIso:
+ {
+ float p = 1.0f/11.0f;
+ _slices.push_back(ColorSlice(p*0,
RGBAColor(Color(RGB(0.528, 0.528, 1.0)), 1)));
+ _slices.push_back(ColorSlice(p*1,
RGBAColor(Color(RGB(0.304, 0.5824, 1.0)), 1)));
+ _slices.push_back(ColorSlice(p*2,
RGBAColor(Color(RGB(0.0, 0.6656, 0.832)), 1)));
+ _slices.push_back(ColorSlice(p*3,
RGBAColor(Color(RGB(0.0, 0.712, 0.5696)), 1)));
+ _slices.push_back(ColorSlice(p*4,
RGBAColor(Color(RGB(0.0, 0.744, 0.2976)), 1)));
+ _slices.push_back(ColorSlice(p*5,
RGBAColor(Color(RGB(0.0, 0.76, 0.0)), 1)));
+ _slices.push_back(ColorSlice(p*6,
RGBAColor(Color(RGB(0.304, 0.76, 0.0)), 1)));
+ _slices.push_back(ColorSlice(p*7,
RGBAColor(Color(RGB(0.5504, 0.688, 0.0)), 1)));
+ _slices.push_back(ColorSlice(p*8,
RGBAColor(Color(RGB(0.68, 0.624, 0.0)), 1)));
+ _slices.push_back(ColorSlice(p*9,
RGBAColor(Color(RGB(0.752, 0.6016, 0.0)), 1)));
+ _slices.push_back(ColorSlice(p*10,
RGBAColor(Color(RGB(1.0, 0.5008, 0.168)), 1)));
+ _slices.push_back(ColorSlice(p*11,
RGBAColor(Color(RGB(1.0, 0.424, 0.424)), 1)));
+ break;
+ }
+ default:
+ cerr<<"RGBAColorMap::fillColor(type="<<type
+ <<"): Invalid type, using
gray scale\n";
+ fillColor(GreyScale);
+ }
+ }
+
+ RGBAColorMap::RGBAColorMap(vector<Color>& colors, float* positions,
float* opacities, int numSlices)
+ {
+ _colorScaled = false;
+ _numSlices = numSlices;
+ _min =0;
+ _max = 1;
+ _invRange = 0;
+ oldTInc = 0.01f;
+ for(size_t i = 0; i < colors.size(); i++)
+ {
+ if (opacities)
+ _slices.push_back(ColorSlice(positions[i],
RGBAColor(colors[i], opacities[i]) ));
+ else
+ _slices.push_back(ColorSlice(positions[i],
RGBAColor(colors[i], 1.0f) ));
+ }
+ std::sort(_slices.begin(), _slices.end(), lessThan);
+ computeHash();
+ }
+
+ RGBAColorMap::RGBAColorMap(vector<ColorSlice>& colors, int numSlices)
+ {
+ _colorScaled = false;
+ _numSlices = numSlices;
+ _min =0;
+ _max = 1;
+ _invRange = 0;
+ oldTInc = 0.01f;
+ //for(int i = 0; i < colors.size(); i++)
+ //{
+ _slices = colors;
+ //}
+ std::sort(_slices.begin(), _slices.end(), lessThan);
+ computeHash();
+ }
+
+ RGBAColor RGBAColorMap::GetColor(float v) // get the interpolated
color at value v (0-1)
+ {
+ if (_slices.size() == 0)
+ return RGBAColor(Color(RGB(0,0,0)), 1.0f);
+ int start = 0;
+ int end = _slices.size()-1;
+
+ if (v < _slices[0].value)
+ return _slices[0].color;
+ if (v > _slices[end].value)
+ return _slices[end].color;
+ //binary search over values
+ int middle;
+ while(end > start+1)
+ {
+ middle = (start+end)/2;
+ if (v < _slices[middle].value)
+ end = middle;
+ else if (v > _slices[middle].value)
+ start = middle;
+ else
+ return _slices[middle].color;
+ }
+ float interp = (v - _slices[start].value) /
(_slices[end].value - _slices[start].value);
+ return RGBAColor(
Color(_slices[start].color.color*(1.0f-interp) +
_slices[end].color.color*interp),
+
_slices[start].color.a*(1.0f-interp) + _slices[end].color.a*interp);
+
+ }
+
+ RGBAColor RGBAColorMap::GetColorAbs(float v)
+ {
+ return GetColor((v - _min)/_invRange);
+ }
+
+ RGBAColor RGBAColorMap::GetPreInterpolatedColor(float v1, float v2,
float d)
+ {
+ /*int i1 = v1*_numSlices;
+ int i2 = v2*_numSlices;
+ int i3 = d*_numSlices;
+ if (i1 >= _numSlices)
+ i1 = _numSlices-1;
+ if (i2 >= _numSlices)
+ i2 = _numSlices-1;
+ if (i3 >= _numSlices)
+ i3 = _numSlices-1;
+ return RGBAColor(_colors[i1][i2][i3],
_opacities[i1][i2][i3]);*/
+ return RGBAColor(0,0,0,1);
+ }
+
+ void RGBAColorMap::SetColors(vector<ColorSlice>& colors)
+ {
+ _slices = colors;
+ std::sort(_slices.begin(), _slices.end(), lessThan);
+ if (_colorScaled)
+ {
+ oldTInc = 0.01;
+ scaleAlphas(tInc);
+ }
+ }
+
+ void RGBAColorMap::SetMinMax(float min, float max)
+ {
+ _min = min;
+ _max = max;
+ _invRange = 1.0/(max - min);
+ }
+
+ void RGBAColorMap::BuildSlices()
+ {
+// Array1<Color>* colors = new Array1<Color>(slices.size());
+// sTColors.set_results_ptr(colors);
+// Array1<float>* opacities = new Array1<float>(slices.size());
+// sTOpacities.set_results_ptr(opacities);
+ //
+// for (int i = 0; i < slices.size(); i++)
+// {
+// sTColors[i] = slices.color.colorl
+ //
+// }
+ }
+
+ void RGBAColorMap::PreInterpolate()
+ {
+ /* cout << "PreInterpolating values... ";
+ _colors = new Color**[_numSlices];
+ _opacities = new float**[_numSlices];
+ for(int i = 0; i < _numSlices; i++)
+ {
+ _colors[i] = new Color*[_numSlices];
+ _opacities[i] = new float*[_numSlices];
+ for(int j = 0; j < _numSlices; j++)
+ {
+ _colors[i][j] = new Color[_numSlices];
+ _opacities[i][j] = new float[_numSlices];
+ }
+ }
+ float step = 1.0/_numSlices;
+ for (int sf = 0; sf < _numSlices; sf++)
+ {
+ for(int sb = 0; sb < _numSlices; sb++)
+ {
+ for(int d = 0; d < _numSlices; d++)
+ {
+ float sff =
float(sf)/float(_numSlices);
+ float sbf =
float(sb)/float(_numSlices);
+ float df = float(d)/float(_numSlices);
+ float Tsf = 0.0f, Tsb = 0.0f;
+ float sumOpacity = 0.0f;
+ Color sumColor(RGB(0,0,0));
+ int minS, maxS;
+ if (sb <= sf)
+ {
+ minS = sb;
+ maxS = sf;
+ }
+ else
+ {
+ minS = sf;
+ maxS = sb;
+ }
+ RGBAColor color(Color(), 1.0);
+ for (int i = minS; i <= maxS; i++)
+ {
+ color =
GetColor(float(i)/float(_numSlices));
+ sumOpacity += color.a*step;
+ sumColor += color.color*step;
+ }
+ if (sb == sf)
+ {
+ _colors[sb][sf][d] =
color.color;
+ _opacities[sb][sf][d] =
color.a;
+ }
+ else
+ {
+ _colors[sb][sf][d] =
sumColor*df/fabs(sbf-sff);
+ _opacities[sb][sf][d] = 1.0f
- pow(2.718281828, - df/fabs(sbf - sff)*sumOpacity);
+ }
+ }
+ }
+ }
+ cout << "done.\n"; */
+ }
+
+
+ void RGBAColorMap::computeHash()
+ {
+
+ unsigned long long new_course_hash = 0;
+ int nindex = 63; // 64 - 1
+
+ for (int a_index = 0; a_index < int(_slices.size()-1);
a_index++) {
+ // This code looks for segments where either the start or the end
+ // is non zeoro. When this happens indices are produces which
+ // round down on the start and round up on the end. This can
+ // cause some overlap at adjacent non zero segments, but this is
+ // OK as we are only turning bits on.
+ float val = _slices[a_index].color.a;
+ float next_val = _slices[a_index+1].color.a;
+ if (val != 0 || next_val != 0) {
+ int start, end;
+ start = (int)(_slices[a_index].value *
nindex);
+ end = (int)ceilf(_slices[a_index+1].value *
nindex);
+ for (int i = start; i <= end; i++)
+ // Turn on the bits.
+ new_course_hash |= 1ULL << i;
+ }
+ }
+
+ course_hash = new_course_hash;
+ cout << "course hash: ";
+ size_t num_bits = sizeof(course_hash)*8;
+ for(size_t i = 0; i < num_bits; i++)
+ if (course_hash & ( 1ULL << i ))
+ cout << 1;
+ else
+ cout << 0;
+ cout << "\n";
+
+ }
+
+
+//scale alpha values by t value
+// Preconditions:
+// assuming that the data in alpha_transform and alpha_stripes is already
+// allocated.
+// alpha_list.size() >= 2
+// alpha_list[0].x == 0
+// alpha_list[alpha_list.size()-1].x == 1
+
+ void RGBAColorMap::scaleAlphas(float t) {
+ _colorScaled = true;
+ cout << "t, oldt: " << t << " " << oldTInc << endl;
+ // the ratio of values as explained in rescale_alphas
+ float d2_div_d1 = t/oldTInc;
+
+ // loop over the alpha values and fill in the array
+ for (int a_index = 0; a_index < int(_slices.size());
a_index++) {
+ // the start has already been computed, so you need to figure out
+ // the end.
+ /*end = (int)(_slices[a_index+1].value * nindex);
+ float val = _slices[a_index].color.a;
+ float val_inc = (_slices[a_index+1].val - val) / (end
- start);
+ for (int i = start; i <= end; i++) {
+ // cout << "val = "<<val<<", ";
+ alpha_stripes[i] = val;
+ // apply the alpha scaling.
+ //
+ // Here we have to be careful with powf. It doesn't like values
+ // that equal 0. I tried just comparing val to 1, but that
+ // didn't work I'm assuming because of floating point precision
+ // problems. I tried comparing 1-val to zero, but that failed
+ // for the same reasons. I finally tried comparing to some
+ // epsilon. That is what you see below.
+ float one_minus_val = 1.0f - val;
+ if (one_minus_val >= 1e-6f)
+ alpha_transform[i] = 1 - powf(one_minus_val,
d2_div_d1);
+ else
+ alpha_transform[i] = 1;
+ // cout <<"alpha_transform[i="<<i<<"] =
"<<alpha_transform[i]<<"\n";
+ val += val_inc;
+ }
+ start = end;*/
+ float val = _slices[a_index].color.a;
+ float one_minus_val = 1.0f - val;
+ if (one_minus_val >= 1e-6f)
+ _slices[a_index].color.a = 1 -
powf(one_minus_val, d2_div_d1);
+ else
+ _slices[a_index].color.a = 1;
+ }
+ oldTInc = tInc;
+ tInc = t;
+ }
+
+
+
+float Volume::resolutionFactor = 5.0f;
+
+void computeMinMax(double& min, double& max, GridArray3<double>& grid)
+{
+ unsigned int i, total = grid.getNx() * grid.getNy() * grid.getNz();
+
+ min = max = grid[0];
+ for(i=1; i<total; i++)
+ {
+ if (min > grid[i]) min = grid[i];
+ else if (max < grid[i]) max = grid[i];
+ }
+}
+
+//template<class T>
+Volume::Volume(string dataNrrd, RGBAColorMap* colorMap, const Vector&
minBound, const Vector& maxBound, double cellStepSize, double metersPerUnit,
int depth, float forceDataMin, float forceDataMax)
+ :_colorMap(colorMap)
+ {
+ _data = new GridArray3<double>();
+ loadNrrd(dataNrrd, _data);
+ //_emittedLight = new GridArray3<double>();
+ //loadNrrd(emittedLightNrrdFile, _emittedLight);
+ _minBound = minBound;
+ _maxBound = maxBound;
+
+ Vector diag = maxBound - minBound;
+
+ _cellSize = diag*Vector( 1.0 / (double)(_data->getNx()-1),
+ 1.0 / (double)(_data->getNy()-1),
+ 1.0 / (double)(_data->getNz()-1));
+
+ //_stepSize = cellStepSize*_cellSize.length()/sqrt(3.0);
+ _stepSize = cellStepSize;
+
+ _metersPerUnit = metersPerUnit;
+ _maxDistance = (_maxBound - _minBound).length();
+
+
+ double min = 0,max = 0;
+ //_data->getMinMax(&min, &max);
+ min = FLT_MAX;
+ max = -FLT_MAX;
+ //computeMinMax(min, max, *_data);
+ unsigned int size = _data->getNx()*_data->getNy()*_data->getNz();
+ for (int i = 0; i < (int)size; i++)
+ {
+ if ((*_data)[i] < min) min = (*_data)[i];
+ if ((*_data)[i] > max) max = (*_data)[i];
+ }
+ if (forceDataMin != -FLT_MAX || forceDataMax != -FLT_MAX)
+ {
+ min = forceDataMin;
+ max = forceDataMax;
+ }
+ _dataMin = min;
+ _dataMax = max;
+ _colorScalar = 1.0f/(_dataMax - _dataMin);
+ cout << "datamin/max: " << min << " " << max << endl;
+
+ _depth = depth;
+ if (_depth <= 0)
+ _depth=1;
+ _datadiag = diag;
+ _nx = _data->getNx();
+ _ny = _data->getNy();
+ _nz = _data->getNz();
+ _sdiag = _datadiag/Vector(_nx-1,_ny-1,_nz-1);
+ diag = (_maxBound - _minBound);
+ inv_diag = Vector(1.0f/diag.x(), 1.0f/diag.y(), 1.0f/diag.z());
+
+ // Compute all the grid stuff
+ _xsize=new int[depth];
+ _ysize=new int[depth];
+ _zsize=new int[depth];
+ int tx=_nx-1;
+ int ty=_ny-1;
+ int tz=_nz-1;
+ for(int i=depth-1;i>=0;i--){
+ int nx=(int)(pow(tx, 1./(i+1))+.9);
+ tx=(tx+nx-1)/nx;
+ _xsize[depth-i-1]=nx;
+ int ny=(int)(pow(ty, 1./(i+1))+.9);
+ ty=(ty+ny-1)/ny;
+ _ysize[depth-i-1]=ny;
+ int nz=(int)(pow(tz, 1./(i+1))+.9);
+ tz=(tz+nz-1)/nz;
+ _zsize[depth-i-1]=nz;
+ }
+ _ixsize=new double[depth];
+ _iysize=new double[depth];
+ _izsize=new double[depth];
+ cerr << "Calculating depths...\n";
+ for(int i=0;i<depth;i++){
+ cerr << "xsize=" << _xsize[i] << ", ysize=" << _ysize[i] <<
", zsize=" << _zsize[i] << '\n';
+ _ixsize[i]=1./_xsize[i];
+ _iysize[i]=1./_ysize[i];
+ _izsize[i]=1./_zsize[i];
+ }
+ cerr << "X: ";
+ tx=1;
+ for(int i=depth-1;i>=0;i--){
+ cerr << _xsize[i] << ' ';
+ tx*=_xsize[i];
+ }
+ cerr << "(" << tx << ")\n";
+ if(tx<_nx-1){
+ cerr << "TX TOO SMALL!\n";
+ exit(1);
+ }
+ cerr << "Y: ";
+ ty=1;
+ for(int i=depth-1;i>=0;i--){
+ cerr << _ysize[i] << ' ';
+ ty*=_ysize[i];
+ }
+ cerr << "(" << ty << ")\n";
+ if(ty<_ny-1){
+ cerr << "TY TOO SMALL!\n";
+ exit(1);
+ }
+ cerr << "Z: ";
+ tz=1;
+ for(int i=depth-1;i>=0;i--){
+ cerr << _zsize[i] << ' ';
+ tz*=_zsize[i];
+ }
+ cerr << "(" << tz << ")\n";
+ if(tz<_nz-1){
+ cerr << "TZ TOO SMALL!\n";
+ exit(1);
+ }
+ _hierdiag=_datadiag*Vector(tx,ty,tz)/Vector(_nx-1,_ny-1,_nz-1);
+ _ihierdiag=Vector(1.,1.,1.)/_hierdiag;
+
+ if(depth==1){
+ macrocells=0;
+ } else {
+ macrocells=new GridArray3<VMCell>[depth+1];
+ int xs=1;
+ int ys=1;
+ int zs=1;
+ for(int d=depth-1;d>=1;d--){
+ xs*=_xsize[d];
+ ys*=_ysize[d];
+ zs*=_zsize[d];
+ macrocells[d].resize(xs, ys, zs);
+ cerr << "Depth " << d << ": " << xs << "x" << ys <<
"x" << zs << '\n';
+ }
+ cerr << "Building hierarchy\n";
+ VMCell top;
+ calc_mcell(depth-1,0,0,0,top);
+ cerr << "done\n";
+ }
+ //cerr << "**************************************************\n";
+ //print(cerr);
+ //cerr << "**************************************************\n";
+ }
+
+ void Volume::calc_mcell(int depth, int startx, int starty, int startz,
VMCell& mcell)
+ {
+ int endx=startx+_xsize[depth];
+ int endy=starty+_ysize[depth];
+ int endz=startz+_zsize[depth];
+ int nx = _nx;
+ int ny = _ny;
+ int nz = _nz;
+ if(endx>nx-1)
+ endx=nx-1;
+ if(endy>ny-1)
+ endy=ny-1;
+ if(endz>nz-1)
+ endz=nz-1;
+ if(startx>=endx || starty>=endy || startz>=endz){
+ /* This cell won't get used... */
+ return;
+ }
+ if(depth==0){
+ // We are at the data level. Loop over each voxel and compute the
+ // mcell for this group of voxels.
+ GridArray3<double>& data = *_data;
+ float data_min = _dataMin;
+ float data_max = _dataMax;
+ for(int ix=startx;ix<endx;ix++){
+ for(int iy=starty;iy<endy;iy++){
+ for(int iz=startz;iz<endz;iz++){
+ float rhos[8];
+ rhos[0]=data(ix, iy, iz);
+ rhos[1]=data(ix, iy, iz+1);
+ rhos[2]=data(ix, iy+1, iz);
+ rhos[3]=data(ix, iy+1, iz+1);
+ rhos[4]=data(ix+1, iy, iz);
+ rhos[5]=data(ix+1, iy, iz+1);
+ rhos[6]=data(ix+1, iy+1, iz);
+ rhos[7]=data(ix+1, iy+1, iz+1);
+ float minr=rhos[0];
+ float maxr=rhos[0];
+ for(int i=1;i<8;i++){
+ if(rhos[i]<minr)
+ minr=rhos[i];
+ if(rhos[i]>maxr)
+ maxr=rhos[i];
+ }
+ // Figure out what bits to turn on running from min to max.
+ mcell.turn_on_bits(minr, maxr,
data_min, data_max);
+ }
+ }
+ }
+ } else {
+ int nxl=_xsize[depth-1];
+ int nyl=_ysize[depth-1];
+ int nzl=_zsize[depth-1];
+ GridArray3<VMCell>& mcells=macrocells[depth];
+ for(int x=startx;x<endx;x++){
+ for(int y=starty;y<endy;y++){
+ for(int z=startz;z<endz;z++){
+ // Compute the mcell for this block and store it in tmp
+ VMCell tmp;
+ calc_mcell(depth-1, x*nxl, y*nyl,
z*nzl, tmp);
+ // Stash it away
+ mcells(x,y,z)=tmp;
+ // Now aggregate all the mcells created for this depth by
+ // doing a bitwise or.
+ mcell |= tmp;
+ }
+ }
+ }
+ }
+ /* if (1) //(depth == (_depth-1))
+ {
+ unsigned long long course_hash = mcell.course_hash;
+ cout << "course hash2: ";
+size_t num_bits = sizeof(course_hash)*8;
+ for(size_t i = 0; i < num_bits; i++)
+ if (course_hash & ( 1ULL << i ))
+ cout << 1;
+ else
+ cout << 0;
+ cout << "\n";
+
+
+ }*/
+
+ }
+
+// template<class T>
+ Volume::~Volume()
+ {
+
+ }
+
+// template<class T>
+ void Volume::preprocess(const PreprocessContext& context)
+ {
+
+ }
+
+// template<class T>
+ void Volume::loadNrrd(string file, GridArray3<double> * grid)
+ {
+ cout << "READING NRRD: " << file << " ";
+ Nrrd *new_nrrd = nrrdNew();
+
////////////////////////////////////////////////////////////////////////////
+ // Attempt to open the nrrd
+ if (nrrdLoad( new_nrrd, file.c_str(), 0 )) {
+ char *reason = biffGetDone( NRRD );
+ std::cout << "WARNING Loading Nrrd Failed: " << reason << std::endl;
+ exit(__LINE__);
+ }
+
+ // Check to make sure the nrrd is the proper dimensions.
+ if (new_nrrd->dim != 3) {
+ std::cout << "WARNING Nrrd must three dimension RGB" << std::endl;
+ exit(__LINE__);
+ }
+
+ // Check that the nrrd is the correct type.
+ // if (new_nrrd->type != nrrdTypeFloat) {
+ // std::cout << "WARNING Nrrd must be type
nrrdTypeDouble" << std::endl;
+ // exit(__LINE__);
+ // }
+
+ // Check that the nrrd could be a RGBA volume.
+ // if (new_nrrd->axis[0].size != 3) {
+ // std::cout << "WARNING Nrrd first axis must be RGB
size 3." << std::endl;
+ // exit(__LINE__);
+ // }
+
+ //Determine scaling for the volume.
+ Vector size( new_nrrd->axis[0].size * new_nrrd->axis[0].spacing,
+ new_nrrd->axis[1].size * new_nrrd->axis[1].spacing,
+ new_nrrd->axis[2].size * new_nrrd->axis[2].spacing );
+ size = Vector(new_nrrd->axis[0].size, new_nrrd->axis[1].size,
new_nrrd->axis[2].size);
+ cout << "size: " << size.x() << " " << size.y() << " " << size.z() <<
endl;
+ // Rescale.
+ float max_dim = max( size.x(), max( size.y(), size.z() ) );
+ cout << "resizing\n";
+ size /= max_dim;
+
+ //_minBound = -size/2.0;
+ //_maxBound = size/2.0;
+cout << "resizing grid\n";
+ grid->resize(new_nrrd->axis[0].size, new_nrrd->axis[1].size,
new_nrrd->axis[2].size);
+cout << "copying nrrd data\n";
+ for(int i=0; i<grid->getNz(); i++)
+ for(int j=0; j<grid->getNy(); j++)
+ for(int k=0; k<grid->getNx(); k++)
+ {
+ //(*grid)(k,j,i) = 0.0;
+ //(*grid)(k,j,i) = (*dataPtr)++;
+ //void* vd = (void *) (&(idata[((i*ny + j)*nx + k) *
formatSize]));
+ if (new_nrrd->type == nrrdTypeFloat)
+ (*grid)(k,j,i) =
((float*)new_nrrd->data)[(int)(((i*grid->getNy() + j)*grid->getNx() + k))];
+ else if (new_nrrd->type == nrrdTypeChar)
+ (*grid)(k,j,i) =
((char*)new_nrrd->data)[(int)(((i*grid->getNy() + j)*grid->getNx() + k))];
+ else if (new_nrrd->type == nrrdTypeDouble)
+ (*grid)(k,j,i) =
((double*)new_nrrd->data)[(int)(((i*grid->getNy() + j)*grid->getNx() + k))];
+ else if (new_nrrd->type == nrrdTypeUChar)
+ (*grid)(k,j,i) = ((unsigned
char*)new_nrrd->data)[(int)(((i*grid->getNy() + j)*grid->getNx() + k))];
+ else if (new_nrrd->type == nrrdTypeShort)
+ (*grid)(k,j,i) =
((short*)new_nrrd->data)[(int)(((i*grid->getNy() + j)*grid->getNx() + k))];
+ else if (new_nrrd->type == nrrdTypeUShort)
+ (*grid)(k,j,i) = ((unsigned
short*)new_nrrd->data)[(int)(((i*grid->getNy() + j)*grid->getNx() + k))];
+ else if (new_nrrd->type == nrrdTypeInt)
+ (*grid)(k,j,i) = ((int*)new_nrrd->data)[(int)(((i*grid->getNy()
+ j)*grid->getNx() + k))];
+ else if (new_nrrd->type == nrrdTypeUInt)
+ (*grid)(k,j,i) = ((unsigned
int*)new_nrrd->data)[(int)(((i*grid->getNy() + j)*grid->getNx() + k))];
+ // if (new_nrrd->type == nrrdTypeLong)
+ // (*grid)(k,j,i) =
((long*)new_nrrd->data)[(int)(((i*grid->getNy() + j)*grid->getNx() + k))];
+ // if (new_nrrd->type == nrrdTypeULong)
+ // (*grid)(k,j,i) = ((unsigned
long*)new_nrrd->data)[(int)(((i*grid->getNy() + j)*grid->getNx() + k))];
+ }
+ cout << "done\n";
+ }
+
+
+#define RAY_TERMINATION_THRESHOLD 0.9
+
+ void Volume::isect(int depth, double t_sample,
+ double dtdx, double dtdy, double dtdz,
+ double next_x, double next_y, double next_z,
+ int ix, int iy, int iz,
+ int startx, int starty, int startz,
+ const Vector& cellcorner, const Vector&
celldir,
+ HVIsectContext &isctx) const
+ {
+
+ int cx=_xsize[depth];
+ int cy=_ysize[depth];
+ int cz=_zsize[depth];
+
+ GridArray3<double>& data = *_data;
+
+ if(depth==0){
+ int nx = _nx;
+ int ny = _ny;
+ int nz = _nz;
+ for(;;){
+ int gx=startx+ix;
+ int gy=starty+iy;
+ int gz=startz+iz;
+
+ // t is our t_sample
+
+ // If we have valid samples
+ if(gx<nx-1 && gy<ny-1 && gz<nz-1){
+
+ float rhos[8];
+ rhos[0]=data(gx, gy, gz);
+ rhos[1]=data(gx, gy, gz+1);
+ rhos[2]=data(gx, gy+1, gz);
+ rhos[3]=data(gx, gy+1, gz+1);
+ rhos[4]=data(gx+1, gy, gz);
+ rhos[5]=data(gx+1, gy, gz+1);
+ rhos[6]=data(gx+1, gy+1, gz);
+ rhos[7]=data(gx+1, gy+1, gz+1);
+
+
////////////////////////////////////////////////////////////
+ // get the weights
+
+ Vector weights =
cellcorner+celldir*t_sample;
+ double x_weight_high = weights.x()-ix;
+ double y_weight_high = weights.y()-iy;
+ double z_weight_high = weights.z()-iz;
+
+
+ double lz1, lz2, lz3, lz4, ly1, ly2, value;
+ lz1 = rhos[0] * (1 - z_weight_high) +
rhos[1] * z_weight_high;
+ lz2 = rhos[2] * (1 - z_weight_high) +
rhos[3] * z_weight_high;
+ lz3 = rhos[4] * (1 - z_weight_high) +
rhos[5] * z_weight_high;
+ lz4 = rhos[6] * (1 - z_weight_high) +
rhos[7] * z_weight_high;
+
+ ly1 = lz1 * (1 - y_weight_high) + lz2 *
y_weight_high;
+ ly2 = lz3 * (1 - y_weight_high) + lz4 *
y_weight_high;
+
+ value = ly1 * (1 - x_weight_high) + ly2 *
x_weight_high;
+ value = (value - _dataMin)*_colorScalar;
+
+ //float alpha_factor =
this->dpy->lookup_alpha(value) * (1-isctx.alpha);
+ RGBAColor color =
_colorMap->GetColor(value);
+ float alpha_factor =
color.a*(1-isctx.alpha);
+ if (alpha_factor > 0.001) {
+ // the point is contributing, so compute the color
+
+ /*Light*
light=isctx.cx->scene->light(0);
+ if (light->isOn()) {
+
+ // compute the gradient
+ Vector gradient;
+ double dx = ly2 - ly1;
+
+ double dy, dy1, dy2;
+ dy1 = lz2 - lz1;
+ dy2 = lz4 - lz3;
+ dy = dy1 * (1 -
x_weight_high) + dy2 * x_weight_high;
+
+ double dz, dz1, dz2, dz3,
dz4, dzly1, dzly2;
+ dz1 = rhos[1] - rhos[0];
+ dz2 = rhos[3] - rhos[2];
+ dz3 = rhos[5] - rhos[4];
+ dz4 = rhos[7] - rhos[6];
+ dzly1 = dz1 * (1 -
y_weight_high) + dz2 * y_weight_high;
+ dzly2 = dz3 * (1 -
y_weight_high) + dz4 * y_weight_high;
+ dz = dzly1 * (1 -
x_weight_high) + dzly2 * x_weight_high;
+ double length2 =
dx*dx+dy*dy+dz*dz;
+ if (length2){
+ // this lets the compiler use a special 1/sqrt() operation
+ double ilength2 =
1/sqrt(length2);
+ gradient =
Vector(dx*ilength2, dy*ilength2,dz*ilength2);
+ } else {
+ gradient =
Vector(0,0,0);
+ }
+
+ Vector light_dir;
+ Point current_p =
isctx.ray.origin() + isctx.ray.direction()*t_sample - this->min.vector();
+ light_dir =
light->get_pos()-current_p;
+
+ Color temp =
color(gradient, isctx.ray.direction(),
+
light_dir.normal(),
+
*(this->dpy->lookup_color(value)),
+
light->get_color());
+ isctx.total += temp *
alpha_factor;
+ } else {*/
+ isctx.total +=
color.color*(alpha_factor);
+ //}
+ isctx.alpha += alpha_factor;
+ }
+
+ }
+
+ // Update the new position
+ t_sample += isctx.t_inc;
+
+ // If we have overstepped the limit of the ray
+ if(t_sample >= isctx.t_max)
+ break;
+
+ // Check to see if we leave the level or not
+ bool break_forloop = false;
+ while (t_sample > next_x) {
+ // Step in x...
+ next_x+=dtdx;
+ ix+=isctx.dix_dx;
+ if(ix<0 || ix>=cx) {
+ break_forloop = true;
+ break;
+ }
+ }
+ while (t_sample > next_y) {
+ next_y+=dtdy;
+ iy+=isctx.diy_dy;
+ if(iy<0 || iy>=cy) {
+ break_forloop = true;
+ break;
+ }
+ }
+ while (t_sample > next_z) {
+ next_z+=dtdz;
+ iz+=isctx.diz_dz;
+ if(iz<0 || iz>=cz) {
+ break_forloop = true;
+ break;
+ }
+ }
+
+ if (break_forloop)
+ break;
+
+ // This does early ray termination when we don't have anymore
+ // color to collect.
+ if (isctx.alpha >= RAY_TERMINATION_THRESHOLD)
+ break;
+ }
+ } else {
+ GridArray3<VMCell>& mcells=macrocells[depth];
+ for(;;){
+ int gx=startx+ix;
+ int gy=starty+iy;
+ int gz=startz+iz;
+ bool hit = mcells(gx,gy,gz) &
isctx.transfunct;
+ if(hit){
+ // Do this cell...
+ int new_cx=_xsize[depth-1];
+ int new_cy=_ysize[depth-1];
+ int new_cz=_zsize[depth-1];
+ int
new_ix=(int)((cellcorner.x()+t_sample*celldir.x()-ix)*new_cx);
+ int
new_iy=(int)((cellcorner.y()+t_sample*celldir.y()-iy)*new_cy);
+ int
new_iz=(int)((cellcorner.z()+t_sample*celldir.z()-iz)*new_cz);
+ if(new_ix<0)
+ new_ix=0;
+ else if(new_ix>=new_cx)
+ new_ix=new_cx-1;
+ if(new_iy<0)
+ new_iy=0;
+ else if(new_iy>=new_cy)
+ new_iy=new_cy-1;
+ if(new_iz<0)
+ new_iz=0;
+ else if(new_iz>=new_cz)
+ new_iz=new_cz-1;
+
+ double new_dtdx=dtdx*_ixsize[depth-1];
+ double new_dtdy=dtdy*_iysize[depth-1];
+ double new_dtdz=dtdz*_izsize[depth-1];
+ const Vector dir(isctx.ray.direction());
+ double new_next_x;
+ if(dir.x() > 0){
+
new_next_x=next_x-dtdx+new_dtdx*(new_ix+1);
+ } else {
+ new_next_x=next_x-new_ix*new_dtdx;
+ }
+ double new_next_y;
+ if(dir.y() > 0){
+
new_next_y=next_y-dtdy+new_dtdy*(new_iy+1);
+ } else {
+ new_next_y=next_y-new_iy*new_dtdy;
+ }
+ double new_next_z;
+ if(dir.z() > 0){
+
new_next_z=next_z-dtdz+new_dtdz*(new_iz+1);
+ } else {
+ new_next_z=next_z-new_iz*new_dtdz;
+ }
+ int new_startx=gx*new_cx;
+ int new_starty=gy*new_cy;
+ int new_startz=gz*new_cz;
+ Vector cellsize(new_cx, new_cy, new_cz);
+ isect(depth-1, t_sample,
+ new_dtdx, new_dtdy,
new_dtdz,
+ new_next_x, new_next_y,
new_next_z,
+ new_ix, new_iy, new_iz,
+ new_startx, new_starty,
new_startz,
+ (cellcorner-Vector(ix, iy,
iz))*cellsize, celldir*cellsize,
+ isctx);
+ }
+
+ // We need to determine where the next sample is. We do this
+ // using the closest crossing in x/y/z. This will be the next
+ // sample point.
+ double closest;
+ if(next_x < next_y && next_x < next_z){
+ // next_x is the closest
+ closest = next_x;
+ } else if(next_y < next_z){
+ closest = next_y;
+ } else {
+ closest = next_z;
+ }
+
+ double step = ceil((closest -
t_sample)*isctx.t_inc_inv);
+ t_sample += isctx.t_inc * step;
+
+ if(t_sample >= isctx.t_max)
+ break;
+
+ // Now that we have the next sample point, we need to determine
+ // which cell it ended up in. There are cases (grazing corners
+ // for example) which will make the next sample not be in the
+ // next cell. Because this case can happen, this code will try
+ // to determine which cell the sample will live.
+ //
+ // Perhaps there is a way to use cellcorner and cellsize to get
+ // ix/y/z. The result would have to be cast to an int, and then
+ // next_x/y/z would have to be updated as needed. (next_x +=
+ // dtdx * (newix - ix). The advantage of something like this
+ // would be the lack of branches, but it does use casts.
+ bool break_forloop = false;
+ while (t_sample > next_x) {
+ // Step in x...
+ next_x+=dtdx;
+ ix+=isctx.dix_dx;
+ if(ix<0 || ix>=cx) {
+ break_forloop = true;
+ break;
+ }
+ }
+ while (t_sample > next_y) {
+ next_y+=dtdy;
+ iy+=isctx.diy_dy;
+ if(iy<0 || iy>=cy) {
+ break_forloop = true;
+ break;
+ }
+ }
+ while (t_sample > next_z) {
+ next_z+=dtdz;
+ iz+=isctx.diz_dz;
+ if(iz<0 || iz>=cz) {
+ break_forloop = true;
+ break;
+ }
+ }
+ if (break_forloop)
+ break;
+
+ if (isctx.alpha >= RAY_TERMINATION_THRESHOLD)
+ break;
+ }
+ }
+ }
+
+//#define CD_SSE
+#define RAY_TERMINATION 0.9
+//template<class T>
+ void Volume::shade(const RenderContext & context, RayPacket& rays) const
+ {
+ bool BITMASK_VERSION = true;
+if (BITMASK_VERSION)
+{
+ rays.normalizeDirections();
+ double t_inc = _stepSize;
+
+ RayPacketData rpData1;
+ RayPacket lRays1(rpData1, RayPacket::SquarePacket, rays.begin(),
rays.end(), rays.getDepth()+1,
+
RayPacket::NormalizedDirections);
+ float tMins[rays.end()];
+ float tMaxs[rays.end()];
+ float alphas[rays.end()];
+ Color totals[rays.end()];
+ for(int rayIndex = rays.begin(); rayIndex < rays.end(); rayIndex++)
+ {
+ lRays1.setRay(rayIndex, Ray(rays.getOrigin(rayIndex)+
rays.getMinT(rayIndex)*rays.getDirection(rayIndex),
rays.getDirection(rayIndex)));
+ tMins[rayIndex] = rays.getMinT(rayIndex);
+ alphas[rayIndex] = 0;
+ }
+ lRays1.resetHits();
+ context.scene->getObject()->intersect(context, lRays1);
+ for(int i = rays.begin(); i < rays.end(); i++)
+ {
+ tMaxs[i] = lRays1.getMinT(i);
+ if (lRays1.wasHit(i))
+ tMaxs[i] = lRays1.getMinT(i)+ tMins[i];
+ else
+ tMaxs[i] = -1;
+ totals[i] = Color(RGB(0,0,0));
+ }
+
+
+ for(int rayIndex = rays.begin(); rayIndex < rays.end(); rayIndex++)
+ {
+ Ray ray = rays.getRay(rayIndex);
+ float t_min = tMins[rayIndex];
+ float t_max = tMaxs[rayIndex];
+
+ const Vector dir(ray.direction());
+ const Vector orig(ray.origin());
+ int dix_dx;
+ int ddx;
+ if(dir.x() > 0){
+ dix_dx=1;
+ ddx=1;
+ } else {
+ dix_dx=-1;
+ ddx=0;
+ }
+ int diy_dy;
+ int ddy;
+ if(dir.y() > 0){
+ diy_dy=1;
+ ddy=1;
+ } else {
+ diy_dy=-1;
+ ddy=0;
+ }
+ int diz_dz;
+ int ddz;
+ if(dir.z() > 0){
+ diz_dz=1;
+ ddz=1;
+ } else {
+ diz_dz=-1;
+ ddz=0;
+ }
+
+ Vector start_p(orig+dir*t_min);
+ Vector s((start_p-_minBound)*_ihierdiag);
+ int cx=_xsize[_depth-1];
+ int cy=_ysize[_depth-1];
+ int cz=_zsize[_depth-1];
+ int ix=(int)(s.x()*cx);
+ int iy=(int)(s.y()*cy);
+ int iz=(int)(s.z()*cz);
+ if(ix>=cx)
+ ix--;
+ if(iy>=cy)
+ iy--;
+ if(iz>=cz)
+ iz--;
+ if(ix<0)
+ ix++;
+ if(iy<0)
+ iy++;
+ if(iz<0)
+ iz++;
+
+ double next_x, next_y, next_z;
+ double dtdx, dtdy, dtdz;
+
+ double icx=_ixsize[_depth-1];
+ double x=_minBound.x()+_hierdiag.x()*double(ix+ddx)*icx;
+ double xinv_dir=1./dir.x();
+ next_x=(x-orig.x())*xinv_dir;
+ dtdx=dix_dx*_hierdiag.x()*icx*xinv_dir;
+
+ double icy=_iysize[_depth-1];
+ double y=_minBound.y()+_hierdiag.y()*double(iy+ddy)*icy;
+ double yinv_dir=1./dir.y();
+ next_y=(y-orig.y())*yinv_dir;
+ dtdy=diy_dy*_hierdiag.y()*icy*yinv_dir;
+
+ double icz=_izsize[_depth-1];
+ double z=_minBound.z()+_hierdiag.z()*double(iz+ddz)*icz;
+ double zinv_dir=1./dir.z();
+ next_z=(z-orig.z())*zinv_dir;
+ dtdz=diz_dz*_hierdiag.z()*icz*zinv_dir;
+
+ Vector cellsize(cx,cy,cz);
+ // cellcorner and celldir can be used to get the location in terms
+ // of the metacell in index space.
+ //
+ // For example if you wanted to get the location at time t (world
+ // space units) in terms of indexspace you would do the following
+ // computation:
+ //
+ // Vector pos = cellcorner + celldir * t + Vector(startx, starty,
startz);
+ //
+ // If you wanted to get how far you are inside a given cell you
+ // could use the following code:
+ //
+ // Vector weights = cellcorner + celldir * t - Vector(ix, iy, iz);
+ Vector cellcorner((orig-_minBound)*_ihierdiag*cellsize);
+ Vector celldir(dir*_ihierdiag*cellsize);
+
+ HVIsectContext isctx;
+ isctx.total = totals[rayIndex];
+ isctx.alpha = alphas[rayIndex];
+ isctx.dix_dx = dix_dx;
+ isctx.diy_dy = diy_dy;
+ isctx.diz_dz = diz_dz;
+ isctx.transfunct.course_hash = _colorMap->course_hash;
+ isctx.t_inc = t_inc;
+ isctx.t_min = t_min;
+ isctx.t_max = t_max;
+ isctx.t_inc_inv = 1/isctx.t_inc;
+ isctx.ray = ray;
+
+ isect(_depth-1, t_min, dtdx, dtdy, dtdz, next_x, next_y,
next_z,
+ ix, iy, iz, 0, 0, 0,
+ cellcorner, celldir,
+ isctx);
+
+ alphas[rayIndex] = isctx.alpha;
+ totals[rayIndex] = isctx.total;
+ }
+
+ for(int i = rays.begin(); i < rays.end(); i++)
+ {
+ if (alphas[i] < RAY_TERMINATION)
+ {
+ Ray ray;
+ ray = rays.getRay(i);
+ if (lRays1.getHitMaterial(i) == this || tMaxs[i] ==
-1)
+ {
+ lRays1.setRay(i, Ray(ray.origin() +
ray.direction()*tMaxs[i], ray.direction()));
+ if (tMaxs[i] == -1)
+ lRays1.setRay(i,
Ray(ray.origin()+ray.direction()*(tMins[i]+0.0001f), ray.direction()));
+ }
+ else
+ lRays1.setRay(i,
Ray(ray.origin()+ray.direction()*(tMins[i]), ray.direction()));
+ }
+ }
+ bool depth = (rays.getDepth() <
context.scene->getRenderParameters().maxDepth);
+ context.renderer->traceRays(context, lRays1);
+ for(int i = rays.begin(); i < rays.end(); i++)
+ {
+ Color bgColor(RGB(0,0,0));
+ if (depth)
+ bgColor = lRays1.getColor(i);
+ totals[i] += bgColor*(1.-alphas[i]);
+ rays.setColor(i, totals[i]);
+ }
+ return;
+}
+
+ rays.normalizeDirections();
+ double t_inc = _stepSize;
+ //t_inc = 0.003;
+ Vector diag = (_maxBound - _minBound);
+ Vector inv_diag (1.0f/diag.x(), 1.0f/diag.y(), 1.0f/diag.z());
+ GridArray3<double>& data = *_data;
+ int nx = data.getNx();
+ int ny = data.getNy();
+ int nz = data.getNz();
+ RayPacketData rpData1;
+ RayPacket lRays1(rpData1, RayPacket::SquarePacket, rays.begin(),
rays.end(), rays.getDepth()+1,
+ RayPacket::NormalizedDirections);
+ float tMins[rays.end()];
+ float tMaxs[rays.end()];
+ float alphas[rays.end()];
+ Color totals[rays.end()];
+ for(int rayIndex = rays.begin(); rayIndex < rays.end(); rayIndex++)
+ {
+ lRays1.setRay(rayIndex, Ray(rays.getOrigin(rayIndex)+
rays.getMinT(rayIndex)*rays.getDirection(rayIndex),
rays.getDirection(rayIndex)));
+ tMins[rayIndex] = rays.getMinT(rayIndex);
+ alphas[rayIndex] = 0;
+ }
+ lRays1.resetHits();
+ context.scene->getObject()->intersect(context, lRays1);
+ for(int i = rays.begin(); i < rays.end(); i++)
+ {
+ tMaxs[i] = lRays1.getMinT(i);
+ if (lRays1.wasHit(i))
+ tMaxs[i] = lRays1.getMinT(i)+ tMins[i];
+ else
+ tMaxs[i] = -1;
+ totals[i] = Color(RGB(0,0,0));
+ }
+
+
+ for(int rayIndex = rays.begin(); rayIndex < rays.end(); rayIndex++)
+ {
+ Ray ray = rays.getRay(rayIndex);
+
+ for(double t = tMins[rayIndex]; t < tMaxs[rayIndex]; t+=
t_inc)
+ {
+ if (alphas[rayIndex] < RAY_TERMINATION)
+ {
+ Vector current_p =
rays.getOrigin(rayIndex)+rays.getDirection(rayIndex)*t-_minBound;
+ ////////////////////////////////////////////////////////////
+ // interpolate the point
+
+ // get the indices and weights for the indicies
+ double norm = current_p.x() * inv_diag.x();
+ double step = norm * (nx - 1);
+ int x_low = SCIRun::Clamp((int)step, 0, nx-2);
+ float x_weight_high = step - x_low;
+
+ norm = current_p.y() * inv_diag.y();
+ step = norm * (ny - 1);
+ int y_low = SCIRun::Clamp((int)step, 0, ny-2);
+ float y_weight_high = step - y_low;
+
+ norm = current_p.z() * inv_diag.z();
+ step = norm * (nz - 1);
+ int z_low = SCIRun::Clamp((int)step, 0, nz-2);
+ float z_weight_high = step - z_low;
+
+
+ double a,b,c,d,e,f,g,h;
+ a = data(x_low, y_low, z_low);
+ b = data(x_low, y_low, z_low+1);
+ c = data(x_low, y_low+1, z_low);
+ d = data(x_low, y_low+1, z_low+1);
+ e = data(x_low+1, y_low, z_low);
+ f = data(x_low+1, y_low, z_low+1);
+ g = data(x_low+1, y_low+1, z_low);
+ h = data(x_low+1, y_low+1, z_low+1);
+
+ float lz1, lz2, lz3, lz4, ly1, ly2, value;
+ lz1 = a * (1 - z_weight_high) + b * z_weight_high;
+ lz2 = c * (1 - z_weight_high) + d * z_weight_high;
+ lz3 = e * (1 - z_weight_high) + f * z_weight_high;
+ lz4 = g * (1 - z_weight_high) + h * z_weight_high;
+
+ ly1 = lz1 * (1 - y_weight_high) + lz2 * y_weight_high;
+ ly2 = lz3 * (1 - y_weight_high) + lz4 * y_weight_high;
+
+ value = ly1 * (1 - x_weight_high) + ly2 *
x_weight_high;
+ //nmormalize
+ value = (value - _dataMin)*_colorScalar;
+ RGBAColor color = _colorMap->GetColor(value);
+ float
alpha_factor = color.a*(1.0-alphas[rayIndex]);
+ if (alpha_factor >0.001)
+ {
+ //TODO: lighting?
+ totals[rayIndex] +=
color.color*(alpha_factor);
+ alphas[rayIndex] += alpha_factor;
+ }
+ else
+ continue;
+ }
+ }
+
+ }
+
+ for(int i = rays.begin(); i < rays.end(); i++)
+ {
+ if (alphas[i] < RAY_TERMINATION)
+ {
+ Ray ray;
+ ray = rays.getRay(i);
+ if (lRays1.getHitMaterial(i) == this || tMaxs[i] ==
-1)
+ {
+ lRays1.setRay(i, Ray(ray.origin() +
ray.direction()*tMaxs[i], ray.direction()));
+ if (tMaxs[i] == -1)
+ lRays1.setRay(i,
Ray(ray.origin()+ray.direction()*(tMins[i]+0.0001f), ray.direction()));
+ }
+ else
+ lRays1.setRay(i,
Ray(ray.origin()+ray.direction()*(tMins[i]), ray.direction()));
+ }
+ }
+ bool depth = (rays.getDepth() <
context.scene->getRenderParameters().maxDepth);
+ context.renderer->traceRays(context, lRays1);
+ for(int i = rays.begin(); i < rays.end(); i++)
+ {
+ Color bgColor(RGB(0,0,0));
+ if (depth)
+ bgColor = lRays1.getColor(i);
+ totals[i] += bgColor*(1.-alphas[i]);
+ rays.setColor(i, totals[i]);
+ }
+ return;
+
+
+ #ifndef CD_SSE
+ for( int rayIndex = rays.begin(); rayIndex < rays.end()-3;
rayIndex+=4)
+ {
+ float minTsF[4];
+ float maxTsF[4];
+ float currTsF[4];
+ float alphas[4] = {0,0,0,0};
+
+ Color accumColors[4];
+ float accumTransmittances[4];
+
+ Vector steps[4];
+ Vector latticePositions[4];
+ RayPacketData rpData;
+ RayPacket lRays(rpData, RayPacket::SquarePacket, 0, 4,
rays.getDepth()+1,
+ RayPacket::NormalizedDirections);
+
+ for(int i = 0; i < 4; i++)
+ {
+ Vector temp(_cellSize);
+ Ray ray = rays.getRay(rayIndex+i);
+ steps[i] = ray.direction()*_stepSize/temp;
+ accumColors[i] = Color(RGB(0,0,0));
+ accumTransmittances[i] = 1.0f;
+ minTsF[i] = rays.getMinT(rayIndex+i);
+ currTsF[i] = 0.0f;
+ Ray lRay;
+ lRay.set(ray.origin() + minTsF[i]*ray.direction(),
ray.direction());
+ lRays.setRay(i, lRay);
+ }
+ lRays.resetHits();
+ context.scene->getObject()->intersect(context, lRays);
+ for(int i = 0; i < 4; i++)
+ {
+ if (lRays.wasHit(i))
+ maxTsF[i] = lRays.getMinT(i);
+ else
+ maxTsF[i] = 0.0;
+ Vector hitPosition = lRays.getOrigin(i);
+ latticePositions[i] = ((hitPosition - _minBound) /
_cellSize);
+ }
+ while (currTsF[0] < maxTsF[0] || currTsF[1] < maxTsF[1] ||
currTsF[2] < maxTsF[2] || currTsF[3] < maxTsF[3])
+ {
+
+ for(int i = 0; i < 4; i++)
+ {
+ float hitLatticef[3] = {latticePositions[i].x(),
latticePositions[i].y(), latticePositions[i].z() };
+ int hitLattice[3] =
{static_cast<int>(floor(hitLatticef[0])),
static_cast<int>(floor(hitLatticef[1])),
static_cast<int>(floor(hitLatticef[2]))};
+ //check bounds
+ if (currTsF[i] <= maxTsF[i] && 0 <= hitLattice[0] &&
hitLattice[0] < (int)_data->getNx() - 2 &&
+ 0 <= hitLattice[1] && hitLattice[1] <
(int)_data->getNy() - 2 &&
+ 0 <= hitLattice[2] && hitLattice[2] <
(int)_data->getNz() - 2 )
+ {
+ float wx, wy, wz, wx1, wy1, wz1, absCoeff =
0.0;
+ Color sourceColor;
+ double tIntegration;
+
+ wx = hitLatticef[0] - (float)hitLattice[0];
+ wy = hitLatticef[1] - (float)hitLattice[1];
+ wz = hitLatticef[2] - (float)hitLattice[2];
+ wx1 = 1.0f-wx;
+ wy1 = 1.0f-wy;
+ wz1 = 1.0f-wz;
+
+ tIntegration = maxTsF[i]-currTsF[i];
+ if (tIntegration > _stepSize)
+ tIntegration = _stepSize;
+ if (_data)
+ absCoeff = wx1 * ( wy1 * ( wz1 *
(*_data)(hitLattice[0] , hitLattice[1] , hitLattice[2] ) +
+ wz *
(*_data)(hitLattice[0] , hitLattice[1] , hitLattice[2]+1) ) +
+ wy * ( wz1 *
(*_data)(hitLattice[0] , hitLattice[1]+1, hitLattice[2] ) +
+ wz *
(*_data)(hitLattice[0] , hitLattice[1]+1, hitLattice[2]+1) ) ) +
+ wx * ( wy1 * ( wz1 *
(*_data)(hitLattice[0]+1, hitLattice[1] , hitLattice[2] ) +
+ wz *
(*_data)(hitLattice[0]+1, hitLattice[1] , hitLattice[2]+1) ) +
+ wy * ( wz1 *
(*_data)(hitLattice[0]+1, hitLattice[1]+1, hitLattice[2] ) +
+ wz *
(*_data)(hitLattice[0]+1, hitLattice[1]+1, hitLattice[2]+1) ) ) ;
+ else
+ absCoeff = 0.0f;
+ float val = (absCoeff -
_dataMin)/(_dataMax-_dataMin);
+ if (val > 1.0)
+ cout << "OVER" << val
<< endl;
+ if (val < 0.0)
+ cout << "UNDER" <<
val << " " << absCoeff << endl;
+
+ RGBAColor color = _colorMap->GetColor(val);
+ sourceColor = color.color;
+ float alpha_factor = color.a*(1.0-alphas[i]);
+ if (alpha_factor > 0.001)
+ accumColors[i] += sourceColor*alpha_factor;
+ alphas[i] += alpha_factor;
+
+
+ }
+ latticePositions[i] += steps[i];
+ currTsF[i] += _stepSize;
+ }
+ }
+ lRays.resetHits();
+ context.renderer->traceRays(context, lRays);
+ for(int i = 0; i < 4; i++)
+ {
+ rays.setColor(rayIndex+i, accumColors[i]);
+ //}
+ }
+ }
+ #else
+ //SSE
+ static sse_t absCoeffNx = set4((int)_data->getNx() - 1);
+ static sse_t absCoeffNy = set4((int)_data->getNy() - 1);
+ static sse_t absCoeffNz = set4((int)_data->getNz() - 1);
+ static sse_t zeros = set4(0);
+ static sse_t ones = set4(1);
+ static Vector _cellSizeInv(1.0f/_cellSize.x(),
1.0f/_cellSize.y(), 1.0f/_cellSize.z());
+ for( int rayIndex = rays.begin(); rayIndex < rays.end()-3;
rayIndex+=4)
+ {
+ Ray rays4[4] = { rays.getRay(rayIndex),
rays.getRay(rayIndex+1), rays.getRay(rayIndex+2), rays.getRay(rayIndex+3) };
+ sse_t dir_x = set44(rays4[0].direction().x(),
rays4[1].direction().x(),
+ rays4[2].direction().x()
,rays4[3].direction().x());
+ sse_t dir_y = set44(rays4[0].direction().y(),
rays4[1].direction().y(),
+ rays4[2].direction().y()
,rays4[3].direction().y());
+ sse_t dir_z = set44(rays4[0].direction().z(),
rays4[1].direction().z(),
+ rays4[2].direction().z()
,rays4[3].direction().z());
+
+ sse_t steps_x = set4(_cellSizeInv.x());
+ sse_t steps_y = set4(_cellSizeInv.y());
+ sse_t steps_z = set4(_cellSizeInv.z());
+ sse_t stepSize4 = set4(_stepSize);
+ steps_x = mul4(steps_x, dir_x); steps_x = mul4(steps_x,
stepSize4);
+ steps_y = mul4(steps_y, dir_y); steps_y = mul4(steps_y,
stepSize4);
+ steps_z = mul4(steps_z, dir_z); steps_z = mul4(steps_z,
stepSize4);
+
+ float minTsF[4] = {0.1, 0.1, 0.1, 0.1};
+ float maxTsF[4] = {0,0,0,0};
+ float currTsF[4] = {minTsF[0], minTsF[1], minTsF[2],
minTsF[3]};
+ float alphas[4] = {0,0,0,0};
+ sse_t currTsF4 = set44(minTsF[0], minTsF[1], minTsF[2],
minTsF[3]);
+
+ Color accumColors[4];
+ float accumTransmittances[4];
+
+ RayPacketData rpData;
+ RayPacket lRays(rpData, RayPacket::SquarePacket, 0, 4,
rays.getDepth()+1, RayPacket::NormalizedDirections);
+
+ for(int i = 0; i < 4; i++)
+ {
+ //cout << steps[i] << endl;
+ accumColors[i] = Color(RGB(0,0,0));
+ accumTransmittances[i] = 1.0;
+ Ray ray = rays.getRay(rayIndex+i);
+ minTsF[i] = rays.getMinT(rayIndex+i);
+ //currTsF[i] = 0.0f;
+ Ray lRay;
+ lRay.set(ray.origin() + minTsF[i]*ray.direction(),
ray.direction());
+ lRays.setRay(i, lRay);
+ }
+ sse_t minTsF4 = set44(minTsF[0], minTsF[1], minTsF[2],
minTsF[3]);
+ lRays.resetHits();
+ context.scene->getObject()->intersect(context, lRays);
+ sse_t origins_x = set44(lRays.getOrigin(0).x(),
lRays.getOrigin(1).x(), lRays.getOrigin(2).x(), lRays.getOrigin(3).x());
+ sse_t origins_y = set44(lRays.getOrigin(0).y(),
lRays.getOrigin(1).y(), lRays.getOrigin(2).y(), lRays.getOrigin(3).y());
+ sse_t origins_z = set44(lRays.getOrigin(0).z(),
lRays.getOrigin(1).z(), lRays.getOrigin(2).z(), lRays.getOrigin(3).z());
+ sse_t minBound_x = set44(_minBound.x(), _minBound.x(),
_minBound.x(), _minBound.x());
+ sse_t minBound_y = set44(_minBound.y(), _minBound.y(),
_minBound.y(), _minBound.y());
+ sse_t minBound_z = set44(_minBound.z(), _minBound.z(),
_minBound.z(), _minBound.z());
+ sse_t cellSizeInv_x = set44(_cellSizeInv.x(),
_cellSizeInv.x(), _cellSizeInv.x(), _cellSizeInv.x());
+ sse_t cellSizeInv_y = set44(_cellSizeInv.y(),
_cellSizeInv.y(), _cellSizeInv.y(), _cellSizeInv.y());
+ sse_t cellSizeInv_z = set44(_cellSizeInv.z(),
_cellSizeInv.z(), _cellSizeInv.z(), _cellSizeInv.z());
+ for(int i = 0; i < 4; i++)
+ {
+ if (lRays.wasHit(i))
+ maxTsF[i] = lRays.getMinT(i);
+ else
+ maxTsF[i] = 0.0;
+ }
+ sse_t maxTsF4 = set44(maxTsF[0], maxTsF[1], maxTsF[2],
maxTsF[3]);
+ sse_t latticePositions_x = latticePositions_x =
sub4(origins_x, minBound_x);
+ latticePositions_x = mul4(latticePositions_x, cellSizeInv_x);
+ sse_t latticePositions_y = latticePositions_y =
sub4(origins_y, minBound_y);
+ latticePositions_y = mul4(latticePositions_y, cellSizeInv_y);
+ sse_t latticePositions_z = latticePositions_z =
sub4(origins_z, minBound_z);
+ latticePositions_z = mul4(latticePositions_z, cellSizeInv_z);
+ while ((bool)getmask4(cmp4_lt(currTsF4, maxTsF4)))
+ {
+ sse_union hitLatticef_x; hitLatticef_x.sse =
latticePositions_x;
+ sse_union hitLatticef_y; hitLatticef_y.sse =
latticePositions_y;
+ sse_union hitLatticef_z; hitLatticef_z.sse =
latticePositions_z;
+ int hitLattice_x3[4] = {hitLatticef_x.f[0],
hitLatticef_x.f[1], hitLatticef_x.f[2], hitLatticef_x.f[3]};
+ int hitLattice_y3[4] = {hitLatticef_y.f[0],
hitLatticef_y.f[1], hitLatticef_y.f[2], hitLatticef_y.f[3]};
+ int hitLattice_z3[4] = {hitLatticef_z.f[0],
hitLatticef_z.f[1], hitLatticef_z.f[2], hitLatticef_z.f[3]};
+
+ sse_t hitLattice_x = hitLatticef_x.sse;
cast_f2i(hitLattice_x);
+ sse_t hitLattice_y = hitLatticef_y.sse;
cast_f2i(hitLattice_y);
+ sse_t hitLattice_z = hitLatticef_z.sse;
+ cast_f2i(hitLattice_z);
+
+ sse_t absorptionCoeffs[8];
+ float absorptionCoeffsTemp[4][8];
+ //CHECK BOUNDS
+ //LOAD COEFFS
+ float coeffsTemp[4] = {0,0,0,0};
+ sse_union bounds;
+ bounds.sse = cmp4_le(currTsF4, maxTsF4);
+ bounds.sse = and4(bounds.sse, cmp4_le(zeros,
hitLattice_x));
+ bounds.sse = and4(bounds.sse, cmp4_lt(hitLattice_x,
absCoeffNx));
+ bounds.sse = and4(bounds.sse, cmp4_le(zeros,
hitLattice_y));
+ bounds.sse = and4(bounds.sse, cmp4_lt(hitLattice_y,
absCoeffNy));
+ bounds.sse = and4(bounds.sse, cmp4_le(zeros,
hitLattice_z));
+ bounds.sse = and4(bounds.sse, cmp4_lt(hitLattice_z,
absCoeffNz));
+ for(int i = 0; i < 4; i++)
+ {
+ if (bounds.f[i])
+ {
+ absorptionCoeffsTemp[i][0] =
(*_data)(hitLattice_x3[i] , hitLattice_y3[i] , hitLattice_z3[i] );
+ absorptionCoeffsTemp[i][1] =
(*_data)(hitLattice_x3[i] , hitLattice_y3[i] , hitLattice_z3[i]+1 );
+ absorptionCoeffsTemp[i][2] =
(*_data)(hitLattice_x3[i] , hitLattice_y3[i]+1 , hitLattice_z3[i] );
+ absorptionCoeffsTemp[i][3] =
(*_data)(hitLattice_x3[i] , hitLattice_y3[i]+1 , hitLattice_z3[i] +1 );
+ absorptionCoeffsTemp[i][4] =
(*_data)(hitLattice_x3[i]+1 , hitLattice_y3[i] , hitLattice_z3[i] );
+ absorptionCoeffsTemp[i][5] =
(*_data)(hitLattice_x3[i]+1 , hitLattice_y3[i] , hitLattice_z3[i]+1 );
+ absorptionCoeffsTemp[i][6] =
(*_data)(hitLattice_x3[i]+1 , hitLattice_y3[i]+1 , hitLattice_z3[i] );
+ absorptionCoeffsTemp[i][7] =
(*_data)(hitLattice_x3[i]+1 , hitLattice_y3[i]+1 , hitLattice_z3[i]+1 );
+ }
+ else
+ {
+ absorptionCoeffsTemp[i][0] = 0;
+ absorptionCoeffsTemp[i][1] = 0;
+ absorptionCoeffsTemp[i][2] = 0;
+ absorptionCoeffsTemp[i][3] = 0;
+ absorptionCoeffsTemp[i][4] = 0;
+ absorptionCoeffsTemp[i][5] = 0;
+ absorptionCoeffsTemp[i][6] = 0;
+ absorptionCoeffsTemp[i][7] = 0;
+ }
+ }
+ for(int i = 0; i <8;i++)
+ absorptionCoeffs[i] =
set44(absorptionCoeffsTemp[0][i], absorptionCoeffsTemp[1][i],
absorptionCoeffsTemp[2][i], absorptionCoeffsTemp[3][i]);
+ sse_t wxs = sub4(hitLatticef_x.sse, hitLattice_x);
+ sse_t wys = sub4(hitLatticef_y.sse, hitLattice_y);
+ sse_t wzs = sub4(hitLatticef_z.sse, hitLattice_z);
+
+ sse_t wx1s = sub4(ones, wxs);
+ sse_t wy1s = sub4(ones, wys);
+ sse_t wz1s = sub4(ones, wzs);
+
+ sse_t c[15];
+ c[0] = mul4(wz1s, absorptionCoeffs[0]);
+ c[1] = mul4(wzs, absorptionCoeffs[1]);
+ c[2] = add4(c[0], c[1]); c[2] = mul4(c[2], wy1s);
+
+ c[3] = mul4(wz1s, absorptionCoeffs[2]);
+ c[4] = mul4(wzs, absorptionCoeffs[3]);
+ c[5] = add4(c[3], c[4]); c[5] = mul4(c[5], wys);
+
+ c[6] = mul4(wz1s, absorptionCoeffs[4]);
+ c[7] = mul4(wzs, absorptionCoeffs[5]);
+ c[8] = add4(c[6], c[7]); c[8] = mul4(c[8], wy1s);
+
+ c[9] = mul4(wz1s, absorptionCoeffs[6]);
+ c[10] = mul4(wzs, absorptionCoeffs[7]);
+ c[11] = add4(c[9], c[10]); c[11] = mul4(c[11], wys);
+
+ c[12] = add4(c[2], c[5]); c[12] = mul4(c[12], wx1s);
+ c[13] = add4(c[8], c[11]); c[13] = mul4(c[13], wxs);
+ sse_union absCoeffs; absCoeffs.sse =
absorptionCoeffs[0];//add4(c[12], c[13]);
+
+
+
+ for(int i = 0; i < 4; i++)
+ {
+ if (!bounds.f[i])
+ continue;
+ if (alphas[i] > 0.9f)
+ continue;
+ if (absCoeffs.f[i] == 1.0) //TODO: FIX THIS
CHEAP HACK THAT SOLVES A BANDING ISSUE!
+ continue;
+ if (absCoeffs.f[i] == 0) //TODO: also get
rid of this
+ continue;
+ float val = (absCoeffs.f[i] -
_dataMin)/(_dataMax - _dataMin);
+ if (val > 1.0)
+ cout << "OVER" << val << endl;
+ if (val < 0.0)
+ cout << "UNDER" << val << " " <<
absCoeffs.f[i] << endl;
+ RGBAColor color = _colorMap->GetColor(val);
+ Color sourceColor = color.color;
+ float alpha_factor = color.a * (1.0 -
alphas[i]);
+ if (alpha_factor > 0.001)
+ {
+ accumColors[i] +=
sourceColor*alpha_factor;
+ }
+ alphas[i] += alpha_factor;
+ }
+
+ currTsF4 = add4(currTsF4, stepSize4);
+ latticePositions_x = add4(latticePositions_x,
steps_x);
+ latticePositions_y = add4(latticePositions_y,
steps_y);
+ latticePositions_z = add4(latticePositions_z,
steps_z);
+ }
+ lRays.resetHits();
+ context.renderer->traceRays(context, lRays);
+ for(int i = 0; i < 4; i++)
+ {
+ accumColors[i] += lRays.getColor(i)*(1.0 - alphas[i]);
+ rays.setColor(rayIndex+i, accumColors[i]);
+ }
+ }
+
+#endif
+ }
+
+//template<class T>
+ void Volume::computeHistogram(int numBuckets, int* histValues)
+ {
+ float dataMin = _dataMin;
+ float dataMax = _dataMax;
+ int nx = _data->getNx()-1;
+ int ny = _data->getNy()-1;
+ int nz = _data->getNz()-1;
+ float scale = (numBuckets-1)/(dataMax-dataMin);
+ for(int i =0; i<numBuckets;i++)
+ histValues[i] = 0;
+ for(int ix=0;ix<nx;ix++){
+ for(int iy=0;iy<ny;iy++){
+ for(int iz=0;iz<nz;iz++){
+ double p000=(*_data)(ix,iy,iz);
+ double p001=(*_data)(ix,iy,iz+1);
+ double p010=(*_data)(ix,iy+1,iz);
+ double p011=(*_data)(ix,iy+1,iz+1);
+ double p100=(*_data)(ix+1,iy,iz);
+ double p101=(*_data)(ix+1,iy,iz+1);
+ double p110=(*_data)(ix+1,iy+1,iz);
+ double p111=(*_data)(ix+1,iy+1,iz+1);
+ double min=std::min(std::min(std::min(p000,
p001), std::min(p010, p011)), std::min(std::min(p100, p101), std::min(p110,
p111)));
+ double max=std::max(std::max(std::max(p000,
p001), std::max(p010, p011)), std::max(std::max(p100, p101), std::max(p110,
p111)));
+ int nmin=(int)((min-dataMin)*scale);
+ int nmax=(int)((max-dataMin)*scale+.999999);
+ if(nmax>=numBuckets)
+ nmax=numBuckets-1;
+ if(nmin<0)
+ nmin=0;
+ for(int i=nmin;i<nmax;i++){
+ histValues[i]++;
+ }
+ }
+ }
+ }
+ }
+
+//template<class T>
+ void Volume::attenuateShadows(const RenderContext& context, RayPacket&
shadowRays) const
+ {
+
+ }
+
+//template<class T>
+float Volume::getValue(int x, int y, int z)
+{
+ return data[x + Nx*(y + Ny*z)];
+}
+
+
+
+} // namespace Manta
Added: trunk/Model/Materials/Volume.h
==============================================================================
--- (empty file)
+++ trunk/Model/Materials/Volume.h Tue Oct 30 02:53:28 2007
@@ -0,0 +1,274 @@
+/*
+ * Volume.h
+ * Author: Carson brownlee (brownleeATcs.utah.edu)
+ * Description: basic unoptomized volume renderer
+ *
+ *
+*/
+
+#ifndef VOLUME_H
+#define VOLUME_H
+
+#include <Interface/Context.h>
+#include <Interface/Material.h>
+#include <Model/Materials/OpaqueShadower.h>
+#include <Core/Geometry/Vector.h>
+#include <Model/Groups/GriddedGroup.h>
+#include <teem/nrrd.h>
+#include <Core/Thread/Barrier.h>
+#include <SCIRun/Core/Thread/Mutex.h>
+#include <Core/Util/AlignedAllocator.h>
+#include <assert.h>
+//#include "SIMD.hxx"
+
+namespace Manta
+{
+ struct RGBAColor
+ {
+ RGBAColor(Color c, float opacity) : color(c), a(opacity) {}
+ RGBAColor(float r, float g, float b, float alpha) :
color(Color(RGB(r,g,b))), a(alpha){ }
+ Color color;
+ float a; //opacity
+ };
+
+ struct ColorSlice
+ {
+ ColorSlice() : value(0),
color(RGBAColor(Color(RGBColor(0,0,0)), 1.0)) {}
+ ColorSlice(float v, RGBAColor c) : value(v), color(c) {}
+ float value; //0-1 position of slice
+ RGBAColor color;
+ };
+
+ class RGBAColorMap
+ {
+ public:
+ RGBAColorMap(unsigned int
type=RGBAColorMap::InvRainbowIso, int numSlices = 256);
+ RGBAColorMap(vector<Color>& colors, float* positions,
float* opacities = NULL, int numSlices = 256);
+ RGBAColorMap(vector<ColorSlice>& colors, int
numSlices = 256);
+ virtual ~RGBAColorMap();
+
+ void fillColor(unsigned int type);
+ virtual RGBAColor GetColor(float v); // get the
interpolated color at value v (0-1)
+ RGBAColor GetColorAbs(float v); //don't use now
+ RGBAColor GetPreInterpolatedColor(float s1, float s2,
float d);
+ void SetColors(vector<ColorSlice>& colors);
+ vector<ColorSlice> GetColors() { return _slices; }
+ ColorSlice GetSlice(int i) { return _slices.at(i); }
+ int GetNumSlices() { return _slices.size(); }
+ void SetMinMax(float min, float max);
+ void BuildSlices();
+ void PreInterpolate();
+ void computeHash();
+ void scaleAlphas(float t);
+
+ float tInc; // t increment value for volume renderers
+ float oldTInc;
+ bool _colorScaled;
+ unsigned long long course_hash;
+ vector<ColorSlice> _slices;
+ int _numSlices;
+ //float*** _opacities;
+ //Color*** _colors;
+ //ScalarTransform1D<float, Color> sTColors;
+ //ScalarTransform1D<float, float> sTOpacities;
+
+ enum {
+ InvRainbowIso=0,
+ InvRainbow,
+ Rainbow,
+ InvGreyScale,
+ InvBlackBody,
+ BlackBody,
+ GreyScale,
+ Unknown
+ };
+
+ protected:
+ float _min, _max, _invRange;
+ };
+
+ struct MANTA_ALIGN(16) Box4: public AlignedAllocator<Box4>
+ {
+ sse_t min, max;
+ sse_t diameter() const { return sub4(max,min); }
+ };
+ //std::ostream &operator<< (std::ostream &o, const Box4 &v) {
+ // o << "(" << v.min << ".." << v.max << ")";
+ // return o;
+ //}
+
+struct VMCell
+{
+// Need to make sure we have a 64 bit thing
+ unsigned long long course_hash;
+ // The max of unsigned long long is ULLONG_MAX
+ VMCell(): course_hash(0) {}
+
+ void turn_on_bits(float min, float max, float data_min, float
data_max) {
+ // We know that we have 64 bits, so figure out where min and max map
+ // into [0..63].
+ int min_index = (int)((min-data_min)/(data_max-data_min)*63);
+ /* T max_index_prep =
((max-data_min)/(data_max-data_min)*63); */
+ /* int max_index = max_index_prep-(int)max_index_prep>0?
*/
+ /* (int)max_index_prep+1:(int)max_index_prep; */
+ int max_index =
(int)ceil(double((max-data_min)/(data_max-data_min)*63));
+ // Do some checks
+
+ // This will handle clamping, I hope.
+ if (min_index > 63)
+ min_index = 63;
+ if (max_index < 0)
+ max_index = 0;
+ if (min_index < 0)
+ min_index = 0;
+ if (max_index > 63)
+ max_index = 63;
+ for (int i = min_index; i <= max_index; i++)
+ course_hash |= 1ULL << i;
+ }
+ inline VMCell& operator |= (const VMCell& v) {
+ course_hash |= v.course_hash;
+ return *this;
+ }
+ inline bool operator & (const VMCell& v) {
+ return (course_hash & v.course_hash) != 0;
+ }
+ void print(bool print_endl = true) {
+ for( int i = 0; i < 64; i++) {
+ unsigned long long bit= course_hash & (1ULL << i);
+ if (bit)
+ cout << "1";
+ else
+ cout << "0";
+ }
+ if (print_endl) cout << endl;
+ }
+};
+const int packetSize = 32; //TODO: look this up
+const int numPacklets = 32/4;
+struct RayPacketInfo {
+ sse_union minTs[numPacklets];
+ sse_union maxTs[numPacklets];
+ sse_union currTs[numPacklets];
+ sse_union alphas[numPacklets];
+ Color accumColors[numPacklets][4];
+ Color bgColor[numPacklets][4];
+ sse_union latticePositions_x[numPacklets];
+ sse_union latticePositions_y[numPacklets];
+ sse_union latticePositions_z[numPacklets];
+ sse_union steps_x[numPacklets]; //lattice steps per sample
+ sse_union steps_y[numPacklets];
+ sse_union steps_z[numPacklets];
+ sse_union dir_x[numPacklets];
+ sse_union dir_y[numPacklets];
+ sse_union dir_z[numPacklets];
+ sse_t steps_kt[numPacklets]; // stepping t values along slices
+
+};
+
+//template<class T>
+class Volume : public Material, public AlignedAllocator<Volume>
+{
+ protected:
+ class HVIsectContext {
+ public:
+ // These parameters could be modified and
hold accumulating state
+ Color total;
+ float alpha;
+ // These parameters should not change
+ int dix_dx, diy_dy, diz_dz;
+ VMCell transfunct;
+ double t_inc;
+ double t_min;
+ double t_max;
+ double t_inc_inv;
+ Ray ray;
+ //RenderContext *cx;
+ };
+public:
+ Box4 bounds;
+ sse_t sse_M, sse_N;
+ sse_t sse_one_over_N;
+ sse_t diam, // diameter of grid
+ scale, // 1 / diam
+ scaleM, // scale * M
+ scaleN; // scale * N == N / diam == number of cells per unit distance
+ //length of cell in k direction divided by u or v direction cell length.
+ sse_t pcDuStretch[3];
+ sse_t pcDvStretch[3];
+
+ int N[3]; // grid resolution
+ int M[3]; //N - 1
+ int oldN[3];
+
+ int N_mc[3]; // grid resolution
+ int M_mc[3]; //N - 1
+
+ static float resolutionFactor;
+
+ float splitLength; //above this length we split packets.
+
+ struct CellData {
+ VMCell cell;
+ };
+
+ vector <CellData> cellVector;
+
+ static void newFrame();
+ typedef unsigned long long mcT;
+ vector<mcT> cellVector_mc; // color hashes
+
+ // double cellStepSize, double metersPerUnit);
+ Volume(string dataNrrd, RGBAColorMap* colorMap, const Vector&
minBound, const Vector& maxBound, double cellStepSize, double metersPerUnit,
int depth, float forceDataMin = -FLT_MAX, float forceDataMax = -FLT_MAX);
+ void loadNrrd(string nrrd, GridArray3<double> * grid);
+ virtual ~Volume();
+
+ void getMinMax(double* min, double* max){ *min = _dataMin; *max =
_dataMax; }
+ void computeHistogram(int numBuckets, int* histValues);
+ virtual void preprocess(const PreprocessContext& context);
+ virtual void shade(const RenderContext & context, RayPacket& rays)
const;
+ virtual void attenuateShadows(const RenderContext& context,
RayPacket& shadowRays) const;
+ GridArray3<double>* _data;
+ float getValue(int x, int y, int z);
+ GridArray3<VMCell>* macrocells;
+ void calc_mcell(int depth, int ix, int iy, int iz, VMCell& mcell);
+ void parallel_calc_mcell(int cell);
+ void isect(int depth, double t_sample,
+ double dtdx, double dtdy, double dtdz,
+ double next_x, double next_y, double next_z,
+ int ix, int iy, int iz,
+ int startx, int starty, int startz,
+ const Vector& cellcorner, const Vector&
celldir,
+ HVIsectContext &isctx) const;
+ inline void SamplePacklet(RayPacketInfo& packetInfo, int packlet,
sse_t maxTs) const;
+ Vector diag;
+ Vector inv_diag;
+protected:
+ int _nx,_ny,_nz;
+ Nrrd* nrrd;
+ float* data;
+ int Nx, Ny, Nz; // dimensions
+ RGBAColorMap* _colorMap;
+ Vector _minBound;
+ Vector _maxBound;
+ Vector _datadiag;
+ Vector _sdiag;
+ Vector _hierdiag;
+ Vector _ihierdiag;
+ int _depth;
+ int* _xsize,*_ysize,*_zsize;
+ double* _ixsize, *_iysize, *_izsize;
+ double _metersPerUnit;
+ Vector _cellSize;
+ double _stepSize;
+ float _dataMin;
+ float _dataMax;
+ float _colorScalar;
+ float _maxDistance;
+
+ int done_count;
+};
+
+}
+
+#endif
Modified: trunk/scenes/CMakeLists.txt
==============================================================================
--- trunk/scenes/CMakeLists.txt (original)
+++ trunk/scenes/CMakeLists.txt Tue Oct 30 02:53:28 2007
@@ -54,6 +54,14 @@
TARGET_LINK_LIBRARIES(scene_fence ${MANTA_SCENE_LINK})
ENDIF(SCENE_FENCE)
+IF(FOUND_TEEM)
+SET(SCENE_VOLUMETEST 0 CACHE BOOL "volume test")
+IF(SCENE_VOLUMETEST)
+ ADD_LIBRARY(scene_volumeTest volumeTest.cc)
+ TARGET_LINK_LIBRARIES(scene_volumeTest ${MANTA_SCENE_LINK})
+ENDIF(SCENE_VOLUMETEST)
+ENDIF(FOUND_TEEM)
+
# octree isosurface
SET(SCENE_OCTISOVOL 0 CACHE BOOL "octree isosurface")
IF(SCENE_OCTISOVOL)
Added: trunk/scenes/volumeTest.cc
==============================================================================
--- (empty file)
+++ trunk/scenes/volumeTest.cc Tue Oct 30 02:53:28 2007
@@ -0,0 +1,268 @@
+/* File: volumeTest.cc
+ * Author: Carson Brownlee (brownleeATcs.utah.edu)
+ * Description: test for basic single ray volume renderer
+ * example usage: bin/manta -scene "lib/libscene_volumeTest.so (-i
/home/collab/brownlee/Data-Explosion/vol/temp_CC_M02_0496.nrrd -dataMinMax
300 2100)" -np 4
+
+ *
+*/
+
+
+#include <Model/Materials/Volume.h>
+
+#include <Core/Color/ColorDB.h>
+#include <Core/Color/Color.h>
+#include <Core/Geometry/Vector.h>
+#include <Core/Exceptions/IllegalArgument.h>
+#include <Core/Util/Args.h>
+#include <Interface/Context.h>
+#include <Interface/LightSet.h>
+#include <Interface/MantaInterface.h>
+#include <Interface/Scene.h>
+#include <Model/Readers/PlyReader.h>
+#include <Model/AmbientLights/ArcAmbient.h>
+#include <Model/Backgrounds/LinearBackground.h>
+#include <Model/Backgrounds/ConstantBackground.h>
+#include <Model/Groups/Group.h>
+#include <Model/Lights/PointLight.h>
+#include <Model/Materials/Lambertian.h>
+#include <Model/Materials/MetalMaterial.h>
+#include <Model/Materials/Phong.h>
+#include <Model/Materials/Checker.h>
+#include <Model/Materials/Flat.h>
+#include <Model/Primitives/Parallelogram.h>
+#include <Model/Primitives/Sphere.h>
+#include <Model/Primitives/Cube.h>
+#include <Model/Primitives/Disk.h>
+#include <Model/Primitives/Parallelogram.h>
+#include <Model/Primitives/SuperEllipsoid.h>
+#include <Model/Primitives/Heightfield.h>
+#include <Model/Primitives/Hemisphere.h>
+#include <Model/MiscObjects/Intersection.h>
+#include <Model/MiscObjects/Difference.h>
+#include <Model/Textures/CheckerTexture.h>
+#include <Model/Textures/MarbleTexture.h>
+//#include <Model/Textures/NormalTexture.h>
+#include <Model/Textures/WoodTexture.h>
+#include <Model/Textures/OakTexture.h>
+#include <Model/TexCoordMappers/UniformMapper.h>
+#include <Model/TexCoordMappers/SphericalMapper.h>
+#include <Model/TexCoordMappers/LinearMapper.h>
+#include <Core/Geometry/AffineTransform.h>
+#include <Core/Util/NotFinished.h>
+#include <Model/Groups/GriddedGroup.h>
+#include <Model/Instances/InstanceRT.h>
+#include <Model/Instances/InstanceT.h>
+#include <Model/Instances/InstanceST.h>
+#include <Model/Instances/InstanceRST.h>
+#include <Model/Instances/Instance.h>
+#include <Model/Cameras/PinholeCamera.h>
+#include <Core/Color/RegularColorMap.h>
+#include <Model/Readers/ParticleNRRD.h>
+
+#include <Engine/Factory/Factory.h>
+
+#include <Model/MiscObjects/CuttingPlane.h>
+
+#include <sgi_stl_warnings_off.h>
+#include <iostream>
+#include <sgi_stl_warnings_on.h>
+
+#include <math.h>
+#include <string>
+#include <vector>
+#include <fstream>
+
+using namespace Manta;
+using namespace std;
+
+extern "C"
+Scene* make_scene(const ReadContext& context, const vector<string>& args)
+{
+ double t_inc = 0.00125;
+ double x,y,z;
+ double dataMin = -FLT_MAX;
+ double dataMax = -FLT_MAX;
+ Vector minBound(-0.001, -0.101, -0.001);
+ Vector maxBound(0.101, 0.201, 0.101);
+ string cMapType = "";
+ string cMapFile = "";
+ int depth = 2;
+ string volFile = "";
+ RGBAColorMap* cMap = NULL;
+ size_t argc = args.size();
+ vector<ColorSlice> slices;
+ cout << "argc: " << argc << endl;
+ for(size_t i = 0; i < argc; ++i) {
+ string arg = args[i];
+ if(arg == "-t")
+ {
+ if(!getDoubleArg(i, args, t_inc))
+ throw IllegalArgument("scene volumeTest -t",i,args);
+ }
+ else if(arg == "-dataMinMax")
+ {
+ if(!(getDoubleArg(i, args, dataMin) && getDoubleArg(i, args,
dataMax)))
+ throw IllegalArgument("scene volumeTest -dataMinMax",i,args);
+ cout << "datamin/max: " << dataMin << " " << dataMax << endl;
+ }
+ else if(arg == "-minBound")
+ {
+ if(!(getDoubleArg(i, args, x) && !getDoubleArg(i, args, y) &&
!getDoubleArg(i, args, z)))
+ throw IllegalArgument("scene volumeTest -minBound",i,args);
+ minBound = Vector(x,y,z);
+ }
+ else if(arg == "-maxBound")
+ {
+ if(!(getDoubleArg(i, args, x) && !getDoubleArg(i, args, y) &&
!getDoubleArg(i, args, z)))
+ throw IllegalArgument("scene volumeTest -maxBound",i,args);
+ maxBound = Vector(x,y,z);
+ }
+ else if(arg == "-i")
+ {
+ if(!getStringArg(i, args, volFile))
+ throw IllegalArgument("scene volumeTest -i", i, args);
+ }
+ else if (arg == "-depth") {
+ if (!getIntArg(i, args, depth))
+ throw IllegalArgument("scene volumeTest -depth", i, args);
+ cout << "i: " << i << endl;
+ }
+ else if (arg == "-cMapFile")
+ {
+ if(!getStringArg(i, args, cMapFile))
+ throw IllegalArgument("scene volumeTest -cMapFile",i,args);
+ ifstream in(cMapFile.c_str());
+ float pos,r,g,b,a;
+ while (in >> pos && in >> r && in >> g && in >> b && in >> a)
+ slices.push_back(ColorSlice(pos, RGBAColor(r,g,b,a)));
+ in.close();
+ cMap = new RGBAColorMap(slices);
+ }
+ else if (arg == "-cMapType")
+ {
+ if(!getStringArg(i, args, cMapType))
+ throw IllegalArgument("scene volumeTest -cMapType",i,args);
+ if(cMapType == "Rainbow")
+ cMap = new RGBAColorMap(RGBAColorMap::Rainbow);
+ else if(cMapType == "InvRainbow")
+ cMap = new RGBAColorMap(RGBAColorMap::InvRainbow);
+ else if(cMapType == "BlackBody")
+ cMap = new RGBAColorMap(RGBAColorMap::BlackBody);
+ else if(cMapType == "InvBlackBody")
+ cMap = new RGBAColorMap(RGBAColorMap::InvBlackBody);
+ else if(cMapType == "GreyScale")
+ cMap = new RGBAColorMap(RGBAColorMap::GreyScale);
+ else if(cMapType == "InvGreyScale")
+ cMap = new RGBAColorMap(RGBAColorMap::InvGreyScale);
+ else if(cMapType == "InvRainbowIso")
+ cMap = new RGBAColorMap(RGBAColorMap::InvRainbowIso);
+ else
+ throw IllegalArgument("scene volumeTest -cMapType",i,args);
+
+ } else {
+ cerr<<"Valid options for scene pcgt:\n";
+ cerr<<" -depth <int> depth of macrocells (1-5 is good, 2
default)\n";
+ cerr<<" -i <string> input filename\n";
+ cerr<<" -minBound <float> <float> <float> minimum corner of bounding
box\n";
+ cerr<<" -maxBound <float> <float> <float> maximum corner of bounding
box\n";
+ cerr<<" -t <float> t increment value for volume\n";
+ cerr<<" -cMapType <string> type of colormap:
Rainbow/InvRainbow/BlackBody/InvBlackBody/GreyScale/InvRainbowISo/InvGreyScale\n";
+ cerr<<" -cMapFile <string> file containing colormap info, in format
of position(0-1) r g b a\n";
+ throw IllegalArgument("scene volumeTest", i, args);
+ }
+ }
+
+ Group* group = new Group();
+
+ float div = 1.0/255.0;
+ float a = 1.0;
+
+/*
+ slices.push_back(ColorSlice(0.0f, RGBAColor(Color(RGB(0, 0, 0)) * div,
0*a)));
+ slices.push_back(ColorSlice(0.109804f, RGBAColor(Color(RGB(52, 0, 0)) *
div, 0*a)));
+ slices.push_back(ColorSlice(0.2f, RGBAColor(Color(RGB(102, 2, 0)) * div,
0.1*a)));
+ slices.push_back(ColorSlice(0.328571f, RGBAColor(Color(RGB(153, 18, 0)) *
div, 0.216667*a)));
+ slices.push_back(ColorSlice(0.4f, RGBAColor(Color(RGB(200, 41, 0)) * div,
0.23*a)));
+ slices.push_back(ColorSlice(0.5f, RGBAColor(Color(RGB(230, 71, 0)) * div,
0.27*a)));
+ slices.push_back(ColorSlice(0.618367f, RGBAColor(Color(RGB(255, 120, 0)) *
div, 0.3375*a)));
+ slices.push_back(ColorSlice(0.68f, RGBAColor(Color(RGB(255, 163, 20)) *
div, 0.35*a)));
+ slices.push_back(ColorSlice(0.72f, RGBAColor(Color(RGB(255, 204, 55)) *
div, 0.37*a)));
+ slices.push_back(ColorSlice(0.79f, RGBAColor(Color(RGB(255, 228, 80)) *
div, 0.39*a)));
+ slices.push_back(ColorSlice(0.85f, RGBAColor(Color(RGB(255, 247, 120)) *
div, 0.43*a)));
+ slices.push_back(ColorSlice(0.92f, RGBAColor(Color(RGB(255, 255, 180)) *
div, 0.47*a)));
+ slices.push_back(ColorSlice(1.0f, RGBAColor(Color(RGB(255, 255, 255)) *
div, 0.5*a)));
+*/
+if(cMap == NULL)
+{
+ slices.push_back(ColorSlice(0.0f, RGBAColor(Color(RGB(0, 0, 0)) * div,
0*a)));
+ slices.push_back(ColorSlice(0.11f, RGBAColor(Color(RGB(52, 0, 0)) * div,
0*a)));
+ slices.push_back(ColorSlice(0.4f, RGBAColor(Color(RGB(200, 41, 0)) * div,
0.23*a)));
+ slices.push_back(ColorSlice(0.5f, RGBAColor(Color(RGB(230, 71, 0)) * div,
0.27*a)));
+ slices.push_back(ColorSlice(0.618367f, RGBAColor(Color(RGB(255, 120, 0)) *
div, 0.3375*a)));
+ slices.push_back(ColorSlice(0.68f, RGBAColor(Color(RGB(255, 163, 20)) *
div, 0.35*a)));
+
+//slices.push_back(ColorSlice(0.5f, RGBAColor(Color(RGB(255, 163, 20)) *
div, 0*a)));
+//slices.push_back(ColorSlice(0.51f, RGBAColor(Color(RGB(255, 163, 20)) *
div, 0.35*a)));
+
+
+ slices.push_back(ColorSlice(0.72f, RGBAColor(Color(RGB(255, 204, 55)) *
div, 0.37*a)));
+ slices.push_back(ColorSlice(0.79f, RGBAColor(Color(RGB(255, 228, 80)) *
div, 0.39*a)));
+ slices.push_back(ColorSlice(0.85f, RGBAColor(Color(RGB(255, 247, 120)) *
div, 0.43*a)));
+ slices.push_back(ColorSlice(0.92f, RGBAColor(Color(RGB(255, 255, 180)) *
div, 0.47*a)));
+ slices.push_back(ColorSlice(1.0f, RGBAColor(Color(RGB(255, 255, 255)) *
div, 0.5*a)));
+
+ /*Array1<Color> matls;
+ Array1<AlphaPos> alphas;
+ matls.add(Color(RGB(0, 0, 0)) * div);
+ matls.add(Color(RGB(52, 0, 0)) * div);
+ matls.add(Color(RGB(102, 2, 0)) * div);
+ matls.add(Color(RGB(153, 18, 0)) * div);
+ matls.add(Color(RGB(200, 41, 0)) * div);
+ matls.add(Color(RGB(230, 71, 0)) * div);
+ matls.add(Color(RGB(255, 120, 0)) * div);
+ matls.add(Color(RGB(255, 163, 20)) * div);
+ matls.add(Color(RGB(255, 204, 55)) * div);
+ matls.add(Color(RGB(255, 228, 80)) * div);
+ matls.add(Color(RGB(255, 247, 120)) * div);
+ matls.add(Color(RGB(255, 255, 180)) * div);
+ matls.add(Color(RGB(255, 255, 255)) * div);
+
+
+ alphas.add(AlphaPos(0 , 0)); // pos 0
+ alphas.add(AlphaPos(0.109804, 0));
+ alphas.add(AlphaPos(0.328571, 0.216667));
+ alphas.add(AlphaPos(0.618367, 0.3375));
+ alphas.add(AlphaPos(1, 0.5));
+ CDColorMap2* cmap2 = new CDColorMap2(matls, alphas, 256, t_inc);
+ */
+ cMap = new RGBAColorMap(slices);
+};
+
+
+cMap->scaleAlphas(t_inc);
+Volume* mat = new Volume(volFile, cMap, minBound, maxBound, t_inc, 10.0,
depth, dataMin, dataMax);
+ //double min, max;
+ //mat->getMinMax(&min, &max);
+ //cmap2->setMinMax(min, max);
+ Primitive* prim = new Cube(mat, minBound - Vector(0.001, 0.001,
0.001), maxBound + Vector(0.001, 0.001, 0.001));
+ group->add(prim);
+
+
+ Scene* scene = new Scene();
+ scene->setBackground(new LinearBackground(Color(RGB(0.2, 0.4, 0.9)),
+ Color(RGB(0.0,0.0,0.0)),
+ Vector(0,1,0)));
+ scene->setObject(group);
+
+ LightSet* lights = new LightSet();
+ lights->add(new PointLight(Vector(5,5,8), Color(RGB(1,1,1))*2));
+ Color cup(RGB(0.1, 0.3, 0.8));
+ Color cdown(RGB(0.82, 0.62, 0.62));
+ Vector up(0,1,0);
+ lights->setAmbientLight(new ArcAmbient(cup, cdown, up));
+ scene->setLights(lights);
+
+
+ return scene;
+}
+
- [Manta] r1792 - in trunk: Model/Materials scenes, brownlee, 10/30/2007
Archive powered by MHonArc 2.6.16.
