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

[Zachary Pincus]

Hi Sachin (and any others interested in texture calculations),

I just checked into the ITK CVS a modification to the 
ScalarImageTextureCalculator class that should allow for approximately 
four-fold faster texture calculations.

By default, texure features are computed for each spatial direction and 
then averaged afterward, so it is possible to access the standard 
deviations of the texture features. These values give a clue as to 
texture anisotropy. However, doing this is much more work, because it 
involved computing one GLCM for each offset given. To compute a single 
GLCM for all of the offsets, call FastCalculationsOn(). If this is 
called, then the texture standard deviations will not be computed (and 
will be set to zero), but texture computation will be much faster.

So, if you don't use the texture standard deviations, just set 
FastCalculationsOn() and you should get similar (not identical; the 
procedure to ensure rotation invariance is different, sort of like the 
difference between the average of a sum and the sum of averages!) 
texture results faster.

Zach Pincus

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


On Mar 3, 2005, at 11:54 PM, Zachary Pincus wrote:

> I think that the glcm calculations are inherently expensive: they 
> require several passes through a 2D histogram in order to calculate 
> various features.
>
> I don't know of any particular shortcuts around this. My suggestion 
> would be to run a profiler on the code to see what the slow parts are. 
> (Is it scanning across the image? Iterating through the histogram? 
> Actually computing the features?) Once you know the slow parts, you 
> can try to speed them up.
>
> For example, you could modify the texture coefficient calculator class 
> to only calculate the single texture value you require, if that value 
> needs fewer passes through the histogram than the calculator currently 
> takes. Or modifying the code to use a sparse frequency container for 
> the histogram might be a good idea.
>
> Anyhow, you're basically in optimization-land at this point, so all 
> the general advice applies. I don't know of any particularly better 
> way to compute the GLCM values specifically.
>
> You might also think about down-sampling your images first. Do you 
> really need texture calculated at *every* pixel? Would interpolation 
> work instead? If not, then there can be some savings there. It's worth 
> thinking about alternate solutions.
>
> Zach
>
>
>
> On Mar 3, 2005, at 11:09 PM, Sachin Jambawalikar wrote:
>
>> Hi  Zach,
>> thanks  for your prompt reply.
>> I made the necessary changes you asked for. But still the  Glcm takes
>> about  360 second
>>  for one slice with  the radius of neighboorhood 12. I was wondering
>> if you  know of a faster implementation of glcm or have come across
>> some references  for fast glcm matrix computation, I had implemented a
>>  C code for (without using the itk library) glcm maps computation
>> which  computed   glcm maps for window 5x5  and histogram bins 256 in
>> about  780 second ,so I feel that the glcm computation with  your
>> class  is slightly faster
>> or round about the same ( for your class  I am using radius 12 ie
>> window of 25x25). Is there any scope  for improving the speed .
>>
>> I'm trying to compute slice by slice glcm feature maps for a 3D
>> volume.And my volumes are 512x512x50 so it would take me like 5 hours
>> for a volume.
>>
>> this is the   out put and the changed code:
>>           Probe Tag    Starts    Stops           Time
>>           Slice glcm           1            1           360.968
>> computed0thslice
>>           Probe Tag    Starts    Stops           Time
>>           Slice glcm           1            1           364.125
>> computed1thslice
>>           Probe Tag    Starts    Stops           Time
>>           Slice glcm           1            1           370.906
>> computed2thslice
>>           Probe Tag    Starts    Stops           Time
>>           Slice glcm           1            1           371.359
>> computed3thslice
>>
>>
>>
>>
>> ComputeNeighboorhoodConvolution(ImagePointer2D inimage)
>> 	{
>> 	
>> //typedef  PixelType short;
>>
>> 	 itk::TimeProbesCollectorBase collector;
>>
>> typedef itk::Image<double,2> InternalPixelType;
>> InternalPixelType::Pointer output = InternalPixelType::New();
>> typedef itk::ImageRegionIterator< InternalPixelType>        
>> Iterator2DType;
>>
>> output->SetRegions(inimage->GetRequestedRegion());
>> output->Allocate();
>> output->FillBuffer(0.0);
>>
>> ImageType2D::Pointer BBOut = ImageType2D::New();
>> BBOut=inimage;
>>
>> NeighborhoodIteratorType::RadiusType radius;
>> 	
>>
>> itk::NeighborhoodInnerProduct<ImageType2D> innerProduct;
>> 	radius.Fill(12);
>> 	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;
>> double result=0.0;
>>
>> ROIFilterType::Pointer roifilter=ROIFilterType::New();
>>  	roifilter->SetInput(inimage);
>> 			fit=faceList.begin();
>> 		it = NeighborhoodIteratorType( radius,inimage, *fit );
>> 		out = Iterator2DType( output, *fit );
>>    collector.Start("Slice glcm");
>>  texturecal->SetNumberOfBinsPerAxis(16);	
>> 	 for (it.GoToBegin(), out.GoToBegin(); ! it.IsAtEnd(); ++it, ++out)
>> 			{
>> 		
>>    	BBOut->SetRequestedRegion(it.GetBoundingBoxAsImageRegion());
>> 	     texturecal->SetInput(BBOut);
>> 			  texturecal->Compute();
>> 		 	FeatureValueVectorPointer tmp= texturecal->GetFeatureMeans();
>> 				result=static_cast<double>(tmp->GetElement(0));
>> 				out.Set(result);
>> 			}
>> 	
>> 		collector.Stop("Slice glcm");
>> 		collector.Report();
>>
>> 		
>> 	return inimage;
>> 	}
>>
>>
>>
>>
>> Regards
>> --Sachin
>>
>