OTS/OTSCPP/OTSImagePro/OTSImageProcess.cpp

#pragma once
#include "stdafx.h"
#include <opencv2/core/core.hpp>  
#include <opencv2/highgui/highgui.hpp> 
#include <opencv2/opencv.hpp>
#include "OTSImageProcess.h"
#include "OTSImageProcessParam.h"
#include <OTSFieldData.h>
#include "OTSMorphology.h"
#include "../OTSLog/COTSUtilityDllFunExport.h"
#include "FieldMgr.h"
#include "BaseFunction.h"
#include "GreyPeak.h"


namespace OTSIMGPROC
{
	using namespace cv;
	using namespace std;

	const int nBlackColor = 255;
	namespace
	{
	
		cv::Mat GetMatDataFromParticles(COTSParticleList parts, int width, int height)
		{
			int rows, cols;
			cols = width;
			rows = height;
			BYTE* pPixel = new BYTE[width * height];
			Mat cvcopyImg = Mat(rows, cols, CV_8UC1, pPixel);
			cvcopyImg = Scalar::all(0);
			for (auto p : parts)
			{
				auto area = p->GetPixelArea();
				if (area > 100)
				{
					auto fea = p->GetFeature();
					auto segs = fea->GetSegmentsList();
					for (auto seg : segs)
					{
						auto row = seg->GetHeight();
						for (int i = 0; i < seg->GetLength(); i++)
						{
							auto col = seg->GetStart() + i;
							cvcopyImg.at<char>(row, col) = 255;
						}

					}
			    }
			}



			return cvcopyImg;
		}
	
	
	}
	

	

	

	COTSImageProcess::COTSImageProcess(COTSImageProcessParamPtr a_pImageProcessParam)
	{
		m_imageProcessParam = a_pImageProcessParam;
	}


	COTSImageProcess::~COTSImageProcess()
	{
	}

	BOOL COTSImageProcess::RemoveBGByCVconnectivities(CBSEImgPtr inBSEImg, double a_pixelSize,  COTSFieldDataPtr m_pFieldData)
	{

		ASSERT(m_pFieldData);

		ASSERT(inBSEImg);

		ASSERT(m_imageProcessParam);


		int nWidthImg = inBSEImg->GetWidth();
		int nHeightImg = inBSEImg->GetHeight();
		m_pFieldData->Width = nWidthImg;
		m_pFieldData->Height = nHeightImg;
		long nImgSize = nWidthImg * nHeightImg;


		BYTE* pSrcImg = inBSEImg->GetImageDataPointer();

		BYTE* pTempImg = new BYTE[nImgSize];

		CRect r = CRect(0, 0, nWidthImg, nHeightImg);
		CBSEImgPtr imgNoBGBinary = CBSEImgPtr(new CBSEImg(r));
		
		BinaryProcess(inBSEImg, m_imageProcessParam, imgNoBGBinary/*, nNumParticle*/);




		BYTE* pPixel = imgNoBGBinary->GetImageDataPointer();




		long nEigenGrayStart = m_imageProcessParam->GetParticleGray().GetStart();
		long nEigenGrayEnd = m_imageProcessParam->GetParticleGray().GetEnd();


		

			// get the area image	

			Mat cvcopyImg = Mat(nHeightImg, nWidthImg, CV_8UC1, pPixel);
		
			Mat labels = Mat::zeros(cvcopyImg.size(), CV_32S);
			Mat  stats, centroids;
			int number = connectedComponentsWithStats(cvcopyImg, labels, stats, centroids, 8, CV_32S);
			
			double rMin = m_imageProcessParam->GetIncArea().GetStart()/2.0;
			double rMax = m_imageProcessParam->GetIncArea().GetEnd()/2.0;
			double partAreaMin = rMin * rMin * 3.14159;
			double partAreaMax = rMax * rMax * 3.14159;
		

			COTSParticleList listParticleOut;
			for (size_t i = 1; i < number; i++)
			{
				int center_x = centroids.at<double>(i, 0);
				  int center_y = centroids.at<double>(i, 1);
				  //矩形边框
				  int x = stats.at<int>(i, CC_STAT_LEFT);
				  int y = stats.at<int>(i, CC_STAT_TOP);
				  int w = stats.at<int>(i, CC_STAT_WIDTH);
				  int h = stats.at<int>(i, CC_STAT_HEIGHT);
				  int area = stats.at<int>(i, CC_STAT_AREA);

				  double actualArea = area * a_pixelSize * a_pixelSize;

				  if (actualArea >= partAreaMin && actualArea < partAreaMax)
				  {
				 
					Rect rectMax = Rect(x, y, w, h);
					
				
					
					Mat  rectROI = labels(rectMax).clone();
					Mat imageROI = Mat::zeros(rectMax.size(), cvcopyImg.type());

					//exclude the point which  intersect into this bounding box but is not in this group.
					int label = i;
					for (int row = 0; row < rectROI.rows; row++)
					{
						for (int col = 0; col < rectROI.cols; col++)
						{
							int v = rectROI.at<int>(row, col);

							if (v == label)
							{
								imageROI.at<uchar>(row, col) = 255;//set the value to 255,so we won't consider other pixel when we find segment and feature in this ROI.
							}
						

						}

					}
					COTSParticleList roiParts;
					
					if (GetOneParticleFromROI(rectMax.x, rectMax.y, rectMax.width, rectMax.height, imageROI.data, roiParts))
					{
						if (roiParts.size() > 0)
						{
							COTSParticlePtr roiPart = roiParts[0];//we will find only one part in the roi.
							roiPart->SetXRayPos(CPoint(center_x, center_y));
							CRect r = CRect(x, y, x + w, y + h);
							roiPart->SetParticleRect(r);
							roiPart->SetActualArea(actualArea);
							roiPart->SetPixelArea(area);
							listParticleOut.push_back(roiPart);
						}


					}
					
					
				  }

			}
			int  nTagId;
			
			COTSParticleList listParticleSearched;
			for (auto pParticle : listParticleOut)
			{
				COTSFeaturePtr pFeature = pParticle->GetFeature();
				COTSSegmentsList listSegment = pFeature->GetSegmentsList();

				long nPixelNum = 0;
				long eigenPixelNum = 0;
				long nPixelAll = 0;
				int nStartS = 0;
				int nHeightS = 0;
				int nLengthS = 0;
				for (auto pSegment : listSegment)
				{
	
					// get particle average gray
					nStartS = pSegment->GetStart();
					nHeightS = pSegment->GetHeight();
					nLengthS = pSegment->GetLength();
					nPixelNum += (long)nLengthS;


					for (unsigned int i = 0; i < nLengthS; i++)
					{
			

						long nValueTemp = (long)*(pSrcImg + nHeightS * nWidthImg + nStartS + i);
						nPixelAll += nValueTemp;
						if (nValueTemp >= nEigenGrayStart && nValueTemp <= nEigenGrayEnd)
						{
							eigenPixelNum += 1;
						}
					}
				}
				double r = m_imageProcessParam->GetIncArea().GetStart() / 2;
				double minPartArea = 3.14159 * r * r;
				if (eigenPixelNum* a_pixelSize*a_pixelSize > minPartArea)
				{
					BYTE nAveGray = (BYTE)(nPixelAll / nPixelNum);
					pParticle->SetAveGray(nAveGray);
					auto fieldOTSRect = m_pFieldData->GetOTSRect();
					CPoint leftTop = fieldOTSRect.GetTopLeft();

					CRect rectInSinglefld = pParticle->GetParticleRect();
					CPoint OTSLeftTop = CPoint(leftTop.x + rectInSinglefld.left * a_pixelSize, leftTop.y - rectInSinglefld.top * a_pixelSize);
					CPoint OTSRightBottom = CPoint(leftTop.x + rectInSinglefld.right * a_pixelSize, leftTop.y - rectInSinglefld.bottom * a_pixelSize);

					COTSRect recInOTSCord = COTSRect(OTSLeftTop, OTSRightBottom);
					pParticle->SetOTSRect(recInOTSCord);
					listParticleSearched.push_back(pParticle);
				}
				

			}
			m_pFieldData->SetParticleList(listParticleSearched);
		


		delete[]pTempImg;


		return TRUE;
	}
	BOOL COTSImageProcess::GetParticlesBySpecialGrayRange(CBSEImgPtr a_pBSEImg, CIntRangePtr a_grayRange,CDoubleRangePtr a_diameterRange,double a_pixelSize, COTSFieldDataPtr m_pFieldData)
	{
		ASSERT(m_pFieldData);

		ASSERT(a_pBSEImg);

		ASSERT(a_grayRange);


		int nWidthImg = a_pBSEImg->GetWidth();
		int nHeightImg = a_pBSEImg->GetHeight();
		m_pFieldData->Width = nWidthImg;
		m_pFieldData->Height = nHeightImg;
		long nImgSize = nWidthImg * nHeightImg;


		BYTE* pSrcImg = a_pBSEImg->GetImageDataPointer();

		BYTE* pTempImg = new BYTE[nImgSize];
		CRect r = CRect(0, 0, nWidthImg, nHeightImg);

		CBSEImgPtr imgNoBGBinary = CBSEImgPtr(new CBSEImg(r));
		
		long nNumParticle = 0;

		GetSpecialGrayRangeImage(a_pBSEImg, *a_grayRange, imgNoBGBinary, nNumParticle);



		BYTE* pPixel = imgNoBGBinary->GetImageDataPointer();



		if (nNumParticle == 0)
		{

			COTSParticleList listParticleEmpty;
			listParticleEmpty.clear();
			m_pFieldData->SetParticleList(listParticleEmpty);

		
		}
		else
		{

			// get the area image


			Mat cvcopyImg = Mat(nHeightImg, nWidthImg, CV_8UC1, pPixel);
		
		
			Mat labels = Mat::zeros(cvcopyImg.size(), CV_32S);
			Mat  stats, centroids;
			int number = connectedComponentsWithStats(cvcopyImg, labels, stats, centroids, 8, CV_32S);

			double rStart = a_diameterRange->GetStart() / 2.0;
			double rEnd = a_diameterRange->GetEnd() / 2.0;
			double areaStart = rStart * rStart * 3.14159;
			double areaEnd = rEnd * rEnd * 3.14159;

			COTSParticleList listParticleOut;
			for (size_t i = 1; i < number; i++)
			{
				int center_x = centroids.at<double>(i, 0);
				int center_y = centroids.at<double>(i, 1);
				//矩形边框
				int x = stats.at<int>(i, CC_STAT_LEFT);
				int y = stats.at<int>(i, CC_STAT_TOP);
				int w = stats.at<int>(i, CC_STAT_WIDTH);
				int h = stats.at<int>(i, CC_STAT_HEIGHT);
				int area = stats.at<int>(i, CC_STAT_AREA);

				double actualArea = area * a_pixelSize * a_pixelSize;

				

				if (actualArea >= areaStart && actualArea < areaEnd)
				{
					Rect rectMax = Rect(x, y, w, h);



					Mat  rectROI = labels(rectMax).clone();
					Mat imageROI = Mat::zeros(rectMax.size(), cvcopyImg.type());

					//exclude the point which  intersect into this bounding box but is not in this group.
					int label = i;
					for (int row = 0; row < rectROI.rows; row++)
					{
						for (int col = 0; col < rectROI.cols; col++)
						{
							int v = rectROI.at<int>(row, col);

							if (v == label)
							{
								imageROI.at<uchar>(row, col) = 255;
							}


						}

					}

					COTSParticleList roiParts;
					if (!GetOneParticleFromROI(rectMax.x, rectMax.y, rectMax.width, rectMax.height, imageROI.data, roiParts))
					{
						continue;
					}

					if (roiParts.size() > 0)
					{
						COTSParticlePtr roiPart = roiParts[0];
						roiPart->SetXRayPos(CPoint(center_x, center_y));
						CRect r = CRect(x, y, x + w, y + h);
						roiPart->SetParticleRect(r);
						roiPart->SetActualArea(actualArea);
						roiPart->SetPixelArea(area);
						listParticleOut.push_back(roiPart);

					}
				
				}
				

			}



			// form a image only have particles on
			//COTSSegmentsList listImage;
			
			for (auto pParticle : listParticleOut)
			{
				int area = pParticle->GetActualArea();
				double pActualArea = area ;
			
					COTSFeaturePtr pFeature = pParticle->GetFeature();
					COTSSegmentsList listSegment = pFeature->GetSegmentsList();

					long nPixelNum = 0;
					long nPixelAll = 0;
					int nStartS = 0;
					int nHeightS = 0;
					int nLengthS = 0;
					for (auto pSegment : listSegment)
					{

						// get particle average gray
						nStartS = pSegment->GetStart();
						nHeightS = pSegment->GetHeight();
						nLengthS = pSegment->GetLength();
						nPixelNum += (long)nLengthS;


						for (unsigned int i = 0; i < nLengthS; i++)
						{
							long nValueTemp = (long)*(pSrcImg + nHeightS * nWidthImg + nStartS + i);
							nPixelAll += nValueTemp;
						}
					}
					BYTE nAveGray = (BYTE)(nPixelAll / nPixelNum);
					pParticle->SetAveGray(nAveGray);
					auto fieldOTSRect = m_pFieldData->GetOTSRect();
					CPoint leftTop = fieldOTSRect.GetTopLeft();

					CRect rectInSinglefld = pParticle->GetParticleRect();
					CPoint OTSLeftTop = CPoint(leftTop.x + rectInSinglefld.left * a_pixelSize, leftTop.y - rectInSinglefld.top * a_pixelSize);
					CPoint OTSRightBottom = CPoint(leftTop.x + rectInSinglefld.right * a_pixelSize, leftTop.y - rectInSinglefld.bottom * a_pixelSize);

					COTSRect recInOTSCord = COTSRect(OTSLeftTop, OTSRightBottom);
					pParticle->SetOTSRect(recInOTSCord);

			}

			m_pFieldData->SetParticleList(listParticleOut);
		}


		delete[]pTempImg;


		return TRUE;
	}

