gimage.cpp

一个由Mike Gashler完成的机器学习方面的includes neural net, naive bayesian classifier, decision tree, KNN, a genet

			pixOld = m_pPixels[dst + xx];
			m_pPixels[dst + xx] = gRGB((a * gRed(pix) + (256 - a) * gRed(pixOld)) >> 8,
										(a * gGreen(pix) + (256 - a) * gGreen(pixOld)) >> 8,
										(a * gBlue(pix) + (256 - a) * gBlue(pixOld)) >> 8);
		}
		dst += m_nWidth;
		src += pSource->m_nWidth;
	}
}

void GImage::MoveLight(double dRadians, float fAmount)
{
	float dx = (float)cos(dRadians);
	float dy = (float)sin(dRadians);
	int nDelta;
	GColor c1, c2;
	int x, y;
	for(y = 0; y < m_nHeight; y++)
	{
		for(x = 0; x < m_nWidth; x++)
		{
			c1 = GetPixel(x, y);
			c2 = InterpolatePixel(x + dx, y + dy);
			nDelta = (int)gGray(c1) - gGray(c2);
			SetPixel(x, y, gARGB(
				gAlpha(c1),
				ClipChan(gRed(c1) + (int)(nDelta * fAmount)),
				ClipChan(gGreen(c1) + (int)(nDelta * fAmount)),
				ClipChan(gBlue(c1) + (int)(nDelta * fAmount))));
		}
	}
}

GImage* GImage::Munge(int nStyle, float fExtent)
{
	GImage* pMunged = new GImage();
	pMunged->SetSize(m_nWidth, m_nHeight);
	int x, y;
	GColor col;
	double d;
	switch(nStyle)
	{
		case 0: // particle-blur (pick a random pixel in the proximity)
			for(y = 0; y < m_nHeight; y++)
			{
				for(x = 0; x < m_nWidth; x++)
				{
					col = SafeGetPixel(
										(int)(x + fExtent * m_nWidth * (GBits::GetRandomDouble() - .5)),
										(int)(y + fExtent * m_nHeight * (GBits::GetRandomDouble() - .5))
									);
					pMunged->SetPixel(x, y, col);
				}
			}
			break;

		case 1: // shadow threshold (throw out all pixels below a certain percent of the total brighness)
			{
				fExtent = 1 - fExtent;
				fExtent *= fExtent;
				fExtent *= fExtent;
				fExtent = 1 - fExtent;

				// Create the histogram data
				unsigned int pnHistData[257];
				memset(pnHistData, '\0', sizeof(int) * 257);
				unsigned int nSize = m_nWidth * m_nHeight;
				unsigned int nPos;
				unsigned int nGray;
				for(nPos = 0; nPos < nSize; nPos++)
					pnHistData[gGray(m_pPixels[nPos]) >> 8]++;

				// Turn it into cumulative histogram data
				int n;
				for(n = 1; n < 256; n++)
					pnHistData[n] += pnHistData[n - 1];

				// Find the cut-off
				unsigned int nCutOff = (unsigned int)(fExtent * pnHistData[255]);
				for(n = 0; n < 256 && pnHistData[n] < nCutOff; n++)
				{
				}

				// Copy all the data above the threshold
				for(y = 0; y < m_nHeight; y++)
				{
					for(x = 0; x < m_nWidth; x++)
					{
						col = GetPixel(x, y);
						nGray = gGray(col) >> 8;
						if(nGray > (unsigned int)n)
							pMunged->SetPixel(x, y, col);
					}
				}
			}
			break;

		case 2: // waves
			for(y = 0; y < m_nHeight; y++)
			{
				for(x = 0; x < m_nWidth; x++)
				{
					col = SafeGetPixel(
										(int)(x + fExtent * (m_nWidth / 2) * cos((double)x * 16 / m_nWidth)),
										(int)(y + fExtent * (m_nHeight / 2) * sin((double)y * 16 / m_nHeight))
									);
					pMunged->SetPixel(x, y, col);
				}
			}
			break;

		case 3: // waved in or out of the middle
				fExtent = 1 - fExtent;
				fExtent *= fExtent;
				fExtent = 1 - fExtent;
			for(y = 0; y < m_nHeight; y++)
			{
				for(x = 0; x < m_nWidth; x++)
				{
					d = atan2((double)(y - m_nHeight / 2), (double)(x - m_nWidth / 2));
					d = fExtent * cos(d * 5);
					col = SafeGetPixel(
										(int)((1 - d) * x + d * (m_nWidth / 2)),
										(int)((1 - d) * y + d * (m_nHeight / 2))
									);
					pMunged->SetPixel(x, y, col);
				}
			}
			break;
	}
	return pMunged;
}

