I'm trying to convert the code for using the maximum entropy thresholding from this matlab code:
%**************************************************************************
%**************************************************************************
%
% maxentropie is a function for thresholding using Maximum Entropy
%
%
% input = I ==> Image in gray level
% output =
% I1 ==> binary image
% threshold ==> the threshold choosen by maxentropie
%
% F.Gargouri
%
%
%**************************************************************************
%**************************************************************************
function [threshold I1]=maxentropie(I)
[n,m]=size(I);
h=imhist(I);
%normalize the histogram ==> hn(k)=h(k)/(n*m) ==> k in [1 256]
hn=h/(n*m);
%Cumulative distribution function
c(1) = hn(1);
for l=2:256
c(l)=c(l-1)+hn(l);
end
hl = zeros(1,256);
hh = zeros(1,256);
for t=1:256
%low range entropy
cl=double(c(t));
if cl>0
for i=1:t
if hn(i)>0
hl(t) = hl(t)- (hn(i)/cl)*log(hn(i)/cl);
end
end
end
%high range entropy
ch=double(1.0-cl); %constraint cl+ch=1
if ch>0
for i=t+1:256
if hn(i)>0
hh(t) = hh(t)- (hn(i)/ch)*log(hn(i)/ch);
end
end
end
end
% choose best threshold
h_max =hl(1)+hh(1)
threshold = 0;
entropie(1)=h_max;
for t=2:256
entropie(t)=hl(t)+hh(t);
if entropie(t)>h_max
h_max=entropie(t);
threshold=t-1;
end
end
% Display
I1 = zeros(size(I));
I1(I<threshold) = 0;
I1(I>threshold) = 255;
%imshow(I1)
end
The problem is that I'm getting floating point excpetion error in the code, and I cannot understand why
This is my implementation:
#include <iostream>
#include <opencv2/highgui/highgui.hpp>
#include <opencv2/imgproc/imgproc.hpp>
#include <math.h>
using namespace cv;
using namespace std;
int main(){
cout.setf(std::ios_base::fixed, std::ios_base::floatfield);
cout.precision(9);
Mat old_image=imread("2.png",CV_LOAD_IMAGE_GRAYSCALE);
double minval, maxval;
minMaxLoc(old_image,&minval, &maxval);
cout<<minval<<" "<<maxval<<endl;
Mat image;
old_image.convertTo(image, CV_8UC1, 255.0/(maxval-minval), -minval*255.0/(maxval-minval));
minMaxLoc(image,&minval, &maxval);
cout<<minval<<" "<<maxval;
int k=0;
imshow("im",image);
waitKey(0);
for(int y=0; y<image.rows;y++){
for(int x=0; x<image.cols;x++){
if((int) image.at<uchar>(y,x)==0){
k++;
}
}
}
cout<<k<<endl<<endl;
int i, l, j, t;
int histSize = 256;
float range[] = { 0, 255 };
const float *ranges[] = { range };
Mat hist, histogram, c, ctmp, hl, hh, hhtmp, entropy;
calcHist( &image, 1, 0, Mat(), hist, 1, &histSize, ranges, true, false );
for( int h = 1; h < histSize; h++){
histogram.push_back(hist.at<float>(h,0));
cout<<histogram.rows<<endl;
cout<<histogram.row(h-1)<<endl;
cout<<hist.row(h)<<endl;
}
histogram=histogram/(image.rows*image.cols-hist.at<float>(0,0));
//cumulative distribution function
float cl,ch;
ctmp.push_back(histogram.row(0));
c.push_back(histogram.row(0));
cout<<c.row(0)<<endl;
for(l=1;l<255;l++){
c.push_back(ctmp.at<float>(0)+histogram.at<float>(l));
ctmp.push_back(c.row(l));
cout<<c.at<float>(l)<<endl;
//c.row(l)=c.row(l-1)+histogram.row(l);
}
Mat hltmp= Mat::zeros(1,256,CV_8U);
// THE PROBLEM IS IN THIS TWO FOR CYCLES
for(t=0;t<255;t++){
//low range entropy
cl=c.at<float>(t);
if(cl>0){
for(i=0;i<=t;i++){
if(histogram.at<float>(t)>0){
printf("here\n");
hl.push_back(hltmp.at<float>(0)-(histogram.at<float> (i)/cl)*log(histogram.at<float>(i)/cl));
printf("here\n");
cout<<hl.at<float>(i);
printf("here\n");
hltmp.push_back(hl.row(t));
printf("here\n");
}
}
}
printf("here\n");
//high range entropy
ch=1.0-cl;
if(ch>0){
for(i=t+1;i<255;i++){
if(histogram.at<float>(i)>0){
hh.push_back(hh.at<float>(t)-(histogram.at<float> (i)/ch)*log(histogram.at<float>(i)/ch));
}
}
}
}
//choose the best threshold
float h_max=hl.at<float>(0,0)+hh.at<float>(0,0);
float threshold=0;
entropy.at<float>(0,0)=h_max;
for(t=1;t<255;t++){
entropy.at<float>(t)=hl.at<float>(t)+hh.at<float>(t);
if(entropy.at<float>(t)>h_max){
h_max=entropy.at<float>(t);
threshold=t-1;
}
cout<<threshold<<endl;
}
//display
Mat I1= Mat::zeros(image.rows,image.cols,CV_8UC1);
for(int y=0; y<image.rows;y++){
for(int x=0; x<image.cols;x++){
if((int) image.at<uchar>(y,x)<threshold){
I1.at<uchar>(y,x)=0;
}
else{
I1.at<uchar>(y,x)=255;
}
}
}
imshow("image",I1);
waitKey(0);*/
return 0;
}
Your problem is that you're reading float elements from a CV_8U (aka uchar) Mat.
Mat hltmp = Mat::zeros(1, 256, CV_8U);
...
hltmp.at<float>(0)
You should learn how to use a debugger, and you'll find out these problems very soon.
Since you over-complicated things in your implementation, made some errors, and the code is cluttered from debug prints, I propose the one below instead of punctually correct your (not many, but mainly conceptual) errors. You can see that, if written properly, there is almost a 1:1 conversion from Matlab to OpenCV.
#include <opencv2/opencv.hpp>
#include <iostream>
using namespace std;
using namespace cv;
uchar maxentropie(const Mat1b& src, Mat1b& dst)
{
// Histogram
Mat1d hist(1, 256, 0.0);
for (int r=0; r<src.rows; ++r)
for (int c=0; c<src.cols; ++c)
hist(src(r,c))++;
// Normalize
hist /= double(src.rows * src.cols);
// Cumulative histogram
Mat1d cumhist(1, 256, 0.0);
float sum = 0;
for (int i = 0; i < 256; ++i)
{
sum += hist(i);
cumhist(i) = sum;
}
Mat1d hl(1, 256, 0.0);
Mat1d hh(1, 256, 0.0);
for (int t = 0; t < 256; ++t)
{
// low range entropy
double cl = cumhist(t);
if (cl > 0)
{
for (int i = 0; i <= t; ++i)
{
if (hist(i) > 0)
{
hl(t) = hl(t) - (hist(i) / cl) * log(hist(i) / cl);
}
}
}
// high range entropy
double ch = 1.0 - cl; // constraint cl + ch = 1
if (ch > 0)
{
for (int i = t+1; i < 256; ++i)
{
if (hist(i) > 0)
{
hh(t) = hh(t) - (hist(i) / ch) * log(hist(i) / ch);
}
}
}
}
// choose best threshold
Mat1d entropie(1, 256, 0.0);
double h_max = hl(0) + hh(0);
uchar threshold = 0;
entropie(0) = h_max;
for (int t = 1; t < 256; ++t)
{
entropie(t) = hl(t) + hh(t);
if (entropie(t) > h_max)
{
h_max = entropie(t);
threshold = uchar(t);
}
}
// Create output image
dst = src > threshold;
return threshold;
}
int main()
{
Mat1b img = imread("path_to_image", IMREAD_GRAYSCALE);
Mat1b res;
uchar th = maxentropie(img, res);
imshow("Original", img);
imshow("Result", res);
waitKey();
return 0;
}