	BOOL COTSImageProcess::SplitFieldImageIntoMatricsParticle(CBSEImgPtr fieldImg, double a_PixelSize, COTSFieldDataPtr outFldData)
	{
		int nWidthImg = fieldImg->GetWidth();
		int nHeightImg = fieldImg->GetHeight();
		outFldData->Width = nWidthImg;
		outFldData->Height = nHeightImg;

		CRect r = CRect(0, 0, nWidthImg, nHeightImg);
		CBSEImgPtr imgNoBGBinary = CBSEImgPtr(new CBSEImg(r));
	

		double xrayStep = m_imageProcessParam->GetXrayStep();
		xrayStep = xrayStep / a_PixelSize;
		auto rect = fieldImg->GetImageRect();


		std::vector<CPoint> matrixPs;

		int colnum = ceil((double)rect.Width() / xrayStep + 0.5);
		if (colnum % 2 == 0) colnum += 1;//let the number to be an odd number.Then we can make the middle point in the center of the particle exactly.
		int rownum = ceil((double)rect.Height() / xrayStep + 0.5);
		if (rownum % 2 == 0) rownum += 1;

		CPoint theFirst = CPoint(rect.left - (colnum * xrayStep - rect.Width()) / 2 + xrayStep / 2, rect.top - (rownum * xrayStep - rect.Height()) / 2 + xrayStep / 2);


		for (int i = 0; i < rownum; i++)
		{
			for (int j = 0; j < colnum; j++)
			{
				double x = (double)theFirst.x + (double)j * xrayStep;
				double y = (double)theFirst.y + (double)i * xrayStep;
				CPoint thePoint = CPoint(x, y);


				matrixPs.push_back(thePoint);


			}
		}		
		int i = 0;
		COTSParticleList matrixParts;
		for (auto point : matrixPs)
		{
			COTSParticlePtr part = COTSParticlePtr(new COTSParticle());


			COTSSegmentsList segs;
			for (int i = 0; i < xrayStep; i++)
			{
				COTSSegmentPtr seg = COTSSegmentPtr(new COTSSegment());
				seg->SetStart(point.x - xrayStep / 2);
				seg->SetLength(xrayStep);
				seg->SetHeight(point.y - xrayStep / 2 + i);

				//auto originalSegs = segsOnTheSameHeight[seg->GetHeight()];
				int currentH = seg->GetHeight();
				int segStart = seg->GetStart();
				int segEnd = seg->GetEnd();

				if (segStart < 0) segStart = 0;
				if (segEnd >= rect.Width()) segEnd = rect.Width()-1;
				if (currentH >= rect.Height()) continue;
				if (currentH < 0) continue;
				seg->SetStart(segStart);
				seg->SetEnd(segEnd);
				segs.push_back(seg);
			
			}


			if (segs.size() > 0)
			{
				
				part->GetFeature()->SetSegmentsList(segs);
				part->SetActualArea(xrayStep * xrayStep);
				part->CalXRayPos();
				part->SetFieldId(outFldData->GetId());
				part->SetAnalysisId(i);
				matrixParts.push_back(part);
				i++;
			}

		}
		outFldData->SetParticleList(matrixParts);
		

		return true;
		//-----------
	}
	
