static float measure_facepair(EditVert *v1, EditVert *v2, EditVert *v3, EditVert *v4, float limit){
	
	/*gives a 'weight' to a pair of triangles that join an edge to decide how good a join they would make*/
	float	measure = 0.0, noA1[3], noA2[3], noB1[3], noB2[3], normalADiff, normalBDiff,
			edgeVec1[3], edgeVec2[3], edgeVec3[3], edgeVec4[3], diff,
			minarea, maxarea, areaA, areaB;
	
	/*First Test: Normal difference*/
	CalcNormFloat(v1->co, v2->co, v3->co, noA1);
	CalcNormFloat(v1->co, v3->co, v4->co, noA2);
	
	if(noA1[0] == noA2[0] && noA1[1] == noA2[1] && noA1[2] == noA2[2]) normalADiff = 0.0;
	else normalADiff = VecAngle2(noA1, noA2);
		//if(!normalADiff) normalADiff = 179;
	CalcNormFloat(v2->co, v3->co, v4->co, noB1);
	CalcNormFloat(v4->co, v1->co, v2->co, noB2);
	
	if(noB1[0] == noB2[0] && noB1[1] == noB2[1] && noB1[2] == noB2[2]) normalBDiff = 0.0;
	else normalBDiff = VecAngle2(noB1, noB2);
		//if(!normalBDiff) normalBDiff = 179;
	
	measure += (normalADiff/360) + (normalBDiff/360);
	if(measure > limit) return measure;
	
	/*Second test: Colinearity*/
	VecSubf(edgeVec1, v1->co, v2->co);
	VecSubf(edgeVec2, v2->co, v3->co);
	VecSubf(edgeVec3, v3->co, v4->co);
	VecSubf(edgeVec4, v4->co, v1->co);
	
	diff = 0.0;
	
	diff = (
		fabs(VecAngle2(edgeVec1, edgeVec2) - 90) +
		fabs(VecAngle2(edgeVec2, edgeVec3) - 90) + 
		fabs(VecAngle2(edgeVec3, edgeVec4) - 90) + 
		fabs(VecAngle2(edgeVec4, edgeVec1) - 90)) / 360;
	//if(!diff) return 0.0;
	
	measure +=  diff;
	if(measure > limit) return measure;

	/*Third test: Concavity*/
	areaA = AreaT3Dfl(v1->co, v2->co, v3->co) + AreaT3Dfl(v1->co, v3->co, v4->co);
	areaB = AreaT3Dfl(v2->co, v3->co, v4->co) + AreaT3Dfl(v4->co, v1->co, v2->co);
	
	if(areaA <= areaB) minarea = areaA;
	else minarea = areaB;
	
	if(areaA >= areaB) maxarea = areaA;
	else maxarea = areaB;
	
	if(!maxarea) measure += 1;
	else measure += (1 - (minarea / maxarea));

	return measure;
}

#define T2QUV_LIMIT 0.005
#define T2QCOL_LIMIT 3
/*Bitflags for delimters on face compare*/
#define T2QCMPUV  1
#define T2QCMPVCOL 2
static int compareFaceAttribs(EditFace *f1, EditFace *f2, int flag) 
{
	/*Test to see if two faces with custom data match within a certain limit*/
	EditVert *va, *vb;
	int i, j;
	MTFace *tf1, *tf2;
	unsigned int *col1, *col2;

	tf1 = CustomData_em_get(&G.editMesh->fdata, f1->data, CD_MTFACE);
	tf2 = CustomData_em_get(&G.editMesh->fdata, f2->data, CD_MTFACE);

	col1 = CustomData_em_get(&G.editMesh->fdata, f1->data, CD_MCOL);
	col2 = CustomData_em_get(&G.editMesh->fdata, f2->data, CD_MCOL);
	
	if(flag & T2QCMPUV){
		if(tf1 == NULL || tf2 == NULL) return 1;
		else if(tf1->tpage == NULL && tf2->tpage == NULL)
			return 1;
		else if(tf1->tpage != tf2->tpage)
			return 0;
	}
	if(flag & T2QCMPVCOL)
		if(!col1 || !col2) return 1;

	i = 0;
	va = f1->v1;
	while(i < 4 && va != NULL){
		j = 0;
		vb = f2->v1;
		while(j < 4 && vb != NULL){
			if(va == vb){
				if(flag & T2QCMPUV){
					if( !(tf1->uv[i][0] + T2QUV_LIMIT > tf2->uv[j][0] && tf1->uv[i][0] - T2QUV_LIMIT < tf2->uv[j][0] &&
						tf1->uv[i][1] + T2QUV_LIMIT > tf2->uv[j][1] && tf1->uv[i][1] - T2QUV_LIMIT < tf2->uv[j][1]))
						return 0;
				}
				if(flag & T2QCMPVCOL){
					char *v1, *v2;
					v1 = (char *)&(col1[i]);
					v2 = (char *)&(col2[j]);
		
					if(!(v1[1] + T2QCOL_LIMIT > v2[1] && v1[1] - T2QCOL_LIMIT < v2[1] &&
						v1[2] + T2QCOL_LIMIT > v2[2] && v1[2] - T2QCOL_LIMIT < v2[2] &&
						v1[3] + T2QCOL_LIMIT > v2[3] && v1[3] - T2QCOL_LIMIT < v2[3]))
						return 0;
				}
			}
			j++;
			vb++;
		}
		i++;
		va++;
	}
	return 1;
}	
	