void GImage::FillTriangle(int x1, int y1, int x2, int y2, int x3, int y3, GColor c)
{
	// Get y1 on top
	int t;
	if(y2 < y1)
	{
		t = y2;
		y2 = y1;
		y1 = t;
		t = x2;
		x2 = x1;
		x1 = t;
	}
	if(y3 < y1)
	{
		t = y3;
		y3 = y1;
		y1 = t;
		t = x3;
		x3 = x1;
		x1 = t;
	}

	// Get y2 in the middle
	if(y3 < y2)
	{
		t = y3;
		y3 = y2;
		y2 = t;
		t = x3;
		x3 = x2;
		x2 = t;
	}

	// Compute step sizes
	float fx1 = (float).5 + x1;
	float fx2 = (float).5 + x1;
	float dx1, dx2;
	if(y1 == y2)
	{
		fx1 = (float).5 + x2;
		dx1 = 0;
	}
	else
		dx1 = (float)(x2 - x1) / (y2 - y1);
	if(y1 == y3)
	{
		fx2 = (float).5 + x3;
		dx2 = 0;
	}
	else
		dx2 = (float)(x3 - x1) / (y3 - y1);

	// Draw the first half
	int x, xMax, y;
	for(y = y1; y <= y2; y++)
	{
		if(fx1 < fx2)
		{
			x = (int)fx1;
			xMax = (int)fx2;
		}
		else
		{
			x = (int)fx2;
			xMax = (int)fx1;
		}
		for( ; x <= xMax; x++)
			SetPixel(x, y, c);
		fx1 += dx1;
		fx2 += dx2;
	}

	// Draw the second half
	fx1 = (float).5 + x2;
	if(y2 == y3)
		dx1 = 0;
	else
		dx1 = (float)(x3 - x2) / (y3 - y2);
	fx1 += dx1;
	for( ; y <= y3; y++)
	{
		if(fx1 < fx2)
		{
			x = (int)fx1;
			xMax = (int)fx2;
		}
		else
		{
			x = (int)fx2;
			xMax = (int)fx1;
		}
		for( ; x <= xMax; x++)
			SetPixel(x, y, c);
		fx1 += dx1;
		fx2 += dx2;
	}
}

void GImage::Stretch(int nXStart, int nYStart, int nXEnd, int nYEnd)
{
	GImage tmp;
	tmp.SetSize(GetWidth(), GetHeight());
	GColor col;
	float d;
	double dSize = sqrt((double)((nXEnd - nXStart) * (nXEnd - nXStart) + (nYEnd - nYStart) * (nYEnd - nYStart)));
	if(dSize == 0)
		dSize = .01;
	int x, y;
	for(y = 0; y < (int)tmp.GetHeight(); y++)
	{
		for(x = 0; x < (int)tmp.GetWidth(); x++)
		{
			d = (float)(exp(-(((x - nXEnd) * (x - nXEnd) + (y - nYEnd) * (y - nYEnd)) / (dSize * dSize))));
			col = InterpolatePixel((float)x - d * (nXEnd - nXStart), (float)y - d * (nYEnd - nYStart));
			tmp.SetPixel(x, y, col);
		}
	}
	SwapData(&tmp);
}

void GImage::MakeCaptcha(const char* szText)
{
	int nWidth = MeasureHardTextWidth(32, szText, 1);
	SetSize(nWidth + 32, 64);
	Clear(0xffeeeecc);
	GRect r(16, 16, nWidth + 16, 32);
	DrawHardText(&r, szText, 0xff444400, 1);
	GImage* pTmp = Munge(0, (float).01);
	SwapData(pTmp);
	delete(pTmp);
	int i;
	for(i = 0; i < 3; i++)
		DrawLine(rand() % (nWidth + 32), rand() % 64, rand() % (nWidth + 32), rand() % 64, 0xff444400);
	pTmp = Munge(3, (float).05);
	SwapData(pTmp);
	delete(pTmp);
	for(i = 0; i < 3; i++)
	{
		int x = rand() % (nWidth + 32);
		int y = rand() % 64;
		int dx = rand() % 20 - 10;
		int dy = rand() % 20 - 10;
		Stretch(x, y, x + dx, y + dy);
	}
	for(i = 0; i < 2; i++)
		DrawLine(rand() % (nWidth + 32), rand() % 64, rand() % (nWidth + 32), rand() % 64, 0xff444400);
	pTmp = Munge(0, (float).01);
	SwapData(pTmp);
	delete(pTmp);
}