	BOOL COTSImageProcess::SplitRawImageIntoParticlesByWaterShed(CBSEImgPtr fieldImg, double a_PixelSize, COTSFieldDataPtr outFldData)
	{
		int nWidthImg = fieldImg->GetWidth();
		int nHeightImg = fieldImg->GetHeight();
		outFldData->Width = nWidthImg;
		outFldData->Height = nHeightImg;

		CRect r = CRect(0, 0, nWidthImg, nHeightImg);
		CBSEImgPtr imgNoBGBinary = CBSEImgPtr(new CBSEImg(r));
		
		COTSParticleList listMainParticle;
		BinaryProcess(fieldImg, m_imageProcessParam, imgNoBGBinary);

		
		CDoubleRange areaRng = m_imageProcessParam->GetIncArea();
		GetParticlesFromBinaryImage(imgNoBGBinary, areaRng, a_PixelSize, listMainParticle);

		outFldData->SetParticleList(listMainParticle);
		

		
		CBSEImgPtr noBgImg = CBSEImgPtr(new CBSEImg(fieldImg->GetImageRect()));
		RemoveBackGround(fieldImg, noBgImg);

		BlurImage(noBgImg);

		std::vector<CPoint> seeds;
		seeds.clear();
		FindSeedsByGrayScale(noBgImg, a_PixelSize, seeds);

		Mat matImg = GetMatDataFromBseImg(noBgImg);

		Mat marks(matImg.size(), CV_32S);

		marks = Scalar::all(0);

		for (int index = 0; index < seeds.size(); index++)
		{

			auto p = seeds[index];

			marks.at<INT32>(p.y, p.x) = index + 10;

		}

		Mat imageGray3;

	
		cv::cvtColor(matImg, imageGray3, cv::COLOR_GRAY2RGB);//灰度转换 opencv4
		cv::watershed(imageGray3, marks);


		auto matBinNobg = GetMatDataFromParticles(listMainParticle, nWidthImg, nHeightImg);

		//the first "long" is the index of particle,the second "long" is the row number of the segment, the third "long" is the sequence of the pixel.
		std::map<long, std::map<long, std::set<long>>> pixelMap;

		for (int i = 0; i < marks.rows; i++)
		{
			for (int j = 0; j < marks.cols; j++)
			{
				if (matBinNobg.at<char>(i, j) == 0)
				{
					marks.at<int>(i, j) = 0;
					continue;
				}

				long index = marks.at<int>(i, j);

				if (index == -1)//if this pixel is on the border,then it will be absorbed by it's neighbor.
				{
					if (i + 1 < marks.rows && marks.at<int>(i + 1, j) != -1 && matBinNobg.at<char>(i + 1, j) != 0)
					{
						index = marks.at<int>(i + 1, j);
					}
					else if (j + 1 < marks.cols && marks.at<int>(i, j + 1) != -1 && matBinNobg.at<char>(i, j + 1) != 0)
					{
						index = marks.at<int>(i, j + 1);
					}
					else if (i > 0 && marks.at<int>(i - 1, j) != -1 && matBinNobg.at<char>(i - 1, j) != 0)
					{
						index = marks.at<int>(i - 1, j);
					}
					else if (j > 0 && marks.at<int>(i, j - 1) != -1 && matBinNobg.at<char>(i, j - 1) != 0)
					{
						index = marks.at<int>(i, j - 1);
					}
					marks.at<int>(i, j) = index;
				}

				pixelMap[index][i].insert(j);

			}
		}

		COTSParticleList listParticleOut;
		for (auto it = pixelMap.begin(); it != pixelMap.end(); it++)
		{
			if (it->first == -1) continue;
			if (it->first == 0) continue;
			auto pixelsegs = it->second;
			COTSParticlePtr newGrain = COTSParticlePtr(new COTSParticle());
			COTSFeaturePtr newFea = COTSFeaturePtr(new COTSFeature());
			for (auto pixelseg : pixelsegs)
			{
				auto start = (*pixelseg.second.begin());
				for (auto pix = pixelseg.second.begin(); pix != pixelseg.second.end(); pix++)
				{
					auto curNode = pix;
					++curNode;
					if ((curNode) != pixelseg.second.end())
					{
						if (*curNode - *pix > 1)
						{
							COTSSegmentPtr  seg = COTSSegmentPtr(new COTSSegment());

							seg->SetStart(start);
							seg->SetEnd(*pix);
							seg->SetHeight(pixelseg.first);
							newFea->AddSegment(seg);

							start = *curNode;
						}
					}
					else
					{
						COTSSegmentPtr  seg = COTSSegmentPtr(new COTSSegment());

						seg->SetStart(start);
						seg->SetEnd(*pix);
						seg->SetHeight(pixelseg.first);
						newFea->AddSegment(seg);
					}

				}

			}
			newGrain->SetFeature(newFea);


			newGrain->CalXRayPos();

			CPoint pos = newGrain->GetXRayPos();

			for (auto mp : listMainParticle)
			{
			
				auto fieldOTSRect = outFldData->GetOTSRect();
				CPoint leftTop = fieldOTSRect.GetTopLeft();

				CRect rectInSinglefld = mp->GetParticleRect();
				CPoint OTSLeftTop = CPoint(leftTop.x + rectInSinglefld.left * a_PixelSize, leftTop.y - rectInSinglefld.top * a_PixelSize);
				CPoint OTSRightBottom = CPoint(leftTop.x + rectInSinglefld.right * a_PixelSize, leftTop.y - rectInSinglefld.bottom * a_PixelSize);

				COTSRect recInOTSCord = COTSRect(OTSLeftTop, OTSRightBottom);
				mp->SetOTSRect(recInOTSCord);
				mp->SetFieldId(outFldData->GetId());
				if (mp->IfContain(pos))
				{
					CRect rectInSinglefld = newGrain->GetParticleRect();
					CPoint OTSLeftTop = CPoint(leftTop.x + rectInSinglefld.left * a_PixelSize, leftTop.y - rectInSinglefld.top * a_PixelSize);
					CPoint OTSRightBottom = CPoint(leftTop.x + rectInSinglefld.right * a_PixelSize, leftTop.y - rectInSinglefld.bottom * a_PixelSize);

					COTSRect recInOTSCord = COTSRect(OTSLeftTop, OTSRightBottom);
					newGrain->SetFieldId(mp->GetFieldId());
					newGrain->SetOTSRect(recInOTSCord);
					mp->AddSubParticle(newGrain);
					break;
				}

			}



		}

		return true;
	}
	
