Image Alignment and Feature Locating Using Shape Contexts

Description:

This deliverable explores the feasibility of using the idea of shape contexts to provide a robust means of quantifying and representing the shape of an object. This representation is then used along with several other techniques to support image alignment and feature location of an image. The idea is to manually select features on a template image, and then align a test image to the template image thereby locating the features in question on the test image.

Shape context is characterized by the spatial relationship between a point and all other points in a shape. Specifically, it is defined by a set of points sampled from the internal or external edges of a shape. The edges may be obtained by using an edge detection algorithm of which there are several to choose from, depending on the application. For example, a simpler edge detector applies convolution with Prewitt or Sobel gradient operators to filter through differentiation and neighborhood averaging. A more advanced edge detector such as the Boie-Cox uses efficient filters to perform both detection and segmentation followed by noise filtering. This application uses the Boie-Cox edge detector. The picture below shows an example of an edge-detected face.

100 Dollar Bill

The sampling of points may be random, which does not guarantee uniform spacing between sampled points, which is preferred, although not critical. Then for a given sampled point, a shape context descriptor is defined by determining the set of vectors from this point to all other sampled points on the shape. Specifically, the shape context for a point is a log-polar histogram that sorts all vectors for a given point by relative distance and angular orientation. Therefore, for N sampled points, each point will have N - 1 vectors from that point to all other points. And the corresponding log-polar histogram that has m radial and n angular separations will have m * n bins. The log-polar plot can be visualized as a series of concentric circles enclosing a number of wedge bins. The wedge bins are uniform in angular spacing, but vary logarithmically in the radial direction. In other words, in Cartesian space, the inner circles are closer together than the outer circles. But in log space, the circles are spaced uniformly in the radial direction. The figure below taken from [Belongie01] illustrates this. The second picture shows a flattened shape context histogram with the relative darkness of the bins indicating point density.

Shape Context

The advantage of using shape contexts is that several geometric invariances are inherent with this method. Specifically, invariance to translation is clearly built in since all vectors are between points on a shape. Invariance to scaling may be achieved by normalizing all vectors with respect to the mean between all n * (n - 1) pairs of points. Rotational invariance may be achieved by using the tangent vector at every point as the axis to measure angular orientation. This method is more robust to angular variations in images as opposed to using a fixed axis (e.g. x-axis). Given the relative difficulty and computational expense of determining the order of points to arrive at the tangent angle, I have come up with a different method for rotational invariance. This method involves determining the average vector heading from a point to all other points and using this average vector as the axis to measure the angle from for each point. The advantage of this method is that it is easily obtained and is relatively robust to local position variations.

Qualitatively, the shape context for each point gives a relatively precise description of the relative position of a point to all other points. Thus, the shape context can be used as an effective measurement of shape similarity for a corresponding point on another shape to be compared. This implies that there must be a point-to-point correspondence between points on two shapes. To obtain point-to-point correspondence, the log-polar histogram representing the shape context for each point on a shape, may be used to calculate the measurement cost between points on two shapes. Specifically, the chi-squared distance metric is used. Given two shapes each with n sampled points, a cost matrix of size n by n is established where each matrix element is the chi-squared cost between a point on one shape with a point on another shape. The cost equation is given below, where Cij represents the cost between points pi and pj on different shapes. hi(k) represents the kth bin of the normalized histograms at pi and pj. In essence, we are comparing the similarity between histograms at two points. Note that the denominator may go to zero, in which case we just assign a zero to the ratio for a particular k.

Cost Equation

Therefore, the matrix represents the cost between all possible pairs of points. Although two shapes need not have the same number of sampled points, it is easier to work with a balanced cost matrix. For an unbalanced matrix, dummy elements with a constant threshold cost may be added to represent the extra points. This method may also be used even when shapes do have the same number of points by eliminating outliers that do not meet the threshold cost. This application assumes both shapes have the same number of sampled points.

The cost matrix represents exactly, an instance of the weighted bipartite matching problem. This problem may be cast in the more familiar form of the machine scheduling problem, where the goal is to assign one task to every machine so that every task will be tended to. The cost elements of the matrix measure the effectiveness of a machine to perform a specific task, while the objective value aims to maximize the value for all machine and task pairs. In modeling the bipartite matching problem with respect to shape contexts, we wish to minimize the total cost for all point pairs. The bipartite matching problem may be solved in O(n^3) by using the Hungarian algorithm. This application uses a more efficient algorithm provided by [Jonker87] in solving the linear assignment problem.

Having solved the point-to-point correspondence problem, we must now align the shapes so that we may correspond features on one image to another image. [Belongie01] suggests using the thin-plate-spline method of point-set alignment. This technique attempts to model coordinate transformations from one set of points to another by using a weighted combination of thin plate splines centered about each control point that allows a mapping function to be interpolated through these points exactly. The figure below taken from [Donato02] illustrates this.

Thin plate warpage

Implementation wise, the thin-plate-spline is relatively straightforward. It involves solving the following matrix equation:

Thin Plate Equation

Given U(r) = r² * log(r) which is the fundamental solution to the biharmonic equation, and matrix K of size p by p, matrix element Kij is equal to U(|| rij ||), where rij is the radial distance between points pi and pj; i and j are from 0 to p - 1 for p points. P is a p (row) by 3 (col) matrix with all 1's for the first column, x coordinates for the second column, and y coordinates for the third column. Therefore, (Pi1, Pi2) is the (x, y) coordinate of the ith point. O is a 3 by 3 matrix with all 0's. v is a column vector of size p where the entries are the relative distances between corresponding sample and control points in the x-direction. It should be noted that two thin plate splines are needed for a two-dimensional coordinate transformation. Therefore, there will also be a v column vector of size p for the y-direction. Hence, when solving the matrix equation, v will actually be a p by 2 matrix. Also, o is a 3 by 2 matrix with all 0's.

The resultant matrices w and a will be of sizes p by 2 and 3 by 2, respectively. To perform the actual mapping, the following equation is used:

Thin Plate Spline Mapping Function

The first column of w and a are used to determine the new mapping in the x-direction, and second column for the y-direction. A final point to note is that matrices K and P are determined with respect to the sample points to be mapped to the control points.

[Donato02] notes additional improvements to the thin-plate-spline method such as smoothing to handle noise in the point-sets, and computational efficiency. Specifically, an additional term may be added to the K matrix. This added parameter is the regularization parameter "lambda" which is added to the K by lambda * I (identity matrix). By increasing lambda, the amount of smoothing also increases until a least squares affine model is reached. Given that the matrix inversion can be relatively expensive for calculating the mapping function coefficients, simple subsampling may be employed that reduces the size of the K matrix by substituting it with a smaller matrix A with randomly selected values from K. Then interpolating the rest of the points after the mapping process. These additional improvements are not explored in this application.

Operations of the Image Alignment/Feature Locate Program

As detailed above, the code implementation of this application first performs edge detection to retrieve the edge data of images. These edges are then sampled for both test and target point-sets for two images to be aligned. Point-to-point correspondence is performed between the two point-sets to ensure a minimum overall matching between point-sets. The thin-plate-spline interpolation is then performed on the point-sets to provide a closed form mapping function from the test to target point-sets. The figure below shows a screen shot of the application program.

Application Program

Running the Program

The gui has four views, the control, template, test, and feature locate views. The control view allows the user to load and display images for the template and test views. The idea is to match and align the template image to the test image. To run the program, first load images and then perform edge detection on both template and test images. Then press the "Map" button to apply thin-plate-spline mapping. Afterwards, the user can click and create boxes in the template view to select features to map to the test view, which are then shown in the feature locate view. The "Locate Feature(s)" button takes all the user-selected points from the template view, maps it to the test view, groups them into convex hulls according to user selections, and then displays the centroids of each group in the feature locate view.

