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"


namespace OTSIMGPROC
{
	using namespace cv;
	using namespace std;

	const int nBlackColor = 255;

	

	

	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));
		
		RemoveBackGround(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;
	}
	
	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, CIntRangePtr 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::RemoveBackGround(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)
		{
			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
		{
		
			/*auto range = CalBackground(a_pImgIn);
			nBGStart = range->GetStart();

			nBGEnd = range->GetEnd();*/


			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);
		
		//foundedPixelNum = nNumParticle;
		//foundedPixelNum = 100;
		
		
		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::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;
	}
	
	

	
}