	CIntRangePtr COTSImageProcess::CalBackground(CBSEImgPtr m_pBSEImg)
	{
		auto ranges = CalcuGrayLevelRange(m_pBSEImg);

		return ranges[0];
	}
	std::vector<CIntRangePtr> COTSImageProcess::CalcuGrayLevelRange(CBSEImgPtr m_pBSEImg)
	{
		CIntRangePtr pBackground = CIntRangePtr(new CIntRange());
		WORD originChartData[MAXBYTE];
		WORD firstSmoothChart[MAXBYTE];
		WORD secondSmooth[MAXBYTE];
		//1. get chart data
		m_pBSEImg->SetChartData();
	
		memcpy(originChartData, m_pBSEImg->GetBSEChart(), sizeof(WORD) * MAXBYTE);
		originChartData[0] = 0;
		originChartData[254] = 0;

		linearSmooth5(originChartData, firstSmoothChart, MAXBYTE);
		linearSmooth5(firstSmoothChart, secondSmooth, MAXBYTE);
		
		
		//2. get down edge		
		int nLengthEdge = MAXBYTE + 2;

		WORD n_aBSEChart[MAXBYTE + 2];

		memset(n_aBSEChart, 0, sizeof(WORD) * nLengthEdge);

		std::map<long, std::vector<int>> peakMap;// hold all the peaks in this spectrum which are sorted by there area.
		std::vector<int> currentUpSeries;
	
		std::vector<int> currentPeakSeries;

		// make sure the wave begin with up edge and end with down edge
		n_aBSEChart[0] = 0;
		n_aBSEChart[nLengthEdge - 1] = 0;

		memcpy(&n_aBSEChart[1], &secondSmooth, sizeof(WORD) * MAXBYTE);


		int nLengthCom = MAXBYTE + 1;

		// up edge
		for (int i = 0; i < nLengthCom; i++)
		{
			if (n_aBSEChart[i] <= n_aBSEChart[i + 1])//this is a upward edge
			{
				
				if (currentPeakSeries.size() > 0)
				{
					int seriesSize = currentPeakSeries.size();
					long area = 0;
					for (int i = 0; i < seriesSize; i++)
					{
						area = area + n_aBSEChart[currentPeakSeries[i]];
					}
					peakMap[area] = currentPeakSeries;
					currentPeakSeries.clear();
				}

				currentUpSeries.push_back(i + 1);// save all the continuous up edge

			}
			else//this is a downward edge
			{
				// encounter a downward edge means upward edge series end, 
			

				if (currentUpSeries.size() > 0)
				{
					
					currentPeakSeries = currentUpSeries;
					currentUpSeries.clear();
				}

				currentPeakSeries.push_back(i + 1);

			}
		}
		if (currentPeakSeries.size() > 0)
		{
			int seriesSize = currentPeakSeries.size();
			long area = 0;
			for (int i = 0; i < seriesSize; i++)
			{
				area = area + n_aBSEChart[currentPeakSeries[i]];
			}
			peakMap[area] = currentPeakSeries;
			currentPeakSeries.clear();

		}
		std::vector<CIntRangePtr> ranges;
		std::map<long, std::vector<int>>::reverse_iterator it;
		for (it=peakMap.rbegin();it!=peakMap.rend();it++)
		{
			CIntRangePtr pRange = CIntRangePtr(new CIntRange());

			pRange->SetStart(it->second[0]);
			pRange->SetEnd(it->second[it->second.size()-1]);
			ranges.push_back(pRange);

		}
	


		return ranges;
	
	}
	void COTSImageProcess::GetSpecialGrayRangeImage(CBSEImgPtr a_pImgIn, CIntRange a_SpecialGrayRange, CBSEImgPtr a_pBinImgOut, long& foundedPixelNum)
	{
		// the background  pixel will be 0,and the other part will be 255.
		ASSERT(a_pImgIn);


	
		int nWidthImg = a_pImgIn->GetWidth();
		int nHeightImg = a_pImgIn->GetHeight();

		long nImgSize = nWidthImg * nHeightImg;

		BYTE* pTempImg = new BYTE[nImgSize];

		BYTE* pSrcImg = a_pImgIn->GetImageDataPointer();


		BYTE* pPixel = new byte[nImgSize];



		long nBGStart;
		long nBGEnd;
	
		long nNumParticle = 0;
	
			nBGStart = a_SpecialGrayRange.GetStart();
			nBGEnd = a_SpecialGrayRange.GetEnd();
			

			// delete background 
			for (unsigned int i = 0; i < nImgSize; i++)
			{
				if (pSrcImg[i] >= nBGStart && pSrcImg[i] <= nBGEnd)
				{
					pPixel[i] = 255;
					nNumParticle++;
				}
				else
				{
					pPixel[i] = 0;
					
				}
		
			}
		
				BErode3(pPixel, pTempImg, 5, nHeightImg, nWidthImg);
				BDilate3(pTempImg, pPixel, 5, nHeightImg, nWidthImg);
			
			
			
		a_pBinImgOut->SetImageData(pPixel, nWidthImg, nHeightImg);

		foundedPixelNum = nNumParticle;
		delete[] pTempImg;

		return;
	}
	void COTSImageProcess::BinaryProcess(CBSEImgPtr a_pImgIn, COTSImageProcessParamPtr a_pImageProcessParam, CBSEImgPtr a_pBinImgOut/*,long& foundedPixelNum*/)
	{
		// the background  pixel will be 0,and the other part will be 255.
		ASSERT(a_pImgIn);


		ASSERT(a_pImageProcessParam);
		int nWidthImg = a_pImgIn->GetWidth();
		int nHeightImg = a_pImgIn->GetHeight();

		long nImgSize = nWidthImg * nHeightImg;

	

		BYTE* pSrcImg = a_pImgIn->GetImageDataPointer();
	

	    BYTE* pPixel= new BYTE[nImgSize];
		Mat srcImgMat = GetMatDataFromBseImg(a_pImgIn);
		Mat rstMat;
		

		long nBGStart;
		long nBGEnd;
		long nPartStart;
		long nPartEnd;
		long nNumParticle = 0;
		if (a_pImageProcessParam->GetBGRemoveType() == OTS_BGREMOVE_TYPE::MANUAL || a_pImageProcessParam->GetBGRemoveType() == OTS_BGREMOVE_TYPE::WaterShed)
		{
			nBGStart = a_pImageProcessParam->GetBGGray().GetStart();
			nBGEnd = a_pImageProcessParam->GetBGGray().GetEnd();
			nPartStart = a_pImageProcessParam->GetParticleGray().GetStart();
			nPartEnd = a_pImageProcessParam->GetParticleGray().GetEnd();

			
			
			
			//CVRemoveBG(srcImgMat, rstMat, nBGStart, nBGEnd);
			RemoveBG_old(srcImgMat, rstMat, nBGStart, nBGEnd, nNumParticle);

			pPixel = rstMat.data;
		

			
		}
		else
		{
		
		

			switch (a_pImageProcessParam->GetAutoBGRemoveType())
			{
			case OTS_AUTOBGREMOVE_TYPE::DOWNWARD:
				
				//RemoveBG_old(srcImgMat, rstMat, 0, nBGEnd, nNumParticle);
				//CVRemoveBG(srcImgMat, rstMat, 0, nBGEnd, nNumParticle);
				AutoRemove_background_OTS(srcImgMat, rstMat, 1);
				break;
			case OTS_AUTOBGREMOVE_TYPE::UPWARD:
				
				//RemoveBG_old(srcImgMat, rstMat, nBGStart, 255, nNumParticle);
				//CVRemoveBG(srcImgMat, rstMat, nBGStart, 255, nNumParticle);
				AutoRemove_background_OTS(srcImgMat, rstMat, 0);
				break;
			case OTS_AUTOBGREMOVE_TYPE::MIDDLE:
				
				//RemoveBG_old(srcImgMat, rstMat, nBGStart, nBGEnd, nNumParticle);
				//CVRemoveBG(srcImgMat, rstMat, nBGStart, nBGEnd, nNumParticle);
				AutoRemove_background_OTS(srcImgMat, rstMat, 2);
				break;
			default:
				break;
			}
			pPixel = rstMat.data;
			
		
		}
		a_pBinImgOut->SetImageData(pPixel,nWidthImg,nHeightImg);
		
		
		
		
		return ;
	}
	void COTSImageProcess::RemoveBackGround(CBSEImgPtr a_pImgIn,  CBSEImgPtr a_pBinImgOut)
	{
		// the background  pixel will be 0,and the other part will be 255.
		ASSERT(a_pImgIn);


		ASSERT(m_imageProcessParam);
		int nWidthImg = a_pImgIn->GetWidth();
		int nHeightImg = a_pImgIn->GetHeight();

		long nImgSize = nWidthImg * nHeightImg;

		BYTE* pTempImg = new BYTE[nImgSize];

		BYTE* pSrcImg = a_pImgIn->GetImageDataPointer();


		BYTE* pPixel = new byte[nImgSize];



		long nBGStart;
		long nBGEnd;
		long nPartStart;
		long nPartEnd;
		long nNumParticle = 0;
		if (m_imageProcessParam->GetBGRemoveType() == OTS_BGREMOVE_TYPE::MANUAL|| m_imageProcessParam->GetBGRemoveType() == OTS_BGREMOVE_TYPE::WaterShed)
		{
			nBGStart = m_imageProcessParam->GetBGGray().GetStart();
			nBGEnd = m_imageProcessParam->GetBGGray().GetEnd();
			nPartStart = m_imageProcessParam->GetParticleGray().GetStart();
			nPartEnd = m_imageProcessParam->GetParticleGray().GetEnd();

			// delete background 
			for (unsigned int i = 0; i < nImgSize; i++)
			{
				if (pSrcImg[i] >= nBGStart && pSrcImg[i] <= nBGEnd)
				{
					pPixel[i] = 0;
				}
				else
				{
					
						pPixel[i] = pSrcImg[i];
						nNumParticle++;
					

				}
				if (pSrcImg[i]<nPartStart || pSrcImg[i]>nPartEnd)
				{
					pPixel[i] = 0;
				}
			}
			




		}
		
		a_pBinImgOut->SetImageData(pPixel, nWidthImg, nHeightImg);

		
		delete[] pTempImg;

		return;
	}
	BOOL COTSImageProcess::GetParticles(long left, long top, long a_nWidth, long a_nHeight, const BYTE* a_pPixel, COTSParticleList& a_listParticles)
	{
		ASSERT(a_pPixel);
		if (!a_pPixel)
		{
			return FALSE;
		}

		//a_listParticles.clear();
		COTSParticleList findedParts;
		COTSSegmentsList listSegment;
		listSegment.clear();


		//1. get segment line by line	
		if (!GetSegmentList(left, top, a_nWidth, a_nHeight, a_pPixel, listSegment))
		{
			return FALSE;
		}
		if ((int)listSegment.size() == 0)
		{

			return FALSE;
		}


		//2. save the temp feature
		COTSFeatureList listFeature;
		listFeature.clear();
		if (!GetFeatureList(listSegment, listFeature))//get every feature for all the particle,the complete feature.
		{
			return FALSE;
		}


		if ((int)listFeature.size() == 0)
		{
			return FALSE;
		}

		/*COTSParticleList listParticles;
		listParticles.clear();*/
		if (!ChangeFeaturelist(listFeature, findedParts))
		{
			return FALSE;
		}
		for (auto f : findedParts)
		{
			a_listParticles.push_back(f);
		}

		return TRUE;
	}
	BOOL COTSImageProcess::GetOneParticleFromROI(long left, long top, long a_nWidth, long a_nHeight, const BYTE* a_pPixel, COTSParticleList& a_listParticles)
	{
		ASSERT(a_pPixel);
		if (!a_pPixel)
		{
			return FALSE;
		}

		//a_listParticles.clear();
		COTSParticleList findedParts;
		COTSSegmentsList listSegment;
		listSegment.clear();


		//1. get segment line by line	
		if (!GetSegmentList(left, top, a_nWidth, a_nHeight, a_pPixel, listSegment))
		{
			return FALSE;
		}
		if ((int)listSegment.size() == 0)
		{

			return FALSE;
		}


		//2. save the temp feature
		COTSFeatureList listFeature;
		listFeature.clear();
		COTSFeaturePtr fea = COTSFeaturePtr(new COTSFeature());
		
		fea->SetSegmentsList(listSegment);
		
		listFeature.push_back(fea);

		if ((int)listFeature.size() == 0)
		{
			return FALSE;
		}

	
		if (!ChangeFeaturelist(listFeature, findedParts))
		{
			return FALSE;
		}
		for (auto f : findedParts)
		{
			a_listParticles.push_back(f);
		}

		return TRUE;
	}
	BOOL COTSImageProcess::GetSegmentList(long left, long top, long a_nWidth, long a_nHeight, const BYTE* a_pPixel, COTSSegmentsList& a_listSegments)
	{
		ASSERT(a_pPixel);
		

		long nImgSize = a_nWidth * a_nHeight;

		a_listSegments.clear();

		//1. get segment line by line
		long nLine, nm, nn;
		long nStart = 0, nLength = 0;
		for (nLine = 0; nLine < a_nHeight; nLine++)
		{
			for (nm = 0; nm < a_nWidth; nm += (nLength + 1))
			{
				nLength = 0;
				// get start
				if (*(a_pPixel + nLine * a_nWidth + nm) != 0)
				{
					nStart = nm;
					nLength++;
					//get length
					for (nn = nm + 1; nn < a_nWidth; nn++)
					{
						// check if segment is over, break
						if (nLength != 0)
						{
							if (*(a_pPixel + nLine * a_nWidth + nn) == 0)
								break;
						}

						if (*(a_pPixel + nLine * a_nWidth + nn) != 0)
						{
							nLength++;
						}
					}

					// generate segment
					COTSSegmentPtr pSegment = COTSSegmentPtr(new COTSSegment(nLine + top, nStart + left, nLength));
					a_listSegments.push_back(pSegment);

				}
				else
				{
					continue;
				}
			}
		}

		if ((int)a_listSegments.size() == 0)
		{
			//LogErrorTrace(__FILE__, __LINE__, _T("no particle is found."));
			return FALSE;
		}

		return TRUE;
	}
	BOOL COTSImageProcess::GetFeatureList(COTSSegmentsList& a_listSegments, COTSFeatureList& a_listFeatures)
	{
		COTSSegmentsList listSegmentNew;
		std::map<long, COTSSegmentsList > mapOneLineSegments;
		for each (auto s in a_listSegments)
		{
			mapOneLineSegments[s->GetHeight()].push_back(s);//sorting all the segments base on the line number.
		}
		std::map<long, COTSSegmentsList >::iterator lineItr = mapOneLineSegments.begin();//find the highest line
		while (lineItr != mapOneLineSegments.end())
		{
			for (auto s = lineItr->second.begin(); s < lineItr->second.end(); )//find  one segment of this line.
			{
				COTSSegmentPtr bottomSeg = *s;
				listSegmentNew.clear();
				listSegmentNew.push_back(*s);
				s = lineItr->second.erase(s);
				std::map<long, COTSSegmentsList >::iterator tempItr = lineItr;
				tempItr++;
				for (; tempItr != mapOneLineSegments.end(); tempItr++)//find all other lines of segments
				{
					if (tempItr->first - bottomSeg->GetHeight() > 1)
					{
						break;
					}

					for (auto nextLineSegment = tempItr->second.begin(); nextLineSegment < tempItr->second.end();)//find next line's all segments
					{
						if (((*nextLineSegment)->GetStart() - (bottomSeg->GetStart() + bottomSeg->GetLength())) > 1)
						{
							break;
						}

						if (bottomSeg->UpDownConection(**nextLineSegment))
						{
							listSegmentNew.push_back(*nextLineSegment);
							bottomSeg = *nextLineSegment;
							nextLineSegment = tempItr->second.erase(nextLineSegment);
							break;
						}

						if (tempItr->second.size() > 0)
						{
							nextLineSegment++;
						}
						else
						{
							break;
						}
					}

				}
				COTSFeaturePtr pFeature = COTSFeaturePtr(new COTSFeature());
				pFeature->SetSegmentsList(listSegmentNew);
				//check if this new feature is connected with other found feature.
				COTSSegmentPtr topSeg = listSegmentNew[0];//find the toppest segment of this new feature.
				COTSSegmentPtr bottomSegment = listSegmentNew[listSegmentNew.size() - 1];//find the lowest segment of this new feature.

				bool haveMerged = false;
				for each (auto f in a_listFeatures)
				{
					for (auto seg : f->GetSegmentsList())
					{
						if (bottomSegment->UpDownConection(*seg) || topSeg->UpDownConection(*seg))
						{
							COTSSegmentsList segs = f->GetSegmentsList();
							for (auto s : listSegmentNew)
							{
								segs.push_back(s);

							}

							f->SetSegmentsList(segs);
							haveMerged = true;
							break;
						}
					}
					if (haveMerged)
					{
						break;
					}

				}
				if (!haveMerged)
				{
					a_listFeatures.push_back(pFeature);
				}

				if (lineItr->second.size() == 0)
				{
					break;
				}
			}
			lineItr++;
		}
		return true;
	}

