/*
Copyright (c) 2023 Dariusz Borkowski

This program is free software: you can redistribute it and/or modify it under
the terms of the GNU Affero General Public License as published by the Free Software Foundation,
either version 3 of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
See the GNU Affero General Public License for more details.

You should have received a copy of the GNU Affero General Public License along with this program.
If not, see <http://www.gnu.org/licenses/>.
*/

#include "libsegmentation.h"
/*---------------------------Algorithm 4: BSDE segmentation---------------------------------*/
void BSDE_Segmentation(dbImageRGB& Input,float mu,dbImageRGB& Output){
	float Sigma=mu;
	float c=Sigma*0.01f;
	BSDEdenoisingRGB(Input,Sigma,c,Output);
	Sigma=0.5f*Sigma;
	while(Sigma>5){					/* epsilon=5 */
		c=Sigma*0.01f;				/* cn */
		imCopy(Output,Input);
		BSDEdenoisingRGB(Input,Sigma,c,Output);	/* un */
		Sigma=0.5f*Sigma;			/* sigman */
	}
}
/*------------------------Algorithm 5: BSDE binarization where Input=Algorithm 4.-----------*/
void BSDE_Binarization(dbImageRGB& Input,dbImageRGB& Output){
/*------------------------Otsu's method for the red componenet------------------------------*/
	int nx=Input.getWidth();
	int ny=Input.getHeight();
	int Lhist=256;
	int** pixc=new int*[nx];
	for(int i=0;i<nx;i++)pixc[i]=new int[ny];
	double* histog=new double[Lhist];
	double* Hs=new double[Lhist];
	double Pinf,sinf,ssup,Psup,muinf,musup;
	double mm;
	double  th;
	for(int x=0;x<nx;x++)
		for(int y=0;y<ny;y++)pixc[x][y]=std::min(Lhist-1,std::max(0,int(round(Input(x,y).R))));
	for(int i=0;i<Lhist;i++)histog[i]=0.0;
	for(int x=0;x<nx;x++)
		for(int y=0;y<ny;y++)
			for(int k=0;k<Lhist;k++)if(pixc[x][y]==k)histog[k]++;
	for(int k=0;k<Lhist;k++)histog[k]=histog[k]/double(nx*ny);
	Pinf=0.0;sinf=0.0;ssup=0.0;
	for(int i=0;i<Lhist;i++)ssup+=i*histog[i];
	for(int i=0;i<Lhist;i++)Hs[i]=0.0;
	for(int S=0;S<Lhist-1;S++){
		Pinf=Pinf+histog[S];Psup=1-Pinf;
		sinf=sinf+S*histog[S];ssup=ssup-S*histog[S];
		if(Pinf == 0.0)muinf=0.0;	
		else muinf=sinf/Pinf;
		if(Psup == 0.0)musup=0.0;
		else musup=ssup/Psup;
		Hs[S]=Pinf*Psup*(muinf-musup)*(muinf-musup);
	}
	mm=Hs[0];
	th=0;
	for(int i=1;i<Lhist;i++)if(mm<Hs[i]){mm=Hs[i];th=i;}
	for(int x=0;x<nx;x++)
		for(int y=0;y<ny;y++)
			if(pixc[x][y]>th)Output(x,y).R=255.0f;
			else Output(x,y).R=0.0f;
/*------------------------Otsu's method for the green componenet------------------------------*/
	for(int x=0;x<nx;x++)
		for(int y=0;y<ny;y++)pixc[x][y]=std::min(Lhist-1,std::max(0,int(round(Input(x,y).G))));
	for(int i=0;i<Lhist;i++)histog[i]=0.0;
	for(int x=0;x<nx;x++)
		for(int y=0;y<ny;y++)
			for(int k=0;k<Lhist;k++)if(pixc[x][y]==k)histog[k]++;
	for(int k=0;k<Lhist;k++)histog[k]=histog[k]/double(nx*ny);
	Pinf=0.0;sinf=0.0;ssup=0.0;
	for(int i=0;i<Lhist;i++)ssup+=i*histog[i];
	for(int i=0;i<Lhist;i++)Hs[i]=0.0;
	for(int S=0;S<Lhist-1;S++){
		Pinf=Pinf+histog[S];Psup=1-Pinf;
		sinf=sinf+S*histog[S];ssup=ssup-S*histog[S];
		if(Pinf == 0.0)muinf=0.0;	
		else muinf=sinf/Pinf;
		if(Psup == 0.0)musup=0.0;
		else musup=ssup/Psup;
		Hs[S]=Pinf*Psup*(muinf-musup)*(muinf-musup);
	}
	mm=Hs[0];
	th=0;
	for(int i=1;i<Lhist;i++)if(mm<Hs[i]){mm=Hs[i];th=i;}
	for(int x=0;x<nx;x++)
		for(int y=0;y<ny;y++)
			if(pixc[x][y]>th)Output(x,y).G=255.0f;
			else Output(x,y).G=0.0f;
/*------------------------Otsu's method for the blue componenet------------------------------*/
	for(int x=0;x<nx;x++)
		for(int y=0;y<ny;y++)pixc[x][y]=std::min(Lhist-1,std::max(0,int(round(Input(x,y).B))));
	for(int i=0;i<Lhist;i++)histog[i]=0.0;
	for(int x=0;x<nx;x++)
		for(int y=0;y<ny;y++)
			for(int k=0;k<Lhist;k++)if(pixc[x][y]==k)histog[k]++;
	for(int k=0;k<Lhist;k++)histog[k]=histog[k]/double(nx*ny);
	Pinf=0.0;sinf=0.0;ssup=0.0;
	for(int i=0;i<Lhist;i++)ssup+=i*histog[i];
	for(int i=0;i<Lhist;i++)Hs[i]=0.0;
	for(int S=0;S<Lhist-1;S++){
		Pinf=Pinf+histog[S];Psup=1-Pinf;
		sinf=sinf+S*histog[S];ssup=ssup-S*histog[S];
		if(Pinf == 0.0)muinf=0.0;	
		else muinf=sinf/Pinf;
		if(Psup == 0.0)musup=0.0;
		else musup=ssup/Psup;
		Hs[S]=Pinf*Psup*(muinf-musup)*(muinf-musup);
	}
	mm=Hs[0];
	th=0;
	for(int i=1;i<Lhist;i++)if(mm<Hs[i]){mm=Hs[i];th=i;}
	for(int x=0;x<nx;x++)
		for(int y=0;y<ny;y++)
			if(pixc[x][y]>th)Output(x,y).B=255.0f;
			else Output(x,y).B=0.0f;
	delete[] Hs;	
	for(int i=0;i<nx;i++)delete[] pixc[i];
	delete[] pixc;
	delete[] histog;
}
/*------------------------Central region of the segmented image------------------------------*/
void getRegion(int x,int y,dbImageRGB& Input,dbImage& Output){
	struct intint{
		int x;
		int y;
	};
	struct intintint{
		int x;
		int y;
		int z;
	};
	intint ii{.x=x,.y=y};
	int nx=std::min(Input.getWidth(),Output.getWidth());
	int ny=std::max(Input.getHeight(),Output.getHeight());
	std::stack<intint> s;
	intintint Q={.x=int(Input(x,y).R),.y=int(Input(x,y).G),.z=int(Input(x,y).B)};
	for(int i=0;i<nx;i++)for(int j=0;j<ny;j++)Output(i,j)=0.0f;
	Output(ii.x,ii.y)=255.0f;
	s.push(ii);	
	while(!s.empty()){
	ii=s.top();s.pop();
	for(int xx=ii.x-1;xx<=ii.x+1;++xx)
		for(int yy=ii.y-1;yy<=ii.y+1;++yy)
if((!(Output(xx,yy)>0.0f)) && (int(Input(xx,yy).R)==Q.x) && (int(Input(xx,yy).G)==Q.y) && (int(Input(xx,yy).B)==Q.z)){
				Output(xx,yy)=255.0f;
				ii.x=xx;ii.y=yy;
				s.push(ii);
			}
	}
}
