OTS/OTSCPP/OTSImagePro/GBImgPropCal.cpp

#pragma once

#include "math.h"
#include "stdafx.h"
#include "GBImgPropCal.h"
#include <opencv2/core/core.hpp>  
#include <opencv2/highgui/highgui.hpp> 
#include <opencv2/opencv.hpp>
#include "OTSImageProcess.h"
// calculate image properties
namespace OTSIMGPROC {

	using namespace OTSDATA;
	using namespace cv;
	
	using namespace std;
	//make matrix filled with 255 
	const int nBlackColor = 255;
	// construct
	CGBImgPropCal::CGBImgPropCal()
	{
		Init();
	}

	// destruct
	CGBImgPropCal::~CGBImgPropCal()
	{
		Clear();
	}
	
	// get particle ferret diameter
	BOOL CGBImgPropCal::GetParticleFTD(COTSParticlePtr a_pOTSPart, double a_PixelSize , double &dPartFTD, double &dMinWidth, double &dMaxLength, double &dRatio)
	{
		// safety check
		ASSERT(a_pOTSPart);
		

		// 1. convert particle data to image data
		// get the segment list
		COTSSegmentsList listSegment = a_pOTSPart->GetFeature()->GetSegmentsList();
		// get rectangle of the particle
		CRect rect = a_pOTSPart->GetParticleRect();

		// calculate image data size, expand 1 pixel at the edge
		int nWidth = rect.Width() + 1 + nExpand_Size*2;
		int nHeight = rect.Height() + 1 + nExpand_Size*2;

		// define the image data
		Mat cvImageData = Mat::ones(nHeight, nWidth, CV_8U) * nWhiteColor;
		// go through segment list
		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(cvImageData, Point(nStart, nHeight), Point(nEnd, nHeight), Scalar(nBlackColor), nThickness, nLineType);
		}
		

		// 2. draw the contours
		Mat cvcopyImg = cvImageData.clone();
		//cvImageData.release();
		// contours
		threshold(cvcopyImg, cvcopyImg, nWhiteColor, 100, CV_THRESH_BINARY_INV);
		Canny(cvcopyImg, cvcopyImg, 20, 20 * 2, 3);

		//  find contours	
		vector<vector<Point>> contours;
		findContours(cvcopyImg, contours, CV_RETR_EXTERNAL, CV_CHAIN_APPROX_NONE);
		//Mat cvcopyImg1 = cvcopyImg.clone();
		//drawContours(cvcopyImg1, contours, -1, Scalar(nWhiteColor), 1, 8);
		//cvcopyImg = cvcopyImg - cvcopyImg1;
		//findContours(cvcopyImg, contours, CV_RETR_EXTERNAL, CV_CHAIN_APPROX_NONE);
		
		vector<Point> listEdge;
		for (auto listcont : contours)
		{
			for (auto pt : listcont)
			{
				auto itr = std::find_if(listEdge.begin(), listEdge.end(), [pt](Point p) 
				{ 
					return (p.x == pt.x)&&(p.y == pt.y);
				});
				if (itr == listEdge.end())
				{
					listEdge.push_back(pt);
				}
			}
		}

		double nx = 0, ny = 0;
		for (auto pt : listEdge)
		{
			nx = nx + pt.x*pt.x;
			ny = ny + pt.y*pt.y;
		}

		nx = sqrt(nx / (double)listEdge.size());
		ny = sqrt(ny / (double)listEdge.size());


		// 3. calculate the center point
		// find the center point
		/*Moments edgeMoments = moments(listEdge, true);
		Point2f CenterPoint = Point2f(float(edgeMoments.m10 / edgeMoments.m00), float(edgeMoments.m01 / edgeMoments.m00));
		Point ptCenter = (Point)CenterPoint;
*/
		Point ptCenter = Point((int)nx, (int)ny);
		// 4. sort contours into 4 quadrant

		vector<CPoint> QuaOnelist;
		vector<CPoint> QuaTwolist;
		vector<CPoint> QuaThreelist;
		vector<CPoint> QuaFourlist;
		vector<CPoint> Xaxislist;
		vector<CPoint> XaxisUplist;
		vector<CPoint> XaxisDownlist;
		vector<CPoint> Yaxislist;
		vector<CPoint> YaxisRightlist;
		vector<CPoint> YaxisLeftlist;
		
		int Xmax = 0;
		int Ymax = 0;

		// coordinate transformation
		Point ptPosition;