	BOOL COTSImageProcess::ChangeFeaturelist(COTSFeatureList& a_listFeatures, COTSParticleList& a_listParticle)
	{
	
		for (auto pFeature : a_listFeatures)
		{
			COTSParticlePtr pParticle = COTSParticlePtr(new COTSParticle());
			pParticle->SetFeature(pFeature);
	
			a_listParticle.push_back(pParticle);
		}

		if ((int)a_listParticle.size() == 0)
		{
			return FALSE;
		}

		return TRUE;
	}

	BOOL COTSImageProcess::GetParticlesFromBinaryImage(CBSEImgPtr m_pBSEImg, CDoubleRange a_diameterRange, double a_pixelSize, COTSParticleList& listParticleOut)
	{
		int height, width;
		height = m_pBSEImg->GetHeight();
		width = m_pBSEImg->GetWidth();

		BYTE* pPixel = m_pBSEImg->GetImageDataPointer();
		Mat nobgMat = Mat(height, width, CV_8UC1, pPixel);


		Mat labels = Mat::zeros(nobgMat.size(), CV_32S);
		Mat  stats, centroids;
		int number = connectedComponentsWithStats(nobgMat, labels, stats, centroids, 8, CV_32S);

		double rStart = a_diameterRange.GetStart()/2.0;
		double rEnd = a_diameterRange.GetEnd()/2.0;
		double areaStart =  rStart* rStart * 3.14159;
		double areaEnd =  rEnd* rEnd * 3.14159;

		//COTSParticleList listParticleOut;
		for (size_t i = 1; i < number; i++)
		{
			int center_x = centroids.at<double>(i, 0);
			int center_y = centroids.at<double>(i, 1);
			//矩形边框
			int x = stats.at<int>(i, CC_STAT_LEFT);
			int y = stats.at<int>(i, CC_STAT_TOP);
			int w = stats.at<int>(i, CC_STAT_WIDTH);
			int h = stats.at<int>(i, CC_STAT_HEIGHT);
			int area = stats.at<int>(i, CC_STAT_AREA);

			double actualArea = area * a_pixelSize * a_pixelSize;



			if (actualArea >= areaStart && actualArea < areaEnd)
			{
				Rect rectMax = Rect(x, y, w, h);

				Mat  rectROI = labels(rectMax).clone();
				Mat imageROI = Mat::zeros(rectMax.size(), nobgMat.type());

				//exclude the point which  intersect into this bounding box but is not in this group.
				int label = i;
				for (int row = 0; row < rectROI.rows; row++)
				{
					for (int col = 0; col < rectROI.cols; col++)
					{
						int v = rectROI.at<int>(row, col);

						if (v == label)
						{
							imageROI.at<uchar>(row, col) = 255;
						}


					}

				}

				COTSParticleList roiParts;
				if (!GetOneParticleFromROI(rectMax.x, rectMax.y, rectMax.width, rectMax.height, imageROI.data, roiParts))
				{
					continue;
				}

				if (roiParts.size() > 0)
				{
					COTSParticlePtr roiPart = roiParts[0];
					roiPart->SetXRayPos(CPoint(center_x, center_y));
					CRect r = CRect(x, y, x + w, y + h);
					roiPart->SetParticleRect(r);
					roiPart->SetActualArea(actualArea);
					roiPart->SetPixelArea(area);
					listParticleOut.push_back(roiPart);

				}



			}


		}

		return TRUE;
	}
	void  COTSImageProcess::findPeakAndValley(const vector<int>& v, vector<int>& peakPositions, vector<int>& valleyPositions)
	{
		vector<int> diff_v(v.size() - 1, 0);
		// 计算V的一阶差分和符号函数trend
		for (vector<int>::size_type i = 0; i != diff_v.size(); i++)
		{
			if (v[i + 1] - v[i] > 0)
				diff_v[i] = 1;
			else if (v[i + 1] - v[i] < 0)
				diff_v[i] = -1;
			else
				diff_v[i] = 0;
		}
		// 对Trend作了一个遍历
		for (int i = diff_v.size() - 1; i >= 0; i--)
		{
			if (diff_v[i] == 0 && i == diff_v.size() - 1)
			{
				diff_v[i] = 1;
			}
			else if (diff_v[i] == 0)
			{
				if (diff_v[i + 1] >= 0)
					diff_v[i] = 1;
				else
					diff_v[i] = -1;
			}
		}

		for (vector<int>::size_type i = 0; i != diff_v.size() - 1; i++)
		{
			if (diff_v[i + 1] - diff_v[i] == -2)
				peakPositions.push_back(i + 1);
			if (diff_v[i + 1] - diff_v[i] == 2)
				valleyPositions.push_back(i + 1);
		}
	}
	std::vector<CIntRange> COTSImageProcess::GetValidGreyLevelRanges(CBSEImgPtr pBSEImg)
	{


		WORD originChartData[MAXBYTE];


		//1. get chart data
		pBSEImg->SetChartData();

		memcpy(originChartData, pBSEImg->GetBSEChart(), sizeof(WORD) * MAXBYTE);
		originChartData[0] = 0;

		std::vector<int> greyChart;
		for (int i = 0; i < MAXBYTE; i++)
		{
			greyChart.push_back((int)originChartData[i]);
		}
		std::vector<int> peaks;
		std::vector<int> valleys;
		findPeakAndValley(greyChart, peaks, valleys);

		if ((peaks.size() - valleys.size()) == 1)
		{
			valleys.push_back(254);
		}




		std::map<int, MyGreyPeak*> peakMap;// hold all the peaks in this spectrum which are sorted by there area.

		MyGreyPeak* lastRng = NULL;
		MyGreyPeak* firstRng = NULL;

		for (int i = 0; i < peaks.size(); i++)
		{
			MyGreyPeak* gr = new MyGreyPeak();
			if (i == 0)
			{
				gr->start = 0;
			}
			else
			{
				gr->start = valleys[i - 1];
			}

			gr->end = valleys[i];
			for (int j = gr->start; j < gr->end; j++)
			{
				gr->area += originChartData[j];
			}

			gr->peakPos = peaks[i];
			gr->peakValue = originChartData[peaks[i]];
			gr->valleyValue = valleys[i];
			if (firstRng == NULL)
			{
				firstRng = gr;
				lastRng = gr;
			}
			else
			{
				gr->preRng = lastRng;
				lastRng->nextRng = gr;
			}

			lastRng = gr;

			peakMap[gr->start] = gr;

		}
		int validGreyLevelNum = 0;//statistic the none zero grey level number
		long partArea = 0;
		for (int i = 0; i < 255; i++)
		{
			if (originChartData[i] != 0)
			{
				validGreyLevelNum += 1;
				partArea += originChartData[i];
			}
		}


		auto curRange = firstRng;
		int maxPeakValue = curRange->peakValue;
		while (curRange != NULL)
		{
			if (curRange->peakValue > maxPeakValue)
			{
				maxPeakValue = curRange->peakValue;
			}
			curRange = curRange->nextRng;
		}
		curRange = firstRng;
		while (curRange != NULL)
		{
			curRange->maxPeakvalue = maxPeakValue;
			curRange->allNoneZeroGreyLevelNum = validGreyLevelNum;
			curRange = curRange->nextRng;
		}
		MyGreyPeak* firstRng1=NULL;

		MyGreyPeak::MergeGreyPeaks(firstRng,firstRng1);

		std::vector<CIntRange> ranges;
		auto curRange1 = firstRng1;
		while (curRange1 != NULL)
		{

			CIntRange pRange;
			pRange.SetStart(curRange1->start);
			pRange.SetEnd(curRange1->end);
			ranges.push_back(pRange);
			curRange1 = curRange1->nextRng;
		}

		return ranges;

	}
	BOOL COTSImageProcess::FindSeedsByGrayScale(CBSEImgPtr fieldImg, double a_PixelSize, std::vector<CPoint>& OutSeeds)
	{
		ASSERT(fieldImg);


		ASSERT(m_imageProcessParam);
		int nWidthImg = fieldImg->GetWidth();
		int nHeightImg = fieldImg->GetHeight();

		long nImgSize = nWidthImg * nHeightImg;

		std::vector<CIntRange> validPeaks = GetValidGreyLevelRanges(fieldImg);



		CRect r = CRect(0, 0, nWidthImg, nHeightImg);

		CBSEImgPtr imgNoBGBinary = CBSEImgPtr(new CBSEImg(r));

		for (int i = 0; i < validPeaks.size(); i++)
		{
			CIntRange pRange = validPeaks[i];
			if (pRange.GetStart() == 0)//prevent the 0 grey pixels being recognized as particle.
			{
				pRange.SetStart(3);
			}



			long nNumParticle = 0;

			GetSpecialGrayRangeImage(fieldImg, pRange, imgNoBGBinary, nNumParticle);



			if (nNumParticle == 0)
			{

				continue;


			}
			else
			{

				// get the area image


				Mat matImg = GetMatDataFromBseImg(imgNoBGBinary);
				Mat labels = Mat::zeros(matImg.size(), CV_32S);
				Mat  stats, centroids;
				int number = connectedComponentsWithStats(matImg, labels, stats, centroids, 8, CV_32S);
				auto incsize = m_imageProcessParam->GetIncArea();
				double rStart = incsize.GetStart() / 2.0;
				double rEnd = incsize.GetEnd() / 2.0;
				double pixelarea = a_PixelSize * a_PixelSize;
				double areaStart = rStart* rStart * 3.14159*pixelarea;
				double areaEnd =  rEnd* rEnd * 3.14159*pixelarea;
				/*double areaStart = 0;
				double areaEnd = 1000;*/


				for (size_t i = 1; i < number; i++)
				{
					int center_x = centroids.at<double>(i, 0);
					int center_y = centroids.at<double>(i, 1);

					int area = stats.at<int>(i, CC_STAT_AREA);

					double actualArea = area * a_PixelSize * a_PixelSize;



					if (actualArea >= areaStart && actualArea < areaEnd)
					{


						CPoint seed = CPoint(center_x, center_y);
						OutSeeds.push_back(seed);
					}
				}



			}

		}




		return true;
	}
	
