text: ﻿#include "LBP.h"
  
  LBP::LBP(LBPParam* param) : mode(param->mode)
  {
  	switch (mode)
  	{
  	case Sqr_R1P8: r = 1, p = 8; break;
  	case Sqr_R2P8: r = 2, p = 8; break;
  	case Sqr_R2P16: r = 2, p = 16; break;
  	case Sqr_R3P8: r = 3, p = 8; break;
  	case Sqr_R3P24: r = 3, p = 24; break;
  	case Sqr_R4P8: r = 4, p = 8; break;
  	case Sqr_R5P8: r = 5, p = 8; break;
  	case Sqr_R6P8: r = 6, p = 8; break;
  	case Sqr_R7P8: r = 7, p = 8; break;
  	case Sqr_R8P8: r = 8, p = 8; break;
  	case Sqr_R9P8: r = 9, p = 8; break;
  	case Cir_R1P8: r = 1, p = 8; break;
  	case Cir_R2P8: r = 2, p = 8; break;
  	case Cir_R2P16: r = 2, p = 16; break;
  	case Cir_R3P8: r = 3, p = 8; break;
  	case Cir_R3P24: r = 3, p = 24; break;
  	case Cir_R4P8: r = 4, p = 8; break;
  	case Cir_R5P8: r = 5, p = 8; break;
  	case Cir_R6P8: r = 6, p = 8; break;
  	case Cir_R7P8: r = 7, p = 8; break;
  	case Cir_R8P8: r = 8, p = 8; break;
  	case Cir_R9P8: r = 9, p = 8; break;
  
  	case Cir_R4P24: r = 4, p = 24; break;
  	case Cir_R5P24: r = 5, p = 24; break;
  	case Cir_R6P24: r = 6, p = 24; break;
  	case Cir_R7P24: r = 7, p = 24; break;
  	case Cir_R8P24: r = 8, p = 24; break;
  	case Cir_R9P24: r = 9, p = 24; break;
  
  	default:break;
  	}
  	PN = 1 << p;
  
  	if (mode < Cir_R1P8)
  	{
  		relIntCoord = new Point[p]{};
  
  		if (p == 8)
  		{
  			int index[8][2] = { {r,0}, {r,1} , {0, r}, {-r,r}, {-r, 0}, {-r, -r}, {0, -r}, {r, -r} };
  			for (size_t i = 0; i < 8; i++)
  			{
  				relIntCoord[i] = Point(index[i][0], index[i][1]);
  			}
  		}
  		else if (p == 16)
  		{
  
  		}
  		else if (p == 24)
  		{
  
  		}
  		//int16_t index[8][2] = { {1,0}, {-1, -1}, {0, 1}, {-1,1}, {0, -1}, {1,1} , {-1, 0}, {1, -1} };
  		//int16_t index[8][2] = { {0, -1}, {1,1} , {1,0}, {-1, -1}, {0, 1}, {-1,1} , {-1, 0}, {1, -1} };
  		
  	}
  	else if (mode >= Cir_R1P8)
  	{
  		relFloatCoord = new Point2d[p]{};
  		relIntCoord = new Point[p]{};
  		for (size_t i = 0; i < p; ++i) 
  		{
  			relFloatCoord[i] = Point2d(r*cos(i * 2 * 3.14159 / p), r * sin(i * 2 * 3.14159 / p));
  		}
  		relIntCoord[0] = Point(r, 0);
  		relIntCoord[p/4] = Point(0, r);
  		relIntCoord[p/2] = Point(-r, 0);
  		relIntCoord[p*3/4] = Point(0, -r);
  	}
  }
  
  void LBP::run(const Mat & src, Mat & dst) const
  {
  	dst = Mat(src.rows - 2 * r, src.cols - 2 * r, CV_32SC1);
  	if (mode < Cir_R1P8)
  	{
  		for (int i = r; i < src.rows - r; ++i)
  		{
  			for (int j = r; j < src.cols - r; ++j)
  			{
  				uint32_t lbpCode = 0;
  				int16_t tmp = 0;
  
  				uint8_t c = src.at<uint8_t>(i, j);
  				for (size_t t = 0; t < p; ++t)
  				{
  					tmp = src.at<uint8_t>(i + relIntCoord[t].y, j + relIntCoord[t].x) - c;
  					lbpCode |= (tmp > 0 ? 1 : 0) << t;
  				}
  				dst.at<int>(i-r, j-r) = lbpCode;
  			}
  		}
  	}
  	else
  	{
  		for (int i = r; i < src.rows - r; ++i)
  		{
  			for (int j = r; j < src.cols - r; ++j)
  			{
  				uint32_t lbpCode = 0;
  				int16_t tmp = 0;
  
  				uint8_t c = src.at<uint8_t>(i, j);
  
  				for (size_t t = 0; t < p; ++t)
  				{
  					if (t == 0 || t == p/4 || t == p /2 || t == p*3/4 )
  					{
  						tmp = src.at<uint8_t>(i + relIntCoord[t].y, j + relIntCoord[t].x) - c;
  					}
  					else
  					{
  						tmp = biLinearInterpolation(src, Point2d(j + relFloatCoord[t].x, i + relFloatCoord[t].y)) - c;
  					}
  					
  					lbpCode |= (tmp > 0 ? 1 : 0) << t;
  				}
  				dst.at<int>(i-r, j-r) = lbpCode;
  			}
  		}
  	}
  }
  
  void LBP::myImgConvert(const Mat& src, Mat& dst)
  {
  	if (p == 8)
  	{
  		dst = Mat(src.rows , src.cols, CV_8UC1);
  		for (int i = 0; i < src.rows; ++i)
  		{
  			int32_t const* pSrc = src.ptr<int32_t>(i);
  			uint8_t * pDst = dst.ptr<uint8_t>(i);
  			for (int j = 0; j < src.cols; ++j)
  			{
  				pDst[j] = pSrc[j];
  			}
  		}
  	}
  	else if (p == 16)
  	{
  		dst = Mat(src.rows, src.cols, CV_16UC1);
  		for (int i = 0; i < src.rows; ++i)
  		{
  			int32_t const* pSrc = src.ptr<int32_t>(i);
  			uint16_t * pDst = dst.ptr<uint16_t>(i);
  			for (int j = 0; j < src.cols; ++j)
  			{
  				pDst[j] = pSrc[j];
  			}
  		}
  	}
  	else
  	{
  		cout << "waiting for realizing!\n";
  	}
  }
  
  MRLBP::MRLBP(MRLBPParam* param_):LBP(param_)
  {
  	lbm = new uint32_t[p]{ 0 };
  	lbw = new uint32_t[p]{ 0 };
  	lbs = new bool[p]{ 0 };
  
  	param = param_;
  
  	if (param->ratio < 1)
  	{
  		if (param->lbpCodeHistgram == nullptr)
  		{
  			param->lbpCodeHistgram = new uint32_t[1 << param->enob]{ 0 };
  		}
  	}
  	else
  	{
  		PN = 1 << param->enob;
  		param->PN = this->PN;
  		if (param->lookupTable == nullptr)
  		{
  			param->lookupTable = new uint32_t[PN];
  			for (size_t i = 0; i < PN; i++)
  			{
  				param->lookupTable[i] = i;
  			}
  		}
  	}
  }
  
  void MRLBP::run(const Mat & src, Mat & dst) const
  {
  	dst = Mat(src.rows - 2 * r, src.cols - 2 * r, CV_32SC1);
  
  	if (mode < Cir_R1P8)
  	{
  		for (int i = r; i < src.rows - r; ++i)
  		{
  			for (int j = r; j < src.cols - r; ++j)
  			{
  				uint8_t c = src.at<uint8_t>(i, j);
  				int16_t tmp = 0;
  
  				for (uint8_t t = 0; t < p; t++)
  				{
  					tmp = src.at<uint8_t>(i + relIntCoord[t].y, j + relIntCoord[t].x) - c;
  					lbm[t] = abs(tmp);
  					lbs[t] = (tmp > 0) ? 1 : 0;
  				}
  				uint32_t lbpCode = 0;
  				getBinWeight(lbm, lbw, p);
  				for (uint32_t t = 0; t < p; ++t)
  				{
  					lbpCode += (lbs[t] << lbw[t]);
  				}
  				lbpCode = lbpCode >> (p - param->enob);
  				dst.at<int>(i - r, j - r) = param->lookupTable[lbpCode];
  			}
  		}
  	}
  	else
  	{
  		for (int i = r; i < src.rows - r; ++i)
  		{
  			for (int j = r; j < src.cols - r; ++j)
  			{
  				uint8_t c = src.at<uint8_t>(i, j);
  				int16_t tmp = 0;
  
  				for (size_t t = 0; t < p; ++t)
  				{
  					if (t == 0 || t == p / 4 || t == p / 2 || t == p * 3 / 4)
  					{
  						tmp = src.at<uint8_t>(i + relIntCoord[t].y, j + relIntCoord[t].x) - c;
  					}
  					else
  					{
  						tmp = biLinearInterpolation(src, Point2d(j + relFloatCoord[t].x, i + relFloatCoord[t].y)) - c;
  					}
  					lbm[t] = abs(tmp);
  					lbs[t] = (tmp > 0) ? 1 : 0;
  				}
  				uint32_t lbpCode = 0;
  				getBinWeight(lbm, lbw, p);
  				for (uint32_t t = 0; t < p; ++t)
  				{
  					lbpCode += (lbs[t] << lbw[t]);
  				}
  				lbpCode = lbpCode >> (p - param->enob);
  				dst.at<int>(i - r, j - r) = param->lookupTable[lbpCode];
  			}
  		}
  	}
  }
  
  void MRLBP::preTrainStep0(const cv::Mat & src)
  {
  	int16_t tmp = 0;
  	if (mode < Cir_R1P8)
  	{
  		for (int i = r; i < src.rows - r; ++i)
  		{
  			for (int j = r; j < src.cols - r; ++j)
  			{
  				uint8_t c = src.at<uint8_t>(i, j);
  
  				for (uint8_t t = 0; t < p; t++)
  				{
  					tmp = src.at<uint8_t>(i + relIntCoord[t].y, j + relIntCoord[t].x) - c;
  					lbm[t] = abs(tmp);
  					lbs[t] = (tmp > 0) ? 1 : 0;
  				}
  				uint32_t lbpCode = 0;
  				getBinWeight(lbm, lbw, p);
  				for (uint32_t t = 0; t < p; ++t)
  				{
  					lbpCode += (lbs[t] << lbw[t]);
  				}
  				lbpCode = lbpCode >> (p - param->enob);
  				threadLock.lock();
  				param->lbpCodeHistgram[lbpCode]++;
  				threadLock.unlock();
  			}
  		}
  	}
  	else
  	{
  		for (int i = r; i < src.rows - r; ++i)
  		{
  			for (int j = r; j < src.cols - r; ++j)
  			{
  				uint8_t c = src.at<uint8_t>(i, j);
  
  				for (uint8_t t = 0; t < p; t++)
  				{
  					if (t == 0 || t == p / 4 || t == p / 2 || t == p * 3 / 4)
  					{
  						tmp = src.at<uint8_t>(i + relIntCoord[t].y, j + relIntCoord[t].x) - c;
  					}
  					else
  					{
  						tmp = biLinearInterpolation(src, Point2d(j + relFloatCoord[t].x, i + relFloatCoord[t].y)) - c;
  					}
  					lbm[t] = abs(tmp);
  					lbs[t] = (tmp > 0) ? 1 : 0;
  				}
  				uint32_t lbpCode = 0;
  				getBinWeight(lbm, lbw, p);
  				for (uint32_t t = 0; t < p; ++t)
  				{
  					lbpCode += (lbs[t] << lbw[t]);
  				}
  				lbpCode = lbpCode >> (p - param->enob);
  				threadLock.lock();
  				param->lbpCodeHistgram[lbpCode]++;
  				threadLock.unlock();
  			}
  		}
  	}
  }
  
  
  void MRLBP::preTrainStep1()
  {
  	threadLock.lock();
  	if (param->lookupTable == nullptr)
  	{
  		uint32_t N = 1 << (param->enob);
  		double* histgramNorm = new double[N] {0};
  		uint64_t histgramAdd = 0;
  		for (uint32_t i = 0; i < N; i++)
  		{
  			histgramAdd += param->lbpCodeHistgram[i];
  		}
  		for (uint32_t i = 0; i < N; i++)
  		{
  			histgramNorm[i] = param->lbpCodeHistgram[i] / static_cast<double>(histgramAdd);
  		}
  
  		std::vector<int32_t> effPattern; //effective pattern
  		effPattern.reserve(N);
  
  		double ratioAdd = 0;
  		double max;
  
  		for (uint32_t i = 0; i < N; i++)//外层循环累加概率 判断是否大于给定ratio，然后跳出循环
  		{
  			effPattern.push_back(0);
  			max = 0;
  			for (uint32_t j = 0; j < N; j++)//内层循环，依次找到最大概率，第二大的概率，，，，不能直接用sort排序，因为排序之后找不到对应的索引值
  			{
  				if (max < histgramNorm[j])
  				{
  					effPattern[i] = j;
  					max = histgramNorm[j];
  				}
  			}
  			ratioAdd += histgramNorm[effPattern[i]];
  			histgramNorm[effPattern[i]] = 0;//清零，下次循环还是找最大值
  			if (ratioAdd >= param->ratio)
  				break;
  		}
  		delete[] histgramNorm;
  
  		PN = effPattern.size();
  		std::sort(effPattern.begin(), effPattern.end(), std::less<int>());
  
  		param->PN = this->PN;
  		param->lookupTable = new uint32_t[N]{};
  		
  		//uint32_t margin = param->margin;
  		for (uint32_t i = 0; i < N; i++)//查表法，此处填表
  		{
  			uint32_t /*diff1, diff2, */diff, minDiff = 2 * N;
  			param->lookupTable[i] = PN-1;
  			for (uint32_t j = 0; j < effPattern.size(); j++)
  			{
  				diff = abs(effPattern[j] - i);
  				if (diff < minDiff)
  				{
  					minDiff = diff;
  					param->lookupTable[i] = j;
  					continue;
  				}
  			}
  		}
  	}
  	threadLock.unlock();
  }
  
  
  void GNP::run(const Mat & src, Mat & dst) const
  {
  	uint32_t maxGrayLevel = 256;
  	uint32_t*  cHist = new uint32_t[maxGrayLevel]{ 0 };
  	uint8_t* cLookupTable = new uint8_t[maxGrayLevel]{ 0 };
  	for (size_t i = 0; i < src.rows; i++)
  	{
  		const uint8_t* pData = src.ptr<uint8_t>(i);
  		for (size_t j = 0; j < src.cols; j++)
  		{
  			cHist[pData[j]]++;
  		}
  	}
  
  	uint32_t pointsNum = src.rows*src.cols;
  
  	uint8_t* median = new uint8_t[this->mode]{ 0 };
  	uint32_t tmp1 = 0;
  	for (uint16_t i = 0, j = 1; i < maxGrayLevel; i++)
  	{
  		tmp1 += cHist[i];
  		if (tmp1 > pointsNum * j / double(PN))
  		{
  			median[j - 1] = i;
  			if (j == PN)
  				break;
  			j++;
  		}
  	}
  	for (uint16_t i = 0, j = 0; i < maxGrayLevel; i++)
  	{
  		if (j < this->mode)
  		{
  			if (i <= median[j])
  			{
  				cLookupTable[i] = j;
  			}
  			else
  			{
  				j++;
  				cLookupTable[i] = j;
  			}
  		}
  		else
  		{
  			cLookupTable[i] = j;
  		}
  	}
  
  	dst = Mat(src.rows, src.cols, CV_8UC1);
  	for (size_t i = 0; i < src.rows; i++)
  	{
  		const uint8_t* pSrcData = src.ptr<uint8_t>(i);
  		uint8_t* pDstData = dst.ptr<uint8_t>(i);
  		for (size_t j = 0; j < src.cols; j++)
  		{
  			pDstData[j] = cLookupTable[pSrcData[j]];
  		}
  	}
  	delete[] median;
  	delete[] cHist;
  	delete[] cLookupTable;
  }
  