		for (auto pt : listEdge)
		{
			// for touch the contour
			if (pt.x > Xmax)
				Xmax = pt.x;
			if (pt.y > Ymax)
				Ymax = pt.y;

			ptPosition.x = pt.x - ptCenter.x;
			ptPosition.y = ptCenter.y - pt.y;
			
			if (ptPosition.x >0 && ptPosition.y>0)
			{
				QuaOnelist.push_back(CPoint(ptPosition.x,ptPosition.y));
			}
			else if (ptPosition.x <0 && ptPosition.y>0)
			{
				QuaTwolist.push_back(CPoint(ptPosition.x, ptPosition.y));
			}
			else if (ptPosition.x <0 && ptPosition.y<0)
			{
				QuaThreelist.push_back(CPoint(ptPosition.x, ptPosition.y));
			}
			else if (ptPosition.x >0 && ptPosition.y<0)
			{
				QuaFourlist.push_back(CPoint(ptPosition.x, ptPosition.y));
			}
			else if (pt.x == ptCenter.x)
			{
				Yaxislist.push_back(CPoint(ptPosition.x, ptPosition.y));
			}
			else if (pt.y == ptCenter.y)
			{
				Xaxislist.push_back(CPoint(ptPosition.x, ptPosition.y));
			}

			if (ptPosition.y == 1)
			{
				XaxisUplist.push_back(CPoint(ptPosition.x, ptPosition.y));
			}
			else if (ptPosition.y == -1)
			{
				XaxisDownlist.push_back(CPoint(ptPosition.x, ptPosition.y));
			}

			if (ptPosition.x == 1)
			{
				YaxisRightlist.push_back(CPoint(ptPosition.x, ptPosition.y));
			}
			else if (ptPosition.x == -1)
			{
				YaxisLeftlist.push_back(CPoint(ptPosition.x, ptPosition.y));
			}
		}

		// get ferret diameter
		vector<double> FltDiameter;

		if ((Yaxislist.size() == 2) && Xmax >= ptCenter.x)
		{
			double d = (double)abs(Yaxislist[0].y - Yaxislist[1].y);
			FltDiameter.push_back(d);
		}
		else if (Yaxislist.size() > 2)
		{
			double d = GetDisFromRLNeigbor(Yaxislist, YaxisRightlist, YaxisLeftlist);
			FltDiameter.push_back(d);
		}		

		if ((Xaxislist.size() == 2) && Ymax >= ptCenter.y)
		{
			double d = (double)abs(Xaxislist[0].x - Xaxislist[1].x);
			FltDiameter.push_back(d);
		}
		else if (Xaxislist.size() > 2)
		{
			double d = GetDisFromUDNeigbor(Xaxislist, XaxisUplist, XaxisDownlist);
			FltDiameter.push_back(d);
		}		

		if ((int)QuaOnelist.size() <= ANGLE_MAX || (int)QuaThreelist.size() <= ANGLE_MAX)
		{
			if (!GetDisFrom2list(QuaOnelist, QuaThreelist, FltDiameter))
			{
				return FALSE;
			}
		}
		else 
		{
			if (!GetDisFrom13list(QuaOnelist, QuaThreelist, FltDiameter))
			{
				return FALSE;
			}
		}
		
		if ((int)QuaTwolist.size() <= ANGLE_MAX || (int)QuaFourlist.size() <= ANGLE_MAX)
		{			
			if (!GetDisFrom2list(QuaTwolist, QuaFourlist, FltDiameter))
			{
				return FALSE;
			}
		}
		else
		{
			if (!GetDisFrom24list(QuaTwolist,QuaFourlist, FltDiameter))
			{
				return FALSE;
			}
		}

		double dSum = 0;
		double dMax = 0;
		double dMin = 0;

		double a_dAveDia = 0;
		GetSuperParticleSize(a_pOTSPart, dMin, dMax, a_dAveDia);

		for (auto d : FltDiameter)
		{
			dSum = dSum + d;
			if (d > dMax)
				dMax = d;

			if (d < dMin && d >0)
				dMin = d;
		}
		
		//calculate the distance
		dPartFTD = a_PixelSize * dSum / (int)FltDiameter.size();

		if ((dSum == 0) || (int)FltDiameter.size() == 0)
		{
			dPartFTD = a_PixelSize * a_dAveDia;
		}

