[Insight-users] Generating coocurence feature maps for 2D images

[Zachary Pincus]

Hi,

First, you should note that the texture calculator class uses only the 
requested region set in the image for its calculations. So instead of 
mucking about by dividing up the input image into regions with the ROI 
filter, just set the input image's requested region to the appropriate 
region.

Second, note that you only want the first entry in the face list : 
that's the region where the neighborhood won't go out of bounds. Other 
regions in the face list are regions where the neighborhood will be 
partially out of bounds. It is also *precisely* these regions which 
cause the ROI filter to throw exceptions. Recall that throwing 
exceptions is VERY slow. For EACH pixel in EACH face list entry 
(excluding the first face), you cause an exception to be thrown. The 
moral of this story is to only iterate over the region in the first 
face. (In fact, why did you use the face list calculator in the first 
place, if you were just going to iterate over all the regions it gave 
you?)

Third, note that while you ask the texture calculator to calculate all 
of the texture values, you only use one value. It would be faster to 
explicitly request only the texture feature values you actually plan on 
using.

Fourth, note that you have not set the number of histogram bins for the 
coocurrence histogram. This number defaults to 256. Calculating 
features over a 256x256 coocurrence matrix is slow. Try setting the 
number of bins to something like 16 or 32. Alternately, if you must use 
256 bins, you would do well to make your own modification to the 
texture calculator class so that it could optionally use a histogram 
with a sparse frequency container instead of a dense container. Doing 
the texture calculations on a sparse histogram would be much faster. 
(I've been planning on making these changes for a while; I just haven't 
gotten around to it.)

So, I suspect that if you don't use the ROI filter, only iterate over 
the first entry in the face list, only request that the calculator 
compute the feature you're interested in, and tun down the number of 
bins in the historgram, you will get much faster computation.

Zach Pincus

Department of Biochemistry and Program in Biomedical Informatics
Stanford University School of Medicine

On Mar 3, 2005, at 4:13 PM, Sachin Jambawalikar wrote:

> Hi all,
>
> Especially Luis, Zachary Pincus and Glenn Pierce
>
> I'm trying to use the  itkScalarImageTextureCalculator class
>  to compute texture features for an image slice.
> However instead of finding  one mean  feature value(energy entropy
> etc) for the 2D image I  am trying to find  the same for a 5x5 window
> . Similar to convolution but here the glcm and the features are
> computed for  each of the windowed subimages and the center pixel is
> replaced by the feature value. However this  seems to be
> computationally expensive  for my  image of 512*512  . Can any one
> suggest a faster way of computing the cooccurence feature maps or tell
> me what I'm doing wrong in my code
>
>
> ImagePointer2D ComputeNeighboorhoodConvolution(ImagePointer2D inimage)
> {
> typedef itk::Image<double,2> InternalPixelType;
> InternalPixelType::Pointer output = InternalPixelType::New();
> typedef itk::ImageRegionIterator< InternalPixelType>        
> Iterator2DType;
>
> output->SetRegions(inimage->GetRequestedRegion());
> output->Allocate();
>
> ImageType2D::Pointer BBOut = ImageType2D::New();
>
> NeighborhoodIteratorType::RadiusType radius;
> radius.Fill(2);
> typedef itk::NeighborhoodAlgorithm::ImageBoundaryFacesCalculator<
> ImageType2D> FaceCalculatorType;
> 	FaceCalculatorType faceCalculator;
> 	FaceCalculatorType::FaceListType faceList;
> 	faceList = faceCalculator(inimage, output->GetRequestedRegion(),
> 		radius);
> 	FaceCalculatorType::FaceListType::iterator fit;
>
> 	Iterator2DType out;
>   NeighborhoodIteratorType it;
> 	
> typedef itk::Statistics::ScalarImageTextureCalculator<ImageType2D>
> ScalarImageTextureCalculatorType;
> ScalarImageTextureCalculatorType::Pointer
> texturecal=ScalarImageTextureCalculatorType::New();
> typedef itk::RegionOfInterestImageFilter< ImageType2D,ImageType2D >
> ROIFilterType;
>
> typedef VectorContainer<unsigned char, double>      FeatureValueVector;
> typedef typename FeatureValueVector::Pointer        
> FeatureValueVectorPointer;
>
> ROIFilterType::Pointer roifilter=ROIFilterType::New();
>
> roifilter->SetInput(inimage);
> 	for ( fit=faceList.begin(); fit != faceList.end(); ++fit)
> 		{
> 		it = NeighborhoodIteratorType( radius,inimage, *fit );
> 		out = Iterator2DType( output, *fit );
> 		for (it.GoToBegin(), out.GoToBegin(); ! it.IsAtEnd(); ++it, ++out)
> 			{
> 		
> 		
> double result=0.0;
> bool BoundryFlag=false;			
>                                 try
> 				{
> 	
> 				roifilter->SetRegionOfInterest(it.GetBoundingBoxAsImageRegion());
> 				roifilter->Update();
> 				BBOut=roifilter->GetOutput();
>                                }
> 			catch( itk::ExceptionObject & err )
> 				{
> 				std::cout << "ExceptionObject caught !" << std::endl;
> 				std::cout << err << std::endl;
> 				BoundryFlag=true;
> 			        result=0.0;
> 				}
> 			
> 			if(BoundryFlag==false)
> 				{	
>                                 texturecal->SetInput(BBOut);
> 				 texturecal->SetNumberOfBinsPerAxis(16);
> 			        texturecal->Compute();
> 				FeatureValueVectorPointer tmp= texturecal->GetFeatureMeans();
> 				result=static_cast<double>(tmp->GetElement(0));
> 				}
> 		
> 			out.Set(result);
> 			}
> 	
> 			std::cout<<"faceList"<< "done out of="<<faceList.size()<<std::endl;
> 		}
> 		
> 	return inimage;// for jus checking time taken  in reality outimgae
> should be recasted
>                               //  and returned
> 	}
>
>
>
> Regards
> --Sachin
> _______________________________________________
> Insight-users mailing list
> Insight-users at itk.org
> http://www.itk.org/mailman/listinfo/insight-users
>