text: /**
   * LBP.cpp
   * Implements the local binary pattern (LBP) texture descriptors
   *
   *  Created on: Jan 25, 2013
   *      Author: Navid Nourani-Vatani
   *      Email: Navid.Nourani-Vatani@sydney.edu.au
   *
   *  The methods implemented here are inspired by the Matlab code available
   *  from web site of the University of Oulu:
   *  	http://www.cse.oulu.fi/CMV/Downloads/LBPMatlab
   *  You should cite the appropriate publications when using this code.
   *
   *  This program is free software: you can redistribute it and/or modify
   *  it under the terms of the GNU 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 General Public License for more details.
   *
   *  You should have received a copy of the GNU General Public License
   *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
   *
   */
  
  #include "lbp/LBP.hpp"
  
  using namespace lbp;
  
  /**
   * Constructors
   */
  LBP::LBP(void)
  	: type(LBP_MAPPING_NONE), samples(0), num(0), fftIn(NULL), fftOut(NULL), fftPlan(NULL),
  	fftN(0), fftHermN(0) {
  }
  LBP::LBP(unsigned int _samples, MappingType _type)
  	: type(_type), samples(_samples), num(0), fftIn(NULL), fftOut(NULL),
  	fftPlan(NULL), fftN(0), fftHermN(0) {
  	generateMapping();
  
  	if (type == LBP_MAPPING_HF) {
  		initHF();
  	}
  }
  LBP::~LBP() {
  	if (fftIn != NULL) {
  		fftw_destroy_plan(fftPlan);
  		fftw_free(fftIn);
  		delete[] fftOut;
  	}
  }
  
  /** ******************************************************************
   *
   * Mapping part
   *
   */
  LBP & LBP::generateMapping() {
  	return generateMapping(this->samples, this->type);
  }
  
  LBP & LBP::generateMapping(unsigned int samples, MappingType type) {
  	this->orbits.clear();
  	this->table.clear();
  	this->num = 0;
  	this->type = type;
  	this->samples = samples;
  
  	int newMax = 0; //number of patterns in the resulting LBP code
  	int index = 0;
  
  	if (type == LBP_MAPPING_NONE) {
  		newMax = (int)pow(2., (int)samples);
  		for (unsigned int i = 0; i < newMax; i++) {
  			table.push_back(i);
  		}
  	}
  	else if (type == LBP_MAPPING_U2) {
  		// Uniform 2
  		newMax = samples * (samples - 1) + 3;
  
  		for (unsigned int i = 0; i < pow(2., (int)(samples)); i++) {
  
  			// Rotate left
  			//			unsigned int bg = ((i & (1 << (samples - 1))) >> (samples - 1)); // bitget(i,samples)
  			//			unsigned int bs = (i << 1) & ((int) pow( 2., (int) samples ) - 1); // bitshift( i, 1, samples )
  			//			unsigned int j = (bs + bg) & ((int) pow( 2., (int) samples ) - 1); // bitset( bs, 1, bg )
  			unsigned int j = rotateLeft(i, samples);
  
  			int numt = NumberOfSetBits(i ^ j); // sum(bitget(bitxor(i,j),1:samples));
  			//number of 1->0 and 0->1 transitions
  			//in binary string
  			//x is equal to the
  			//number of 1-bits in
  			//XOR(x,Rotate left(x))
  
  			if (numt <= 2) {
  				table.push_back(index);
  				index = index + 1;
  			}
  			else {
  				table.push_back(newMax - 1);
  			}
  		}
  	}
  	else if (type == LBP_MAPPING_RI) {
  		long N = (int)pow(2., (int)samples);
  		// Rotation Invariant
  		int * tmpMap = new int[N];
  		memset((void *)tmpMap, -1, N);
  
  		for (unsigned long i = 0; i < N; i++) {
  			tmpMap[i] = -1;
  
  			unsigned long rm = i;
  			unsigned long r = i;
  			for (int j = 1; j <= samples - 1; j++) {
  				r = rotateLeft(r, samples);
  				if (r < rm)
  					rm = r;
  			}
  			if (tmpMap[rm] < 0) {
  				tmpMap[rm] = newMax;
  				newMax = newMax + 1;
  			}
  			table.push_back(tmpMap[rm]);
  		}
  	}
  	else if (type == LBP_MAPPING_RIU2) {
  		// Rotation invariant uniform 2
  		newMax = samples + 2;
  		for (unsigned int i = 0; i <= pow(2., (int)samples) - 1; i++) {
  			unsigned int j = rotateLeft(i, samples); //bitset( bitshift( i, 1, samples ), 1, bitget( i, samples ) ); // rotate left
  			unsigned int numt = NumberOfSetBits(i ^ j); //sum(bitget(bitxor(i,j),1:samples));
  			if (numt <= 2)
  				table.push_back(NumberOfSetBits(i));
  			else
  				table.push_back(samples + 1);
  		}
  	}
  	else if (type == LBP_MAPPING_HF) {
  		// Histogram Fourier
  		newMax = samples * (samples - 1) + 3;
  		table.push_back(newMax - 3);
  
  		for (unsigned int i = 1; i <= pow(2., (int)samples) - 2; i++) {
  			unsigned int j = rotateLeft(i, samples); // bitset(bitshift(i,1,samples),1,bitget(i,samples)); %rotate left
  			unsigned int numt = NumberOfSetBits(i ^ j); // sum(bitget(bitxor(i,j),1:samples)); %number of 1->0 and 0->1 transitions
  
  			if (numt == 2) { // Uniform pattern
  				unsigned int n = NumberOfSetBits(i); // sum(bitget(i,1:samples)); %Number of 1-bits
  
  				unsigned int bc = j ^ ((unsigned int)pow(2., (int)samples) - 1);
  				unsigned int ba = bc & i;
  				unsigned int f = trailingZeroInd(ba) + 1; //find(bitget(bitand(i,bitcmp(j,samples)),1:samples)); // Rotation index of the bit pattern
  				unsigned int r = ((int)floor(n / 2.) + f) % samples;
  				index = (n - 1) * samples + r;
  				table.push_back(index);
  			}
  			else { // Non-uniform
  				table.push_back(newMax - 1);
  			}
  		}
  		table.push_back(newMax - 2);
  
  		vector<int> o;
  		for (int i = 1; i <= samples - 1; i++) {
  			o.clear();
  			for (int j = ((i - 1) * samples); j <= (i * samples - 1); j++) {
  				o.push_back(j);
  			}
  			orbits.push_back(o);
  		}
  		o.clear();
  		o.push_back(newMax - 3);
  		orbits.push_back(o);
  		o[0] = newMax - 2;
  		orbits.push_back(o);
  		o[0] = newMax - 1;
  		orbits.push_back(o);
  	}
  	else {
  		cerr << "Unknown mapping!" << endl;
  		exit(1);
  	}
  
  	this->num = newMax;
  
  	return *this;
  }
  
  /**
   *
   */
  
  bool LBP::saveMapping(string fileName) {
  	ofstream ofs(fileName.c_str(), ios::out);
  	if (!ofs) {
  		cerr << "File \'" << fileName << "\' could not be opened" << endl;
  		return false;
  	}
  
  	ofs << "LBPMapping" << endl;
  	ofs << "version " << 1 << endl;
  	ofs << "type " << MappingTypeStr[type] << endl;
  	ofs << "samples " << samples << endl;
  	ofs << "maxnum " << num << endl;
  	ofs << "table ";
  	for (int i = 0; i < table.size(); i++) {
  		ofs << table[i] << " ";
  	}
  	ofs << endl;
  	if (type == LBP_MAPPING_HF) {
  		ofs << "orbits ";
  		for (int i = 0; i < orbits.size(); i++) {
  			for (int j = 0; j < (orbits[i]).size(); j++) {
  				ofs << orbits[i][j] << " ";
  			}
  			ofs << "-1 ";
  		}
  		ofs << endl;
  	}
  
  
  	return true;
  }
  bool LBP::loadMapping(string fileName) {
  	ifstream ifs(fileName.c_str(), ios::in);
  	if (!ifs) {
  		cerr << "File \'" << fileName << "\' could not be opened" << endl;
  		return false;
  	}
  
  
  
  	string s; int i;
  	// Get file type
  	ifs >> s;
  	if (s.compare("LBPMapping")) {
  		cerr << fileName << " is not a LBPMapping file." << endl;
  		return false;
  	}
  
  	// Get verion
  	ifs >> s >> i;
  
  	// Get mapping type
  	ifs >> s >> s;
  	this->type = strToType(s);
  
  	// Get samples
  	ifs >> s >> this->samples;
  
  	// Get maxnum
  	ifs >> s >> this->num;
  
  	// Get table
  	ifs >> s;
  	this->table.clear();
  	for (int j = 0; j < pow(2., (double)samples); j++) {
  		ifs >> i;
  		table.push_back(i);
  	}
  
  	if (type != LBP_MAPPING_HF) {
  		return true;
  	}
  	// Get orbits for HF
  	this->orbits.clear();
  	ifs >> s;
  	vector<int> o;
  	while (ifs >> i) {
  		if (i < 0) { // -1 are used as separators
  			orbits.push_back(o);
  			o.clear();
  			continue;
  		}
  		o.push_back(i);
  	}
  
  	return true;
  }
  
  /** ******************************************************************
   *
   * Descriptor part
   *
   */
  LBP & LBP::calcLBP(Mat d_img, double radius, bool borderCopy) {
  
  	//	clock_t startTime, endTime, sT, eT;
  	//	vector<double> times;
  	//	double minVal, maxVal;
  	//	namedWindow( "lbp", 0 );
  	//	Mat dummy( 300, 260, CV_8UC1);
  
  	// Make sure the image has Double precision version
  	if (d_img.type() < CV_64F) {
  		d_img.convertTo(d_img, CV_64F);
  	}
  
  	// Make a copy of the image border the same size as the radius
  	if (borderCopy) {
  		Mat tmp(d_img.rows + 2 * radius, d_img.cols + 2 * radius, CV_64F);
  		copyMakeBorder(d_img, tmp, radius, radius, radius, radius, BORDER_WRAP, Scalar(0));
  		d_img = tmp.clone();
  	}
  
  	const unsigned int nSamples = samples;
  	double** spoints = new double*[nSamples];
  	double a = 2 * M_PI / samples;
  	double miny = +INT_MAX;
  	double maxy = -INT_MAX;
  	double minx = +INT_MAX;
  	double maxx = -INT_MAX;
  	for (int i = 0; i < samples; i++) {
  		spoints[i] = new double[2];
  		spoints[i][0] = +radius * cos(double(i * a));
  		spoints[i][1] = -radius * sin(double(i * a));
  
  		minx = (spoints[i][0] < minx ? spoints[i][0] : minx);
  		maxx = (spoints[i][0] > maxx ? spoints[i][0] : maxx);
  		miny = (spoints[i][1] < miny ? spoints[i][1] : miny);
  		maxy = (spoints[i][1] > maxy ? spoints[i][1] : maxy);
  
  	}
  
  	// Determine the dimensions of the input image.
  	int xsize = d_img.cols;
  	int ysize = d_img.rows;
  
  	// Block size, each LBP code is computed within a block of size bsizey*bsizex
  	int bsizex = ceil(max(maxx, 0.)) - floor(min(minx, 0.)) + 1;
  	int bsizey = ceil(max(maxy, 0.)) - floor(min(miny, 0.)) + 1;
  
  	// Minimum allowed size for the input image depends
  	// on the radius of the used LBP operator.
  	if (xsize < bsizex || ysize < bsizey) {
  		cerr << "Too small input image. Should be at least (2*radius+1) x (2*radius+1)" << endl;
  		return *this;
  	}
  
  	// Coordinates of origin (0,0) in the block
  	int origx = 1 - floor(min(minx, 0.)) - 1;
  	int origy = 1 - floor(min(miny, 0.)) - 1;
  
  	// Calculate dx and dy;
  	int dx = xsize - bsizex + 1;
  	int dy = ysize - bsizey + 1;
  
  	// Fill the center pixel matrix C.
  	// d_C is a shallow copie. But that's OK because we are not changing the values
  	//	but only comparing to N
  	Mat d_C(d_img, Rect(origx, origy, dx, dy));
  
  	// Initialize the result matrix with zeros.
  	Mat result(dy, dx, CV_64FC1);
  	result = Mat::zeros(dy, dx, CV_64FC1);
  	Mat D(dy, dx, CV_64FC1);
  	Mat N(dy, dx, CV_64FC1);
  
  	// Compute the LBP code image
  	//	startTime = clock();
  	for (int i = 0; i < samples; i++) {
  		double x = spoints[i][0] + origx;
  		double y = spoints[i][1] + origy;
  		// Calculate floors, ceils and rounds for the x and y.
  		int fy = floor(y);
  		int cy = ceil(y);
  		int ry = round(y);
  		int fx = floor(x);
  		int cx = ceil(x);
  		int rx = round(x);
  
  		// Check if interpolation is needed.
  		if ((fabs(x - rx) < 1e-6) && (fabs(y - ry) < 1e-6)) {
  			// Interpolation is not needed, use original data types
  			// N is a shallow copy. But that's OK because we are not changing the value
  			//	but only comparing to C
  			Mat N(d_img, Rect(rx, ry, dx, dy));
  			compare(N, d_C, D, CMP_GE); // D = (N >= C);
  		}
  		else {
  			// Interpolation needed, use double type images
  			double tx = x - fx;
  			double ty = y - fy;
  
  			// Calculate the interpolation weights.
  			double w1 = (1 - tx) * (1 - ty);
  			double w2 = tx * (1 - ty);
  			double w3 = (1 - tx) * ty;
  			double w4 = tx * ty;
  
  			// Compute interpolated pixel values
  			//			N = w1 * d_img( Rect( fx, fy, dx, dy ) ) + w2 * d_img( Rect( cx, fy, dx, dy ) )
  			//						+ w3 * d_img( Rect( fx, cy, dx, dy ) )
  			//						+ w4 * d_img( Rect( cx, cy, dx, dy ) );
  			// The below operations are about 20% faster than the above
  			addWeighted(d_img(Rect(fx, fy, dx, dy)), w1,
  				d_img(Rect(cx, fy, dx, dy)), w2, 0, N);
  			addWeighted(d_img(Rect(fx, cy, dx, dy)), w3, N, 1, 0, N);
  			addWeighted(d_img(Rect(cx, cy, dx, dy)), w4, N, 1, 0, N);
  
  			compare(N, d_C, D, CMP_GE); // D = (N >= C);
  		}
  		delete[] spoints[i];
  
  		// Update the result matrix.
  		double v = pow(2., i) / 255.; // Divide by 255 because D is 0/255 rather than 0/1
  		D.convertTo(D, CV_64F);
  		result = result + (v * D);
  
  	}
  	delete[] spoints;
  	//	endTime = clock();
  	//	times.push_back( (endTime - startTime) );
  	//	cout << "lbp calc took " << times.back() << " cycles" << endl;
  
  	//	startTime = clock();
  	// Apply mapping if it is defined
  	if (type != LBP_MAPPING_NONE) {
  		MatIterator_<double> it = result.begin<double>(), it_end = result.end<
  			double>();
  		for (; it != it_end; ++it) {
  			*it = table[(*it)];
  		}
  	}
  	result.convertTo(result, CV_8U);
  	//	endTime = clock();
  	//	times.push_back( (endTime - startTime) );
  	//	cout << "mapping took " << times.back() << " cycles" << endl;
  
  	// Store the final result
  	lbpImage = result.clone();
  #if 0
  	for (int j = 0; j < lbpImage.rows; j++) {
  		for (int i = 0; i < lbpImage.cols; i++) {
  			//cout.width(3);
  			cout << int(lbpImage.at<unsigned char>(j, i)) << " ";
  		}
  		cout << endl;
  	}
  #endif
  	return *this;
  }
  
  bool LBP::saveLBPImage(string fileName) {
  	return cv::imwrite(fileName, this->lbpImage);
  }
  /** ******************************************************************
   *
   * Histogram part
   *
   */
  LBP & LBP::calcHist(void) {
  	return calcHist(&lbpImage);
  }
  LBP & LBP::calcHist(Mat mask) {
  	return calcHist(&lbpImage, &mask);
  }
  LBP & LBP::calcHist(Mat * lbpImg, Mat * mask) {
  
  	if (lbpImg == NULL) {
  		lbpImg = &(this->lbpImage);
  	}
  
  	int histSize = num;
  	float range[] = { 0, num };
  	const float* histRange = { range };
  	if (mask == NULL) {
  		cv::calcHist(lbpImg, 1, 0, Mat(), // do not use mask
  			hist, 1, &histSize, &histRange, true, // the histogram is uniform
  			false // do not accumulate
  		);
  	}
  	else {
  		cv::calcHist(lbpImg, 1, 0, *mask, // use mask
  			hist, 1, &histSize, &histRange, true, // the histogram is uniform
  			false // do not accumulate
  		);
  	}
  	return *this;
  }
  
  vector<double> LBP::getHist(bool norm) {
  	vector<double> h(hist.rows);
  	Scalar sum(1);
  
  	// normalization value
  	if (norm || type == LBP_MAPPING_HF) {
  		sum = cv::sum(hist);
  	}
  
  	for (int i = 0; i < hist.rows; i++) {
  		h[i] = hist.at<float>(i) / sum[0];
  	}
  
  	if (type == LBP_MAPPING_HF) {
  		h = constructHF(h);
  	}
  
  	return h;
  }
  
  void LBP::initHF(void) {
  	// All the vectors in the orbit have at least the same size as the first one.
  	//	only the last 3 are off size 1 which are not converted anyway
  	fftN = this->orbits[0].size();
  	// Since the input data are real we take advantage of Hermetian redundancy. This
  	// 	gives us a speed up
  	fftHermN = floor(fftN / 2) + 1;
  
  	// If the size of this fft array is different from previous
  	if (fftN != fftN && fftIn != NULL) {
  		fftw_free(fftIn);
  		delete[] fftOut;
  		fftIn = NULL;
  		fftOut = NULL;
  	}
  	// If the in/out arrays are uninitialized
  	if (fftIn == NULL) {
  		fftIn = (double *)fftw_malloc(sizeof(double) * fftN);
  		fftOut = new complex<double>[fftHermN];
  	}
  	// Setup the fft plan
  	fftPlan = fftw_plan_dft_r2c_1d(
  		fftN, fftIn, reinterpret_cast<fftw_complex *>(fftOut), FFTW_ESTIMATE);
  
  }
  
  vector<double> LBP::constructHF(vector<double> h) {
  	if (this->type != LBP_MAPPING_HF) {
  		cerr << "The mapping type must be " << MappingTypeStr[LBP_MAPPING_HF] << endl;
  		return h;
  	}
  
  	initHF();
  
  	// Size of the output vector
  	//	int FVLEN = (samples - 1) * (floor( samples / 2 ) + 1) + 3;
  	//	hf.reserve( FVLEN );
  	hf.clear();
  
  	for (int j = 0; j < this->orbits.size(); j++) {
  		if (orbits[j].size() > 1) {
  			// transfer in the data
  			for (int i = 0; i < fftN; i++) {
  				fftIn[i] = h[orbits[j][i]];
  			}
  
  			fftw_execute(fftPlan);
  
  			// read out the data
  			for (int i = 0; i < fftHermN; i++) {
  				hf.push_back(abs(fftOut[i]));
  			}
  		}
  		else {
  			hf.push_back(h[orbits[j][0]]);
  		}
  	}
  
  	return hf;
  }
  
  /** ******************************************************************
   *
   * Others part
   *
   */
  //std::string LBP::toString(bool verbose) const {
  //	string s = "LBP = \n";
  //	s += "\t    type: " + MappingTypeStr[type] + "\n";
  //	if (verbose) {
  //		s += "\t   table: [";
  //		for (int i = 0; i < table.size(); i++)
  //			s += SSTR(table[i]) + (i < table.size() - 1 ? ", " : "");
  //		s += "]\n";
  //	}
  //	else {
  //		s += "\t   table: [1x" + SSTR(pow(2., (int)this->samples)) + (string) "]\n";
  //	}
  //	s += "\t samples: " + SSTR(this->samples) + (string) "\n";
  //	s += "\t     num: " + SSTR(this->num) + (string) "\n";
  //	if (this->type == LBP_MAPPING_HF) {
  //		if (verbose) {
  //			s += "\t  orbits: {";
  //			for (int i = 0; i < orbits.size(); i++) {
  //				s += "{";
  //				for (int j = 0; j < orbits[i].size(); j++) {
  //					s += SSTR(orbits[i][j]) + (j < orbits[i].size() - 1 ? ", " : "");
  //				}
  //				s += (string)"}" + (i < orbits.size() - 1 ? ", " : "");
  //			}
  //			s += "}\n";
  //		}
  //		else {
  //			s += "\t  orbits: {" + SSTR(this->orbits.size()) + (string) "x1}\n";
  //		}
  //	}
  //	return s;
  //}