		dMaxLength = a_PixelSize * dMax;
		dMinWidth = a_PixelSize * dMin;
		if (dMin != 0)
		{ 
			dRatio = dMax / dMin;
		}
		//cvcopyImg.release();
		return TRUE;
	}
	
	// initialization
	void CGBImgPropCal::Init()
	{
	}

	// clear
	void CGBImgPropCal::Clear()
	{
	}

	double CGBImgPropCal::GetDisFrom3(double a_s1, double a_s2, double a_s3)
	{
		double d = 0;

		if (a_s1*a_s2 > 0)
		{
			d = (double)abs(a_s2 - a_s1);
		}
		else if (a_s2*a_s3 > 0)
		{
			d = (double)abs(a_s2 - a_s3);
		}
		else if (a_s2*a_s3 > 0)
		{
			d = (double)abs(a_s3 - a_s1);
		}

		return d;
	}

	BOOL CGBImgPropCal::GetDisFrom2list(std::vector<CPoint> a_list1, std::vector<CPoint> a_list2, std::vector<double>& a_listFlt)
	{
	
		if ((int)a_list1.size() < (int)a_list2.size())
		{
			if (a_list1.empty())
			{
				for (auto pt : a_list2)
				{
					double d = sqrt(pt.x*pt.x + pt.y *pt.y);
					a_listFlt.push_back(d);
				}
			}
			else
			{
				for (auto pt : a_list1)
				{
					double r1 = sqrt(pt.x*pt.x + pt.y *pt.y);

					double k = (double)pt.y / (double)pt.x;
					double kmin = 100;

					for (auto pot : a_list2)
					{
						double kThree = (double)pot.y / (double)pot.x;
						if (abs(kThree - k) < kmin)
						{
							kmin = abs(kThree - k);
						}
					}

					double r2;
					for (auto pot : a_list2)
					{
						double kThree = (double)pot.y / (double)pot.x;
						if (abs(kThree - k) == kmin)
						{
							r2 = sqrt(pot.x * pot.x + pot.y *pot.y);
							break;
						}
					}
					a_listFlt.push_back((r1 + r2));
				}
			}
		}
		else
		{
			if (a_list2.empty())
			{
				for (auto pt : a_list1)
				{
					double d = sqrt(pt.x*pt.x + pt.y *pt.y);
					a_listFlt.push_back(d);
				}
			}
			else
			{
				for (auto pt : a_list2)
				{
					double r1 = sqrt(pt.x*pt.x + pt.y *pt.y);

					double k = (double)pt.y / (double)pt.x;
					double kmin = 100;

					for (auto pot : a_list1)
					{
						double kThree = (double)pot.y / (double)pot.x;
						if (abs(kThree - k) < kmin)
						{
							kmin = abs(kThree - k);
						}
					}

					double r2;
					for (auto pot : a_list1)
					{
						double kThree = (double)pot.y / (double)pot.x;
						if (abs(kThree - k) == kmin)
						{
							r2 = sqrt(pot.x * pot.x + pot.y *pot.y);
							break;
						}
					}


					a_listFlt.push_back((r1 + r2));
				}
			}
		}

		return TRUE;
	}

	BOOL CGBImgPropCal::GetDisFrom4list(std::vector<CPoint> a_list1, std::vector<CPoint> a_list2,
		std::vector<CPoint> a_list3, std::vector<CPoint> a_list4, std::vector<double>& a_listFlt)
	{
		if ((int)a_list1.size() < ANGLE_MAX ||
			(int)a_list2.size() < ANGLE_MAX ||
			(int)a_list3.size() < ANGLE_MAX ||
			(int)a_list4.size() < ANGLE_MAX)
		{
			return FALSE;
		}
		//there should be enough radius in one and two quadrant
		for (int i = 0; i < ANGLE_MAX; i++)
		{
			// one and three
			double k = TAN_VALUES[i];
			double r1 = 0;
			double r2 = 0;

			r1 = GetRadiusAsK(a_list1, k);
			r2 = GetRadiusAsK(a_list3, k);

			if(r1!=0 && r2!=0)
				a_listFlt.push_back((r1 + r2));
			// two and four
			k = -TAN_VALUES[i];
			r1 = GetRadiusAsK(a_list2, k);
			r2 = GetRadiusAsK(a_list4, k);

			if (r1 != 0 && r2 != 0)
				a_listFlt.push_back((r1 + r2));
			
		}

		return TRUE;
	}

	BOOL CGBImgPropCal::GetDisFrom13list(std::vector<CPoint> a_list1,
		std::vector<CPoint> a_list3,std::vector<double>& a_listFlt)
	{
		if ((int)a_list1.size() < ANGLE_MAX ||			
			(int)a_list3.size() < ANGLE_MAX )
		{
			return FALSE;
		}
		//there should be enough radius in one and two quadrant
		for (int i = 0; i < ANGLE_MAX; i++)
		{
			// one and three
			double k = TAN_VALUES[i];
			double r1 = 0;
			double r2 = 0;

			r1 = GetRadiusAsK(a_list1, k);
			r2 = GetRadiusAsK(a_list3, k);

			if (r1 != 0 && r2 != 0)
				a_listFlt.push_back((r1 + r2));		

		}

		return TRUE;
	}

	BOOL CGBImgPropCal::GetDisFrom24list(std::vector<CPoint> a_list2,
		std::vector<CPoint> a_list4, std::vector<double>& a_listFlt)
	{
		if ((int)a_list2.size() < ANGLE_MAX ||
			(int)a_list4.size() < ANGLE_MAX)
		{
			return FALSE;
		}
		//there should be enough radius in one and two quadrant
		for (int i = 0; i < ANGLE_MAX; i++)
		{
			
			// two and four
			double k = -TAN_VALUES[i];
			double r1 = GetRadiusAsK(a_list2, k);
			double r2 = GetRadiusAsK(a_list4, k);

			if (r1 != 0 && r2 != 0)
				a_listFlt.push_back((r1 + r2));

		}

		return TRUE;
	}

	BOOL CGBImgPropCal::GetSuperParticleSize(COTSParticlePtr a_pOTSPart, double &a_dMinWidth, double &a_dMaxLength, double &a_dAveDia)
	{
		// safety check
		ASSERT(a_pOTSPart);
		

		// 1. convert particle data to image data
		// get the segment list
		COTSSegmentsList listSegment = a_pOTSPart->GetFeature()->GetSegmentsList();
		// get rectangle of the particle
		CRect rect = a_pOTSPart->GetParticleRect();

		//get the particle rectangle size
		a_dMinWidth = rect.Width();
		a_dMaxLength = rect.Height();
		a_dAveDia = (a_dMinWidth + a_dMaxLength) / 2;

		if (a_dMinWidth > a_dMaxLength)
		{
			a_dMaxLength = a_dMinWidth;
		}

		return TRUE;
	}

	// note: this not include touch 
	double CGBImgPropCal::GetDisFromRLNeigbor(std::vector<CPoint> a_Yaxislist, std::vector<CPoint> a_YaxisRightlist, std::vector<CPoint> a_YaxisLeftlist)
	{
		double d = 0;

		int nXMax,nXMin;
		GetMaxMin(a_Yaxislist, AXIAS_DIRECTION::XAX, nXMax, nXMin);
		int nXRMax, nXRMin;
		GetMaxMin(a_YaxisRightlist, AXIAS_DIRECTION::XAX, nXRMax, nXRMin);
		int nXLMax, nXLMin;
		GetMaxMin(a_YaxisLeftlist, AXIAS_DIRECTION::XAX, nXLMax, nXLMin);

		if (nXMax*nXMin<0 && (nXLMax*nXLMin < 0 || nXRMax*nXLMin < 0))
		{
			d = abs(nXMax - nXMin);
		}

		return d;
		
	}
	
	double CGBImgPropCal::GetDisFromUDNeigbor(std::vector<CPoint> a_Xaxislist, std::vector<CPoint> a_XaxisUplist, std::vector<CPoint> a_XaxisDownlist)
	{
		double d = 0;

		int nYMax, nYMin;
		GetMaxMin(a_Xaxislist, AXIAS_DIRECTION::YAX, nYMax, nYMin);
		int nYUMax, nYUMin;
		GetMaxMin(a_XaxisUplist, AXIAS_DIRECTION::YAX, nYUMax, nYUMin);
		int nYDMax, nYDMin;
		GetMaxMin(a_XaxisDownlist, AXIAS_DIRECTION::YAX, nYDMax, nYDMin);

		if (nYMax*nYMin<0 && (nYUMax*nYUMin < 0 || nYDMax*nYDMin < 0))
		{
			d = abs(nYMax - nYMin);
		}

		return d;
	}

	BOOL CGBImgPropCal::GetMaxMin(std::vector<CPoint> a_Xaxislist, AXIAS_DIRECTION a_nDirection, int& a_nMax, int& a_nMin)
	{
		int nMax = 0;
		int nMin = 100;

		if (a_nDirection == AXIAS_DIRECTION::YAX)
		{			
			for (auto pt : a_Xaxislist)
			{
				if (pt.x > nMax)
					nMax = pt.x;
				if (pt.x < nMin)
					nMin = pt.x;
			}			
		}
		else if (a_nDirection == AXIAS_DIRECTION::XAX)
		{
			for (auto pt : a_Xaxislist)
			{
				if (pt.y > nMax)
					nMax = pt.y;
				if (pt.y < nMin)
					nMin = pt.y;
			}
		}

		a_nMax = nMax;
		a_nMin = nMin;

		return TRUE;
	}

	double CGBImgPropCal::GetRadiusAsK(std::vector<CPoint> a_Ptlist, double a_k)
	{
		
		double r1 = 0;
		
		vector<CPoint> ptLengthList;
		for (auto pt : a_Ptlist)
		{
			if (abs(pt.x * a_k - pt.y) < VALUE_MIN)
			{
				ptLengthList.push_back(pt);
			}
		}

		if (ptLengthList.empty())
			return 0;

		double rAve = 0;
		for (auto len : ptLengthList)
		{
			double r = sqrt(len.x * len.x + len.y *len.y);
			rAve += r;
		}

		rAve = rAve / (double)ptLengthList.size();

		return rAve;
	}
}