[ITK-users] Histogram for a double precision array

[SDobashi]

Hello, everyone,

I want to calculate a histogram for a double precision C-array (or
std::vector<double>).

As a simple example, I tried to calculate a histogram for the data

double data[] = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9};

for the range from 5 to 10 in 10 bins, and I wrote the code listed in the
last part of this post.

This code worked well for "unsined char" but not for "double".

Could anyone suggest me why this could happen?
Please tell me a proper way to calculate a histogram of the double precision
array.

Thanks a lot in advance,.

Kind regards,

# Results for "unsigned char" data and this is what I expected:

Frequency = [ 1,
0,
1,
0,
1,
0,
1,
0,
1,
0 ]
Frequency of 0 : (5 to 5.5) = 1
Frequency of 1 : (5.5 to 6) = 0
Frequency of 2 : (6 to 6.5) = 1
Frequency of 3 : (6.5 to 7) = 0
Frequency of 4 : (7 to 7.5) = 1
Frequency of 5 : (7.5 to 8) = 0
Frequency of 6 : (8 to 8.5) = 1
Frequency of 7 : (8.5 to 9) = 0
Frequency of 8 : (9 to 9.5) = 1
Frequency of 9 : (9.5 to 10) = 0
Total count 5

# Results for "double" precision data:

Frequency = [ 1,
1,
1,
1,
1,
1,
1,
1,
1,
1 ]
Frequency of 0 : (0 to 0.9009) = 1
Frequency of 1 : (0.9009 to 1.8018) = 1
Frequency of 2 : (1.8018 to 2.7027) = 1
Frequency of 3 : (2.7027 to 3.6036) = 1
Frequency of 4 : (3.6036 to 4.5045) = 1
Frequency of 5 : (4.5045 to 5.4054) = 1
Frequency of 6 : (5.4054 to 6.3063) = 1
Frequency of 7 : (6.3063 to 7.2072) = 1
Frequency of 8 : (7.2072 to 8.1081) = 1
Frequency of 9 : (8.1081 to 9.009) = 1
Total count 10


# Code I wrote:

#include "itkImage.h"
#include "itkImportImageFilter.h"
#include "itkImageToHistogramFilter.h"
 
using namespace std;
 
int main(int argc, char* argv[])
{
    //typedef unsigned char PixelType;
    typedef double PixelType;
    const int Dimension = 1;
 
    const int numElems = 10;
    PixelType data[numElems] = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9};
 
    typedef itk::Image<PixelType, Dimension> ImageType;
    typedef itk::ImportImageFilter<PixelType, Dimension> ImportFilterType;
 
    ImportFilterType::Pointer importFilter = ImportFilterType::New();
 
    ImportFilterType::SizeType size;
    size[0] = numElems;
 
    ImportFilterType::IndexType start;
    start[0] = 0;
 
    ImportFilterType::RegionType region;
    region.SetIndex(start);
    region.SetSize(size);
 
    importFilter->SetRegion(region);
 
    importFilter->SetImportPointer(data, numElems, false);
 
    importFilter->Update();
 
    ImageType::Pointer img = importFilter->GetOutput();
 
    const unsigned int binsPerDimension = numElems;
 
    const unsigned int MeasurementVectorSize = 1;
 
    typedef itk::Statistics::ImageToHistogramFilter<ImageType>
ImageToHistogramFilterType;
 
    ImageToHistogramFilterType::HistogramType::MeasurementVectorType
lowerBound(binsPerDimension);
    lowerBound.Fill(5);
 
    ImageToHistogramFilterType::HistogramType::MeasurementVectorType
upperBound(binsPerDimension);
    upperBound.Fill(10);
 
    ImageToHistogramFilterType::HistogramType::SizeType
histSize(MeasurementVectorSize);
    histSize.Fill(binsPerDimension);
 
    ImageToHistogramFilterType::Pointer imageToHistogramFilter =
ImageToHistogramFilterType::New();
    imageToHistogramFilter->SetInput(img);
    imageToHistogramFilter->SetHistogramBinMinimum(lowerBound);
    imageToHistogramFilter->SetHistogramBinMaximum(upperBound);
    imageToHistogramFilter->SetHistogramSize(histSize);
 
    try {
        imageToHistogramFilter->Update();
    }
    catch (itk::ExceptionObject& error) {
        std::cerr << "Error: " << error << std::endl;
        return EXIT_FAILURE;
    }
 
    ImageToHistogramFilterType::HistogramType* histogram =
imageToHistogramFilter->GetOutput();
 
    std::cout << "Frequency = [ ";
    for (unsigned int i = 0; i < histogram->GetSize()[0]; ++i) {
        std::cout << histogram->GetFrequency(i);
 
        if (i != histogram->GetSize()[0] - 1) {
            std::cout << "," << std::endl;
        }
    }
    std::cout << " ]" << std::endl;
 
    for (unsigned int i = 0; i < histogram->GetSize()[0]; i++) {
        std::cout << "Frequency of " << i << " : ("
        << histogram->GetBinMin(0, i) << " to " << histogram->GetBinMax(0,
i)
        << ") = " << histogram->GetFrequency(i) << std::endl;
    }
 
    std::cout << "Total count " << histogram->GetTotalFrequency() <<
std::endl;
}




--
View this message in context: http://itk-insight-users.2283740.n2.nabble.com/Histogram-for-a-double-precision-array-tp7589566.html
Sent from the ITK Insight Users mailing list archive at Nabble.com.