	 BOOL COTSImageProcess::CalcuParticleImagePropertes(COTSParticlePtr a_pOTSPart, double a_PixelSize)
	{	
		//--------- convert this particle data to image data,construct an image only with this particle.------
		const int nExpand_Size = 3;
		const int nWhiteColor = 0;
		const int nThickness = 1;
		// lineType Type of the line
		const int nLineType = 8;
		// get rectangle of the particle
		CRect rect = a_pOTSPart->GetParticleRect();
		if (a_pOTSPart->GetActualArea() < 30 * a_PixelSize)// the particle is too small that openCV can't calculate a width value of it. Then we take the upright rect of the particle as it's minArea rect.
		{
			double w = 0, h = 0;
			w = (double)rect.Width()*a_PixelSize;
			h = (double)rect.Height()*a_PixelSize;
			a_pOTSPart->SetDMax(MAX(w, h));
			a_pOTSPart->SetDMin(MIN(w, h));
			a_pOTSPart->SetDMean((w + h) / 2);
			a_pOTSPart->SetFeretDiameter((w + h) / 2);
			a_pOTSPart->SetDElong(MAX(w, h));
			a_pOTSPart->SetPerimeter((w+h)*2);
			a_pOTSPart->SetDPerp(MIN(w, h));
			a_pOTSPart->SetDInscr(MIN(w, h));
			a_pOTSPart->SetMinWidth(w);
			a_pOTSPart->SetMinHeight(h);
			return true;
		}
		// calculate the particle image data size, expand 3 pixel at the edge
		Mat particleImage = Mat::zeros(rect.Height() + nExpand_Size , rect.Width() + nExpand_Size , CV_8U);
		// get the segment list
		COTSSegmentsList listSegment = a_pOTSPart->GetFeature()->GetSegmentsList();
		for (auto pSegment : listSegment)
		{
			int nStart = pSegment->GetStart() - rect.left + nExpand_Size;
			int nEnd = pSegment->GetStart() + pSegment->GetLength() - rect.left - 1 + nExpand_Size;
			int nHeight = pSegment->GetHeight() - rect.top + nExpand_Size;
			line(particleImage, Point(nStart, nHeight), Point(nEnd, nHeight), Scalar(nBlackColor), nThickness, nLineType);
		}		
	//--------abstract the contour of the particle.

	
		vector<vector<Point>>contours;
	
		findContours(particleImage, contours, CV_RETR_EXTERNAL, CV_CHAIN_APPROX_NONE);
		if (contours.size()==0)// the particle is too odd that openCV can't find a contour of it. Then we take the upright rect of the particle as it's minArea rect.
		{
			double w = 0, h = 0;
			w = (double)rect.Width()*a_PixelSize;
			h = (double)rect.Height()*a_PixelSize;
			a_pOTSPart->SetDMax(MAX(w, h));
			a_pOTSPart->SetDMin(MIN(w, h));
			a_pOTSPart->SetDMean((w + h) / 2);
			a_pOTSPart->SetFeretDiameter((w + h) / 2);
			a_pOTSPart->SetDElong(MAX(w, h));
			a_pOTSPart->SetPerimeter((w + h) * 2);
			a_pOTSPart->SetDPerp(MIN(w, h));
			a_pOTSPart->SetDInscr(MIN(w, h));
			a_pOTSPart->SetMinWidth(w);
			a_pOTSPart->SetMinHeight(h);
			return true;
		}
		int imaxcontour = 0, imax = 0;
		for (unsigned int i = 0; i < contours.size(); i++) {                  

			int itmp = contourArea(contours[i]);

			if (imaxcontour < itmp) {

				imax = i;

				imaxcontour = itmp;

			}
		}

		vector<Point > listEdge = contours[imax];
		
		vector<vector<Point>>Outcontours;
		Outcontours.push_back(listEdge);

		//---------calculate the minimium  rectangle 

		auto rRect = cv::minAreaRect(listEdge);
		Point2f p[4];
		rRect.points(p);
		int D_MIN =getDistance(p[0], p[1]);
		int D_MinRecLen = 0;//minareaRect's length(the longger side).
		for (int j = 0; j <= 2; j++)
		{
			//line(cvContourImg, p[j], p[(j + 1) % 4], Scalar(100, 100, 0), 2);
			int d = getDistance(p[j], p[j + 1]);
			if (d < D_MIN)
			{
				D_MIN = d;
			}
			if (d > D_MinRecLen)
			{
				D_MinRecLen = d;
			}
		}

		a_pOTSPart->SetDMin(D_MIN*a_PixelSize);

		a_pOTSPart->SetMinWidth(rRect.size.width);
		a_pOTSPart->SetMinHeight(rRect.size.height);
		//----------calculate the perimeter
		double d = arcLength(listEdge, true);
		a_pOTSPart->SetPerimeter(d*a_PixelSize);
		//-----------calculate the Max Diameter. Find the min enclosing circle first ,then find the two farthest circle connected point.
		Point2f center; float radius;
		minEnclosingCircle(listEdge, center, radius);
		//circle(cvContourImg, center, radius, Scalar(100), 2);
		std::vector <Point> outContour = listEdge;
		
		std::vector <Point> rst;
		for (unsigned int k = 0; k < outContour.size(); k++)
		{
			Point p = outContour[k];
			double d = sqrt(pow((p.x - center.x), 2) + pow((p.y - center.y), 2));
			if (fabs(d - radius) < 0.01)
			{
				rst.push_back(p);
			}
		}
		double D_MAX = 0;
		Point lineDmax[2];
		for (unsigned int m = 0; m < rst.size(); m++)
		{
			Point p = rst[m];
			for (unsigned int n = m + 1; n < rst.size(); n++)
			{
				Point p1 = rst[n];
				double d = sqrt(powf((p.x - p1.x), 2) + powf((p.y - p1.y), 2));
				if (d > D_MAX)
				{
					D_MAX = d;
					lineDmax[0] = p;
					lineDmax[1] = p1;
				}
			}
		}
		a_pOTSPart->SetDMax(D_MAX*a_PixelSize);
		
		//--------calculate the D_PERP property using the D_MAX's two endpoints.
		std::vector<Point> curve1;
		std::vector<Point> curve2;
		for (unsigned int i = 0; i < outContour.size(); i++)
		{
			Point pt = outContour[i];
			bool start = false;
			int clockwise = Side(lineDmax[0], lineDmax[1], pt);// devide these points into two group ,separate into the two sides.
			if (clockwise > 0)
			{
				curve1.push_back(pt);
			}
			else
			{
				curve2.push_back(pt);
			}
		}

		double d_perp1 = 0, d_perp2 = 0;
		for (unsigned int i = 0; i < curve1.size(); i++)
		{
			double d = getDist_P2L(curve1[i], lineDmax[0], lineDmax[1]);
			if (d > d_perp1)
			{
				d_perp1 = d;
			}
		}
		for (unsigned int i = 0; i < curve2.size(); i++)
		{
			double d = getDist_P2L(curve2[i], lineDmax[0], lineDmax[1]);
			if (d > d_perp2)
			{
				d_perp2 = d;
			}
		}
		a_pOTSPart->SetDPerp((d_perp1 + d_perp2)*a_PixelSize);
		//----------find the  diameter of max inscribed circle
		int r;
		Point inscribeCirclecenter;
		FindInnerCircleInContour(outContour, inscribeCirclecenter, r);
		//--------------------------------------------------------calculate the xraypos !
		
		CRect rec = a_pOTSPart->GetParticleRect();
	
		a_pOTSPart->SetXRayPos(CPoint(inscribeCirclecenter.x - nExpand_Size + rec.left - 1, inscribeCirclecenter.y - nExpand_Size + rec.top - 1));
	

		a_pOTSPart->SetDInscr(r * 2 * a_PixelSize);
		//---------------calculate the image other caracater: length/width  realArea/minRectangeArea etc. we can use these propertes to do forward process.
		double minRectArea = D_MIN * D_MinRecLen*a_PixelSize*a_PixelSize;//最小外接矩形面积
		double fillRatio = a_pOTSPart->GetActualArea() / minRectArea;//实际面积与最小外接矩形面积比,that's the fill rate.
		double lengthWidthRatio;
		lengthWidthRatio = (double)D_MinRecLen / D_MIN;//长宽比

	
		 //decide if this shape is a strip shape :if the lenthWidthRatio>2 then it is. if the lengthWidthRatio<2 and  the areaRatio<0.5 then it is.
		bool isStripShape = false;
		double curveLength = 0;
		double D_MEAN=0;
		Moments mu;
		mu = moments(listEdge, false);
		int nx = mu.m10 / mu.m00;
		int ny = mu.m01 / mu.m00;
		//circle(cvcopyImg, Point(nx, ny), 1, (255), 1);
		Point ptCenter = Point((int)nx, (int)ny);
		if (pointPolygonTest(listEdge, ptCenter, false) != 1)// the center point doesn't contain in the contour, we think it as curve shape.
		{
			isStripShape = true;
		}
		/*if (lengthWidthRatio >= 2 )// in PartA software this is true,but IncA because of the GB definition the everage feret diameter is always the mean value of all the chord.
		{
			
			isStripShape = true;
		}*/
		
		if (fillRatio <= 0.4)// only when the fill rate is very low,we think it as a curve shape,then we choose the mean width as the feret diameter.
		{
			isStripShape = true;
		}
		
		if (isStripShape)
		{

			curveLength = a_pOTSPart->GetPerimeter()/2 - a_pOTSPart->GetDInscr()/2;// thinking this particle as a strip rectangle.the width is the max inscribe circle diameter/2.
			if (curveLength < D_MAX)
			{
				curveLength = D_MAX;
			}
			if (curveLength < MIN_DOUBLE_VALUE || a_pOTSPart->GetActualArea()<MIN_DOUBLE_VALUE)
			{
				D_MEAN = 0;
			}
			else
			{
				D_MEAN = a_pOTSPart->GetActualArea() / curveLength;
			}
			
			a_pOTSPart->SetDMean(D_MEAN*a_PixelSize);
			a_pOTSPart->SetFeretDiameter(D_MEAN*a_PixelSize);
			a_pOTSPart->SetDElong (curveLength*a_PixelSize);
		}
		else//it's a ball shape particle
		{
			curveLength = D_MAX;
			double ftd = 0, maxD = 0, minD = 0, dratio = 0;
			GetParticleAverageChord(outContour, a_PixelSize, ftd);
			
			a_pOTSPart->SetDMean(ftd);
			a_pOTSPart->SetFeretDiameter(ftd);
			a_pOTSPart->SetDElong(curveLength*a_PixelSize);
		}

		
		if (a_pOTSPart->GetAspectRatio() > 2)
		{
			double angle;
			if (rRect.size.width > rRect.size.height) // w > h
			{
				angle = abs(rRect.angle);
			}
			else
			{
				angle = 90.0 + abs(rRect.angle);
			}
			a_pOTSPart->SetOrientation(angle);
		}
		else
		{
			a_pOTSPart->SetOrientation(0);
		}
		
	return true;
		
	}
	
	
	
