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"
#include <unordered_map>


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;
		GetMatrixPointsFromRect(rect, xrayStep, matrixPs);
			
		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);
		std::vector<CPoint> matrixseeds;
		matrixseeds.clear();
		if (seeds.size() <50)
		{
			auto step = areaRng.GetStart() * 4;
			GetMatrixPointsFromRect(fieldImg->GetImageRect(), step / a_PixelSize, matrixseeds);
			seeds.insert(seeds.end(), matrixseeds.begin(), matrixseeds.end());
		}
		
		

		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 newsubpart = 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);
					}

				}

			}
			newsubpart->SetFeature(newFea);


			newsubpart->CalXRayPos();
			newsubpart->CalActualArea(a_PixelSize);

			CPoint pos = newsubpart->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 = newsubpart->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);
					newsubpart->SetFieldId(mp->GetFieldId());
					newsubpart->SetOTSRect(recInOTSCord);
					mp->AddSubParticle(newsubpart);
					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->InitChartData();
	
		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);
	
		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;
		}
		
		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->InitChartData();

		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 chart 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::GetMatrixPointsFromRect(CRect rect, double xrayStep, 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);

				if (rect.PtInRect(thePoint))
				{
					matrixPs.push_back(thePoint);	
				}
				


			}
		}
		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 = 4;
		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() < 20 * 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/2;
			int nEnd = nStart + pSegment->GetLength() ; 
			int nHeight = pSegment->GetHeight() - rect.top+ nExpand_Size/2;
			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->InitChartData();
		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)
	{
		// Lightweight wrapper used only in this function (no global/static map).
		struct BorderPart
		{
			explicit BorderPart(const COTSParticlePtr& p) : myPart(p) {}
			COTSParticlePtr LockMyPart() const { return myPart; }
			COTSParticlePtr LockHeaderParticle() const { return headerParticle.lock(); }

			COTSParticlePtr myPart;
			std::weak_ptr<COTSParticle> headerParticle; // points to representative particle if merged
		};

		// Map lives only for this invocation -> no leaked persistent state between calls.
		std::unordered_map<COTSParticle*, std::shared_ptr<BorderPart>> allborderPart;
		allborderPart.reserve(1024);

			// Helper: convert COTSParticleList -> vector<BorderPart shared_ptr>
	auto ConvertPartToBorderPart = [&](const COTSParticleList& parts) -> std::vector<std::shared_ptr<BorderPart>>
	{
		std::vector<std::shared_ptr<BorderPart>> borderParts;
		borderParts.reserve(parts.size());
		for (auto p : parts)
		{
			COTSParticle* raw = p.get();
			auto it = allborderPart.find(raw);
			if (it == allborderPart.end())
			{
				auto bp = std::make_shared<BorderPart>(p);
				borderParts.push_back(bp);
				allborderPart[raw] = bp;
			}
			else
			{
				borderParts.push_back(it->second);
			}
		}
		return borderParts;
	};

	// Helper: get header particle pointer from raw COTSParticle*
	auto GetHeaderParticlePtr = [&](COTSParticle* raw) -> COTSParticlePtr
	{
		if (!raw) return nullptr;
		auto it = allborderPart.find(raw);
		if (it == allborderPart.end()) return nullptr;
		return it->second->LockHeaderParticle();
	};

	// Use RAII to avoid leaks
	std::unique_ptr<CFieldMgr> FldMgr(new CFieldMgr(scanFieldSize, ResolutionSize));

	std::map<COTSParticle*, COTSParticleList> mapMergeParticles; // hold boundary-connected particles
	std::map<COTSParticle*, COTSSegmentsList> mapMergedSegments;  // kept for compatibility with existing logic (unused below)

	// Iterate fields and check 4 neighbor directions
	for (auto centerfld : allFields)
	{
		// LEFT neighbor
		if (auto leftFld = FldMgr->FindNeighborField(allFields, centerfld, SORTING_DIRECTION::LEFT))
		{
			auto lParts = centerfld->GetLeftBorderedBigParticles();
			auto rParts = leftFld->GetRightBorderedBigParticles();
			auto leftParts = ConvertPartToBorderPart(lParts);
			auto rightParts = ConvertPartToBorderPart(rParts);

			for (auto& leftp : leftParts)
			{
				for (auto& rightp : rightParts)
				{
					auto left_sp = leftp->LockMyPart();
					auto right_sp = rightp->LockMyPart();
					if (!left_sp || !right_sp) continue;

					if (left_sp->IsConnected(right_sp.get(), centerfld->Width, centerfld->Height, (int)SORTING_DIRECTION::LEFT))
					{
						auto left_hdr = leftp->LockHeaderParticle();
						auto right_hdr = rightp->LockHeaderParticle();
						if (left_hdr)
						{
							if (!right_hdr)
							{
								rightp->headerParticle = leftp->headerParticle;
								mapMergeParticles[left_hdr.get()].push_back(right_sp);
							}
						}
						else
						{
							if (right_hdr)
							{
								leftp->headerParticle = right_hdr;
								mapMergeParticles[right_hdr.get()].push_back(left_sp);
							}
							else
							{
								leftp->headerParticle = left_sp;
								rightp->headerParticle = left_sp;
								mapMergeParticles[left_sp.get()].push_back(right_sp);
							}
						}
					}
				}
			}
		}

		// UP neighbor
		if (auto upFld = FldMgr->FindNeighborField(allFields, centerfld, SORTING_DIRECTION::UP))
		{
			auto topBorderParts = centerfld->GetTopBorderedBigParticles();
			auto bottomBorderParts = upFld->GetBottomBorderedBigParticles();
			auto upParts = ConvertPartToBorderPart(topBorderParts);
			auto downParts = ConvertPartToBorderPart(bottomBorderParts);

			for (auto& upprt : upParts)
			{
				for (auto& downprt : downParts)
				{
					auto up_sp = upprt->LockMyPart();
					auto down_sp = downprt->LockMyPart();
					if (!up_sp || !down_sp) continue;

					if (up_sp->IsConnected(down_sp.get(), centerfld->Width, centerfld->Height, (int)SORTING_DIRECTION::UP))
					{
						auto up_hdr = upprt->LockHeaderParticle();
						auto down_hdr = downprt->LockHeaderParticle();
						if (up_hdr)
						{
							if (!down_hdr)
							{
								downprt->headerParticle = up_hdr;
								mapMergeParticles[up_hdr.get()].push_back(down_sp);
							}
						}
						else
						{
							if (down_hdr)
							{
								upprt->headerParticle = down_hdr;
								mapMergeParticles[down_hdr.get()].push_back(up_sp);
							}
							else
							{
								upprt->headerParticle = up_sp;
								downprt->headerParticle = up_sp;
								mapMergeParticles[up_sp.get()].push_back(down_sp);
							}
						}
					}
				}
			}
		}

		// DOWN neighbor
		if (auto downFld = FldMgr->FindNeighborField(allFields, centerfld, SORTING_DIRECTION::DOWN))
		{
			auto bottomParts = centerfld->GetBottomBorderedBigParticles();
			auto topParts = downFld->GetTopBorderedBigParticles();
			auto downParts = ConvertPartToBorderPart(bottomParts);
			auto upParts = ConvertPartToBorderPart(topParts);

			for (auto& downprt : downParts)
			{
				for (auto& upprt : upParts)
				{
					auto down_sp = downprt->LockMyPart();
					auto up_sp = upprt->LockMyPart();
					if (!down_sp || !up_sp) continue;

					if (down_sp->IsConnected(up_sp.get(), centerfld->Width, centerfld->Height, (int)SORTING_DIRECTION::DOWN))
					{
						auto down_hdr = downprt->LockHeaderParticle();
						auto up_hdr = upprt->LockHeaderParticle();
						if (down_hdr)
						{
							if (!up_hdr)
							{
								upprt->headerParticle = down_hdr;
								mapMergeParticles[down_hdr.get()].push_back(up_sp);
							}
						}
						else
						{
							if (up_hdr)
							{
								downprt->headerParticle = up_hdr;
								mapMergeParticles[up_hdr.get()].push_back(down_sp);
							}
							else
							{
								downprt->headerParticle = down_sp;
								upprt->headerParticle = down_sp;
								mapMergeParticles[down_sp.get()].push_back(up_sp);
							}
						}
					}
				}
			}
		}

		// RIGHT neighbor
		if (auto rightFld = FldMgr->FindNeighborField(allFields, centerfld, SORTING_DIRECTION::RIGHT))
		{
			auto rParts = centerfld->GetRightBorderedBigParticles();
			auto lParts = rightFld->GetLeftBorderedBigParticles();
			auto rightParts = ConvertPartToBorderPart(rParts);
			auto leftParts = ConvertPartToBorderPart(lParts);

			for (auto& rightprt : rightParts)
			{
				for (auto& leftprt : leftParts)
				{
					auto right_sp = rightprt->LockMyPart();
					auto left_sp = leftprt->LockMyPart();
					if (!right_sp || !left_sp) continue;

					if (right_sp->IsConnected(left_sp.get(), centerfld->Width, centerfld->Height, (int)SORTING_DIRECTION::RIGHT))
					{
						auto right_hdr = rightprt->LockHeaderParticle();
						auto left_hdr = leftprt->LockHeaderParticle();
						if (right_hdr)
						{
							if (!left_hdr)
							{
								leftprt->headerParticle = right_hdr;
								mapMergeParticles[right_hdr.get()].push_back(left_sp);
							}
						}
						else
						{
							if (left_hdr)
							{
								rightprt->headerParticle = left_hdr;
								mapMergeParticles[left_hdr.get()].push_back(right_sp);
							}
							else
							{
								rightprt->headerParticle = right_sp;
								leftprt->headerParticle = right_sp;
								mapMergeParticles[right_sp.get()].push_back(left_sp);
							}
						}
					}
				}
			}
		}
	}

	// Merge groups into new particles
	static int partTagId=1000; // preserved behavior (was static previously). If you want thread-safety use std::atomic<int>.
	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->SetFieldSEMPos(p->GetFieldSEMPos());
		
		// build parts string and total area
		std::string partsStr = std::to_string(p->GetFieldId()) + ":" + std::to_string(p->GetAnalysisId());
		double allPartArea = p->GetActualArea();

		for (auto other : pair.second)
		{
			partsStr += "," + std::to_string(other->GetFieldId()) + ":" + std::to_string(other->GetAnalysisId());
			allPartArea += other->GetActualArea();
		}

		// collect sub-parts (use header particle for the primary entry)
		COTSParticleList allSubParts;
		allSubParts.push_back(GetHeaderParticlePtr(p));
		for (auto other : pair.second)
		{
			allSubParts.push_back(other);
		}

		// compute merged segments in global OTS coords
		for (auto subp : allSubParts)
		{
			int fid = subp->GetFieldId();
			CPoint myFldPos;
			for (auto f : allFields)
			{
				if (f->GetId() == fid)
				{
					myFldPos = f->GetPosition();
					break;
				}
			}
			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);
				newseg->SetLength(s->GetLength());
				newSegs.push_back(newseg);
			}
		}

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

		// merge element chemistry weighted by area
		std::map<std::string, std::vector<EleAreaPercentage>> mapEleData;
		CPosXrayPtr pXray1 = p->GetXrayInfo();
		if (pXray1)
		{
			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)
			{
				for (auto eledata : otherXray->GetElementQuantifyData())
				{
					mapEleData[eledata->GetName().GetBuffer()].push_back(EleAreaPercentage(other->GetActualArea() / allPartArea, eledata));
				}
			}
		}

		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);
		partTagId++;
		xray->SetIndex(partTagId);
		newPart->SetXrayInfo(xray);
		newPart->SetConnectedParticlesSequentialString(partsStr);
		newPart->SetSubParticles(allSubParts);
		newPart->SetActualArea(allPartArea);
		newPart->SetOTSScanFieldPos(p->GetOTSScanFieldPos());	
		newPart->SetFieldSEMPos(p->GetFieldSEMPos());
		newPart->SetParticleId(partTagId);
		newPart->SetAnalysisId(partTagId);
		std::string name = p->GetTypeName();
		newPart->SetTypeName(name);
		newPart->SetFieldId(p->GetFieldId());
		newPart->SetColor(p->GetColor());
		mergedParts.push_back(newPart);
	}

	return true;
	}

	BOOL COTSImageProcess::CalcuBrightPhaseDarkPhaseGrayByOtsu(CBSEImgPtr a_pImgIn, int& bright, int& dark)
	{
		Mat srcImgMat = GetMatDataFromBseImg(a_pImgIn);
		Mat brightPhaseMat, darkPhaseMat;			
		double thresholdValue = 0;
		thresholdValue =cv::threshold(srcImgMat, brightPhaseMat, thresholdValue, 255, THRESH_TOZERO|cv::THRESH_OTSU);
		cv::threshold(srcImgMat, darkPhaseMat, thresholdValue, 255, THRESH_TOZERO_INV);
		/*cv::Scalar mean1, std1;
		cv::Scalar mean2, std2;
		
		cv::meanStdDev(brightPhaseMat, mean1, std1);
		cv::meanStdDev(darkPhaseMat, mean2, std2);
		bright = (int)mean1[0];
		dark = (int)mean2[0];*/	
		int imagesize = a_pImgIn->GetWidth() * a_pImgIn->GetHeight();

		BYTE* briData = brightPhaseMat.data;
		BYTE* darkData = srcImgMat.data;

		double sum = 0;
		int num = 0;
		for (int i = 0;  i < imagesize;i++)
		{
			if (briData[i] > 0)
			{
				num += 1;
				sum += briData[i];
			}
		}
		bright =(int)( sum / num);
		sum = 0;
		num = 0;
		for (int i = 0;  i < imagesize;i++)
		{
			if (darkData[i] < thresholdValue)
			{
				num += 1;
				sum += darkData[i];
			}
		}
		dark =(int)( sum / num);
		
		return TRUE;
	}
	
	

	
}