void GImage::MakeGaussianKernel(int nWidth, float fDepth)
{
	GAssert(nWidth >= 1, "out of range");
	GAssert(fDepth >= 1 && fDepth <= 255, "out of range");
	double dRadius = sqrt(MAX(0.0, -2.0 * log(1.0 / (double)fDepth)));
	SetSize(nWidth, nWidth);
	double dCenter = (double)(nWidth - 1) / 2.0;
	int x, y, v;
	double r;
	float val;
	for(y = 0; y < nWidth; y++)
	{
		for(x = 0; x < nWidth; x++)
		{
			r = sqrt((dCenter - x) * (dCenter - x) + (dCenter - y) * (dCenter - y));
			val = fDepth * (float)GMath::gaussian(r * dRadius / dCenter);
			v = ClipChan((int)(val + (float).5));
			SetPixel(x, y, gARGB(0xff, v, v, v));
		}
	}
}

void GImage::HighPassFilter(double dExtent)
{
	dExtent *= (2.0 * MAX(GetWidth(), GetHeight()));

	// Convert to the Fourier domain
	int nChannelWidth, nChannelHeight;
	struct ComplexNumber* pArray = GFourier::ImageToFFTArray(this, &nChannelWidth, &nChannelHeight);
	ArrayHolder<struct ComplexNumber*> hArray(pArray);
	int nChannelSize = nChannelWidth * nChannelHeight;

	// Filter out the low frequency data
	double d, fac;
	int radius = (int)dExtent;
	int x, y, i, nStart;
	for(y = 0; y < radius && y < nChannelHeight; y++)
	{
		for(x = 0; x < radius && x < nChannelWidth; x++)
		{
			d = sqrt((double)(x * x + y * y));
			if((int)d >= radius)
				break;
			fac = GMath::smoothedIdentity(d / dExtent, 2.0);
			nStart = 0;
			for(i = 0; i < 3; i++)
			{
				pArray[nStart + y * nChannelWidth + x].real *= fac;
				pArray[nStart + y * nChannelWidth + x].imag *= fac;
				pArray[nStart + y * nChannelWidth + nChannelWidth - 1 - x].real *= fac;
				pArray[nStart + y * nChannelWidth + nChannelWidth - 1 - x].imag *= fac;
				pArray[nStart + (nChannelHeight - 1 - y) * nChannelWidth + x].real *= fac;
				pArray[nStart + (nChannelHeight - 1 - y) * nChannelWidth + x].imag *= fac;
				pArray[nStart + (nChannelHeight - 1 - y) * nChannelWidth + nChannelWidth - 1 - x].real *= fac;
				pArray[nStart + (nChannelHeight - 1 - y) * nChannelWidth + nChannelWidth - 1 - x].imag *= fac;
				nStart += nChannelSize;
			}
		}
	}

	// Convert back to the spatial domain
	GFourier::FFTArrayToImage(pArray, nChannelWidth, nChannelHeight, this, true);
}

class GImageThresholdFinder : public GRealVectorCritic
{
protected:
	unsigned int* m_pHistData;

public:
	GImageThresholdFinder(unsigned int* pHistData)
		: GRealVectorCritic(6)
	{
		m_pHistData = pHistData;
	}

