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

Revision 8951, 81.9 kB (checked in by robert, 6 years ago)

Added setting of various filters and setResizeNonPowerOfTwoHint(false) on normals maps etc.

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