[Insight-users] co-occurrence matrix

[Glenn Pierce]

--=-y5KyOvD2dlTzj5owqB0i
Content-Type: text/plain
Content-Transfer-Encoding: 7bit

Ok

I have a first version of a co occurrence filter.
I sure there are many things I can improve.

The main problem is I have stored each pair in a ListSample (each pair
as a new MeasurementVector).
However, I need to increase the frequency of duplicate pairs. I wasn't
sure how to efficiently find if a pair have already been stored and then
increment the frequency of that pair.
see line 76 of itkImageToCooccurenceList.txx


I had to use quite a few typename constructs to remove compiler warnings
of 'implicit typename' if anyone wonders about them.


Also I inherited from ImageToListAdaptor
but I am storing values in a ListSample that is a member of this
class. This is different from the other adaptors as it is storing
values. I don't know whether this is suitable for you?

Any suggestions of improving this filter. Even better itk method names
etc would be great.

I eventually am going to use this to produced an co-occurrence image in
vtk. Does anyone know how you can produce a grey level image with
different intensity values for the pixels in vtk? I probably should ask
on the vtk list -sorry.


I have attached the filter and a small driver program.

Thanks for the help.


Glenn


On Fri, 2004-05-14 at 11:12 -0400, Luis Ibanez wrote:
> Hi Glenn,
> 
> Thanks for the clarification.
> 
> Well...
> 
>      You probably will have to write that
>      ImageToListAdaptor class after alll    :-)
> 
> 
> I would suggest you to write it as a
> 
>          ImageToListFilter
> 
> 
> and take advantage of the Neighborhood
> Iterators.  You probably can get very
> close by looking at the code used for
> the Morphological filters.
> 
> Please read the Iterators chapter on
> the SoftwareGuide. That will help you
> find your way.
> 
> It will be nice let the option to the
> user for defining what neighbors to
> use as cliques for the coorcurrence
> matrix.
> 
> Your filter will populate a List of
> 2D vector values having all the pairs
> of cliques from your images.
> 
> 
> Please let us know if you need any
> help implementing this filter.
> 
> and... of course, we will be happy
> to get it back into the toolkit
> once you are done   :-)
> 



--=-y5KyOvD2dlTzj5owqB0i
Content-Disposition: attachment; filename=itkImageToCooccurenceList.h
Content-Type: text/x-chdr; name=itkImageToCooccurenceList.h; charset=iso-8859-1
Content-Transfer-Encoding: 7bit

#ifndef __itkImageToCooccurenceList_h
#define __itkImageToCooccurenceList_h

#include <typeinfo>

#include "itkImage.h"
#include "itkPixelTraits.h"
#include "itkImageToListAdaptor.h"
#include "itkSmartPointer.h"
#include "itkImageRegionIterator.h"
#include "itkShapedNeighborhoodIterator.h"
#include "itkNeighborhoodIterator.h"
#include "itkNeighborhoodAlgorithm.h"
#include "itkConstantBoundaryCondition.h"
#include "itkListSample.h"
#include "itkFixedArray.h"
#include "itkMacro.h"

#include <vector>
#include <algorithm>
#include <iostream>

namespace itk{ 
namespace Statistics{

template < class TImage >
class ITK_EXPORT ImageToCooccurenceList
  : public ImageToListAdaptor< 
  TImage, 
  FixedArray< typename TImage::PixelType, 2 > >
{
public:
  typedef TImage ImageType;

  typedef FixedArray< typename TImage::PixelType, 2 > MeasurementVectorType ;

  typedef typename itk::Statistics::ListSample< MeasurementVectorType > SampleType ;

  /** Standard class typedefs */
  typedef ImageToCooccurenceList Self;
  typedef ImageToListAdaptor< TImage, MeasurementVectorType > Superclass;
  typedef SmartPointer< Self > Pointer;
  typedef SmartPointer<const Self> ConstPointer;

  /** Neighborhood iterator type. */
  typedef itk::ShapedNeighborhoodIterator< TImage , ConstantBoundaryCondition<TImage> > ShapedNeighborhoodIteratorType;

  /** Offset type used for Neighborhoods **/
  typedef typename ShapedNeighborhoodIteratorType::OffsetType OffsetType;
  typedef typename std::vector<OffsetType> OffsetTable;
  
  void UseNeighbor(OffsetType offset);

  /** Triggers the Computation of the co-occurence matrix */
  void Compute();
 
  /** Run-time type information (and related methods). */
  itkTypeMacro(ImageToCooccurenceList, ListSampleBase);
  
  /** Method for creation through the object factory. */
  itkNewMacro(Self);
  
  /** the number of components in a measurement vector */
  itkStaticConstMacro(MeasurementVectorSize, unsigned int, 2);

  /** Superclass typedefs for Measurement vector, measurement, 
   * Instance Identifier, frequency, size, size element value */
  typedef typename Superclass::FrequencyType FrequencyType ;
  typedef typename Superclass::MeasurementType MeasurementType ;
  typedef typename Superclass::InstanceIdentifier InstanceIdentifier ;
  typedef MeasurementVectorType ValueType ;

protected:
  ImageToCooccurenceList();
  virtual ~ImageToCooccurenceList() {}
  void PrintSelf(std::ostream& os, Indent indent) const;  

private:
  ImageToCooccurenceList(const Self&) ; //purposely not implemented
  void operator=(const Self&) ; //purposely not implemented
  OffsetTable m_OffsetTable;
  typename SampleType::Pointer sample;
} ; // end of class ScalarImageToListAdaptor

} // end of namespace Statistics
} // end of namespace itk