	void COTSImageProcess::ImshowImage(CBSEImgPtr img)
	{
		BYTE* data = img->GetImageDataPointer();
		//Mat cvImg;
		cv::Size s;
		s.width = img->GetImageSize().cx;
		s.height = img->GetImageSize().cy;
		Mat cvImg=Mat::zeros(s, CV_8U);
		cvImg.data = data;
	
		cv::imshow("dd", cvImg);

	    cv::waitKey();
	
	}
	void COTSImageProcess::ImshowChartData(CBSEImgPtr img)
	{
		img->SetChartData();
		WORD* data = img->GetBSEChart();
		//Mat cvImg;
		cv::Size s;
		s.width = 255;
		s.height = 100;
		Mat cvImg = Mat::zeros(s, CV_8U);
		//cvImg.data = data;
		WORD nBSEChart[MAXBYTE];
		//1. get chart data
	
		linearSmooth5(data, nBSEChart, MAXBYTE);
		for (int i=1;i<255;i++)
		{
			
			line(cvImg, Point(i, 100-nBSEChart[i]), Point(i+1, 100-nBSEChart[i+1]), Scalar(nBlackColor), 1, 8);
		}

		cv::imshow("chart", cvImg);

		cv::waitKey();
	}
	

	BOOL COTSImageProcess::MergeBigBoundaryParticles(COTSFieldDataList allFields,double pixelSize,int scanFieldSize, CSize ResolutionSize, COTSParticleList& mergedParts)
	{
		class BorderPart
		{
			typedef std::shared_ptr<BorderPart>  CBorderPartPtr;
			BorderPart(COTSParticlePtr p)
			{
				myPart = p;
				headerParticle = NULL;
			}
		public:
			COTSParticlePtr myPart;
			COTSParticle* headerParticle;//used to merge particles ,if this particle has been merged then this pointer will point to the first particle of these merged particles or else it's NULL.
		
			
			static std::vector<CBorderPartPtr> ConvertPartToBorderPart(COTSParticleList parts)
			{
				static std::map<COTSParticle*, CBorderPartPtr> allborderPart;
				std::vector<CBorderPartPtr> borderParts;

				for (auto p : parts)
				{
					if (allborderPart.find(p.get()) == allborderPart.end())
					{
						auto borderp = CBorderPartPtr(new BorderPart(p));
						borderParts.push_back(borderp);
						allborderPart[p.get()] = borderp;
					}
					else
					{
						borderParts.push_back(allborderPart[p.get()]);
					}
					
				}
				return borderParts;
			
			}

		};

		auto FldMgr = new CFieldMgr(scanFieldSize, ResolutionSize);
	

		std::map<COTSParticle*, COTSParticleList> mapMergeParticles;//hold up all the boundary connected particles. the pair's first is also the member of these particles.
		std::map<COTSParticle*, COTSSegmentsList> mapMergedSegments;//hold up all the segment's corresponding clone in the connected particles.

		for (auto centerfld : allFields)
		{
			// find neighbor field on the left.
			auto leftFld = FldMgr->FindNeighborField(allFields, centerfld, SORTING_DIRECTION::LEFT);
			if (leftFld != nullptr)
			{
				auto lParts = centerfld->GetLeftBorderedBigParticles();
				auto rParts = leftFld->GetRightBorderedBigParticles();
				auto leftParts = BorderPart::ConvertPartToBorderPart(lParts);
				auto rightParts = BorderPart::ConvertPartToBorderPart(rParts);

				for (auto leftp : leftParts)
				{
					for (auto rightp : rightParts)
					{
						if (leftp->myPart->IsConnected(rightp->myPart.get(), centerfld->Width, centerfld->Height, (int)SORTING_DIRECTION::LEFT))
						{
							if (leftp->headerParticle != NULL)
							{
								if (rightp->headerParticle == NULL)
								{
									rightp->headerParticle = leftp->headerParticle;
									mapMergeParticles[leftp->headerParticle].push_back(rightp->myPart);
								}
							}
							else
							{
								if (rightp->headerParticle != NULL)
								{
									leftp->headerParticle = rightp->myPart.get();

									mapMergeParticles[rightp->myPart.get()].push_back(leftp->myPart);
								}
								else
								{
									leftp->headerParticle = leftp->myPart.get();
									rightp->headerParticle = leftp->myPart.get();
									mapMergeParticles[leftp->myPart.get()].push_back(rightp->myPart);
								}

							}

						}
					}
				}

			}
			//find neighbor field on the upward
			auto upFld = FldMgr->FindNeighborField(allFields, centerfld, SORTING_DIRECTION::UP);
			if (upFld != nullptr)
			{
				auto topBorderParts = centerfld->GetTopBorderedBigParticles();
				auto bottomBorderParts = upFld->GetBottomBorderedBigParticles();
				auto upParts = BorderPart::ConvertPartToBorderPart(topBorderParts);
				auto downParts = BorderPart::ConvertPartToBorderPart(bottomBorderParts);

				for (auto upprt : upParts)
				{
					for (auto downprt : downParts)
					{
						if (upprt->myPart->IsConnected(downprt->myPart.get(), centerfld->Width, centerfld->Height, (int)SORTING_DIRECTION::UP))
						{

							if (upprt->headerParticle != NULL)
							{
								if (downprt->headerParticle == NULL)
								{
									downprt->headerParticle = upprt->headerParticle;
									mapMergeParticles[upprt->headerParticle].push_back(downprt->myPart);
								}

							}
							else
							{
								if (downprt->headerParticle != NULL)
								{
									upprt->headerParticle = downprt->headerParticle;
									mapMergeParticles[downprt->myPart.get()].push_back(upprt->myPart);
								}
								else
								{
									upprt->headerParticle = upprt->myPart.get();
									downprt->headerParticle = upprt->myPart.get();
									mapMergeParticles[upprt->myPart.get()].push_back(downprt->myPart);
								}

							}

						}
					}
				}

			}
			//find neighbor field on the downward.
			auto downFld = FldMgr->FindNeighborField(allFields, centerfld,SORTING_DIRECTION::DOWN);
			if (downFld != nullptr)
			{
				auto bottomParts = centerfld->GetBottomBorderedBigParticles();
				auto topParts = downFld->GetTopBorderedBigParticles();
				auto downParts = BorderPart::ConvertPartToBorderPart(bottomParts);
				auto upParts= BorderPart::ConvertPartToBorderPart(topParts);

				for (auto downprt : downParts)
				{
					for (auto upprt : upParts)
					{
						if (downprt->myPart->IsConnected(upprt->myPart.get(), centerfld->Width, centerfld->Height, (int)SORTING_DIRECTION::DOWN))
						{

							if (downprt->headerParticle != NULL)
							{
								if (upprt->headerParticle == NULL)
								{
									upprt->headerParticle = downprt->headerParticle;
									mapMergeParticles[downprt->headerParticle].push_back(upprt->myPart);
								}

							}
							else
							{
								if (upprt->headerParticle != NULL)
								{
									downprt->headerParticle = upprt->headerParticle;
									mapMergeParticles[upprt->headerParticle].push_back(downprt->myPart);
								}
								else
								{
									downprt->headerParticle = downprt->myPart.get();
									upprt->headerParticle = downprt->myPart.get();
									mapMergeParticles[downprt->myPart.get()].push_back(upprt->myPart);
								}

							}

						}
					}
				}

			}
			//find neighbor field on the right.
			auto rightFld = FldMgr->FindNeighborField(allFields, centerfld, SORTING_DIRECTION::RIGHT);
			if (rightFld != nullptr)
			{
				auto rParts = centerfld->GetRightBorderedBigParticles();
				auto lParts = rightFld->GetLeftBorderedBigParticles();
				auto rightParts = BorderPart::ConvertPartToBorderPart(rParts);
				auto leftParts = BorderPart::ConvertPartToBorderPart(lParts);

				for (auto rightprt : rightParts)
				{
					for (auto leftprt : leftParts)
					{
						if (rightprt->myPart->IsConnected(leftprt->myPart.get(), centerfld->Width, centerfld->Height, (int)SORTING_DIRECTION::RIGHT))
						{

							if (rightprt->headerParticle != NULL)
							{
								if (leftprt->headerParticle == NULL)
								{
									leftprt->headerParticle = rightprt->headerParticle;
									mapMergeParticles[rightprt->headerParticle].push_back(leftprt->myPart);
								}

							}
							else
							{
								if (leftprt->headerParticle != NULL)
								{
									rightprt->headerParticle = leftprt->headerParticle;
									mapMergeParticles[leftprt->headerParticle].push_back(rightprt->myPart);
								}
								else
								{
									rightprt->headerParticle = rightprt->myPart.get();
									leftprt->headerParticle = rightprt->myPart.get();
									mapMergeParticles[rightprt->myPart.get()].push_back(leftprt->myPart);
								}

							}

						}
					}
				}

			}
		}
		/*for (auto particle : mapMergeParticles)
		{
			
		}*/
		static int partTagId;
		for (auto pair : mapMergeParticles)
		{
			struct EleAreaPercentage
			{
				EleAreaPercentage(double p, CElementChemistryPtr e)
				{
					areaPercentage = p;
					eleData = e;
				}
				double areaPercentage;
				CElementChemistryPtr eleData;
			};

			auto newPart = COTSParticlePtr(new COTSParticle());
			COTSSegmentsList newSegs;
			auto p = pair.first;
			newPart->SetSEMPos(p->GetSEMPos());
			
			//firstly,we sum up all the merged particles's area and get the represent string.
			std::string partsStr = std::to_string(p->GetFieldId()) + ":" + std::to_string(p->GetAnalysisId());
			double allPartArea = p->GetActualArea();//Get the first particle's area.


			for (auto other : pair.second)// Get the total area of all these particles for  the use of ele calcu.
			{
				partsStr += "," + std::to_string(other->GetFieldId()) + ":" + std::to_string(other->GetAnalysisId());//Get the subparticles string such as "1:1,2:1" etc.
				allPartArea += other->GetActualArea();//Get other particle's area

			}
			// calculate all the new segment's position.
			std::vector <COTSParticle*> allSubParts;
			allSubParts.push_back(p);
			for (auto other : pair.second)// Get the total area of all these particles for  the use of ele calcu.
			{
				allSubParts.push_back(other.get());
			}
			for (auto subp : allSubParts)
			{
				int fid = subp->GetFieldId();
				CPoint myFldPos;
				for (auto f : allFields)//find  this particle's filed.
				{
					if (f->GetId() == fid)
					{
						myFldPos = f->GetPosition();
					}
				}
				int fldWidth = allFields[0]->Width;
				int fldHeight = allFields[0]->Height;
				CPoint fldLeftUpPos = CPoint(myFldPos.x + fldWidth / 2 , myFldPos.y + fldHeight / 2 );
				for (auto s : subp->GetFeature()->GetSegmentsList())
				{
					COTSSegmentPtr newseg = COTSSegmentPtr(new COTSSegment());
					newseg->SetStart(s->GetStart()  + fldLeftUpPos.x);
					newseg->SetHeight((0 - s->GetHeight()) + fldLeftUpPos.y);//the coordinate system of segment in a field is different with the OTS coordinate system.OTS system's y axis is upward positive ,yet the field is downward positive.
					newseg->SetLength(s->GetLength());
					newSegs.push_back(newseg);
				}

			}

			COTSFeaturePtr newFeature = COTSFeaturePtr(new COTSFeature());
			newFeature->SetSegmentsList(newSegs);
			newPart->SetFeature(newFeature);
			newPart->CalCoverRectFromSegment();



			//second, we get all the element data and their  area percentage .
			std::map<std::string, std::vector<EleAreaPercentage>> mapEleData;

			CPosXrayPtr pXray1 = p->GetXrayInfo();
			if (pXray1 != nullptr)
			{
				for (auto ele : pXray1->GetElementQuantifyData())
				{
					mapEleData[ele->GetName().GetBuffer()].push_back(EleAreaPercentage(p->GetActualArea() / allPartArea, ele));
				}
			}
			

			for (auto other : pair.second)
			{
				auto otherXray = other->GetXrayInfo();
				if (otherXray != nullptr)
				{
					for (auto eledata : otherXray->GetElementQuantifyData())
					{
						mapEleData[eledata->GetName().GetBuffer()].push_back(EleAreaPercentage(other->GetActualArea() / allPartArea, eledata));

					}
				}
				
			}
			// third,we calculate all the element's new percentage data and get a new element chemistry list.

			CElementChemistriesList newCheList;
			for (auto eledata : mapEleData)
			{
				CElementChemistryPtr newEleche = CElementChemistryPtr(new CElementChemistry());
				newEleche->SetName(CString(eledata.first.c_str()));
				double newPercentage = 0;
				for (auto d : eledata.second)
				{
					newPercentage += d.areaPercentage * d.eleData->GetPercentage();
				}
				newEleche->SetPercentage(newPercentage);
				newCheList.push_back(newEleche);
			}

			CPosXrayPtr xray(new CPosXray());
			xray->SetElementQuantifyData(newCheList);
			newPart->SetXrayInfo(xray);
			newPart->SetConnectedParticlesSequentialString(partsStr);
			newPart->SetActualArea(allPartArea);
			partTagId++;
			newPart->SetParticleId(partTagId);
			newPart->SetAnalysisId(partTagId);
			std::string name = p->GetClassifyName();
			newPart->SetClassifyName(name);
			newPart->SetColor(p->GetColor());
			mergedParts.push_back(newPart);
		}
		return true;
	}
	
	

	
}