static int fplcmp(const void *v1, const void *v2)
{
	const EditEdge *e1= *((EditEdge**)v1), *e2=*((EditEdge**)v2);
	
	if( e1->crease > e2->crease) return 1;
	else if( e1->crease < e2->crease) return -1;
	
	return 0;
}

/*Bitflags for edges.*/
#define T2QDELETE	1
#define T2QCOMPLEX	2
#define T2QSIMPLE	4
void join_triangles(void)
{
	EditMesh *em=G.editMesh;
	EditVert *v1, *v2, *v3, *v4, *eve;
	EditEdge *eed, **edsortblock = NULL, **edb = NULL;
	EditFace *efa;
	EVPTuple *efaar = NULL;
	EVPtr *efaa = NULL;
	float *creases = NULL;
	float limit = G.scene->toolsettings->select_thresh * 10; /*Used to set tolerance*/
	int i, ok, totedge=0, totseledge=0, complexedges, vindex[4];
	
	/*test for multi-resolution data*/
	if(multires_test()) return;
	
	/*if we take a long time on very dense meshes we want waitcursor to display*/
	waitcursor(1);
	
	totseledge = count_selected_edges(em->edges.first);
	if(totseledge==0) return;
	
	/*abusing crease value to store weights for edge pairs. Nasty*/
	for(eed=em->edges.first; eed; eed=eed->next) totedge++;
	if(totedge) creases = MEM_callocN(sizeof(float) * totedge, "Join Triangles Crease Array"); 
	for(eed=em->edges.first, i = 0; eed; eed=eed->next, i++){
		creases[i] = eed->crease; 
		eed->crease = 0.0;
	}
	
	/*clear temp flags*/
	for(eve=em->verts.first; eve; eve=eve->next) eve->f1 = eve->f2 = 0;
	for(eed=em->edges.first; eed; eed=eed->next) eed->f2 = eed->f1 = 0;
	for(efa=em->edges.first; efa; efa=efa->next) efa->f1 = 0;

	/*For every 2 manifold edge, create pointers to its two faces.*/
	efaar= (EVPTuple *) MEM_callocN(totseledge * sizeof(EVPTuple), "Tri2Quad");
	ok = collect_quadedges(efaar, em->edges.first, em->faces.first);
	
	if(ok){
		/*clear tmp.l flag and store number of faces that are selected and coincident to current face here.*/  
		for(efa=em->faces.first; efa; efa=efa->next) efa->tmp.l = 0;
		for(eed=em->edges.first; eed; eed=eed->next){
			if(eed->tmp.p){
				efaa= (EVPtr *) eed->tmp.p;
				if(efaa[0]) efaa[0]->tmp.l++;
				if(efaa[1]) efaa[1]->tmp.l++;
			}
		}
		complexedges = 0;
		for(eed=em->edges.first; eed; eed=eed->next){
			if(eed->tmp.p){
				efaa= (EVPtr *) eed->tmp.p;
				if(efaa[0] && efaa[1]){
					v1 = v2 = v3 = v4 = NULL;
					givequadverts(efaa[0], efaa[1], &v1, &v2, &v3, &v4, vindex);
					if(v1 && v2 && v3 && v4){
						/*test if simple island first. This mimics 2.42 behaviour and the tests are more restrictive.*/
						if(efaa[0]->tmp.l == 1 && efaa[1]->tmp.l == 1){
							eed->f1 |= T2QSIMPLE;
							if( convex(v1->co, v2->co, v3->co, v4->co) ) {
								/*mark joined faces/edges for deletion*/
								efaa[0]->f1 = 1;
								efaa[1]->f1 = 1;
								eed->f1 |= T2QDELETE;
							}
						}
						else{ /*its part of a 'complex' island, tests are different*/
							if(efaa[0]->mat_nr != efaa[1]->mat_nr); /*do nothing*/
							else if(eed->sharp); /*do nothing*/
							else{
								ok = compareFaceAttribs(efaa[0], efaa[1], T2QCMPUV|T2QCMPVCOL);
								if(ok){
									measure = measure_facepair(v1, v2, v3, v4, limit);
									if(measure < limit){
										complexedges++;
										eed->f1 |= T2QCOMPLEX;
										eed->crease = measure; /*we dont mark edges for join YET*/
									}
								}
							}
						}
					}
				}
			}
		}
		
		/*Quicksort the complex edges and mark the most appropriate ones for deletion*/
		if(complexedges){
			edsortblock = edb = MEM_callocN(sizeof(EditEdge*) * complexedges, "Face Pairs quicksort Array");
			for(eed = em->edges.first; eed; eed=eed->next){
				if(eed->f1 & T2QCOMPLEX){
					*edb = eed;
					edb++;
				}
			}
			qsort(edsortblock, complexedges, sizeof(EditEdge*), fplcmp);
			/*now go through and mark the edges who get the highest weighting*/
			for(edb=edsortblock, i=0; i < complexedges; edb++, i++){ 
				efaa = (EVPtr *)((*edb)->tmp.p);
				if( !(efaa[0]->f1) && !(efaa[1]->f1) ){
					v1 = v2 = v3 = v4 = NULL;
					givequadverts(efaa[0], efaa[1], &v1, &v2, &v3, &v4, vindex);
					if(v1 && v2 && v3 && v4){
						/*mark elements for deletion later*/
						efaa[0]->f1 = 1;
						efaa[1]->f1 = 1;
						(*edb)->f1 |= T2QDELETE;
					}
				}
			}
		}
		/*finally go through all edges marked for deletion and join their faces*/ 
		for(eed=em->edges.first; eed; eed=eed->next){
			if(eed->f1 & T2QDELETE){
				efaa= (EVPtr *)eed->tmp.p;
				v1 = v2 = v3 = v4 = NULL;
				givequadverts(efaa[0], efaa[1], &v1, &v2, &v3, &v4, vindex);
				if(exist_face(v1, v2, v3, v4)==0) {
					efa = EM_face_from_faces(efaa[0], efaa[1], vindex[0], vindex[1], 4+vindex[2], 4+vindex[3]);
					EM_select_face(efa,1);
				}
			}
		}
	}
	
	/*free data and cleanup*/
	if(creases){
		for(eed=em->edges.first, i = 0; eed; eed=eed->next, i++) eed->crease = creases[i]; 
		MEM_freeN(creases);
	}
	for(eed=em->edges.first; eed; eed=eed->next){
		if(eed->f1 & T2QDELETE) eed->f1 = 1;
		else eed->f1 = 0;
	}
	free_tagged_facelist(em->faces.first);
	free_tagged_edgelist(em->edges.first);
	if(efaar) MEM_freeN(efaar);
	if(edsortblock) MEM_freeN(edsortblock);
		
	EM_selectmode_flush();
	countall();
	allqueue(REDRAWVIEW3D, 0);
	DAG_object_flush_update(G.scene, G.obedit, OB_RECALC_DATA);
	#ifdef WITH_VERSE
	if(G.editMesh->vnode)
		sync_all_versefaces_with_editfaces((VNode*)G.editMesh->vnode);
	#endif
	waitcursor(0);
	BIF_undo_push("Convert Triangles to Quads");
}