#include "ZImage.h"
#include "ZLogger.h"
#include <iostream>
#include <fstream>
#include <atlimage.h>
#include <regex>

#define ANGLE(b, w) (b < 0 ? (b + w) : (b >= w ? (b - w) : b))

void quicksort(int* a, double* b, int l, int r)
{
	int left, right, i, j, tmp;
	double pivot, x, y, z;
	int taskl[1000];
	int taskr[1000];
	int taskIndex = 0;
	taskl[taskIndex] = l;
	taskr[taskIndex] = r;
	while (taskIndex >= 0) {
		left = taskl[taskIndex];
		right = taskr[taskIndex];
		taskIndex--;

		if (left < right) {
			i = left;
			j = right;
			x = b[a[i]];
			y = b[a[i + (j - i) / 2]];
			z = b[a[j]];
			if (x < y)
				if (y < z) pivot = y;
				else if (z < x) pivot = x;
				else pivot = z;
			else if (z < y) pivot = y;
			else if (x < z) pivot = x;
			else pivot = z;

			while (true) {
				while (b[a[i]] < pivot) i++;
				while (pivot < b[a[j]]) j--;
				if (i >= j) break;
				tmp = a[i]; a[i] = a[j]; a[j] = tmp;
				i++;
				j--;
			}

			taskIndex++;
			taskl[taskIndex] = j + 1;
			taskr[taskIndex] = right;

			taskIndex++;
			taskl[taskIndex] = left;
			taskr[taskIndex] = i - 1;
		}
	}
}

void quicksort(int* a, int* b, int l, int r)
{
	int left, right, i, j, tmp;
	int pivot, x, y, z;
	int taskl[1000];
	int taskr[1000];
	int taskIndex = 0;
	taskl[taskIndex] = l;
	taskr[taskIndex] = r;
	while (taskIndex >= 0) {
		left = taskl[taskIndex];
		right = taskr[taskIndex];
		taskIndex--;

		if (left < right) {
			i = left;
			j = right;
			x = b[a[i]];
			y = b[a[i + (j - i) / 2]];
			z = b[a[j]];
			if (x < y)
				if (y < z) pivot = y;
				else if (z < x) pivot = x;
				else pivot = z;
			else if (z < y) pivot = y;
			else if (x < z) pivot = x;
			else pivot = z;

			while (true) {
				while (b[a[i]] < pivot) i++;
				while (pivot < b[a[j]]) j--;
				if (i >= j) break;
				tmp = a[i]; a[i] = a[j]; a[j] = tmp;
				i++;
				j--;
			}

			taskIndex++;
			taskl[taskIndex] = j + 1;
			taskr[taskIndex] = right;

			taskIndex++;
			taskl[taskIndex] = left;
			taskr[taskIndex] = i - 1;
		}
	}
}

struct Frame {
	int Frame;
	int CellCount;
	double CellSize;
	double CellCircularity;
	int EukaryoteCount;
	double EukaryoteSize;
	double EukaryoteCircularity;
	double NucleiSize;
	double NucleiPart;
	double NucleoliSize;
	double NucleoliPart;
	double LinkCount;
	double LinkLength;
	int OrganelleCount;
	int CellWithOrganelleCount;
	double OrganelleTotalSize;
	double OrganelleMaxSize;
};

struct Cell {
	int MinX;
	int MinY;
	int MaxX;
	int MaxY;
	double MinDelta;
	double MaxDelta;

	int QueStart;
	int Size;
	int Line;
	bool IsOuter;

	int OrganelleCount;
	double OrganelleTotalSize;
	double OrganelleMaxSize;

	bool IsChecked;

	int NucleiCenter;
	double NucleiHume;
	double NucleiSize;
	double NucleoliSize;
	double NucleiLengthes[180];
	double NucleoliLengthes[180];

	int LinkCount;
	int LinkFrom;
	int Id;
	double LinkLength;
	double Aspect;
};

struct Part {
	int Id;
	int Position;
	int Size;
	bool IsCandidate;
	double Hume;
	int Center;
	double Lengthes[180];
	double Area;
};

struct View {
	double NucleiLength;
	double NucleoliLength;
	double NucleiMoving;
	double NucleoliMoving;
	double Sin;
	double Cos;
};

struct Link {
	int Frame;
	int Id;
	int From;
	int To;
};

