#include "stdafx.h"
#include "threshold.h"
using namespace std;
Log_value::Log_value()
{
is_get=false;
value=0;
}
Log_value threshold(Photoinfo* info, DS_cube * cube, int smoothingoption, int lower, int upper, int threshold_method)
//smoothingoption=0; 1 and 255
{
register unsigned short int i;
unsigned short int iinc;
Log_value result;
result.is_get = false;
unsigned long int* hist=(unsigned long int*)malloc(sizeof(unsigned long int)*256);
if(!hist)
{
write_to_log("Threshold: out of memory for hist");
return result;
}
//initialise histogram
for(i=0;i < 256;i++)
hist[i]=0;
//determine how many slices to analyse-> faster process
iinc = (unsigned short)(1 + (info->sizeofx*info->sizeofy*info->sizeofz) / (100*200*200)); //200^3pixels, arbitrary value
iinc = 1;
{
char log_txt[512];
sprintf(log_txt, "Threshold: #slice step =%d", iinc);
write_to_log(log_txt);
}
if (!cube)
{//allocate memory for a slice if file input type
DS_slice * slice = new DS_slice(info->sizeofy, info->sizeofx);
for(i = 0; i < info->sizeofz; i = i + iinc)
{
//read from file into memory
slice->read_from_file(info, i);
//fill histogram
for (register unsigned short aa = 0; aa < slice->wY; aa++)
for (register unsigned short bb = 0; bb < slice->wX; bb++)
++hist[slice->data[aa][bb]];
}
//free slice
delete slice;
}
else
{//RAW filters performed
for(i = 0; i < info->sizeofz; i = i+iinc)
for (register unsigned short aa = 0; aa < cube->wY; aa++)
for (register unsigned short bb = 0; bb < cube->wX; bb++)
++hist[cube->data[i][aa][bb]];
}
if(smoothingoption==3)
{
unsigned long int* hist3=(unsigned long int*)malloc(sizeof(unsigned long int)*256);
if(!hist3)
{
write_to_log("Threshold: out of memory for hist3");
return result;
}
hist3[0]=hist[0];
for(i=1;i<=254;i++)
hist3[i]=(hist[i-1]+hist[i]+hist[i+1])/3;
hist3[255]=hist[255];
for(i=0;i<=255;i++)
hist[i]=hist3[i];
if (hist3)
free(hist3);
}
if(smoothingoption==5)
{
unsigned long int* hist5=(unsigned long int*)malloc(sizeof(unsigned long int)*256);
if(!hist5)
{
write_to_log("Threshold: out of memory for hist5");
return result;
}
hist5[0]=hist[0];
hist5[1]=(hist[0]+hist[1]+hist[2])/3;
for(i=2;i<=253;i++)
hist5[i]=(hist[i-2]+hist[i-1]+hist[i]+hist[i+1]+hist[i+2])/5;
hist5[255]=hist[255];
hist5[254]=(hist[255]+hist[254]+hist[253])/3;
for(i=0;i<=255;i++)
hist[i]=hist5[i];
if (hist5)
free(hist5);
}
//
//get_threshold using minimum error
//
unsigned short int thresh_min_error=0;
if ( (threshold_method == 1) || (threshold_method == 3) )
{
List<Val_instance>* thresholdlist=new List<Val_instance>;
short iterationsFlag=0,imaginaryFlag=0;
register unsigned short int iterations=0;
short T, Tl = (short)lower, Tu = (short)upper;
float P1=0,P2=0;
double u1=0,u2=0,u11=0,u22=0;
double s1=0,s2=0,s11=0,s22=0;
double a=0,b=0,c1=0,d=0,D=0;
double t,t0,r1=0,r2=0,diff;
while(Tl<=Tu)
{
T=Tl;
t0=(double)T;
do
{
P1=P2=1;
u11=u22=s11=s22=0;
iterations++;
if(iterations>10000)
{
iterationsFlag=1;
break;
}
for(i=0;i<=T;i++)
P1=P1+hist[i];
for(i=(T+1);i<=255;i++)
P2=P2+hist[i];
for(i=0;i<=T;i++)
u11 += i*hist[i];
for(i=(T+1); i<=255; i++)
u22 += i*hist[i];
u1=(u11/P1);
u2=(u22/P2);
for(i=0;i<=T;i++)
s11=s11+(i-u1)*(i-u1)*hist[i];
for(i=(T+1);i<=255;i++)
s22=s22+(i-u2)*(i-u2)*hist[i];
s1=sqrt(s11/P1);
s2=sqrt(s22/P2);
if ((s1 ==0) || (s2 ==0))
break;
a=(1/(s1*s1)-1/(s2*s2));
b=-2*(u1/(s1*s1)-u2/(s2*s2));
c1=(u1*u1)/(s1*s1)-(u2*u2)/(s2*s2);
c1=c1+2*(log(s1)-log(s2))-2*(log(P1)-log(P2));
d=(b*b-4*a*c1);
if(d < 0)
{
imaginaryFlag = 1;
break;
}
D=sqrt(d);
r1=(-b+D)/(2*a);
r2=(-b-D)/(2*a);
if(r1>u1 && r1<u2)
t=r1;
else
t=r2;
if(t0>t)
diff=t0-t;
else
diff=t-t0;
t0 = t;
T = (short)t;
//if ((t-(double)T) > 0.5) //round to closest int
// T ++; //no, because [>] threshold in binarise.cpp
} while (diff >= 0.01);
if(iterationsFlag==1)
{
char log_txt[512];
sprintf(log_txt, "Threshold: iterations > 10000 for Tl=%d", Tl);
write_to_log(log_txt);
iterationsFlag=0;
}
else if(imaginaryFlag==1)
{
char log_txt[512];
sprintf(log_txt, "Threshold: b2-4ac<0 for Tl=%d", Tl);
write_to_log(log_txt);
imaginaryFlag=0;
}
else //fill list
{
if ((T>=0)&&(T<=255))
{
bool is_add=false;
for(i=1; i<=thresholdlist->getNumber(); i++)
if(thresholdlist->getWhich(i)->value==T)
{
thresholdlist->getWhich(i)->instance++;
is_add=true;
break;
}
if(!is_add)
{
Val_instance* newthreshold=new Val_instance;
newthreshold->value=T;
newthreshold->instance=1;
thresholdlist->Add(*newthreshold);
}
}
}
Tl++;
iterations=0;
}
if (thresholdlist->getNumber() > 0)
{
result.is_get=true;
int mostcommon=0;
for(i=1; i<=thresholdlist->getNumber() ;i++) //ignore value >upper and <lower
if ((thresholdlist->getWhich(i)->instance > mostcommon)&&(thresholdlist->getWhich(i)->value <=upper)&&(thresholdlist->getWhich(i)->value >=lower))
{
mostcommon = thresholdlist->getWhich(i)->instance;
thresh_min_error = (unsigned short)thresholdlist->getWhich(i)->value;
}
}
delete thresholdlist;
{
char log_txt[512];
sprintf(log_txt, "Threshold: based on minimum error thresholding: %d", thresh_min_error);
write_to_log(log_txt);
}
}
//get the threshold based on the composite average (ImageJ)
//The automatic thresholding function used by Image/Adjust/Threshold
//divides the image into objects and background. It does this by taking a test threshold
//and computing the average of the pixels at or below the threshold and pixels above.
//It then computes the average of those two, increments the threshold, and repeats the process.
//Incrementing stops when the threshold is larger than the composite average.
//That is, threshold = (average background + average objects)/2
unsigned short int thresh_comp_ave = 0;
if ( (threshold_method == 2) || (threshold_method == 3) )
{
for (thresh_comp_ave=1; thresh_comp_ave<=255; thresh_comp_ave++)
{
unsigned long long av1, av2, s1, s2;
av1=0;
s1=0;
for(i=0; i<thresh_comp_ave; i++)
{
av1 += i*hist[i];
s1 += hist[i];
}
if (!s1)
continue;
av2=0;
s2=0;
for(i=thresh_comp_ave; i<256; i++)
{
av2 += i*hist[i];
s2 += hist[i];
}
if (!s2)
continue;
//double comp_av = round( (av1/s1 + av2/s2) * .5);
//unsigned short int T=(int)comp_av;
//if ((comp_av-(double)T) > 0.5) //round to closest int
// T ++;
//if (thresh_comp_ave > T)
// break;
unsigned short comp_av = (unsigned short)( (double)((double)av1/(double)s1 + (double)av2/(double)s2) * .5);
if (thresh_comp_ave > comp_av)
break;
}
{
char log_txt[512];
sprintf(log_txt, "Threshold: based on the composite average: %d", thresh_comp_ave);
write_to_log(log_txt);
}
}
//free histogram
if (hist)
free(hist);
if (threshold_method == 1) //minimum error
{
result.is_get=true;
result.value = thresh_min_error;
}
else if (threshold_method == 2) //composite average
{
result.is_get=true;
result.value = thresh_comp_ave;
}
else if ( threshold_method == 3 ) //average
{
double r = (double)((double)(thresh_min_error + thresh_comp_ave)/(double)2);
int R = int(r);
//if ((r-(double)R) > 0.5) //round to closest int?
// R ++; //no, because [>] threshold in binarise.cpp
result.is_get=true;
result.value = R;
}
{//write final value out to the log
char log_txt[512];
sprintf(log_txt, "Threshold: final = %d", result.value);
write_to_log(log_txt);
}
return result;
}