-
sparker
thiago
Manta Interactive Ray Tracer Development Mailing List
Text archives Help
- From: brownlee@sci.utah.edu
- To: manta@sci.utah.edu
- Subject: [Manta] r1827 - in trunk: Model/Materials Model/Readers scenes
- Date: Tue, 6 Nov 2007 15:33:57 -0700 (MST)
Author: brownlee
Date: Tue Nov 6 15:33:55 2007
New Revision: 1827
Added:
trunk/Model/Readers/VolumeNRRD-stub.cc
trunk/Model/Readers/VolumeNRRD.cc
trunk/Model/Readers/VolumeNRRD.h
Modified:
trunk/Model/Materials/CMakeLists.txt
trunk/Model/Materials/Volume.cc
trunk/Model/Materials/Volume.h
trunk/Model/Readers/CMakeLists.txt
trunk/scenes/volumeTest.cc
Log:
templated, cleaned Volume Code, seperated it a bit from Nrrd
Modified: trunk/Model/Materials/CMakeLists.txt
==============================================================================
--- trunk/Model/Materials/CMakeLists.txt (original)
+++ trunk/Model/Materials/CMakeLists.txt Tue Nov 6 15:33:55 2007
@@ -33,16 +33,11 @@
Materials/Transparent.h
)
-IF(FOUND_TEEM AND MANTA_SSE)
+IF(MANTA_SSE)
SET (Manta_Materials_SRCS
${Manta_Materials_SRCS}
Materials/Volume.h
Materials/Volume.cc
- )
- INCLUDE_DIRECTORIES(${TEEM_INCLUDE_DIRS})
- LINK_DIRECTORIES (${TEEM_LIBRARY_DIRS})
- SET(Manta_Model_extra_libs ${Manta_Model_extra_libs}
- ${TEEM_LIBRARY}
- )
-ENDIF(FOUND_TEEM AND MANTA_SSE)
+ )
+ENDIF(MANTA_SSE)
Modified: trunk/Model/Materials/Volume.cc
==============================================================================
--- trunk/Model/Materials/Volume.cc (original)
+++ trunk/Model/Materials/Volume.cc Tue Nov 6 15:33:55 2007
@@ -1,1706 +1,399 @@
-#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)
+
+ bool lessThan(const ColorSlice& left, const ColorSlice& right)
{
- _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";
+ return left.value <= right.value;
}
-
- void Volume::calc_mcell(int depth, int startx, int starty, int startz,
VMCell& mcell)
+
+ RGBAColorMap::RGBAColorMap(unsigned int type, int numSlices)
+ :oldTInc(0.01), _colorScaled(false)
{
- 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";
+ 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);
+ }
}
-// template<class T>
- Volume::~Volume()
+ 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();
}
-
-// template<class T>
- void Volume::preprocess(const PreprocessContext& context)
+
+ 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();
}
-// template<class T>
- void Volume::loadNrrd(string file, GridArray3<double> * grid)
+ RGBAColor RGBAColorMap::GetColor(float v) // get the interpolated color
at value v (0-1)
{
- 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__);
- }
+ 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);
+
+ }
- // 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";
+ RGBAColor RGBAColorMap::GetColorAbs(float v)
+ {
+ return GetColor((v - _min)/_invRange);
}
-
-
-#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
+
+ 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);
+ }
- 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;
- }
+ 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);
}
-//#define CD_SSE
-#define RAY_TERMINATION 0.9
-//template<class T>
- void Volume::shade(const RenderContext & context, RayPacket& rays) const
+ void RGBAColorMap::BuildSlices()
{
- bool BITMASK_VERSION = true;
-if (BITMASK_VERSION)
-{
- rays.normalizeDirections();
- double t_inc = _stepSize;
+ // 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
+ //
+ // }
+ }
- 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));
- }
+ 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"; */
+ }
- 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]);
- }
- }
+ void RGBAColorMap::computeHash()
+ {
-#endif
- }
+ unsigned long long new_course_hash = 0;
+ int nindex = 63; // 64 - 1
-//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]++;
- }
- }
+ 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;
}
}
- }
-//template<class T>
- void Volume::attenuateShadows(const RenderContext& context, RayPacket&
shadowRays) const
- {
+ 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";
}
-//template<class T>
-float Volume::getValue(int x, int y, int z)
-{
- return data[x + Nx*(y + Ny*z)];
-}
+ //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;
+ }
+} //namespace manta
-} // namespace Manta
Modified: trunk/Model/Materials/Volume.h
==============================================================================
--- trunk/Model/Materials/Volume.h (original)
+++ trunk/Model/Materials/Volume.h Tue Nov 6 15:33:55 2007
@@ -4,7 +4,7 @@
* Description: basic unoptomized volume renderer
*
*
-*/
+ */
#ifndef VOLUME_H
#define VOLUME_H
@@ -14,261 +14,1486 @@
#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 <Model/Readers/VolumeNRRD.h>
+#include <Interface/Scene.h>
+#include <Interface/Renderer.h>
+#include <Interface/RayPacket.h>
+#include <cassert>
+#include <float.h>
+
+#include <SCIRun/Core/Thread/Time.h>
+#include <SCIRun/Core/Thread/Thread.h>
+#include <Core/Geometry/BBox.h>
+
#include <assert.h>
//#include "SIMD.hxx"
namespace Manta
{
- struct RGBAColor
+ 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<T> >
+ {
+ 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
+
+ // float cellStepSize, float metersPerUnit);
+ Volume(string dataNrrd, RGBAColorMap* colorMap, const Vector&
minBound, const Vector& maxBound, double cellStepSize, double metersPerUnit,
int depth, double forceDataMin = -FLT_MAX, double forceDataMax = -FLT_MAX);
+ 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<T>* _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;
+ double* 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;
+ float _metersPerUnit;
+ Vector _cellSize;
+ float _stepSize;
+ float _dataMin;
+ float _dataMax;
+ float _colorScalar;
+ float _maxDistance;
+
+ int done_count;
+ };
+
+ template<class T>
+ Volume<T>::Volume(string dataNrrd, RGBAColorMap* colorMap, const Vector&
minBound, const Vector& maxBound, double cellStepSize, double metersPerUnit,
int depth, double forceDataMin, double forceDataMax)
+ :_colorMap(colorMap)
+ {
+ //_data = new GridArray3<double>();
+ _data = loadNRRDToGrid<T>(dataNrrd);
+ //_emittedLight = new GridArray3<float>();
+ //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();
+
+
+ float 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++)
{
- 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 ((*_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";
+ }
+
+ template<class T>
+ Volume<T>::~Volume()
+ {
+
+ }
+
+
+ template<class T>
+ float Volume<T>::resolutionFactor = 5.0f;
+
+ template<class T>
+ void computeMinMax(float& min, float& max, GridArray3<T>& 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>
+ void Volume<T>::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<T>& 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>
+ void Volume<T>::preprocess(const PreprocessContext& context)
+ {
+
+ }
+
+
+#define RAY_TERMINATION_THRESHOLD 0.9
+ template<class T>
+ void Volume<T>::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<T>& 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;
+ float dx = ly2 - ly1;
+
+ float dy, dy1, dy2;
+ dy1 = lz2 - lz1;
+ dy2 = lz4 - lz3;
+ dy = dy1 * (1 - x_weight_high) + dy2 * x_weight_high;
+
+ float 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;
+ float length2 = dx*dx+dy*dy+dz*dz;
+ if (length2){
+ // this lets the compiler use a special 1/sqrt() operation
+ float 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<T>::shade(const RenderContext & context, RayPacket& rays)
const
+ {
+ bool BITMASK_VERSION = true;
+ if (BITMASK_VERSION)
+ {
+ rays.normalizeDirections();
+ float 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();
+ float 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<T>& 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(float 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;
+
+
+ float 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;
+ float 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;
}
- if (print_endl) cout << endl;
+ }
+ lRays.resetHits();
+ context.renderer->traceRays(context, lRays);
+ for(int i = 0; i < 4; i++)
+ {
+ rays.setColor(rayIndex+i, accumColors[i]);
+ //}
+ }
}
-};
-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
-
-};
+#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]};
-//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;
+ 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<T>::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++){
+ float p000=(*_data)(ix,iy,iz);
+ float p001=(*_data)(ix,iy,iz+1);
+ float p010=(*_data)(ix,iy+1,iz);
+ float p011=(*_data)(ix,iy+1,iz+1);
+ float p100=(*_data)(ix+1,iy,iz);
+ float p101=(*_data)(ix+1,iy,iz+1);
+ float p110=(*_data)(ix+1,iy+1,iz);
+ float p111=(*_data)(ix+1,iy+1,iz+1);
+ float min=std::min(std::min(std::min(p000, p001), std::min(p010,
p011)), std::min(std::min(p100, p101), std::min(p110, p111)));
+ float 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<T>::attenuateShadows(const RenderContext& context,
RayPacket& shadowRays) const
+ {
- int done_count;
-};
+ }
-}
+ template<class T>
+ float Volume<T>::getValue(int x, int y, int z)
+ {
+ return data[x + Nx*(y + Ny*z)];
+ }
+}//namespace manta
#endif
Modified: trunk/Model/Readers/CMakeLists.txt
==============================================================================
--- trunk/Model/Readers/CMakeLists.txt (original)
+++ trunk/Model/Readers/CMakeLists.txt Tue Nov 6 15:33:55 2007
@@ -15,6 +15,23 @@
INCLUDE_DIRECTORIES(/usr/include/malloc/)
ENDIF (APPLE)
+IF(FOUND_TEEM)
+ SET(Manta_Readers_SRCS
+ ${Manta_Readers_SRCS}
+ Readers/VolumeNRRD.h
+ Readers/VolumeNRRD.cc
+ )
+ INCLUDE_DIRECTORIES(${TEEM_INCLUDE_DIRS})
+ LINK_DIRECTORIES (${TEEM_LIBRARY_DIRS})
+ELSE (FOUND_TEEM)
+ SET(Manta_Readers_SRCS
+ ${Manta_Readers_SRCS}
+ Readers/VolumeNRRD.h
+ Readers/VolumeNRRD-stub.cc
+ )
+ENDIF(FOUND_TEEM)
+
+
# Reader for NRRD particle data
IF(BUILD_NRRDPARTICLES)
SET(Manta_Readers_SRCS
Added: trunk/Model/Readers/VolumeNRRD-stub.cc
==============================================================================
--- (empty file)
+++ trunk/Model/Readers/VolumeNRRD-stub.cc Tue Nov 6 15:33:55 2007
@@ -0,0 +1,10 @@
+#include <Model/Readers/VolumeNRRD.h>
+
+using namespace Manta;
+
+template<class T>
+GridArray3<T>* loadNRRDToGrid(string file)
+{
+ cerr << "loadNRRDToGrid: Error: TEEM not supported\n";
+ return NULL;
+}
Added: trunk/Model/Readers/VolumeNRRD.cc
==============================================================================
--- (empty file)
+++ trunk/Model/Readers/VolumeNRRD.cc Tue Nov 6 15:33:55 2007
@@ -0,0 +1,85 @@
+#include <Model/Readers/VolumeNRRD.h>
+#include <teem/nrrd.h>
+
+using namespace Manta;
+
+template<class T>
+GridArray3<T>* loadNRRDToGrid(string file)
+{
+ GridArray3<T>* grid = new GridArray3<T>();
+ 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";
+ return grid;
+}
+
Added: trunk/Model/Readers/VolumeNRRD.h
==============================================================================
--- (empty file)
+++ trunk/Model/Readers/VolumeNRRD.h Tue Nov 6 15:33:55 2007
@@ -0,0 +1,20 @@
+/*
+ *Author: Carson Brownlee (brownleeATcs.utah.edu)
+ *Description: simple loader of Nrrd files into GridArray3
+ *
+*/
+
+#ifndef VOLUMENRRD_H
+#define VOLUMENRRD_H
+
+#include <Model/Groups/GriddedGroup.h>
+#include <iostream>
+using namespace std;
+
+namespace Manta
+{
+template<class T>
+GridArray3<T>* loadNRRDToGrid(string file);
+}
+
+#endif
Modified: trunk/scenes/volumeTest.cc
==============================================================================
--- trunk/scenes/volumeTest.cc (original)
+++ trunk/scenes/volumeTest.cc Tue Nov 6 15:33:55 2007
@@ -18,7 +18,6 @@
#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>
@@ -57,7 +56,7 @@
#include <Model/Instances/Instance.h>
#include <Model/Cameras/PinholeCamera.h>
#include <Core/Color/RegularColorMap.h>
-#include <Model/Readers/ParticleNRRD.h>
+//#include <Model/Readers/ParticleNRRD.h>
#include <Engine/Factory/Factory.h>
@@ -240,7 +239,7 @@
cMap->scaleAlphas(t_inc);
-Volume* mat = new Volume(volFile, cMap, minBound, maxBound, t_inc, 10.0,
depth, dataMin, dataMax);
+Volume<float>* mat = new Volume<float>(volFile, cMap, minBound, maxBound,
t_inc, 10.0, depth, dataMin, dataMax);
//double min, max;
//mat->getMinMax(&min, &max);
//cmap2->setMinMax(min, max);
- [Manta] r1827 - in trunk: Model/Materials Model/Readers scenes, brownlee, 11/06/2007
Archive powered by MHonArc 2.6.16.