The "Reduce Outliers & Update" button attempts to reduce the number of outliers during grouping. This is done by removing any points of a group that are greater than a fixed number of standard deviations from the mean radial distance between the group centroid and a point. This is a fairly simple and inflexible technique, as there are more sophistocated methods for outlier removal (e.g. http://research.microsoft.com/~jdunagan/outliers-2002.pdf). But it produces reasonable results for this initial attempt at feature location. I have more work to do in tweaking the shape context and thin-plate-spline algorithms to obtain better results. In doing so, that may eliminate the need to post-process the data by removing outliers.

The "Apply Hysteresis" checkbox allows post-processing of edge detection where image pixels that should be, and have not been suppressed, are done so. Thinning is as it sounds, that is, edges are thinned to its minimum, i.e. single pixel width edges.

Code Implementation

The program was written and compiled in Visual C++ .Net. Several resources were drawn from to implement the application program. Boie-Cox edge detection code was borrowed from http://www.ee.ucl.ac.uk/~icox and adapted to work for this application. Code for solving the linear assignment problem was borrowed from http://www.magiclogic.com/assignment.html. Support for C++ matrix operations was borrowed and adapted from the Matrix TCL Lite package downloaded from http://www.techsoftpl.com (the C++ Boost libraries could have been another alternative). Source code for GUI controls and directory manipulation in MFC was borrowed from [Giannini01]. And code for solving the convex hull for a set of points was borrowed from http://softsurfer.com/Archive/algorithm_0109/algorithm_0109.htm.

Image File Format

The application program only reads and displays .pgm image files, given that this format is easy to parse. The format is such that intensity values are represented in ASCII format. The .pgm image format comes in two versions. A "P2" and "P5" version. The "P2" version is represented by ASCII number strings (e.g. 101). The "P5" version is represented by ASCII in 8 or 16 bit binary. The application program reads both versions. The decision to use .pgm format was based mostly on the requirement of the edge detection code to process grayscale images. Also, extra color information would be superfluous for simple feature location. And lastly, given current time constraints, the decision was made not to delve into common image formats such as .bmp, .jpg, and .gif, etc. More flexible image file support may be added in the future.

Downloadable files

For the source code and executable, download image_align.zip. Also, for .pgm images to test with, download pgm_images.zip.

References

[Belongie01] S. Belongie, J. Malik, and J. Puzicha. Shape matching and object recognition using shape contexts. Technical Report UCB//CSD-00-1128, UC Berkeley, January 2001.

[Bookstein89] F. L. Bookstein. Principal warps: thin-plate splines and decomposition of deformations. IEEE Trans. Pattern Analysis and Machine Intelligence, June 1989.

[Donato02] Gianluca Donato and Serge Belongie. Approximation Methods for Thin Plate Spline Mappings and Principal Warps. http://www-cse.ucsd.edu/~sjb/pami_tps.pdf, 2002.

[Elonen03] Jarno Elonen. Thin Plate Spline editor - an example program in C++. http://elonen.iki.fi/code/tpsdemo, 2003.

[Giannini01] Mario Giannini and Jim Keogh.Windows Programming: Programmer's Notebook. Prentice Hall, 2001.

[Jonker87] R. Jonker and A. Volgenant. "A Shortest Augmenting Path Algorithm for Dense and Sparse Linear Assignment Problems", Computing 38, 325-340, 1987.

[Russell03] Stuart Russell and Peter Norvig. Artificial Intelligence: A Modern Approach. Pearson Education, Inc., 2003.

[Seul00] Michael Seul, Lawrence O'Gorman, and Michael J. Sammon. Practical Algorithms for Image Analysis: Description, Examples, and Code. Cambridge University Press, 2000.

Source code for the Image Alignment Program:

/***************************************************************************

PROJECT:			CS 297 Deliverable 3

FILE:				controlview.cpp

PURPOSE:			Class implementation for forms view where functionality
					for controls such as buttons, edit boxes, etc. are
					implemented.  The CControlView class defines the
					control view for users to interact with.

COMPILER / TOOL:	Visual C++ .NET

PROGRAMMER:			Wallun Chan

START DATE:			4/1/05

***************************************************************************/

#include "stdafx.h"
#include "utility.h"
#include "imagealign.h"
#include "imagealigndoc.h"
#include "edgefinder.h"
#include "ControlView.h"
#include "leftimageview.h"
#include "rightimageview.h"
#include "featureview.h"
#include "directory.h"
#include "convexhull.h"


// CControlView

IMPLEMENT_DYNCREATE(CControlView, CFormView)

CControlView::CControlView() : CFormView(CControlView::IDD)
// Default constructor
{
}

CControlView::~CControlView()
// Destructor
{
}

void CControlView::OnInitialUpdate()
/* PURPOSE:		Initializes GUI controls to their proper states.
*/
{
	CFormView::OnInitialUpdate();
	SetDlgItemText(IDC_SAMPLED_EDIT, "300");
	SetDlgItemText(IDC_ANGULAR_EDIT, "10");
	SetDlgItemText(IDC_RADIAL_EDIT, "8");
	SetDlgItemText(IDC_TEMPLATE_SIGMA_EDIT, "1.5");
	SetDlgItemText(IDC_TEMPLATE_THRESHOLD_EDIT, "12");
	SetDlgItemText(IDC_TEST_SIGMA_EDIT, "1.5");
	SetDlgItemText(IDC_TEST_THRESHOLD_EDIT, "12");
	GetDlgItem(IDC_TEMPLATE_HYS_CHECK)->EnableWindow(FALSE);
	GetDlgItem(IDC_TEMPLATE_THIN_CHECK)->EnableWindow(FALSE);
	GetDlgItem(IDC_TEMPLATE_SHOW_IMAGE_CHECK)->EnableWindow(FALSE);
	GetDlgItem(IDC_TEMPLATE_SHOW_EDGES_CHECK)->EnableWindow(FALSE);
	GetDlgItem(IDC_TEMPLATE_SAMPLED_CHECK)->EnableWindow(FALSE);
	GetDlgItem(IDC_TEMPLATE_DETECT_BUTTON)->EnableWindow(FALSE);
	((CButton*)GetDlgItem(IDC_TEMPLATE_HYS_CHECK))->SetCheck(BST_CHECKED);
	((CButton*)GetDlgItem(IDC_TEMPLATE_THIN_CHECK))->SetCheck(BST_CHECKED);
	((CButton*)GetDlgItem(IDC_TEMPLATE_SHOW_IMAGE_CHECK))->SetCheck(BST_CHECKED);
	((CButton*)GetDlgItem(IDC_TEMPLATE_SHOW_EDGES_CHECK))->SetCheck(BST_CHECKED);
	((CButton*)GetDlgItem(IDC_TEMPLATE_SAMPLED_CHECK))->SetCheck(BST_CHECKED);
	GetDlgItem(IDC_TEST_HYS_CHECK)->EnableWindow(FALSE);
	GetDlgItem(IDC_TEST_THIN_CHECK)->EnableWindow(FALSE);
	GetDlgItem(IDC_TEST_SHOW_IMAGE_CHECK)->EnableWindow(FALSE);
	GetDlgItem(IDC_TEST_SHOW_EDGES_CHECK)->EnableWindow(FALSE);
	GetDlgItem(IDC_TEST_SAMPLED_CHECK)->EnableWindow(FALSE);
	GetDlgItem(IDC_TEST_DETECT_BUTTON)->EnableWindow(FALSE);
	((CButton*)GetDlgItem(IDC_TEST_HYS_CHECK))->SetCheck(BST_CHECKED);
	((CButton*)GetDlgItem(IDC_TEST_THIN_CHECK))->SetCheck(BST_CHECKED);
	((CButton*)GetDlgItem(IDC_TEST_SHOW_IMAGE_CHECK))->SetCheck(BST_CHECKED);
	((CButton*)GetDlgItem(IDC_TEST_SHOW_EDGES_CHECK))->SetCheck(BST_CHECKED);
	((CButton*)GetDlgItem(IDC_TEST_SAMPLED_CHECK))->SetCheck(BST_CHECKED);
	((CButton*)GetDlgItem(IDC_POINT_GROUP_CHECK))->SetCheck(BST_CHECKED);
}

void CControlView::DoDataExchange(CDataExchange* pDX)
// Standard windows initialization
{
	CFormView::DoDataExchange(pDX);
}

BEGIN_MESSAGE_MAP(CControlView, CFormView)
	ON_BN_CLICKED(IDC_TEMPLATE_LOAD_BUTTON, OnBnClickedLoadLeftButton)
	ON_BN_CLICKED(IDC_TEST_LOAD_BUTTON, OnBnClickedLoadRightButton)
	ON_BN_CLICKED(IDC_TEMPLATE_CLEAR_BUTTON, OnBnClickedTemplateClearButton)
	ON_BN_CLICKED(IDC_TEMPLATE_DETECT_BUTTON, OnBnClickedTemplateDetectButton)
	ON_BN_CLICKED(IDC_TEST_CLEAR_BUTTON, OnBnClickedTestClearButton)
	ON_BN_CLICKED(IDC_TEST_DETECT_BUTTON, OnBnClickedTestDetectButton)
	ON_BN_CLICKED(IDC_ALIGN_BUTTON, OnBnClickedAlignButton)
	ON_BN_CLICKED(IDC_TEMPLATE_SHOW_IMAGE_CHECK, OnBnClickedTemplateShowImageCheck)
	ON_BN_CLICKED(IDC_TEMPLATE_SHOW_EDGES_CHECK, OnBnClickedTemplateShowEdgesCheck)
	ON_BN_CLICKED(IDC_TEMPLATE_SAMPLED_CHECK, OnBnClickedTemplateSampledCheck)
	ON_BN_CLICKED(IDC_TEMPLATE_HYS_CHECK, OnBnClickedTemplateHysCheck)
	ON_BN_CLICKED(IDC_TEMPLATE_THIN_CHECK, OnBnClickedTemplateThinCheck)
	ON_BN_CLICKED(IDC_TEST_HYS_CHECK, OnBnClickedTestHysCheck)
	ON_BN_CLICKED(IDC_TEST_THIN_CHECK, OnBnClickedTestThinCheck)
	ON_BN_CLICKED(IDC_TEST_SHOW_IMAGE_CHECK, OnBnClickedTestShowImageCheck)
	ON_BN_CLICKED(IDC_TEST_SHOW_EDGES_CHECK, OnBnClickedTestShowEdgesCheck)
	ON_BN_CLICKED(IDC_TEST_SAMPLED_CHECK, OnBnClickedTestSampledCheck)
	ON_BN_CLICKED(IDC_ERASE_BOXES_BUTTON, OnBnClickedEraseBoxesButton)
	ON_BN_CLICKED(IDC_LOCATE_BUTTON, OnBnClickedLocateButton)
	ON_BN_CLICKED(IDC_FEATURE_RESET_BUTTON, OnBnClickedFeatureResetButton)
	ON_BN_CLICKED(IDC_POINT_GROUP_CHECK, OnBnClickedPointGroupCheck)
	ON_BN_CLICKED(ID_SAVE_TEMPLATE_BUTTON, OnBnClickedSaveTemplateButton)
	ON_BN_CLICKED(ID_SAVE_TEST_BUTTON, OnBnClickedSaveTestButton)
	ON_BN_CLICKED(ID_SAVE_FEATURE_BUTTON, OnBnClickedSaveFeatureButton)
	ON_BN_CLICKED(IDC_LOCATE_REMOVE_OUTLIERS_BUTTON, OnBnClickedLocateRemoveOutliersButton)
END_MESSAGE_MAP()


// CControlView diagnostics

#ifdef _DEBUG
void CControlView::AssertValid() const
{
	CFormView::AssertValid();
}

void CControlView::Dump(CDumpContext& dc) const
{
	CFormView::Dump(dc);
}
#endif //_DEBUG

void CControlView::updateTemplateViewState(BOOL updateSample)
/* PURPOSE:		Checks the state of the GUI and based on user
				selections (e.g. check boxes, etc.) will update
				template view to reflect latest state.

   RECEIVES:	Boolean to determine whether pixels should
				be sampled from edge data of template image.
*/
{
	BYTE* image;
	BYTE* edgeMap;
	CString sigmaStr;
	CString thresholdStr;
	CString sampleStr;

	// Get user input values
	GetDlgItemText(IDC_TEMPLATE_SIGMA_EDIT, sigmaStr);
	GetDlgItemText(IDC_TEMPLATE_THRESHOLD_EDIT, thresholdStr);
	GetDlgItemText(IDC_SAMPLED_EDIT, sampleStr);
	double sigma = atof(sigmaStr);
	int threshold = atoi(thresholdStr);
	int oldThreshold = threshold;
	int samples = atoi(sampleStr);

	// Check to see if checkboxes are checked, or not.
	int hysCheck = ((CButton*)GetDlgItem(IDC_TEMPLATE_HYS_CHECK))->GetCheck();
	int thinCheck =	((CButton*)GetDlgItem(IDC_TEMPLATE_THIN_CHECK))->GetCheck();
	int imageCheck = ((CButton*)GetDlgItem(IDC_TEMPLATE_SHOW_IMAGE_CHECK))->GetCheck();
	int edgeCheck = ((CButton*)GetDlgItem(IDC_TEMPLATE_SHOW_EDGES_CHECK))->GetCheck();
	int sampleCheck = ((CButton*)GetDlgItem(IDC_TEMPLATE_SAMPLED_CHECK))->GetCheck();

	CimagealignDoc* doc = (CimagealignDoc*)GetDocument();

	// Release allocated memory for template view image
	delete [] doc->templateImage;

	// Get edge detected template image
	image = doc->templateDetector->get_image();

	/* Get edge detected data based on user selections of
	   hysteresis, thinning, etc. for template image view. */
	edgeMap = getEdgeMap(doc->templateDetector, hysCheck, thinCheck, sigma, &threshold);

	/* If a zero for threshold is passed in, the edge detection code will
	   estimate a best effort threshold to use. So update new threshold
	   for threshold if this is the case. */
	if (oldThreshold == 0)
	{	thresholdStr.Format("%d", threshold);
		SetDlgItemText(IDC_TEMPLATE_THRESHOLD_EDIT, thresholdStr);
		CString editStr;
		editStr.Format("Template view threshold estimated to be %s.",
			thresholdStr);
		MessageBox(editStr, "Threshold estimation ...", MB_OK | MB_ICONWARNING);
	}

	// If updateSample is true, then (re)sample pixels from edges.
	if (updateSample)
	{	doc->templateSampledPts.clear();
		getSampledPoints(doc->templateSampledPts, edgeMap,
			doc->templateWidth, doc->templateHeight, samples);
	}

	/* Based on user selections, this function call generates the proper image
	   to be displayed in the template image view. */
	setDisplayImage(doc->templateImage, image, edgeMap, doc->templateSampledPts,
		doc->templateWidth, doc->templateHeight, imageCheck, edgeCheck, sampleCheck);

	// Redraw template view
	CLeftImageView* view = (CLeftImageView*)doc->getView(RUNTIME_CLASS(CLeftImageView));
	view->Invalidate();
	view->UpdateWindow();
}

void CControlView::updateTestViewState(BOOL updateSample)
/* PURPOSE:		Checks the state of the GUI and based on user
				selections (e.g. check boxes, etc.) will update
				test view to reflect latest state.

   RECEIVES:	Boolean to determine whether pixels should
				be sampled from edge data of test image.
*/
{
	BYTE* image;
	BYTE* edgeMap;
	CString sigmaStr;
	CString thresholdStr;
	CString sampleStr;

	// Get user selections and input in GUI controls
	GetDlgItemText(IDC_TEST_SIGMA_EDIT, sigmaStr);
	GetDlgItemText(IDC_TEST_THRESHOLD_EDIT, thresholdStr);
	GetDlgItemText(IDC_SAMPLED_EDIT, sampleStr);
	double sigma = atof(sigmaStr);
	int threshold = atoi(thresholdStr);
	int oldThreshold = threshold;
	int samples = atoi(sampleStr);

	// Check to see if checkboxes are checked, or not.
	int hysCheck = ((CButton*)GetDlgItem(IDC_TEST_HYS_CHECK))->GetCheck();
	int thinCheck =	((CButton*)GetDlgItem(IDC_TEST_THIN_CHECK))->GetCheck();
	int imageCheck = ((CButton*)GetDlgItem(IDC_TEST_SHOW_IMAGE_CHECK))->GetCheck();
	int edgeCheck = ((CButton*)GetDlgItem(IDC_TEST_SHOW_EDGES_CHECK))->GetCheck();
	int sampleCheck = ((CButton*)GetDlgItem(IDC_TEST_SAMPLED_CHECK))->GetCheck();

	CimagealignDoc* doc = (CimagealignDoc*)GetDocument();

	// Release any memory for test image view
	delete [] doc->testImage;

	// Get edge detected test image
	image = doc->testDetector->get_image();

	/* Get edge detected data based on user selections of hysteresis,
	   thinning, etc for test image view. */
	edgeMap = getEdgeMap(doc->testDetector, hysCheck, thinCheck, sigma, &threshold);

	/* If a zero for threshold is passed in, the edge detection code will
	   estimate a best effort threshold to use. So update new threshold
	   for threshold if this is the case. */
	if (oldThreshold == 0)
	{	thresholdStr.Format("%d", threshold);
		SetDlgItemText(IDC_TEST_THRESHOLD_EDIT, thresholdStr);
		CString editStr;
		editStr.Format("Test view threshold estimated to be %s.",
			thresholdStr);
		MessageBox(editStr, "Threshold estimation ...", MB_OK | MB_ICONWARNING);
	}

	// If updateSample is true, then (re)sample pixels from edges.
	if (updateSample)
	{	doc->testSampledPts.clear();
		getSampledPoints(doc->testSampledPts, edgeMap, doc->testWidth, doc->testHeight, samples);
	}

	/* Based on user selections, this function call generates the proper image
	   to be displayed in the test image view. */
	setDisplayImage(doc->testImage, image, edgeMap, doc->testSampledPts,
		doc->testWidth, doc->testHeight, imageCheck, edgeCheck, sampleCheck);

	// Redraw template view
	CRightImageView* view = (CRightImageView*)doc->getView(RUNTIME_CLASS(CRightImageView));
	view->Invalidate();
	view->UpdateWindow();
}


// CControlView message handlers

void CControlView::OnBnClickedLoadLeftButton()
/* PURPOSE:		Allows user to select a .pgm image to be displayed
				in template image view.
*/
{
	// Open file display and selection dialog
	CFileDialog dialog(1, 0, 0, OFN_HIDEREADONLY | OFN_OVERWRITEPROMPT, pgmFilter);
	if (dialog.DoModal() != IDOK)
	{	return;
	}

	int width;
	int height;

	/* Read .pgm image for image dimensions and return a byte string
	   containing the image data. */
	BYTE* bytes = pgmByteRead(dialog.GetPathName(), width, height);

	CimagealignDoc* doc = (CimagealignDoc*)GetDocument();

	// Release previously allocated memory
	doc->templateDetector->freeImage();
	delete [] doc->templateImage;
	doc->templateSampledPts.clear();
	doc->boxes.clear();

	// Initialize edge detector with image to be edge detected
	doc->templateDetector->setImage(bytes, width, height);

	/* Set document data so that the template view will update to the
	   latest image. */
	doc->templateImage = (BYTE*)rgbFormat(bytes, width * height);
	doc->templateWidth = width;
	doc->templateHeight = height;

	// Set border title
	doc->tempFileName.Format("   (Current image: %s)",
		GetSubPath(dialog.GetPathName()));

	/* Flag to check if template image is mapped to
	   test image. */
	doc->isMapped = false;

	// Redraw template image view
	CLeftImageView* view = (CLeftImageView*)doc->
		getView(RUNTIME_CLASS(CLeftImageView));
	view->Invalidate();
	view->UpdateWindow();

	// Resize scroll length to fit image
	CSize sizeTotal;
	sizeTotal.cx = (long)(SCROLL_FACTOR * width);
	sizeTotal.cy = (long)(SCROLL_FACTOR * height);
	view->SetScrollSizes(MM_TEXT, sizeTotal);

	// Enable the proper GUI controls
	GetDlgItem(IDC_TEMPLATE_DETECT_BUTTON)->EnableWindow();
	GetDlgItem(IDC_TEMPLATE_HYS_CHECK)->EnableWindow();
	GetDlgItem(IDC_TEMPLATE_THIN_CHECK)->EnableWindow();
	GetDlgItem(IDC_TEMPLATE_SHOW_IMAGE_CHECK)->EnableWindow();
	GetDlgItem(IDC_TEMPLATE_SHOW_EDGES_CHECK)->EnableWindow();
	GetDlgItem(IDC_TEMPLATE_SAMPLED_CHECK)->EnableWindow();
}

void CControlView::OnBnClickedLoadRightButton()
/* PURPOSE:		Allows user to select a .pgm image to be displayed
				in test view.
*/
{
	// Open file display and selection dialog
	CFileDialog dialog(1, 0, 0, OFN_HIDEREADONLY | OFN_OVERWRITEPROMPT, pgmFilter);
	if (dialog.DoModal() != IDOK)
	{	return;
	}

	int width;
	int height;

	/* Read .pgm image for image dimensions and return a byte string
	   containing the image data. */
	BYTE* bytes = pgmByteRead(dialog.GetPathName(), width, height);

	CimagealignDoc* doc = (CimagealignDoc*)GetDocument();

	// Release previously allocated memory
	doc->testDetector->freeImage();
	delete [] doc->testImage;
	doc->testSampledPts.clear();

	// Initialize edge detector with image to be edge detected
	doc->testDetector->setImage(bytes, width, height);

	/* Set document data so that the test view will update to the
	   latest image. */
	doc->testImage = (BYTE*)rgbFormat(bytes, width * height);
	doc->testWidth = width;
	doc->testHeight = height;

	// Set border title
	doc->testFileName.Format("   (Current image: %s)",
		GetSubPath(dialog.GetPathName()));

	/* Flag to check if template image is mapped to
	   test image. */
	doc->isMapped = false;

	// Redraw template image view
	CLeftImageView* view = (CLeftImageView*)doc->getView(RUNTIME_CLASS(CRightImageView));
	view->Invalidate();
	view->UpdateWindow();

	// Resize scroll length to fit image
	CSize sizeTotal;
	sizeTotal.cx = (long)(1.25 * width);
	sizeTotal.cy = (long)(1.25 * height);
	view->SetScrollSizes(MM_TEXT, sizeTotal);

	// Enable the proper GUI controls
	GetDlgItem(IDC_TEST_DETECT_BUTTON)->EnableWindow();
	GetDlgItem(IDC_TEST_HYS_CHECK)->EnableWindow();
	GetDlgItem(IDC_TEST_THIN_CHECK)->EnableWindow();
	GetDlgItem(IDC_TEST_SHOW_IMAGE_CHECK)->EnableWindow();
	GetDlgItem(IDC_TEST_SHOW_EDGES_CHECK)->EnableWindow();
	GetDlgItem(IDC_TEST_SAMPLED_CHECK)->EnableWindow();
}

void CControlView::OnBnClickedTemplateClearButton()
/* PURPOSE:		Clears template view display and releases memory
				allocated for image that was displayed.
*/
{
	CimagealignDoc* doc = (CimagealignDoc*)GetDocument();

	// Free allocated memory for template image view
	doc->templateDetector->freeImage();
	delete [] doc->templateImage;
	doc->templateImage = 0;
	doc->templateSampledPts.clear();
	doc->boxes.clear();

	// Reset border title
	doc->tempFileName.Empty();

	/* Flag to check if template image is mapped to
	   test image. */
	doc->isMapped = false;

	// Update template window
	CLeftImageView* view = (CLeftImageView*)doc->getView(RUNTIME_CLASS(CLeftImageView));
	view->Invalidate();
	view->UpdateWindow();

	// Reset scroll length
	view->SetScrollSizes(MM_TEXT, CSize(DEFAULT_SCROLL, DEFAULT_SCROLL));

	// Reset GUI controls
	GetDlgItem(IDC_TEMPLATE_DETECT_BUTTON)->EnableWindow(FALSE);
	GetDlgItem(IDC_TEMPLATE_HYS_CHECK)->EnableWindow(FALSE);
	GetDlgItem(IDC_TEMPLATE_THIN_CHECK)->EnableWindow(FALSE);
	GetDlgItem(IDC_TEMPLATE_SHOW_IMAGE_CHECK)->EnableWindow(FALSE);
	GetDlgItem(IDC_TEMPLATE_SHOW_EDGES_CHECK)->EnableWindow(FALSE);
	GetDlgItem(IDC_TEMPLATE_SAMPLED_CHECK)->EnableWindow(FALSE);
	((CButton*)GetDlgItem(IDC_TEMPLATE_HYS_CHECK))->SetCheck(BST_CHECKED);
	((CButton*)GetDlgItem(IDC_TEMPLATE_THIN_CHECK))->SetCheck(BST_CHECKED);
	((CButton*)GetDlgItem(IDC_TEMPLATE_SHOW_IMAGE_CHECK))->SetCheck(BST_CHECKED);
	((CButton*)GetDlgItem(IDC_TEMPLATE_SHOW_EDGES_CHECK))->SetCheck(BST_CHECKED);
	((CButton*)GetDlgItem(IDC_TEMPLATE_SAMPLED_CHECK))->SetCheck(BST_CHECKED);
}

void CControlView::OnBnClickedTemplateDetectButton()
/* PURPOSE:		Invokes edge detection and the display of edge
				detection data for template image view.
*/
{
	// Update template view to reflect latest image
	updateTemplateViewState();

	// Set flag to indicate that images need to be mapped
	((CimagealignDoc*)GetDocument())->isMapped = false;
}

void CControlView::OnBnClickedTestClearButton()
/* PURPOSE:		Clears test view display and releases memory
				allocated for image that was displayed.
*/
{
	CimagealignDoc* doc = (CimagealignDoc*)GetDocument();

	// Free allocated memory for template image view
	doc->testDetector->freeImage();
	delete [] doc->testImage;
	doc->testImage = 0;
	doc->testSampledPts.clear();

	// Reset border title
	doc->testFileName.Empty();

	/* Flag to check if template image is mapped to
	   test image. */
	doc->isMapped = false;

	// Update template window
	CRightImageView* view = (CRightImageView*)doc->getView(RUNTIME_CLASS(CRightImageView));
	view->Invalidate();
	view->UpdateWindow();

	// Reset scroll length
	view->SetScrollSizes(MM_TEXT, CSize(DEFAULT_SCROLL, DEFAULT_SCROLL));

	// Reset GUI controls
	GetDlgItem(IDC_TEST_DETECT_BUTTON)->EnableWindow(FALSE);
	GetDlgItem(IDC_TEST_HYS_CHECK)->EnableWindow(FALSE);
	GetDlgItem(IDC_TEST_THIN_CHECK)->EnableWindow(FALSE);
	GetDlgItem(IDC_TEST_SHOW_IMAGE_CHECK)->EnableWindow(FALSE);
	GetDlgItem(IDC_TEST_SHOW_EDGES_CHECK)->EnableWindow(FALSE);
	GetDlgItem(IDC_TEST_SAMPLED_CHECK)->EnableWindow(FALSE);
	((CButton*)GetDlgItem(IDC_TEST_HYS_CHECK))->SetCheck(BST_CHECKED);
	((CButton*)GetDlgItem(IDC_TEST_THIN_CHECK))->SetCheck(BST_CHECKED);
	((CButton*)GetDlgItem(IDC_TEST_SHOW_IMAGE_CHECK))->SetCheck(BST_CHECKED);
	((CButton*)GetDlgItem(IDC_TEST_SHOW_EDGES_CHECK))->SetCheck(BST_CHECKED);
	((CButton*)GetDlgItem(IDC_TEST_SAMPLED_CHECK))->SetCheck(BST_CHECKED);
}

void CControlView::OnBnClickedTestDetectButton()
/* PURPOSE:		Invokes edge detection and the display of edge
				detection data for test image view.
*/
{
	// Update test view to reflect latest image
	updateTestViewState();

	/* Flag to check if template image is mapped to
	   test image. */
	((CimagealignDoc*)GetDocument())->isMapped = false;
}

void CControlView::OnBnClickedAlignButton()
/* PURPOSE:		This function takes the sampled points of edge data
				from the template image and aligns it with the sampled
				points of the test image.  Based on this alignment,
				the rest of the edge data is mapped over to the
				test image.  The mapped image is then displayed.
*/
{
	CimagealignDoc* doc = (CimagealignDoc*)GetDocument();

	// Check to see if images have been loaded before aligning
	if (doc->templateImage == 0)
	{	MessageBox("Template image not loaded.",
			"Load image ...", MB_OK | MB_ICONWARNING);
		return;
	}
	if (doc->testImage == 0)
	{	MessageBox("Test image not loaded.",
			"Load image ...", MB_OK | MB_ICONWARNING);
		return;
	}
	if (doc->templateSampledPts.size() == 0)
	{	MessageBox("Template image not edge detected.",
			"Load image ...", MB_OK | MB_ICONWARNING);
		return;
	}
	if (doc->testSampledPts.size() == 0)
	{	MessageBox("Test image not edge detected.",
			"Load image ...", MB_OK | MB_ICONWARNING);
		return;
	}

	// Get user selected and input data from GUI.
	CString sampNumStr;
	CString angBinStr;
	CString radialBinStr;
	GetDlgItemText(IDC_SAMPLED_EDIT, sampNumStr);
	GetDlgItemText(IDC_ANGULAR_EDIT, angBinStr);
	GetDlgItemText(IDC_RADIAL_EDIT, radialBinStr);
	int samples = atoi(sampNumStr);
	int angBinNum = atoi(angBinStr);
	int radialBinNum = atoi(radialBinStr);

	Map tempContext;
	Map testContext;

	// Make copies of sampled points
	VecIntPr tempPts = doc->templateSampledPts;
	VecIntPr testPts = doc->testSampledPts;

	/* Generate shape context data for sampled pixels from edges
	   of template and test images. */
	setShapeContextPair(tempContext, testContext,
		tempPts, testPts, radialBinNum, angBinNum);

	// Update progress bar
	SendDlgItemMessage(IDC_PROGRESS, PBM_SETPOS, 25, 0);

	// Find point-to-point correspondence between sampled points
	matchPoints(tempContext, testContext, tempPts, testPts);

	// Update progress bar
	SendDlgItemMessage(IDC_PROGRESS, PBM_SETPOS, 75, 0);

	// Map template sampled edge points to test image space
	doc->tps.setPoints(testPts, tempPts);

	// Flag to indicate that template and test images have been mapped
	doc->isMapped = true;

	// Update progress bar
	SendDlgItemMessage(IDC_PROGRESS, PBM_SETPOS, 100, 0);

	MessageBox("Template image mapped to test image.",
		"Images are mapped ...", MB_OK | MB_ICONWARNING);

	// Reset progress bar
	SendDlgItemMessage(IDC_PROGRESS, PBM_SETPOS, 0, 0);
}

void CControlView::OnBnClickedTemplateShowImageCheck()
/* PURPOSE:		To show or not show loaded image in template
				image view.
*/
{
	/* Update template view to reflect latest image
	   without point sampling. */
	updateTemplateViewState(FALSE);
}

void CControlView::OnBnClickedTemplateShowEdgesCheck()
/* PURPOSE:		To show or not show edge data of image in template
				image view.
*/
{
	/* Update template view to reflect latest image
	   without point sampling. */
	updateTemplateViewState(FALSE);
}

void CControlView::OnBnClickedTemplateSampledCheck()
/* PURPOSE:		To show or not show sampled pixels of edges of image
				in template image view.
*/
{
	/* Update template view to reflect latest image
	   without point sampling. */
	updateTemplateViewState(FALSE);
}

void CControlView::OnBnClickedTemplateHysCheck()
/* PURPOSE:		To show or not show hysteresis applied edge data of
				image in template image view.
*/
{
	/* Update template view to reflect latest image
	   without point sampling. */
	updateTemplateViewState(FALSE);
}

void CControlView::OnBnClickedTemplateThinCheck()
/* PURPOSE:		To show or not show thinned edge data of
				image in template image view.
*/
{
	/* Update template view to reflect latest image
	   without point sampling. */
	updateTemplateViewState(FALSE);
}

void CControlView::OnBnClickedTestHysCheck()
/* PURPOSE:		To show or not show hysteresis applied edge data of
				image in test image view.
*/
{
	/* Update test view to reflect latest image
	   without point sampling. */
	updateTestViewState(FALSE);
}

void CControlView::OnBnClickedTestThinCheck()
/* PURPOSE:		To show or not show thinned edge data of
				image in test image view.
*/
{
	/* Update test view to reflect latest image
	   without point sampling. */
	updateTestViewState(FALSE);
}

void CControlView::OnBnClickedTestShowImageCheck()
/* PURPOSE:		To show or not show loaded image in test
				image view.
*/
{
	/* Update test view to reflect latest image
	   without point sampling. */
	updateTestViewState(FALSE);
}

void CControlView::OnBnClickedTestShowEdgesCheck()
/* PURPOSE:		To show or not show edge data of image in test
				image view.
*/
{
	/* Update test view to reflect latest image
	   without point sampling. */
	updateTestViewState(FALSE);
}

void CControlView::OnBnClickedTestSampledCheck()
/* PURPOSE:		To show or not show sampled pixels of edges of image
				in test image view.
*/
{
	/* Update test view to reflect latest image
	   without point sampling. */
	updateTestViewState(FALSE);
}

void CControlView::OnBnClickedEraseBoxesButton()
/* PURPOSE:		Erases rectangles that are drawn in template
				image view.
*/
{
	CimagealignDoc* doc = (CimagealignDoc*)GetDocument();
	doc->boxes.clear();
	doc->ptGroup.clear();
	CLeftImageView* view = (CLeftImageView*)doc->
		getView(RUNTIME_CLASS(CLeftImageView));
	view->Invalidate();
	view->UpdateWindow();
}

void CControlView::featureLocate(BOOL showOutliers)
{
	CimagealignDoc* doc = (CimagealignDoc*)GetDocument();

	/* Used to store coordinates of all user selected points in
	   template view. */
	VecIntPr mapped;

	doc->ptGroup.clear();
	doc->centroids.clear();

	/* Goes through all points contained within rectangles
	   drawn by user and stores them in a list to be displayed
	   in feature locate view. */
	unsigned n;
	for (n = 0; n < doc->boxes.size(); n++)
	{	int left = doc->boxes[n].first.first;
		int top	= doc->boxes[n].first.second;
		int right = doc->boxes[n].second.first;
		int bottom = doc->boxes[n].second.second;

		/* Used to store groups of user selected sampled points in
		   template image.  Difference between this and "mapped" is
		   that "mapped" is used to show the mapped points, whereas
		   "pts" is used to draw the convex hull around the points
		   stored in "map". */
		VecIntPr pts;

		// Checks for points within user-specified rectangle
		unsigned i;
		for (i = 0; i < doc->templateSampledPts.size(); i++)
		{	IntPr pt = doc->templateSampledPts[i];
			if (pt.first > left && pt.first < right &&
				pt.second > top - BORDER_HEIGHT &&
				pt.second < bottom - BORDER_HEIGHT)
			{	pts.push_back(pt);
			}
		}
		doc->tps.mapPoints(pts);

		// Show all points, including extreme outliers
		if (showOutliers)
		{	// Store centroids of groups of points
			doc->centroids.push_back(centroid(pts));

			/* Store gift-wrapped points to be displayed in feature
			   locate view. */
			doc->ptGroup.push_back(convexHull(pts));
			mapped.insert(mapped.end(), pts.begin(), pts.end());
		}
		/* Remove outliers that are not within a specified number
		   of standard deviations from the mean. */
		else
		{	VecIntPr adjPts = removeOutliers(pts);
			doc->centroids.push_back(centroid(adjPts));
			doc->ptGroup.push_back(convexHull(adjPts));
			mapped.insert(mapped.end(), adjPts.begin(), adjPts.end());
		}
	}

	doc->featureWidth = doc->testWidth;
	doc->featureHeight = doc->testHeight;
	delete [] doc->featureImage;

	/* Overlay mapped points onto test image edge detected data and
	   save to document data to be displayed in feature locate view. */
	doc->featureImage = (BYTE*)overlay(doc->testDetector->get_edge_map(),
		doc->testWidth, doc->testHeight, mapped, BLUE, YLW);
}

void CControlView::OnBnClickedLocateButton()
/* PURPOSE:		Based on rectangles drawn in the template image
				by user, this function maps the sampled points
				to the test image space, that are contained
				within these rectangles and displays them
				in the feature locate view.
*/
{
	CimagealignDoc* doc = (CimagealignDoc*)GetDocument();
	if (!doc->isMapped)
	{	MessageBox("Images have not been mapped.",
			"Map image ...", MB_OK | MB_ICONWARNING);
		return;
	}

	/* Locate user-selected features on template image to
	   be displayed in feature locate view; show outliers. */
	featureLocate(TRUE);

	// Update feature locate view
	CFeatureView* view = (CFeatureView*)doc->getView(RUNTIME_CLASS(CFeatureView));
	view->Invalidate();
	view->UpdateWindow();

	// Resize scroll length to fit mapped image
	CSize sizeTotal;
	sizeTotal.cx = (long)(SCROLL_FACTOR * doc->featureWidth);
	sizeTotal.cy = (long)(SCROLL_FACTOR * doc->featureHeight);
	view->SetScrollSizes(MM_TEXT, sizeTotal);
}

void CControlView::OnBnClickedLocateRemoveOutliersButton()
{
	CimagealignDoc* doc = (CimagealignDoc*)GetDocument();
	if (!doc->isMapped)
	{	MessageBox("Images have not been mapped.",
			"Map image ...", MB_OK | MB_ICONWARNING);
		return;
	}

	/* Locate user-selected features on template image to
	   be displayed in feature locate view; remove outliers. */
	featureLocate(FALSE);

	// Update feature locate view
	CFeatureView* view = (CFeatureView*)doc->getView(RUNTIME_CLASS(CFeatureView));
	view->Invalidate();
	view->UpdateWindow();

	// Resize scroll length to fit mapped image
	CSize sizeTotal;
	sizeTotal.cx = (long)(SCROLL_FACTOR * doc->featureWidth);
	sizeTotal.cy = (long)(SCROLL_FACTOR * doc->featureHeight);
	view->SetScrollSizes(MM_TEXT, sizeTotal);
}

void CControlView::OnBnClickedFeatureResetButton()
/* PURPOSE:		Clears display and releases memory allocated for
				image displayed in feature locate view.
*/
{
	CimagealignDoc* doc = (CimagealignDoc*)GetDocument();
	delete [] doc->featureImage;
	doc->featureImage = 0;
	doc->ptGroup.clear();
	CFeatureView* view = (CFeatureView*)doc->getView(RUNTIME_CLASS(CFeatureView));
	view->Invalidate();
	view->UpdateWindow();
	view->SetScrollSizes(MM_TEXT, CSize(DEFAULT_SCROLL, DEFAULT_SCROLL));
}

void CControlView::OnBnClickedPointGroupCheck()
/* PURPOSE:		Based on user selection, this function controls the display
				of the grouping of selected points in the template image
				view that are mapped over to the test view shown in
				the feature locate view.
*/
{
	CimagealignDoc* doc = (CimagealignDoc*)GetDocument();
	if (doc->featureImage == 0)
	{	MessageBox("Feature(s) not located.",
			"No features ...", MB_OK | MB_ICONWARNING);
		return;
	}

	/* If group checkbox is checked, show convex hull of mapped points
	   in feature locate view. */
	doc->showGrouping = ((CButton*)GetDlgItem(IDC_POINT_GROUP_CHECK))->GetCheck();

	// Update feature locate view
	CFeatureView* view = (CFeatureView*)doc->getView(RUNTIME_CLASS(CFeatureView));
	view->Invalidate();
	view->UpdateWindow();
}

void CControlView::OnBnClickedSaveTemplateButton()
/* PURPOSE:		Allows user to save displayed image in template
				image view as .bmp image file.
*/
{
	CimagealignDoc* doc = (CimagealignDoc*)GetDocument();
	if (doc->templateImage == 0)
	{	MessageBox("Template image has not been loaded.",
			"Load template image ...", MB_OK | MB_ICONWARNING);
		return;
	}
	CFileDialog dialog(0, 0, 0, OFN_HIDEREADONLY | OFN_OVERWRITEPROMPT, bmpFilter);
	if (dialog.DoModal() != IDOK)
	{	return;
	}
	imageOutput(dialog.GetPathName(), doc->templateImage,
		doc->templateWidth, doc->templateHeight);
}

void CControlView::OnBnClickedSaveTestButton()
/* PURPOSE:		Allows user to save displayed image in test
				image view as .bmp image file.
*/
{
	CimagealignDoc* doc = (CimagealignDoc*)GetDocument();
	if (doc->testImage == 0)
	{	MessageBox("Test image has not been loaded.",
			"Load test image ...", MB_OK | MB_ICONWARNING);
		return;
	}
	CFileDialog dialog(0, 0, 0, OFN_HIDEREADONLY | OFN_OVERWRITEPROMPT, bmpFilter);
	if (dialog.DoModal() != IDOK)
	{	return;
	}
	imageOutput(dialog.GetPathName(), doc->testImage,
		doc->testWidth, doc->testHeight);
}

void CControlView::OnBnClickedSaveFeatureButton()
/* PURPOSE:		Allows user to save displayed image in feature
				locate image view as .bmp image file.
*/
{
	CimagealignDoc* doc = (CimagealignDoc*)GetDocument();
	if (doc->featureImage == 0)
	{	MessageBox("Features have not been located.",
			"Features not located ...", MB_OK | MB_ICONWARNING);
		return;
	}
	CFileDialog dialog(0, 0, 0, OFN_HIDEREADONLY | OFN_OVERWRITEPROMPT, bmpFilter);
	if (dialog.DoModal() != IDOK)
	{	return;
	}
	imageOutput(dialog.GetPathName(), doc->featureImage,
		doc->featureWidth, doc->featureHeight);
}





/***************************************************************************

PROJECT:			CS 297 Deliverable 3

FILE:				imagealigndoc.h

PURPOSE:			Document class definition to store all important
					data used by the template, test, map, and feature
					locate views.

COMPILER / TOOL:	Visual C++ .NET

PROGRAMMER:			Wallun Chan

START DATE:			4/1/05

***************************************************************************/

#pragma once

#include "edgefinder.h"
#include "utility.h"
#include "tps.h"

class CimagealignDoc : public CDocument
{
protected: // create from serialization only
	CimagealignDoc();
	DECLARE_DYNCREATE(CimagealignDoc)

// Attributes
public:
	// Detector for template view
	EdgeDetector* templateDetector;

	// Detector for test view
	EdgeDetector* testDetector;

	// Thin-plate-spline algorithm
	TPSSolver tps;

	// Points sampled from edges of template image
	VecIntPr templateSampledPts;

	// Points sampled from edges of test image
	VecIntPr testSampledPts;

	// Centroids of mapped user-selected points
	VecIntPr centroids;

	// Coordinates of user-specified selection rectangle
	Boxes boxes;

	//
	PtGroup ptGroup;

	/* Used to track the upper-left and lower-right points
	   of a user-specified selection rectangle. */
	IntPr boxBoundary[2];

	// Pixel data of template image
	BYTE* templateImage;

	// Pixel data of test image
	BYTE* testImage;

	// Pixel data of feature located image
	BYTE* featureImage;

	// Pixel dimensions of template image
	int templateWidth;
	int templateHeight;

	// Pixel dimensions of test image
	int testWidth;
	int testHeight;

	// Pixel dimensions of feature located image
	int featureWidth;
	int featureHeight;

	/* Used to track the number of user clicks inside template
	   view area. */
	int pointCount;

	// Checks to see template and test images are mapped
	BOOL isMapped;

	// Checks to see if mapped points should be convex hulled
	BOOL showGrouping;

	// Checks to see if outliers should be removed from mapped points
	BOOL removeOutliers;

	// Border title text of template view
	CString tempBorderTitle;

	// Border title text of template view
	CString testBorderTitle;

	// File name of current image displayed in template view
	CString tempFileName;

	// File name of current image displayed in test view
	CString testFileName;

// Operations
public:

// Overrides
public:
	virtual BOOL OnNewDocument();
	CView* getView(CRuntimeClass* runtimeClass);

// Implementation
public:
	virtual ~CimagealignDoc();
#ifdef _DEBUG
	virtual void AssertValid() const;
	virtual void Dump(CDumpContext& dc) const;
#endif

protected:

// Generated message map functions
protected:
	DECLARE_MESSAGE_MAP()
};





/***************************************************************************

PROJECT:			CS 297 Deliverable 3

FILE:				imagealigndoc.cpp

PURPOSE:			Document class implementation to store data
					used by the template, test, map, and feature
					locate views.

COMPILER / TOOL:	Visual C++ .NET

PROGRAMMER:			Wallun Chan

START DATE:			4/1/05

***************************************************************************/

#include "stdafx.h"
#include "imagealign.h"
#include "imagealignDoc.h"

#ifdef _DEBUG
#define new DEBUG_NEW
#endif


// CimagealignDoc

IMPLEMENT_DYNCREATE(CimagealignDoc, CDocument)

BEGIN_MESSAGE_MAP(CimagealignDoc, CDocument)
END_MESSAGE_MAP()


// CimagealignDoc construction / destruction

CimagealignDoc::CimagealignDoc() : templateImage(0),
testImage(0), featureImage(0), templateWidth(0),
templateHeight(0), testWidth(0), testHeight(0),
featureWidth(0), featureHeight(0), pointCount(0),
isMapped(FALSE), removeOutliers(FALSE),
showGrouping(TRUE)
// Default constructor
{
	templateDetector = new EdgeDetector();
	testDetector = new EdgeDetector();
}

CimagealignDoc::~CimagealignDoc()
// Destructor to clear all document data
{
	delete templateDetector;
	delete testDetector;
	delete [] templateImage;
	delete [] testImage;
	templateSampledPts.clear();
	testSampledPts.clear();
	templateWidth = 0;
	templateHeight = 0;
	testWidth = 0;
	testHeight = 0;
}

BOOL CimagealignDoc::OnNewDocument()
/* PURPOSE:		Resets entire application.
*/
{
	if (!CDocument::OnNewDocument())
		return FALSE;
	templateSampledPts.clear();
	testSampledPts.clear();
	boxes.clear();
	ptGroup.clear();
	delete [] templateImage;
	delete [] testImage;
	delete [] featureImage;
	templateImage = 0;
	testImage = 0;
	featureImage = 0;
	isMapped = false;
	UpdateAllViews(0);
	return TRUE;
}


// CimagealignDoc diagnostics

#ifdef _DEBUG
void CimagealignDoc::AssertValid() const
{
	CDocument::AssertValid();
}

void CimagealignDoc::Dump(CDumpContext& dc) const
{
	CDocument::Dump(dc);
}
#endif //_DEBUG


CView* CimagealignDoc::getView(CRuntimeClass* runtimeClass)
/* PURPOSE:		Searches for and returns specific view attached
				to parent window (e.g. template view, test view,
				map view, or feature locate view.

   RECEIVES:	Run-time class of view to retrieve.

   RETURNS:		Specific view attached to parent window.
*/
{
	POSITION pos = GetFirstViewPosition();
	while (pos != NULL)
	{	CView* pView = GetNextView(pos);
		if (pView->GetRuntimeClass() == runtimeClass)
		{	return pView;
		}
	}
	return NULL;
}


// CimagealignDoc commands





/***************************************************************************

PROJECT:			CS 297 Deliverable 3

FILE:				leftimageview.cpp

PURPOSE:			Class implementation for template view where user
					may select and display an image to be mapped to
					the test image.

COMPILER / TOOL:	Visual C++ .NET

PROGRAMMER:			Wallun Chan

START DATE:			4/1/05

***************************************************************************/

#include "stdafx.h"
#include "imagealign.h"
#include "LeftImageView.h"
#include "imagealigndoc.h"
#include "edgefinder.h"


// CLeftImageView

IMPLEMENT_DYNCREATE(CLeftImageView, CScrollView)

CLeftImageView::CLeftImageView()
// Default constructor
{
}

CLeftImageView::~CLeftImageView()
// Destructor
{
}


BEGIN_MESSAGE_MAP(CLeftImageView, CScrollView)
	ON_WM_LBUTTONDOWN()
END_MESSAGE_MAP()


// CLeftImageView drawing

void CLeftImageView::OnInitialUpdate()
/* PURPOSE:		Initializes view during startup
*/
{
	CScrollView::OnInitialUpdate();

	// Default scroll initialization
	CSize sizeTotal;
	sizeTotal.cx = sizeTotal.cy = DEFAULT_SCROLL;
	SetScrollSizes(MM_TEXT, sizeTotal);

	// Set border title
	((CimagealignDoc*)GetDocument())->tempBorderTitle =
		"     Template Image View";
}

void CLeftImageView::OnDraw(CDC* pDC)
/* PURPOSE:		Called when client window needs to be redrawn.

   RECEIVES:	Device context to client window for drawing.

   REMARKS:		Draws image border, paints background black,
				draws current image with edge detected data
				and sampled points, and draws user selected
				rectangles around sampled points.
*/
{
	CRect rect;
	CRect rcClient;

	CimagealignDoc* doc = (CimagealignDoc*)GetDocument();

	// Get client window area in pixel dimensions
	GetClientRect(rcClient);

	// Get system defined colors
	COLORREF crLight = GetSysColor(COLOR_BTNHIGHLIGHT);
	COLORREF crShadow = GetSysColor(COLOR_BTNSHADOW);
	COLORREF crBtnFace = GetSysColor(COLOR_BTNFACE);
	pDC->SetBkMode(TRANSPARENT);
	CGdiObject *pOldFont = pDC->SelectStockObject(ANSI_VAR_FONT);
	rect = rcClient;

	// Draw title border and set title text
	rect.bottom = rect.top + BORDER_HEIGHT;
	rect.right = BORDER_LIMIT;
	pDC->FillSolidRect(rect, crBtnFace);
	pDC->Draw3dRect(rect, crLight, crShadow);
	pDC->DrawText(doc->tempBorderTitle + doc->tempFileName,
		rect, DT_SINGLELINE | DT_VCENTER);

	// Fill client area as black
	rect.top = rect.bottom;
	rect.bottom = BORDER_LIMIT;
	pDC->FillSolidRect(rect, BLK);
	pDC->SelectObject(pOldFont);

	int width = doc->templateWidth;
	int height = doc->templateHeight;

	// Create device compatible bitmap
	CBitmap bitmap;
	bitmap.CreateCompatibleBitmap(pDC, width, height);

	/* Create compatible device context as storage to
	   blit image data from. */
	CDC cdc;
	cdc.CreateCompatibleDC(pDC);

	// Initialize bitmap with template image data and dimensions
	bitmap.SetBitmapBits(width * height * sizeof(DWORD), doc->templateImage);
	CBitmap* oldBitmap = cdc.SelectObject(&bitmap);

	// Display bitmap in client window
	pDC->BitBlt(0, BORDER_HEIGHT, width, height, &cdc, 0, 0, SRCCOPY);

	CBrush brush(CYAN);
	pDC->SelectObject(&brush);

	// Draw rectangle around user selected sampled points
	int i;
	for (i = 0; i < (int)doc->boxes.size(); i++)
	{	Box box = doc->boxes[i];
		IntPr pt1 = box.first;
		IntPr pt2 = box.second;
		CRect cRect(pt1.first, pt1.second, pt2.first, pt2.second);
		pDC->FrameRect(&cRect, &brush);
	}

	pDC->SelectObject(oldBitmap);
}

void CLeftImageView::OnLButtonDown(UINT nFlags, CPoint point)
/* PURPOSE:		Allows user to draw rectangles around sampled
				points in the template view that are then
				mapped over to the test view.  Points
				inside rectangle are then high-lighted
				in feature locate view.

   RECEIVES:	Coordinates of clicked cursor when in client window.
*/
{
	CimagealignDoc* doc = (CimagealignDoc*)GetDocument();
	if (doc->templateImage == 0)
	{	return;
	}

	/* This index keeps track of whether user has clicked twice in
	   client window to draw rectangle. */
	int index = doc->pointCount % 2;

	// Gets current position of scroll bar in pixel dimension
	CPoint scrollPos = GetScrollPosition();

	int width = doc->templateWidth;
	int height = doc->templateHeight;

	/* Add length scrolled to current position of cursor
	   in client window. */
	int posX = point.x + scrollPos.x;
	int posY = point.y + scrollPos.y;

	// If cursor is outside image dimensions, alert user.
	if (posX >= width || posY >= height || posY <= BORDER_HEIGHT)
	{	MessageBox("Cursor outside of image.",
			"Outside image ...", MB_OK | MB_ICONWARNING);
		return;
	}

	// Store user-clicked position in client window
	doc->boxBoundary[index].first = posX;
	doc->boxBoundary[index].second = posY;

	/* This converts point coordinates so that user may draw
	   rectangle in any orientation (e.g. from left to right,
	   right to left, top to bottom, or bottom to top). */
	if (index == 1)
	{	int x1 = doc->boxBoundary[0].first;
		int y1 = doc->boxBoundary[0].second;
		int x2 = doc->boxBoundary[1].first;
		int y2 = doc->boxBoundary[1].second;
		int minX = __min(x1, x2);
		int minY = __min(y1, y2);
		int maxX = __max(x1, x2);
		int maxY = __max(y1, y2);

		// Store coordinates of rectangles in vector
		doc->boxes.push_back(Box(IntPr(minX, minY), IntPr(maxX, maxY)));
		Invalidate();
		UpdateWindow();
	}
	doc->pointCount++;
}


// CLeftImageView diagnostics

#ifdef _DEBUG
void CLeftImageView::AssertValid() const
{
	CScrollView::AssertValid();
}

void CLeftImageView::Dump(CDumpContext& dc) const
{
	CScrollView::Dump(dc);
}
#endif //_DEBUG


// CLeftImageView message handlers





/***************************************************************************

PROJECT:			CS 297 Deliverable 3

FILE:				rightimageview.cpp

PURPOSE:			Class implementation for test image view where user
					may select and display an image that provides
					control points for the template image to be
					mapped to.

COMPILER / TOOL:	Visual C++ .NET

PROGRAMMER:			Wallun Chan

START DATE:			4/1/05

***************************************************************************/

#include "stdafx.h"
#include "imagealign.h"
#include "RightImageView.h"
#include "imagealigndoc.h"


// CRightImageView

IMPLEMENT_DYNCREATE(CRightImageView, CScrollView)

CRightImageView::CRightImageView()
// Default constructor
{
}

CRightImageView::~CRightImageView()
// Destructor
{
}


BEGIN_MESSAGE_MAP(CRightImageView, CScrollView)
END_MESSAGE_MAP()


// CRightImageView drawing

void CRightImageView::OnInitialUpdate()
{
	CScrollView::OnInitialUpdate();

	// Default scroll initialization
	CSize sizeTotal;
	sizeTotal.cx = sizeTotal.cy = DEFAULT_SCROLL;
	SetScrollSizes(MM_TEXT, sizeTotal);

	// Set border title
	((CimagealignDoc*)GetDocument())->testBorderTitle = "     Test Image View";
}

void CRightImageView::OnDraw(CDC* pDC)
{
	CRect rect;
	CRect rcClient;

	CimagealignDoc* doc = (CimagealignDoc*)GetDocument();

	// Get client window area in pixel dimensions
	GetClientRect(rcClient);

	// Get system defined colors
	COLORREF crLight = GetSysColor(COLOR_BTNHIGHLIGHT);
	COLORREF crShadow = GetSysColor(COLOR_BTNSHADOW);
	COLORREF crBtnFace = GetSysColor(COLOR_BTNFACE);
	pDC->SetBkMode(TRANSPARENT);
	CGdiObject *pOldFont = pDC->SelectStockObject(ANSI_VAR_FONT);
	rect = rcClient;

	// Draw title border and set title text
	rect.bottom = rect.top + BORDER_HEIGHT;
	rect.right = BORDER_LIMIT;
	pDC->FillSolidRect(rect, crBtnFace);
	pDC->Draw3dRect(rect, crLight, crShadow);
	pDC->DrawText(doc->testBorderTitle + doc->testFileName,
		rect, DT_SINGLELINE | DT_VCENTER);

	// Fill client area as black
	rect.top = rect.bottom;
	rect.bottom = BORDER_LIMIT;
	pDC->FillSolidRect(rect, BLK);
	pDC->SelectObject(pOldFont);

	int width = doc->testWidth;
	int height = doc->testHeight;

	if (doc->testImage == 0)
		return;

	// Create device compatible bitmap
	CBitmap bitmap;
	bitmap.CreateCompatibleBitmap(pDC, width, height);

	/* Create compatible device context as storage to
	   blit image data from. */
	CDC cdc;
	cdc.CreateCompatibleDC(pDC);

	// Initialize bitmap with test image data and dimensions
	bitmap.SetBitmapBits(width * height * sizeof(DWORD), doc->testImage);

	CBitmap* oldBitmap = cdc.SelectObject(&bitmap);

	// Display bitmap in client window
	pDC->BitBlt(0, BORDER_HEIGHT, width, height, &cdc, 0, 0, SRCCOPY);

	pDC->SelectObject(oldBitmap);
}


// CRightImageView diagnostics

#ifdef _DEBUG
void CRightImageView::AssertValid() const
{
	CScrollView::AssertValid();
}

void CRightImageView::Dump(CDumpContext& dc) const
{
	CScrollView::Dump(dc);
}
#endif //_DEBUG


// CRightImageView message handlers





/***************************************************************************

PROJECT:			CS 297 Deliverable 3

FILE:				featureview.cpp

PURPOSE:			Class implementation for view that displays
					user selected sampled points that are mapped
					from template view to test view.

COMPILER / TOOL:	Visual C++ .NET

PROGRAMMER:			Wallun Chan

START DATE:			4/1/05

***************************************************************************/

#include "stdafx.h"
#include "imagealign.h"
#include "FeatureView.h"
#include "imagealigndoc.h"
#include "utility.h"


// CFeatureView

IMPLEMENT_DYNCREATE(CFeatureView, CScrollView)

CFeatureView::CFeatureView()
// Default constructor
{
}

CFeatureView::~CFeatureView()
// Destructor
{
}


BEGIN_MESSAGE_MAP(CFeatureView, CScrollView)
END_MESSAGE_MAP()


// CFeatureView drawing

void CFeatureView::OnInitialUpdate()
/* PURPOSE:		Initializes view during startup
*/
{
	CScrollView::OnInitialUpdate();

	// Scroll length initialization
	CSize sizeTotal;
	sizeTotal.cx = sizeTotal.cy = DEFAULT_SCROLL;
	SetScrollSizes(MM_TEXT, sizeTotal);
}

void CFeatureView::OnDraw(CDC* pDC)
/* PURPOSE:		Called when client window needs to be redrawn.

   RECEIVES:	Device context to client window for drawing.

   REMARKS:		Draws image border, paints background black,
				draws current image with user selected points
				mapped from template view, and draws convex
				hull around them.
*/
{
	CRect rect;
	CRect rcClient;

	CimagealignDoc* doc = (CimagealignDoc*)GetDocument();

	// Get client window area in pixel dimensions
	GetClientRect(rcClient);

	// Get system defined colors
	COLORREF crLight = GetSysColor(COLOR_BTNHIGHLIGHT);
	COLORREF crShadow = GetSysColor(COLOR_BTNSHADOW);
	COLORREF crBtnFace = GetSysColor(COLOR_BTNFACE);
	pDC->SetBkMode(TRANSPARENT);
	CGdiObject *pOldFont = pDC->SelectStockObject(ANSI_VAR_FONT);
	rect = rcClient;

	// Draw title border and set title text
	rect.bottom = rect.top + BORDER_HEIGHT;
	rect.right = BORDER_LIMIT;
	pDC->FillSolidRect(rect, crBtnFace);
	pDC->Draw3dRect(rect, crLight, crShadow);
	CString borderTitle = "     Feature Locate View - ";
	borderTitle += "( Mapped user-selected points from template view )";
	pDC->DrawText(borderTitle, rect, DT_SINGLELINE | DT_VCENTER);

	// Fill client area as black
	rect.top = rect.bottom;
	rect.bottom = BORDER_LIMIT;
	pDC->FillSolidRect(rect, BLK);
	pDC->SelectObject(pOldFont);

	int width = doc->testWidth;
	int height = doc->testHeight;

	if (doc->featureImage == 0)
		return;

	// Create device compatible bitmap
	CBitmap bitmap;
	bitmap.CreateCompatibleBitmap(pDC, width, height);

	/* Create compatible device context as storage to
	   blit image data from. */
	CDC cdc;
	cdc.CreateCompatibleDC(pDC);

	// Initialize bitmap with feature located image data and dimensions
	bitmap.SetBitmapBits(width * height * sizeof(DWORD), doc->featureImage);
	CBitmap* oldBitmap = cdc.SelectObject(&bitmap);

	// Blit image to client window
	pDC->BitBlt(0, BORDER_HEIGHT, width, height, &cdc, 0, 0, SRCCOPY);

	// Display Cross hairs and convex hulls in white
	pDC->SelectStockObject(WHITE_PEN);

	unsigned i;
	unsigned j;

	// Draw cross-hairs representing centroids of groups of points
	for (i = 0; i < (int)doc->centroids.size(); i++)
	{	IntPr pt = doc->centroids[i];
		pDC->MoveTo(pt.first, pt.second + BORDER_HEIGHT);
		pDC->LineTo(pt.first + CROSS_HAIR_SIZE, pt.second +
			CROSS_HAIR_SIZE + BORDER_HEIGHT);
		pDC->MoveTo(pt.first, pt.second + BORDER_HEIGHT);
		pDC->LineTo(pt.first - CROSS_HAIR_SIZE, pt.second +
			CROSS_HAIR_SIZE + BORDER_HEIGHT);
		pDC->MoveTo(pt.first, pt.second + BORDER_HEIGHT);
		pDC->LineTo(pt.first + CROSS_HAIR_SIZE, pt.second -
			CROSS_HAIR_SIZE + BORDER_HEIGHT);
		pDC->MoveTo(pt.first, pt.second + BORDER_HEIGHT);
		pDC->LineTo(pt.first - CROSS_HAIR_SIZE, pt.second -
			CROSS_HAIR_SIZE + BORDER_HEIGHT);
	}

	if (doc->showGrouping)
	{	// Draw convex hull around user selected sampled points
		for (i = 0; i < doc->ptGroup.size(); i++)
		{	VecIntPr group = doc->ptGroup[i];
			CPoint* pts = new CPoint[group.size()];
			for (j = 0; j < group.size(); j++)
			{	pts[j].x = group[j].first;
				pts[j].y = group[j].second + BORDER_HEIGHT;
			}
			pDC->Polyline(pts, (int)group.size());
			delete [] pts;
		}
	}

	pDC->SelectObject(oldBitmap);
}


// CFeatureView diagnostics

#ifdef _DEBUG
void CFeatureView::AssertValid() const
{
	CScrollView::AssertValid();
}

void CFeatureView::Dump(CDumpContext& dc) const
{
	CScrollView::Dump(dc);
}
#endif //_DEBUG


// CFeatureView message handlers





/***************************************************************************

PROJECT:			CS 297 Deliverable 3

FILE:				mainfrm.h

PURPOSE:			Class definition for frame window to attach
					splitter windows and child views to.

COMPILER / TOOL:	Visual C++ .NET

PROGRAMMER:			Wallun Chan

START DATE:			4/1/05

***************************************************************************/

#pragma once

#include "splitoverride.h"

class CMainFrame : public CFrameWnd
{

protected: // create from serialization only
	CMainFrame();
	DECLARE_DYNCREATE(CMainFrame)

// Attributes
public:
	CSplitOverride m_wndSplitterLR;
	CSplitterWnd m_wndSplitterUD;
	CSplitterWnd m_wndSplitterPIC;

// Operations
public:

// Overrides
public:
	virtual BOOL OnCreateClient(LPCREATESTRUCT lpcs, CCreateContext* pContext);
	virtual BOOL PreCreateWindow(CREATESTRUCT& cs);

// Implementation
public:
	virtual ~CMainFrame();
#ifdef _DEBUG
	virtual void AssertValid() const;
	virtual void Dump(CDumpContext& dc) const;
#endif

protected:  // control bar embedded members
	CStatusBar  m_wndStatusBar;
	CToolBar    m_wndToolBar;

// Generated message map functions
protected:
	afx_msg int OnCreate(LPCREATESTRUCT lpCreateStruct);
	DECLARE_MESSAGE_MAP()
};





/***************************************************************************

PROJECT:			CS 297 Deliverable 3

FILE:				mainfrm.cpp

PURPOSE:			Class implementation for frame window to
					attach splitter windows and child views to.

COMPILER / TOOL:	Visual C++ .NET

PROGRAMMER:			Wallun Chan

START DATE:			4/1/05

***************************************************************************/

#include "stdafx.h"
#include "imagealign.h"
#include "leftimageview.h"
#include "rightimageview.h"
#include "controlview.h"
#include "featureview.h"

#include "MainFrm.h"

#ifdef _DEBUG
#define new DEBUG_NEW
#endif


// CMainFrame

IMPLEMENT_DYNCREATE(CMainFrame, CFrameWnd)

BEGIN_MESSAGE_MAP(CMainFrame, CFrameWnd)
	ON_WM_CREATE()
END_MESSAGE_MAP()

static UINT indicators[] =
{
	ID_SEPARATOR,           // status line indicator
	ID_INDICATOR_CAPS,
	ID_INDICATOR_NUM,
	ID_INDICATOR_SCRL,
};


// CMainFrame construction/destruction

CMainFrame::CMainFrame()
// Default constructor
{
}

CMainFrame::~CMainFrame()
// Destructor
{
}


int CMainFrame::OnCreate(LPCREATESTRUCT lpCreateStruct)
/* PURPOSE:		Initializes and aligns toolbar.
*/
{
	if (CFrameWnd::OnCreate(lpCreateStruct) == -1)
		return -1;

	if (!m_wndToolBar.CreateEx(this, TBSTYLE_FLAT, WS_CHILD | WS_VISIBLE | CBRS_TOP
		| CBRS_GRIPPER | CBRS_TOOLTIPS | CBRS_FLYBY | CBRS_SIZE_DYNAMIC) ||
		!m_wndToolBar.LoadToolBar(IDR_MAINFRAME))
	{
		TRACE0("Failed to create toolbar\n");
		return -1;      // fail to create
	}

	if (!m_wndStatusBar.Create(this) ||
		!m_wndStatusBar.SetIndicators(indicators,
		  sizeof(indicators)/sizeof(UINT)))
	{
		TRACE0("Failed to create status bar\n");
		return -1;      // fail to create
	}
	// TODO: Delete these three lines if you don't want the toolbar to be dockable
	m_wndToolBar.EnableDocking(CBRS_ALIGN_ANY);
	EnableDocking(CBRS_ALIGN_ANY);
	DockControlBar(&m_wndToolBar);
	SetTitle("Image Align / Feature Locate Program");

	return 0;
}

BOOL CMainFrame::OnCreateClient(LPCREATESTRUCT lpcs, CCreateContext* pContext)
/* PURPOSE:		Creates and initializes the splitter windows and child views.
*/
{
	if (!m_wndSplitterLR.CreateStatic(this, 1, 2))
	{	TRACE0("Failed to create left/right splitter window.\n");
		return FALSE;
	}
	if (!m_wndSplitterLR.CreateView(0, 0, RUNTIME_CLASS(CControlView), CSize(), pContext))
	{	TRACE0("Failed to create control view.\n");
		return FALSE;
	}
	if (!m_wndSplitterUD.CreateStatic(&m_wndSplitterLR, 2, 1,
		WS_CHILD | WS_VISIBLE, m_wndSplitterLR.IdFromRowCol(0, 1)))
	{	TRACE0("Failed to create up/down splitter window.\n");
		return FALSE;
	}
	if (!m_wndSplitterPIC.CreateStatic(&m_wndSplitterUD, 1, 2,
		WS_CHILD | WS_VISIBLE, m_wndSplitterUD.IdFromRowCol(0, 0)))
	{	TRACE0("Failed to create picture splitter window.\n");
		return FALSE;
	}
	if (!m_wndSplitterPIC.CreateView(0, 0, RUNTIME_CLASS(CLeftImageView), CSize(), pContext))
	{	TRACE0("Failed to create template view.\n");
		return FALSE;
	}
	if (!m_wndSplitterPIC.CreateView(0, 1, RUNTIME_CLASS(CRightImageView), CSize(), pContext))
	{	TRACE0("Failed to create test view.\n");
		return FALSE;
	}
	if (!m_wndSplitterUD.CreateView(1, 0, RUNTIME_CLASS(CFeatureView), CSize(), pContext))
	{	TRACE0("Failed to create feature view.\n");
		return FALSE;
	}
	m_wndSplitterLR.SetColumnInfo(0, 300, 0);
	m_wndSplitterUD.SetRowInfo(0, 450, 0);
	m_wndSplitterPIC.SetColumnInfo(0, 550, 0);
	return TRUE;
}

BOOL CMainFrame::PreCreateWindow(CREATESTRUCT& cs)
// Windows initialization routine
{
	if( !CFrameWnd::PreCreateWindow(cs) )
		return FALSE;
	// TODO: Modify the Window class or styles here by modifying
	//  the CREATESTRUCT cs

	return TRUE;
}


// CMainFrame diagnostics

#ifdef _DEBUG
void CMainFrame::AssertValid() const
{
	CFrameWnd::AssertValid();
}

void CMainFrame::Dump(CDumpContext& dc) const
{
	CFrameWnd::Dump(dc);
}

#endif //_DEBUG


// CMainFrame message handlers





/***************************************************************************

PROJECT:			CS 297 Deliverable 3

FILE:				utility.h

PURPOSE:			Provides various data structures and function
					declarations to support random number generation,
					mapping points to shape context histograms, storing
					point coordinates, setting up the shape context for
					a set of points defining a shape, obtaining the edges
					from an image via edge detection, overlaying points
					and edges to an image, sampling pixels from edges,
					and parsing a .pgm image file.

COMPILER / TOOL:	Visual C++ .NET

PROGRAMMER:			Wallun Chan

START DATE:			4/1/05

***************************************************************************/

#ifndef _UTILITY_H_
#define _UTILITY_H_

#include "edgefinder.h"
#include "matrx.h"
#include <time.h>

class TPSSolver;

// Really small number
const double SMALL = 1.0E-6;

// Really big number
const double LARGE = 1E9;

// Approximation for pi constant
const double PI = 3.141592653589793238462643383279;

// Size of cross-hairs used in displaying group centroids
const int CROSS_HAIR_SIZE = 7;

// Colors to be used for drawing in device contexts
const COLORREF WHT = RGB(255, 255, 255);
const COLORREF RED = RGB(255, 0, 0);
const COLORREF GRN = RGB(0, 255, 0);
const COLORREF BLUE = RGB(0, 0, 255);
const COLORREF MAG = RGB(255, 0, 255);
const COLORREF CYAN = RGB(0, 255, 255);
const COLORREF YLW = RGB(255, 255, 0);
const COLORREF BLK = RGB(0, 0, 0);
const COLORREF GRAY = RGB(100, 100, 100);

// Border height of gray title bar in each splitter window
const int BORDER_HEIGHT = 25;

/* Used to paint entire client area black for each splitter
   window. */
const int BORDER_LIMIT = 1500;

// Default scroll length
const int DEFAULT_SCROLL = 100;

// Used to size scroll length accordingly
const double SCROLL_FACTOR = 1.25;

// Used to filter for .pgm files in "open" and "save" dialogs
static char pgmFilter[] = "PGM Files (*.pgm)|*.pgm||";
static char bmpFilter[] = "BMP Files (*.bmp)|*.bmp||";

// Integer tuple
typedef pair<int, int> IntPr;

// Tuple to store mapping between a point and its histogram
typedef pair<IntPr, MatInt*> PrMat;

// Used to map a point with its shape context histogram
typedef map<IntPr, MatInt*> Map;

typedef Map::const_iterator MapIter;
typedef vector<int> VecInt;
typedef vector<double> VecD;
typedef vector<IntPr> VecIntPr;

// Used to locate the endpoints of a rectangle
typedef pair<IntPr, IntPr> Box;

// Stores a number of rectangle coordinates
typedef vector<Box> Boxes;

// Stores a number of point groupings
typedef vector<VecIntPr> PtGroup;

typedef pair<IntPr, double> PrDouble;

/* Used to store average vector heading for each point
   on a shape. */
typedef map<IntPr, double> DirectionMap;
typedef DirectionMap::iterator directionIter;

// Random number generator
class Random
{
public:
	Random();
	double randDouble();
	int randRange(int min, int max);
};

IntPr centroid(const VecIntPr& pts);

VecIntPr removeOutliers(const VecIntPr& pts, int threshold = 0);

void minMaxMeanRadDist(const VecIntPr& pts, double& min,
double& max, double& mean);

void setShapeContextPair(Map& tempContext, Map& testContext,
const VecIntPr& tempSampledPts, const VecIntPr& testSampledPts,
int radialBinNum, int angBinNum);

BOOL getSampledPoints(VecIntPr& sampledPts, BYTE* edgeMap,
int width, int height, int sampleNum);

void matchPoints(const Map& tempShape, const Map& testShape,
const VecIntPr& tempPts, VecIntPr& testPts);

VecIntPr getPointList(BYTE* image, int width, int height);

BYTE* getEdgeMap(EdgeDetector* detector, int hysCheck,
int thinCheck, double sigma, int* threshold);

void setDisplayImage(BYTE*& displayImage, BYTE* image, BYTE* edgeMap,
VecIntPr& sampledPts, int width, int height, int imageCheck,
int edgeCheck, int sampleCheck);

DWORD* rgbFormat(BYTE* image, int byteCount);

BYTE* overlay(BYTE* image, int width, int height, VecIntPr& pts,
COLORREF color1, COLORREF color2);

BYTE* rgbSaveFormat(BYTE* image, int width, int height);

BYTE* pgmByteRead(const char* file, int& width, int& height);

void imageOutput(const char* file, BYTE* image, int width, int height);

#endif





/***************************************************************************

PROJECT:			CS 297 Deliverable 3

FILE:				utility.cpp

PURPOSE:			Implements various functions to support random number
					generation, mapping points to shape context histograms,
					storing point coordinates, setting up the shape context
					for a set of points defining a shape, obtaining the edges
					from an image via edge detection, overlaying points and
					edges to an image, sampling pixels from edges, and parsing
					a .pgm image file.

COMPILER / TOOL:	Visual C++ .NET

PROGRAMMER:			Wallun Chan

START DATE:			4/1/05

***************************************************************************/

#include "stdafx.h"
#include "utility.h"
#include "lap.h"
#include "tps.h"

const char* DELIMITERS = " \t\r\n";
const char* P5_FORMAT = "P5";

// Global random number generator
Random random;

Random::Random()
// Initialize random number generator
{
	srand((unsigned)time(0));
}

double Random::randDouble()
/* PURPOSE:		To obtain random double.

   RETURNS:		Random double n in 0.0 < n < 1.0.
*/
{
	return (double)rand() / RAND_MAX;
}

int Random::randRange(int min, int max)
/* PURPOSE:		To obtain number between min and max.

   RECEIVES:	min and max integers specifying range
				of random integer to return

   RETURNS:		Random integer from min to max, exclusive
				of max.
*/
{
	return min + (int)(randDouble() * (max - min));
}

void meanStdDevRadDist(const VecIntPr& pts, const IntPr& center,
double& mean, double& stdDev)
/* PURPOSE:		To calculate the centroid, and mean and standard
				deviation of the radial distances between the
				group centroid and all points.

   RECEIVES:	Points with (x, y) coordinates, centroid coordinate,
				and mean and standard deviation of radial distance
				between group centroid to be calculated.
*/
{
	/* Calculate mean radial distance between centroid and
	   points. */
	mean = 0.0;
	int i;
	int n = (int)pts.size();
	for (i = 0; i < n; i++)
	{	double difX = pts[i].first - center.first;
		double difY = pts[i].second - center.second;
		mean += sqrt(difX * difX + difY * difY);
	}
	mean /= n;

	/* Calculate standard deviation of radial distance
	   between centroid and points. */
	stdDev = 0.0;
	for (i = 0; i < n; i++)
	{	double difX = pts[i].first - center.first;
		double difY = pts[i].second - center.second;
		double dist = sqrt(difX * difX + difY * difY);
		double dif = mean - dist;
		stdDev += (dif * dif);
	}
	stdDev = sqrt(stdDev / (n - 1));
}

IntPr centroid(const VecIntPr& pts)
/* PURPOSE:		Calculate centroid of point set.

   RECEIVES:	Points with (x, y) coordinates.

   RETURNS:		Centroid in (x, y) coordinate.
*/
{
	int sumX = 0;
	int sumY = 0;
	int n = (int)pts.size();
	int i;
	for (i = 0; i < n; i++)
	{	sumX += pts[i].first;
		sumY += pts[i].second;
	}
	sumX = (int)((double)sumX / n);
	sumY = (int)((double)sumY / n);
	return IntPr(sumX, sumY);
}

VecIntPr removeOutliers(const VecIntPr& pts, int threshold)
/* PURPOSE:		Takes a set of points and attempts to remove
				outliers with respect to the mean and standard
				deviation of the radial distances between the
				group centroid and all points.

   RECEIVES:	Points with (x, y) coordinates and factor
				specifying the number of standard deviations
				away from the mean.

   RETURNS:		Original set of points with all outliers with
				radial distances greater than the specified number
				of standard deviations from the mean are removed.
*/
{
	double meanRadDist;
	double stdDevRadDist;
	IntPr center = centroid(pts);
	meanStdDevRadDist(pts, center, meanRadDist, stdDevRadDist);
	double metric = meanRadDist + stdDevRadDist * threshold;
	VecIntPr adjPts;
	int i;
	for (i = 0; i < (int)pts.size(); i++)
	{	double difX = pts[i].first - center.first;
		double difY = pts[i].second - center.second;
		double dist = sqrt(difX * difX + difY * difY);
		if (dist < metric)
		{	adjPts.push_back(pts[i]);
		}
	}
	return adjPts;
}

void minMaxMeanRadDist(const VecIntPr& pts, double& min,
double& max, double& mean)
/* PURPOSE:		Calculates the min, max, and mean radial
				distances between all pairs of points on
				a shape.

   RECEIVES:	Points in (x, y) coordinates, and min, max,
				and mean radial distances to be calculated.
*/
{
	min = LARGE;
	max = 0.0;
	mean = 0.0;
	int i;
	int j;
	int num = (int)pts.size();
	for (i = 0; i < num; i++)
	{	IntPr pt1 = pts[i];
		for (j = 0; j < num; j++)
		{	if (i != j)
			{	IntPr pt2 = pts[j];
				double xDif = pt1.first - pt2.first;
				double yDif = pt1.second - pt2.second;
				double dist = sqrt(xDif * xDif + yDif * yDif);
				if (dist < min)
				{	min = dist;
				}
				if (dist > max)
				{	max = dist;
				}
				mean += dist;
			}
		}
	}
	mean /= (num * (num - 1));
}

void shapeContextParams(const VecIntPr& tempPts, const VecIntPr& testPts,
int radialBinNum, int angBinNum, double& radialBinSize, double& angBinSize,
double& scale, double& tempMean, double& testMean)
/* PURPOSE:		Sets the size of the log-polar histogram
				in the angular and radial directions.

   RECEIVES:	Sampled points for the template and test images, number
				of radial bins, number of angular bins, sizes of the
				radial bin and angular bin to be calculated, scale
				factor to be multiplied with the template image to make
				it the same size as the test image, and mean radial
				distances of the template and test shapes to be
				calculated.
*/
{
	double tempMin;
	double tempMax;
	double testMin;
	double testMax;
	minMaxMeanRadDist(tempPts, tempMin, tempMax, tempMean);
	minMaxMeanRadDist(testPts, testMin, testMax, testMean);

	/* Scale is equal to the reciprocal of the minimum
	   normalized radial distance. */
	scale = 1.0 / __min(tempMin / tempMean, testMin / testMean);

	// Maximum normalized radial distance
	double maxDist = __max(tempMax / tempMean, testMax / testMean);

	// Split logarithmic radial distance into equal parts
	radialBinSize = log(scale * maxDist) / radialBinNum;

	// Size of angular bin
	angBinSize = 360.0 / angBinNum;
}

void directionVectors(const VecIntPr& pts, DirectionMap& dirMap)
/* PURPOSE:		Calculates a rotation invariant axis to measure the
				angular orientation for each sampled point for a
				given shape.

   RECEIVES:	Points with (x, y) coordinates and map object that
				associates a point with its rotation invariant
				angular measurement axis.
*/
{
	int i;
	int j;
	for (i = 0; i < (int)pts.size(); i++)
	{	double sumX = 0.0;
		double sumY = 0.0;
		for (j = 0; j < (int)pts.size(); j++)
		{	if (i != j)
			{	sumX += (pts[i].first - pts[j].first);
				sumY += (pts[i].second - pts[j].second);
			}
		}
		double angle = atan2(sumY, sumX) * 180.0 / PI;
		angle = angle < 0.0 ? 360 + angle : angle;
		dirMap.insert(PrDouble(pts[i], angle));
	}
}

void setShapeContext(Map& shapeContext, const VecIntPr& pts,
int radialBinNum, int angBinNum, double radialBinSize,
double angBinSize, double scale, double mean)
/* PURPOSE:		Calculates the shape context for a shape by
				determining the histogram for each sampled
				point of a shape.

   RECEIVES:	Map object that maps point coordinates with
				the associated shape context histogram of
				a point.
*/
{
	DirectionMap dirMap;
	directionVectors(pts, dirMap);
	int i;
	int j;
	for (i = 0; i < (int)pts.size(); i++)
	{	IntPr pt1 = pts[i];

		// Set size for log-polar histogram
		MatInt* histogram = new MatInt(radialBinNum, angBinNum);
		histogram->Null();

		/* Calculate rotation invariant axis relative to horizontal axis
		   to measure angular orientation of point. */
		directionIter dirIter = dirMap.find(pt1);

		// Sort points in their respective log-polar histograms
		for (j = 0; j < (int)pts.size(); j++)
		{	if (i != j)
			{	IntPr pt2 = pts[j];
				double xDif = pt1.first - pt2.first;
				double yDif = pt1.second - pt2.second;

				// Calculate log radial distance
				double logDist = log(scale * sqrt(xDif * xDif + yDif * yDif) / mean);

				// Sort in radial direction
				int radialBinIndex = (int)(logDist / radialBinSize);
				if (radialBinIndex == radialBinNum)
				{	radialBinIndex--;
				}

				// Calculate angular orientation
				double angle = atan2(yDif, xDif) * 180.0 / PI;
				angle = (angle < 0.0) ? 360.0 + angle : angle;
				angle = (angle > dirIter->second) ? angle -
					dirIter->second : 360.0 - (dirIter->second - angle);

				// Sort in angular direction
				int angBinIndex = (int)(angle / angBinSize);
				if (angBinIndex == angBinNum)
				{	angBinIndex--;
				}

				// Place point in bin
				(*histogram)(radialBinIndex, angBinIndex)++;
			}
		}
		// Map point coordinate with log-polar histogram
		shapeContext.insert(PrMat(pt1, histogram));
	}
}

void setShapeContextPair(Map& tempContext, Map& testContext,
const VecIntPr& tempSampledPts, const VecIntPr& testSampledPts,
int radialBinNum, int angBinNum)
/* PURPOSE:		Determine shape context for all sampled points
				of template and test images.

   RECEIVES:	Map objects for template and test images that
				associates point coordinates with corresponding
				log-polar histogram for each sampled point,
				sampled points of template and test images,
				and number of radial and angular bins.
*/
{
	double tempMean;
	double testMean;
	double scale;
	double radialBinSize;
	double angBinSize;

	// Obtain log-polar histogram sizing data
	shapeContextParams(tempSampledPts, testSampledPts, radialBinNum,
		angBinNum, radialBinSize, angBinSize, scale, tempMean, testMean);

	// Sets shape context for sampled points of template image
	setShapeContext(tempContext, tempSampledPts, radialBinNum,
		angBinNum, radialBinSize, angBinSize, scale, tempMean);

	// Sets shape context for sampled points of test image
	setShapeContext(testContext, testSampledPts, radialBinNum,
		angBinNum, radialBinSize, angBinSize, scale, testMean);
}

BOOL getSampledPoints(VecIntPr& sampledPts, BYTE* edgeMap,
int width, int height, int sampleNum)
/* PURPOSE:		Samples points of an edge detected image.

   RECEIVES:	List to store sampled points, edge detected
				image, width and height of edge detected
				image, and number of points to sample.

   RETURNS:		True if sampling was successful, and false
				otherwise.
*/
{
	VecIntPr tmp;
	int i;
	int j;
	for (i = 0; i < height; i++)
	{	for (j = 0; j < width; j++)
		{	if ((unsigned)edgeMap[i * width + j] > 0)
			{	tmp.push_back(IntPr(j, i));
			}
		}
	}
	if (sampleNum > (int)tmp.size())
	{	return FALSE;
	}
	for (i = 0; i < sampleNum; i++)
	{	int num = random.randRange(0, (int)tmp.size());
		sampledPts.push_back(tmp[num]);

		// Don't erase an otherwise empty vector
		if (!tmp.empty() && num < (int)tmp.size())
		{	tmp.erase(tmp.begin() + num);
		}
	}
	return TRUE;
}

double cost(const Map& tempShape, const Map& testShape,
const IntPr& tempPt, const IntPr& testPt)
/* PURPOSE:		Calculates the chi-squared distance between
				points on two shapes.

   RECEIVES:	Shape context mapping of template and test
				images, and coordinates of points to
				calculate the associated chi-squared
				distance.

   RETURNS:		Cost of assigning a particular template
				and test sampled point pair.
*/
{
	MapIter tempIt = tempShape.find(tempPt);
	MapIter testIt = testShape.find(testPt);
	MatInt* tempHistogram = (*tempIt).second;
	MatInt* testHistogram = (*testIt).second;
	double sum = 0.0;
	int i;
	int j;
	for (i = 0; i < (int)tempHistogram->RowNo(); i++)
	{	for (j = 0; j < (int)tempHistogram->ColNo(); j++)
		{	int tempBin = (*tempHistogram)(i, j);
			int testBin = (*testHistogram)(i, j);
			if (tempBin + testBin != 0)
			{	double subtract = tempBin - testBin;
				double add = tempBin + testBin;
				sum += (subtract * subtract / add);
			}
		}
	}
	return sum / 2.0;
}

void matchPoints(const Map& tempShape, const Map& testShape,
const VecIntPr& tempPts, VecIntPr& testPts)
/* PURPOSE:		Takes sampled point-sets from the edge detected template
				and test images, and based on the calculated shape
				contexts for each set, rearranges the order so that
				a point-to-point correspondence between the sets is
				established.

   RECEIVES:	Shape context mapping of template and test images,
				and sampled point-sets of template and test images.

   REMARKS:		The shape context mapping contains an associated
				histogram for each point coordinate.
*/
{
	int dim = (int)tempPts.size();

	// Set shape context cost matrix
	MatInt mat(dim, dim);
	int i;
	int j;
	for (i = 0; i < dim; i++)
	{	for (j = 0; j < dim; j++)
		{	mat(i, j) = (int)cost(tempShape, testShape, tempPts[i], testPts[j]);
		}
	}
	VecInt rowSol(dim);
	VecInt colSol(dim);
	VecInt u(dim);
	VecInt v(dim);

	/* Call algorithm to solve assignment problem to find
	   best matching between points. */
	lap(dim, mat, rowSol, colSol, u, v);

	/* Order sampled test image points with respect to matching found by
	   Hungarian algorithm. */
	VecIntPr tmp;
	for (i = 0; i < dim; i++)
	{	tmp.push_back(testPts[rowSol[i]]);
	}
	testPts = tmp;
}

VecIntPr getPointList(BYTE* image, int width, int height)
/* PURPOSE:		Takes an edge detected image and retrieves the
				point coordinates of non-zero pixels, and stores
				the coordinates in a vector.

   RECEIVES:	Edge detected image, and width and height of image.

   RETURNS:		List with coordinates of edge detected pixels.
*/
{
	VecIntPr pts;
	int i;
	int j;
	for (i = 0; i < height; i++)
	{	for (j = 0; j < width; j++)
		{	if (image[i * width + j] != 0)
			{	pts.push_back(IntPr(j, i));
			}
		}
	}
	return pts;
}

DWORD* rgbFormat(BYTE* image, int byteCount)
/* PURPOSE:		Takes grayscale image data and converts it into
				a form suitable to be displayed in a compatible
				device context.

   RECEIVES:	Grayscale image data and byte count of image data.

   RETURNS:		Grayscale image converted into a four-byte group
				format.

   REMARKS:		Specifically, the first byte of the four-byte group is
				reserved, the three bytes that follow are for red,
				green, and blue values.
*/
{
	DWORD* bytes = new DWORD[byteCount];
	int i;
	for (i = 0; i < byteCount; i++)
	{	DWORD pixel = (DWORD)image[i];
		bytes[i] = pixel << 16 | pixel << 8 | pixel & 0x00FFFFFF;
	}
	return bytes;
}

DWORD* rgbFormat(BYTE* image, int byteCount, COLORREF color)
/* PURPOSE:		Takes edge detected image data and converts it into
				a form suitable to be displayed by a compatible
				device context.

   RECEIVES:	Edge detected image data, byte count of image data,
				and color of edge data to be displayed.

   RETURNS:		Edge detected image converted into a four-byte group
				format.

   REMARKS:		Specifically, the first byte of the four-byte group is
				reserved, the three bytes that follow are for red,
				green, and blue values.
*/
{
	DWORD red = GetRValue(color);
	DWORD green = GetGValue(color);
	DWORD blue = GetBValue(color);
	DWORD* bytes = new DWORD[byteCount];
	int i;
	for (i = 0; i < byteCount; i++)
	{	if ((unsigned)image[i] > 0)
		{	bytes[i] = red << 16 | green << 8 | blue & 0x00FFFFFF;
		}
		else
		{	bytes[i] = 0;
		}
	}
	return bytes;
}

DWORD* rgbFormat(VecIntPr& pts, int width, int height, COLORREF color)
/* PURPOSE:		Takes a set of points and paints them on a blank image.

   RECEIVES:	Points in (x, y) coordinates, width and height of image
				to display points, and color of displayed points.

   RETURNS:		Image with sampled points displayed on a blank image
				in the specified color.
*/
{
	DWORD red = GetRValue(color);
	DWORD green = GetGValue(color);
	DWORD blue = GetBValue(color);
	DWORD* bytes = new DWORD[width * height];
	int i;
	for (i = 0; i < width * height; i++)
	{	bytes[i] = 0;
	}
	for (i = 0; i < (int)pts.size(); i++)
	{	IntPr pt = pts[i];
		int index = pt.second * width + pt.first;
		if (index >= 0 && index < width * height)
		{	bytes[index] = red << 16 | green << 8 | blue & 0x00FFFFFF;
		}
	}
	return bytes;
}

BYTE* rgbSaveFormat(BYTE* image, int width, int height)
/* PURPOSE:		Converts image data into an rgb triplet format for every
				three bytes suitable for storage as a bitmap file.

   RECEIVES:	Grayscale image intensity data, and width and height
				of image to be saved to a file.

   RETURNS:		Byte array suitable for storage as a bitmap file.
*/
{
	long extWidth = width * 3 + width % 4;
	long imageSize = extWidth * height;
	BYTE* bytes = new BYTE[imageSize];
	DWORD* imDat = (DWORD*)image;
	int i;
	int j;
	for (i = 0; i < height; i++)
	{	for (j = 0; j < width; j++)
		{	DWORD pixel = imDat[i * width + j];
			int pIndex = i * extWidth + j * 3;
			bytes[pIndex] = (BYTE)(pixel & 0x000000FF);
			bytes[pIndex + 1] = (BYTE)((pixel & 0x0000FF00) >> 8);
			bytes[pIndex + 2] = (BYTE)((pixel & 0x00FF0000) >> 16);
		}
	}
	return bytes;
}

BYTE* overlay(BYTE* image, BYTE* edgeMap, int width, int height,
VecIntPr& pts, COLORREF color1, COLORREF color2)
/* PURPOSE:		Overlays original grayscale image, edge detected image,
				and sampled points.

   RECEIVES:	Original grayscale image, edge detected image of grayscale
				image, width and height of images, sampled points to be
				imposed, and color of edge data and sampled points to be
				displayed.

   RETURNS:		Image with original grayscale image, edge detected image,
				and sampled points all super-imposed.
*/
{
	DWORD red = GetRValue(color1);
	DWORD green = GetGValue(color1);
	DWORD blue = GetBValue(color1);
	DWORD* overlay = rgbFormat(image, width * height);
	int i;
	for (i = 0; i < width * height; i++)
	{	if ((unsigned)edgeMap[i] > 0)
		{	overlay[i] = red << 16 | green << 8 | blue & 0x00FFFFFF;
		}
	}
	red = GetRValue(color2);
	green = GetGValue(color2);
	blue = GetBValue(color2);
	for (i = 0; i < (int)pts.size(); i++)
	{	IntPr pt = pts[i];
		int index = pt.second * width + pt.first;
		if (index >= 0 && index < width * height)
		{	overlay[index] =
				red << 16 | green << 8 | blue & 0x00FFFFFF;
		}
	}
	return (BYTE*)overlay;
}

BYTE* overlay(BYTE* image, BYTE* edgeMap, int byteCount, COLORREF color)
/* PURPOSE:		Overlays original grayscale image with edge detected image.

   RECEIVES:	Original grayscale image, edge detected image of grayscale
				image, byte count of images, and color of edge data to be
				displayed.

   RETURNS:		Image with original grayscale image and
				imposed edge detected image.
*/
{
	DWORD red = GetRValue(color);
	DWORD green = GetGValue(color);
	DWORD blue = GetBValue(color);
	DWORD* overlay = rgbFormat(image, byteCount);
	int i;
	for (i = 0; i < byteCount; i++)
	{	if ((unsigned)edgeMap[i] > 0)
		{	overlay[i] = red << 16 | green << 8 | blue & 0x00FFFFFF;
		}
	}
	return (BYTE*)overlay;
}

BYTE* overlay(BYTE* image, int width, int height,
VecIntPr& pts, COLORREF color)
/* PURPOSE:		Overlays original grayscale image with
				sampled points.

   RECEIVES:	Original grayscale image, width and height
				of image, sampled points to be imposed on
				original image, and color of sampled points
				to be displayed.

   RETURNS:		Image with original grayscale image and
				imposed sampled points.
*/
{
	DWORD red = GetRValue(color);
	DWORD green = GetGValue(color);
	DWORD blue = GetBValue(color);
	DWORD* overlay = rgbFormat(image, width * height);
	int i;
	for (i = 0; i < (int)pts.size(); i++)
	{	IntPr pt = pts[i];
		int index = pt.second * width + pt.first;
		if (index >= 0 && index < width * height)
		{	overlay[index] = red << 16 | green << 8 | blue & 0x00FFFFFF;
		}
	}
	return (BYTE*)overlay;
}

BYTE* overlay(BYTE* image, int width, int height, VecIntPr& pts,
COLORREF color1, COLORREF color2)
/* PURPOSE:		Overlays an edge detected image with sampled points.

   RECEIVES:	Edge detected image, width and height of image,
				sampled points to be drawn onto edge detected
				data, and colors of edge data and sampled points
				to be displayed.

   RETURNS:		Image with edge detected data and imposed sampled points.
*/
{
	DWORD red = GetRValue(color1);
	DWORD green = GetGValue(color1);
	DWORD blue = GetBValue(color1);
	DWORD* overlay = rgbFormat(image, width * height);
	int i;
	for (i = 0; i < width * height; i++)
	{	if ((unsigned)image[i] > 0)
		{	overlay[i] = red << 16 | green << 8 | blue & 0x00FFFFFF;
		}
	}
	red = GetRValue(color2);
	green = GetGValue(color2);
	blue = GetBValue(color2);
	for (i = 0; i < (int)pts.size(); i++)
	{	IntPr pt = pts[i];
		int index = pt.second * width + pt.first;
		if (index >= 0 && index < width * height)
		{	overlay[index] = red << 16 | green << 8 | blue & 0x00FFFFFF;
		}
	}
	return (BYTE*)overlay;
}

BYTE* getEdgeMap(EdgeDetector* detector, int hysCheck,
int thinCheck, double sigma, int* threshold)
/* PURPOSE:		(Re)generates image data for user-selections
				of hysteresis and thinning, and user inputs
				for sigma and threshold.

   RECEIVES:	Edge detector containing image edge data,
				user-selections of hysteresis and thinning,
				and values for sigma (increases the amount
				of smoothing) and threshold (sensitivity
				to intensity changes).

   RETURNS:		Generated image based on user-selections of
				hysteresis and thin checkboxes, and user
				inputs for sigma and threshold.

   REMARKS:		Hysteresis is post-postprocessing of edge
				detection where image pixels that should
				be and have not been suppressed are done
				so. Thinning is as it sounds, that is,
				edges are thinned to its minimum, i.e.
				single pixel width edges.
*/
{
	BYTE* edgeMap;
	if (hysCheck && thinCheck)
	{
		detector->image_find_edges(sigma, threshold);
		detector->image_hysteresis(*threshold / 3);
		edgeMap = detector->image_thin();
	}
	else if (hysCheck)
	{	detector->image_find_edges(sigma, threshold);
		edgeMap = detector->image_hysteresis(*threshold / 3);
	}
	else if (thinCheck)
	{	detector->image_find_edges(sigma, threshold);
		edgeMap = detector->image_thin();
	}
	else
	{	edgeMap = detector->image_find_edges(sigma, threshold);
	}
	return edgeMap;
}

void setDisplayImage(BYTE*& displayImage, BYTE* image, BYTE* edgeMap,
VecIntPr& sampledPts, int width, int height, int imageCheck,
int edgeCheck, int sampleCheck)
/* PURPOSE:		Depending on user-selections of the checkboxes,
				this function generates the appropriate image
				to be displayed.

   RECEIVES:	Image to generate and display,
*/
{
	// Overlay original image, edge detected image, and sampled points
	if (imageCheck && edgeCheck && sampleCheck)
	{	displayImage = overlay(image, edgeMap, width, height,
			sampledPts, BLUE, YLW);
	}
	// Overlay original and edge detected image
	else if (imageCheck && edgeCheck)
	{	displayImage = overlay(image, edgeMap, width * height, BLUE);
	}
	// Overly edge detected image and sampled points
	else if (edgeCheck && sampleCheck)
	{	displayImage = overlay(edgeMap, width, height,
			sampledPts, BLUE, YLW);
	}
	// Overly original image and sampled points
	else if (imageCheck && sampleCheck)
	{	displayImage = overlay(image, width, height, sampledPts, YLW);
	}
	// Show only original image
	else if (imageCheck)
	{	displayImage = (BYTE*)rgbFormat(image, width * height);
	}
	// Show only edge data
	else if (edgeCheck)
	{	displayImage = (BYTE*)rgbFormat(edgeMap, width * height, BLUE);
	}
	// Show only sampled points
	else if (sampleCheck)
	{	displayImage = (BYTE*)rgbFormat(sampledPts, width, height, YLW);
	}
	// Show nothing
	else
	{	displayImage = 0;
	}
}

BYTE* pgmByteRead(const char* file, int& width, int& height)
/* PURPOSE:		Parses a .pgm image file for image dimensions and pixel
				intensities, and stores values in unsigned char array.

   RECEIVES:	Image file path name, unsigned char array to store image
				values to, and image pixel width and height.

   RETURNS:		unsigned char array with initialized values or 0.

   REMARKS:		1) This function reads both the "P2" and "P5" version
				   of the .pgm format; for P2, each pixel intensity
				   value is represented as an ASCII decimal number
				   string (e.g. 101).  For P5, each pixel is
				   represented in 8 or 16-bit binary.

				2) The number of grayscale bits is 8, giving intensity
				   values from 0 to 255.
*/
{
	int curPos = 0;
	int index = 0;
	int linePos = 0;
	BYTE* bytes;
	CString line;
	CString token;
	CString format;
	CStdioFile fileHandle;
	if (!fileHandle.Open(file, CFile::modeRead | CFile::typeBinary))
	{	return 0;
	}
	while(fileHandle.ReadString(line))
	{	// Ignore comments marked by a '#'
		if (line.Find('#') == -1)
		{	// Read .pgm type (P2 or P5)
			if (linePos == 0)
			{	format = line.Tokenize(DELIMITERS, curPos);
			}
			// Read image dimensions
			else if (linePos == 1)
			{	curPos = 0;
				token = line.Tokenize(DELIMITERS, curPos);
				width = atoi(token);
				token = line.Tokenize(DELIMITERS, curPos);
				height = atoi(token);
				bytes = new BYTE[width * height];
			}
			// Check if it's P5 format, since P2 is default .pgm format
			else if (linePos == 2)
			{	// Read pixel intensities for P5 format
				if (format == P5_FORMAT)
				{	BYTE* longText = new BYTE[width * height];
					int numRead = fileHandle.Read(longText, width * height);
					for (index = 0; index < numRead; index++)
					{	bytes[index] = longText[index];
					}
					delete [] longText;
					break;
				}
			}
			// Read pixel intensities for P2 format
			else
			{	curPos = 0;
				token = line.Tokenize(DELIMITERS, curPos);
				while (token != "")
				{	bytes[index] = atoi(token);
					index++;
					token = line.Tokenize(DELIMITERS, curPos);
				}
			}
			linePos++;
		}
	}
	fileHandle.Close();
	return bytes;
}

void imageOutput(const char* file, BYTE* image, int width, int height)
/* PURPOSE:		Saves a bitmap image

   RECEIVES:	File name of bitmap file, pixel data of image, and pixel
				dimensions of image.

   REMARKS:		This is a rather crude and inflexible way of saving
				exclusively .bmp files.  But I wanted to save
				images with the least amount of code, and this
				is probably it.
*/
{
	BITMAPINFOHEADER bmih;
	bmih.biSize = sizeof(BITMAPINFOHEADER);
	bmih.biBitCount = 24;
	bmih.biPlanes = 1;
	bmih.biCompression = BI_RGB;
	bmih.biWidth = width;
	bmih.biHeight = -height;
	bmih.biSizeImage = (width * 3 + width % 4) * height;

    BITMAPFILEHEADER bmfh;
    int nBitsOffset = sizeof(BITMAPFILEHEADER) + bmih.biSize;
    LONG imSize = bmih.biSizeImage;
    LONG lFileSize = nBitsOffset + imSize;
    bmfh.bfType = 'B' + ('M' << 8);
    bmfh.bfOffBits = nBitsOffset;
    bmfh.bfSize = lFileSize;
    bmfh.bfReserved1 = 0;
	bmfh.bfReserved2 = 0;

	FILE* pFile = fopen(file, "wb");
    fwrite(&bmfh, 1, sizeof(BITMAPFILEHEADER), pFile);
    fwrite(&bmih, 1, sizeof(BITMAPINFOHEADER), pFile);
	fwrite(rgbSaveFormat(image, width, height), 1, imSize, pFile);
    fclose(pFile);
}





/***************************************************************************

PROJECT:			CS 297 Deliverable 3

FILE:				edgefinder.h

PURPOSE:			Function declarations for Boie-Cox edge detector to
					perform edge detection on a grayscale image.

COMPILER / TOOL:	Visual C++ .NET

PROGRAMMER:			Wallun Chan

START DATE:			4/1/05

***************************************************************************/

#ifndef _EDGE_FINDER_H_
#define _EDGE_FINDER_H_

const int EDGE_X = 252;
const int EDGE_Y = 253;
const int EDGE_45 = 254;
const int EDGE_135 = 255;
const int FWD = 0;
const int BWD = 1;
const double diag_scale = 1.414;

class EdgeDetector
{
public:
	EdgeDetector(unsigned char* array = 0, int xdim = 0, int ydim = 0);
	~EdgeDetector();
	void freeImage();
	void setImage(unsigned char* array, int xdim, int ydim);
	unsigned char* image_find_edges(double sigma, int* threshold);
	unsigned char* image_hysteresis(int lo_thld);
	unsigned char* image_thin();
	unsigned char* get_image() const { return an_image; }
	unsigned char* get_edge_map() const { return edge_map; }
	unsigned char* get_edge_map_high() const { return edge_map_hi; }
	int get_x() const { return nx; }
	int get_y() const { return ny; }
protected:
	void image_detx();
	void image_dety();
	void image_det45();
	void image_det135();
	void generate_filter(double sigma);
	void image_convolve(unsigned char* pic);
	void image_neighbor_fwd(int iy, int ix, unsigned char* mapP);
	void image_neighbor_bwd(int iy, int ix, unsigned char* mapP);
	unsigned char* image_edges(int threshold);
	int image_locx(int iy, int ix, int index);
	int image_locy(int iy, int ix, int index);
	int image_loc45(int iy, int ix, int index);
	int image_loc135(int iy, int ix, int index);
	int image_zc_x(int ix, int iy, int i);
	int image_zc_y(int ix, int iy, int i);
	int image_zc_135(int ix, int iy, int i);
	int image_zc_45(int ix, int iy, int i);
	int image_edge_map_lo(int ix, int iy, int index);
	int Fcontour(int iy, int ix);
	int Bcontour(int iy, int ix);
	int image_threshold();
private:
	unsigned char* edge_map;
	unsigned char* edge_map_hi;
	unsigned char* an_image;
	int nx;
	int ny;
	int norm;
	int nf;
	int nf2;
	int orient_flag;
	int lo_threshold;
	int* gaussian;
	int* idx;
	int* idy;
	int* id45;
	int* id135;
	int* filter;
};

#endif





/***************************************************************************

PROJECT:			CS 297 Deliverable 3

FILE:				edgefinder.cpp

PURPOSE:			Implementation for Boie-Cox edge detector to
					perform edge detection on a grayscale image.

COMPILER / TOOL:	Visual C++ .NET

PROGRAMMER:			Wallun Chan

START DATE:			4/1/05

REMARKS:			This code was borrowed and adapted from
					http://www.ee.ucl.ac.uk/~icox/#Anchor--Ima-25529.
					The code was originally written in C, but was
					converted into slightly more object-oriented
					C++ code.

					2) The code only works for gray scale images.

***************************************************************************/

#include "stdafx.h"
#include "edgefinder.h"
#include <math.h>

EdgeDetector::EdgeDetector(BYTE* array, int xdim, int ydim)
/* PURPOSE:		Default constructor to initialize detector
				member variables.

   RECEIVES:	Image data in byte array, width, and
				height of image.
*/
{
	an_image = array;
	gaussian = 0;
	filter = 0;
	idx = 0;
	idy = 0;
	id45 = 0;
	id135 = 0;
	edge_map = 0;
	norm = 0;
	nf = 0;
	nf2 = 0;
	orient_flag = 0;
	lo_threshold = 0;
	edge_map_hi = 0;
	nx = xdim;
	ny = ydim;
}

EdgeDetector::~EdgeDetector()
// Destructor
{
	freeImage();
}

void EdgeDetector::freeImage()
// Release allocated memory
{
	delete [] gaussian;
	delete [] filter;
	delete [] idx;
	delete [] idy;
	delete [] id45;
	delete [] id135;
	delete [] edge_map_hi;
	delete [] an_image;
	gaussian = 0;
	filter = 0;
	idx = 0;
	idy = 0;
	id45 = 0;
	id135 = 0;
	norm = 0;
	nf = 0;
	nf2 = 0;
	orient_flag = 0;
	lo_threshold = 0;
	nx = 0;
	ny = 0;
	edge_map = 0;
	edge_map_hi = 0;
	an_image = 0;
}

void EdgeDetector::setImage(BYTE* array, int xdim, int ydim)
/* PURPOSE:		Initializes detector with image data to be
				edge detected.

   RECEIVES:	Image data in byte array, width, and
				height of image.
*/
{
	freeImage();
	an_image = array;
	nx = xdim;
	ny = ydim;
}

BYTE* EdgeDetector::image_find_edges(double sigma, int* threshold)
/* PURPOSE:		Generates edge data for image.

   RECEIVES:	Double value for sigma that determines the width
				of the Gaussian filter (increasing this value
				provides more smoothing, which may decrease
				edge detail), and integer value for the threshold,
				which determines how sensitive the detector is
				to changes in intensity (increasing this
				value decreases edge content).

   RETURNS:		Edge data of image.

   REMARKS:		If a zero is input for threshold, then a best
				effort threshold is estimated.
*/
{
	generate_filter(sigma);
	image_convolve(an_image);
	image_detx();
	image_dety();
	image_det45();
	image_det135();
	if (*threshold == 0)
	{	*threshold = image_threshold();
	}
	edge_map = image_edges(*threshold);
	return edge_map;
}

void EdgeDetector::generate_filter(double sigma)
{
	int i;
	double two_sigma_sq;
	nf = (int)(6 * sigma);
	nf = (nf % 2 > 0 ? nf : nf + 1);
	nf2 = nf / 2;
	filter = new int[nf];
	two_sigma_sq = (double)(2.0 * sigma * sigma);
	norm = 0;
	for (i = 0; i < nf; i++)
	{	double x = (double)i - nf2;
		filter[i] = (int)(255.0 * exp(-x * x / two_sigma_sq));
		norm += filter[i];
	}
	norm *= norm;
}

void EdgeDetector::image_convolve(BYTE* pic)
{
	int ik;
	int ix;
	int iy;
	int nxX2M1;
	int iklim;
	int ikXnx;
	int nfXnx;
	int nyX2M1Xnx;
	int iyMnf2Xnx;
	int iy0;
	int iy0Xnx;
	int it;
	int result;
	int* filtP;
	int* tmpP;
	int* gaussianP;
	int* temp_image;
	BYTE* tpicP;
	BYTE* picP;
	nfXnx = nf * nx;
	nxX2M1 = nx + nx - 1;
	nyX2M1Xnx = (ny + ny - 1) * nx;
	temp_image = new int[nx * ny];
	gaussian = new int[nx * ny];
	tmpP = &temp_image[0];
	picP = pic;
	for (iy = 0; iy < ny; iy++, tmpP += nx, picP += nx)
	{	for (ix = 0; ix < nf2; ix++)
		{	filtP = filter;
			result = 0;
			ik = ix - nf2;
			iklim = ik + nf;
			for (; ik < iklim; ik++)
			{	result += (*filtP++) * picP[abs(ik)];
			}
			*(tmpP + ix) = result;
		}
	}
	tmpP = &temp_image[0];
	tpicP = pic;
	for (iy = 0; iy < ny; iy++, tmpP += nx, tpicP += nx)
	{	int xlim = nx - nf2 - 1;
		for (ix = nf2; ix < xlim; ix++)
		{	picP = tpicP + ix - nf2;
			filtP = filter;
			result = 0;
			ik = nf;
			while (ik--)
			{	result += (*filtP++) * (*picP++);
			}
			*(tmpP + ix) = result;
		}
	}
	tmpP = &temp_image[0];
	picP = pic;
	for (iy = 0; iy < ny; iy++, tmpP += nx, picP += nx)
	{	for (ix = nx - nf2 - 1; ix < nx; ix++)
		{	filtP = filter;
			result = 0;
			ik = ix - nf2;
			iklim = ik + nf;
			for (; ik < iklim; ik++)
			{	if (ik >= nx)
				{	result += *filtP * picP[nxX2M1 - ik];
				}
				else
				{	result += *filtP * picP[ik];
				}
				filtP++;
			}
			*(tmpP + ix) = result;
		}
	}
	it = -nf2 * nx;
	for (ix = 0; ix < nx; ix++)
	{	gaussianP = &gaussian[ix];
		tmpP = &temp_image[ix];
		iy = 0;
		iyMnf2Xnx = it;
		for (; iy < nf2; iy++, iyMnf2Xnx += nx)
		{	filtP = filter;
			result = 0;
			ik = iyMnf2Xnx;
			iklim = ik + nfXnx;
			for (; ik < iklim; ik += nx)
			{	result += (*filtP++) * tmpP[abs (ik)];
			}
			*(gaussianP) = result / norm;
			gaussianP += nx;
		}
	}
	iy0 = nf2;
	iy0Xnx = iy0 * nx;
	for (ix = 0; ix < nx; ix++)
	{	int* col_tmpP = &temp_image[ix];
		int ylim = ny - nf2 - 1;
		iy = iy0;
		gaussianP = &gaussian[iy0Xnx + ix];
		iyMnf2Xnx = 0;
		for (; iy < ylim; iy++, iyMnf2Xnx += nx)
		{	filtP = filter;
			tmpP = col_tmpP + iyMnf2Xnx;
			result = 0;
			ik = nf;
			while (ik--)
			{	result += (*filtP++) * (*tmpP);
				tmpP += nx;
			}
			*(gaussianP) = result / norm;
			gaussianP += nx;
		}
	}
	iy0 = ny - nf2 - 1;
	iy0Xnx = iy0 * nx;
	it = (iy0 - nf2) * nx;
	for (ix = 0; ix < nx; ix++)
	{	tmpP = &temp_image[ix];
		iy = iy0;
		gaussianP = &gaussian[iy0Xnx + ix];
		iyMnf2Xnx = it;
		for (; iy < ny; iy++, iyMnf2Xnx += nx)
		{	filtP = filter;
			result = 0;
			ik = iy - nf2;
			iklim = ik + nf;
			ikXnx = iyMnf2Xnx;
			for (; ik < iklim; ik++)
			{	if (ik >= ny)
				{	result += *filtP * tmpP[nyX2M1Xnx - ikXnx];
				}
				else
				{	result += *filtP * tmpP[ikXnx];
				}
				filtP++;
				ikXnx += nx;
			}
			*(gaussianP) = result / norm;
			gaussianP += nx;
		}
	}
	delete [] temp_image;
}

int EdgeDetector::image_zc_x(int ix, int iy, int i)
{
	return ((image_locx(iy, ix, i) > 0) !=
		(image_locx(iy, ix + 1, i + 1) > 0)) ? EDGE_X : 0;
}

int EdgeDetector::image_zc_y(int ix, int iy, int i)
{
	return ((image_locy(iy + 1, ix, i + nx) > 0) !=
		(image_locy(iy, ix, i) > 0)) ? EDGE_Y : 0;
}

int EdgeDetector::image_zc_135(int ix, int iy, int i)
{
	return ((image_loc135(iy + 1, ix + 1, i + nx + 1) > 0) !=
		(image_loc135(iy, ix, i) > 0)) ? EDGE_135 : 0;
}

int EdgeDetector::image_zc_45(int ix, int iy, int i)
{
	return ((image_loc45(iy + 1,ix - 1, i + nx - 1) > 0) !=
		(image_loc45(iy, ix, i) > 0)) ? EDGE_45 : 0;
}

unsigned char* EdgeDetector::image_edges(int threshold)
{
	int ix;
	int iy;
	int orient_flag;
	int index;
	int* px = idx;
	int* py = idy;
	int* p45 = id45;
	int* p135 = id135;
	BYTE* mp;
	BYTE* mapp;
	mp = new BYTE[nx * ny];
	mapp = mp;
	for (iy = 0; iy < ny; iy++)
	{	for (ix = 0; ix < nx; ix++)
		{	int sx;
			int sy;
			int s45;
			int s135;
			int big;
			sx = abs(*px);
			px++;
			sy = abs(*py);
			py++;
			s45 = (int)((double)(abs(*p45)) / diag_scale);
			p45++;
			s135 = (int)((double)(abs(*p135)) / diag_scale);
			p135++;
			big = sx;
			orient_flag = EDGE_X;
			if (sy > big)
			{	big = sy;
				orient_flag = EDGE_Y;
			}
			if (s45 > big)
			{	big = s45;
				orient_flag = EDGE_45;
			}
			if (s135 > big)
			{	big = s135;
				orient_flag = EDGE_135;
			}
			orient_flag = orient_flag;
			if (big > threshold)
			{	index = (int)(mapp - mp);
				switch (orient_flag)
				{
					case EDGE_X:
						*mapp++ = image_zc_x(ix, iy, index);
						break;
					case EDGE_Y:
						*mapp++ = image_zc_y(ix, iy, index);
						break;
					case EDGE_45:
						*mapp++ = image_zc_45(ix, iy, index);
						break;
					case EDGE_135:
						*mapp++ = image_zc_135(ix, iy, index);
						break;
					default:
						throw "EdgeDetector::image_edges() - error in case statements.";
				}
			}
			else
			{	*mapp++ = 0;
			}
		}
	}
	return mp;
}

int EdgeDetector::image_edge_map_lo(int ix, int iy, int index)
{
	int big;
	int sx;
	int sy;
	int s45;
	int s135;
	sx = abs(idx[index]);
	sy = abs(idy[index]);
	s45 = (int)((double)(abs(id45[index])) / diag_scale);
	s135 = (int)((double)(abs(id135[index])) / diag_scale);
	big = sx;
	orient_flag = EDGE_X;
	if (sy > big)
	{	big = sy;
		orient_flag = EDGE_Y;
	}
	if (s45 > big)
	{	big = s45;
		orient_flag = EDGE_45;
	}
	if (s135 > big)
	{	big = s135;
		orient_flag = EDGE_135;
	}
	if (big > lo_threshold)
	{	switch (orient_flag)
		{
			case EDGE_X:
				return image_zc_x(ix, iy, index);
			case EDGE_Y:
				return image_zc_y(ix, iy, index);
			case EDGE_45:
				return image_zc_45(ix, iy, index);
			case EDGE_135:
				return image_zc_135(ix, iy, index);
			default:
				throw "EdgeDetector::image_edge_map_lo() - error in case statements.";
		}
	}
	else
	{	return 0;
	}
}

void EdgeDetector::image_detx()
{
	int* dP;
	int* gP;
	int ix;
	int iy;
	int nxM1;
	nxM1 = nx - 1;
	idx = new int[nx * ny];
	dP = idx;
	gP = gaussian;
	for (iy = 0; iy < ny; iy++)
	{	*dP = *(gP + 1) - *(gP);
		dP += nx;
		gP += nx;
	}
	for (iy = 0; iy < ny; iy++)
	{	dP = idx + iy * nx + 1;
		gP = gaussian + iy * nx + 2;
		for (ix = 1; ix < nxM1; ix++)
		{	*dP = *gP - *(gP - 2);
			dP++;
			gP++;
		}
	}
	dP = idx + nxM1;
	gP = gaussian + nxM1;
	for (iy = 0; iy < ny; iy++)
	{	*dP = *gP - *(gP - 1);
		dP += nx;
		gP += nx;
	}
}

void EdgeDetector::image_dety()
{
	int* dP;
	int* gP;
	int ix;
	int iy;
	int nyM1;
	int nxPnx;
	nxPnx = nx + nx;
	nyM1 = ny - 1;
	idy = new int[nx * ny];
	dP = idy;
	gP = gaussian;
	for (ix = 0; ix < nx; ix++)
	{	dP[ix] = gP[ix + nx] - gP[ix];
	}
	for (ix = 0; ix < nx; ix++)
	{	dP = idy + ix;
		gP = gaussian + ix;
		for (iy = 1; iy < nyM1; iy++)
		{	dP += nx;
			*dP = *(gP + nxPnx) - *gP;
			gP += nx;
		}
	}
	dP = idy + (ny - 1) * nx;
	gP = gaussian + (ny - 1) * nx;
	for (ix = 0; ix < nx; ix++)
	{	dP[ix] = gP[ix] - gP[ix - nx];
	}
}

void EdgeDetector::image_det45()
{
	int* dP;
	int* gP;
	int i;
	int ix;
	int iy;
	int index;
	int nxP1;
	int nxM1;
	int nyM1;
	nxM1 = nx - 1;
	nxP1 = nx + 1;
	nyM1 = ny - 1;
	id45 = new int[nx * ny];
	dP = id45 + nx - 2;
	gP = gaussian + nx - 2;
	for (ix = nx - 2; ix > 0; ix--)
	{	*dP = *(gP + nxM1) - *(gP + 1);
		dP--;
		gP--;
	}
	dP = id45;
	gP = gaussian;
	for (i = nx - 2; i > 0; i--)
	{	ix = i;
		iy = 1;
		index = nx + ix;
		while (ix > -1 && iy < nyM1)
		{	dP[index] = gP[index + nxM1] - gP[index - nxM1];
			index += nxM1;
			ix--;
			iy++;
		}
	}
	for (i = 1; i < ny; i++)
	{	ix = nx - 2;
		iy = i;
		index = iy * nx + ix;
		while (ix > 0 && iy < nyM1)
		{	dP[index] = gP[index + nxM1] - gP[index - nxM1];
			index += nxM1;
			ix--;
			iy++;
		}
	}
	dP = id45 + nx + nx;
	gP = gaussian + nx + nx;
	for (iy = 1; iy < nyM1; iy++)
	{	*dP = *(gP + nx) - *(gP - nxM1);
	}
	dP = id45 + nx + nx - 1;
	gP = gaussian + nx + nx - 1;
	for (iy = 1; iy < nyM1; iy++)
	{	*dP = *(gP + nxM1) - *(gP - nx);
		dP += nx;
		gP += nx;
	}
	dP = id45 + (ny - 1) * nx;
	gP = gaussian + (ny - 1) * nx;
	for (ix = 1; ix < nxM1; ix++)
	{	*dP = *(gP - 1) - *(gP - nxM1);
		dP++;
		gP++;
	}
	dP = id45;
	gP = gaussian;
	*dP = *(gP + nx) - *(gP + 1);
	dP = id45 + (ny - 1) * nx;
	gP = gaussian + (ny - 1) * nx;
	*dP = *(gP) - *(gP - nxM1);
	dP[nxM1] = *(gP + nxM1 - 1) - *(gP - 1);
	dP = id45 + nxM1;
	gP = gaussian + nxM1;
	*dP = *(gP + nxM1) - *(gP);
}

void EdgeDetector::image_det135()
{
	int* dP;
	int* gP;
	int i;
	int ix;
	int iy;
	int index;
	int nxP1;
	int nxM1;
	int nyM1;
	int nxM2;
	nxM1 = nx - 1;
	nxP1 = nx + 1;
	nxM2 = nx - 2;
	nyM1 = ny - 1;
	id135 = new int[nx * ny];
	dP = id135 + 1;
	gP = gaussian + 1;
	for (ix = 1; ix < nxM1; ix++)
	{	*dP = *(gP + nxP1) - *(gP - 1);
		dP++;
		gP++;
	}
	dP = id135;
	gP = gaussian;
	for (i = 1; i < nxM1; i++)
	{	ix = i;
		iy = 1;
		index = nx + ix;
		while (ix < nxM1 && iy < nyM1)
		{	dP[index] = gP[index + nxP1] - gP[index - nxP1];
			index += nxP1;
			ix++;
			iy++;
		}
	}
	for (i = 1; i < ny; i++)
	{	ix = 1;
		iy = i;
		index = iy * nx + ix;
		while (ix < nxM1 && iy < nyM1)
		{	dP[index] = gP[index + nxP1] - gP[index - nxP1];
			index += nxP1;
			ix++;
			iy++;
		}
	}
	dP = id135 + nx;
	gP = gaussian + nx;
	for (iy = 1; iy < nyM1; iy++)
	{	*dP = *(gP + nxP1) - *(gP - nx);
		dP += nx;
		gP += nx;
	}
	dP = id135 + nx + nx - 1;
	gP = gaussian + nx + nx - 1;
	for (iy = 1; iy < nyM1; iy++)
	{	*dP = *(gP + nx) - *(gP - nxP1);
		dP += nx;
		gP += nx;
	}
	dP = id135 + (ny - 1) * nx + 1;
	gP = gaussian + (ny - 1) * nx + 1;
	for (ix = 1; ix < nxM1; ix++)
	{	*dP = *(gP + 1) - *(gP - nxP1);
		dP++;
		gP++;
	}
	dP = id45;
	gP = gaussian;
	*dP = *(gP + nxP1) - *(gP);
	dP = id45 + (ny - 1) * nx;
	gP = gaussian + (ny - 1) * nx;
	*dP = *(gP + 1) - *(gP - nx);
	dP[nxM1] = *(gP + nxM1) - *(gP - 2);
	dP = id45 + nxM1;
	gP = gaussian + nxM1;
	*dP = *(gP + nx) - *(gP - 1);
}

BYTE* EdgeDetector::image_hysteresis(int lo_thld)
{
	int iy;
	int ix;
	BYTE* lmapP;
	BYTE* lmap_hiP;
	BYTE* mapP;
	lo_threshold = lo_thld;
	edge_map_hi = edge_map;
	edge_map = new BYTE[nx * ny];
	lmap_hiP = &edge_map_hi[0];
	lmapP = &edge_map[0];
	mapP = lmapP;
	for (ix = nx * ny; --ix >= 0;)
	{	*mapP++ = 0;
	}
	for (iy = 0; iy < ny; iy++, lmap_hiP += nx, lmapP += nx)
	{	mapP = lmapP;
		for (ix = 0; ix < nx; ix++, mapP++)
		{	if (lmap_hiP[ix] != 0)
			{	*mapP = lmap_hiP[ix];
				image_neighbor_fwd(iy, ix, mapP);
				image_neighbor_bwd(iy, ix, mapP);
			}
		}
	}
	return edge_map;
}

unsigned char* EdgeDetector::image_thin()
{
	int ix;
	int iy;
	int nxM1;
	int nxP1;
	int nlinks;
	int npieces;
	int nyM1;
	int b01;
	int b12;
	int b21;
	int b10;
	int p1;
	int p2;
	int p3;
	int p4;
	int b00;
	int b02;
	int b20;
	int b22;
	BYTE* emP;
	BYTE* emPiy;
	nxM1 = nx - 1;
	nxP1 = nx + 1;
	nyM1 = ny - 1;
	emPiy = &edge_map[0];
	for (iy = 1; iy < nyM1; iy++)
	{	emPiy += nx;
		for (ix = 1; ix < nxM1; ix++)
		{	emP = emPiy + ix;
			b01 = *(emP - nx) > 0;
			b12 = *(emP + 1) > 0;
			b21 = *(emP + nx) > 0;
			b10 = *(emP - 1) > 0;
			if ((b01 + b12 + b21 + b10) > 1)
			{	b00 = *(emP - nxP1) > 0;
				b02 = *(emP - nxM1) > 0;
				b20 = *(emP + nxM1) > 0;
				b22 = *(emP + nxP1) > 0;
				p1 = b00 | b01;
				p2 = b02 | b12;
				p3 = b22 | b21;
				p4 = b20 | b10;
				nlinks = b01 & p2;
				nlinks += b12 & p3;
				nlinks += b21 & p4;
				nlinks += b10 & p1;
				npieces = p1 + p2 + p3 + p4;
				if ((npieces - nlinks) < 2)
				{	*emP = 0;
				}
			}
		}
	}
	return edge_map;
}

void EdgeDetector::image_neighbor_fwd(int iy, int ix, BYTE* mapP)
{
	switch (*mapP)
	{
		case 5:
			return;
		case EDGE_Y:
			if ((!Fcontour(iy, ix + 1) && Fcontour(iy, ix + 2)) ||
				(!(Fcontour(iy - 1, ix + 1) || Fcontour(iy + 1, ix + 1)) &&
				(Fcontour(iy - 1, ix + 2) || Fcontour(iy + 1, ix + 2))))
			{
				if (ix + 1 < nx)
				{	*(mapP + 1) = 5;
				}
			}
			break;
		case EDGE_135:
			if ((!Fcontour(iy - 1, ix + 1) && Fcontour(iy - 2, ix + 2)) ||
				(!(Fcontour(iy, ix + 1) || Fcontour(iy - 1, ix)) &&
				(Fcontour(iy - 2, ix + 1) || Fcontour(iy - 1, ix + 2))))
			{
				if ((iy - 1 >= 0) && (ix + 1 < nx))
				{	*(mapP - nx + 1) = 5;
				}
			}
			break;
		case EDGE_X:
			if ((!Fcontour(iy - 1, ix) && Fcontour(iy - 2, ix)) ||
				(!(Fcontour(iy - 1, ix - 1) || Fcontour(iy - 1, ix + 1)) &&
				(Fcontour(iy - 2, ix - 1) || Fcontour(iy - 2, ix + 1))))
			{
				if (iy - 1 >= 0)
				{	*(mapP - nx) = 5;
				}
			}
			break;
		case EDGE_45:
			if ((!Fcontour(iy + 1, ix + 1) && Fcontour(iy + 2, ix + 2)) ||
				(!(Fcontour(iy, ix + 1) || Fcontour(iy + 1, ix)) &&
				(Fcontour(iy + 1, ix + 2) || Fcontour(iy + 2, ix + 1))))
			{
				if ((iy + 1 < ny) && (ix + 1 < nx))
				{	*(mapP + nx + 1) = 5;
				}
			}
			break;
		default:
			throw "EdgeDetector::image_neighbor_fwd() - error in case statements.";
	}
}

void EdgeDetector::image_neighbor_bwd(int iy, int ix, BYTE* mapP)
{
	switch (*mapP)
	{
		case 5:
			return;
		case EDGE_Y:
			if ((!Bcontour(iy, ix - 1) && Bcontour(iy, ix - 2)) ||
				(!(Bcontour(iy - 1, ix - 1) || Bcontour(iy + 1, ix - 1)) &&
				(Bcontour(iy - 1, ix - 2) || Bcontour(iy + 1, ix - 2))))
			{
				if (ix - 1 >= 0)
				{	*(mapP - 1) = 5;
				}
			}
			break;
		case EDGE_135:
			if ((!Bcontour(iy + 1, ix - 1) && Bcontour(iy + 2, ix - 2)) ||
				(!(Bcontour(iy, ix - 1) || Bcontour(iy + 1, ix)) &&
				(Bcontour(iy + 1, ix - 2) || Bcontour(iy + 2, ix - 1))))
			{
				if ((iy + 1 < ny) && (ix - 1 >= 0))
				{	*(mapP + nx - 1) = 5;
				}
			}
			break;
		case EDGE_X:
			if ((!Bcontour(iy + 1, ix) && Bcontour(iy + 2, ix)) ||
				(!(Bcontour(iy + 1, ix - 1) || Bcontour(iy + 1, ix + 1)) &&
				(Bcontour(iy + 2, ix - 1) || Bcontour(iy + 2, ix + 1))))
			{
				if (iy + 1 < ny)
				{	*(mapP + nx) = 5;
				}
			}
			break;
		case EDGE_45:
			if ((!Bcontour(iy - 1, ix - 1) && Bcontour(iy - 2, ix - 2)) ||
				(!(Bcontour(iy, ix - 1) || Bcontour(iy - 1, ix)) &&
				(Bcontour(iy - 1, ix - 2) || Bcontour(iy - 2, ix - 1))))
			{
				if ((iy - 1 >= 0) && (ix - 1 >= 0))
				{	*(mapP - nx - 1) = 5;
				}
			}
			break;
		default:
			throw "EdgeDetector::image_neighbor_bwd() - error in case statements.";
	}
}

int EdgeDetector::Fcontour(int iy, int ix)
{
	int edge_pt;
	int index;
	BYTE* mapP;
	index = iy * nx + ix;
	if (index < 0 || index >= nx*ny)
	{	return 0;
	}
	mapP = &edge_map[index];
	if (*mapP > 0)
	{	return 1;
	}
	else if ((edge_pt = image_edge_map_lo(ix, iy, index)) != 0)
	{	*mapP = edge_pt;
		image_neighbor_fwd(iy, ix, mapP);
		return 1;
	}
	else
	{	return 0;
	}
}

int EdgeDetector::Bcontour(int iy, int ix)
{
	BYTE* mapP;
	int edge_pt;
	int index;
	index = iy * nx + ix;
	if (index < 0 || index >= nx * ny)
	{	return 0;
	}
	mapP = &edge_map[index];
	if (*mapP > 0)
	{	return 1;
	}
	else if ((edge_pt = image_edge_map_lo (ix, iy, index)) != 0)
	{	*mapP = edge_pt;
		image_neighbor_bwd (iy, ix, mapP);
		return 1;
	}
	else
	{	return 0;
	}
}

int EdgeDetector::image_locx(int iy, int ix, int index)
{
	int* imP;
	if (ix < 2 || ix > nx - 2 || iy < 0 || iy > ny - 1)
	{	return 0;
	}
	imP = &idx[index];
	return *imP - *(imP - 1);
}

int EdgeDetector::image_locy(int iy, int ix, int index)
{
	int* imP;
	if (ix < 0 || ix > nx - 1 || iy < 2 || iy > ny - 2)
	{	return 0;
	}
	imP = &idy[index];
	return *imP - *(imP - nx);
}

int EdgeDetector::image_loc45(int iy, int ix, int index)
{
	int* imP;
	if (ix < 2 || ix > nx - 2 || iy < 2 || iy > ny - 2)
	{	return 0;
	}
	imP = &id45[index];
	return *imP - *(imP - nx + 1);
}

int EdgeDetector::image_loc135(int iy, int ix, int index)
{
	int* imP;
	if (ix < 2 || ix > nx - 2 || iy < 2 || iy > ny - 2)
	{	return 0;
	}
	imP = &id135[index];
	return *imP - *(imP - nx - 1);
}

int EdgeDetector::image_threshold()
{
	int ix;
	int iy;
	int* detx;
	int tmp;
	double mad = 0.0;
	detx = idx;
	for (iy = 0; iy < ny; iy++)
	{	for (ix = 0; ix < nx; ix++)
		{	tmp = *detx++;
			mad += abs(tmp);
		}
	}
	return (int)((double)(3 * mad / nx / ny) / 0.8);
}





/*******************************************************************************

PROJECT:			CS 297 Deliverable 3

FILE:				lap.cpp

PURPOSE:			Implementation to solve linear assignment problem
					in minimizing the total cost of matching each
					node in a bipartite graph.

COMPILER / TOOL:	Visual C++ .NET

PROGRAMMER:			Wallun Chan

START DATE:			4/1/05

AUTHOR:				1)	This code was downloaded from
						http://www.magiclogic.com/assignment.html
						and adapted to work for this application.

					2)	Version 1.0 - 4 September 1996
						author: Roy Jonker @ MagicLogic Optimization Inc.
						e-mail: roy_jonker@magiclogic.com

						Code for Linear Assignment Problem, according to

						"A Shortest Augmenting Path Algorithm for Dense and
						Sparse Linear Assignment Problems," Computing 38,
						325-340, 1987

						by

						R. Jonker and A. Volgenant, University of Amsterdam.

*******************************************************************************/

#include "stdafx.h"
#include "lap.h"

const int BIG = 1000000;

double lap(int dim, const MatInt& assigncost, VecInt& rowsol,
VecInt& colsol, VecInt& u, VecInt& v)
// input:
// dim			- problem size
// assigncost	- cost matrix
// rowsol		- column assigned to row in solution
// colsol		- row assigned to column in solution
// u			- dual variables, row reduction numbers
// v			- dual variables, column reduction numbers
{
	boolean unassignedfound;
	int f;
	int h;
	int i;
	int i0;
	int j;
	int j1;
	int j2;
	int k;
	int v2;
	int min;
	int umin;
	int usubmin;
	int imin;
	int endofpath;
	int last;
	int low;
	int up;
	int prvnumfree;
	int freerow;
	int numfree = 0;
	int *pred;
	int *free;
	int *collist;
	int *matches;
	int *d;

	free = new int[dim];       // list of unassigned rows.
	collist = new int[dim];    // list of columns to be scanned in various ways.
	matches = new int[dim];    // counts how many times a row could be assigned.
	d = new int[dim];         // 'cost-distance' in augmenting path calculation.
	pred = new int[dim];       // row-predecessor of column in augmenting/alternating path.

	// init how many times a row will be assigned in the column reduction.
	for (i = 0; i < dim; i++)
	{	matches[i] = 0;
	}

	// COLUMN REDUCTION
	for (j = dim - 1; j >= 0; j--)    // reverse order gives better results.
	{	// find minimum cost over rows.
		min = assigncost(0, j);
		imin = 0;
		for (i = 1; i < dim; i++)
		{	if (assigncost(i, j) < min)
			{	min = assigncost(i, j);
				imin = i;
			}
		}
		v[j] = min;
		if (++matches[imin] == 1)
		{	// init assignment if minimum row assigned for first time.
			rowsol[imin] = j;
			colsol[j] = imin;
		}
		else
		{	colsol[j] = -1;        // row already assigned, column not assigned.
		}
	}

	// REDUCTION TRANSFER
	for (i = 0; i < dim; i++)
	{	if (matches[i] == 0)     // fill list of unassigned 'free' rows.
		{	free[numfree++] = i;
		}
		else
		{	if (matches[i] == 1)   // transfer reduction from rows that are assigned once.
			{	j1 = rowsol[i];
				min = BIG;
				for (j = 0; j < dim; j++)
				{	if (j != j1)
					{	if (assigncost(i, j) - v[j] < min)
						{	min = assigncost(i, j) - v[j];
						}
					}
				}
				v[j1] = v[j1] - min;
			}
		}
	}

	// AUGMENTING ROW REDUCTION
	int loopcnt = 0;           // do-loop to be done twice.
	do
	{	loopcnt++;

		// scan all free rows.
		// in some cases, a free row may be replaced with another one to be scanned next.
		k = 0;
		prvnumfree = numfree;
		numfree = 0;             // start list of rows still free after augmenting row reduction.
		while (k < prvnumfree)
		{	i = free[k];
			k++;

			// find minimum and second minimum reduced cost over columns.
			umin = assigncost(i, 0) - v[0];
			j1 = 0;
			usubmin = BIG;
			for (j = 1; j < dim; j++)
			{	h = assigncost(i, j) - v[j];
				if (h < usubmin)
				{	if (h >= umin)
					{	usubmin = h;
						j2 = j;
					}
					else
					{	usubmin = umin;
						umin = h;
						j2 = j1;
						j1 = j;
					}
				}
			}

			i0 = colsol[j1];
			if (umin < usubmin)
			{	// change the reduction of the minimum column to increase the minimum
				// reduced cost in the row to the subminimum.
				v[j1] = v[j1] - (usubmin - umin);
			}
			else                   // minimum and subminimum equal.
			{	if (i0 >= 0)         // minimum column j1 is assigned.
				{	// swap columns j1 and j2, as j2 may be unassigned.
					j1 = j2;
					i0 = colsol[j2];
				}
			}

			// (re-)assign i to j1, possibly de-assigning an i0.
			rowsol[i] = j1;
			colsol[j1] = i;

			if (i0 >= 0)           // minimum column j1 assigned earlier.
			{	if (umin < usubmin)
				{	// put in current k, and go back to that k.
					// continue augmenting path i - j1 with i0.
					free[--k] = i0;
				}
				else
				{	// no further augmenting reduction possible.
					// store i0 in list of free rows for next phase.
					free[numfree++] = i0;
				}
			}
		}
	} while (loopcnt < 2);       // repeat once.

	// AUGMENT SOLUTION for each free row.
	for (f = 0; f < numfree; f++)
	{	freerow = free[f];       // start row of augmenting path.

		// Dijkstra shortest path algorithm.
		// runs until unassigned column added to shortest path tree.
		for (j = 0; j < dim; j++)
		{	d[j] = assigncost(freerow, j) - v[j];
			pred[j] = freerow;
			collist[j] = j;        // init column list.
		}

		low = 0; // columns in 0..low-1 are ready, now none.
		up = 0;  // columns in low..up-1 are to be scanned for current minimum, now none.
				 // columns in up..dim-1 are to be considered later to find new minimum,
				 // at this stage the list simply contains all columns
		unassignedfound = FALSE;
		do
		{	if (up == low)         // no more columns to be scanned for current minimum.
			{	last = low - 1;

				// scan columns for up..dim-1 to find all indices for which new minimum occurs.
				// store these indices between low..up-1 (increasing up).
				min = d[collist[up++]];
				for (k = up; k < dim; k++)
				{	j = collist[k];
					h = d[j];
					if (h <= min)
					{	if (h < min)     // new minimum.
						{	up = low;      // restart list at index low.
							min = h;
						}
						// new index with same minimum, put on undex up, and extend list.
						collist[k] = collist[up];
						collist[up++] = j;
					}
				}

				// check if any of the minimum columns happens to be unassigned.
				// if so, we have an augmenting path right away.
				for (k = low; k < up; k++)
				{	if (colsol[collist[k]] < 0)
					{	endofpath = collist[k];
						unassignedfound = TRUE;
						break;
					}
				}
			}

			if (!unassignedfound)
			{	// update 'distances' between freerow and all unscanned columns, via next scanned column.
				j1 = collist[low];
				low++;
				i = colsol[j1];
				h = assigncost(i, j1) - v[j1] - min;

				for (k = up; k < dim; k++)
				{	j = collist[k];
					v2 = assigncost(i, j) - v[j] - h;
					if (v2 < d[j])
					{	pred[j] = i;
						if (v2 == min)   // new column found at same minimum value
						{	if (colsol[j] < 0)
							{	// if unassigned, shortest augmenting path is complete.
								endofpath = j;
								unassignedfound = TRUE;
								break;
							}
							// else add to list to be scanned right away.
							else
							{	collist[k] = collist[up];
								collist[up++] = j;
							}
						}
						d[j] = v2;
					}
				}
			}
		} while (!unassignedfound);

		// update column prices.
		for (k = 0; k <= last; k++)
		{	j1 = collist[k];
			v[j1] = v[j1] + d[j1] - min;
		}

		// reset row and column assignments along the alternating path.
		do
		{	i = pred[endofpath];
			colsol[endofpath] = i;
			j1 = endofpath;
			endofpath = rowsol[i];
			rowsol[i] = j1;
		} while (i != freerow);
	}

	// calculate optimal cost.
	int lapcost = 0;
	for (i = 0; i < dim; i++)
	{	j = rowsol[i];
		u[i] = assigncost(i, j) - v[j];
		lapcost = lapcost + assigncost(i, j);
	}

	// free reserved memory.
	delete [] pred;
	delete [] free;
	delete [] collist;
	delete [] matches;
	delete [] d;

	return lapcost;
}





/***************************************************************************

PROJECT:			CS 297 Deliverable 3

FILE:				tps.h

PURPOSE:			Class definition for code that implements the
					thin-plate-spline method of landmark alignment.

COMPILER / TOOL:	Visual C++ .NET

PROGRAMMER:			Wallun Chan

START DATE:			4/1/05

***************************************************************************/

#ifndef _THIN_PLATE_SPLINE_H_
#define _THIN_PLATE_SPLINE_H_

#include "utility.h"

class TPSSolver
{
public:
	TPSSolver();
	~TPSSolver();
	void setPoints(const VecIntPr& trgtPts, const VecIntPr& testPts);
	void mapPoints(VecIntPr& pts);
private:
	double shapeFunction(double radius);
	Matrix lMatrix;
	Matrix wMatrix;
	VecIntPr samples;
};

#endif





/***************************************************************************

PROJECT:			CS 297 Deliverable 3

FILE:				tps.cpp

PURPOSE:			Implementation of the thin-plate-spline method
					for point-set alignment.

COMPILER / TOOL:	Visual C++ .NET

PROGRAMMER:			Wallun Chan

START DATE:			4/1/05

REMARKS:			This code is based on the following:

					F. L. Bookstein. Principal warps:
					thin-plate splines and decomposition of deformations.
					IEEE Trans. Pattern Analysis and Machine Intelligence, June 1989,

					and

					Gianluca Donato and Serge Belongie. Approximation
					Methods for Thin Plate Spline Mappings and
					Principal Warps. October 2002;
					http://www-cse.ucsd.edu/~sjb/pami_tps.pdf.

					and

					Jarno Elonen. Thin Plate Spline editor -
					An example program in C++.
					http://elonen.iki.fi/code/tpsdemo, 2003.

***************************************************************************/

#include "stdafx.h"
#include "tps.h"

TPSSolver::TPSSolver()
// Default constructor
{
}

TPSSolver::~TPSSolver()
// Destructor
{
}

void TPSSolver::setPoints(const VecIntPr& trgtPts, const VecIntPr& testPts)
/* PURPOSE:		This function takes a set of test points and maps it to
				the set of target points.  The function implements the
				thin-plate-spline algorithm that uses the notion of a
				thin plate "tacked onto" a set of control points and uses
				the bending of the surface defined by splines on the control
				points as a way of describing the mapping of one set of points
				to another.  The function solves a matrix formula for a set of
				weights and factors to provide a closed form equation mapping
				the positions of one set of points to another.

   RECEIVES:	Set of test and target points such that the test points are
				mapped onto the target points.
*/
{
	int p = (int)testPts.size();
	lMatrix.SetSize(p + 3, p + 3);
	int i;
	int j;

	// Set K matrix
	for (i = 0; i < p; i++)
	{	IntPr pt1 = testPts[i];
		for (j = 0; j < p; j++)
		{	if (i != j)
			{	IntPr pt2 = testPts[j];
				double diffX = pt1.first - pt2.first;
				double diffY = pt1.second - pt2.second;
				lMatrix(i, j) = shapeFunction(diffX * diffX + diffY * diffY);
			}
			else
			{	lMatrix(i, j) = 0.0;
			}
		}
	}

	// Set P matrix and its transpose
	for (i = 0; i < p; i++)
	{	IntPr pt = testPts[i];
		lMatrix(i, p) = 1.0;
		lMatrix(i, p + 1) = pt.first;
		lMatrix(i, p + 2) = pt.second;
		lMatrix(p, i) = 1.0;
		lMatrix(p + 1, i) = pt.first;
		lMatrix(p + 2, i) = pt.second;
	}

	// Set O matrix
	for (i = p; i < p + 3; i++)
	{	lMatrix(i, p) = 0.0;
		lMatrix(i, p + 1) = 0.0;
		lMatrix(i, p + 2) = 0.0;
	}

	// Set v column vectors
	Matrix yMatrix(p + 3, 2);
    for (i = 0; i < p; i++)
    {	IntPr targPt = trgtPts[i];
		IntPr testPt = testPts[i];
		yMatrix(i, 0) = targPt.first - testPt.first;
		yMatrix(i, 1) = targPt.second - testPt.second;
    }

	// Set o column vectors
    yMatrix(p, 0) = yMatrix(p + 1, 0) = yMatrix(p + 2, 0) = 0.0;
    yMatrix(p, 1) = yMatrix(p + 1, 1) = yMatrix(p + 2, 1) = 0.0;

	// Solve for w and a coefficients
	wMatrix = lMatrix.Solve(yMatrix);
	samples = testPts;
}

void TPSSolver::mapPoints(VecIntPr& pts)
/* PURPOSE:		This function maps a set of points to the test
				image space by using the mapping function
				determined by the function setPoints().

   RECEIVES:	Set of points to be mapped test image space.
*/
{
	int n = (int)pts.size();
	int p = (int)samples.size();
	int i;
	int j;
	for (i = 0; i < n; i++)
	{	double x = pts[i].first;
		double y = pts[i].second;

		// dx = a1x + a2x * x + a3x * y
		// dy = a1y + a2y * x + a3y * y
		double dx = wMatrix(p, 0) + wMatrix(p + 1, 0) * x + wMatrix(p + 2, 0) * y;
		double dy = wMatrix(p, 1) + wMatrix(p + 1, 1) * x + wMatrix(p + 2, 1) * y;
		for (j = 0; j < p; j++)
		{	double diffX = samples[j].first - x;
			double diffY = samples[j].second - y;
			double u = shapeFunction(diffX * diffX + diffY * diffY);
			dx += wMatrix(j, 0) * u;
			dy += wMatrix(j, 1) * u;
		}
		pts[i].first += (int)dx;
		pts[i].second += (int)dy;
	}
}

double TPSSolver::shapeFunction(double rad)
/* PURPOSE:		To calculate the shap

   RECEIVES:	Radial distance between two points.

   REMARKS:		Input is passed in as a square since don't want to
				perform a square root, given that the equation is
				shapeFunction(rad) = rad ^ 2 * log(rad).
*/
{
	// log(r ^ 2) = 2 * log(r)
	return (rad == 0.0) ? 0.0 : rad * log(rad) / 2.0;

}





/*************************************************************************************

PROJECT:			CS 297 Deliverable 3

FILE:				convexhull.cpp

PURPOSE:			Implements the Chain Hull algorithm to calculate the convex
					hull of a set of points.

COMPILER / TOOL:	Visual C++ .NET

PROGRAMMER:			Wallun Chan

START DATE:			4/1/05

REMARKS:			This code was borrowed from
					http://softsurfer.com/Archive/algorithm_0109/algorithm_0109.htm
					and adapted to work with this application.

*************************************************************************************/

#include "stdafx.h"
#include "convexhull.h"

int isLeft(IntPr P0, IntPr P1, IntPr P2)
/* PURPOSE:		Tests if a point is Left|On|Right of an infinite line.

   RECEIVES:	Three points P0, P1, and P2.

   RETURNS:		> 0 for P2 left of the line through P0 and P1.
				= 0 for P2 on the line.
				< 0 for P2 right of the line.
*/
{
    return (P1.first - P0.first) * (P2.second - P0.second) -
		(P2.first - P0.first) * (P1.second - P0.second);
}

bool comparator(IntPr element1, IntPr element2)
/* PURPOSE:		Used to compare the coordinates of two elements.

   RECEIVES:	Two elements to be compared.

   RETURNS:		True if x-coordinate of first element is less
				than the second element, or if the y-coordinates
				are the same, then return true if x-coordinate
				of first element is less than the second element.
				Otherwise, return false.
*/
{
	return (element1.first < element2.first ||
		    (element1.first == element2.first &&
			element1.second < element2.second));
}

VecIntPr convexHull(VecIntPr pts)
{
	// hull will be used as the stack
	VecIntPr hull;
	if (pts.empty())
	{	return hull;
	}
	hull.resize(pts.size());
	sort(pts.begin(), pts.end(), comparator);
	int bot = 0;	// index for bottom and top of the stack
	int top = -1;	// index for top of the stack

	// Get the indices of points with min x-coord and min|max y-coord
	int minmin = 0;
	int minmax;
	int xmin = pts[0].first;
	int n = (int)pts.size();
	int i;
	for (i = 1; i < n; i++)
	{	if (pts[i].first != xmin)
		{	break;
		}
	}
	minmax = i - 1;
	if (minmax == n - 1)	// degenerate case: all x-coords == xmin
	{	hull[++top] = pts[minmin];
		if (pts[minmax].second != pts[minmin].second)	// a nontrivial segment
		{    hull[++top] = pts[minmax];
		}
		hull[++top] = pts[minmin];
		hull.resize(top + 1);	// add polygon endpoint
		return hull;
	}

	// Get the indices of points with max x-coord and min|max y-coord
    int maxmin;
	int maxmax = n - 1;
    int xmax = pts[n - 1].first;
    for (i = n - 2; i >= 0; i--)
	{   if (pts[i].first != xmax)
		{	break;
		}
	}
	maxmin = i + 1;

	// Compute the lower hull on the stack H
	hull[++top] = pts[minmin];	// push minmin point onto stack
	i = minmax;
	while (++i <= maxmin)
	{	// the lower line joins P[minmin] with P[maxmin]
		if (isLeft(pts[minmin], pts[maxmin], pts[i]) >= 0 && i < maxmin)
		{	continue;	// ignore P[i] above or on the lower line
		}
		while (top > 0)		// there are at least 2 points on the stack
		{	// test if P[i] is left of the line at the stack top
			if (isLeft(hull[top - 1], hull[top], pts[i]) > 0)
			{	break;	// P[i] is a new hull vertex
			}
			else
			{	top--;	// pop top point off stack
			}
		}
		hull[++top] = pts[i];	// push P[i] onto stack
	}

	// Next, compute the upper hull on the stack H above the bottom hull
    if (maxmax != maxmin)	// if distinct xmax points
	{    hull[++top] = pts[maxmax];		// push maxmax point onto stack
	}
    bot = top;		// the bottom point of the upper hull stack
    i = maxmin;
	while (--i >= minmax)
	{	// the upper line joins P[maxmax] with P[minmax]
		if (isLeft(pts[maxmax], pts[minmax], pts[i]) >= 0 && i > minmax)
		{    continue;		// ignore P[i] below or on the upper line
		}
		while (top > bot)	// at least 2 points on the upper stack
		{	// test if P[i] is left of the line at the stack top
			if (isLeft(hull[top-1], hull[top], pts[i]) > 0)
			{	break;		// P[i] is a new hull vertex
			}
			else
			{	top--;		// pop top point off stack
			}
		}
		hull[++top] = pts[i];	// push P[i] onto stack
	}
    if (minmax != minmin)
	{    hull[++top] = pts[minmin];	// push joining endpoint onto stack
	}
	hull.resize(top + 1);
    return hull;
}





/***************************************************************************

PROJECT:			CS 297 Deliverable 3

FILE:				directory.cpp

PURPOSE:			Implementation that contains code to support the
					display, querying, and traversal of directory and
					file path structures in MFC.

COMPILER / TOOL:	Visual C++ .NET

PROGRAMMER:			Wallun Chan

START DATE:			4/1/05

COMMENTS:			Code was borrowed and adapted from "Windows
					Programming: Programmer's Notebook", by Mario
					Giannini and Jim Keogh, 2001; code is found in
					PNDirTree.cpp in that book.

***************************************************************************/

#include "stdafx.h"
#include "directory.h"

/*** Functions to display, query, and traverse directories and files in MFC. ***/

int CStringCmp(const void* arg1, const void* arg2)
/* PURPOSE:		Comparison for two CStrings, for qsort call below

   RECEIVES:	Pointers to items to be compared.

   RETURNS:		-1 if arg1 < arg2, 1 if arg1 > arg2, and 0 if arg1 == arg2.

   REMARKS:		1) Assumes that function parameters have the proper
				   have the needed comparison operators.

				2) This code was borrowed from "Windows Programming:
				   Programmer's Notebook, by Mario Giannini and Jim Keogh,
				   2001; code found in PNDirTree.cpp in that book.
*/
{
	return (((CString*)arg1)->CompareNoCase(*(CString*)arg2));
}

void SlashPath(CString& Dest, bool MustHave)
/* PURPOSE:		Ensures that a path either has, or does not have,
				a slash at the end.

   RECEIVES:	CString object to slash or unslash depending on
				whether "MustHave" is true or false, respectively.

   REMARKS:		This code was borrowed from "Windows Programming:
				Programmer's Notebook, by Mario Giannini and Jim
				Keogh, 2001; code found in PNDirTree.cpp in that
				book.
*/
{
	char Now = '\0';
	if (!Dest.IsEmpty())
	{	Now = Dest[Dest.GetLength() - 1];
	}
	if (MustHave && Now != '\\')
	{	Dest += '\\';
	}
	else if (!MustHave && Now == '\\')
	{	Dest.TrimRight('\\');
	}
}

LPCSTR GetSubPath(const char* Str)
/* PURPOSE:		To strip a path down to its file name.

   RECEIVES:	String name of full file path to process.

   RETURNS:		File name with path stripped.

   REMARKS:		This code was borrowed from "Windows Programming:
				Programmer's Notebook, by Mario Giannini and Jim
				Keogh, 2001; code found in PNDirTree.cpp in that
				book.
*/
{
	static CString Path;
	Path = Str;
	SlashPath(Path, false);
	int i = Path.ReverseFind('\\');
	if (i != -1)
	{	Path = Path.Mid(i + 1);
	}
	return Path;
}

int BuildFileList(CStringArray& Dest, const char* Path)
/* PURPOSE:		Builds a list of file names for a path.

   RECEIVES:	CStringArray object to store file names that are in "Path".

   RETURNS:		Number of files in directory specified by "Path".

   REMARKS:		This code was borrowed and adapted from "Windows
				Programming: Programmer's Notebook, by Mario
				Giannini and Jim Keogh, 2001; code found in
				PNDirTree.cpp in that book.
*/
{
	int Ret = BuildListEx(Dest, Path, true);
	int i;
	for (i = 0; i < Ret; i++)
	{	Dest[i].Delete(0, 1);
	}
	return (Ret);
}

int BuildDirList(CStringArray& Dest, const char* Path)
/* PURPOSE:		Builds a list of file names for a path.

   RECEIVES:	CStringArray object to store directory names that are in "Path".

   RETURNS:		Number of directories in directory specified by "Path".

   REMARKS:		This code was borrowed and adapted from "Windows
				Programming: Programmer's Notebook, by Mario
				Giannini and Jim Keogh, 2001; code found in
				PNDirTree.cpp in that book.
*/
{
	int Ret = BuildListEx(Dest, Path, false);
	int i;
	for (i = 0; i < Ret; i++)
	{	Dest[i].Delete(0, 1);
	}
	return (Ret);
}

int BuildListEx(CStringArray& Dest, const char* Path, bool showFiles)
/* PURPOSE:		Builds a list of file and directory names for a path.

   RECEIVES:	CStringArray object to store file and directory names
				that are in "Path", and boolean "showFiles" to determine
				whether file names should also be stored.

   RETURNS:		Total number of files and directories in directory
				specified by "Path".

   REMARKS:		This code was borrowed and adapted from "Windows
				Programming: Programmer's Notebook, by Mario
				Giannini and Jim Keogh, 2001; code found in
				PNDirTree.cpp in that book.
*/
{
	CFileFind Find;
	CString MyPath = Path;
	BOOL Found;
	CWaitCursor Wait;
	SlashPath(MyPath, true);
	MyPath += "*.*";
	Found = Find.FindFile(MyPath);
	while (Found)
	{	Found = Find.FindNextFile();
		if (Find.IsDirectory() && !Find.IsDots())
		{	Dest.Add("A" + Find.GetFilePath());
		}
		if (showFiles && !Find.IsDirectory())
		{	Dest.Add("B" + Find.GetFilePath());
		}
	}
	qsort((void*)Dest.GetData(), Dest.GetSize(), sizeof(CString), CStringCmp);
	return (int)Dest.GetSize();
}

BOOL DemoBrowseForFolder(char* Dest, char* caption, bool showFiles)
/* PURPOSE:		Allows user to select a folder through a dialog GUI where
				directories and file names are displayed in a complete
				directory structure.

   RECEIVES:	Windows app handle hWnd, char string "Dest" to store
				user selected file or directory name, and boolean
				"showFiles" to determine whether files should be
				displayed in dialog GUI.

   RETURNS:		BOOL value determining whether user selected a file,
				or cancelled.

   REMARKS:		This code was borrowed and adapted from "Windows
				Programming: Programmer's Notebook, by Mario
				Giannini and Jim Keogh, 2001; code found in
				PNDirTree.cpp in that book.
*/
{
	BROWSEINFO bi;
	ITEMIDLIST* pItemIDList;
	IMalloc *pMalloc;
	if (CoGetMalloc(1, &pMalloc) != S_OK)
	{	return FALSE;
	}
	char FolderOnly[_MAX_PATH];
	memset(&bi, 0, sizeof(bi));
	bi.hwndOwner = 0;
	bi.lpszTitle = caption;
	bi.pszDisplayName = FolderOnly;
	if (showFiles)
	{	bi.ulFlags = BIF_BROWSEINCLUDEFILES;
	}
	if ((pItemIDList = SHBrowseForFolder(&bi)) != NULL)
	{	SHGetPathFromIDList(pItemIDList, Dest);
		pMalloc->Free(pItemIDList);
		return TRUE;
	}
	return FALSE;
}