root/OpenSceneGraph/trunk/examples/osgvolume/osgvolume.cpp @ 8856

Revision 8856, 77.0 kB (checked in by robert, 6 years ago)

Added optional usage of DCMTK in the dicom plugin

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
Line 
1/* OpenSceneGraph example, osgvolume.
2*
3*  Permission is hereby granted, free of charge, to any person obtaining a copy
4*  of this software and associated documentation files (the "Software"), to deal
5*  in the Software without restriction, including without limitation the rights
6*  to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
7*  copies of the Software, and to permit persons to whom the Software is
8*  furnished to do so, subject to the following conditions:
9*
10*  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
11*  IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
12*  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
13*  AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
14*  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
15*  OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
16*  THE SOFTWARE.
17*/
18
19#include <osg/Node>
20#include <osg/Geometry>
21#include <osg/Notify>
22#include <osg/Texture3D>
23#include <osg/TexGen>
24#include <osg/Geode>
25#include <osg/Billboard>
26#include <osg/PositionAttitudeTransform>
27#include <osg/ClipNode>
28#include <osg/AlphaFunc>
29#include <osg/TexGenNode>
30#include <osg/TexEnv>
31#include <osg/TexEnvCombine>
32#include <osg/Material>
33#include <osg/PrimitiveSet>
34#include <osg/Endian>
35#include <osg/BlendFunc>
36#include <osg/BlendEquation>
37#include <osg/TransferFunction>
38
39#include <osgDB/Registry>
40#include <osgDB/ReadFile>
41#include <osgDB/WriteFile>
42#include <osgDB/FileUtils>
43#include <osgDB/FileNameUtils>
44
45#include <osgGA/EventVisitor>
46
47#include <osgUtil/CullVisitor>
48
49#include <osgViewer/Viewer>
50#include <osgViewer/ViewerEventHandlers>
51
52#include <osg/io_utils>
53
54#include <algorithm>
55#include <iostream>
56
57
58typedef std::vector< osg::ref_ptr<osg::Image> > ImageList;
59
60//  example ReadOperator
61// struct ReadOperator
62// {
63//     inline void luminance(float l) const { rgba(l,l,l,1.0f); }
64//     inline void alpha(float a) const { rgba(1.0f,1.0f,1.0f,a); }
65//     inline void luminance_alpha(float l,float a) const { rgba(l,l,l,a); }
66//     inline void rgb(float r,float g,float b) const { rgba(r,g,b,1.0f); }
67//     inline void rgba(float r,float g,float b,float a) const { std::cout<<"pixel("<<r<<", "<<g<<", "<<b<<", "<<a<<")"<<std::endl; }
68// };
69
70
71template <typename T, class O>   
72void _readRow(unsigned int num, GLenum pixelFormat, T* data,float scale, const O& operation)
73{
74    switch(pixelFormat)
75    {
76        case(GL_LUMINANCE):         { for(unsigned int i=0;i<num;++i) { float l = float(*data++)*scale; operation.luminance(l); } }  break;
77        case(GL_ALPHA):             { for(unsigned int i=0;i<num;++i) { float a = float(*data++)*scale; operation.alpha(a); } }  break;
78        case(GL_LUMINANCE_ALPHA):   { for(unsigned int i=0;i<num;++i) { float l = float(*data++)*scale; float a = float(*data++)*scale; operation.luminance_alpha(l,a); } }  break;
79        case(GL_RGB):               { for(unsigned int i=0;i<num;++i) { float r = float(*data++)*scale; float g = float(*data++)*scale; float b = float(*data++)*scale; operation.rgb(r,g,b); } }  break;
80        case(GL_RGBA):              { for(unsigned int i=0;i<num;++i) { float r = float(*data++)*scale; float g = float(*data++)*scale; float b = float(*data++)*scale; float a = float(*data++)*scale; operation.rgba(r,g,b,a); } }  break;
81        case(GL_BGR):               { for(unsigned int i=0;i<num;++i) { float b = float(*data++)*scale; float g = float(*data++)*scale; float r = float(*data++)*scale; operation.rgb(r,g,b); } }  break;
82        case(GL_BGRA):              { for(unsigned int i=0;i<num;++i) { float b = float(*data++)*scale; float g = float(*data++)*scale; float r = float(*data++)*scale; float a = float(*data++)*scale; operation.rgba(r,g,b,a); } }  break;
83    }
84}
85
86template <class O>   
87void readRow(unsigned int num, GLenum pixelFormat, GLenum dataType, unsigned char* data, const O& operation)
88{
89    switch(dataType)
90    {
91        case(GL_BYTE):              _readRow(num,pixelFormat, (char*)data,            1.0f/128.0f,        operation); break;
92        case(GL_UNSIGNED_BYTE):     _readRow(num,pixelFormat, (unsigned char*)data,   1.0f/255.0f,        operation); break;
93        case(GL_SHORT):             _readRow(num,pixelFormat, (short*) data,          1.0f/32768.0f,      operation); break;
94        case(GL_UNSIGNED_SHORT):    _readRow(num,pixelFormat, (unsigned short*)data,  1.0f/65535.0f,      operation); break;
95        case(GL_INT):               _readRow(num,pixelFormat, (int*) data,            1.0f/2147483648.0f, operation); break;
96        case(GL_UNSIGNED_INT):      _readRow(num,pixelFormat, (unsigned int*) data,   1.0f/4294967295.0f, operation); break;
97        case(GL_FLOAT):             _readRow(num,pixelFormat, (float*) data,          1.0f,               operation); break;
98    }
99}
100
101template <class O>   
102void readImage(osg::Image* image, const O& operation)
103{
104    if (!image) return;
105   
106    for(int r=0;r<image->r();++r)
107    {
108        for(int t=0;t<image->t();++t)
109        {
110            readRow(image->s(), image->getPixelFormat(), image->getDataType(), image->data(0,t,r), operation);
111        }
112    }
113}
114
115//  example ModifyOperator
116// struct ModifyOperator
117// {
118//     inline void luminance(float& l) const {}
119//     inline void alpha(float& a) const {}
120//     inline void luminance_alpha(float& l,float& a) const {}
121//     inline void rgb(float& r,float& g,float& b) const {}
122//     inline void rgba(float& r,float& g,float& b,float& a) const {}
123// };
124
125
126template <typename T, class M>   
127void _modifyRow(unsigned int num, GLenum pixelFormat, T* data,float scale, const M& operation)
128{
129    float inv_scale = 1.0f/scale;
130    switch(pixelFormat)
131    {
132        case(GL_LUMINANCE):         { for(unsigned int i=0;i<num;++i) { float l = float(*data)*scale; operation.luminance(l); *data++ = T(l*inv_scale); } }  break;
133        case(GL_ALPHA):             { for(unsigned int i=0;i<num;++i) { float a = float(*data)*scale; operation.alpha(a); *data++ = T(a*inv_scale); } }  break;
134        case(GL_LUMINANCE_ALPHA):   { for(unsigned int i=0;i<num;++i) { float l = float(*data)*scale; float a = float(*(data+1))*scale; operation.luminance_alpha(l,a); *data++ = T(l*inv_scale); *data++ = T(a*inv_scale); } }  break;
135        case(GL_RGB):               { for(unsigned int i=0;i<num;++i) { float r = float(*data)*scale; float g = float(*(data+1))*scale; float b = float(*(data+2))*scale; operation.rgb(r,g,b); *data++ = T(r*inv_scale); *data++ = T(g*inv_scale); *data++ = T(b*inv_scale); } }  break;
136        case(GL_RGBA):              { for(unsigned int i=0;i<num;++i) { float r = float(*data)*scale; float g = float(*(data+1))*scale; float b = float(*(data+2))*scale; float a = float(*(data+3))*scale; operation.rgba(r,g,b,a); *data++ = T(r*inv_scale); *data++ = T(g*inv_scale); *data++ = T(g*inv_scale); *data++ = T(a*inv_scale); } }  break;
137        case(GL_BGR):               { for(unsigned int i=0;i<num;++i) { float b = float(*data)*scale; float g = float(*(data+1))*scale; float r = float(*(data+2))*scale; operation.rgb(r,g,b); *data++ = T(b*inv_scale); *data++ = T(g*inv_scale); *data++ = T(r*inv_scale); } }  break;
138        case(GL_BGRA):              { for(unsigned int i=0;i<num;++i) { float b = float(*data)*scale; float g = float(*(data+1))*scale; float r = float(*(data+2))*scale; float a = float(*(data+3))*scale; operation.rgba(r,g,b,a); *data++ = T(g*inv_scale); *data++ = T(b*inv_scale); *data++ = T(r*inv_scale); *data++ = T(a*inv_scale); } }  break;
139    }
140}
141
142template <class M>   
143void modifyRow(unsigned int num, GLenum pixelFormat, GLenum dataType, unsigned char* data, const M& operation)
144{
145    switch(dataType)
146    {
147        case(GL_BYTE):              _modifyRow(num,pixelFormat, (char*)data,            1.0f/128.0f,        operation); break;
148        case(GL_UNSIGNED_BYTE):     _modifyRow(num,pixelFormat, (unsigned char*)data,   1.0f/255.0f,        operation); break;
149        case(GL_SHORT):             _modifyRow(num,pixelFormat, (short*) data,          1.0f/32768.0f,      operation); break;
150        case(GL_UNSIGNED_SHORT):    _modifyRow(num,pixelFormat, (unsigned short*)data,  1.0f/65535.0f,      operation); break;
151        case(GL_INT):               _modifyRow(num,pixelFormat, (int*) data,            1.0f/2147483648.0f, operation); break;
152        case(GL_UNSIGNED_INT):      _modifyRow(num,pixelFormat, (unsigned int*) data,   1.0f/4294967295.0f, operation); break;
153        case(GL_FLOAT):             _modifyRow(num,pixelFormat, (float*) data,          1.0f,               operation); break;
154    }
155}
156
157template <class M>   
158void modifyImage(osg::Image* image, const M& operation)
159{
160    if (!image) return;
161   
162    for(int r=0;r<image->r();++r)
163    {
164        for(int t=0;t<image->t();++t)
165        {
166            modifyRow(image->s(), image->getPixelFormat(), image->getDataType(), image->data(0,t,r), operation);
167        }
168    }
169}
170
171struct PassThroughTransformFunction
172{
173    unsigned char operator() (unsigned char c) const { return c; }
174};
175
176
177struct ProcessRow
178{
179    virtual ~ProcessRow() {}
180
181    virtual void operator() (unsigned int num,
182                    GLenum source_pixelFormat, unsigned char* source,
183                    GLenum dest_pixelFormat, unsigned char* dest) const 
184    {
185        switch(source_pixelFormat)
186        {
187        case(GL_LUMINANCE):
188        case(GL_ALPHA):
189            switch(dest_pixelFormat)
190            {
191            case(GL_LUMINANCE):
192            case(GL_ALPHA): A_to_A(num, source, dest); break;
193            case(GL_LUMINANCE_ALPHA): A_to_LA(num, source, dest); break;
194            case(GL_RGB): A_to_RGB(num, source, dest); break;
195            case(GL_RGBA): A_to_RGBA(num, source, dest); break;
196            }
197            break;
198        case(GL_LUMINANCE_ALPHA):
199            switch(dest_pixelFormat)
200            {
201            case(GL_LUMINANCE):
202            case(GL_ALPHA): LA_to_A(num, source, dest); break;
203            case(GL_LUMINANCE_ALPHA): LA_to_LA(num, source, dest); break;
204            case(GL_RGB): LA_to_RGB(num, source, dest); break;
205            case(GL_RGBA): LA_to_RGBA(num, source, dest); break;
206            }
207            break;
208        case(GL_RGB):
209            switch(dest_pixelFormat)
210            {
211            case(GL_LUMINANCE):
212            case(GL_ALPHA): RGB_to_A(num, source, dest); break;
213            case(GL_LUMINANCE_ALPHA): RGB_to_LA(num, source, dest); break;
214            case(GL_RGB): RGB_to_RGB(num, source, dest); break;
215            case(GL_RGBA): RGB_to_RGBA(num, source, dest); break;
216            }
217            break;
218        case(GL_RGBA):
219            switch(dest_pixelFormat)
220            {
221            case(GL_LUMINANCE):
222            case(GL_ALPHA): RGBA_to_A(num, source, dest); break;
223            case(GL_LUMINANCE_ALPHA): RGBA_to_LA(num, source, dest); break;
224            case(GL_RGB): RGBA_to_RGB(num, source, dest); break;
225            case(GL_RGBA): RGBA_to_RGBA(num, source, dest); break;
226            }
227            break;
228        }
229    }
230
231    ///////////////////////////////////////////////////////////////////////////////
232    // alpha sources..   
233    virtual void A_to_A(unsigned int num, unsigned char* source, unsigned char* dest) const
234    {
235        for(unsigned int i=0;i<num;++i)
236        {
237            *dest++ = *source++;
238        }
239    }
240
241    virtual void A_to_LA(unsigned int num, unsigned char* source, unsigned char* dest) const
242    {
243        for(unsigned int i=0;i<num;++i)
244        {
245            *dest++ = *source;
246            *dest++ = *source++;
247        }
248    }
249                   
250    virtual void A_to_RGB(unsigned int num, unsigned char* source, unsigned char* dest) const
251    {
252        for(unsigned int i=0;i<num;++i)
253        {
254            *dest++ = *source;
255            *dest++ = *source;
256            *dest++ = *source++;
257        }
258    }
259
260    virtual void A_to_RGBA(unsigned int num, unsigned char* source, unsigned char* dest) const
261    {
262        for(unsigned int i=0;i<num;++i)
263        {
264            *dest++ = *source;
265            *dest++ = *source;
266            *dest++ = *source;
267            *dest++ = *source++;
268        }
269    }
270
271    ///////////////////////////////////////////////////////////////////////////////
272    // alpha luminance sources..   
273    virtual void LA_to_A(unsigned int num, unsigned char* source, unsigned char* dest) const
274    {
275        for(unsigned int i=0;i<num;++i)
276        {
277            ++source;
278            *dest++ = *source++;
279        }
280    }
281
282    virtual void LA_to_LA(unsigned int num, unsigned char* source, unsigned char* dest) const
283    {
284        for(unsigned int i=0;i<num;++i)
285        {
286            *dest++ = *source++;
287            *dest++ = *source++;
288        }
289    }
290                   
291    virtual void LA_to_RGB(unsigned int num, unsigned char* source, unsigned char* dest) const
292    {
293        for(unsigned int i=0;i<num;++i)
294        {
295            *dest++ = *source;
296            *dest++ = *source;
297            *dest++ = *source;
298            source+=2;
299        }
300    }
301
302    virtual void LA_to_RGBA(unsigned int num, unsigned char* source, unsigned char* dest) const
303    {
304        for(unsigned int i=0;i<num;++i)
305        {
306            *dest++ = *source;
307            *dest++ = *source;
308            *dest++ = *source++;
309            *dest++ = *source++;
310        }
311    }
312
313    ///////////////////////////////////////////////////////////////////////////////
314    // RGB sources..   
315    virtual void RGB_to_A(unsigned int num, unsigned char* source, unsigned char* dest) const
316    {
317        for(unsigned int i=0;i<num;++i)
318        {
319            unsigned char val = *source;
320            *dest++ = val;
321            source += 3;
322        }
323    }
324
325    virtual void RGB_to_LA(unsigned int num, unsigned char* source, unsigned char* dest) const
326    {
327        for(unsigned int i=0;i<num;++i)
328        {
329            unsigned char val = *source;
330            *dest++ = val;
331            *dest++ = val;
332            source += 3;
333        }
334    }
335                   
336    virtual void RGB_to_RGB(unsigned int num, unsigned char* source, unsigned char* dest) const
337    {
338        for(unsigned int i=0;i<num;++i)
339        {
340            *dest++ = *source++;
341            *dest++ = *source++;
342            *dest++ = *source++;
343        }
344    }
345
346    virtual void RGB_to_RGBA(unsigned int num, unsigned char* source, unsigned char* dest) const
347    {
348        for(unsigned int i=0;i<num;++i)
349        {
350            unsigned char val = *source;
351            *dest++ = *source++;
352            *dest++ = *source++;
353            *dest++ = *source++;
354            *dest++ = val;
355        }
356    }
357
358    ///////////////////////////////////////////////////////////////////////////////
359    // RGBA sources..   
360    virtual void RGBA_to_A(unsigned int num, unsigned char* source, unsigned char* dest) const
361    {
362        for(unsigned int i=0;i<num;++i)
363        {
364            source += 3;
365            *dest++ = *source++;
366        }
367    }
368
369    virtual void RGBA_to_LA(unsigned int num, unsigned char* source, unsigned char* dest) const
370    {
371        for(unsigned int i=0;i<num;++i)
372        {
373            unsigned char val = *source;
374            source += 3;
375            *dest++ = val;
376            *dest++ = *source++;
377        }
378    }
379                   
380    virtual void RGBA_to_RGB(unsigned int num, unsigned char* source, unsigned char* dest) const
381    {
382        for(unsigned int i=0;i<num;++i)
383        {
384            *dest++ = *source++;
385            *dest++ = *source++;
386            *dest++ = *source++;
387            ++source;
388        }
389    }
390
391    virtual void RGBA_to_RGBA(unsigned int num, unsigned char* source, unsigned char* dest) const
392    {
393        for(unsigned int i=0;i<num;++i)
394        {
395            *dest++ = *source++;
396            *dest++ = *source++;
397            *dest++ = *source++;
398            *dest++ = *source++;
399        }
400    }
401};
402
403
404void clampToNearestValidPowerOfTwo(int& sizeX, int& sizeY, int& sizeZ, int s_maximumTextureSize, int t_maximumTextureSize, int r_maximumTextureSize)
405{
406    // compute nearest powers of two for each axis.
407    int s_nearestPowerOfTwo = 1;
408    while(s_nearestPowerOfTwo<sizeX && s_nearestPowerOfTwo<s_maximumTextureSize) s_nearestPowerOfTwo*=2;
409
410    int t_nearestPowerOfTwo = 1;
411    while(t_nearestPowerOfTwo<sizeY && t_nearestPowerOfTwo<t_maximumTextureSize) t_nearestPowerOfTwo*=2;
412
413    int r_nearestPowerOfTwo = 1;
414    while(r_nearestPowerOfTwo<sizeZ && r_nearestPowerOfTwo<r_maximumTextureSize) r_nearestPowerOfTwo*=2;
415
416    sizeX = s_nearestPowerOfTwo;
417    sizeY = t_nearestPowerOfTwo;
418    sizeZ = r_nearestPowerOfTwo;
419}
420
421osg::Image* createTexture3D(ImageList& imageList, ProcessRow& processRow,
422            unsigned int numComponentsDesired,
423            int s_maximumTextureSize,
424            int t_maximumTextureSize,
425            int r_maximumTextureSize,
426            bool resizeToPowerOfTwo)
427{
428    int max_s = 0;
429    int max_t = 0;
430    unsigned int max_components = 0;
431    int total_r = 0;
432    ImageList::iterator itr;
433    for(itr=imageList.begin();
434        itr!=imageList.end();
435        ++itr)
436    {
437        osg::Image* image = itr->get();
438        GLenum pixelFormat = image->getPixelFormat();
439        if (pixelFormat==GL_ALPHA ||
440            pixelFormat==GL_INTENSITY ||
441            pixelFormat==GL_LUMINANCE ||
442            pixelFormat==GL_LUMINANCE_ALPHA ||
443            pixelFormat==GL_RGB ||
444            pixelFormat==GL_RGBA)
445        {
446            max_s = osg::maximum(image->s(), max_s);
447            max_t = osg::maximum(image->t(), max_t);
448            max_components = osg::maximum(osg::Image::computeNumComponents(pixelFormat), max_components);
449            total_r += image->r();
450        }
451        else
452        {
453            osg::notify(osg::NOTICE)<<"Image "<<image->getFileName()<<" has unsuitable pixel format"<< std::hex<< pixelFormat << std::dec << std::endl;
454        }
455    }
456   
457    if (numComponentsDesired!=0) max_components = numComponentsDesired;
458   
459    GLenum desiredPixelFormat = 0;
460    switch(max_components)
461    {
462    case(1):
463        osg::notify(osg::NOTICE)<<"desiredPixelFormat = GL_LUMINANCE" << std::endl;
464        desiredPixelFormat = GL_LUMINANCE;
465        break;
466    case(2):
467        osg::notify(osg::NOTICE)<<"desiredPixelFormat = GL_LUMINANCE_ALPHA" << std::endl;
468        desiredPixelFormat = GL_LUMINANCE_ALPHA;
469        break;
470    case(3):
471        osg::notify(osg::NOTICE)<<"desiredPixelFormat = GL_RGB" << std::endl;
472        desiredPixelFormat = GL_RGB;
473        break;
474    case(4):
475        osg::notify(osg::NOTICE)<<"desiredPixelFormat = GL_RGBA" << std::endl;
476        desiredPixelFormat = GL_RGBA;
477        break;
478    }   
479    if (desiredPixelFormat==0) return 0;
480   
481    // compute nearest powers of two for each axis.
482   
483    int s_nearestPowerOfTwo = 1;
484    int t_nearestPowerOfTwo = 1;
485    int r_nearestPowerOfTwo = 1;
486
487    if (resizeToPowerOfTwo)
488    {
489        while(s_nearestPowerOfTwo<max_s && s_nearestPowerOfTwo<s_maximumTextureSize) s_nearestPowerOfTwo*=2;
490        while(t_nearestPowerOfTwo<max_t && t_nearestPowerOfTwo<t_maximumTextureSize) t_nearestPowerOfTwo*=2;
491        while(r_nearestPowerOfTwo<total_r && r_nearestPowerOfTwo<r_maximumTextureSize) r_nearestPowerOfTwo*=2;
492
493        osg::notify(osg::NOTICE)<<"max image width = "<<max_s<<"  nearest power of two = "<<s_nearestPowerOfTwo<<std::endl;
494        osg::notify(osg::NOTICE)<<"max image height = "<<max_t<<"  nearest power of two = "<<t_nearestPowerOfTwo<<std::endl;
495        osg::notify(osg::NOTICE)<<"max image depth = "<<total_r<<"  nearest power of two = "<<r_nearestPowerOfTwo<<std::endl;
496    }
497    else
498    {
499        s_nearestPowerOfTwo = max_s;
500        t_nearestPowerOfTwo = max_t;
501        r_nearestPowerOfTwo = total_r;
502    }
503   
504    // now allocate the 3d texture;
505    osg::ref_ptr<osg::Image> image_3d = new osg::Image;
506    image_3d->allocateImage(s_nearestPowerOfTwo,t_nearestPowerOfTwo,r_nearestPowerOfTwo,
507                            desiredPixelFormat,GL_UNSIGNED_BYTE);
508       
509
510    unsigned int r_offset = (total_r<r_nearestPowerOfTwo) ? r_nearestPowerOfTwo/2 - total_r/2 : 0;
511
512    int curr_dest_r = r_offset;
513
514    // copy across the values from the source images into the image_3d.
515    for(itr=imageList.begin();
516        itr!=imageList.end();
517        ++itr)
518    {
519        osg::Image* image = itr->get();
520        GLenum pixelFormat = image->getPixelFormat();
521        if (pixelFormat==GL_ALPHA ||
522            pixelFormat==GL_LUMINANCE ||
523            pixelFormat==GL_INTENSITY ||
524            pixelFormat==GL_LUMINANCE_ALPHA ||
525            pixelFormat==GL_RGB ||
526            pixelFormat==GL_RGBA)
527        {
528       
529            int num_r = osg::minimum(image->r(), (image_3d->r() - curr_dest_r));
530            int num_t = osg::minimum(image->t(), image_3d->t());
531            int num_s = osg::minimum(image->s(), image_3d->s());
532       
533            unsigned int s_offset_dest = (image->s()<s_nearestPowerOfTwo) ? s_nearestPowerOfTwo/2 - image->s()/2 : 0;
534            unsigned int t_offset_dest = (image->t()<t_nearestPowerOfTwo) ? t_nearestPowerOfTwo/2 - image->t()/2 : 0;
535
536            for(int r=0;r<num_r;++r, ++curr_dest_r)
537            {
538                for(int t=0;t<num_t;++t)
539                {
540                    unsigned char* dest = image_3d->data(s_offset_dest,t+t_offset_dest,curr_dest_r);
541                    unsigned char* source = image->data(0,t,r);
542
543                    processRow(num_s, image->getPixelFormat(), source, image_3d->getPixelFormat(), dest);
544                }
545            }
546        }
547    }
548    return image_3d.release();
549}
550
551
552osg::Image* createNormalMapTexture(osg::Image* image_3d)
553{
554    unsigned int sourcePixelIncrement = 1;
555    unsigned int alphaOffset = 0;
556    switch(image_3d->getPixelFormat())
557    {
558    case(GL_ALPHA):
559    case(GL_LUMINANCE):
560        sourcePixelIncrement = 1;
561        alphaOffset = 0;
562        break;
563    case(GL_LUMINANCE_ALPHA):
564        sourcePixelIncrement = 2;
565        alphaOffset = 1;
566        break;
567    case(GL_RGB):
568        sourcePixelIncrement = 3;
569        alphaOffset = 0;
570        break;
571    case(GL_RGBA):
572        sourcePixelIncrement = 4;
573        alphaOffset = 3;
574        break;
575    default:
576        osg::notify(osg::NOTICE)<<"Source pixel format not support for normal map generation."<<std::endl;
577        return 0;
578    }
579   
580    osg::ref_ptr<osg::Image> normalmap_3d = new osg::Image;
581    normalmap_3d->allocateImage(image_3d->s(),image_3d->t(),image_3d->r(),
582                            GL_RGBA,GL_UNSIGNED_BYTE);
583
584    if (osg::getCpuByteOrder()==osg::LittleEndian) alphaOffset = sourcePixelIncrement-alphaOffset-1;
585
586    for(int r=1;r<image_3d->r()-1;++r)
587    {
588        for(int t=1;t<image_3d->t()-1;++t)
589        {
590            unsigned char* ptr = image_3d->data(1,t,r)+alphaOffset;
591            unsigned char* left = image_3d->data(0,t,r)+alphaOffset;
592            unsigned char* right = image_3d->data(2,t,r)+alphaOffset;
593            unsigned char* above = image_3d->data(1,t+1,r)+alphaOffset;
594            unsigned char* below = image_3d->data(1,t-1,r)+alphaOffset;
595            unsigned char* in = image_3d->data(1,t,r+1)+alphaOffset;
596            unsigned char* out = image_3d->data(1,t,r-1)+alphaOffset;
597
598            unsigned char* destination = (unsigned char*) normalmap_3d->data(1,t,r);
599
600            for(int s=1;s<image_3d->s()-1;++s)
601            {
602
603                osg::Vec3 grad((float)(*left)-(float)(*right),
604                               (float)(*below)-(float)(*above),
605                               (float)(*out) -(float)(*in));
606
607                grad.normalize();
608
609                if (grad.x()==0.0f && grad.y()==0.0f && grad.z()==0.0f)
610                {
611                    grad.set(128.0f,128.0f,128.0f);
612                }
613                else
614                {
615                    grad.x() = osg::clampBetween((grad.x()+1.0f)*128.0f,0.0f,255.0f);
616                    grad.y() = osg::clampBetween((grad.y()+1.0f)*128.0f,0.0f,255.0f);
617                    grad.z() = osg::clampBetween((grad.z()+1.0f)*128.0f,0.0f,255.0f);
618                }
619
620                *(destination++) = (unsigned char)(grad.x()); // scale and bias X.
621                *(destination++) = (unsigned char)(grad.y()); // scale and bias Y.
622                *(destination++) = (unsigned char)(grad.z()); // scale and bias Z.
623
624                *destination++ = *ptr;
625
626                ptr += sourcePixelIncrement;
627                left += sourcePixelIncrement;
628                right += sourcePixelIncrement;
629                above += sourcePixelIncrement;
630                below += sourcePixelIncrement;
631                in += sourcePixelIncrement;
632                out += sourcePixelIncrement;
633            }
634        }
635    }
636   
637    return normalmap_3d.release();
638}
639
640
641
642osg::Node* createCube(float size,float alpha, unsigned int numSlices, float sliceEnd=1.0f)
643{
644
645    // set up the Geometry.
646    osg::Geometry* geom = new osg::Geometry;
647
648    float halfSize = size*0.5f;
649    float y = halfSize;
650    float dy =-size/(float)(numSlices-1)*sliceEnd;
651
652    //y = -halfSize;
653    //dy *= 0.5;
654
655    osg::Vec3Array* coords = new osg::Vec3Array(4*numSlices);
656    geom->setVertexArray(coords);
657    for(unsigned int i=0;i<numSlices;++i, y+=dy)
658    {
659        (*coords)[i*4+0].set(-halfSize,y,halfSize);
660        (*coords)[i*4+1].set(-halfSize,y,-halfSize);
661        (*coords)[i*4+2].set(halfSize,y,-halfSize);
662        (*coords)[i*4+3].set(halfSize,y,halfSize);
663    }
664   
665    osg::Vec3Array* normals = new osg::Vec3Array(1);
666    (*normals)[0].set(0.0f,-1.0f,0.0f);
667    geom->setNormalArray(normals);
668    geom->setNormalBinding(osg::Geometry::BIND_OVERALL);
669
670    osg::Vec4Array* colors = new osg::Vec4Array(1);
671    (*colors)[0].set(1.0f,1.0f,1.0f,alpha);
672    geom->setColorArray(colors);
673    geom->setColorBinding(osg::Geometry::BIND_OVERALL);
674
675    geom->addPrimitiveSet(new osg::DrawArrays(osg::PrimitiveSet::QUADS,0,coords->size()));
676
677    osg::Billboard* billboard = new osg::Billboard;
678    billboard->setMode(osg::Billboard::POINT_ROT_WORLD);
679    billboard->addDrawable(geom);
680    billboard->setPosition(0,osg::Vec3(0.0f,0.0f,0.0f));
681   
682    return billboard;
683}
684
685class FollowMouseCallback : public osgGA::GUIEventHandler, public osg::StateSet::Callback
686{
687    public:
688   
689        FollowMouseCallback(bool shader = false):
690            _shader(shader)
691        {
692            _updateTransparency = false;
693            _updateAlphaCutOff = false;
694            _updateSampleDensity = false;
695        }
696
697        FollowMouseCallback(const FollowMouseCallback&,const osg::CopyOp&) {}
698
699        META_Object(osg,FollowMouseCallback);
700
701        virtual void operator() (osg::StateSet* stateset, osg::NodeVisitor* nv)
702        {
703            if (nv->getVisitorType()==osg::NodeVisitor::EVENT_VISITOR)
704            {
705                osgGA::EventVisitor* ev = dynamic_cast<osgGA::EventVisitor*>(nv);
706                if (ev)
707                {
708                    osgGA::GUIActionAdapter* aa = ev->getActionAdapter();
709                    osgGA::EventQueue::Events& events = ev->getEvents();
710                    for(osgGA::EventQueue::Events::iterator itr=events.begin();
711                        itr!=events.end();
712                        ++itr)
713                    {
714                        handle(*(*itr), *aa, stateset, ev);
715                    }
716                }
717            }
718        }
719       
720        virtual bool handle(const osgGA::GUIEventAdapter& ea,osgGA::GUIActionAdapter&, osg::Object* object, osg::NodeVisitor*)
721        {
722            osg::StateSet* stateset = dynamic_cast<osg::StateSet*>(object);
723            if (!stateset) return false;
724           
725            switch(ea.getEventType())
726            {
727                case(osgGA::GUIEventAdapter::MOVE):
728                case(osgGA::GUIEventAdapter::DRAG):
729                {
730                    float v = (ea.getY()-ea.getYmin())/(ea.getYmax()-ea.getYmin());
731                    if (_shader)
732                    {
733                        osg::Uniform* uniform = 0;
734                        if (_updateTransparency && (uniform = stateset->getUniform("transparency"))) uniform->set(v);
735                        if (_updateAlphaCutOff && (uniform = stateset->getUniform("alphaCutOff"))) uniform->set(v);
736                        if (_updateSampleDensity && (uniform = stateset->getUniform("sampleDensity"))) uniform->set(powf(v,5));
737                    }
738                    else
739                    {                   
740                        if (_updateAlphaCutOff)
741                        {
742                            osg::AlphaFunc* alphaFunc = dynamic_cast<osg::AlphaFunc*>(stateset->getAttribute(osg::StateAttribute::ALPHAFUNC));
743                            if (alphaFunc)
744                            {
745                                alphaFunc->setReferenceValue(v);
746                            }
747                        }
748                       
749                        if (_updateTransparency)
750                        {
751                            osg::Material* material = dynamic_cast<osg::Material*>(stateset->getAttribute(osg::StateAttribute::MATERIAL));
752                            if (material)
753                            {
754                                material->setAlpha(osg::Material::FRONT_AND_BACK,v);
755                            }
756                        }
757                    }
758
759                    break;
760                }
761                case(osgGA::GUIEventAdapter::KEYDOWN):
762                {
763                    if (ea.getKey()=='t') _updateTransparency = true;
764                    if (ea.getKey()=='a') _updateAlphaCutOff = true;
765                    if (ea.getKey()=='d') _updateSampleDensity = true;
766                    break;
767                }
768                case(osgGA::GUIEventAdapter::KEYUP):
769                {
770                    if (ea.getKey()=='t') _updateTransparency = false;
771                    if (ea.getKey()=='a') _updateAlphaCutOff = false;
772                    if (ea.getKey()=='d') _updateSampleDensity = false;
773                    break;
774                }
775                default:
776                    break;
777            }
778            return false;
779        }
780       
781        bool _shader;
782        bool _updateTransparency;
783        bool _updateAlphaCutOff;
784        bool _updateSampleDensity;
785
786};
787
788osg::Node* createShaderModel(osg::ref_ptr<osg::Image>& image_3d, osg::ref_ptr<osg::Image>& /*normalmap_3d*/,
789                       osg::Texture::InternalFormatMode internalFormatMode,
790                       float xSize, float ySize, float zSize,
791                       float /*xMultiplier*/, float /*yMultiplier*/, float /*zMultiplier*/,
792                       unsigned int /*numSlices*/=500, float /*sliceEnd*/=1.0f, float alphaFuncValue=0.02f, bool maximumIntensityProjection = false)
793{
794
795    osg::Group* root = new osg::Group;
796   
797    osg::Geode* geode = new osg::Geode;
798    root->addChild(geode);
799   
800    osg::StateSet* stateset = geode->getOrCreateStateSet();
801   
802    stateset->setEventCallback(new FollowMouseCallback(true));
803   
804    stateset->setMode(GL_ALPHA_TEST,osg::StateAttribute::ON);
805
806    // set up the 3d texture itself,
807    // note, well set the filtering up so that mip mapping is disabled,
808    // gluBuild3DMipsmaps doesn't do a very good job of handled the
809    // imbalanced dimensions of the 256x256x4 texture.
810    osg::Texture3D* texture3D = new osg::Texture3D;
811    texture3D->setResizeNonPowerOfTwoHint(false);
812    texture3D->setFilter(osg::Texture3D::MIN_FILTER,osg::Texture3D::LINEAR);
813    texture3D->setFilter(osg::Texture3D::MAG_FILTER,osg::Texture3D::LINEAR);
814    texture3D->setWrap(osg::Texture3D::WRAP_R,osg::Texture3D::CLAMP);
815    texture3D->setWrap(osg::Texture3D::WRAP_S,osg::Texture3D::CLAMP);
816    texture3D->setWrap(osg::Texture3D::WRAP_T,osg::Texture3D::CLAMP);
817    if (image_3d->getPixelFormat()==GL_ALPHA ||
818        image_3d->getPixelFormat()==GL_LUMINANCE)
819    {
820        texture3D->setInternalFormatMode(osg::Texture3D::USE_USER_DEFINED_FORMAT);
821        texture3D->setInternalFormat(GL_INTENSITY);
822    }
823    else
824    {
825        texture3D->setInternalFormatMode(internalFormatMode);
826    }
827
828    texture3D->setImage(image_3d.get());
829
830    stateset->setTextureAttributeAndModes(0,texture3D,osg::StateAttribute::ON);
831
832    osg::Program* program = new osg::Program;
833    stateset->setAttribute(program);
834
835    // get shaders from source
836    std::string vertexShaderFile = osgDB::findDataFile("volume.vert");
837    if (!vertexShaderFile.empty())
838    {
839        program->addShader(osg::Shader::readShaderFile(osg::Shader::VERTEX, vertexShaderFile));
840    }
841    else
842    {
843        char vertexShaderSource[] =
844            "#version 110\n"
845            "varying vec4 cameraPos;\n"
846            "varying vec4 vertexPos;\n"
847            "varying mat4 texgen;\n"
848            "\n"
849            "void main(void)\n"
850            "{\n"
851            "        gl_Position = ftransform();\n"
852            "\n"
853            "        cameraPos = gl_ModelViewMatrixInverse*vec4(0,0,0,1);\n"
854            "        vertexPos = gl_Vertex;\n"
855            "\n"
856            "        texgen = mat4(gl_ObjectPlaneS[0], \n"
857            "                      gl_ObjectPlaneT[0],\n"
858            "                      gl_ObjectPlaneR[0],\n"
859            "                      gl_ObjectPlaneQ[0]);\n"
860            "}\n";
861
862        osg::Shader* vertex_shader = new osg::Shader(osg::Shader::VERTEX, vertexShaderSource);
863        program->addShader(vertex_shader);
864
865    }
866   
867    std::string fragmentShaderFile = osgDB::findDataFile("volume.frag");
868    if (!fragmentShaderFile.empty())
869    {
870        program->addShader(osg::Shader::readShaderFile(osg::Shader::FRAGMENT, fragmentShaderFile));
871    }
872    else
873    {
874        //////////////////////////////////////////////////////////////////
875        // fragment shader
876        //
877        char fragmentShaderSource[] =
878            "uniform sampler3D baseTexture;\n"
879            "uniform float sampleDensity;\n"
880            "uniform float transparency;\n"
881            "uniform float alphaCutOff;\n"
882            "\n"
883            "varying vec4 cameraPos;\n"
884            "varying vec4 vertexPos;\n"
885            "varying mat4 texgen;\n"
886            "\n"
887            "void main(void)\n"
888            "{ \n"
889            "    vec3 t0 = (texgen * vertexPos).xyz;\n"
890            "    vec3 te = (texgen * cameraPos).xyz;\n"
891            "\n"
892            "    if (te.x>=0.0 && te.x<=1.0 &&\n"
893            "        te.y>=0.0 && te.y<=1.0 &&\n"
894            "        te.z>=0.0 && te.z<=1.0)\n"
895            "    {\n"
896            "        // do nothing... te inside volume\n"
897            "    }\n"
898            "    else\n"
899            "    {\n"
900            "        if (te.x<0.0)\n"
901            "        {\n"
902            "            float r = -te.x / (t0.x-te.x);\n"
903            "            te = te + (t0-te)*r;\n"
904            "        }\n"
905            "\n"
906            "        if (te.x>1.0)\n"
907            "        {\n"
908            "            float r = (1.0-te.x) / (t0.x-te.x);\n"
909            "            te = te + (t0-te)*r;\n"
910            "        }\n"
911            "\n"
912            "        if (te.y<0.0)\n"
913            "        {\n"
914            "            float r = -te.y / (t0.y-te.y);\n"
915            "            te = te + (t0-te)*r;\n"
916            "        }\n"
917            "\n"
918            "        if (te.y>1.0)\n"
919            "        {\n"
920            "            float r = (1.0-te.y) / (t0.y-te.y);\n"
921            "            te = te + (t0-te)*r;\n"
922            "        }\n"
923            "\n"
924            "        if (te.z<0.0)\n"
925            "        {\n"
926            "            float r = -te.z / (t0.z-te.z);\n"
927            "            te = te + (t0-te)*r;\n"
928            "        }\n"
929            "\n"
930            "        if (te.z>1.0)\n"
931            "        {\n"
932            "            float r = (1.0-te.z) / (t0.z-te.z);\n"
933            "            te = te + (t0-te)*r;\n"
934            "        }\n"
935            "    }\n"
936            "\n"
937            "    const float max_iteratrions = 256.0;\n"
938            "    float num_iterations = length(te-t0)/sampleDensity;\n"
939            "    if (num_iterations>max_iteratrions) \n"
940            "    {\n"
941            "        num_iterations = max_iteratrions;\n"
942            "    }\n"
943            "\n"
944            "    vec3 deltaTexCoord=(te-t0)/float(num_iterations-1.0);\n"
945            "    vec3 texcoord = t0;\n"
946            "\n"
947            "    vec4 fragColor = vec4(0.0, 0.0, 0.0, 0.0); \n"
948            "    while(num_iterations>0.0)\n"
949            "    {\n"
950            "        vec4 color = texture3D( baseTexture, texcoord);\n"
951            "        float r = color[3]*transparency;\n"
952            "        if (r>alphaCutOff)\n"
953            "        {\n"
954            "            fragColor.xyz = fragColor.xyz*(1.0-r)+color.xyz*r;\n"
955            "            fragColor.w += r;\n"
956            "        }\n"
957            "        texcoord += deltaTexCoord; \n"
958            "\n"
959            "        --num_iterations;\n"
960            "    }\n"
961            "\n"
962            "    if (fragColor.w>1.0) fragColor.w = 1.0; \n"
963            "    if (fragColor.w==0.0) discard;\n"
964            "    gl_FragColor = fragColor;\n"
965            "}\n";
966
967        osg::Shader* fragment_shader = new osg::Shader(osg::Shader::FRAGMENT, fragmentShaderSource);
968        program->addShader(fragment_shader);
969    }
970
971    osg::Uniform* baseTextureSampler = new osg::Uniform("baseTexture",0);
972    stateset->addUniform(baseTextureSampler);
973
974    osg::Uniform* sampleDensity = new osg::Uniform("sampleDensity", 0.01f);
975    stateset->addUniform(sampleDensity);
976
977    osg::Uniform* transpancy = new osg::Uniform("transparency",0.5f);
978    stateset->addUniform(transpancy);
979
980    osg::Uniform* alphaCutOff = new osg::Uniform("alphaCutOff",alphaFuncValue);
981    stateset->addUniform(alphaCutOff);
982
983    stateset->setMode(GL_CULL_FACE, osg::StateAttribute::ON);
984
985    osg::TexGen* texgen = new osg::TexGen;
986    texgen->setMode(osg::TexGen::OBJECT_LINEAR);
987    texgen->setPlane(osg::TexGen::S, osg::Plane(1.0f/xSize,0.0f,0.0f,0.0f));
988    texgen->setPlane(osg::TexGen::T, osg::Plane(0.0f,1.0f/ySize,0.0f,0.0f));
989    texgen->setPlane(osg::TexGen::R, osg::Plane(0.0f,0.0f,1.0f/zSize,0.0f));
990    texgen->setPlane(osg::TexGen::Q, osg::Plane(0.0f,0.0f,0.0f,1.0f));
991   
992    stateset->setTextureAttributeAndModes(0, texgen, osg::StateAttribute::ON);
993
994    {
995        osg::Geometry* geom = new osg::Geometry;
996
997        osg::Vec3Array* coords = new osg::Vec3Array(8);
998        (*coords)[0].set(0,0,0);
999        (*coords)[1].set(xSize,0,0);
1000        (*coords)[2].set(xSize,ySize,0);
1001        (*coords)[3].set(0,ySize,0);
1002        (*coords)[4].set(0,0,zSize);
1003        (*coords)[5].set(xSize,0,zSize);
1004        (*coords)[6].set(ySize,ySize,zSize);
1005        (*coords)[7].set(0,ySize,zSize);
1006        geom->setVertexArray(coords);
1007
1008        osg::Vec4Array* colours = new osg::Vec4Array(1);
1009        (*colours)[0].set(1.0f,1.0f,1.0,1.0f);
1010        geom->setColorArray(colours);
1011        geom->setColorBinding(osg::Geometry::BIND_OVERALL);
1012
1013        osg::DrawElementsUShort* drawElements = new osg::DrawElementsUShort(GL_QUADS);
1014        // bottom
1015        drawElements->push_back(0);
1016        drawElements->push_back(1);
1017        drawElements->push_back(2);
1018        drawElements->push_back(3);
1019       
1020        // bottom
1021        drawElements->push_back(3);
1022        drawElements->push_back(2);
1023        drawElements->push_back(6);
1024        drawElements->push_back(7);
1025
1026        // left
1027        drawElements->push_back(0);
1028        drawElements->push_back(3);
1029        drawElements->push_back(7);
1030        drawElements->push_back(4);
1031
1032        // right
1033        drawElements->push_back(5);
1034        drawElements->push_back(6);
1035        drawElements->push_back(2);
1036        drawElements->push_back(1);
1037
1038        // front
1039        drawElements->push_back(1);
1040        drawElements->push_back(0);
1041        drawElements->push_back(4);
1042        drawElements->push_back(5);
1043
1044        // top
1045        drawElements->push_back(7);
1046        drawElements->push_back(6);
1047        drawElements->push_back(5);
1048        drawElements->push_back(4);
1049
1050        geom->addPrimitiveSet(drawElements);
1051
1052        geode->addDrawable(geom);
1053
1054    }
1055    return root;
1056}
1057
1058osg::Node* createModel(osg::ref_ptr<osg::Image>& image_3d,
1059                       osg::ref_ptr<osg::Image>& normalmap_3d,
1060                       osg::Texture::InternalFormatMode internalFormatMode,
1061                       float xSize, float ySize, float zSize,
1062                       float xMultiplier, float yMultiplier, float zMultiplier,
1063                       unsigned int numSlices=500, float sliceEnd=1.0f, float alphaFuncValue=0.02f, bool maximumIntensityProjection = false)
1064{
1065    bool two_pass = normalmap_3d.valid() && (image_3d->getPixelFormat()==GL_RGB || image_3d->getPixelFormat()==GL_RGBA);
1066
1067    osg::BoundingBox bb(-xSize*0.5f,-ySize*0.5f,-zSize*0.5f,xSize*0.5f,ySize*0.5f,zSize*0.5f);
1068
1069    float maxAxis = xSize;
1070    if (ySize > maxAxis) maxAxis = ySize;
1071    if (zSize > maxAxis) maxAxis = zSize;
1072
1073    osg::Group* group = new osg::Group;
1074   
1075    osg::TexGenNode* texgenNode_0 = new osg::TexGenNode;
1076    texgenNode_0->setTextureUnit(0);
1077    texgenNode_0->getTexGen()->setMode(osg::TexGen::EYE_LINEAR);
1078    texgenNode_0->getTexGen()->setPlane(osg::TexGen::S, osg::Plane(xMultiplier/xSize,0.0f,0.0f,0.5f));
1079    texgenNode_0->getTexGen()->setPlane(osg::TexGen::T, osg::Plane(0.0f,yMultiplier/ySize,0.0f,0.5f));
1080    texgenNode_0->getTexGen()->setPlane(osg::TexGen::R, osg::Plane(0.0f,0.0f,zMultiplier/zSize,0.5f));
1081   
1082    if (two_pass)
1083    {
1084        osg::TexGenNode* texgenNode_1 = new osg::TexGenNode;
1085        texgenNode_1->setTextureUnit(1);
1086        texgenNode_1->getTexGen()->setMode(osg::TexGen::EYE_LINEAR);
1087        texgenNode_1->getTexGen()->setPlane(osg::TexGen::S, texgenNode_0->getTexGen()->getPlane(osg::TexGen::S));
1088        texgenNode_1->getTexGen()->setPlane(osg::TexGen::T, texgenNode_0->getTexGen()->getPlane(osg::TexGen::T));
1089        texgenNode_1->getTexGen()->setPlane(osg::TexGen::R, texgenNode_0->getTexGen()->getPlane(osg::TexGen::R));
1090
1091        texgenNode_1->addChild(texgenNode_0);
1092
1093        group->addChild(texgenNode_1);
1094    }
1095    else
1096    { 
1097        group->addChild(texgenNode_0);
1098    }
1099
1100    float cubeSize = sqrtf(xSize*xSize+ySize*ySize+zSize*zSize);
1101
1102    osg::ClipNode* clipnode = new osg::ClipNode;
1103    clipnode->addChild(createCube(cubeSize,1.0f, numSlices,sliceEnd));
1104    clipnode->createClipBox(bb);
1105
1106    {
1107        // set up the Geometry to enclose the clip volume to prevent near/far clipping from affecting billboard
1108        osg::Geometry* geom = new osg::Geometry;
1109
1110        osg::Vec3Array* coords = new osg::Vec3Array();
1111        coords->push_back(bb.corner(0));
1112        coords->push_back(bb.corner(1));
1113        coords->push_back(bb.corner(2));
1114        coords->push_back(bb.corner(3));
1115        coords->push_back(bb.corner(4));
1116        coords->push_back(bb.corner(5));
1117        coords->push_back(bb.corner(6));
1118        coords->push_back(bb.corner(7));
1119
1120        geom->setVertexArray(coords);
1121
1122        osg::Vec4Array* colors = new osg::Vec4Array(1);
1123        (*colors)[0].set(1.0f,1.0f,1.0f,1.0f);
1124        geom->setColorArray(colors);
1125        geom->setColorBinding(osg::Geometry::BIND_OVERALL);
1126
1127        geom->addPrimitiveSet(new osg::DrawArrays(osg::PrimitiveSet::POINTS,0,coords->size()));
1128
1129        osg::Geode* geode = new osg::Geode;
1130        geode->addDrawable(geom);
1131       
1132        clipnode->addChild(geode);
1133       
1134    }
1135
1136    texgenNode_0->addChild(clipnode);
1137
1138    osg::StateSet* stateset = texgenNode_0->getOrCreateStateSet();
1139
1140    stateset->setEventCallback(new FollowMouseCallback(false));
1141 
1142    stateset->setMode(GL_LIGHTING,osg::StateAttribute::ON);
1143    stateset->setMode(GL_BLEND,osg::StateAttribute::ON);
1144    stateset->setAttributeAndModes(new osg::AlphaFunc(osg::AlphaFunc::GREATER,alphaFuncValue), osg::StateAttribute::ON);
1145   
1146    osg::Material* material = new osg::Material;
1147    material->setDiffuse(osg::Material::FRONT_AND_BACK,osg::Vec4(1.0f,1.0f,1.0f,1.0f));
1148    stateset->setAttributeAndModes(material);
1149   
1150    if (maximumIntensityProjection)
1151    {
1152        stateset->setAttribute(new osg::BlendFunc(osg::BlendFunc::ONE, osg::BlendFunc::ONE));
1153        stateset->setAttribute(new osg::BlendEquation(osg::BlendEquation::RGBA_MAX));
1154    }
1155   
1156    osg::Vec3 lightDirection(1.0f,-1.0f,1.0f);
1157    lightDirection.normalize();
1158
1159    if (normalmap_3d.valid())
1160    {
1161        if (two_pass)
1162        {
1163
1164            // set up normal texture
1165            osg::Texture3D* bump_texture3D = new osg::Texture3D;
1166            bump_texture3D->setFilter(osg::Texture3D::MIN_FILTER,osg::Texture3D::LINEAR);
1167            bump_texture3D->setFilter(osg::Texture3D::MAG_FILTER,osg::Texture3D::LINEAR);
1168            bump_texture3D->setWrap(osg::Texture3D::WRAP_R,osg::Texture3D::CLAMP);
1169            bump_texture3D->setWrap(osg::Texture3D::WRAP_S,osg::Texture3D::CLAMP);
1170            bump_texture3D->setWrap(osg::Texture3D::WRAP_T,osg::Texture3D::CLAMP);
1171            bump_texture3D->setImage(normalmap_3d.get());
1172
1173            bump_texture3D->setInternalFormatMode(internalFormatMode);
1174
1175            stateset->setTextureAttributeAndModes(0,bump_texture3D,osg::StateAttribute::ON);
1176
1177            osg::TexEnvCombine* tec = new osg::TexEnvCombine;
1178            tec->setConstantColorAsLightDirection(lightDirection);
1179
1180            tec->setCombine_RGB(osg::TexEnvCombine::DOT3_RGB);
1181            tec->setSource0_RGB(osg::TexEnvCombine::CONSTANT);
1182            tec->setOperand0_RGB(osg::TexEnvCombine::SRC_COLOR);
1183            tec->setSource1_RGB(osg::TexEnvCombine::TEXTURE);
1184
1185            tec->setOperand1_RGB(osg::TexEnvCombine::SRC_COLOR);
1186
1187            tec->setCombine_Alpha(osg::TexEnvCombine::REPLACE);
1188            tec->setSource0_Alpha(osg::TexEnvCombine::PRIMARY_COLOR);
1189            tec->setOperand0_Alpha(osg::TexEnvCombine::SRC_ALPHA);
1190            tec->setSource1_Alpha(osg::TexEnvCombine::TEXTURE);
1191            tec->setOperand1_Alpha(osg::TexEnvCombine::SRC_ALPHA);
1192
1193            stateset->setTextureAttributeAndModes(0, tec, osg::StateAttribute::OVERRIDE|osg::StateAttribute::ON);
1194
1195            stateset->setTextureMode(0,GL_TEXTURE_GEN_S,osg::StateAttribute::ON);
1196            stateset->setTextureMode(0,GL_TEXTURE_GEN_T,osg::StateAttribute::ON);
1197            stateset->setTextureMode(0,GL_TEXTURE_GEN_R,osg::StateAttribute::ON);
1198
1199
1200            // set up color texture
1201            osg::Texture3D* texture3D = new osg::Texture3D;
1202            texture3D->setResizeNonPowerOfTwoHint(false);
1203            texture3D->setFilter(osg::Texture3D::MIN_FILTER,osg::Texture3D::LINEAR);
1204            texture3D->setFilter(osg::Texture3D::MAG_FILTER,osg::Texture3D::LINEAR);
1205            texture3D->setWrap(osg::Texture3D::WRAP_R,osg::Texture3D::CLAMP);
1206            texture3D->setWrap(osg::Texture3D::WRAP_S,osg::Texture3D::CLAMP);
1207            texture3D->setWrap(osg::Texture3D::WRAP_T,osg::Texture3D::CLAMP);
1208            if (image_3d->getPixelFormat()==GL_ALPHA ||
1209                image_3d->getPixelFormat()==GL_LUMINANCE)
1210            {
1211                texture3D->setInternalFormatMode(osg::Texture3D::USE_USER_DEFINED_FORMAT);
1212                texture3D->setInternalFormat(GL_INTENSITY);
1213            }
1214            else
1215            {
1216                texture3D->setInternalFormatMode(internalFormatMode);
1217            }
1218            texture3D->setImage(image_3d.get());
1219
1220            stateset->setTextureAttributeAndModes(1,texture3D,osg::StateAttribute::ON);
1221
1222            stateset->setTextureMode(1,GL_TEXTURE_GEN_S,osg::StateAttribute::ON);
1223            stateset->setTextureMode(1,GL_TEXTURE_GEN_T,osg::StateAttribute::ON);
1224            stateset->setTextureMode(1,GL_TEXTURE_GEN_R,osg::StateAttribute::ON);
1225
1226            stateset->setTextureAttributeAndModes(1,new osg::TexEnv(),osg::StateAttribute::ON);
1227
1228        }
1229        else
1230        {
1231            osg::ref_ptr<osg::Image> normalmap_3d = createNormalMapTexture(image_3d.get());
1232            osg::Texture3D* bump_texture3D = new osg::Texture3D;
1233            bump_texture3D->setResizeNonPowerOfTwoHint(false);
1234            bump_texture3D->setFilter(osg::Texture3D::MIN_FILTER,osg::Texture3D::LINEAR);
1235            bump_texture3D->setFilter(osg::Texture3D::MAG_FILTER,osg::Texture3D::LINEAR);
1236            bump_texture3D->setWrap(osg::Texture3D::WRAP_R,osg::Texture3D::CLAMP);
1237            bump_texture3D->setWrap(osg::Texture3D::WRAP_S,osg::Texture3D::CLAMP);
1238            bump_texture3D->setWrap(osg::Texture3D::WRAP_T,osg::Texture3D::CLAMP);
1239            bump_texture3D->setImage(normalmap_3d.get());
1240
1241            bump_texture3D->setInternalFormatMode(internalFormatMode);
1242
1243            stateset->setTextureAttributeAndModes(0,bump_texture3D,osg::StateAttribute::ON);
1244
1245            osg::TexEnvCombine* tec = new osg::TexEnvCombine;
1246            tec->setConstantColorAsLightDirection(lightDirection);
1247
1248            tec->setCombine_RGB(osg::TexEnvCombine::DOT3_RGB);
1249            tec->setSource0_RGB(osg::TexEnvCombine::CONSTANT);
1250            tec->setOperand0_RGB(osg::TexEnvCombine::SRC_COLOR);
1251            tec->setSource1_RGB(osg::TexEnvCombine::TEXTURE);
1252            tec->setOperand1_RGB(osg::TexEnvCombine::SRC_COLOR);
1253
1254            tec->setCombine_Alpha(osg::TexEnvCombine::MODULATE);
1255            tec->setSource0_Alpha(osg::TexEnvCombine::PRIMARY_COLOR);
1256            tec->setOperand0_Alpha(osg::TexEnvCombine::SRC_ALPHA);
1257            tec->setSource1_Alpha(osg::TexEnvCombine::TEXTURE);
1258            tec->setOperand1_Alpha(osg::TexEnvCombine::SRC_ALPHA);
1259
1260            stateset->setTextureAttributeAndModes(0, tec, osg::StateAttribute::OVERRIDE|osg::StateAttribute::ON);
1261
1262            stateset->setTextureMode(0,GL_TEXTURE_GEN_S,osg::StateAttribute::ON);
1263            stateset->setTextureMode(0,GL_TEXTURE_GEN_T,osg::StateAttribute::ON);
1264            stateset->setTextureMode(0,GL_TEXTURE_GEN_R,osg::StateAttribute::ON);
1265
1266            image_3d = normalmap_3d;
1267        }
1268    }
1269    else
1270    {     
1271        // set up the 3d texture itself,
1272        // note, well set the filtering up so that mip mapping is disabled,
1273        // gluBuild3DMipsmaps doesn't do a very good job of handled the
1274        // imbalanced dimensions of the 256x256x4 texture.
1275        osg::Texture3D* texture3D = new osg::Texture3D;
1276        texture3D->setResizeNonPowerOfTwoHint(false);
1277        texture3D->setFilter(osg::Texture3D::MIN_FILTER,osg::Texture3D::LINEAR);
1278        texture3D->setFilter(osg::Texture3D::MAG_FILTER,osg::Texture3D::LINEAR);
1279        texture3D->setWrap(osg::Texture3D::WRAP_R,osg::Texture3D::CLAMP);
1280        texture3D->setWrap(osg::Texture3D::WRAP_S,osg::Texture3D::CLAMP);
1281        texture3D->setWrap(osg::Texture3D::WRAP_T,osg::Texture3D::CLAMP);
1282        if (image_3d->getPixelFormat()==GL_ALPHA ||
1283            image_3d->getPixelFormat()==GL_LUMINANCE)
1284        {
1285            texture3D->setInternalFormatMode(osg::Texture3D::USE_USER_DEFINED_FORMAT);
1286            texture3D->setInternalFormat(GL_INTENSITY);
1287        }
1288        else
1289        {
1290            texture3D->setInternalFormatMode(internalFormatMode);
1291        }
1292
1293        texture3D->setImage(image_3d.get());
1294
1295        stateset->setTextureAttributeAndModes(0,texture3D,osg::StateAttribute::ON);
1296
1297        stateset->setTextureMode(0,GL_TEXTURE_GEN_S,osg::StateAttribute::ON);
1298        stateset->setTextureMode(0,GL_TEXTURE_GEN_T,osg::StateAttribute::ON);
1299        stateset->setTextureMode(0,GL_TEXTURE_GEN_R,osg::StateAttribute::ON);
1300
1301        stateset->setTextureAttributeAndModes(0,new osg::TexEnv(),osg::StateAttribute::ON);
1302    }
1303 
1304    return group;
1305}
1306
1307struct FindRangeOperator
1308{
1309    FindRangeOperator():
1310        _rmin(FLT_MAX),
1311        _rmax(-FLT_MAX),
1312        _gmin(FLT_MAX),
1313        _gmax(-FLT_MAX),
1314        _bmin(FLT_MAX),
1315        _bmax(-FLT_MAX),
1316        _amin(FLT_MAX),
1317        _amax(-FLT_MAX) {}
1318       
1319    mutable float _rmin, _rmax, _gmin, _gmax, _bmin, _bmax, _amin, _amax;
1320
1321    inline void luminance(float l) const { rgb(l,l,l); }
1322    inline void alpha(float a) const { _amin = osg::minimum(a,_amin); _amax = osg::maximum(a,_amax); }
1323    inline void luminance_alpha(float l,float a) const { rgb(l,l,l); alpha(a); }
1324    inline void rgb(float r,float g,float b) const { _rmin = osg::minimum(r,_rmin); _rmax = osg::maximum(r,_rmax); _gmin = osg::minimum(g,_gmin); _gmax = osg::maximum(g,_gmax); _bmin = osg::minimum(b,_bmin); _bmax = osg::maximum(b,_bmax);  }
1325    inline void rgba(float r,float g,float b,float a) const { rgb(r,g,b); alpha(a); }
1326};
1327 
1328struct ScaleOperator
1329{
1330    ScaleOperator():_scale(1.0f) {}
1331    ScaleOperator(float scale):_scale(scale) {}
1332    ScaleOperator(const ScaleOperator& so):_scale(so._scale) {}
1333   
1334    ScaleOperator& operator = (const ScaleOperator& so) { _scale = so._scale; return *this; }
1335
1336    float _scale;
1337
1338    inline void luminance(float& l) const { l*= _scale; }
1339    inline void alpha(float& a) const { a*= _scale; }
1340    inline void luminance_alpha(float& l,float& a) const { l*= _scale; a*= _scale;  }
1341    inline void rgb(float& r,float& g,float& b) const { r*= _scale; g*=_scale; b*=_scale; }
1342    inline void rgba(float& r,float& g,float& b,float& a) const { r*= _scale; g*=_scale; b*=_scale; a*=_scale; }
1343};
1344
1345struct RecordRowOperator
1346{
1347    RecordRowOperator(unsigned int num):_colours(num),_pos(0) {}
1348
1349    mutable std::vector<osg::Vec4>  _colours;
1350    mutable unsigned int            _pos;
1351   
1352    inline void luminance(float l) const { rgba(l,l,l,1.0f); }
1353    inline void alpha(float a) const { rgba(1.0f,1.0f,1.0f,a); }
1354    inline void luminance_alpha(float l,float a) const { rgba(l,l,l,a);  }
1355    inline void rgb(float r,float g,float b) const { rgba(r,g,b,1.0f); }
1356    inline void rgba(float r,float g,float b,float a) const { _colours[_pos++].set(r,g,b,a); }
1357};
1358
1359struct WriteRowOperator
1360{
1361    WriteRowOperator():_pos(0) {}
1362    WriteRowOperator(unsigned int num):_colours(num),_pos(0) {}
1363
1364    std::vector<osg::Vec4>  _colours;
1365    mutable unsigned int    _pos;
1366   
1367    inline void luminance(float& l) const { l = _colours[_pos++].r(); }
1368    inline void alpha(float& a) const { a = _colours[_pos++].a(); }
1369    inline void luminance_alpha(float& l,float& a) const { l = _colours[_pos].r(); a = _colours[_pos++].a(); }
1370    inline void rgb(float& r,float& g,float& b) const { r = _colours[_pos].r(); g = _colours[_pos].g(); b = _colours[_pos].b(); }
1371    inline void rgba(float& r,float& g,float& b,float& a) const {  r = _colours[_pos].r(); g = _colours[_pos].g(); b = _colours[_pos].b(); a = _colours[_pos++].a(); }
1372};
1373
1374osg::Image* readRaw(int sizeX, int sizeY, int sizeZ, int numberBytesPerComponent, int numberOfComponents, const std::string& endian, const std::string& raw_filename)
1375{
1376    std::ifstream fin(raw_filename.c_str(), std::ifstream::binary);
1377    if (!fin) return 0;
1378
1379    GLenum pixelFormat;
1380    switch(numberOfComponents)
1381    {
1382        case 1 : pixelFormat = GL_LUMINANCE; break;
1383        case 2 : pixelFormat = GL_LUMINANCE_ALPHA; break;
1384        case 3 : pixelFormat = GL_RGB; break;
1385        case 4 : pixelFormat = GL_RGBA; break;
1386        default :
1387            osg::notify(osg::NOTICE)<<"Error: numberOfComponents="<<numberOfComponents<<" not supported, only 1,2,3 or 4 are supported."<<std::endl;
1388            return 0;
1389    }
1390
1391   
1392    GLenum dataType;
1393    switch(numberBytesPerComponent)
1394    {
1395        case 1 : dataType = GL_UNSIGNED_BYTE; break;
1396        case 2 : dataType = GL_UNSIGNED_SHORT; break;
1397        case 4 : dataType = GL_UNSIGNED_INT; break;
1398        default :
1399            osg::notify(osg::NOTICE)<<"Error: numberBytesPerComponent="<<numberBytesPerComponent<<" not supported, only 1,2 or 4 are supported."<<std::endl;
1400            return 0;
1401    }
1402   
1403    int s_maximumTextureSize=256, t_maximumTextureSize=256, r_maximumTextureSize=256;
1404   
1405    int sizeS = sizeX;
1406    int sizeT = sizeY;
1407    int sizeR = sizeZ;
1408    clampToNearestValidPowerOfTwo(sizeS, sizeT, sizeR, s_maximumTextureSize, t_maximumTextureSize, r_maximumTextureSize);
1409
1410    osg::ref_ptr<osg::Image> image = new osg::Image;
1411    image->allocateImage(sizeS, sizeT, sizeR, pixelFormat, dataType);
1412   
1413   
1414    bool endianSwap = (osg::getCpuByteOrder()==osg::BigEndian) ? (endian!="big") : (endian=="big");
1415   
1416    unsigned int r_offset = (sizeZ<sizeR) ? sizeR/2 - sizeZ/2 : 0;
1417   
1418    int offset = endianSwap ? numberBytesPerComponent : 0;
1419    int delta = endianSwap ? -1 : 1;
1420    for(int r=0;r<sizeZ;++r)
1421    {
1422        for(int t=0;t<sizeY;++t)
1423        {
1424            char* data = (char*) image->data(0,t,r+r_offset);
1425            for(int s=0;s<sizeX;++s)
1426            {
1427                if (!fin) return 0;
1428               
1429                for(int c=0;c<numberOfComponents;++c)
1430                {
1431                    char* ptr = data+offset;
1432                    for(int b=0;b<numberBytesPerComponent;++b)
1433                    {
1434                        fin.read((char*)ptr, 1);
1435                        ptr += delta;
1436                    }
1437                    data += numberBytesPerComponent;
1438                }
1439            }
1440        }
1441    }
1442
1443
1444    // normalise texture
1445    {
1446        // compute range of values
1447        FindRangeOperator rangeOp;   
1448        readImage(image.get(), rangeOp);
1449        modifyImage(image.get(),ScaleOperator(1.0f/rangeOp._rmax));
1450    }
1451   
1452   
1453    fin.close();
1454
1455    if (dataType!=GL_UNSIGNED_BYTE)
1456    {
1457        // need to convert to ubyte
1458       
1459        osg::ref_ptr<osg::Image> new_image = new osg::Image;
1460        new_image->allocateImage(sizeS, sizeT, sizeR, pixelFormat, GL_UNSIGNED_BYTE);
1461       
1462        RecordRowOperator readOp(sizeS);
1463        WriteRowOperator writeOp;
1464
1465        for(int r=0;r<sizeR;++r)
1466        {
1467            for(int t=0;t<sizeT;++t)
1468            {
1469                // reset the indices to beginning
1470                readOp._pos = 0;
1471                writeOp._pos = 0;
1472           
1473                // read the pixels into readOp's _colour array
1474                readRow(sizeS, pixelFormat, dataType, image->data(0,t,r), readOp);
1475                               
1476                // pass readOp's _colour array contents over to writeOp (note this is just a pointer swap).
1477                writeOp._colours.swap(readOp._colours);
1478               
1479                modifyRow(sizeS, pixelFormat, GL_UNSIGNED_BYTE, new_image->data(0,t,r), writeOp);
1480
1481                // return readOp's _colour array contents back to its rightful owner.
1482                writeOp._colours.swap(readOp._colours);
1483            }
1484        }
1485       
1486        image = new_image;
1487    }
1488   
1489    return image.release();
1490   
1491   
1492}
1493
1494enum ColourSpaceOperation
1495{
1496    NO_COLOUR_SPACE_OPERATION,
1497    MODULATE_ALPHA_BY_LUMINANCE,
1498    MODULATE_ALPHA_BY_COLOUR,
1499    REPLACE_ALPHA_WITH_LUMINACE
1500};
1501
1502struct ModulateAlphaByLuminanceOperator
1503{
1504    ModulateAlphaByLuminanceOperator() {}
1505
1506    inline void luminance(float&) const {}
1507    inline void alpha(float&) const {}
1508    inline void luminance_alpha(float& l,float& a) const { a*= l; }
1509    inline void rgb(float&,float&,float&) const {}
1510    inline void rgba(float& r,float& g,float& b,float& a) const { float l = (r+g+b)*0.3333333; a *= l;}
1511};
1512
1513struct ModulateAlphaByColourOperator
1514{
1515    ModulateAlphaByColourOperator(const osg::Vec4& colour):_colour(colour) { _lum = _colour.length(); }
1516   
1517    osg::Vec4 _colour;
1518    float _lum;
1519
1520    inline void luminance(float&) const {}
1521    inline void alpha(float&) const {}
1522    inline void luminance_alpha(float& l,float& a) const { a*= l*_lum; }
1523    inline void rgb(float&,float&,float&) const {}
1524    inline void rgba(float& r,float& g,float& b,float& a) const { a = (r*_colour.r()+g*_colour.g()+b*_colour.b()+a*_colour.a()); }
1525};
1526
1527struct ReplaceAlphaWithLuminanceOperator
1528{
1529    ReplaceAlphaWithLuminanceOperator() {}
1530
1531    inline void luminance(float&) const {}
1532    inline void alpha(float&) const {}
1533    inline void luminance_alpha(float& l,float& a) const { a= l; }
1534    inline void rgb(float&,float&,float&) const { }
1535    inline void rgba(float& r,float& g,float& b,float& a) const { float l = (r+g+b)*0.3333333; a = l; }
1536};
1537
1538void doColourSpaceConversion(ColourSpaceOperation op, osg::Image* image, osg::Vec4& colour)
1539{
1540    switch(op)
1541    {
1542        case (MODULATE_ALPHA_BY_LUMINANCE):
1543            std::cout<<"doing conversion MODULATE_ALPHA_BY_LUMINANCE"<<std::endl;
1544            modifyImage(image,ModulateAlphaByLuminanceOperator());
1545            break;
1546        case (MODULATE_ALPHA_BY_COLOUR):
1547            std::cout<<"doing conversion MODULATE_ALPHA_BY_COLOUR"<<std::endl;
1548            modifyImage(image,ModulateAlphaByColourOperator(colour));
1549            break;
1550        case (REPLACE_ALPHA_WITH_LUMINACE):
1551            std::cout<<"doing conversion REPLACE_ALPHA_WITH_LUMINACE"<<std::endl;
1552            modifyImage(image,ReplaceAlphaWithLuminanceOperator());
1553            break;
1554        default:
1555            break;
1556    }
1557}
1558
1559
1560struct ApplyTransferFunctionOperator
1561{
1562    ApplyTransferFunctionOperator(osg::TransferFunction1D* tf, unsigned char* data):
1563        _tf(tf),
1564        _data(data) {}
1565       
1566    inline void luminance(float l) const
1567    {
1568        osg::Vec4 c = _tf->getInterpolatedValue(l);
1569        //std::cout<<"l = "<<l<<" c="<<c<<std::endl;
1570        *(_data++) = (unsigned char)(c[0]*255.0f + 0.5f);
1571        *(_data++) = (unsigned char)(c[1]*255.0f + 0.5f);
1572        *(_data++) = (unsigned char)(c[2]*255.0f + 0.5f);
1573        *(_data++) = (unsigned char)(c[3]*255.0f + 0.5f);
1574    }
1575     
1576    inline void alpha(float a) const
1577    {
1578        luminance(a);
1579    }
1580   
1581    inline void luminance_alpha(float l,float a) const
1582    {
1583        luminance(l);
1584    }
1585     
1586    inline void rgb(float r,float g,float b) const
1587    {
1588        luminance((r+g+b)*0.3333333);
1589    }
1590   
1591    inline void rgba(float r,float g,float b,float a) const
1592    {
1593        luminance(a);
1594    }
1595   
1596    mutable osg::ref_ptr<osg::TransferFunction1D> _tf;
1597    mutable unsigned char* _data;
1598};
1599
1600osg::Image* applyTransferFunction(osg::Image* image, osg::TransferFunction1D* transferFunction)
1601{
1602    std::cout<<"Applying transfer function"<<std::endl;
1603    osg::Image* output_image = new osg::Image;
1604    output_image->allocateImage(image->s(),image->t(), image->r(), GL_RGBA, GL_UNSIGNED_BYTE);
1605   
1606    readImage(image,ApplyTransferFunctionOperator(transferFunction, output_image->data()));
1607   
1608    return output_image;
1609}
1610
1611osg::TransferFunction1D* readTransferFunctionFile(const std::string& filename)
1612{
1613    std::string foundFile = osgDB::findDataFile(filename);
1614    if (foundFile.empty())
1615    {
1616        std::cout<<"Error: could not find transfer function file : "<<filename<<std::endl;
1617        return 0;
1618    }
1619   
1620    std::cout<<"Reading transfer function "<<filename<<std::endl;
1621
1622    osg::TransferFunction1D::ValueMap valueMap;
1623    std::ifstream fin(foundFile.c_str());
1624    while(fin)
1625    {
1626        float value, red, green, blue, alpha;
1627        fin >> value >> red >> green >> blue >> alpha;
1628        if (fin)
1629        {
1630            std::cout<<"value = "<<value<<" ("<<red<<", "<<green<<", "<<blue<<", "<<alpha<<")"<<std::endl;
1631            valueMap[value] = osg::Vec4(red,green,blue,alpha);
1632        }
1633    }
1634   
1635    if (valueMap.empty())
1636    {
1637        std::cout<<"Error: No values read from transfer function file: "<<filename<<std::endl;
1638        return 0;
1639    }
1640   
1641    osg::TransferFunction1D* tf = new osg::TransferFunction1D;
1642    tf->assign(valueMap, true);
1643   
1644    return tf;
1645}
1646
1647
1648class TestSupportOperation: public osg::GraphicsOperation
1649{
1650public:
1651
1652    TestSupportOperation():
1653        osg::GraphicsOperation("TestSupportOperation",false),
1654        supported(true),
1655        errorMessage(),
1656        maximumTextureSize(256) {}
1657
1658    virtual void operator () (osg::GraphicsContext* gc)
1659    {
1660        OpenThreads::ScopedLock<OpenThreads::Mutex> lock(mutex);
1661
1662        glGetIntegerv( GL_MAX_3D_TEXTURE_SIZE, &maximumTextureSize );
1663       
1664        osg::notify(osg::NOTICE)<<"Max texture size="<<maximumTextureSize<<std::endl;
1665    }
1666       
1667    OpenThreads::Mutex  mutex;
1668    bool                supported;
1669    std::string         errorMessage;
1670    GLint               maximumTextureSize;
1671};
1672
1673
1674
1675int main( int argc, char **argv )
1676{
1677    // use an ArgumentParser object to manage the program arguments.
1678    osg::ArgumentParser arguments(&argc,argv);
1679   
1680    // set up the usage document, in case we need to print out how to use this program.
1681    arguments.getApplicationUsage()->setDescription(arguments.getApplicationName()+" is the example which demonstrates use of 3D textures.");
1682    arguments.getApplicationUsage()->setCommandLineUsage(arguments.getApplicationName()+" [options] filename ...");
1683    arguments.getApplicationUsage()->addCommandLineOption("-h or --help","Display this information");
1684    arguments.getApplicationUsage()->addCommandLineOption("-n","Create normal map for per voxel lighting.");
1685    arguments.getApplicationUsage()->addCommandLineOption("-s <numSlices>","Number of slices to create.");
1686    arguments.getApplicationUsage()->addCommandLineOption("--images [filenames]","Specify a stack of 2d images to build the 3d volume from.");
1687    arguments.getApplicationUsage()->addCommandLineOption("--shader","Use OpenGL Shading Language.");
1688    arguments.getApplicationUsage()->addCommandLineOption("--mip","Use Maximum Intensity Projection (MIP) filtering.");
1689    arguments.getApplicationUsage()->addCommandLineOption("--xSize <size>","Relative width of rendered brick.");
1690    arguments.getApplicationUsage()->addCommandLineOption("--ySize <size>","Relative length of rendered brick.");
1691    arguments.getApplicationUsage()->addCommandLineOption("--zSize <size>","Relative height of rendered brick.");
1692    arguments.getApplicationUsage()->addCommandLineOption("--xMultiplier <multiplier>","Tex coord x mulitplier.");
1693    arguments.getApplicationUsage()->addCommandLineOption("--yMultiplier <multiplier>","Tex coord y mulitplier.");
1694    arguments.getApplicationUsage()->addCommandLineOption("--zMultiplier <multiplier>","Tex coord z mulitplier.");
1695    arguments.getApplicationUsage()->addCommandLineOption("--clip <ratio>","clip volume as a ratio, 0.0 clip all, 1.0 clip none.");
1696    arguments.getApplicationUsage()->addCommandLineOption("--maxTextureSize <size>","Set the texture maximum resolution in the s,t,r (x,y,z) dimensions.");
1697    arguments.getApplicationUsage()->addCommandLineOption("--s_maxTextureSize <size>","Set the texture maximum resolution in the s (x) dimension.");
1698    arguments.getApplicationUsage()->addCommandLineOption("--t_maxTextureSize <size>","Set the texture maximum resolution in the t (y) dimension.");
1699    arguments.getApplicationUsage()->addCommandLineOption("--r_maxTextureSize <size>","Set the texture maximum resolution in the r (z) dimension.");
1700    arguments.getApplicationUsage()->addCommandLineOption("--compressed","Enable the usage of compressed textures.");
1701    arguments.getApplicationUsage()->addCommandLineOption("--compressed-arb","Enable the usage of OpenGL ARB compressed textures.");
1702    arguments.getApplicationUsage()->addCommandLineOption("--compressed-dxt1","Enable the usage of S3TC DXT1 compressed textures.");
1703    arguments.getApplicationUsage()->addCommandLineOption("--compressed-dxt3","Enable the usage of S3TC DXT3 compressed textures.");
1704    arguments.getApplicationUsage()->addCommandLineOption("--compressed-dxt5","Enable the usage of S3TC DXT5 compressed textures.");
1705    arguments.getApplicationUsage()->addCommandLineOption("--modulate-alpha-by-luminance","For each pixel multiple the alpha value by the luminance.");
1706    arguments.getApplicationUsage()->addCommandLineOption("--replace-alpha-with-luminance","For each pixel mSet the alpha value to the luminance.");
1707    arguments.getApplicationUsage()->addCommandLineOption("--num-components <num>","Set the number of components to in he target image.");
1708//    arguments.getApplicationUsage()->addCommandLineOption("--raw <sizeX> <sizeY> <sizeZ> <numberBytesPerComponent> <numberOfComponents> <endian> <filename>","read a raw image data");
1709
1710    // construct the viewer.
1711    osgViewer::Viewer viewer(arguments);
1712
1713    // add the window size toggle handler
1714    viewer.addEventHandler(new osgViewer::WindowSizeHandler);
1715       
1716    // add the stats handler
1717    viewer.addEventHandler(new osgViewer::StatsHandler);
1718
1719    viewer.getCamera()->setClearColor(osg::Vec4(0.0f,0.0f,0.0f,0.0f));
1720
1721    // if user request help write it out to cout.
1722    if (arguments.read("-h") || arguments.read("--help"))
1723    {
1724        arguments.getApplicationUsage()->write(std::cout);
1725        return 1;
1726    }
1727
1728    std::string outputFile;
1729    while (arguments.read("-o",outputFile)) {}
1730
1731
1732
1733    osg::ref_ptr<osg::TransferFunction1D> transferFunction;
1734    std::string tranferFunctionFile;
1735    while (arguments.read("--tf",tranferFunctionFile))
1736    {
1737        transferFunction = readTransferFunctionFile(tranferFunctionFile);
1738    }
1739
1740    unsigned int numSlices=500;
1741    while (arguments.read("-s",numSlices)) {}
1742   
1743   
1744    float sliceEnd=1.0f;
1745    while (arguments.read("--clip",sliceEnd)) {}
1746
1747    float alphaFunc=0.02f;
1748    while (arguments.read("--alphaFunc",alphaFunc)) {}
1749
1750
1751    bool createNormalMap = false;
1752    while (arguments.read("-n")) createNormalMap=true;
1753   
1754    bool maximumIntensityProjection = false;
1755    while(arguments.read("--mip")) maximumIntensityProjection = true;
1756
1757    float xSize=1.0f, ySize=1.0f, zSize=1.0f;
1758    while (arguments.read("--xSize",xSize)) {}
1759    while (arguments.read("--ySize",ySize)) {}
1760    while (arguments.read("--zSize",zSize)) {}
1761
1762    float xMultiplier=1.0f, yMultiplier=1.0f, zMultiplier=1.0f;
1763    while (arguments.read("--xMultiplier",xMultiplier)) {}
1764    while (arguments.read("--yMultiplier",yMultiplier)) {}
1765    while (arguments.read("--zMultiplier",zMultiplier)) {}
1766
1767    osg::ref_ptr<TestSupportOperation> testSupportOperation = new TestSupportOperation;
1768    viewer.setRealizeOperation(testSupportOperation.get());
1769   
1770    viewer.realize();
1771
1772    int maximumTextureSize = testSupportOperation->maximumTextureSize;
1773    int s_maximumTextureSize = maximumTextureSize;
1774    int t_maximumTextureSize = maximumTextureSize;
1775    int r_maximumTextureSize = maximumTextureSize;
1776    while(arguments.read("--maxTextureSize",maximumTextureSize))
1777    {
1778        s_maximumTextureSize = maximumTextureSize;
1779        t_maximumTextureSize = maximumTextureSize;
1780        r_maximumTextureSize = maximumTextureSize;
1781    }
1782    while(arguments.read("--s_maxTextureSize",s_maximumTextureSize)) {}
1783    while(arguments.read("--t_maxTextureSize",t_maximumTextureSize)) {}
1784    while(arguments.read("--r_maxTextureSize",r_maximumTextureSize)) {}
1785
1786    osg::Texture::InternalFormatMode internalFormatMode = osg::Texture::USE_IMAGE_DATA_FORMAT;
1787    while(arguments.read("--compressed") || arguments.read("--compressed-arb")) { internalFormatMode = osg::Texture::USE_ARB_COMPRESSION; }
1788
1789    while(arguments.read("--compressed-dxt1")) { internalFormatMode = osg::Texture::USE_S3TC_DXT1_COMPRESSION; }
1790    while(arguments.read("--compressed-dxt3")) { internalFormatMode = osg::Texture::USE_S3TC_DXT3_COMPRESSION; }
1791    while(arguments.read("--compressed-dxt5")) { internalFormatMode = osg::Texture::USE_S3TC_DXT5_COMPRESSION; }
1792   
1793   
1794    // set up colour space operation.
1795    ColourSpaceOperation colourSpaceOperation = NO_COLOUR_SPACE_OPERATION;
1796    osg::Vec4 colourModulate(0.25f,0.25f,0.25f,0.25f);
1797    while(arguments.read("--modulate-alpha-by-luminance")) { colourSpaceOperation = MODULATE_ALPHA_BY_LUMINANCE; }
1798    while(arguments.read("--modulate-alpha-by-colour", colourModulate.x(),colourModulate.y(),colourModulate.z(),colourModulate.w() )) { colourSpaceOperation = MODULATE_ALPHA_BY_COLOUR; }
1799    while(arguments.read("--replace-alpha-with-luminance")) { colourSpaceOperation = REPLACE_ALPHA_WITH_LUMINACE; }
1800       
1801    bool resizeToPowerOfTwo = false;
1802   
1803    unsigned int numComponentsDesired = 0;
1804    while(arguments.read("--num-components", numComponentsDesired)) {}
1805
1806    bool useShader = false;
1807    while(arguments.read("--shader")) { useShader = true; }
1808
1809    osg::ref_ptr<osg::Image> image_3d;
1810
1811    int sizeX, sizeY, sizeZ, numberBytesPerComponent, numberOfComponents;
1812    std::string endian, raw_filename;
1813    while (arguments.read("--raw", sizeX, sizeY, sizeZ, numberBytesPerComponent, numberOfComponents, endian, raw_filename))
1814    {
1815        image_3d = readRaw(sizeX, sizeY, sizeZ, numberBytesPerComponent, numberOfComponents, endian, raw_filename);
1816    }
1817
1818    while (arguments.read("--images"))
1819    {
1820        ImageList imageList;
1821        for(int pos=1;pos<arguments.argc() && !arguments.isOption(pos);++pos)
1822        {
1823            // not an option so assume string is a filename.
1824            osg::Image *image = osgDB::readImageFile( arguments[pos]);
1825
1826            if(image)
1827            {
1828                imageList.push_back(image);
1829            }
1830        }
1831       
1832        // pack the textures into a single texture.
1833        ProcessRow processRow;
1834        image_3d = createTexture3D(imageList, processRow, numComponentsDesired, s_maximumTextureSize, t_maximumTextureSize, r_maximumTextureSize, resizeToPowerOfTwo);
1835    }
1836
1837
1838    // any option left unread are converted into errors to write out later.
1839    arguments.reportRemainingOptionsAsUnrecognized();
1840
1841    // report any errors if they have occurred when parsing the program arguments.
1842    if (arguments.errors())
1843    {
1844        arguments.writeErrorMessages(std::cout);
1845        return 1;
1846    }
1847
1848    // assume remaining arguments are file names of textures.
1849    for(int pos=1;pos<arguments.argc() && !image_3d;++pos)
1850    {
1851        if (!arguments.isOption(pos))
1852        {
1853            std::string filename = arguments[pos];
1854            if (osgDB::getLowerCaseFileExtension(filename)=="dicom")
1855            {
1856                image_3d = osgDB::readImageFile( filename );
1857            }
1858            else
1859            {
1860                osgDB::FileType fileType = osgDB::fileType(filename);
1861                if (fileType == osgDB::FILE_NOT_FOUND)
1862                {
1863                    filename = osgDB::findDataFile(filename);
1864                    fileType = osgDB::fileType(filename);
1865                }
1866
1867                if (fileType == osgDB::DIRECTORY)
1868                {
1869                    osgDB::DirectoryContents contents = osgDB::getDirectoryContents(filename);
1870
1871                    std::sort(contents.begin(), contents.end());
1872
1873                    ImageList imageList;
1874                    for(osgDB::DirectoryContents::iterator itr = contents.begin();
1875                        itr != contents.end();
1876                        ++itr)
1877                    {
1878                        std::string localFile = filename + "/" + *itr;
1879                        std::cout<<"contents = "<<localFile<<std::endl;
1880                        if (osgDB::fileType(localFile) == osgDB::REGULAR_FILE)
1881                        {
1882                            // not an option so assume string is a filename.
1883                            osg::Image *image = osgDB::readImageFile(localFile);
1884                            if(image)
1885                            {
1886                                imageList.push_back(image);
1887                            }
1888                        }
1889                    }
1890
1891                    // pack the textures into a single texture.
1892                    ProcessRow processRow;
1893                    image_3d = createTexture3D(imageList, processRow, numComponentsDesired, s_maximumTextureSize, t_maximumTextureSize, r_maximumTextureSize, resizeToPowerOfTwo);
1894
1895                }
1896                else if (fileType == osgDB::REGULAR_FILE)
1897                {
1898                    // not an option so assume string is a filename.
1899                    image_3d = osgDB::readImageFile( filename );
1900                }
1901                else
1902                {
1903                    osg::notify(osg::NOTICE)<<"Error: could not find file: "<<filename<<std::endl;
1904                    return 1;
1905                }
1906            }           
1907        }
1908    }
1909   
1910    if (!image_3d)
1911    {
1912        std::cout<<"No model loaded, please specify and volumetric image file on the command line."<<std::endl;
1913        return 1;
1914    }
1915   
1916    if (colourSpaceOperation!=NO_COLOUR_SPACE_OPERATION)
1917    {
1918        doColourSpaceConversion(colourSpaceOperation, image_3d.get(), colourModulate);
1919    }
1920   
1921    if (transferFunction.valid())
1922    {
1923        image_3d = applyTransferFunction(image_3d.get(), transferFunction.get());
1924    }
1925   
1926    osg::ref_ptr<osg::Image> normalmap_3d = createNormalMap ? createNormalMapTexture(image_3d.get()) : 0;
1927
1928
1929    osg::RefMatrix* matrix = dynamic_cast<osg::RefMatrix*>(image_3d->getUserData());
1930    if (matrix)
1931    {
1932        osg::notify(osg::NOTICE)<<"Image has Matrix = "<<*matrix<<std::endl;
1933        xSize = image_3d->s() * (*matrix)(0,0);
1934        ySize = image_3d->t() * (*matrix)(1,1);
1935        zSize = image_3d->r() * (*matrix)(2,2);
1936    }
1937
1938
1939    // create a model from the images.
1940    osg::Node* rootNode = 0;
1941   
1942    if (useShader)
1943    {
1944        rootNode = createShaderModel(image_3d, normalmap_3d,
1945                               internalFormatMode,
1946                               xSize, ySize, zSize,
1947                               xMultiplier, yMultiplier, zMultiplier,
1948                               numSlices, sliceEnd, alphaFunc, maximumIntensityProjection);
1949    }
1950    else
1951    {
1952        rootNode = createModel(image_3d, normalmap_3d,
1953                               internalFormatMode,
1954                               xSize, ySize, zSize,
1955                               xMultiplier, yMultiplier, zMultiplier,
1956                               numSlices, sliceEnd, alphaFunc, maximumIntensityProjection);
1957    }
1958   
1959    if (!outputFile.empty())
1960    {   
1961        std::string ext = osgDB::getFileExtension(outputFile);
1962        std::string name_no_ext = osgDB::getNameLessExtension(outputFile);
1963        if (ext=="osg")
1964        {
1965            if (image_3d.valid())
1966            {
1967                image_3d->setFileName(name_no_ext + ".dds");           
1968                osgDB::writeImageFile(*image_3d, image_3d->getFileName());
1969            }
1970            if (normalmap_3d.valid())
1971            {
1972                normalmap_3d->setFileName(name_no_ext + "_normalmap.dds");           
1973                osgDB::writeImageFile(*normalmap_3d, normalmap_3d->getFileName());
1974            }
1975           
1976            osgDB::writeNodeFile(*rootNode, outputFile);
1977        }
1978        else if (ext=="ive")
1979        {
1980            osgDB::writeNodeFile(*rootNode, outputFile);       
1981        }
1982        else if (ext=="dds")
1983        {
1984            osgDB::writeImageFile(*image_3d, outputFile);       
1985        }
1986        else
1987        {
1988            std::cout<<"Extension not support for file output, not file written."<<std::endl;
1989        }
1990       
1991        return 0;
1992    }
1993
1994
1995    if (rootNode)
1996    {
1997
1998        // set the scene to render
1999        viewer.setSceneData(rootNode);
2000       
2001        // the the viewers main frame loop
2002        viewer.run();
2003    }   
2004   
2005    return 0;
2006
2007}
Note: See TracBrowser for help on using the browser.