	virtual ~GImageThresholdFinder()
	{
	}

protected:
	virtual double ComputeError(double* pVector)
	{
		double dSumSquaredError = 0;
		double d, error;
		int i;
		for(i = 0; i < 1020; i++)
		{
			d = pVector[0] * GMath::gaussian(((double)i - pVector[1]) / pVector[2]) +
				pVector[3] * GMath::gaussian(((double)i - pVector[4]) / pVector[5]);
			error = m_pHistData[i] - d;
			dSumSquaredError += (error * error);
		}
		return dSumSquaredError;
	}
};
/*
// todo: remove this function
void TestThreshold(unsigned int* pHist, double* pVector)
{
	unsigned int max = 0;
	int i;
	for(i = 0; i < 1020; i++)
		max = MAX(max, pHist[i]);
	GImage image;
	image.SetSize(1020, 1020);

	// Draw histogram
	for(i = 0; i < 1020; i++)
		image.SetPixel(i, pHist[i] * 1019 / max, 0xffff0000);

	// Draw gaussians
	for(i = 0; i < 1020; i++)
	{
		image.SetPixel(i, (int)(pVector[0] * GMath::gaussian(((double)i - pVector[1]) / pVector[2]) * 1019 / max), 0xff00ff00);
		image.SetPixel(i, (int)(pVector[3] * GMath::gaussian(((double)i - pVector[4]) / pVector[5]) * 1019 / max), 0xff0000ff);
	}
	if(!image.SaveBMPFile("c:\\mike\\test.bmp"))
	{
		GAssert(false, "failed to save bmp");
	}
}
*/
int GImage::ComputeOptimalGrayscaleThreshold()
{
	// Create the histogram data
	unsigned int pnHistData[1021];
	memset(pnHistData, '\0', sizeof(unsigned int) * 1021);
	int nSize = m_nWidth * m_nHeight;
	int nPos;
	int nGray;
	for(nPos = 0; nPos < nSize; nPos++)
	{
		nGray = gGray(m_pPixels[nPos]) >> 6;
		GAssert(nGray >= 0 && nGray < 1021, "out of range");
		pnHistData[nGray]++;
	}
	pnHistData[1019] += pnHistData[1020];

	// Fit two gaussians to the histogram
	GImageThresholdFinder critic(pnHistData);
	GMomentumGreedySearch searcher(&critic);
	double initialState[6];
	initialState[0] = (double)(GetWidth() * GetHeight()) / 128;
	initialState[1] = 340;
	initialState[2] = 200;
	initialState[3] = initialState[0];
	initialState[4] = 680;
	initialState[5] = initialState[2];
	searcher.SetState(initialState);
	searcher.SetAllStepSizes(1);

	// Do the search
	int i;
	for(i = 0; i < 5000; i++)
		searcher.Iterate();

	// Use quadratic formula to solve for intersection of gaussians
	double* vec = critic.GetBestYet();
//TestThreshold(pnHistData, vec);
	double vec2squared = vec[2] * vec[2];
	double vec5squared = vec[5] * vec[5];
	if(vec5squared == vec2squared)
		vec5squared -= 1e-9;
	double a = vec2squared - vec5squared;
	double b = 2.0 * (vec5squared * vec[1] - vec2squared * vec[4]);
	double c = vec2squared * vec[4] * vec[4] -
				vec5squared * vec[1] * vec[1] +
				2.0 * vec2squared * vec5squared * (log(vec[0]) - log(vec[3]));
	double discriminant = b * b - 4.0 * a * c;
	if(discriminant < 0)
	{
		//GAssert(false, "The gaussians don't intersect");
		return 32640;
	}
	double d = sqrt(discriminant);
	double answer1 = (-b + d) / (2.0 * a);
	double answer2 = (-b - d) / (2.0 * a);
	if(answer1 > vec[1] && answer1 < vec[4])
		return (int)(answer1 * 64 + .5);
	else if(answer1 > vec[4] && answer1 < vec[1])
		return (int)(answer1 * 64 + .5);
	else if(answer2 > vec[1] && answer2 < vec[4])
		return (int)(answer2 * 64 + .5);
	else if(answer2 > vec[4] && answer2 < vec[1])
		return (int)(answer2 * 64 + .5);
	else
	{
		//GAssert(false, "Neither intersection is between the two peaks");
		return (int)((vec[1] + vec[4]) * 32.0 + .5);
	}
}

void GImage::Threshold(int nGrayscaleValue)
{
	int x, y, g;
	for(y = 0; y < (int)m_nHeight; y++)
	{
		for(x = 0; x < (int)m_nWidth; x++)
		{
			g = gGray(GetPixel(x, y));
			if(g >= nGrayscaleValue)
				SetPixel(x, y, 0xffffffff);
			else
				SetPixel(x, y, 0xff000000);
		}
	}
}

void GImage::MedianFilter(float fRadius)
{
	int x, y, i, j, d2, median, size;
	int r = (int)(fRadius + .9999999);
	GIntArray array(r * r);
	float fRadiusSquared = fRadius * fRadius;
	for(y = 0; y < GetHeight(); y++)
	{
		for(x = 0; x < GetWidth(); x++)
		{
			array.Clear();
			for(j = y - r; j <= y + r; j++)
			{
				if(j < 0)
					continue;
				if(j >= GetHeight())
					break;
				for(i = x - r; i <= x + r; i++)
				{
					if(i < 0)
						continue;
					if(i >= GetWidth())
						break;
					d2 = (x - i) * (x - i) + (y - j) * (y - j);
					if(d2 > (int)fRadiusSquared)
						continue;
					array.AddInt(gGray(GetPixel(i, j)));
				}
			}
			size = array.GetSize();
			GAssert(size > 0, "no data");
			array.SortByValue();
			median = (array.GetInt((size - 1) / 2) + array.GetInt(size / 2)) / 2;
			median = median >> 8;
			GAssert(median >= 0 && median <= 255, "out of range");
			SetPixel(x, y, gARGB(0xff, median, median, median));
		}
	}
}

void GImage::Dialate(GImage* pStructuringElement)
{
	GImage target;
	target.SetSize(GetWidth(), GetHeight());
	int x, y, i, j, r, g, b, u, v;
	GColor c1, c2;
	int dx = pStructuringElement->GetWidth() / 2;
	int dy = pStructuringElement->GetHeight() / 2;
	for(y = 0; y < GetHeight(); y++)
	{
		for(x = 0; x < GetWidth(); x++)
		{
			r = 0;
			g = 0;
			b = 0;
			for(j = 0; j < (int)pStructuringElement->GetHeight(); j++)
			{
				v = y + j -