#ifndef ITK_MANUAL_INSTANTIATION
#include "itkImageToCooccurenceList.txx"
#endif

#endif

--=-y5KyOvD2dlTzj5owqB0i
Content-Disposition: attachment; filename=itkImageToCooccurenceList.txx
Content-Type: text/plain; name=itkImageToCooccurenceList.txx; charset=iso-8859-1
Content-Transfer-Encoding: 7bit

#ifndef _itkImageToCooccurenceList_txx
#define _itkImageToCooccurenceList_txx

namespace itk{ 
namespace Statistics{

template < class TImage >
ImageToCooccurenceList< TImage >
::ImageToCooccurenceList()
{
  sample = SampleType::New();
}

template < class TImage >
void
ImageToCooccurenceList< TImage >
::PrintSelf(std::ostream& os, Indent indent) const
{
  Superclass::PrintSelf(os,indent);
}

template < class TImage >
void
ImageToCooccurenceList< TImage >
::Compute()
{
  typename ShapedNeighborhoodIteratorType::RadiusType radius;
  radius.Fill(1);

  itk::ConstantBoundaryCondition<TImage> boundaryCondition;
  // 0 is valid so I chose -1. Is this valid for all images ?
  boundaryCondition.SetConstant( -1 );

  typedef itk::NeighborhoodAlgorithm::ImageBoundaryFacesCalculator<
                                                ImageType > FaceCalculatorType;
                                                                                                                            
  FaceCalculatorType faceCalculator;
  typename FaceCalculatorType::FaceListType faceList;
  typename FaceCalculatorType::FaceListType::iterator fit;
  typename ShapedNeighborhoodIteratorType::ConstIterator ci;    
  ShapedNeighborhoodIteratorType it;
  typename ImageType::IndexType index;
  typename ImageType::PixelType pixel_intensity, center_pixel_intensity;
  typename OffsetTable::iterator iter;
  typename SampleType::MeasurementVectorType coords;  
                                                                                                                      
  faceList = faceCalculator( GetImage(),
                             GetImage()->GetRequestedRegion(),
                             radius );

  OffsetType center_offset = {{0,0}};

  for ( fit=faceList.begin(); fit != faceList.end(); ++fit) {

    it = ShapedNeighborhoodIteratorType( radius, GetImage(), *fit );
    it.OverrideBoundaryCondition(&boundaryCondition);

    for (iter = m_OffsetTable.begin(); iter != m_OffsetTable.end(); iter++) {
      it.ActivateOffset(*iter);
    }
                                                                                                                            
    for (it.GoToBegin(); !it.IsAtEnd(); ++it) {

        center_pixel_intensity = it.GetPixel(center_offset);

        for (ci = it.Begin(); ci != it.End(); ci++) {

            pixel_intensity = ci.Get();
            std::cout << "pixel value: " << pixel_intensity << "  " << std::endl;

            // We have the intensity values for the center pixel and one of it's neighbours.
            // We can now place these in the SampleList
            coords[0] = center_pixel_intensity;
            coords[1] = pixel_intensity;

            sample->PushBack(coords) ;

        }
    }
   }

}


template < class TImage >
void
ImageToCooccurenceList< TImage >
::UseNeighbor(OffsetType offset)
{
  // Don't add the center pixel
  if(offset[0] == 0 && offset[1] == 0) {
      //std::cout << "HERE " << std::endl;
      return;
  }

  m_OffsetTable.push_back(offset);
}


} // end of namespace Statistics 
} // end of namespace itk

#endif

--=-y5KyOvD2dlTzj5owqB0i
Content-Disposition: attachment; filename=Test.cxx
Content-Type: text/x-c++src; name=Test.cxx; charset=iso-8859-1
Content-Transfer-Encoding: 7bit

#include "itkImage.h"

#include "itkImageFileReader.h"
#include "itkImageToCooccurenceList.h"

int main( int argc, char * argv [] )
{
  typedef float      PixelType;
  const   unsigned int        Dimension = 2;
  typedef itk::Image< PixelType, Dimension >    ImageType;
  typedef itk::ImageFileReader< ImageType >  ReaderType;
  typedef itk::Statistics::ImageToCooccurenceList < ImageType > CooccurenceListType;

  ReaderType::Pointer reader = ReaderType::New();
                                                                                                                            
  reader->SetFileName( "/home/glenn/development/histogram/test.png" );
                                                                                                                            
  try {
      reader->Update();
  }
  catch( itk::ExceptionObject & exp )
  {
      std::cerr << "Exception caught" << std::endl;
      std::cerr << exp << std::endl;
  }

  CooccurenceListType::Pointer list = CooccurenceListType::New();
  list->SetImage(reader->GetOutput());

  CooccurenceListType::OffsetType offset0 = {{-1,0}};
  CooccurenceListType::OffsetType offset1 = {{1,0}};
  CooccurenceListType::OffsetType offset2 = {{0,-1}};
  CooccurenceListType::OffsetType offset3 = {{0,1}};

  list->UseNeighbor(offset0);
  list->UseNeighbor(offset1);
  list->UseNeighbor(offset2);
  list->UseNeighbor(offset3);

  list->Compute();

  return 0;
  
}



--=-y5KyOvD2dlTzj5owqB0i--