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 "../OTSLog/COTSUtilityDllFunExport.h"
#include "FieldMgr.h"
using namespace cv;
using namespace std;

namespace OTSIMGPROC
{   
	namespace 
	{
		
		/***** 求两点间距离*****/
		float getDistance(Point pointO, Point pointA)
		{
			float distance;
			distance = powf((pointO.x - pointA.x), 2) + powf((pointO.y - pointA.y), 2);
			distance = sqrtf(distance);
			return distance;

		}
		/***** 点到直线的距离:P到AB的距离*****/
		//P为线外一点，AB为线段两个端点
		float  getDist_P2L(Point pointP, Point pointA, Point pointB)
		{
			//求直线方程
			int A = 0, B = 0, C = 0;
			A = pointA.y - pointB.y;
			B = pointB.x - pointA.x;
			C = pointA.x*pointB.y - pointA.y*pointB.x;
			//代入点到直线距离公式
			float distance = 0;
			distance = ((float)abs(A*pointP.x + B * pointP.y + C)) / ((float)sqrtf(A*A + B * B));
			return distance;
		}
		int  Side(Point P1, Point P2, Point point)
		{
			/*Point P1 = line.P1;
			Point P2 = line.P2;*/
			return ((P2.y - P1.y) * point.x + (P1.x - P2.x) * point.y + (P2.x*P1.y - P1.x*P2.y));
		}
		void  FindInnerCircleInContour(vector<Point> contour, Point &center, int &radius)
		{
			Rect r = boundingRect(contour);
			int nL = r.x, nR = r.br().x; //轮廓左右边界
			int nT = r.y, nB = r.br().y; //轮廓上下边界

			double dist = 0;
			double maxdist = 0;

			for (int i = nL; i < nR; i++)  //列
			{
				for (int j = nT; j < nB; j++)  //行
				{
					//计算轮廓内部各点到最近轮廓点的距离
					dist = pointPolygonTest(contour, Point(i, j), true);
					if (dist > maxdist)
					{
						//求最大距离，只有轮廓最中心的点才距离最大
						maxdist = dist;
						center = Point(i, j);
					}
				}
			}
			radius = maxdist;  //圆半径
		}
		BOOL GetParticleAverageChord(std::vector<Point>  listEdge, double a_PixelSize, double &dPartFTD)
		{
			// safety check

			double nx = 0, ny = 0;
			Moments mu;
			mu = moments(listEdge, false);
			nx = mu.m10 / mu.m00;
			ny = mu.m01 / mu.m00;
			//circle(cvcopyImg, Point(nx, ny), 1, (255), 1);
			Point ptCenter = Point((int)nx, (int)ny);
			// coordinate transformation
			Point ptPosition;
			int radiusNum = 0;
			// get ferret diameter
			double sumFltDiameter = 0;
			int interval;
			int edgePointNum = listEdge.size();
			if (edgePointNum > 10)
			{
				interval = edgePointNum / 10;//get one line per 10 degree  aproxemately 
			}
			else
			{
				interval = 1;
			}
			for (int i = 0; i < edgePointNum; i++)
			{
				Point pt = listEdge[i];
				ptPosition.x = abs(pt.x - ptCenter.x);
				ptPosition.y = abs(pt.y - ptCenter.y);


				if (i % interval == 0)//calculate one line per 10 point ,so to speed up.don't calculate all the diameter.
				{
					double r1 = sqrt(pow(ptPosition.x, 2) + pow(ptPosition.y, 2));
					sumFltDiameter += r1;
					radiusNum += 1;
					//line(cvImageData, ptCenter, pt, Scalar(nBlackColor), nThickness, nLineType);
				}
			}
			if (radiusNum == 0)
			{
				dPartFTD = 0;
			}
			else
			{
				dPartFTD = a_PixelSize * sumFltDiameter / radiusNum * 2;

			}
			//imshow("feret center", cvImageData);
			return TRUE;
		}
		void linearSmooth5(WORD wordIn[], WORD wordOut[], int N = 255)//smooth algorithm
		{
			double in[256];
			double out[256];
			double smoothCurveData[256];
			for (int i = 0; i < 256; i++)
			{
				in[i] = (double)wordIn[i];
			}

			int i;
			if (N < 5)
			{
				for (i = 0; i <= N - 1; i++)
				{
					out[i] = in[i];
				}
			}
			else
			{
				out[0] = (3.0 * in[0] + 2.0 * in[1] + in[2] - in[4]) / 5.0;
				out[1] = (4.0 * in[0] + 3.0 * in[1] + 2 * in[2] + in[3]) / 10.0;
				for (i = 2; i <= N - 3; i++)
				{
					out[i] = (in[i - 2] + in[i - 1] + in[i] + in[i + 1] + in[i + 2]) / 5.0;
				}
				out[N - 2] = (4.0 * in[N - 1] + 3.0 * in[N - 2] + 2 * in[N - 3] + in[N - 4]) / 10.0;
				out[N - 1] = (3.0 * in[N - 1] + 2.0 * in[N - 2] + in[N - 3] - in[N - 5]) / 5.0;
			}
			for (int i = 0; i < N; i++)
			{
				wordOut[i] = (WORD)out[i];
			}
		}
		void GetMatricsParticlesFromRawParticle(COTSParticlePtr a_pOTSPart,int imageWidth,int imageHeight, double a_PixelSize, int xrayStep, COTSParticleList& matricsParts)
		{
		
			auto originalSegs = a_pOTSPart->GetFeature()->GetSegmentsList();
			std::map<int, COTSSegmentsList> segsOnTheSameHeight;
			for (auto s : originalSegs)
			{
				segsOnTheSameHeight[s->GetHeight()].push_back(s);
			}
			auto rect = a_pOTSPart->GetParticleRect();
			

			std::vector<CPoint> matrixPs;

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

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

				}
			}
			a_pOTSPart->SetXrayMatrixPoints(matrixPs);
			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();
					for (int i = 0; i < originalSegs.size();i++)//judge if the seg is in the original particle scope.
					{
						auto rseg = originalSegs[i];
						
						int rsegStart = rseg->GetStart();
						int rsegEnd = rseg->GetEnd();

						if (segStart > rsegEnd || rsegStart > segEnd)//there's no intersection.is not a valid seg for this original segment.
						{
							
							continue;
						}

						if (segStart>=rsegStart && segEnd <= rsegEnd)//contained in the original segment,is a valid seg.
						{
							segs.push_back(seg);
							continue;
						}						
						else if (segStart>= rsegStart&& segEnd >= rsegEnd)// intersect in the head end.Modify the end of the seg .
						{
							COTSSegmentPtr newseg = COTSSegmentPtr(new COTSSegment());
							newseg->SetStart(segStart);
							newseg->SetEnd(rsegEnd);
							newseg->SetHeight(currentH);
							segs.push_back(newseg);
	
							continue;
						}
						else if (segStart<= rsegStart&& segEnd >= rsegEnd)
						{
							COTSSegmentPtr newseg = COTSSegmentPtr(new COTSSegment());
							newseg->SetStart(rsegStart);
							newseg->SetEnd(rsegEnd);
							newseg->SetHeight(currentH);
					
							segs.push_back(newseg);
							continue;

						}
						else if (segStart<= rsegStart&& rsegEnd >= segEnd)
						{
							COTSSegmentPtr newseg = COTSSegmentPtr(new COTSSegment());
							newseg->SetStart(rsegStart);
							newseg->SetEnd(segEnd);
							newseg->SetHeight(currentH);

							segs.push_back(newseg);
							continue;
						}
						
					}

					
				}
				if (segs.size() > 0)
				{
					part->GetFeature()->SetSegmentsList(segs);
					part->CalXRayPos();
					part->SetFieldId(a_pOTSPart->GetFieldId());
					part->SetAnalysisId(a_pOTSPart->GetAnalysisId());
					matricsParts.push_back(part);
				}
				
			}


		}
	}

	COTSImageProcess::COTSImageProcess()
	{
	}


	COTSImageProcess::~COTSImageProcess()
	{
	}

	// use verticl line of 3 pixel to erode a image
	void COTSImageProcess::BErodeVertical3(LPBYTE source, LPBYTE target, WORD rows, WORD columns)
	{
		WORD x, y, wcounts;

		if (rows <= 2 || columns <= 2)return;

		// top line
		for (x = 0; x < columns; x++)
		{
			*(target + x) = 0;
		}

		// bottom line
		for (x = 0; x < columns; x++)
		{
			*(target + (DWORD)(rows - 1)*columns + x) = 0;
		}

		for (y = 1; y<rows - 1; y++)
		{
			for (x = 0; x<columns; x++)
			{
				if (*(source + (DWORD)y*columns + x) == 0)
				{
					*(target + (DWORD)y*columns + x) = 0;
					continue;
				}

				wcounts = 0;
				if (*(source + (DWORD)(y - 1)*columns + x) == 255)
				{
					wcounts++;
				}
				if (*(source + (DWORD)(y + 1)*columns + x) == 255)
				{
					wcounts++;
				}

				if (wcounts == 2) *(target + (DWORD)y*columns + x) = 255;
				else *(target + (DWORD)y*columns + x) = 0;
			}
		}

	}
	// use left 45 degree line of 3 pixel to erode a image
	void COTSImageProcess::BErodeLeft45Degree3(LPBYTE source, LPBYTE target, WORD rows, WORD columns)
	{
		WORD x, y, wcounts;

		if (rows <= 2 || columns <= 2)return;

		// top line
		for (x = 0; x < columns; x++)
		{
			*(target + x) = 0;
		}

		// bottom line
		for (x = 0; x < columns; x++)
		{
			*(target + (DWORD)(rows - 1)*columns + x) = 0;
		}

		// left line
		for (y = 0; y<rows; y++)
		{
			*(target + (DWORD)y*columns) = 0;
		}

		// right line
		for (y = 0; y<rows; y++)
		{
			*(target + (DWORD)y*columns + columns - 1) = 0;
		}


		for (y = 1; y < rows - 1; y++)
		{
			for (x = 1; x < columns - 1; x++)
			{
				if (*(source + (DWORD)y*columns + x) == 0)
				{
					*(target + (DWORD)y*columns + x) = 0;
					continue;
				}

				wcounts = 0;
				if (*(source + (DWORD)(y - 1)*columns + x - 1) == 255)
				{
					wcounts++;
				}
				if (*(source + (DWORD)(y + 1)*columns + x + 1) == 255)
				{
					wcounts++;
				}

				if (wcounts == 2) *(target + (DWORD)y*columns + x) = 255;
				else *(target + (DWORD)y*columns + x) = 0;
			}

		}
		

	}
	// use horizoontal line of 3 pixel to erode a image
	void COTSImageProcess::BErodeHorizontal3(LPBYTE source, LPBYTE target, WORD rows, WORD columns)
	{
		WORD x, y, wcounts;

		if (rows <= 2 || columns <= 2)return;

		// left line
		for (y = 0; y<rows; y++)
		{
			*(target + (DWORD)y*columns) = 0;
		}

		// right line
		for (y = 0; y<rows; y++)
		{
			*(target + (DWORD)y*columns + columns - 1) = 0;
		}


		for (y = 0; y<rows; y++)
		{
			for (x = 1; x<columns - 1; x++)
			{
				if (*(source + (DWORD)y*columns + x) == 0)
				{
					*(target + (DWORD)y*columns + x) = 0;
					continue;
				}

				wcounts = 0;
				if (*(source + (DWORD)y*columns + x - 1) == 255)
				{
					wcounts++;
				}
				if (*(source + (DWORD)y*columns + x + 1) == 255)
				{
					wcounts++;
				}

				if (wcounts == 2) *(target + (DWORD)y*columns + x) = 255;
				else *(target + (DWORD)y*columns + x) = 0;
			}
		}

	}
	// use right 45 degree line of 3 pixel to erode a image
	void COTSImageProcess::BErodeRight45Degree3(LPBYTE source, LPBYTE target, WORD rows, WORD columns)
	{
		WORD x, y, wcounts;

		if (rows <= 2 || columns <= 2)return;

		// top line
		for (x = 0; x < columns; x++)
		{
			*(target + x) = 0;
		}

		// bottom line
		for (x = 0; x < columns; x++)
		{
			*(target + (DWORD)(rows - 1)*columns + x) = 0;
		}

		// left line
		for (y = 0; y<rows; y++)
		{
			*(target + (DWORD)y*columns) = 0;
		}

		// right line
		for (y = 0; y<rows; y++)
		{
			*(target + (DWORD)y*columns + columns - 1) = 0;
		}


		for (y = 1; y<rows - 1; y++)
		{
			for (x = 1; x<columns - 1; x++)
			{
				if (*(source + (DWORD)y*columns + x) == 0)
				{
					*(target + (DWORD)y*columns + x) = 0;
					continue;
				}

				wcounts = 0;
				if (*(source + (DWORD)(y - 1)*columns + x + 1) == 255)
				{
					wcounts++;
				}
				if (*(source + (DWORD)(y + 1)*columns + x - 1) == 255)
				{
					wcounts++;
				}

				if (wcounts == 2) *(target + (DWORD)y*columns + x) = 255;
				else *(target + (DWORD)y*columns + x) = 0;
			}
		}

	
	}
	void COTSImageProcess::BDilateVertical3(LPBYTE source, LPBYTE target, WORD rows, WORD columns)
	{
		WORD x, y;

		if (rows <= 2 || columns <= 2)return;

		// top line
		for (x = 0; x<columns; x++)
		{
			if (*(source + x) != 0)
			{
				*(target + x) = 0xff;
			}
		}

		// bottom line
		for (x = 0; x<columns; x++)
		{
			if (*(source + (DWORD)(rows - 1)*columns + x) != 0)
			{
				*(target + (DWORD)(rows - 1)*columns + x) = 0xff;
			}
		}

		for (y = 1; y<rows - 1; y++)
		{
			for (x = 1; x<columns - 1; x++)
			{
				if (*(source + (DWORD)y*columns + x) != 0)
				{
					*(target + (DWORD)y*columns + x) = 0xff;
					*(target + (DWORD)(y - 1)*columns + x) = 255;
					*(target + (DWORD)(y + 1)*columns + x) = 255;
				}
				else *(target + (DWORD)y*columns + x) = 0;
			}
		}

	}
	void COTSImageProcess::BDilateLeft45Degree3(LPBYTE source, LPBYTE target, WORD rows, WORD columns)
	{
		WORD x, y;

		if (rows <= 2 || columns <= 2)return;

		// top line
		for (x = 0; x<columns; x++)
		{
			if (*(source + x) != 0)
			{
				*(target + x) = 0xff;
			}
		}

		// bottom line
		for (x = 0; x<columns; x++)
		{
			if (*(source + (DWORD)(rows - 1)*columns + x) != 0)
			{
				*(target + (DWORD)(rows - 1)*columns + x) = 0xff;
			}
		}

		// left line
		for (y = 0; y<rows; y++)
		{
			if (*(source + (DWORD)y*columns) != 0)
			{
				*(target + (DWORD)y*columns) = 0xff;
			}
		}

		// right line
		for (y = 0; y<rows; y++)
		{
			if (*(source + (DWORD)y*columns + columns - 1) != 0)
			{
				*(target + (DWORD)y*columns + columns - 1) = 0xff;
			}
		}

		for (y = 1; y<rows - 1; y++)
		{
			for (x = 1; x<columns - 1; x++)
			{
				if (*(source + (DWORD)y*columns + x) != 0)
				{
					*(target + (DWORD)y*columns + x) = 0xff;
					*(target + (DWORD)(y - 1)*columns + x - 1) = 255;
					*(target + (DWORD)(y + 1)*columns + x + 1) = 255;
				}
				else *(target + (DWORD)y*columns + x) = 0;
			}
		}

	}
	void COTSImageProcess::BDilateHorizontal3(LPBYTE source, LPBYTE target, WORD rows, WORD columns)
	{
		WORD x, y;

		if (rows <= 2 || columns <= 2)return;

		// left line
		for (y = 0; y<rows; y++)
		{
			if (*(source + (DWORD)y*columns) != 0)
			{
				*(target + (DWORD)y*columns) = 0xff;
			}
		}

		// right line
		for (y = 0; y<rows; y++)
		{
			if (*(source + (DWORD)y*columns + columns - 1) != 0)
			{
				*(target + (DWORD)y*columns + columns - 1) = 0xff;
			}
		}

		for (y = 1; y<rows - 1; y++)
		{
			for (x = 1; x<columns - 1; x++)
			{
				if (*(source + (DWORD)y*columns + x) != 0)
				{
					*(target + (DWORD)y*columns + x) = 0xff;
					*(target + (DWORD)y*columns + x - 1) = 255;
					*(target + (DWORD)y*columns + x + 1) = 255;
				}
				else *(target + (DWORD)y*columns + x) = 0;
			}
		}
	}
	void COTSImageProcess::BDilateRight45Degree3(LPBYTE source, LPBYTE target, WORD rows, WORD columns)
	{
		WORD x, y;

		if (rows <= 2 || columns <= 2)return;

		// top line
		for (x = 0; x<columns; x++)
		{
			if (*(source + x) != 0)
			{
				*(target + x) = 0xff;
			}
		}

		// bottom line
		for (x = 0; x<columns; x++)
		{
			if (*(source + (DWORD)(rows - 1)*columns + x) != 0)
			{
				*(target + (DWORD)(rows - 1)*columns + x) = 0xff;
			}
		}

		// left line
		for (y = 0; y<rows; y++)
		{
			if (*(source + (DWORD)y*columns) != 0)
			{
				*(target + (DWORD)y*columns) = 0xff;
			}
		}

		// right line
		for (y = 0; y<rows; y++)
		{
			if (*(source + (DWORD)y*columns + columns - 1) != 0)
			{
				*(target + (DWORD)y*columns + columns - 1) = 0xff;
			}
		}

		for (y = 1; y<rows - 1; y++)
		{
			for (x = 1; x<columns - 1; x++)
			{
				if (*(source + (DWORD)y*columns + x) != 0)
				{
					*(target + (DWORD)y*columns + x) = 0xff;
					*(target + (DWORD)(y - 1)*columns + x + 1) = 255;
					*(target + (DWORD)(y + 1)*columns + x - 1) = 255;
				}
				else *(target + (DWORD)y*columns + x) = 0;
			}
		}

	}
	void COTSImageProcess::BErode3(LPBYTE source, LPBYTE target, WORD wDegree, WORD rows, WORD columns)
	{
		WORD x, y, i, j, wcounts;

		if (rows == 1 || columns == 1)return;

		for (y = 1; y<rows - 1; y++)
		{
			for (x = 1; x<columns - 1; x++)
			{
				if (*(source + (DWORD)y*columns + x) == 0)
				{
					*(target + (DWORD)y*columns + x) = 0;
					continue;
				}
				wcounts = 0;
				for (i = (WORD)(y - 1); i <= (WORD)(y + 1); i++)

				{
					for (j = (WORD)(x - 1); j <= (WORD)(x + 1); j++)
					{
						if (*(source + (DWORD)i*columns + j) == 0)
						{
							wcounts++;
						}
					}
				}
				if (wcounts >= wDegree) *(target + (DWORD)y*columns + x) = 0;
				else *(target + (DWORD)y*columns + x) = 0xff;
			}
		}

		// top line
		for (x = 1; x<columns - (WORD)1; x++)
		{
			if (*(source + x) == 0)
			{
				*(target + x) = 0;
				continue;
			}
			wcounts = 0;
			for (i = 0; i <= 1; i++)

			{
				for (j = (WORD)(x - 1); j <= (WORD)(x + 1); j++)
				{
					if (*(source + (DWORD)i*columns + j) == 0) wcounts++;
				}
			}
			if (wcounts >= wDegree * 5 / 8) *(target + x) = 0;
			else *(target + x) = 0xff;
		}

		// bottom line
		for (x = 1; x<columns - 1; x++)
		{
			if (*(source + (DWORD)(rows - 1)*columns + x) == 0)
			{
				*(target + (DWORD)(rows - 1)*columns + x) = 0;
				continue;
			}
			wcounts = 0;
			for (i = (WORD)(rows - 2); i <= (WORD)(rows - 1); i++)

			{
				for (j = (WORD)(x - 1); j <= (WORD)(x + 1); j++)
				{
					if (*(source + (DWORD)i*columns + j) == 0) wcounts++;
				}
			}
			if (wcounts >= wDegree * 5 / 8) *(target + (DWORD)(rows - 1)*columns + x) = 0;
			else *(target + (DWORD)(rows - 1)*columns + x) = 0xff;
		}

		// left line
		for (y = 1; y<rows - 1; y++)
		{
			if (*(source + (DWORD)y*columns) == 0)
			{
				*(target + (DWORD)y*columns) = 0;
				continue;
			}
			wcounts = 0;
			for (i = (WORD)(y - 1); i <= (WORD)(y + 1); i++)
			{
				for (j = 0; j <= 1; j++)
				{
					if (*(source + (DWORD)i*columns + j) == 0) wcounts++;
				}
			}
			if (wcounts >= wDegree * 5 / 8) *(target + (DWORD)y*columns) = 0;
			else *(target + (DWORD)y*columns) = 0xff;
		}

		// right line
		for (y = 1; y<rows - 1; y++)
		{
			if (*(source + (DWORD)y*columns + columns - 1) == 0)
			{
				*(target + (DWORD)y*columns + columns - 1) = 0;
				continue;
			}
			wcounts = 0;
			for (i = (WORD)(y - 1); i <= (WORD)(y + 1); i++)

			{
				for (j = (WORD)(columns - 2); j <= (WORD)(columns - 1); j++)
				{
					if (*(source + (DWORD)i*columns + j) == 0) wcounts++;
				}
			}
			if (wcounts >= wDegree * 5 / 8) *(target + (DWORD)y*columns + columns - 1) = 0;
			else *(target + (DWORD)y*columns + columns - 1) = 0xff;
		}

		return;
	}
	void COTSImageProcess::BDilate3(LPBYTE source, LPBYTE target, WORD wDegree, WORD rows, WORD columns)
	{
		WORD x, y, i, j, wcounts;

		for (y = 1; y<rows - 1; y++)
		{
			for (x = 1; x<columns - 1; x++)
			{
				if (*(source + (DWORD)y*columns + x) != 0)
				{
					*(target + (DWORD)y*columns + x) = 0xff;
					continue;
				}
				wcounts = 0;
				for (i = (WORD)(y - 1); i <= (WORD)(y + 1); i++)
				{
					for (j = (WORD)(x - 1); j <= (WORD)(x + 1); j++)
					{
						if (*(source + (DWORD)i*columns + j) != 0)	wcounts++;
					}
				}
				if (wcounts >= wDegree) *(target + (DWORD)y*columns + x) = 0xff;
				else *(target + (DWORD)y*columns + x) = 0;
			}
		}

		// top line
		for (x = 1; x<columns - 1; x++)
		{
			if (*(source + x) != 0)
			{
				*(target + x) = 0xff;
				continue;
			}
			wcounts = 0;
			for (i = 0; i <= 1; i++)
			{
				for (j = (WORD)(x - 1); j <= (WORD)(x + 1); j++)
				{
					if (*(source + (DWORD)i*columns + j) != 0) wcounts++;
				}
			}

			if (wcounts >= wDegree * 5 / 8)	// but does not mater, as we have border of 2 now
			{
				*(target + x) = 0xff;
			}
			else { *(target + x) = 0; }
		}

		// bottom line
		for (x = 1; x<columns - 1; x++)
		{
			if (*(source + (DWORD)(rows - 1)*columns + x) != 0)
			{
				*(target + (DWORD)(rows - 1)*columns + x) = 0xff;
				continue;
			}
			wcounts = 0;
			for (i = (WORD)(rows - 2); i <= (WORD)(rows - 1); i++)
			{
				for (j = (WORD)(x - 1); j <= (WORD)(x + 1); j++)
				{
					if (*(source + (DWORD)i*columns + j) != 0) wcounts++;
				}
			}

			if (wcounts > wDegree * 5 / 8)
			{
				*(target + (DWORD)(rows - 1)*columns + x) = 0xff;
			}
			else
			{
				*(target + (DWORD)(rows - 1)*columns + x) = 0;
			}
		}

		// left line
		for (y = 1; y<rows - 1; y++)
		{
			if (*(source + (DWORD)y*columns) != 0)
			{
				*(target + (DWORD)y*columns) = 0xff;
				continue;
			}
			wcounts = 0;
			for (i = (WORD)(y - 1); i <= (WORD)(y + 1); i++)
			{
				for (j = 0; j <= (WORD)1; j++)
				{
					if (*(source + (DWORD)i*columns + j) != 0) wcounts++;
				}
			}
			if (wcounts >= wDegree * 5 / 8)
			{
				*(target + (DWORD)y*columns) = 0xff;
			}
			else
			{
				*(target + (DWORD)y*columns) = 0;
			}
		}

		// right line
		for (y = 1; y<rows - 1; y++)
		{
			if (*(source + (DWORD)y*columns + columns - 1) != 0)
			{
				*(target + (DWORD)y*columns + columns - 1) = 0xff;
				continue;
			}
			wcounts = 0;
			for (i = (WORD)(y - 1); i <= (WORD)(y + 1); i++)
			{
				for (j = (WORD)(columns - 2); j <= (WORD)(columns - 1); j++)
				{
					if (*(source + (DWORD)i*columns + j) != 0) wcounts++;
				}
			}
			if (wcounts >= wDegree * 5 / 8)
			{
				*(target + (DWORD)y*columns + columns - 1) = 0xff;
			}
			else
			{
				*(target + (DWORD)y*columns + columns - 1) = 0;
			}
		}

		// four cornor points treated separately here
		// top-left
		if (*(source) != 0)
		{
			*target = 0xff;
		}
		else
		{
			wcounts = 0;
			if (*(source + 1) != 0) wcounts++;
			if (*(source + columns) != 0) wcounts++;
			if (*(source + columns + 1) != 0) wcounts++;

			//        if (wcounts >= wDegree*3/8)  // this is a bug here - interger division
			if (wcounts * 8 >= wDegree * 3)
			{
				*target = 0xff;
			}
			else
			{
				*target = 0;
			}
		}

		//top-right
		if (*(source + columns - 1) != 0)
		{
			*(target + columns - 1) = 0xff;
		}
		else
		{
			wcounts = 0;
			if (*(source + columns - 2) != 0) wcounts++;
			if (*(source + columns * 2 - 1) != 0) wcounts++;
			if (*(source + columns * 2 - 2) != 0) wcounts++;

			//        if (wcounts >= wDegree*3/8)  // this is a bug here - interger division
			if (wcounts * 8 >= wDegree * 3)
			{
				*(target + columns - 1) = 0xff;
			}
			else
			{
				*(target + columns - 1) = 0;
			}
		}

		//bottom-left
		if (*(source + (DWORD)columns * (rows - 1)) != 0)
		{
			*(target + (DWORD)columns * (rows - 1)) = 0xff;
		}
		else
		{
			wcounts = 0;
			if (*(source + (DWORD)columns * (rows - 1) + 1) != 0) wcounts++;
			if (*(source + (DWORD)columns * (rows - 2)) != 0) wcounts++;
			if (*(source + (DWORD)columns * (rows - 2) + 1) != 0) wcounts++;

			//        if (wcounts >= wDegree*3/8)  // this is a bug here - interger division
			if (wcounts * 8 >= wDegree * 3)
			{
				*(target + (DWORD)columns * (rows - 1)) = 0xff;
			}
			else
			{
				*(target + (DWORD)columns * (rows - 1)) = 0;
			}
		}

		//bottom-right
		if (*(source + (DWORD)columns * rows - 1) != 0)
		{
			*(target + (DWORD)columns * rows - 1) = 0xff;
		}
		else
		{
			wcounts = 0;
			if (*(source + (DWORD)columns * rows - 2) != 0) wcounts++;
			if (*(source + (DWORD)columns * (rows - 1) - 2) != 0) wcounts++;
			if (*(source + (DWORD)columns * (rows - 1) - 1) != 0) wcounts++;

			//        if (wcounts >= wDegree*3/8)  // this is a bug here - interger division
			if (wcounts * 8 >= wDegree * 3)
			{
				*(target + (DWORD)columns * rows - 1) = 0xff;
			}
			else
			{
				*(target + (DWORD)columns * rows - 1) = 0;
			}
		}

		return;
	}

	// ReZoom the picture with re-magnification 
	BOOL COTSImageProcess::ReZoom(CString InPutPath, CString OutPutPath)
	{
		Mat cvSrcImg;
		string strInputPath;
		strInputPath = CStringA(InPutPath);

		// Pictures loop in folder
		std::vector<cv::String> ImageFolder;
		cv::glob(strInputPath, ImageFolder);
		if (ImageFolder.size() == 0)
		{
			return FALSE;
		}

		for (unsigned int nImgNum = 0; nImgNum < ImageFolder.size(); ++nImgNum) {
			cvSrcImg = cv::imread(ImageFolder[nImgNum], CV_LOAD_IMAGE_GRAYSCALE);

			// Image convolution operation	
			//// convolution kernel
			float kernel[] = { -1,  -1 , -1,  -1 , 0, -1,  -1 ,  -1 , -1 };
			cv::Mat ker = cv::Mat(nImage_Size, nImage_Size, CV_32F, &kernel);
			cv::Mat cvDstImg = cv::Mat(cvSrcImg.size(), cvSrcImg.type());
			// anchor of the kernel
			cv::Point anchor(-1, -1);
			cv::filter2D(cvSrcImg, cvDstImg, CV_32F, ker, anchor, delta, cv::THRESH_TRUNC);

			// Maximum Pixel Value
			cvDstImg = abs(cvDstImg);
			double minVal, maxVal;
			minMaxLoc(cvDstImg, &minVal, &maxVal);

			// Grayscale image
			int nReduce;
			Mat onesImg = Mat::ones(cvDstImg.rows, cvDstImg.cols, CV_32F) * (int)minVal;
			absdiff(cvDstImg, onesImg, cvDstImg);
			nReduce = (int)maxVal - minVal;
			cvDstImg = cvDstImg * nBlackColor / nReduce;

			// Output image convert data to int
			cvDstImg.convertTo(cvDstImg, CV_8U);

			// Process the picture to 128 pixels
			resize(cvDstImg, cvDstImg, Size(nPictureSize, nPictureSize));
			threshold(cvDstImg, cvDstImg, nProcessParam, nBlackColor, CV_THRESH_BINARY);

			string strOutPutPath;
			strOutPutPath = CStringA(OutPutPath);

			imwrite(strOutPutPath , cvDstImg);
		}
		return TRUE;
	}
	BOOL COTSImageProcess::RemoveBSEImageBG(CBSEImgPtr m_pBSEImg, COTSImageProcessParamPtr a_pImgProcessParam,COTSFieldDataPtr m_pFieldData)
	{
		ASSERT(m_pFieldData);
		if (!m_pFieldData)
		{
			LogErrorTrace(__FILE__, __LINE__, _T("RemoveBSEImageBG: there is no field data"));
			return FALSE;
		}

		ASSERT(m_pBSEImg);
		if (!m_pBSEImg)
		{
			LogErrorTrace(__FILE__, __LINE__, _T("RemoveBSEImageBG: there is no image data"));
			return FALSE;
		}

		ASSERT(a_pImgProcessParam);
		if (!a_pImgProcessParam)
		{
			LogErrorTrace(__FILE__, __LINE__, _T("RemoveBSEImageBG: there is no image process data"));
			return FALSE;
		}

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

		
		

		BYTE* pSrcImg = m_pBSEImg->GetImageDataPointer();

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

		CBSEImgPtr imgNoBGBinary = CBSEImgPtr(new CBSEImg(r));
		long nNumParticle = 0;
		RemoveBackGround(m_pBSEImg, a_pImgProcessParam, imgNoBGBinary,nNumParticle);

	

		
		BYTE* pPixel = imgNoBGBinary->GetImageDataPointer();
		
	
		

		long nPtStart = a_pImgProcessParam->GetParticleGray().GetStart();
		long nPtEnd = a_pImgProcessParam->GetParticleGray().GetEnd();


		if (nNumParticle == 0)
		{

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

		}
		else
		{

			

			// get the area image
			

			Mat blurImg;
			Mat srcImg = Mat(nHeightImg, nWidthImg, CV_8UC1, pPixel);
			
			//medianBlur(srcImg, blurImg, 3);//smooth the edge


			COTSParticleList listParticleOut;
			if (!GetParticles(0,0,nWidthImg, nHeightImg, srcImg.data, listParticleOut))
			{
				COTSParticleList listParticleEmpty;
				listParticleEmpty.clear();
				m_pFieldData->SetParticleList(listParticleEmpty);
				

			}
			


			// form a image only have particles on
			COTSSegmentsList listImage;
			for (auto pParticle : listParticleOut)
			{
				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)
				{
					// update image list
					COTSSegmentPtr pSegNew = COTSSegmentPtr(new COTSSegment(*pSegment.get()));
					listImage.push_back(pSegNew);


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


					if (nHeightS > nHeightImg)
					{
						LogErrorTrace(__FILE__, __LINE__, _T("seg height is wrong."));
						return FALSE;
					}
					if ((nStartS + nLengthS - 1) > nWidthImg)
					{
						LogErrorTrace(__FILE__, __LINE__, _T("seg starst and length is wrong."));
						return FALSE;
					}

					for (unsigned int i = 0; i < nLengthS; i++)
					{
						if ((nStartS + i) > nWidthImg)
						{
							LogErrorTrace(__FILE__, __LINE__, _T("seg start is wrong."));
							return FALSE;
						}
						else if (nHeightS > nHeightImg)
						{
							LogErrorTrace(__FILE__, __LINE__, _T("seg height is wrong."));
							return FALSE;
						}

						long nValueTemp = (long)*(pSrcImg + nHeightS * nWidthImg + nStartS + i);
						nPixelAll += nValueTemp;
					}
				}
				BYTE nAveGray = (BYTE)(nPixelAll / nPixelNum);
				pParticle->SetAveGray(nAveGray);
				pParticle->SetActualArea(nPixelNum);

			}
			m_pFieldData->SetParticleList(listParticleOut);
		}

		
		delete[]pTempImg;


		return TRUE;
	
	
	}
	BOOL COTSImageProcess::RemoveBGByFindContour(CBSEImgPtr m_pBSEImg, COTSImageProcessParamPtr a_pImageProcessParam, COTSFieldDataPtr m_pFieldData)
	{

		ASSERT(m_pFieldData);

		ASSERT(m_pBSEImg);	

		ASSERT(a_pImageProcessParam);
	

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


		BYTE* pSrcImg = m_pBSEImg->GetImageDataPointer();

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

		CBSEImgPtr imgNoBGBinary = CBSEImgPtr(new CBSEImg(r));
		long nNumParticle = 0;
		RemoveBackGround(m_pBSEImg, a_pImageProcessParam, imgNoBGBinary, nNumParticle);




		BYTE* pPixel = imgNoBGBinary->GetImageDataPointer();




		long nPtStart = a_pImageProcessParam->GetParticleGray().GetStart();
		long nPtEnd = a_pImageProcessParam->GetParticleGray().GetEnd();


		if (nNumParticle == 0)
		{

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

		}
		else
		{

			// get the area image	

			Mat cvcopyImg = Mat(nHeightImg, nWidthImg, CV_8UC1, pPixel);
		

			vector<vector<Point>>contours;
			

			findContours(cvcopyImg, contours, CV_RETR_EXTERNAL, CV_CHAIN_APPROX_NONE);

			COTSParticleList listParticleOut;
			for (size_t i = 0; i < contours.size(); i++)
			{
				
				Rect rectMax = boundingRect(contours[i]);
				Mat  rectROI = cvcopyImg(rectMax).clone();


				//exclude the point which  intersect into this bounding box but is not in this contour.
				for (int nX = 0; nX < rectROI.rows; nX++)
				{
					for (int nY = 0; nY < rectROI.cols; nY++)
					{
						double localPos = pointPolygonTest(contours[i], Point2f(nX + rectMax.x, nY + rectMax.y), false);
						if (localPos == -1)
						{
							rectROI.data[nX, nY] = 0;//set the value to 0,so we won't consider it when we find segment and feature in this ROI.
						}

					}

				}

				GetParticles(rectMax.x, rectMax.y, rectMax.width, rectMax.height, rectROI.data, listParticleOut);
				
			}



			// form a image only have particles on
			//COTSSegmentsList listImage;
			for (auto pParticle : listParticleOut)
			{
				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);
				pParticle->SetActualArea(nPixelNum);

			}
			m_pFieldData->SetParticleList(listParticleOut);
		}


		delete[]pTempImg;


		return TRUE;
	}
	BOOL COTSImageProcess::RemoveBGByCVconnectivities(CBSEImgPtr inBSEImg, COTSImageProcessParamPtr a_pImageProcessParam,double a_pixelSize,  COTSFieldDataPtr m_pFieldData)
	{

		ASSERT(m_pFieldData);

		ASSERT(inBSEImg);

		ASSERT(a_pImageProcessParam);


		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));
		long nNumParticle = 0;
		RemoveBackGround(inBSEImg, a_pImageProcessParam, imgNoBGBinary, nNumParticle);




		BYTE* pPixel = imgNoBGBinary->GetImageDataPointer();




		long nPtStart = a_pImageProcessParam->GetParticleGray().GetStart();
		long nPtEnd = a_pImageProcessParam->GetParticleGray().GetEnd();


		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 blurImg;
			//medianBlur(cvcopyImg, blurImg, 5);//get rid of the noise point.
			Mat labels = Mat::zeros(cvcopyImg.size(), CV_32S);
			Mat  stats, centroids;
			int number = connectedComponentsWithStats(cvcopyImg, labels, stats, centroids, 8, CV_32S);
			
			double rMin = a_pImageProcessParam->GetIncArea().GetStart()/2.0;
			double rMax = a_pImageProcessParam->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;
			

			for (auto pParticle : listParticleOut)
			{
				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);
			


			}
			m_pFieldData->SetParticleList(listParticleOut);
		}


		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);


			}

			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);
		linearSmooth5(secondSmooth, secondSmooth, MAXBYTE);
		linearSmooth5(secondSmooth, secondSmooth, MAXBYTE);
		linearSmooth5(secondSmooth, 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;
					
				}
		
			}
		

			//Mat cvcopyImg = Mat(nHeightImg, nWidthImg, CV_8UC1, pPixel);// use the medianblur method to achieve the same effect as open morphology(errod and dialate).
			//pPixel = cvcopyImg.data;
			
				COTSImageProcess::BErode3(pPixel, pTempImg, 5, nHeightImg, nWidthImg);
				COTSImageProcess::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* 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 (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();

			// delete background 
			for (unsigned int i = 0; i < nImgSize; i++)
			{
				if (pSrcImg[i] >= nBGStart && pSrcImg[i] <= nBGEnd)
				{
					pPixel[i] = 0;
				}
				else
				{
					pPixel[i] = 255;
					nNumParticle++;
				}
				if (pSrcImg[i]<nPartStart || pSrcImg[i]>nPartEnd)
				{
					pPixel[i] = 0;
				}
			}
			int errodDilateParam = a_pImageProcessParam->GetErrodDilateParam();
			if (errodDilateParam > 0)
			{
				COTSImageProcess::BErode3(pPixel, pTempImg, errodDilateParam, nHeightImg, nWidthImg);
				COTSImageProcess::BDilate3(pTempImg, pPixel, errodDilateParam, nHeightImg, nWidthImg);
			}
			
			//Mat cvcopyImg = Mat(nHeightImg, nWidthImg, CV_8UC1, pPixel);// use the medianblur method to achieve the same effect as open morphology(errod and dialate).
			//medianBlur(cvcopyImg, cvcopyImg, 5);
			//pPixel = cvcopyImg.data;

			
		}
		else
		{
		
			auto range = CalBackground(a_pImgIn);
			nBGStart = range->GetStart();

			nBGEnd = range->GetEnd();


			switch (a_pImageProcessParam->GetAutoBGRemoveType())
			{
			case OTS_AUTOBGREMOVE_TYPE::DOWNWARD:
				for (unsigned int i = 0; i < nImgSize; i++)
				{
					if (pSrcImg[i] <= nBGEnd)
					{
						pPixel[i] = 0;
					}
					else
					{
						pPixel[i] = 255;
						nNumParticle++;
					}
				}
				break;
			case OTS_AUTOBGREMOVE_TYPE::UPWARD:
				for (unsigned int i = 0; i < nImgSize; i++)
				{
					if (pSrcImg[i] >= nBGStart)
					{
						pPixel[i] = 0;
					}
					else
					{
						pPixel[i] = 255;
						nNumParticle++;
					}
				}
				break;
			case OTS_AUTOBGREMOVE_TYPE::MIDDLE:
				for (unsigned int i = 0; i < nImgSize; i++)
				{
					if (pSrcImg[i] >= nBGStart && pSrcImg[i] <= nBGEnd)
					{
						pPixel[i] = 0;
					}
					else
					{
						pPixel[i] = 255;
						nNumParticle++;
					}
				}
				break;
			default:
				break;
			}
			int errodDilateParam = a_pImageProcessParam->GetErrodDilateParam();
			if (errodDilateParam > 0)
			{
				COTSImageProcess::BErode3(pPixel, pTempImg, errodDilateParam, nHeightImg, nWidthImg);
				COTSImageProcess::BDilate3(pTempImg, pPixel, errodDilateParam, nHeightImg, nWidthImg);
			}
			/*Mat cvcopyImg = Mat(nHeightImg, nWidthImg, CV_8UC1, pPixel);
			medianBlur(cvcopyImg, cvcopyImg, 5);
			pPixel = cvcopyImg.data;*/
		}
		a_pBinImgOut->SetImageData(pPixel,nWidthImg,nHeightImg);
		
		foundedPixelNum = nNumParticle;
		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::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));
			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.
		Mat cvcopyImg;
		medianBlur(particleImage, cvcopyImg, 5);//smooth the edge
	
		vector<vector<Point>>contours;
	
		findContours(cvcopyImg, 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));
			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 lenger 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->SetOrientation(rRect.angle);
		//----------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);
		}
	return true;
		
	}
	BOOL COTSImageProcess::SplitRawParticleIntoMatricsParticle(COTSParticlePtr a_pOTSPart,int imageWidth,int imageHeight, double a_PixelSize, double a_XrayStep)
	{
		//--------- 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));
			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.
		Mat cvcopyImg;
		medianBlur(particleImage, cvcopyImg, 5);//smooth the edge
		vector<vector<Point>>contours;

		findContours(cvcopyImg, 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));
			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];
		if (a_XrayStep > 0)
		{
			COTSParticleList matricsParts;
			int xrayStep = a_XrayStep;// *a_PixelSize;
			GetMatricsParticlesFromRawParticle(a_pOTSPart, imageWidth,imageHeight,a_PixelSize, xrayStep, matricsParts);
			a_pOTSPart->SetSubParticles(matricsParts);
		}
		//-----------
	}
	BOOL COTSImageProcess::SplitRawParticleIntoGreyScaleParticle(COTSParticlePtr a_pOTSPart, double a_PixelSize ,CBSEImgPtr fieldImg)
	{
		//--------- 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() < 5 * 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));
			return true;
		}
		
		// calculate the particle image data size, expand 3 pixel at the edge
	
		CBSEImgPtr onePartImg = CBSEImgPtr(new CBSEImg(CRect(0,0, fieldImg->GetWidth(), fieldImg->GetHeight())));
		// get the segment list
		COTSSegmentsList listSegment = a_pOTSPart->GetFeature()->GetSegmentsList();
		for (auto pSegment : listSegment)
		{
			for (int i = 0; i < pSegment->GetLength(); i++)
			{
				int x = pSegment->GetStart() + i;
				int y = pSegment->GetHeight();
				int bseValue = fieldImg->GetBSEValue(x,y);
				onePartImg->SetBSEValue(x, y,bseValue);

			}
		}
		/*ImshowImage(onePartImg);
		ImshowChartData(onePartImg);*/
		
		std::vector<CIntRangePtr> rngs = CalcuGrayLevelRange(onePartImg);
		
		

		COTSFieldDataPtr partData = COTSFieldDataPtr(new COTSFieldData());
		std::map<int, std::vector<COTSParticlePtr>> partAreaMap;
		for (int i = 0; i < rngs.size(); i++)
		{
			partAreaMap.clear();
			GetParticlesBySpecialGrayRange(onePartImg, rngs[i], CDoubleRangePtr(new CDoubleRange(0, 2000)), a_PixelSize, partData);

		
			
			for (auto p : partData->GetParticleList())//sorting and filtering
			{
				/*if (p->GetActualArea() > 50)
				{*/
					partAreaMap[p->GetPixelArea()].push_back(p);
				//}
				
			}
			if(partAreaMap.size()>0)
			{
				auto theBiggestPart = partAreaMap.rbegin()->second[0];
					theBiggestPart->CalXRayPos();
					std::map<int, std::vector<COTSParticlePtr>>::reverse_iterator it;
					auto partsegs = theBiggestPart->GetFeature()->GetSegmentsList();
					it = partAreaMap.rbegin()++;
					for (; it != partAreaMap.rend(); it++)
					{
						for (auto sameAreaP : it->second)
						{
							auto segs = sameAreaP->GetFeature()->GetSegmentsList();
							for (auto s : segs)
							{
								partsegs.push_back(s);
							}
						}

					}
				theBiggestPart->GetFeature()->SetSegmentsList(partsegs, true);
				theBiggestPart->SetFieldId(a_pOTSPart->GetFieldId());
				theBiggestPart->SetAnalysisId(a_pOTSPart->GetAnalysisId());
				a_pOTSPart->AddSubParticle(theBiggestPart);
			
			
			}
			

			/*for (auto p : partData->GetParticleList())
			{
				a_pOTSPart->AddSubParticle(p); 
			}*/
		}
		



		return 0;
	}
	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)
	{
		COTSSegmentsList boarderSegs;

		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 leftParts = centerfld->GetLeftBorderedBigParticles();
				auto rightParts = leftFld->GetRightBorderedBigParticles();

				for (auto leftp : leftParts)
				{
					for (auto rightp : rightParts)
					{
						if (leftp->IsConnected(rightp.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);
								}
							}
							else
							{
								if (rightp->headerParticle != NULL)
								{
									leftp->headerParticle = rightp.get();

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

							}

						}
					}
				}

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

				for (auto upprt : upParts)
				{
					for (auto downprt : downParts)
					{
						if (upprt->IsConnected(downprt.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);
								}

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

							}

						}
					}
				}

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

				for (auto downprt : downParts)
				{
					for (auto upprt : upParts)
					{
						if (downprt->IsConnected(upprt.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);
								}

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

							}

						}
					}
				}

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

				for (auto rightprt : rightParts)
				{
					for (auto leftprt : leftParts)
					{
						if (rightprt->IsConnected(leftprt.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);
								}

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

							}

						}
					}
				}

			}
		}
		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->SetAbsolutePos(p->GetAbsolutPos());
			
			//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->CalCoverRect();



			//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->SetSubParticlesName(partsStr);
			newPart->SetActualArea(allPartArea);
			partTagId++;
			newPart->SetParticleId(partTagId);
			newPart->SetAnalysisId(partTagId);
			std::string name = p->TypeName();
			newPart->TypeName(name);
			newPart->TypeColor(p->TypeColor());
			mergedParts.push_back(newPart);
		}
		return true;
	}
	
		
	
}