Visione Artificiale - Algoritmo di Otsu

#include <opencv2\opencv.hpp>
#include <map>
#include <fstream>
#include <iterator>
#include <iostream>


using namespace cv;
using namespace std;

int main(){
Mat1b imm = imread("../images/Immagini/Tesi00.bmp", 0);
ofstream his, hisn;
his.open("Istogramma.txt");
hisn.open("Istogramma normalizzato.txt");

//Istogramma
map<uchar, int> histo;
map<uchar, int>::iterator im;

for(int i=0;i<=255;++i)
histo.insert(pair<uchar, int>(i,0));

for(int r=0;r<imm.rows;++r){
for(int c=0;c<imm.cols;++c){
histo[imm(r,c)]++;
}
}

for(im=histo.begin();im!=histo.end();++im)
his << (int)im->first << "\t" << im->second << "\n";
his.close();
 
//Istogramma normalizzato
map<uchar, float> histo_n;
map<uchar, float>::iterator imn;

for(int i=0;i<=255;++i)
histo_n.insert(pair<uchar, float>(i,0.0));

for(int i=0;i<=255;++i)
histo_n[i] = (float)histo[i]/(imm.rows*imm.cols);

for(imn=histo_n.begin();imn!=histo_n.end();++imn)
hisn << (int)imn->first << "\t" << imn->second << "\n";
hisn.close();

//Otsu
int soglia = 0;
float min;
int flag = 0;

for(int t=0;t<=255;++t){
float q1=0.0, q2=0.0;
float m1=0.0, m2=0.0;
float s1=0.0, s2=0.0;
float costo=0.0;//per calcolo del minore

//Probabilità a priori
for(int i=0;i<=t;++i)
q1 = q1 + histo_n[i];
for(int j=t+1;j<=255;++j)
q2 = q2 + histo_n[j];

if(q1!=0 && q2!=0){
//Media
for(int i=0;i<=t;++i)
m1 = m1 + (i*histo_n[i])/q1;
for(int j=t+1;j<=255;++j)
m2 = m2 + (j*histo_n[j]/q2);

//Varianza
for(int i=0;i<=t;++i)
s1 = s1 + ((i-m1)*(i-m1))*histo_n[i]/q1;
for(int j=t+1;j<=255;++j)
s2 = s2 + ((j-m2)*(j-m2))*histo_n[j]/q2;

//Costo
costo = q1*s1 + q2*s2;
cout << costo << "\t" << t << "\n";
//inizializza minimo al primo valore calcolato
if(flag==0){
min = costo;
flag = 1;
}
//aggiorno la soglia se trovato valore migliore
if(costo < min){
soglia = t;
min = costo;
}
}
}

//Soglia con valore calcolato
Mat1b out(imm.size());

cout << "Valore soglia selezionato: " << soglia << "\n";
for(int r=0;r<imm.rows;++r)
for(int c=0;c<imm.cols;++c){
if(imm(r,c)>soglia)
out(r,c) = 255;
else
out(r,c) = 0;
}
imshow("Soglia", out);
waitKey();
}