www.gusucode.com > VC++遥感原理与数字摄影测量实习源码程序 > VC++遥感原理与数字摄影测量实习源码程序\code\Moravec\Moravec.cpp

    // Moravec.cpp : Defines the entry point for the console application.
// Download by http://www.NewXing.com

#include "stdafx.h"
#include "Moravec.h"
#include <vector>
#include <math.h>

#ifdef _DEBUG
#define new DEBUG_NEW
#undef THIS_FILE
static char THIS_FILE[] = __FILE__;
#endif

#define FLT_MAX         3.402823466e+38F        /* max value */

/////////////////////////////////////////////////////////////////////////////
// The one and only application object

CWinApp theApp;

using namespace std;

struct FEATUREPOINT
{
	int x;
	int y;		//坐标
	float IV;	//兴趣值
};

vector<FEATUREPOINT> RemoveReplicative(vector<FEATUREPOINT> v)
{
	vector<FEATUREPOINT> ret;

	vector<FEATUREPOINT>::iterator iter = v.begin();

	ret.clear();

	ret.push_back(*iter);

	vector<FEATUREPOINT>::iterator iter2;

	int i(0);

	BOOL bEXIST(FALSE);

	for (iter = v.begin(); iter != v.end(); ++iter)
	{
		bEXIST = false;

		for (iter2 = ret.begin(); iter2 != ret.end(); ++iter2)
		{
			int a = iter->x;
			int aa = iter2->x;
			int b = iter->y;
			int bb = iter2->y;

			if (a == aa && b == bb)
			{
				bEXIST = TRUE;	//存在
				break;
			}
		}

		if (bEXIST == FALSE)
		{
			ret.push_back(*iter);
			++i;
		}
	}

	cout << "i = " << i << endl;

	return ret;
}

void SaveBand(int width, int height, int byteCount, int biBitCount, LPBYTE pBits, CString SavePath)
{
	BITMAPFILEHEADER bmfh;
	BITMAPINFOHEADER bmih;
	
	bmfh.bfType = 0x4d42;        // 0x42 = "B" 0x4d = "M"	
	bmfh.bfReserved1 = 0;
	bmfh.bfReserved2 = 0;
	
	if (biBitCount == 8)
	{
		bmfh.bfOffBits = sizeof(BITMAPFILEHEADER) + sizeof(BITMAPINFOHEADER) + 256*4;
		bmfh.bfSize = sizeof(BITMAPFILEHEADER) + sizeof(BITMAPINFOHEADER) + 256*4 + byteCount;
	}
	else if (biBitCount == 24)
	{
		bmfh.bfOffBits = sizeof(BITMAPFILEHEADER) + sizeof(BITMAPINFOHEADER);
		bmfh.bfSize = sizeof(BITMAPFILEHEADER) + sizeof(BITMAPINFOHEADER) + byteCount;
	}
	else
	{
		MessageBox(NULL, "8bit or 24bit accepted", "错误", MB_ICONERROR);
		return;
	}
	
	
	bmih.biBitCount = biBitCount;
	bmih.biWidth = width;
	bmih.biHeight = height;
	bmih.biSize = 40;
	bmih.biPlanes = 1;
	bmih.biCompression = BI_RGB;
	bmih.biSizeImage = bmfh.bfSize;
	bmih.biClrImportant = 0;
	bmih.biClrUsed = 0;
	bmih.biXPelsPerMeter = 0;
	bmih.biYPelsPerMeter = 0;
	
	CFile f;
	
	//保存位图
	if (f.Open(SavePath, CFile::modeCreate | CFile::modeWrite))
	{
		f.Write(&bmfh, sizeof(bmfh));
		
		f.Write(&bmih, sizeof(bmih));
		
		if (biBitCount == 8)
		{
			RGBQUAD rgb[256];
			int i(0);
			for (; i < 256; ++i)
			{
				rgb[i].rgbBlue = i;
				rgb[i].rgbGreen = i;
				rgb[i].rgbRed = i;
				rgb[i].rgbReserved = 0;
			}

			//修改颜色表255对应的颜色
			rgb[255].rgbBlue = 0;
			rgb[255].rgbGreen = 0;
			rgb[255].rgbRed = 255;
			f.Write(rgb, sizeof(RGBQUAD)*256);
		}
		
		f.Write(pBits, byteCount);
		
		f.Close();
	}
	else
	{
		MessageBox(NULL, "保存图像失败", "错误", MB_ICONERROR);
		return;
	}
}

