mitk: mitkColourImageProcessor.cpp Source File

Go to the documentation of this file.
00001 /*=========================================================================
00002 
00003 Program:   Medical Imaging & Interaction Toolkit
00004 Language:  C++
00005 Date:      $Date: 2009-05-28 17:19:30 +0200 (Do, 28 Mai 2009) $
00006 Version:   $Revision: 17495 $ 
00007  
00008 Copyright (c) German Cancer Research Center, Division of Medical and
00009 Biological Informatics. All rights reserved.
00010 See MITKCopyright.txt or https://www.mitk.org/copyright.html for details.
00011 
00012 This software is distributed WITHOUT ANY WARRANTY; without even
00013 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
00014 PURPOSE.  See the above copyright notices for more information.
00015 
00016 =========================================================================*/
00017 
00018 #include "mitkColourImageProcessor.h"
00019 #include <mitkImage.h>
00020 #include <itkImage.h>
00021 #include <mitkImageCast.h>
00022 
00023 
00024 
00025 namespace mitk
00026 {
00027 
00028 mitkColourImageProcessor::mitkColourImageProcessor()
00029 {
00030 }
00031 
00032 mitkColourImageProcessor::~mitkColourImageProcessor()
00033 {
00034 }
00035 
00036 template <class T> class ScalarToRGBAConverter
00037 {
00038   const T *dataPtr;
00039   unsigned char *tmpPtr;
00040   int sizeX;
00041   int sizeY;
00042   int sizeZ;
00043   int sizeXY;
00044   int sizeXm1;
00045   int sizeYm1;
00046   int sizeZm1;
00047   mitk::TransferFunction::Pointer tf;
00048   
00049   public:
00050 
00051   ScalarToRGBAConverter( const T *_dataPtr,unsigned char *_tmpPtr,int _sizeX,int _sizeY,int _sizeZ,mitk::TransferFunction::Pointer _tf )
00052   {
00053     dataPtr=_dataPtr;
00054     tmpPtr=_tmpPtr;
00055     sizeX=_sizeX;
00056     sizeY=_sizeY;
00057     sizeZ=_sizeZ;
00058 
00059     sizeXY=sizeX*sizeY;
00060     sizeXm1=sizeX-1;
00061     sizeYm1=sizeY-1;
00062     sizeZm1=sizeZ-1;
00063 
00064     tf=_tf;
00065   }
00066 
00067   inline float sample(int x,int y,int z)
00068   {
00069     return float(dataPtr[ x + y * sizeX + z * sizeXY ]);
00070   }
00071   
00072   inline int clamp(int x)
00073   {
00074     if(x<0) x=0; else if(x>255) x=255;
00075     return x;
00076   }
00077   
00078   inline void write(int x,int y,int z,float grayValue,float gx,float gy,float gz)
00079   {
00080   
00081  /*       
00082     gx /= aspect[0];
00083     gy /= aspect[1];
00084     gz /= aspect[2];
00085  */
00086     // Compute the gradient magnitude
00087 
00088     float t = sqrtf( gx*gx + gy*gy + gz*gz );
00089 
00090     int doff = x + y * sizeX + z * sizeXY;
00091     
00092     vtkPiecewiseFunction* opacityTransferFunction = tf->GetScalarOpacityFunction();
00093     vtkPiecewiseFunction* gradientTransferFunction = tf->GetGradientOpacityFunction();
00094     vtkColorTransferFunction* colorTransferFunction = tf->GetColorTransferFunction();
00095 
00096     double rgb[3];
00097     colorTransferFunction->GetColor( double(grayValue), rgb);
00098     
00099     double opacity= opacityTransferFunction->GetValue( double(grayValue) );
00100     
00101     opacity *= gradientTransferFunction->GetValue( double(0.5f*t) );
00102     
00103     tmpPtr[doff*4+0] = int( rgb[0]*255 + 0.5 );
00104     tmpPtr[doff*4+1] = int( rgb[1]*255 + 0.5 );
00105     tmpPtr[doff*4+2] = int( rgb[2]*255 + 0.5 );
00106     tmpPtr[doff*4+3] = int( opacity*255 + 0.5 );
00107   }
00108   
00109   
00110   inline void compute(int x,int y,int z)
00111   {
00112     float grayValue = sample(x,y,z);
00113     float gx,gy,gz;
00114 
00115     gx = sample(x+1,y,z) - sample(x-1,y,z);
00116     gy = sample(x,y+1,z) - sample(x,y-1,z);
00117     gz = sample(x,y,z+1) - sample(x,y,z-1);
00118     
00119     write( x, y, z, grayValue, gx, gy, gz );
00120  
00121   }
00122 
00123   inline void computeClamp(int x,int y,int z)
00124   {
00125     float grayValue = sample(x,y,z);
00126     float gx,gy,gz;
00127 
00128     if(x==0)            gx = 2.0f * ( sample(x+1,y,z) - grayValue       );
00129     else if(x==sizeXm1) gx = 2.0f * ( grayValue       - sample(x-1,y,z) );
00130     else                gx =          sample(x+1,y,z) - sample(x-1,y,z);
00131     
00132     if(y==0)            gy = 2.0f * ( sample(x,y+1,z) - grayValue       );
00133     else if(y==sizeYm1) gy = 2.0f * ( grayValue       - sample(x,y-1,z) );
00134     else                gy =          sample(x,y+1,z) - sample(x,y-1,z);
00135     
00136     if(z==0)            gz = 2.0f * ( sample(x,y,z+1) - grayValue       );
00137     else if(z==sizeZm1) gz = 2.0f * ( grayValue       - sample(x,y,z-1) );
00138     else                gz =          sample(x,y,z+1) - sample(x,y,z-1);
00139     
00140     write( x, y, z, grayValue, gx, gy, gz );
00141   }
00142 
00143   inline void compute1D(int y,int z)
00144   {
00145     int x;
00146     
00147     x=0;
00148     computeClamp(x,y,z);
00149     x++;
00150     
00151     while(x<sizeX-1)
00152     {
00153       compute(x,y,z);     
00154       x++;
00155     }
00156     
00157     if(x<sizeX)
00158     {
00159       computeClamp(x,y,z);
00160       x++;
00161     }
00162   }      
00163 
00164   inline void computeClamp1D(int y,int z)
00165   {
00166     int x;
00167     
00168     x=0;
00169     
00170     while(x<sizeX)
00171     {
00172       computeClamp(x,y,z);     
00173       x++;
00174     }
00175   }      
00176         
00177   inline void computeClamp2D(int z)
00178   {
00179     int y;
00180 
00181     y=0;
00182 
00183     while(y<sizeY)
00184     {
00185       computeClamp1D(y,z);     
00186       y++;
00187     }
00188   }
00189   
00190   inline void compute2D(int z)
00191   {
00192     int y;
00193 
00194     y=0;
00195     computeClamp1D(y,z);     
00196     y++;
00197     
00198     while(y<sizeY-1)
00199     {
00200       compute1D(y,z);     
00201       y++;
00202     }
00203     
00204     if(y<sizeY)
00205     {
00206       computeClamp1D(y,z);
00207       y++;
00208     }
00209   }
00210   
00211   inline void fillSlices()
00212   {
00213     int z;
00214 
00215     for(z=0;z<sizeZ;z++)
00216     {
00217       if(z==0 || z==sizeZ-1)
00218         computeClamp2D(z);
00219        else
00220         compute2D(z);
00221     }
00222   }    
00223 };
00224 
00225 
00226 
00227 template<class TType> mitk::Image::Pointer mitkColourImageProcessor::ScalarToRGBA( itk::Image<TType, 3>* input  , mitk::TransferFunction::Pointer tf)
00228 { 
00229   
00230   const TType *inputData=input->GetBufferPointer();
00231 
00232   typename itk::Image<TType, 3>::SizeType ioSize = input->GetLargestPossibleRegion().GetSize();
00233 
00234   MITK_INFO << "size input image: " << ioSize[0] << ", " << ioSize[1] << ", " << ioSize[2];
00235   MITK_INFO << "size voxel: " << ioSize[0] * ioSize[1] * ioSize[2];
00236   
00237   int voxel = ioSize[0] * ioSize[1] * ioSize[2];
00238   
00239   unsigned char* RGBABuffer = new unsigned char[4*voxel];
00240   
00241   // Convert volume  
00242   {
00243     ScalarToRGBAConverter<TType> strc(inputData,RGBABuffer,ioSize[0],ioSize[1],ioSize[2],tf);
00244     strc.fillSlices();
00245   }
00246     
00247   // Create MITK Image out of the raw data
00248   {
00249     mitk::Image::Pointer image = mitk::Image::New();
00250 
00251     unsigned int dimensions[ 3 ];
00252     dimensions[ 0 ] = ioSize[0];
00253     dimensions[ 1 ] = ioSize[1];
00254     dimensions[ 2 ] = ioSize[2];
00255 
00256     mitk::PixelType pixelType( typeid(RGBAPixel), 1, itk::ImageIOBase::RGBA );
00257     image->Initialize( pixelType, 3, dimensions );
00258     image->SetImportChannel( RGBABuffer, 0, Image::ManageMemory );
00259 
00260     return image;
00261   }
00262 }
00263 
00264  mitk::Image::Pointer mitkColourImageProcessor::convertToRGBAImage( mitk::Image::Pointer mitkInput , mitk::TransferFunction::Pointer tf )
00265 {
00266   MITK_INFO << "convertToRGBAImage" ;
00267   
00268   mitk::Image::Pointer mitkResult= mitk::Image::New();
00269  
00270   if (*mitkInput->GetPixelType().GetTypeId() == typeid(unsigned char))
00271   {
00272     //cast to itkImage
00273     itk::Image< unsigned char, 3 >::Pointer itkInput;
00274     mitk::CastToItkImage(mitkInput,itkInput);
00275     mitkResult = ScalarToRGBA<unsigned char>(itkInput, tf);
00276   }
00277   else if (*mitkInput->GetPixelType().GetTypeId() == typeid(short))
00278   {
00279     //cast to itkImage
00280     itk::Image< short, 3 >::Pointer itkInput;
00281     mitk::CastToItkImage(mitkInput,itkInput);
00282     mitkResult = ScalarToRGBA<short>(itkInput, tf);
00283   }
00284   else
00285   {
00286     MITK_ERROR << "unsupported pixel type";
00287     return NULL;
00288   }
00289 
00290 
00291    mitkResult->SetSpacing( mitkInput->GetGeometry()->GetSpacing() );
00292   
00293   return mitkResult;
00294 }
00295 
00296 
00297 
00298 
00299 
00300 
00301 
00302 
00303 
00304 
00305 
00306 
00307 template <class T,class B> 
00308 class ScalarBinaryToRGBAConverter
00309 {
00310   const T *dataPtr;
00311   const B *data2Ptr;
00312   unsigned char *tmpPtr;
00313   int sizeX;
00314   int sizeY;
00315   int sizeZ;
00316   int sizeXY;
00317   int sizeXm1;
00318   int sizeYm1;
00319   int sizeZm1;
00320   mitk::TransferFunction::Pointer tf;
00321   
00322   public:
00323 
00324   ScalarBinaryToRGBAConverter( const T *_dataPtr,const B *_data2Ptr,unsigned char *_tmpPtr,int _sizeX,int _sizeY,int _sizeZ,mitk::TransferFunction::Pointer _tf )
00325   {
00326     dataPtr=_dataPtr;
00327     data2Ptr=_data2Ptr;
00328     tmpPtr=_tmpPtr;
00329     sizeX=_sizeX;
00330     sizeY=_sizeY;
00331     sizeZ=_sizeZ;
00332 
00333     sizeXY=sizeX*sizeY;
00334     sizeXm1=sizeX-1;
00335     sizeYm1=sizeY-1;
00336     sizeZm1=sizeZ-1;
00337 
00338     tf=_tf;
00339   }
00340 
00341   inline float sample(int x,int y,int z)
00342   {
00343     return float(dataPtr[ x + y * sizeX + z * sizeXY ]);
00344   }
00345   
00346   inline bool sampleBinary(int x,int y,int z)
00347   {
00348     return data2Ptr[ x + y * sizeX + z * sizeXY ];
00349   }
00350   
00351   inline int clamp(int x)
00352   {
00353     if(x<0) x=0; else if(x>255) x=255;
00354     return x;
00355   }
00356   
00357   inline void write(int x,int y,int z,float grayValue,float gx,float gy,float gz)
00358   {
00359     if(sampleBinary(x,y,z))
00360     {
00361       
00362    /*       
00363       gx /= aspect[0];
00364       gy /= aspect[1];
00365       gz /= aspect[2];
00366    */
00367       // Compute the gradient magnitude
00368 
00369       float t = sqrtf( gx*gx + gy*gy + gz*gz );
00370 
00371       int doff = x + y * sizeX + z * sizeXY;
00372       
00373       vtkPiecewiseFunction* opacityTransferFunction = tf->GetScalarOpacityFunction();
00374       vtkPiecewiseFunction* gradientTransferFunction = tf->GetGradientOpacityFunction();
00375       vtkColorTransferFunction* colorTransferFunction = tf->GetColorTransferFunction();
00376 
00377       double rgb[3];
00378       colorTransferFunction->GetColor( double(grayValue), rgb);
00379       
00380       double opacity= opacityTransferFunction->GetValue( double(grayValue) );
00381       
00382       opacity *= gradientTransferFunction->GetValue( double(0.5f*t) );
00383       
00384       tmpPtr[doff*4+0] = int( rgb[0]*255 + 0.5 );
00385       tmpPtr[doff*4+1] = int( rgb[1]*255 + 0.5 );
00386       tmpPtr[doff*4+2] = int( rgb[2]*255 + 0.5 );
00387       tmpPtr[doff*4+3] = int( opacity*255 + 0.5 );
00388     }
00389     else
00390     {
00391       int doff = x + y * sizeX + z * sizeXY;
00392       tmpPtr[doff*4+0] = 0;
00393       tmpPtr[doff*4+1] = 0;
00394       tmpPtr[doff*4+2] = 0;
00395       tmpPtr[doff*4+3] = 0;
00396     }
00397   }
00398   
00399   
00400   inline void compute(int x,int y,int z)
00401   {
00402     float grayValue = sample(x,y,z);
00403     float gx,gy,gz;
00404 
00405     gx = sample(x+1,y,z) - sample(x-1,y,z);
00406     gy = sample(x,y+1,z) - sample(x,y-1,z);
00407     gz = sample(x,y,z+1) - sample(x,y,z-1);
00408     
00409     write( x, y, z, grayValue, gx, gy, gz );
00410  
00411   }
00412 
00413   inline void computeClamp(int x,int y,int z)
00414   {
00415     float grayValue = sample(x,y,z);
00416     float gx,gy,gz;
00417 
00418     if(x==0)            gx = 2.0f * ( sample(x+1,y,z) - grayValue       );
00419     else if(x==sizeXm1) gx = 2.0f * ( grayValue       - sample(x-1,y,z) );
00420     else                gx =          sample(x+1,y,z) - sample(x-1,y,z);
00421     
00422     if(y==0)            gy = 2.0f * ( sample(x,y+1,z) - grayValue       );
00423     else if(y==sizeYm1) gy = 2.0f * ( grayValue       - sample(x,y-1,z) );
00424     else                gy =          sample(x,y+1,z) - sample(x,y-1,z);
00425     
00426     if(z==0)            gz = 2.0f * ( sample(x,y,z+1) - grayValue       );
00427     else if(z==sizeZm1) gz = 2.0f * ( grayValue       - sample(x,y,z-1) );
00428     else                gz =          sample(x,y,z+1) - sample(x,y,z-1);
00429     
00430     write( x, y, z, grayValue, gx, gy, gz );
00431   }
00432 
00433   inline void compute1D(int y,int z)
00434   {
00435     int x;
00436     
00437     x=0;
00438     computeClamp(x,y,z);
00439     x++;
00440     
00441     while(x<sizeX-1)
00442     {
00443       compute(x,y,z);     
00444       x++;
00445     }
00446     
00447     if(x<sizeX)
00448     {
00449       computeClamp(x,y,z);
00450       x++;
00451     }
00452   }      
00453 
00454   inline void computeClamp1D(int y,int z)
00455   {
00456     int x;
00457     
00458     x=0;
00459     
00460     while(x<sizeX)
00461     {
00462       computeClamp(x,y,z);     
00463       x++;
00464     }
00465   }      
00466         
00467   inline void computeClamp2D(int z)
00468   {
00469     int y;
00470 
00471     y=0;
00472 
00473     while(y<sizeY)
00474     {
00475       computeClamp1D(y,z);     
00476       y++;
00477     }
00478   }
00479   
00480   inline void compute2D(int z)
00481   {
00482     int y;
00483 
00484     y=0;
00485     computeClamp1D(y,z);     
00486     y++;
00487     
00488     while(y<sizeY-1)
00489     {
00490       compute1D(y,z);     
00491       y++;
00492     }
00493     
00494     if(y<sizeY)
00495     {
00496       computeClamp1D(y,z);
00497       y++;
00498     }
00499   }
00500   
00501   inline void fillSlices()
00502   {
00503     int z;
00504 
00505     for(z=0;z<sizeZ;z++)
00506     {
00507       if(z==0 || z==sizeZ-1)
00508         computeClamp2D(z);
00509        else
00510         compute2D(z);
00511     }
00512   }    
00513 };
00514 
00515 
00516 
00517 template<class TType,class BType>   
00518 mitk::Image::Pointer mitkColourImageProcessor::ScalarAndBinaryToRGBA(itk::Image<TType, 3>* input  ,itk::Image<BType, 3>* input2  , mitk::TransferFunction::Pointer tf)
00519 {
00520   const TType *inputData=input->GetBufferPointer();
00521   const BType *input2Data=input2->GetBufferPointer();
00522 
00523   typename itk::Image<TType, 3>::SizeType ioSize = input->GetLargestPossibleRegion().GetSize();
00524 
00525   MITK_INFO << "size input image: " << ioSize[0] << ", " << ioSize[1] << ", " << ioSize[2];
00526   MITK_INFO << "size voxel: " << ioSize[0] * ioSize[1] * ioSize[2];
00527   
00528   int voxel= ioSize[0] * ioSize[1] * ioSize[2];
00529   
00530   unsigned char* RGBABuffer = new unsigned char[4*voxel];
00531   
00532   //for(int i=0;i<voxel;i++)
00533 
00534   // Convert volume  
00535   {
00536     ScalarBinaryToRGBAConverter<TType,BType> strc(inputData,input2Data,RGBABuffer,ioSize[0],ioSize[1],ioSize[2],tf);
00537     strc.fillSlices();
00538   }
00539     
00540   // Create MITK Image out of the raw data
00541   {
00542     mitk::Image::Pointer image = mitk::Image::New();
00543 
00544     unsigned int dimensions[ 3 ];
00545     dimensions[ 0 ] = ioSize[0];
00546     dimensions[ 1 ] = ioSize[1];
00547     dimensions[ 2 ] = ioSize[2];
00548 
00549     mitk::PixelType pixelType( typeid(RGBAPixel), 1, itk::ImageIOBase::RGBA );
00550     image->Initialize( pixelType, 3, dimensions );
00551     image->SetImportChannel( RGBABuffer, 0, Image::ManageMemory );
00552 
00553     return image;
00554   }
00555 }
00556 
00557  mitk::Image::Pointer mitkColourImageProcessor::convertWithBinaryToRGBAImage( mitk::Image::Pointer input1 ,mitk::Image::Pointer input2  , mitk::TransferFunction::Pointer tf )
00558 {
00559   MITK_INFO << "convertWithBinaryToRGBAImage" ;
00560     
00561   itk::Image< short, 3 >::Pointer inputCT;
00562   itk::Image< unsigned char, 3 >::Pointer inputBinary;
00563  
00564   if (*input1->GetPixelType().GetTypeId() == typeid(unsigned char) && *input2->GetPixelType().GetTypeId() == typeid(short))
00565   {
00566     mitk::CastToItkImage(input1,inputBinary);
00567     mitk::CastToItkImage(input2,inputCT);
00568   }
00569   else if (*input1->GetPixelType().GetTypeId() == typeid(short) && *input2->GetPixelType().GetTypeId() == typeid(unsigned char))
00570   {
00571     mitk::CastToItkImage(input1,inputCT);
00572     mitk::CastToItkImage(input2,inputBinary);
00573   }
00574   else
00575   {
00576     MITK_ERROR << "unsupported pixel type";
00577     return 0;
00578   }
00579   
00580   return ScalarAndBinaryToRGBA<short,unsigned char>(inputCT,inputBinary, tf);
00581 }
00582 
00584 
00585 
00586 
00587 
00588 template <class T,class B> 
00589 class ScalarBinaryColorToRGBAConverter
00590 {
00591   const T *dataPtr;
00592   const B *data2Ptr;
00593   unsigned char *tmpPtr;
00594   int sizeX;
00595   int sizeY;
00596   int sizeZ;
00597   int sizeXY;
00598   int sizeXm1;
00599   int sizeYm1;
00600   int sizeZm1;
00601   mitk::TransferFunction::Pointer tf;
00602   int *color;
00603   
00604   public:
00605 
00606   ScalarBinaryColorToRGBAConverter( const T *_dataPtr,const B *_data2Ptr,unsigned char *_tmpPtr,int _sizeX,int _sizeY,int _sizeZ,mitk::TransferFunction::Pointer _tf ,int *_color)
00607   {
00608     dataPtr=_dataPtr;
00609     data2Ptr=_data2Ptr;
00610     tmpPtr=_tmpPtr;
00611     sizeX=_sizeX;
00612     sizeY=_sizeY;
00613     sizeZ=_sizeZ;
00614 
00615     sizeXY=sizeX*sizeY;
00616     sizeXm1=sizeX-1;
00617     sizeYm1=sizeY-1;
00618     sizeZm1=sizeZ-1;
00619 
00620     tf=_tf;
00621 
00622     color = _color;
00623   }
00624 
00625   inline float sample(int x,int y,int z)
00626   {
00627     return float(dataPtr[ x + y * sizeX + z * sizeXY ]);
00628   }
00629   
00630   inline bool sampleBinary(int x,int y,int z)
00631   {
00632     return data2Ptr[ x + y * sizeX + z * sizeXY ];
00633   }
00634   
00635   inline int clamp(int x)
00636   {
00637     if(x<0) x=0; else if(x>255) x=255;
00638     return x;
00639   }
00640   
00641   inline void write(int x,int y,int z,float grayValue,float gx,float gy,float gz)
00642   {
00643     if(sampleBinary(x,y,z))
00644     {
00645       
00646    /*       
00647       gx /= aspect[0];
00648       gy /= aspect[1];
00649       gz /= aspect[2];
00650    */
00651       // Compute the gradient magnitude
00652 
00653       float t = sqrtf( gx*gx + gy*gy + gz*gz );
00654 
00655       int doff = x + y * sizeX + z * sizeXY;
00656       
00657       vtkPiecewiseFunction* opacityTransferFunction = tf->GetScalarOpacityFunction();
00658       vtkPiecewiseFunction* gradientTransferFunction = tf->GetGradientOpacityFunction();
00659       vtkColorTransferFunction* colorTransferFunction = tf->GetColorTransferFunction();
00660 
00661       double rgb[3];
00662       colorTransferFunction->GetColor( double(grayValue), rgb);
00663       
00664       double opacity= opacityTransferFunction->GetValue( double(grayValue) );
00665       
00666       opacity *= gradientTransferFunction->GetValue( double(0.5f*t) );
00667       
00668       tmpPtr[doff*4+0] = int( rgb[0]*255 + 0.5 );
00669       tmpPtr[doff*4+1] = int( rgb[1]*255 + 0.5 );
00670       tmpPtr[doff*4+2] = int( rgb[2]*255 + 0.5 );
00671       tmpPtr[doff*4+3] = int( opacity*255 + 0.5 );
00672 
00673       tmpPtr[doff*4+0] = color[0];
00674       tmpPtr[doff*4+1] = color[1];
00675       tmpPtr[doff*4+2] = color[2];
00676     }
00677     else
00678     {
00679       int doff = x + y * sizeX + z * sizeXY;
00680       tmpPtr[doff*4+0] = 0;
00681       tmpPtr[doff*4+1] = 0;
00682       tmpPtr[doff*4+2] = 0;
00683       tmpPtr[doff*4+3] = 0;
00684     }
00685   }
00686   
00687   
00688   inline void compute(int x,int y,int z)
00689   {
00690     float grayValue = sample(x,y,z);
00691     float gx,gy,gz;
00692 
00693     gx = sample(x+1,y,z) - sample(x-1,y,z);
00694     gy = sample(x,y+1,z) - sample(x,y-1,z);
00695     gz = sample(x,y,z+1) - sample(x,y,z-1);
00696     
00697     write( x, y, z, grayValue, gx, gy, gz );
00698  
00699   }
00700 
00701   inline void computeClamp(int x,int y,int z)
00702   {
00703     float grayValue = sample(x,y,z);
00704     float gx,gy,gz;
00705 
00706     if(x==0)            gx = 2.0f * ( sample(x+1,y,z) - grayValue       );
00707     else if(x==sizeXm1) gx = 2.0f * ( grayValue       - sample(x-1,y,z) );
00708     else                gx =          sample(x+1,y,z) - sample(x-1,y,z);
00709     
00710     if(y==0)            gy = 2.0f * ( sample(x,y+1,z) - grayValue       );
00711     else if(y==sizeYm1) gy = 2.0f * ( grayValue       - sample(x,y-1,z) );
00712     else                gy =          sample(x,y+1,z) - sample(x,y-1,z);
00713     
00714     if(z==0)            gz = 2.0f * ( sample(x,y,z+1) - grayValue       );
00715     else if(z==sizeZm1) gz = 2.0f * ( grayValue       - sample(x,y,z-1) );
00716     else                gz =          sample(x,y,z+1) - sample(x,y,z-1);
00717     
00718     write( x, y, z, grayValue, gx, gy, gz );
00719   }
00720 
00721   inline void compute1D(int y,int z)
00722   {
00723     int x;
00724     
00725     x=0;
00726     computeClamp(x,y,z);
00727     x++;
00728     
00729     while(x<sizeX-1)
00730     {
00731       compute(x,y,z);     
00732       x++;
00733     }
00734     
00735     if(x<sizeX)
00736     {
00737       computeClamp(x,y,z);
00738       x++;
00739     }
00740   }      
00741 
00742   inline void computeClamp1D(int y,int z)
00743   {
00744     int x;
00745     
00746     x=0;
00747     
00748     while(x<sizeX)
00749     {
00750       computeClamp(x,y,z);     
00751       x++;
00752     }
00753   }      
00754         
00755   inline void computeClamp2D(int z)
00756   {
00757     int y;
00758 
00759     y=0;
00760 
00761     while(y<sizeY)
00762     {
00763       computeClamp1D(y,z);     
00764       y++;
00765     }
00766   }
00767   
00768   inline void compute2D(int z)
00769   {
00770     int y;
00771 
00772     y=0;
00773     computeClamp1D(y,z);     
00774     y++;
00775     
00776     while(y<sizeY-1)
00777     {
00778       compute1D(y,z);     
00779       y++;
00780     }
00781     
00782     if(y<sizeY)
00783     {
00784       computeClamp1D(y,z);
00785       y++;
00786     }
00787   }
00788   
00789   inline void fillSlices()
00790   {
00791     int z;
00792 
00793     for(z=0;z<sizeZ;z++)
00794     {
00795       if(z==0 || z==sizeZ-1)
00796         computeClamp2D(z);
00797        else
00798         compute2D(z);
00799     }
00800   }    
00801 };
00802 
00803 
00804 template<class TType,class BType>   
00805 mitk::Image::Pointer mitkColourImageProcessor::ScalarAndBinaryAndColorToRGBA(itk::Image<TType, 3>* input  ,itk::Image<BType, 3>* input2  , mitk::TransferFunction::Pointer tf, int * color)
00806 {
00807   const TType *inputData=input->GetBufferPointer();
00808   const BType *input2Data=input2->GetBufferPointer();
00809 
00810   typename itk::Image<TType, 3>::SizeType ioSize = input->GetLargestPossibleRegion().GetSize();
00811 
00812   MITK_INFO << "size input image: " << ioSize[0] << ", " << ioSize[1] << ", " << ioSize[2];
00813   MITK_INFO << "size voxel: " << ioSize[0] * ioSize[1] * ioSize[2];
00814   
00815   int voxel= ioSize[0] * ioSize[1] * ioSize[2];
00816   
00817   unsigned char* RGBABuffer = new unsigned char[4*voxel];
00818   
00819   //for(int i=0;i<voxel;i++)
00820 
00821   // Convert volume  
00822   {
00823     ScalarBinaryColorToRGBAConverter<TType,BType> strc(inputData,input2Data,RGBABuffer,ioSize[0],ioSize[1],ioSize[2],tf,color);
00824     strc.fillSlices();
00825   }
00826     
00827   // Create MITK Image out of the raw data
00828   {
00829     mitk::Image::Pointer image = mitk::Image::New();
00830 
00831     unsigned int dimensions[ 3 ];
00832     dimensions[ 0 ] = ioSize[0];
00833     dimensions[ 1 ] = ioSize[1];
00834     dimensions[ 2 ] = ioSize[2];
00835 
00836     mitk::PixelType pixelType( typeid(RGBAPixel), 1, itk::ImageIOBase::RGBA );
00837     image->Initialize( pixelType, 3, dimensions );
00838     image->SetImportChannel( RGBABuffer, 0, Image::ManageMemory );
00839 
00840     return image;
00841   }
00842 }
00843 
00844  mitk::Image::Pointer mitkColourImageProcessor::convertWithBinaryAndColorToRGBAImage( mitk::Image::Pointer input1 ,mitk::Image::Pointer input2  , mitk::TransferFunction::Pointer tf , int * color)
00845 {
00846   MITK_INFO << "convertWithBinaryToRGBAImage" ;
00847     
00848   itk::Image< short, 3 >::Pointer inputCT;
00849   itk::Image< unsigned char, 3 >::Pointer inputBinary;
00850  
00851   if (*input1->GetPixelType().GetTypeId() == typeid(unsigned char) && *input2->GetPixelType().GetTypeId() == typeid(short))
00852   {
00853     mitk::CastToItkImage(input1,inputBinary);
00854     mitk::CastToItkImage(input2,inputCT);
00855   }
00856   else if (*input1->GetPixelType().GetTypeId() == typeid(short) && *input2->GetPixelType().GetTypeId() == typeid(unsigned char))
00857   {
00858     mitk::CastToItkImage(input1,inputCT);
00859     mitk::CastToItkImage(input2,inputBinary);
00860   }
00861   else
00862   {
00863     MITK_ERROR << "unsupported pixel type";
00864     return 0;
00865   }
00866   
00867   return ScalarAndBinaryAndColorToRGBA<short,unsigned char>(inputCT,inputBinary, tf,color);
00868 }
00869 
00870 static  inline int clamp(int x)
00871   {
00872     if(x<0) x=0; else if(x>255) x=255;
00873     return x;
00874   }
00875   
00876 
00877 mitk::Image::Pointer mitkColourImageProcessor::CombineRGBAImage( unsigned char* input  ,unsigned char* input2, int sizeX,int sizeY,int sizeZ )
00878 { 
00879   int voxel= sizeX*sizeY*sizeZ;
00880   
00881   unsigned char* RGBABuffer = new unsigned char[4*voxel];
00882   
00883   // Convert volume  
00884   {
00885     for( int r=0;r<voxel;r++)
00886     {
00887       int rgbaInput1[4];
00888       int rgbaInput2[4];
00889       
00890       rgbaInput1[0]=input[r*4+0];
00891       rgbaInput1[1]=input[r*4+1];
00892       rgbaInput1[2]=input[r*4+2];
00893       rgbaInput1[3]=input[r*4+3];
00894       
00895       rgbaInput2[0]=input2[r*4+0];
00896       rgbaInput2[1]=input2[r*4+1];
00897       rgbaInput2[2]=input2[r*4+2];
00898       rgbaInput2[3]=input2[r*4+3];
00899 
00900       int result[4];
00901 
00902       /*
00903       
00904       float fac1 = rgbaInput1[3]/255.0f;
00905       float fac2 = rgbaInput2[3]/255.0f;
00906 
00907       float sum = fac1+fac2;
00908       
00909       float cor = 1.0f;
00910       
00911       if(sum > 1.0f)
00912         cor = 1.0f/sum;
00913         
00914       fac1 *= cor;
00915       fac2 *= cor;
00916       
00917       result[0]= clamp(int(fac1 * rgbaInput1[0] + fac2 * rgbaInput2[0] + 0.5f));
00918       result[1]= clamp(int(fac1 * rgbaInput1[1] + fac2 * rgbaInput2[1] + 0.5f));
00919       result[2]= clamp(int(fac1 * rgbaInput1[2] + fac2 * rgbaInput2[2] + 0.5f));
00920       result[3]= clamp(int(fac1 * rgbaInput1[3] + fac2 * rgbaInput2[3] + 0.5f));
00921       */
00922      
00923       if( rgbaInput1[3] )
00924       {
00925         result[0]= rgbaInput1[0];
00926         result[1]= rgbaInput1[1];
00927         result[2]= rgbaInput1[2];
00928         result[3]= rgbaInput1[3];      
00929       }
00930       else
00931       {
00932         result[0]= rgbaInput2[0];
00933         result[1]= rgbaInput2[1];
00934         result[2]= rgbaInput2[2];
00935         result[3]= rgbaInput2[3];
00936       }
00937       
00938       
00939       RGBABuffer[r*4+0]= result[0];
00940       RGBABuffer[r*4+1]= result[1];
00941       RGBABuffer[r*4+2]= result[2];
00942       RGBABuffer[r*4+3]= result[3];
00943     }
00944     
00945   }
00946   
00947   // Create MITK Image out of the raw data
00948   {
00949     mitk::Image::Pointer image = mitk::Image::New();
00950 
00951     unsigned int dimensions[ 3 ];
00952     dimensions[ 0 ] = sizeX;
00953     dimensions[ 1 ] = sizeY;
00954     dimensions[ 2 ] = sizeZ;
00955 
00956     mitk::PixelType pixelType( typeid(RGBAPixel), 1, itk::ImageIOBase::RGBA );
00957     image->Initialize( pixelType, 3, dimensions );
00958     image->SetImportChannel( RGBABuffer, 0, Image::ManageMemory );
00959 
00960     return image;
00961   }
00962 }
00963 
00964 mitk::Image::Pointer mitkColourImageProcessor::combineRGBAImage( mitk::Image::Pointer input1 , mitk::Image::Pointer input2)
00965 {
00966   
00967   RGBAImage::Pointer itk1,itk2;
00968 
00969   unsigned char *data1=(unsigned char *)input1->GetData();
00970   unsigned char *data2=(unsigned char *)input2->GetData();
00971 
00972   unsigned int *dim = input1->GetDimensions();
00973    
00974   return CombineRGBAImage(data1,data2,dim[0],dim[1],dim[2]);
00975 }
00976 
00977 
00978 }//end namespace mitk
00979 
00980