int main(int argc, char* argv[])
{
	// System
	int logLevel = 0;
	const int loopCount = 1;
	const int left = 0, top = 0;
	char filename[256];
	char tempname[16];

	// Arguments
	int width = 960, height = 720;
	const char* imagesPath = "";
	const char* outputPath = "";
	int frameStart = 1, frameCount = 5, frameLength = 4;
	int stackStart = 1, stackCount = 16, stackLength = 3, stackBlur = 4;
	bool outputResultImage = false;
	double thresthold_ditect_cell = 0.1;
	double thresthold_nuclei_hume = 0.6;
	double thresthold_nuclei_valley = 0.1;
	int thresthold_enough_cell = 800;
	int thresthold_minimun_cell = 200;
	double circler_distance = 16.0;
	int circler_width = 30;
	int circler_cutter = 3;
	int imageZoom = 1;

	int mode_type = 0;
	for (int i = 0; i < argc; i++) {
		if (mode_type == 3) width = std::atoi(argv[i]);
		if (mode_type == 4) height = std::atoi(argv[i]);
		if (mode_type == 7) imagesPath = argv[i];
		if (mode_type == 8) outputPath = argv[i];
		if (mode_type == 10) frameStart = std::atoi(argv[i]);
		if (mode_type == 11) frameCount = std::atoi(argv[i]);
		if (mode_type == 12) frameLength = std::atoi(argv[i]);
		if (mode_type == 20) stackStart = std::atoi(argv[i]);
		if (mode_type == 21) stackCount = std::atoi(argv[i]);
		if (mode_type == 22) stackLength = std::atoi(argv[i]);
		if (mode_type == 30) thresthold_nuclei_hume = std::atof(argv[i]);
		if (mode_type == 31) thresthold_ditect_cell = std::atof(argv[i]);
		if (mode_type == 32) thresthold_nuclei_valley = std::atof(argv[i]);
		if (mode_type == 40) stackBlur = std::atoi(argv[i]);
		if (mode_type == 0) {
			if (strcmp(argv[i], "-w") == 0) mode_type = 3;
			if (strcmp(argv[i], "-h") == 0) mode_type = 4;
			if (strcmp(argv[i], "-p") == 0) mode_type = 7;
			if (strcmp(argv[i], "-o") == 0) mode_type = 8;
			if (strcmp(argv[i], "-fs") == 0) mode_type = 10;
			if (strcmp(argv[i], "-c") == 0) mode_type = 11;
			if (strcmp(argv[i], "-fl") == 0) mode_type = 12;
			if (strcmp(argv[i], "-ss") == 0) mode_type = 20;
			if (strcmp(argv[i], "-sc") == 0) mode_type = 21;
			if (strcmp(argv[i], "-sl") == 0) mode_type = 22;
			if (strcmp(argv[i], "-th") == 0) mode_type = 30;
			if (strcmp(argv[i], "-tc") == 0) mode_type = 31;
			if (strcmp(argv[i], "-tv") == 0) mode_type = 32;
			if (strcmp(argv[i], "-i") == 0) outputResultImage = true;
			if (strcmp(argv[i], "-b") == 0) mode_type = 40;
		} else mode_type = 0;
	}

	// Outputs
	sprintf_s(filename, "%s\\_cells.csv", outputPath);
	std::ofstream output1(filename);
	sprintf_s(filename, "%s\\_average.csv", outputPath);
	std::ofstream output2(filename);
	ZImage image(width, height, imageZoom);
	ZLogger logger(logLevel);

	output1 << "Index" << ", ";
	output1 << "Frame" << ", ";
	output1 << "Id" << ", ";
	output1 << "Has Nuclei" << ", ";
	output1 << "Cell Size" << ", ";
	output1 << "Cell Circularity" << ", ";
	output1 << "Move" << ", ";
	output1 << "Min X" << ", ";
	output1 << "Min Y" << ", ";
	output1 << "Max X" << ", ";
	output1 << "Max Y" << ", ";
	output1 << "Aspect" << ", ";
	output1 << "Nuclei Size" << ", ";
	output1 << "Nuclei Value" << ", ";
	output1 << "Nuclei / Cell" << ", ";
	output1 << "Nucleoi Size" << ", ";
	output1 << "Nucleoli / Nuclei" << ", ";
	output1 << "Vacuole Count" << ", ";
	output1 << "Vacuole Size" << ", ";
	output1 << "Vacuole Max" << ", ";
	output1 << std::endl;

	output2 << "Index" << ", ";
	output2 << "Frame" << ", ";
	output2 << "Cell Count (All)" << ", ";
	output2 << "Cell Size (All)" << ", ";
	output2 << "Circularity (All)" << ", ";
	output2 << "Move" << ", ";
	output2 << "Cell Count (Nuclei)" << ", ";
	output2 << "Cell Size (Nuclei)" << ", ";
	output2 << "Circularity (Nuclei)" << ", ";
	output2 << "Nuclei Size" << ", ";
	output2 << "Nuclei / Cell" << ", ";
	output2 << "Nucleoli Size" << ", ";
	output2 << "Nucleoli / Nuclei" << ", ";
	output2 << "Vacuole Count" << ", ";
	output2 << "Vacuole Size" << ", ";
	output2 << "Vacuole Max" << ", ";
	output2 << std::endl; 

	// Memory
	unsigned char* sourceImage = (unsigned char*)malloc(sizeof(unsigned char) * stackCount * width * height);
	double* deltaMap = (double*)malloc(sizeof(double) * width * height);
	int* sortMap = (int*)malloc(sizeof(int) * width * height);
	int* tempMap = (int*)malloc(sizeof(int) * width * height);
	int* taskData = (int*)malloc(sizeof(int) * width * height);
	View* views = (View*)malloc(sizeof(View) * 181);

	const int resultCount = frameCount + 1;
	Frame* frames = (Frame*)malloc(sizeof(Frame) * resultCount);

	int* cellsMap = (int*)malloc(sizeof(int) * width * height);
	int* cellsQue = (int*)malloc(sizeof(int) * width * height);
	Cell* cellsData = (Cell*)malloc(sizeof(Cell) * 1000);

	int cellsResultCount = 0;
	int* cellsResultMap = (int*)malloc(sizeof(int) * width * height);
	Cell* cellsResultData = (Cell*)malloc(sizeof(Cell) * 1000);

	int* partsMap = (int*)malloc(sizeof(int) * width * height);
	int* partsQue = (int*)malloc(sizeof(int) * width * height);
	Part* partsData = (Part*)malloc(sizeof(Part) * 100000);

	Link* links = (Link*)malloc(sizeof(Link) * 500000);

	if (sourceImage != NULL && deltaMap != NULL && sortMap != NULL && cellsMap != NULL && cellsQue != NULL && cellsData != NULL && cellsResultData != NULL && cellsResultMap != NULL && partsMap != NULL && partsData != NULL && partsQue != NULL && tempMap != NULL && taskData != NULL && views != NULL && frames != NULL && links != NULL) {

		int wholeTotalId = 1;
		int wholeLinkId = 0;
		int targets[4];

		const int size = width * height;
		const int deltas[] { -1, 1, -width, width, -1 - width, -1 + width, 1 - width, 1 + width, };
		const int carving = 360 / circler_cutter;
		const int moving = circler_width / circler_cutter;
		for (int a = 0; a < carving + 1; a++) {
			views[a].Sin = sin(a * 360.0 / carving / 180.0 * 3.1415927);
			views[a].Cos = cos(a * 360.0 / carving / 180.0 * 3.1415927);
		}

		int viewMax = 10;
		int viewCount = 0;
		for (int x = -viewMax; x <= viewMax; x++) {
			for (int y = -viewMax; y <= viewMax; y++) {
				sortMap[viewCount] = viewCount;
				deltaMap[viewCount++] = sqrt(x * x + y * y);
			}
		}
		quicksort(sortMap, deltaMap, 0, viewCount - 1);
		char rdx[500];
		char rdy[500];
		for (int i = 0; i < viewCount; i++) {
			int position = sortMap[i];
			if (deltaMap[position] > viewMax) {
				viewCount = i;
				break;
			}
			rdx[i] = position % (viewMax * 2 + 1) - viewMax;
			rdy[i] = position / (viewMax * 2 + 1) - viewMax;
		}
		
		// Analyze frames
		for (int imageIndex = frameStart; imageIndex < frameStart + frameCount; imageIndex++)
		{
			int resultIndex = imageIndex - frameStart + 1;
			if (logger.infomation("Frame Start", true)) continue;

			// Load image
			for (int i = 0; i < stackCount; i++) {
				sprintf_s(tempname, "%%0%dd", frameLength);
				sprintf_s(tempname, tempname, imageIndex);
				sprintf_s(filename, "%s", std::regex_replace(imagesPath, std::regex("%0"), tempname).c_str());
				sprintf_s(tempname, "%%0%dd", stackLength);
				sprintf_s(tempname, tempname, i + stackStart);
				sprintf_s(filename, "%s", std::regex_replace(filename, std::regex("%1"), tempname).c_str());
				CImage originalImage;
				if (originalImage.Load(CA2W(filename)) == S_OK) {
					for (int y = 0; y < height; y++) {
						memcpy(sourceImage + static_cast<__int64>(i) * size + static_cast<__int64>(y) * width, (unsigned char*)originalImage.GetBits() + (static_cast<__int64>(y) + top) * originalImage.GetPitch() + left, width);
					}
					originalImage.Destroy();
				}
				else logger.error("Image Load Failed.");
			}
			if (logger.infomation("Image Loaded", true)) continue;

			// Calculate Z stack images
			int tempLeft = 0, tempRight = stackCount - 1, tempCenter = -1;
			for (int loop = 0; loop < loopCount; loop++) {

				// Center frame
				int distance = 0;
				for (int j = tempLeft; j <= tempRight; j++) {
					int mean = 0, offset = j * size;
					for (int i = 0; i < size; i++) {
						int center = sourceImage[offset + i];
						if (i % width > 0) mean += (center - sourceImage[offset + i - 1]) * (center - sourceImage[offset + i - 1]);
						if (i / width > 0) mean += (center - sourceImage[offset + i - width]) * (center - sourceImage[offset + i - width]);
					}
					if (distance < mean) {
						distance = mean;
						tempCenter = j;
					}
				}
				if (logger.step(10)) break;

				// Fix Z stack
				if (stackBlur > 0) {
					double sourceMin = 255, sourceMax = 0;
					for (int i = 0; i < size; i++) {
						unsigned char value = sourceImage[tempCenter * size + i];
						if (sourceMin > value) sourceMin = value;
						if (sourceMax < value) sourceMax = value;
					}

					double distance = DBL_MAX;
					for (int l = 0; l <= stackBlur; l++) {
						for (int dr = 0; dr <= stackBlur - l; dr++) {

							int r = stackCount - 1 - dr, count = r - l + 1;
							double ax = (0.5 * r * (1.0 + r) - 0.5 * (1.0 + l) * l) / count;
							double peek = 0;
							for (int i = 0; i < size; i++) {
								double ay = 0;
								for (int j = l; j <= r; j++) ay += sourceImage[j * size + i];
								ay /= count;

								double sxx = 0, sxy = 0;
								for (int j = l; j <= r; j++) {
									sxx += (j - ax) * (j - ax);
									sxy += (sourceImage[j * size + i] - ay) * (j - ax);
								}
								double a = sxy / sxx;

								deltaMap[i] = a;
								if (peek < abs(a)) peek = abs(a);
							}

							double totalMean1 = 0;
							for (int i = 0; i < size; i++) {
								double temp1 = (sourceImage[tempCenter * size + i] - sourceMin) / (sourceMax - sourceMin) - (1 - deltaMap[i] / peek) / 2.0;
								totalMean1 += temp1 * temp1;
								deltaMap[i] = temp1 * temp1;
							}

							double totalMean3 = 0;
							for (int i = 0; i < size; i++) {
								double center = deltaMap[i];
								if (i % width > 0) totalMean3 += (center - deltaMap[i - 1]) * (center - deltaMap[i - 1]);
								if (i / width > 0) totalMean3 += (center - deltaMap[i - width]) * (center - deltaMap[i - width]);
							}

							if (distance > totalMean3) {
								distance = totalMean3;
								tempLeft = l;
								tempRight = r;
							}
						}
					}
				}
				if (logger.step(11)) break;

				// Calculate delta
				int count = tempRight - tempLeft + 1;
				double peek = 0, ax = (0.5 * tempRight * (1.0 + tempRight) - 0.5 * (1.0 + tempLeft) * tempLeft) / count;
				for (int i = 0; i < size; i++) {
					double ay = 0;
					for (int j = tempLeft; j <= tempRight; j++) ay += sourceImage[j * size + i];
					ay /= count;

					double sxx = 0, sxy = 0;
					for (int j = tempLeft; j <= tempRight; j++) {
						sxx += (j - ax) * (j - ax);
						sxy += (sourceImage[j * size + i] - ay) * (j - ax);
					}
					double a = sxy / sxx;

					deltaMap[i] = a;
					if (peek < abs(a)) peek = abs(a);
				}
				for (int i = 0; i < size; i++) deltaMap[i] /= peek;
				if (logger.step(12)) break;
			}
			if (logger.infomation("[10] Stack Selected", true, 10)) continue;
			if (logger.infomation("[11] Stack Searched", true, 11)) continue;
			if (logger.infomation("[12] Delta Calculated", true, 12)) continue;

			// Detect all areas
			for (int loop = 0; loop < loopCount; loop++) {
				for (int i = 0; i < size; i++) sortMap[i] = i;
				quicksort(sortMap, deltaMap, 0, size - 1);
				if (logger.step(40)) break;

				memset(tempMap, 0, sizeof(int) * size);

				int centerPosition;
				for (centerPosition = 0; centerPosition < size; centerPosition++) if (deltaMap[sortMap[centerPosition]] > 0) break;
				for (int i = size - 1; i >= centerPosition; i--) {
					int position = sortMap[i];
					if (deltaMap[position] >= thresthold_ditect_cell
						|| (position % width > 0 && tempMap[position - 1] == 1)
						|| (position % width < width - 1 && tempMap[position + 1] == 1)
						|| (position >= width && tempMap[position - width] == 1)
						|| (position < size - width - 1 && tempMap[position + width] == 1)) tempMap[position] = 1;
				}
				for (int i = 0; i < centerPosition; i++) {
					int position = sortMap[i];
					if (deltaMap[position] <= -thresthold_ditect_cell
						|| (position % width > 0 && tempMap[position - 1] == 2)
						|| (position % width < width - 1 && tempMap[position + 1] == 2)
						|| (position >= width && tempMap[position - width] == 2)
						|| (position < size - width - 1 && tempMap[position + width] == 2)) tempMap[position] = 2;
				}

				for (int i = 0; i < size; i++) {
					int position = sortMap[i];
					if (tempMap[position] != 0) continue;
					if (deltaMap[position] >= 0) {
						if ((position % width > 0 && tempMap[position - 1] == 1) ||
							(position % width < width - 1 && tempMap[position + 1] == 1) ||
							(position >= width && tempMap[position - width] == 1) ||
							(position < size - width - 1 && tempMap[position + width] == 1)) tempMap[position] = 3;
					}
					else {
						if ((position % width > 0 && tempMap[position - 1] == 2) ||
							(position % width < width - 1 && tempMap[position + 1] == 2) ||
							(position >= width && tempMap[position - width] == 2) ||
							(position < size - width - 1 && tempMap[position + width] == 2)) tempMap[position] = 3;
					}
				}

				for (int i = 0; i < size; i++) {
					if (tempMap[i] != 0) continue;
					bool touched = false;
					int taskIndex = 0, taskCount = 0;
					tempMap[i] = -1;
					taskData[taskCount++] = i;
					while (taskIndex < taskCount)
					{
						int position = taskData[taskIndex++];
						int targetCount = 0;
						if (position % width > 0) targets[targetCount++] = position - 1;
						if (position % width < width - 1) targets[targetCount++] = position + 1;
						if (position >= width) targets[targetCount++] = position - width;
						if (position < size - width - 1) targets[targetCount++] = position + width;
						if (targetCount != 4) touched = true;
						for (int j = 0; j < targetCount; j++) {
							int targetPosition = targets[j];
							if (tempMap[targetPosition] != 0) continue;
							tempMap[targetPosition] = -1;
							taskData[taskCount++] = targetPosition;
						}
					}
					if (!touched) for (int j = 0; j < taskCount; j++) tempMap[taskData[j]] = 3;
				}

				for (int i = 0; i < size; i++) {
					if (tempMap[i] != 3 || deltaMap[i] < 0) continue;
					bool touched = false;
					int taskIndex = 0, taskCount = 0;
					tempMap[i] = 4;
					taskData[taskCount++] = i;
					while (taskIndex < taskCount)
					{
						int position = taskData[taskIndex++];
						int targetCount = 0;
						if (position % width > 0) targets[targetCount++] = position - 1;
						if (position % width < width - 1) targets[targetCount++] = position + 1;
						if (position >= width) targets[targetCount++] = position - width;
						if (position < size - width - 1) targets[targetCount++] = position + width;
						for (int j = 0; j < targetCount; j++) {
							int targetPosition = targets[j];
							if (tempMap[targetPosition] == 1) touched = true;
							if (tempMap[targetPosition] != 3 || deltaMap[targetPosition] < 0) continue;
							tempMap[targetPosition] = 4;
							taskData[taskCount++] = targetPosition;
						}
					}
					if (!touched) for (int j = 0; j < taskCount; j++) tempMap[taskData[j]] = 5;
				}
				if (logger.step(41)) break;
			}
			if (logger.infomation("[40] Sorted", true, 40)) continue;
			if (logger.infomation("[41] Colored", true, 41)) continue;

			// Detect cells
			int tempCount = 0;
			int cellsCount = 1;
			int cellsQueCount = 0;
			for (int loop = 0; loop < loopCount; loop++) {
				for (int i = 0; i < size; i++) cellsMap[i] = (tempMap[i] == 1 || tempMap[i] == 4) ? 1 : 0;

				while (true)
				{
					int count = 0;
					for (int i = 0; i < size; i++) {
						if (cellsMap[i] != 1) continue;
						int x = i % width;
						int y = i / width;
						if ((0 < x && x < width - 1 && cellsMap[i - 1] != 1 && cellsMap[i + 1] != 1) ||
							(0 < y && y < height - 1 && cellsMap[i - width] != 1 && cellsMap[i + width] != 1) ||
							((0 < x && x < width - 1 && 0 < y && y < height - 1) && (
								(cellsMap[i - 1 - width] != 1 && cellsMap[i + 1 + width] != 1) ||
								(cellsMap[i + 1 - width] != 1 && cellsMap[i - 1 + width] != 1)))
							) 
						{
							cellsMap[i] = 2;
							count++;
						}
					}
					if (count == 0) break;
				}

				for (int j = 3; j <= 4; j++) {
					for (int i = 0; i < size; i++) {
						if (cellsMap[i] == 0 || cellsMap[i] == j) continue;
						int x = i % width;
						int y = i / width;

						if (x > 0 && cellsMap[i - 1] == 0) cellsMap[i - 1] = j;
						if (x < width - 1 && cellsMap[i + 1] == 0) cellsMap[i + 1] = j;
						if (y > 0 && cellsMap[i - width] == 0) cellsMap[i - width] = j;
						if (y < height - 1 && cellsMap[i + width] == 0) cellsMap[i + width] = j;
					}
				}

				memset(tempMap, 0, sizeof(int) * size);
				partsData[0].Id = 0;
				partsData[0].Size = 0;
				for (int i = size - 1; i >= 0; i--) {
					int position = sortMap[i];
					if (cellsMap[position] == 0) continue;

					int targetCount = 0;
					if (position % width > 0) targets[targetCount++] = position - 1;
					if (position % width < width - 1) targets[targetCount++] = position + 1;
					if (position >= width) targets[targetCount++] = position - width;
					if (position < size - width - 1) targets[targetCount++] = position + width;

					int cellId = INT_MAX;
					for (int j = 0; j < targetCount; j++) {
						int targetPosition = targets[j];
						int targetId = tempMap[targetPosition];
						while (targetId != partsData[targetId].Id) targetId = partsData[targetId].Id;
						tempMap[targetPosition] = targetId;
						if (targetId != 0 && cellId > targetId) cellId = targetId;
					}

					if (cellId == INT_MAX) {
						if (cellsMap[position] == 1 || cellsMap[position] == 2) {
							tempMap[position] = ++tempCount;
							partsData[tempCount].Id = tempCount;
							partsData[tempCount].Size = 1;
						}
					}
					else {
						tempMap[position] = cellId;
						partsData[cellId].Size++;

						for (int j = 0; j < targetCount; j++) {
							int targetId = tempMap[targets[j]];
							if (targetId != 0 && targetId != cellId) {
								if (cellsMap[position] > 3 || cellsMap[targets[j]] > 3) continue;
								if ((cellsMap[position] != 1 || cellsMap[targets[j]] != 1) && partsData[cellId].Size >= 100 && partsData[targetId].Size >= 100) continue;
								partsData[cellId].Size += partsData[targetId].Size;
								partsData[targetId].Id = cellId;
								partsData[targetId].Size = 0;
							}
						}
					}
				}

				for (int i = 0; i < size; i++) {
					int id = tempMap[i];
					while (id != partsData[id].Id) id = partsData[id].Id;
					tempMap[i] = partsData[id].Size < 100 ? 0 : id;
				}

				for (int i = size - 1; i >= 0; i--) {
					int position = sortMap[i];

					int partId = tempMap[position];
					if (partId == 0) continue;
					while (partId != partsData[partId].Id) partId = partsData[partId].Id;
					tempMap[position] = partId;

					int targetCount = 0;
					if (position % width > 0) targets[targetCount++] = position - 1;
					if (position % width < width - 1) targets[targetCount++] = position + 1;
					if (position >= width) targets[targetCount++] = position - width;
					if (position < size - width - 1) targets[targetCount++] = position + width;
					for (int j = 0; j < targetCount; j++) {
						int targetId = tempMap[targets[j]];
						if (targetId == 0) continue;
						while (targetId != partsData[targetId].Id) targetId = partsData[targetId].Id;
						tempMap[targets[j]] = targetId;
						if (targetId != partId) {
							if (partsData[partId].Size >= thresthold_enough_cell && partsData[targetId].Size >= thresthold_enough_cell) continue;
							partsData[partId].Size += partsData[targetId].Size;
							partsData[targetId].Id = partId;
							partsData[targetId].Size = 0;
						}
					}
				}

				partsData[0].Position = 0;
				for (int i = 1; i <= tempCount; i++) {
					if (partsData[i].Size >= thresthold_minimun_cell) {
						cellsData[cellsCount].Size = 0;
						cellsData[cellsCount].Id = cellsCount;
						cellsData[cellsCount].LinkCount = 0;
						cellsData[cellsCount].LinkFrom = 0;
						partsData[i].Position = cellsCount++;
					} else partsData[i].Position = 0;
				}
				for (int i = 0; i < size; i++) {
					int id = tempMap[i];
					while (id != partsData[id].Id) id = partsData[id].Id;
					tempMap[i] = partsData[id].Position;
				}
				if (logger.step(42)) break;

				// Link
				if (cellsResultCount > 0) {
					memset(partsQue, 0, sizeof(int) * size);

					for (int i = 1; i < cellsResultCount; i++) {
						cellsResultData[i].LinkCount = 0;
						cellsResultData[i].LinkFrom = 0;
					}

					int* map1 = tempMap;
					int* map2 = cellsResultMap;
					for (int i = 0; i < size; i++) {
						int x = i % width;
						int y = i / width;
						if (map1[i] > 0) {
							for (int j = 0; j < viewCount; j++) {
								int dx = x + rdx[j];
								int dy = y + rdy[j];
								if (dx < 0 || width <= dx || dy < 0 || height <= dy) break;
								if (map2[dy * width + dx] > 0) {
									partsQue[map2[dy * width + dx] * cellsCount + map1[i]]++;
									break;
								}
							}
						}
						if (map2[i] > 0) {
							for (int j = 0; j < viewCount; j++) {
								int dx = x + rdx[j];
								int dy = y + rdy[j];
								if (dx < 0 || width <= dx || dy < 0 || height <= dy) break;
								if (map1[dy * width + dx] > 0) {
									partsQue[map2[i] * cellsCount + map1[dy * width + dx]]++;
									break;
								}
							}
						}
					}
					for (int i = 0; i < cellsResultCount * cellsCount; i++) {
						if (partsQue[i] > 0) {
							int id1 = i % cellsCount;
							int id2 = i / cellsCount;
							cellsData[id1].LinkCount++;
							cellsResultData[id2].LinkCount++;
						}
					}
					for (int i = 0; i < cellsResultCount * cellsCount; i++) {
						if (partsQue[i] > 0) {
							int id1 = i % cellsCount;
							int id2 = i / cellsCount;
							if (cellsData[id1].LinkCount == 1 || cellsResultData[id2].LinkCount == 1) {
								cellsData[id1].LinkFrom = id2;
								cellsResultData[id2].LinkFrom = id1;
								partsQue[i] = -1;
							}
						}
					}
					for (int i = 0; i < cellsResultCount * cellsCount; i++) {
						if (partsQue[i] > 0) {
							int id1 = i % cellsCount;
							int id2 = i / cellsCount;
							if (cellsData[id1].LinkFrom > 0 && cellsResultData[id2].LinkFrom > 0) {
								partsQue[i] = -2;
							}
						}
					}

					int linkCount = 0;
					for (int i = 0; i < cellsResultCount * cellsCount; i++) {
						if (partsQue[i] > 0) {
							partsMap[linkCount++] = i;
						}
					}
					quicksort(partsMap, partsQue, 0, linkCount - 1);

					for (int i = linkCount - 1; i >= 0; i--) {
						int id1 = partsMap[i] % cellsCount;
						int id2 = partsMap[i] / cellsCount;
						if (cellsData[id1].LinkFrom == 0 && cellsResultData[id2].LinkFrom == 0) {
							cellsData[id1].LinkFrom = id2;
							cellsResultData[id2].LinkFrom = id1;
							partsQue[partsMap[i]] = -1;
						}
					}

					for (int i = linkCount - 1; i >= 0; i--) {
						int id1 = partsMap[i] % cellsCount;
						int id2 = partsMap[i] / cellsCount;
						if (cellsData[id1].LinkFrom == 0 || cellsResultData[id2].LinkFrom == 0) {
							cellsData[id1].LinkFrom = id2;
							cellsResultData[id2].LinkFrom = id1;
							partsQue[partsMap[i]] = -1;
						}
					}
				}
				if (logger.step(45)) break;
				
				// Join
				for (int i = size - 1; i >= 0; i--) {
					int position = sortMap[i];
					int cellId = tempMap[position];
					if (cellId == 0) continue;
					while (cellId != cellsData[cellId].Id) cellId = cellsData[cellId].Id;
					tempMap[position] = cellId;
					if (cellsData[cellId].LinkFrom == 0) continue;

					int targetCount = 0;
					if (position % width > 0) targets[targetCount++] = position - 1;
					if (position % width < width - 1) targets[targetCount++] = position + 1;
					if (position >= width) targets[targetCount++] = position - width;
					if (position < size - width - 1) targets[targetCount++] = position + width;
					for (int j = 0; j < targetCount; j++) {
						int targetId = tempMap[targets[j]];
						if (targetId == 0) continue;
						while (targetId != cellsData[targetId].Id) targetId = cellsData[targetId].Id;
						tempMap[targets[j]] = targetId;
						if (targetId != cellId && cellsData[cellId].LinkFrom == cellsData[targetId].LinkFrom) cellsData[targetId].Id = cellId;
					}
				}
				for (int i = 0; i < size; i++) {
					int id = tempMap[i];
					if (id == 0) continue;
					while (id != cellsData[id].Id) id = cellsData[id].Id;
					tempMap[i] = id;
				}
				if (cellsResultCount > 0) {
					for (int i = 1; i < cellsCount; i++) {
						cellsData[i].LinkCount = 0;
						cellsData[i].LinkFrom = 0;
					}
					for (int i = 1; i < cellsResultCount; i++) {
						cellsResultData[i].LinkCount = 0;
						cellsResultData[i].LinkFrom = 0;
					}
					for (int i = 0; i < cellsResultCount * cellsCount; i++) {
						if (partsQue[i] != -1) continue;
						int id1 = i % cellsCount;
						int id2 = i / cellsCount;
						if (cellsData[id1].Id != id1 || cellsResultData[id2].Size == 0) continue;
						cellsData[id1].LinkCount++;
						cellsResultData[id2].LinkCount++;

						links[wholeLinkId].Frame = resultIndex;
						links[wholeLinkId].From = id2;
						links[wholeLinkId].To = id1;
						wholeLinkId++;
					}
					for (int i = 0; i < cellsResultCount * cellsCount; i++) {
						if (partsQue[i] != -1) continue;
						int id1 = i % cellsCount;
						int id2 = i / cellsCount;
						if (cellsData[id1].LinkCount == 1 && cellsResultData[id2].LinkCount == 1) {
							cellsData[id1].LinkFrom = id2;
						}
					}
				}

				// Fill
				for (int i = 0; i < size; i++) {
					if (tempMap[i] != 0) continue;
					int taskIndex = 0, taskCount = 0;
					tempMap[i] = -1;
					taskData[taskCount++] = i;
					int targetId = -1;
					while (taskIndex < taskCount)
					{
						int position = taskData[taskIndex++];
						int targetCount = 0;
						if (position % width > 0) targets[targetCount++] = position - 1;
						if (position % width < width - 1) targets[targetCount++] = position + 1;
						if (position >= width) targets[targetCount++] = position - width;
						if (position < size - width - 1) targets[targetCount++] = position + width;
						for (int j = 0; j < targetCount; j++) {
							int targetPosition = targets[j];
							if (tempMap[targetPosition] > 0) {
								if (targetId == -1)targetId = tempMap[targetPosition];
								if (targetId > 0 && targetId != tempMap[targetPosition])targetId = -2;
							}
							if (tempMap[targetPosition] != 0) continue;
							tempMap[targetPosition] = -1;
							taskData[taskCount++] = targetPosition;
						}
					}
					if (targetId > 0) for (int j = 0; j < taskCount; j++) tempMap[taskData[j]] = targetId;
				}

				while (true)
				{
					int count = 0;
					for (int i = 0; i < size; i++) {
						int id = tempMap[i];
						if (id <= 0) continue;
						int x = i % width;
						int y = i / width;
						if ((0 < x && x < width - 1 && tempMap[i - 1] != id && tempMap[i + 1] != id) || 
							(0 < y && y < height - 1 && tempMap[i - width] != id && tempMap[i + width] != id) ||
							((0 < x && x < width - 1 && 0 < y && y < height - 1) && (
								(tempMap[i - 1 - width] != id && tempMap[i + 1 + width] != id) ||
								(tempMap[i + 1 - width] != id && tempMap[i - 1 + width] != id)))
							)
						{
							tempMap[i] = -1;
							count++;
						}
					}
					if (count == 0) break;
				}

				// Set
				for (int i = 0; i < size; i++) {
					int id = tempMap[i];
					cellsMap[i] = tempMap[i];
					if (id > 0) cellsData[id].Size++;
				}
				for (int i = 1; i < cellsCount; i++) {
					cellsData[i].Id = cellsData[i].LinkFrom == 0 ? wholeTotalId++ : cellsResultData[cellsData[i].LinkFrom].Id;
					cellsData[i].MinX = width;
					cellsData[i].MinY = height;
					cellsData[i].MaxX = 0;
					cellsData[i].MaxY = 0;
					cellsData[i].MinDelta = 1.0;
					cellsData[i].MaxDelta = -1.0;
					cellsData[i].QueStart = cellsQueCount;
					cellsData[i].IsOuter = cellsData[i].Size == 0;
					cellsQueCount += cellsData[i].Size;
					cellsData[i].Size = 0;
					cellsData[i].Line = 0;
					cellsData[i].NucleiSize = 0;
					cellsData[i].NucleoliSize = 0;
					cellsData[i].NucleiHume = 0.0;
					cellsData[i].LinkLength = 0.0;
					cellsData[i].OrganelleCount = 0;
					cellsData[i].OrganelleTotalSize = 0.0;
					cellsData[i].OrganelleMaxSize = 0.0;
				}
				for (int i = size - 1; i >= 0; i--) {
					int position = sortMap[i];
					int id = cellsMap[position];
					if (id == -1) continue;

					int x = position % width;
					int y = position / width;
					if (cellsData[id].MinX > x) cellsData[id].MinX = x;
					if (cellsData[id].MinY > y) cellsData[id].MinY = y;
					if (cellsData[id].MaxX < x) cellsData[id].MaxX = x;
					if (cellsData[id].MaxY < y) cellsData[id].MaxY = y;
					double delta = deltaMap[position];
					if (cellsData[id].MinDelta > delta) cellsData[id].MinDelta = delta;
					if (cellsData[id].MaxDelta < delta) cellsData[id].MaxDelta = delta;
					cellsQue[cellsData[id].QueStart + cellsData[id].Size++] = position;

					if (position % width <= 0) cellsData[id].IsOuter = true; else if (cellsMap[position - 1] != id) cellsData[id].Line++;
					if (position % width >= width - 1) cellsData[id].IsOuter = true; else if (cellsMap[position + 1] != id) cellsData[id].Line++;
					if (position < width)cellsData[id].IsOuter = true; else if (cellsMap[position - width] != id) cellsData[id].Line++;
					if (position >= size - width - 1) cellsData[id].IsOuter = true; else if (cellsMap[position + width] != id) cellsData[id].Line++;
				}
				if (logger.step(43)) break;

				// Bounding Box
				for (int cellId = 1; cellId < cellsCount; cellId++) {
					for (int k = 0; k < 60; k++) {
						views[k].NucleiLength = 0.0;
						views[k].NucleoliLength = 0.0;
					}
					double cx = ((double)cellsData[cellId].MaxX + (double)cellsData[cellId].MinX) / 2.0;
					double cy = ((double)cellsData[cellId].MaxY + (double)cellsData[cellId].MinY) / 2.0;
					for (int i = cellsData[cellId].QueStart; i < cellsData[cellId].QueStart + cellsData[cellId].Size; i++) {
						int position = cellsQue[i];
						int x = position % width;
						int y = position / width;
						double angle1 = atan2(cy - y, cx - x);
						double length = sqrt((cx - x) * (cx - x) + (cy - y) * (cy - y));
						for (int k = 0; k < 60; k++) {
							double value = length * sin(angle1 - k * 3.0 / 180.0 * 3.1415927);
							if (value < views[k].NucleoliLength) views[k].NucleoliLength = value;
							if (value > views[k].NucleiLength) views[k].NucleiLength = value;
						}

						deltaMap[position] = (deltaMap[position] - cellsData[cellId].MinDelta) / (cellsData[cellId].MaxDelta - cellsData[cellId].MinDelta);
					}
					double minLength = DBL_MAX;
					double maxLength = 0.0;
					for (int k = 0; k < 60; k++) {
						double length = views[k].NucleiLength - views[k].NucleoliLength;
						if (length < minLength) minLength = length;
						if (length > maxLength) maxLength = length;
					}
					cellsData[cellId].Aspect = maxLength / minLength;
				}
				if (logger.step(44)) break;
			}
			if (logger.infomation("[42] Celled", true, 42)) continue;
			if (logger.infomation("[43] Filled", true, 43)) continue;
			if (logger.infomation("[44] Adjasted", true, 44)) continue;
			if (logger.infomation("[45] Linked", true, 45)) continue;

			// Detect nucleus
			int partsRedCount = 0;
			for (int loop = 0; loop < loopCount; loop++) {
				memset(partsMap, 0, sizeof(int) * size);
				for (int queIndex = 0; queIndex < cellsQueCount; queIndex++) {
					int position = cellsQue[queIndex];
					int cellId = cellsMap[position];
					if (cellsData[cellId].IsOuter) continue;

					int partId = INT_MAX;
					bool touched = false;
					for (int j = 0; j < 4; j++) {
						int targetPosition = position + deltas[j];
						int targetId = partsMap[targetPosition];
						if (cellsMap[targetPosition] != cellId) {
							touched = true;
							continue;
						}
						if (targetId == 0) continue;
						if (partId != INT_MAX && partId != targetId) partId = -1;
						if (partId > targetId) partId = targetId;
					}

					if (partId == -1) continue;
					if (partId == INT_MAX) {
						partId = ++partsRedCount;
						partsData[partId].Position = position;
						partsData[partId].Size = 0;
						partsData[partId].IsCandidate = true;
					}
					partsMap[position] = partId;
					partsData[partId].Size++;
					if (touched) partsData[partId].IsCandidate = false;
				}

				int queCount = 0;
				for (int queIndex = 0; queIndex < cellsQueCount; queIndex++) {
					int position = cellsQue[queIndex];
					int partId = partsMap[position];
					if (partId == 0 || partsData[partId].Size <= 1 || !partsData[partId].IsCandidate) continue;
					partsQue[queCount++] = position;
				}

				for (int position = 0; position < size; position++) {
					tempMap[position] = -1;
					int cellId = cellsMap[position];
					if (cellId > 0 && !cellsData[cellId].IsOuter) {
						for (int j = 0; j < 8; j++) if (cellsMap[position + deltas[j]] != cellId) tempMap[position] = -2;
					}
				}
				if (logger.step(20)) break;

				for (int queIndex = 0; queIndex < queCount; queIndex++) {
					int position = partsQue[queIndex];
					int cellId = cellsMap[position];
					int partId = partsMap[position];
					int centerPosition = partsData[partId].Position;
					double separator = deltaMap[position];

					bool overed = false;
					for (int a = 0; a < carving; a++) {
						int tx = views[a].Cos < 0 ? -1 : 1;
						int ty = views[a].Sin < 0 ? -1 : 1;
						bool inNucleoli = true;
						double power = circler_distance;
						double remainX = 1.0, remainY = 1.0;
						int point = centerPosition;
						double miniValue = deltaMap[centerPosition];
						tempMap[point] = queIndex;

						while (true) {
							double energyX = tx * remainX / views[a].Cos;
							double energyY = ty * remainY / views[a].Sin;
							if (energyX <= energyY && energyX <= power) {
								point += tx;
								remainX = 1.0;
								remainY -= ty * energyX * views[a].Sin;
								power -= energyX;
							}
							else if (energyX > energyY && energyY <= power) {
								point += ty * width;
								remainX -= tx * energyY * views[a].Cos;
								remainY = 1.0;
								power -= energyY;
							}
							else break;

							if (tempMap[point] == -2) {
								overed = true;
								break;
							}
							if (inNucleoli && (deltaMap[point] < separator || deltaMap[point] > miniValue + 0.1)) {
								views[a].NucleoliLength = circler_distance - power;
								inNucleoli = false;
							}
							if (!inNucleoli && (deltaMap[point] < separator - 0.2 || deltaMap[point] > miniValue + thresthold_nuclei_valley)) {
								views[a].NucleiLength = circler_distance - power;
								break;
							}

							tempMap[point] = queIndex + (inNucleoli ? 0 : (queCount + 10));
							if (miniValue > deltaMap[point]) miniValue = deltaMap[point];
						}
					}
					if (overed) continue;

					double sizes[2] = { 0.0, 0.0 };
					double lines[2] = { 0.0, 0.0 };
					double circularities[2];
					for (int flag = 0; flag < 2; flag++) {
						double lx, ly, length = 0.0;
						for (int b = -moving; b <= moving; b++) length += flag == 0 ? views[ANGLE(b, carving)].NucleoliLength : views[ANGLE(b, carving)].NucleiLength;
						for (int a = 0; a <= carving; a++) {
							length -= flag == 0 ? views[ANGLE(a - (moving + 1), carving)].NucleoliLength : views[ANGLE(a - (moving + 1), carving)].NucleiLength;
							length += flag == 0 ? views[ANGLE(a + moving, carving)].NucleoliLength : views[ANGLE(a + moving, carving)].NucleiLength;
							if (flag == 0) {
								views[a].NucleoliMoving = length / (moving * 2.0 + 1.0);
							}
							else {
								views[a].NucleiMoving = length / (moving * 2.0 + 1.0);
							}
							double px = views[a].Cos * length / (moving * 2.0 + 1.0);
							double py = views[a].Sin * length / (moving * 2.0 + 1.0);
							if (a > 0) {
								sizes[flag] += abs(px * ly - lx * py) / 2.0;
								lines[flag] += sqrt((px - lx) * (px - lx) + (py - ly) * (py - ly));
							}
							lx = px;
							ly = py;
						}
						circularities[flag] = 4 * 3.1415927 * sizes[flag] / lines[flag] / lines[flag];
					}

					double colors[2] = { 0.0, 0.0 };
					int counts[2] = { 0, 0 };
					for (int j = cellsData[cellId].QueStart; j < cellsData[cellId].QueStart + cellsData[cellId].Size; j++) {
						int position = cellsQue[j];
						if (tempMap[position] == queIndex) {
							colors[0] += deltaMap[position];
							counts[0]++;
						}
						else if (tempMap[position] == queIndex + queCount + 10) {
							colors[1] += deltaMap[position];
							counts[1]++;
						}
					}
					if (counts[1] == 0) continue;

					double p1 = min(1, max(0, ((1.0 - sizes[0] / sizes[1]) - 0.3) / 0.4));
					double p2 = min(1, max(0, (circularities[0] - 0.85) / 0.15));
					double p3 = min(1, max(0, (circularities[1] - 0.6) / 0.25));
					double p4 = min(1, max(0, (sqrt(sizes[1]) + 12) / 24.0));
					double p5 = min(1, max(0, (colors[0] / counts[0] - colors[1] / counts[1] + 0.2) / 0.4));
					double p6 = min(1, max(0, (sqrt(sizes[0]) + 5.5) / 11.0));
					double hume = p1 * p2 * p3 * p4 * p5 * p6;

					if (cellsData[cellId].NucleiHume < hume) {
						cellsData[cellId].NucleiHume = hume;
						cellsData[cellId].NucleiSize = sizes[1];
						cellsData[cellId].NucleoliSize = sizes[0];
						cellsData[cellId].NucleiCenter = centerPosition;
						for (int a = 0; a <= carving; a++) {
							cellsData[cellId].NucleoliLengthes[a] = views[a].NucleoliMoving;
							cellsData[cellId].NucleiLengthes[a] = views[a].NucleiMoving;
						}
					}
				}
				if (logger.step(21)) break;
			}
			if (logger.infomation("[20] Analysed", true, 20)) continue;
			if (logger.infomation("[21] Analysed", true, 21)) continue;

			// Detect vacuoles
			int partsBlueCount = 0;
			for (int loop = 0; loop < loopCount; loop++) {
				memset(partsMap, 0, sizeof(int) * size);
				for (int queIndex = cellsQueCount - 1; queIndex >= 0; queIndex--) {
					int position = cellsQue[queIndex];
					int cellId = cellsMap[position];
					if (cellsData[cellId].IsOuter) continue;

					int partId = INT_MAX;
					bool touched = false;
					for (int j = 0; j < 4; j++) {
						int targetPosition = position + deltas[j];
						int targetId = partsMap[targetPosition];
						if (cellsMap[targetPosition] != cellId) {
							touched = true;
							continue;
						}
						if (targetId == 0) continue;
						if (partId != INT_MAX && partId != targetId) partId = -1;
						if (partId > targetId) partId = targetId;
					}

					if (partId == -1) continue;
					if (partId == INT_MAX) {
						partId = ++partsBlueCount;
						partsData[partId].Position = position;
						partsData[partId].Size = 0;
						partsData[partId].Hume = 0;
						partsData[partId].IsCandidate = true;
					}
					partsMap[position] = partId;
					partsData[partId].Size++;
					if (touched) partsData[partId].IsCandidate = false;
				}

				int queCount = 0;
				for (int queIndex = cellsQueCount - 1; queIndex >= 0; queIndex--) {
					int position = cellsQue[queIndex];
					int partId = partsMap[position];
					if (partId == 0 || partsData[partId].Size <= 1 || !partsData[partId].IsCandidate) continue;
					partsQue[queCount++] = position;
				}

				for (int position = 0; position < size; position++) {
					tempMap[position] = -1;
					int cellId = cellsMap[position];
					if (cellId > 0 && !cellsData[cellId].IsOuter) {
						for (int j = 0; j < 8; j++) if (cellsMap[position + deltas[j]] != cellId) tempMap[position] = -2;
					}
				}
				if (logger.step(30)) break;

				for (int queIndex = 0; queIndex < queCount; queIndex++) {
					int position = partsQue[queIndex];
					int cellId = cellsMap[position];
					int partId = partsMap[position];
					int centerPosition = partsData[partId].Position;
					double centerValue = deltaMap[centerPosition];
					double separator = deltaMap[position];

					bool overed = false;
					for (int a = 0; a < carving; a++) {
						int tx = views[a].Cos < 0 ? -1 : 1;
						int ty = views[a].Sin < 0 ? -1 : 1;
						double power = circler_distance;
						double remainX = 1.0, remainY = 1.0;
						int point = centerPosition;
						double maxiValue = centerValue;

						while (true) {
							double energyX = tx * remainX / views[a].Cos;
							double energyY = ty * remainY / views[a].Sin;
							if (energyX <= energyY && energyX <= power) {
								point += tx;
								remainX = 1.0;
								remainY -= ty * energyX * views[a].Sin;
								power -= energyX;
							}
							else if (energyX > energyY && energyY <= power) {
								point += ty * width;
								remainX -= tx * energyY * views[a].Cos;
								remainY = 1.0;
								power -= energyY;
							}
							else break;

							if (tempMap[point] == -2 || deltaMap[point] < centerValue) {
								overed = true;
								break;
							}
							if (deltaMap[point] > separator || deltaMap[point] < maxiValue - 0.1) {
								views[a].NucleoliLength = circler_distance - power;
								break;
							}
							if (maxiValue < deltaMap[point]) maxiValue = deltaMap[point];
						}
					}
					if (overed) continue;

					double lx, ly, length = 0.0, size = 0.0, line = 0.0;
					for (int b = -moving; b <= moving; b++) length += views[ANGLE(b, carving)].NucleoliLength;
					for (int a = 0; a <= carving; a++) {
						length -= views[ANGLE(a - (moving + 1), carving)].NucleoliLength;
						length += views[ANGLE(a + moving, carving)].NucleoliLength;
						views[a].NucleoliMoving = length / (moving * 2.0 + 1.0);
						double px = views[a].Cos * length / (moving * 2.0 + 1.0);
						double py = views[a].Sin * length / (moving * 2.0 + 1.0);
						if (a > 0) {
							size += abs(px * ly - lx * py) / 2.0;
							line += sqrt((px - lx) * (px - lx) + (py - ly) * (py - ly));
						}
						lx = px;
						ly = py;
					}
					double circularity = 4 * 3.1415927 * size / line / line;

					double x1 = width, y1 = height, x2 = 0.0, y2 = 0.0;
					for (int a = 0; a < carving; a++) {
						double x = centerPosition % width + views[a].Cos * views[a].NucleoliMoving;
						double y = centerPosition / width + views[a].Sin * views[a].NucleoliMoving;
						if (x1 > x) x1 = x;
						if (y1 > y) y1 = y;
						if (x2 < x) x2 = x;
						if (y2 < y) y2 = y;
					}
					double cx = 0.5 * (x1 + x2);
					double cy = 0.5 * (y1 + y2);

					double radius = 0;
					for (int a = 0; a < carving; a++) {
						double x = centerPosition % width + views[a].Cos * views[a].NucleoliMoving;
						double y = centerPosition / width + views[a].Sin * views[a].NucleoliMoving;
						radius += sqrt((cx - x) * (cx - x) + (cy - y) * (cy - y));
					}
					radius /= carving;

					double ringValue = 0.0;
					for (int a = 0; a < carving; a++) {
						double x = cx + views[a].Cos * radius;
						double y = cy + views[a].Sin * radius;
						ringValue += deltaMap[(int)(y + 0.5) * width + (int)(x + 0.5)];
					}
					ringValue /= carving;

					double hume = max(0, radius - 2.0) * max(0, circularity - 0.95) * max(0, (ringValue - centerValue) - 0.2);
					if (partsData[partId].Hume < hume) {
						partsData[partId].Hume = hume;
						partsData[partId].Center = centerPosition;
						partsData[partId].Area = radius * radius * 3.1415927;
						for (int a = 0; a < carving; a++) {
							partsData[partId].Lengthes[a] = views[a].NucleoliMoving;
						}
					}
				}
				if (logger.step(31)) break;

				// Exclude nuclei
				memset(tempMap, 0, sizeof(int)* size);
				for (int cellId = 1; cellId < cellsCount; cellId++) {
					Cell cell = cellsData[cellId];
					if (cell.NucleiHume < thresthold_nuclei_hume) continue;
					for (int a = 0; a <= carving; a++) {
						int tx = views[a].Cos < 0 ? -1 : 1;
						int ty = views[a].Sin < 0 ? -1 : 1;
						double power = cell.NucleiLengthes[a];
						double remainX = 1.0, remainY = 1.0;
						int point = cell.NucleiCenter;
						while (true) {
							double energyX = tx * remainX / views[a].Cos;
							double energyY = ty * remainY / views[a].Sin;
							if (energyX <= energyY) {
								point += tx;
								remainX = 1.0;
								remainY -= ty * energyX * views[a].Sin;
								power -= energyX;
							}
							else if (energyX > energyY) {
								point += ty * width;
								remainX -= tx * energyY * views[a].Cos;
								remainY = 1.0;
								power -= energyY;
							}
							tempMap[point] = 1;
							if (power < 0) break;
						}
						tempMap[cell.NucleiCenter] = 1;
					}
				}
				for (int partId = 1; partId <= partsBlueCount; partId++) if (tempMap[partsData[partId].Position] != 0) partsData[partId].IsCandidate = false;
				if (logger.step(32)) break;
			}
			if (logger.infomation("[30] Analysed", true, 30)) continue;
			if (logger.infomation("[31] Analysed", true, 31)) continue;
			if (logger.infomation("[32] Analysed", true, 32)) continue;

			// Output Data
			for (int loop = 0; loop < loopCount; loop++) {
				for (int partId = 1; partId <= partsBlueCount; partId++) {
					Part part = partsData[partId];
					if (part.Hume == 0 || !part.IsCandidate) continue;
					int cellId = cellsMap[part.Position];
					cellsData[cellId].OrganelleCount++;
					cellsData[cellId].OrganelleTotalSize += part.Area;
					if (cellsData[cellId].OrganelleMaxSize < part.Area) cellsData[cellId].OrganelleMaxSize = part.Area;
				}

				frames[resultIndex].Frame = imageIndex;
				frames[resultIndex].CellCount = 0;
				frames[resultIndex].CellSize = 0;
				frames[resultIndex].CellCircularity = 0;
				frames[resultIndex].EukaryoteCount = 0;
				frames[resultIndex].EukaryoteSize = 0;
				frames[resultIndex].EukaryoteCircularity = 0;
				frames[resultIndex].NucleiSize = 0;
				frames[resultIndex].NucleiPart = 0;
				frames[resultIndex].NucleoliSize = 0;
				frames[resultIndex].NucleoliPart = 0;
				frames[resultIndex].LinkCount = 0;
				frames[resultIndex].LinkLength = 0.0;
				frames[resultIndex].OrganelleCount = 0;
				frames[resultIndex].OrganelleTotalSize = 0.0;
				frames[resultIndex].OrganelleMaxSize = 0.0;
				frames[resultIndex].CellWithOrganelleCount = 0;

				for (int i = 0; i < wholeLinkId; i++) {
					if (links[i].Frame != resultIndex) continue;
					int id2 = links[i].To;
					int id1 = links[i].From;
					int x2 = (int)(cellsData[id2].MaxX / 2.0 + cellsData[id2].MinX / 2.0);
					int y2 = (int)(cellsData[id2].MaxY / 2.0 + cellsData[id2].MinY / 2.0);
					int x1 = (int)(cellsResultData[id1].MaxX / 2.0 + cellsResultData[id1].MinX / 2.0);
					int y1 = (int)(cellsResultData[id1].MaxY / 2.0 + cellsResultData[id1].MinY / 2.0);
					double length = sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2));
					links[i].Id = cellsData[id2].Id;
					links[i].To = y2 * width + x2;
					links[i].From = y1 * width + x1;
					cellsData[id2].LinkLength = length;
					frames[resultIndex].LinkLength += length;
					frames[resultIndex].LinkCount++;
				}

				for (int cellId = 1; cellId < cellsCount; cellId++) {
					Cell cell = cellsData[cellId];
					if (cell.IsOuter) continue;
					frames[resultIndex].CellCount++;
					frames[resultIndex].CellSize += cell.Size;
					frames[resultIndex].CellCircularity += 4 * 3.1415927 * cell.Size / cell.Line / cell.Line;
					frames[resultIndex].OrganelleCount += cell.OrganelleCount;
					frames[resultIndex].OrganelleTotalSize += cell.OrganelleTotalSize;
					frames[resultIndex].OrganelleMaxSize += cell.OrganelleMaxSize;
					if (cell.OrganelleCount > 0)frames[resultIndex].CellWithOrganelleCount++;

					if (cell.NucleiHume >= thresthold_nuclei_hume) {
						frames[resultIndex].EukaryoteCount++;
						frames[resultIndex].EukaryoteSize += cell.Size;
						frames[resultIndex].EukaryoteCircularity += 4 * 3.1415927 * cell.Size / cell.Line / cell.Line;
						frames[resultIndex].NucleiSize += cell.NucleiSize;
						frames[resultIndex].NucleiPart += cell.NucleiSize / cell.Size;
						frames[resultIndex].NucleoliSize += cell.NucleoliSize;
						frames[resultIndex].NucleoliPart += cell.NucleoliSize / cell.NucleiSize;
					}

					output1 << resultIndex << ", ";
					output1 << imageIndex << ", ";
					output1 << cell.Id << ", ";
					output1 << (cell.NucleiHume >= thresthold_nuclei_hume ? 1 : 0) << ", ";
					output1 << cell.Size << ", ";
					output1 << 4 * 3.1415927 * cell.Size / cell.Line / cell.Line << ", ";
					output1 << cell.LinkLength << ", ";
					output1 << cell.MinX << ", ";
					output1 << cell.MinY << ", ";
					output1 << cell.MaxX << ", ";
					output1 << cell.MaxY << ", ";
					output1 << cell.Aspect << ", ";
					if (cell.NucleiHume > 0) {
						output1 << cell.NucleiSize << ", ";
						output1 << cell.NucleiHume << ", ";
						output1 << cell.NucleiSize / cell.Size << ", ";
						output1 << cell.NucleoliSize << ", ";
						output1 << cell.NucleoliSize / cell.NucleiSize << ", ";
					}
					else {
						output1 << ", ";
						output1 << ", ";
						output1 << ", ";
						output1 << ", ";
						output1 << ", ";
					}
					output1 << cell.OrganelleCount << ", ";
					output1 << (cell.OrganelleCount == 0 ? 0 : (cell.OrganelleTotalSize / cell.OrganelleCount)) << ", ";
					output1 << cell.OrganelleMaxSize << ", ";
					output1 << std::endl;
				}

				output2 << resultIndex << ", ";
				output2 << imageIndex << ", ";
				output2 << frames[resultIndex].CellCount << ", ";
				output2 << frames[resultIndex].CellSize / frames[resultIndex].CellCount << ", ";
				output2 << frames[resultIndex].CellCircularity / frames[resultIndex].CellCount << ", ";
				output2 << (frames[resultIndex].LinkCount > 0 ? frames[resultIndex].LinkLength / frames[resultIndex].LinkCount : 0) << ", ";
				output2 << frames[resultIndex].EukaryoteCount << ", ";
				if (frames[resultIndex].EukaryoteCount > 0) {
					output2 << frames[resultIndex].EukaryoteSize / frames[resultIndex].EukaryoteCount << ", ";
					output2 << frames[resultIndex].EukaryoteCircularity / frames[resultIndex].EukaryoteCount << ", ";
					output2 << frames[resultIndex].NucleiSize / frames[resultIndex].EukaryoteCount << ", ";
					output2 << frames[resultIndex].NucleiPart / frames[resultIndex].EukaryoteCount << ", ";
					output2 << frames[resultIndex].NucleoliSize / frames[resultIndex].EukaryoteCount << ", ";
					output2 << frames[resultIndex].NucleoliPart / frames[resultIndex].EukaryoteCount << ", ";
				}
				else {
					output2 << ", ";
					output2 << ", ";
					output2 << ", ";
					output2 << ", ";
					output2 << ", ";
					output2 << ", ";
				}
				output2 << (double)frames[resultIndex].OrganelleCount / frames[resultIndex].CellCount << ", ";
				output2 << frames[resultIndex].OrganelleTotalSize / frames[resultIndex].OrganelleCount << ", ";
				output2 << frames[resultIndex].OrganelleMaxSize / frames[resultIndex].CellWithOrganelleCount << ", ";
				output2 << std::endl;
			}

			// Output Image
			if (outputResultImage) {
				image.clearEffect();
				image.drawGrayscaleImage(sourceImage, size * (tempCenter));
				for (int cellId = 1; cellId < cellsCount; cellId++) cellsData[cellId].IsChecked = false;
				for (int i = 0; i < size; i++) {
					int id = cellsMap[i];
					int x = i % width;
					int y = i / width;
					if (id > 0 && cellsData[id].Size > 0 && !cellsData[id].IsChecked) {
						cellsData[id].IsChecked = true;
						image.setNumber(x, max(0, y - 1), cellsData[id].Id, 0, 0, 0);
					}
				}
				for (int cellId = 1; cellId < cellsCount; cellId++) {
					Cell cell = cellsData[cellId];
					if (cell.NucleiHume <= 0) continue;
					for (int a = 0; a <= carving; a++) {
						double px = cell.NucleiCenter % width + views[a].Cos * cell.NucleiLengthes[a];
						double py = cell.NucleiCenter / width + views[a].Sin * cell.NucleiLengthes[a];
						if (cell.NucleiHume >= thresthold_nuclei_hume) {
							image.effectPixel((int)((px + 0.5) * imageZoom), (int)((py + 0.5) * imageZoom), 0, 0, 255);
						}
					}
				}
				for (int partId = 1; partId <= partsBlueCount; partId++) {
					if (partsData[partId].Hume == 0 || !partsData[partId].IsCandidate) continue;
					for (int a = 0; a < carving; a++) {
						double px = partsData[partId].Center % width + views[a].Cos * partsData[partId].Lengthes[a];
						double py = partsData[partId].Center / width + views[a].Sin * partsData[partId].Lengthes[a];
						image.effectPixel((int)((px + 0.5) * imageZoom), (int)((py + 0.5) * imageZoom), 0, 255, 0);
					}
				}
				for (int i = 0; i < size; i++) {
					int id = cellsMap[i];
					int x = i % width, y = i / width;
					if (id > 0) image.effectGrid(x, y, x > 0 && id != cellsMap[y * width + x - 1], y > 0 && id != cellsMap[(y - 1) * width + x], x < width - 1 && id != cellsMap[y * width + x + 1], y < height - 1 && id != cellsMap[(y + 1) * width + x], 255, cellsData[id].IsOuter ? 255 : 0, cellsData[id].IsOuter ? 255 : 0);
				}
				for (int i = 0; i < wholeLinkId; i++) {
					if (links[i].Frame >= resultIndex - 4) image.effectLine(links[i].From, links[i].To, 0, 0, 255);
				}
				sprintf_s(filename, "%s\\result_%04d.bmp", outputPath, resultIndex);
				if (image.save(filename))logger.warning("Image Save Failed.");
				if (logger.infomation("Image Saved", false)) continue;
			}

			// Data for next frame
			int* tempCellsMap = cellsResultMap;
			cellsResultMap = cellsMap;
			cellsMap = tempCellsMap;

			Cell* tempCellsData = cellsResultData;
			cellsResultData = cellsData;
			cellsData = tempCellsData;

			cellsResultCount = cellsCount;

			// Finish
			sprintf_s(filename, "Frame %d Finished", imageIndex);
			logger.restart(filename);
		}
	}
	else logger.error("Memory Allocation Failed.");

	// Release
	image.close();
	logger.finish();
	output1.close();
	output2.close();
	free(sourceImage);
	free(deltaMap);
	free(sortMap);
	free(tempMap);
	free(taskData);
	free(views);
	free(frames);
	free(links);
	free(cellsMap);
	free(cellsQue);
	free(cellsData);
	free(cellsResultMap);
	free(cellsResultData);
	free(partsMap);
	free(partsQue);
	free(partsData);

	return 0;
}