int _tmain(int argc, TCHAR* argv[], TCHAR* envp[])
{
	if (argc != 9)
	{
		::MessageBox(NULL,"启动参数不正确,请从GUI.exe或命令行参数启动","",MB_ICONASTERISK);
		return -1;
	}
	system("Mode con: COLS=50 LINES=15");
	system("title Moravec_Operator");
	int nRetCode = 0;

	// initialize MFC and print and error on failure
	if (!AfxWinInit(::GetModuleHandle(NULL), NULL, ::GetCommandLine(), 0))
	{
		// TODO: change error code to suit your needs
		cerr << _T("Fatal Error: MFC initialization failed") << endl;
		nRetCode = 1;
		return nRetCode;
	}

	CString srcFile,dstFile,txtFile,threshold,window1,window2,flag1,flag2;	//启动参数
	srcFile = argv[1];
	dstFile = argv[2];
	txtFile = argv[3];
	threshold = argv[4];	//阈值
	window1 = argv[5];
	window2 = argv[6];
	flag1 = argv[7];		//判断是否在创建结果文件后显示
	flag2 = argv[8];

	int m_threshold = atoi(threshold);
	int m_window1 = atoi(window1);
	int m_window2 = atoi(window2);
	int m_flag1 = atoi(flag1);
	int m_flag2 = atoi(flag2);

	if (GetFileAttributes(srcFile) == -1)
	{
		//the module does not exist
		CString err;
		err.Format("文件:%s不存在", srcFile);		
		::MessageBox(NULL,err, "Error", MB_ICONERROR);
		return -1;
	}

	FILE* pSrcFile = NULL;
	pSrcFile = fopen(srcFile,"rb");
	if(pSrcFile == NULL)
	{
		printf("原始图像打开失败\n");
		nRetCode = 1;
		return nRetCode;
	}

	BITMAPFILEHEADER bmfh;
	BITMAPINFOHEADER bmih;

	fread(&bmfh, sizeof(BITMAPFILEHEADER), 1, pSrcFile);
	fread(&bmih, sizeof(BITMAPINFOHEADER), 1, pSrcFile);

	if(bmih.biBitCount != 8)
	{
		printf("只支持8位位图\n");
		nRetCode = 1;
		return nRetCode;
	}

	DWORD width = bmih.biWidth;
	DWORD height = bmih.biHeight;
	DWORD _width = (width*8 + 31)/32*4;	//saved width in files
	DWORD pixelCount = _width*height;

	fseek(pSrcFile, 256*sizeof(RGBQUAD), SEEK_CUR);

	BYTE* pSrcBits = new BYTE[pixelCount];

	fread(pSrcBits,sizeof(BYTE),pixelCount,pSrcFile);

	fclose(pSrcFile);

	//开始计算
	int i(0), j(0), m(0), n(0);
	
	int k = m_window1 / 2;

	float v1(0), v2(0), v3(0), v4(0);
	float IV(0);

	float* pInterest = new float[width*height];
	ZeroMemory(pInterest, sizeof(float)*width*height);

	int Count1(0);

	for (i = k; i < height - k; ++i)
	{
		for (j = k; j < width - k; ++j)
		{
			//计算一个点
			v1 = v2 = v3 = v4 = 0.0;
			for (m = 0; m < m_window1-1; ++m)
			{
				v1 += (pSrcBits[(i-k+m)*_width + (j)] - pSrcBits[(i-k+1+m)*_width + (j)])*
					  (pSrcBits[(i-k+m)*_width + (j)] - pSrcBits[(i-k+1+m)*_width + (j)]);

				v2 += (pSrcBits[(i-k+m)*_width + (j-k+m)] - pSrcBits[(i-k+1+m)*_width + (j-k+m+1)])*
					  (pSrcBits[(i-k+m)*_width + (j-k+m)] - pSrcBits[(i-k+1+m)*_width + (j-k+m+1)]);

				v3 += (pSrcBits[(i)*_width + (j-k+m)] - pSrcBits[(i)*_width + (j-k+1+m)])*
					  (pSrcBits[(i)*_width + (j-k+m)] - pSrcBits[(i)*_width + (j-k+1+m)]);

				v4 += (pSrcBits[(i-k+m)*_width + (j+k-m)] - pSrcBits[(i-k+m+1)*_width + (j+k-m-1)])*
					  (pSrcBits[(i-k+m)*_width + (j+k-m)] - pSrcBits[(i-k+m+1)*_width + (j+k-m-1)]);
			}

			IV = min(v1, v2);
			IV = min(IV, v3);
			IV = min(IV, v4);

			if (IV > m_threshold)
			{
				//成为候选点
				pInterest[i*width+j] = IV;
				++Count1;
			}
		}
	}
	
//	cout << "候选点数: " << Count1 << endl;

	//候选点列表
/*	FILE* tempFile1 = fopen("d:\\temp1.txt", "wt");
	
	if (!tempFile1)
		return -1;
	
	for (i = 0; i < height; ++i)
	{
		for (j = 0; j < width; ++j)
		{
			if (pInterest[i*width+height] == 0)
				continue;

			fprintf(tempFile1, "%d\t%d\t%.20f\n", i, j, pInterest[i*width+j]);
		}		
	}
	fclose(tempFile1);*/
	//候选点列表

	bool* bMatrix = new bool[width*height];
	memset(bMatrix, 1, sizeof(bool)*width*height);
//	ZeroMemory(bMatrix, sizeof(bool)*width*height);

	vector<FEATUREPOINT> FeaturePoint;
	FEATUREPOINT tempFP;
	ZeroMemory(&tempFP, sizeof(FEATUREPOINT));

	int Count2(0);

	int k2 = m_window2 / 2;
	
	float max2(0);
	int tempX(0), tempY(0);
	bool b(false);
	
	for(i = k2; i < height - k2; i += m_window2)
	{
		for(j = k2; j < width - k2; j += m_window2)
		{
			max2 = 0;
			b = false;
			for (int m = i-k2; m < i+k2+1; ++m)
			{
				for (int n = j-k2; n<j+k2+1; ++n)
				{
					if(pInterest[m*width + n]>max2)
					{
						max2 = pInterest[m*width + n];
						tempY = m;
						tempX = n;
						b = true;
					}
				}
			}
			if (b)
			{
				bMatrix[tempY*width+tempX] = 1;

				tempFP.IV = pInterest[tempY*width+tempX];
				tempFP.x = tempX;
				tempFP.y = tempY;
				FeaturePoint.push_back(tempFP);

				++Count2;	//尚未剔除重复点
			}
		}
	}

	cout << "未剔除重复: " << FeaturePoint.size() << endl;


	vector<FEATUREPOINT> FeaturePoint2;
	FeaturePoint2 = RemoveReplicative(FeaturePoint);
//	cout << "特征点数: " << Count2 << endl;

	cout << "剔除重复: " << FeaturePoint2.size() << endl;

	//保存结果位图

	vector<FEATUREPOINT>::iterator iter;

	for (iter = FeaturePoint2.begin(); iter != FeaturePoint2.end(); ++iter)
	{
		int i = iter->y;
		int j = iter->x;

		pSrcBits[i*_width+j] = 255;
		
		pSrcBits[i*_width+j+1] = 255;
		pSrcBits[i*_width+j-1] = 255;
		pSrcBits[(i+1)*_width+j] = 255;
		pSrcBits[(i-1)*_width+j] = 255;
		
		pSrcBits[i*_width+j+2] = 255;
		pSrcBits[i*_width+j-2] = 255;
		pSrcBits[(i+2)*_width+j] = 255;
		pSrcBits[(i-2)*_width+j] = 255;
	}
/*	for (i = 0; i < height; ++i)
	{
		for (j = 0; j < width; ++j)
		{
			if (bMatrix[i*width+j])
			{
				pSrcBits[i*_width+j] = 255;

				pSrcBits[i*_width+j+1] = 255;
				pSrcBits[i*_width+j-1] = 255;
				pSrcBits[(i+1)*_width+j] = 255;
				pSrcBits[(i-1)*_width+j] = 255;

				pSrcBits[i*_width+j+2] = 255;
				pSrcBits[i*_width+j-2] = 255;
				pSrcBits[(i+2)*_width+j] = 255;
				pSrcBits[(i-2)*_width+j] = 255;
			}
		}
	}*/

	SaveBand(width, height, _width*height, 8, pSrcBits, dstFile);

	//保存特征点列表
	FILE*  pTxtFile= fopen(txtFile, "wt");
	
	if (!pTxtFile)
		return -1;

	fprintf(pTxtFile, "#\tMoravec算子\n#\t源图像: %s\n#\t结果图: %s\n", srcFile, dstFile);
	fprintf(pTxtFile, "#\t阈值: %d\n", m_threshold);
	fprintf(pTxtFile, "#\t窗口1: %d×%d\n#\t窗口2: %d×%d\n", m_window1, m_window1, m_window2, m_window2);
	fprintf(pTxtFile, "#\t兴趣点数: %d\n", Count1);
	fprintf(pTxtFile, "#\t特征点数: %d\n\n", FeaturePoint2.size());
	fprintf(pTxtFile, "坐标X\t坐标Y\t兴趣值\n");

//	vector<FEATUREPOINT>::iterator iter;

	for (iter = FeaturePoint2.begin(); iter != FeaturePoint2.end(); ++ iter)
		fprintf(pTxtFile, "%d\t%d\t%.1f\n", iter->x, iter->y, iter->IV);

	fclose(pTxtFile);
	
	//clean up
	delete[] pInterest;
	delete[] bMatrix;
	delete[] pSrcBits;

	pInterest = NULL;
	bMatrix = NULL;
	pSrcBits = NULL;

	if (m_flag1 == 1)
	{
		ShellExecute(NULL, "open", dstFile, NULL, NULL, SW_SHOWNORMAL);
	}

	if (m_flag2 == 1)
	{
		ShellExecute(NULL, "open", txtFile, NULL, NULL, SW_SHOWNORMAL);
	}

	return nRetCode;
}