static void get_strand_normal(Material *ma, const float surfnor[3], float surfdist, float nor[3])
{
	float cross[3], nstrand[3], vnor[3], blend;

	if (!((ma->mode & MA_STR_SURFDIFF) || (ma->strand_surfnor > 0.0f)))
		return;

	if (ma->mode & MA_STR_SURFDIFF) {
		cross_v3_v3v3(cross, surfnor, nor);
		cross_v3_v3v3(nstrand, nor, cross);

		blend = dot_v3v3(nstrand, surfnor);
		CLAMP(blend, 0.0f, 1.0f);

		interp_v3_v3v3(vnor, nstrand, surfnor, blend);
		normalize_v3(vnor);
	}
	else {
		copy_v3_v3(vnor, nor);
	}

	if (ma->strand_surfnor > 0.0f) {
		if (ma->strand_surfnor > surfdist) {
			blend = (ma->strand_surfnor - surfdist) / ma->strand_surfnor;
			interp_v3_v3v3(vnor, vnor, surfnor, blend);
			normalize_v3(vnor);
		}
	}

	copy_v3_v3(nor, vnor);
}

/* our own triangle intersection, so we can fully control the epsilons and
 * prevent corner case from going wrong*/
static int meshdeform_tri_intersect(const float orig[3], const float end[3], const float vert0[3],
                                    const float vert1[3], const float vert2[3],
                                    float r_isectco[3], float r_uvw[3])
{
	float edge1[3], edge2[3], tvec[3], pvec[3], qvec[3];
	float det, inv_det, u, v, dir[3], isectdir[3];

	sub_v3_v3v3(dir, end, orig);

	/* find vectors for two edges sharing vert0 */
	sub_v3_v3v3(edge1, vert1, vert0);
	sub_v3_v3v3(edge2, vert2, vert0);

	/* begin calculating determinant - also used to calculate U parameter */
	cross_v3_v3v3(pvec, dir, edge2);

	/* if determinant is near zero, ray lies in plane of triangle */
	det = dot_v3v3(edge1, pvec);

	if (UNLIKELY(det == 0.0f)) {
		return 0;
	}

	inv_det = 1.0f / det;

	/* calculate distance from vert0 to ray origin */
	sub_v3_v3v3(tvec, orig, vert0);

	/* calculate U parameter and test bounds */
	u = dot_v3v3(tvec, pvec) * inv_det;
	if (u < -EPSILON || u > 1.0f + EPSILON)
		return 0;

	/* prepare to test V parameter */
	cross_v3_v3v3(qvec, tvec, edge1);

	/* calculate V parameter and test bounds */
	v = dot_v3v3(dir, qvec) * inv_det;
	if (v < -EPSILON || u + v > 1.0f + EPSILON)
		return 0;

	r_isectco[0] = (1.0f - u - v) * vert0[0] + u * vert1[0] + v * vert2[0];
	r_isectco[1] = (1.0f - u - v) * vert0[1] + u * vert1[1] + v * vert2[1];
	r_isectco[2] = (1.0f - u - v) * vert0[2] + u * vert1[2] + v * vert2[2];

	r_uvw[0] = 1.0f - u - v;
	r_uvw[1] = u;
	r_uvw[2] = v;

	/* check if it is within the length of the line segment */
	sub_v3_v3v3(isectdir, r_isectco, orig);

	if (dot_v3v3(dir, isectdir) < -EPSILON)
		return 0;
	
	if (dot_v3v3(dir, dir) + EPSILON < dot_v3v3(isectdir, isectdir))
		return 0;

	return 1;
}

/**
 * Add randomness to stroke
 * \param gps: Stroke data
 * \param brush: Brush data
 */
void gp_randomize_stroke(bGPDstroke *gps, bGPDbrush *brush)
{
	bGPDspoint *pt1, *pt2, *pt3;
	float v1[3];
	float v2[3];
	if (gps->totpoints < 3) {
		return;
	}

	/* get two vectors using 3 points */
	pt1 = &gps->points[0];
	pt2 = &gps->points[1];
	pt3 = &gps->points[(int)(gps->totpoints * 0.75)];

	sub_v3_v3v3(v1, &pt2->x, &pt1->x);
	sub_v3_v3v3(v2, &pt3->x, &pt2->x);
	normalize_v3(v1);
	normalize_v3(v2);

	/* get normal vector to plane created by two vectors */
	float normal[3];
	cross_v3_v3v3(normal, v1, v2);
	normalize_v3(normal);

	/* get orthogonal vector to plane to rotate random effect */
	float ortho[3];
	cross_v3_v3v3(ortho, v1, normal);
	normalize_v3(ortho);

	/* Read all points and apply shift vector (first and last point not modified) */
	for (int i = 1; i < gps->totpoints - 1; ++i) {
		bGPDspoint *pt = &gps->points[i];
		/* get vector with shift (apply a division because random is too sensitive */
		const float fac = BLI_frand() * (brush->draw_random_sub / 10.0f);
		float svec[3];
		copy_v3_v3(svec, ortho);
		if (BLI_frand() > 0.5f) {
			mul_v3_fl(svec, -fac);
		}
		else {
			mul_v3_fl(svec, fac);
		}

		/* apply shift */
		add_v3_v3(&pt->x, svec);
	}

}

float angle_signed_on_axis_v3v3v3_v3(const float v1[3], const float v2[3], const float v3[3], const float axis[3])
{
	float v1_proj[3], v2_proj[3], tproj[3];
	float angle;

	sub_v3_v3v3(v1_proj, v1, v2);
	sub_v3_v3v3(v2_proj, v3, v2);

	/* project the vectors onto the axis */
	project_v3_v3v3(tproj, v1_proj, axis);
	sub_v3_v3(v1_proj, tproj);

	project_v3_v3v3(tproj, v2_proj, axis);
	sub_v3_v3(v2_proj, tproj);

	angle = angle_v3v3(v1_proj, v2_proj);

	/* calculate the sign (reuse 'tproj') */
	cross_v3_v3v3(tproj, v2_proj, v1_proj);
	if (dot_v3v3(tproj, axis) < 0.0f) {
		angle = ((float)(M_PI * 2.0)) - angle;
	}

	return angle;
}

/* adjust bone roll to align Z axis with vector
 * vec is in local space and is normalized
 */
float ED_rollBoneToVector(EditBone *bone, const float align_axis[3], const short axis_only)
{
	float mat[3][3], nor[3];

	sub_v3_v3v3(nor, bone->tail, bone->head);
	vec_roll_to_mat3(nor, 0.0f, mat);
	
	/* check the bone isn't aligned with the axis */
	if (!is_zero_v3(align_axis) && angle_v3v3(align_axis, mat[2]) > FLT_EPSILON) {
		float vec[3], align_axis_proj[3], roll;
		
		/* project the new_up_axis along the normal */
		project_v3_v3v3(vec, align_axis, nor);
		sub_v3_v3v3(align_axis_proj, align_axis, vec);
		
		if (axis_only) {
			if (angle_v3v3(align_axis_proj, mat[2]) > (float)(M_PI / 2.0)) {
				negate_v3(align_axis_proj);
			}
		}
		
		roll = angle_v3v3(align_axis_proj, mat[2]);
		
		cross_v3_v3v3(vec, mat[2], align_axis_proj);
		
		if (dot_v3v3(vec, nor) < 0) {
			roll = -roll;
		}
		
		return roll;
	}

	return 0.0f;
}

/**
 * finalize after accumulation.
 */
static void calc_tangent_ortho(float ts[3][3])
{
	float v_tan_a[3], v_tan_b[3];
	float t_vec_a[3], t_vec_b[3];

	normalize_v3(ts[2]);

	copy_v3_v3(v_tan_a, ts[0]);
	copy_v3_v3(v_tan_b, ts[1]);

	cross_v3_v3v3(ts[1], ts[2], v_tan_a);
	mul_v3_fl(ts[1], dot_v3v3(ts[1], v_tan_b) < 0.0f ? -1.0f : 1.0f);

	/* orthognalise tangent */
	mul_v3_v3fl(t_vec_a, ts[2], dot_v3v3(ts[2], v_tan_a));
	sub_v3_v3v3(ts[0], v_tan_a, t_vec_a);

	/* orthognalise bitangent */
	mul_v3_v3fl(t_vec_a, ts[2], dot_v3v3(ts[2], ts[1]));
	mul_v3_v3fl(t_vec_b, ts[0], dot_v3v3(ts[0], ts[1]) / dot_v3v3(v_tan_a, v_tan_a));
	sub_v3_v3(ts[1], t_vec_a);
	sub_v3_v3(ts[1], t_vec_b);

	normalize_v3(ts[0]);
	normalize_v3(ts[1]);
}

/* adjust bone roll to align Z axis with vector
 * vec is in local space and is normalized
 */
float ED_armature_ebone_roll_to_vector(const EditBone *bone, const float align_axis[3], const bool axis_only)
{
	float mat[3][3], nor[3];
	float vec[3], align_axis_proj[3], roll = 0.0f;

	BLI_ASSERT_UNIT_V3(align_axis);

	sub_v3_v3v3(nor, bone->tail, bone->head);

	/* If tail == head or the bone is aligned with the axis... */
	if (normalize_v3(nor) <= FLT_EPSILON || (fabsf(dot_v3v3(align_axis, nor)) >= (1.0f - FLT_EPSILON))) {
		return roll;
	}

	vec_roll_to_mat3_normalized(nor, 0.0f, mat);

	/* project the new_up_axis along the normal */
	project_v3_v3v3_normalized(vec, align_axis, nor);
	sub_v3_v3v3(align_axis_proj, align_axis, vec);

	if (axis_only) {
		if (angle_v3v3(align_axis_proj, mat[2]) > (float)(M_PI_2)) {
			negate_v3(align_axis_proj);
		}
	}

	roll = angle_v3v3(align_axis_proj, mat[2]);

	cross_v3_v3v3(vec, mat[2], align_axis_proj);

	if (dot_v3v3(vec, nor) < 0.0f) {
		return -roll;
	}
	return roll;
}

/**
 * Return true if the 2x axis are both aligned when projected into the view.
 * In this case, we can't usefully project the cursor onto the plane.
 */
static bool isPlaneProjectionViewAligned(const TransInfo *t)
{
  const float eps = 0.001f;
  const float *constraint_vector[2];
  int n = 0;
  for (int i = 0; i < 3; i++) {
    if (t->con.mode & (CON_AXIS0 << i)) {
      constraint_vector[n++] = t->con.mtx[i];
      if (n == 2) {
        break;
      }
    }
  }
  BLI_assert(n == 2);

  float view_to_plane[3], plane_normal[3];

  getViewVector(t, t->center_global, view_to_plane);

  cross_v3_v3v3(plane_normal, constraint_vector[0], constraint_vector[1]);
  normalize_v3(plane_normal);

  float factor = dot_v3v3(plane_normal, view_to_plane);
  return fabsf(factor) < eps;
}

static void testAxialSymmetry(BGraph *graph, BNode *root_node, BNode *node1, BNode *node2, BArc *arc1, BArc *arc2, float axis[3], float limit, int group)
{
	const float limit_sq = limit * limit;
	float nor[3], vec[3], p[3];

	sub_v3_v3v3(p, node1->p, root_node->p);
	cross_v3_v3v3(nor, p, axis);

	sub_v3_v3v3(p, root_node->p, node2->p);
	cross_v3_v3v3(vec, p, axis);
	add_v3_v3(vec, nor);
	
	cross_v3_v3v3(nor, vec, axis);
	
	if (fabsf(nor[0]) > fabsf(nor[1]) && fabsf(nor[0]) > fabsf(nor[2]) && nor[0] < 0) {
		negate_v3(nor);
	}
	else if (fabsf(nor[1]) > fabsf(nor[0]) && fabsf(nor[1]) > fabsf(nor[2]) && nor[1] < 0) {
		negate_v3(nor);
	}
	else if (fabsf(nor[2]) > fabsf(nor[1]) && fabsf(nor[2]) > fabsf(nor[0]) && nor[2] < 0) {
		negate_v3(nor);
	}
	
	/* mirror node2 along axis */
	copy_v3_v3(p, node2->p);
	BLI_mirrorAlongAxis(p, root_node->p, nor);
	
	/* check if it's within limit before continuing */
	if (len_squared_v3v3(node1->p, p) <= limit_sq) {
		/* mark node as symmetric physically */
		copy_v3_v3(root_node->symmetry_axis, nor);
		root_node->symmetry_flag |= SYM_PHYSICAL;
		root_node->symmetry_flag |= SYM_AXIAL;

		/* flag side on arcs */
		flagAxialSymmetry(root_node, node1, arc1, group);
		flagAxialSymmetry(root_node, node2, arc2, group);
		
		if (graph->axial_symmetry) {
			graph->axial_symmetry(root_node, node1, node2, arc1, arc2);
		}
	}
	else {
		/* NOT SYMMETRIC */
	}
}

static bool convex(const float p0[3], const float up[3], const float a[3], const float b[3])
{
	/* Vec3 va = a-p0, vb = b-p0; */
	float va[3], vb[3], tmp[3];
	sub_v3_v3v3(va, a, p0);
	sub_v3_v3v3(vb, b, p0);
	cross_v3_v3v3(tmp, va, vb);
	return dot_v3v3(up, tmp) >= 0;
}

static int convex(float *p0, float *up, float *a, float *b)
{
	// Vec3 va = a-p0, vb = b-p0;
	float va[3], vb[3], tmp[3];
	sub_v3_v3v3(va, a, p0);
	sub_v3_v3v3(vb, b, p0);
	cross_v3_v3v3(tmp, va, vb);
	return dot_v3v3(up, tmp) >= 0;
}

static bool depth_read_normal(
        const ViewContext *vc, const bglMats *mats, const int mval[2],
        float r_normal[3])
{
	/* pixels surrounding */
	bool  depths_valid[9] = {false};
	float coords[9][3] = {{0}};

	ARegion *ar = vc->ar;
	const ViewDepths *depths = vc->rv3d->depths;

	for (int x = 0, i = 0; x < 2; x++) {
		for (int y = 0; y < 2; y++) {
			const int mval_ofs[2] = {mval[0] + (x - 1), mval[1] + (y - 1)};

			const double depth = (double)depth_read_zbuf(vc, mval_ofs[0], mval_ofs[1]);
			if ((depth > depths->depth_range[0]) && (depth < depths->depth_range[1])) {
				if (depth_unproject(ar, mats, mval_ofs, depth, coords[i])) {
					depths_valid[i] = true;
				}
			}
			i++;
		}
	}

	const int edges[2][6][2] = {
	    /* x edges */
	    {{0, 1}, {1, 2},
	     {3, 4}, {4, 5},
	     {6, 7}, {7, 8}},
	    /* y edges */
	    {{0, 3}, {3, 6},
	     {1, 4}, {4, 7},
	     {2, 5}, {5, 8}},
	};

	float cross[2][3] = {{0.0f}};

	for (int i = 0; i < 6; i++) {
		for (int axis = 0; axis < 2; axis++) {
			if (depths_valid[edges[axis][i][0]] && depths_valid[edges[axis][i][1]]) {
				float delta[3];
				sub_v3_v3v3(delta, coords[edges[axis][i][0]], coords[edges[axis][i][1]]);
				add_v3_v3(cross[axis], delta);
			}
		}
	}

	cross_v3_v3v3(r_normal, cross[0], cross[1]);

	if (normalize_v3(r_normal) != 0.0f) {
		return true;
	}
	else {
		return false;
	}
}

static void flush_pixel(const MResolvePixelData *data, const int x, const int y)
{
	float st[2] = {(x + 0.5f) / data->w, (y + 0.5f) / data->h};
	float *st0, *st1, *st2;
	float *tang0, *tang1, *tang2;
	float no0[3], no1[3], no2[3];
	float fUV[2], from_tang[3][3], to_tang[3][3];
	float u, v, w, sign;
	int r;

	const int i0 = data->i0;
	const int i1 = data->i1;
	const int i2 = data->i2;

	st0 = data->mtface[data->face_index].uv[i0];
	st1 = data->mtface[data->face_index].uv[i1];
	st2 = data->mtface[data->face_index].uv[i2];

	multiresbake_get_normal(data, no0, data->face_index, i0);   /* can optimize these 3 into one call */
	multiresbake_get_normal(data, no1, data->face_index, i1);
	multiresbake_get_normal(data, no2, data->face_index, i2);

	resolve_tri_uv(fUV, st, st0, st1, st2);

	u = fUV[0];
	v = fUV[1];
	w = 1 - u - v;

	if (data->pvtangent) {
		tang0 = data->pvtangent + data->face_index * 16 + i0 * 4;
		tang1 = data->pvtangent + data->face_index * 16 + i1 * 4;
		tang2 = data->pvtangent + data->face_index * 16 + i2 * 4;

		/* the sign is the same at all face vertices for any non degenerate face.
		 * Just in case we clamp the interpolated value though. */
		sign = (tang0[3] * u + tang1[3] * v + tang2[3] * w) < 0 ? (-1.0f) : 1.0f;

		/* this sequence of math is designed specifically as is with great care
		 * to be compatible with our shader. Please don't change without good reason. */
		for (r = 0; r < 3; r++) {
			from_tang[0][r] = tang0[r] * u + tang1[r] * v + tang2[r] * w;
			from_tang[2][r] = no0[r] * u + no1[r] * v + no2[r] * w;
		}

		cross_v3_v3v3(from_tang[1], from_tang[2], from_tang[0]);  /* B = sign * cross(N, T)  */
		mul_v3_fl(from_tang[1], sign);
		invert_m3_m3(to_tang, from_tang);
	}
	else {
		zero_m3(to_tang);
	}

	data->pass_data(data->lores_dm, data->hires_dm, data->bake_data,
	                data->ibuf, data->face_index, data->lvl, st, to_tang, x, y);
}

/**
 * accumulate edge-vectors from all polys.
 */
static void calc_tangent_loop_accum(const float v_dir_prev[3],
                                    const float v_dir_next[3],
                                    float r_tspace[3][3])
{
  add_v3_v3v3(r_tspace[1], v_dir_prev, v_dir_next);

  if (compare_v3v3(v_dir_prev, v_dir_next, FLT_EPSILON * 10.0f) == false) {
    const float weight = fabsf(acosf(dot_v3v3(v_dir_next, v_dir_prev)));
    float nor[3];

    cross_v3_v3v3(nor, v_dir_prev, v_dir_next);
    normalize_v3(nor);

    cross_v3_v3v3(r_tspace[0], r_tspace[1], nor);

    mul_v3_fl(nor, weight);
    /* accumulate weighted normals */
    add_v3_v3(r_tspace[2], nor);
  }
}

/* builds an X and a Y axis from the given Z axis */
static void build_coordinate_frame(float axisX[3], float axisY[3], const float axisZ[3])
{
	const float faX = fabsf(axisZ[0]);
	const float faY = fabsf(axisZ[1]);
	const float faZ = fabsf(axisZ[2]);

	if (faX <= faY && faX <= faZ) {
		const float len = sqrtf(axisZ[1] * axisZ[1] + axisZ[2] * axisZ[2]);
		axisY[0] = 0; axisY[1] = axisZ[2] / len; axisY[2] = -axisZ[1] / len;
		cross_v3_v3v3(axisX, axisY, axisZ);
	}
	else if (faY <= faZ) {
		const float len = sqrtf(axisZ[0] * axisZ[0] + axisZ[2] * axisZ[2]);
		axisX[0] = axisZ[2] / len; axisX[1] = 0; axisX[2] = -axisZ[0] / len;
		cross_v3_v3v3(axisY, axisZ, axisX);
	}
	else {
		const float len = sqrtf(axisZ[0] * axisZ[0] + axisZ[1] * axisZ[1]);
		axisX[0] = axisZ[1] / len; axisX[1] = -axisZ[0] / len; axisX[2] = 0;
		cross_v3_v3v3(axisY, axisZ, axisX);
	}
}

static void do_kink_spiral_deform(ParticleKey *state, const float dir[3], const float kink[3],
                                  float time, float freq, float shape, float amplitude,
                                  const float spiral_start[3])
{
	float result[3];

	CLAMP(time, 0.f, 1.f);

	copy_v3_v3(result, state->co);

	{
		/* Creates a logarithmic spiral:
		 *   r(theta) = a * exp(b * theta)
		 *
		 * The "density" parameter b is defined by the shape parameter
		 * and goes up to the Golden Spiral for 1.0
		 * https://en.wikipedia.org/wiki/Golden_spiral
		 */
		const float b = shape * (1.0f + sqrtf(5.0f)) / (float)M_PI * 0.25f;
		/* angle of the spiral against the curve (rotated opposite to make a smooth transition) */
		const float start_angle = ((b != 0.0f) ? atanf(1.0f / b) :
		                           (float)-M_PI_2) + (b > 0.0f ? -(float)M_PI_2 : (float)M_PI_2);

		float spiral_axis[3], rot[3][3];
		float vec[3];

		float theta = freq * time * 2.0f * (float)M_PI;
		float radius = amplitude * expf(b * theta);

		/* a bit more intuitive than using negative frequency for this */
		if (amplitude < 0.0f)
			theta = -theta;

		cross_v3_v3v3(spiral_axis, dir, kink);
		normalize_v3(spiral_axis);

		mul_v3_v3fl(vec, kink, -radius);

		axis_angle_normalized_to_mat3(rot, spiral_axis, theta);
		mul_m3_v3(rot, vec);

		madd_v3_v3fl(vec, kink, amplitude);

		axis_angle_normalized_to_mat3(rot, spiral_axis, -start_angle);
		mul_m3_v3(rot, vec);

		add_v3_v3v3(result, spiral_start, vec);
	}

	copy_v3_v3(state->co, result);
}

static float BME_bevel_get_angle(BME_Mesh *UNUSED(bm), BME_Edge *e, BME_Vert *v) {
	BME_Vert *v1, *v2;
	BME_Loop *l1, *l2;
	float vec1[3], vec2[3], vec3[3], vec4[3];

	l1 = e->loop;
	l2 = e->loop->radial.next->data;
	if (l1->v == v) {
		v1 = l1->prev->v;
		v2 = l1->next->v;
	}
	else {
		v1 = l1->next->next->v;
		v2 = l1->v;
	}
	VECSUB(vec1,v1->co,v->co);
	VECSUB(vec2,v2->co,v->co);
	cross_v3_v3v3(vec3,vec1,vec2);

	l1 = l2;
	if (l1->v == v) {
		v1 = l1->prev->v;
		v2 = l1->next->v;
	}
	else {
		v1 = l1->next->next->v;
		v2 = l1->v;
	}
	VECSUB(vec1,v1->co,v->co);
	VECSUB(vec2,v2->co,v->co);
	cross_v3_v3v3(vec4,vec2,vec1);

	normalize_v3(vec3);
	normalize_v3(vec4);

	return dot_v3v3(vec3,vec4);
}

static float RotationBetween(TransInfo *t, const float p1[3], const float p2[3])
{
	float angle, start[3], end[3];

	sub_v3_v3v3(start, p1, t->center_global);
	sub_v3_v3v3(end,   p2, t->center_global);
		
	// Angle around a constraint axis (error prone, will need debug)
	if (t->con.applyRot != NULL && (t->con.mode & CON_APPLY)) {
		float axis[3], tmp[3];
		
		t->con.applyRot(t, NULL, axis, NULL);

		project_v3_v3v3(tmp, end, axis);
		sub_v3_v3v3(end, end, tmp);
		
		project_v3_v3v3(tmp, start, axis);
		sub_v3_v3v3(start, start, tmp);
		
		normalize_v3(end);
		normalize_v3(start);
		
		cross_v3_v3v3(tmp, start, end);
		
		if (dot_v3v3(tmp, axis) < 0.0f)
			angle = -acosf(dot_v3v3(start, end));
		else
			angle = acosf(dot_v3v3(start, end));
	}
	else {
		float mtx[3][3];
		
		copy_m3_m4(mtx, t->viewmat);

		mul_m3_v3(mtx, end);
		mul_m3_v3(mtx, start);
		
		angle = atan2f(start[1], start[0]) - atan2f(end[1], end[0]);
	}
	
	if (angle > (float)M_PI) {
		angle = angle - 2 * (float)M_PI;
	}
	else if (angle < -((float)M_PI)) {
		angle = 2.0f * (float)M_PI + angle;
	}
	
	return angle;
}

/* "Projects" a vector perpendicular to vec2 against vec1, such that
 * the projected vec1 + vec2 has a min distance of 1 from the "edge" defined by vec2.
 * note: the direction, is_forward, is used in conjunction with up_vec to determine
 * whether this is a convex or concave corner. If it is a concave corner, it will
 * be projected "backwards." If vec1 is before vec2, is_forward should be 0 (we are projecting backwards).
 * vec1 is the vector to project onto (expected to be normalized)
 * vec2 is the direction of projection (pointing away from vec1)
 * up_vec is used for orientation (expected to be normalized)
 * returns the length of the projected vector that lies along vec1 */
static float BME_bevel_project_vec(float *vec1, float *vec2, float *up_vec, int is_forward, BME_TransData_Head *UNUSED(td)) {
	float factor, vec3[3], tmp[3],c1,c2;

	cross_v3_v3v3(tmp,vec1,vec2);
	normalize_v3(tmp);
	factor = dot_v3v3(up_vec,tmp);
	if ((factor > 0 && is_forward) || (factor < 0 && !is_forward)) {
		cross_v3_v3v3(vec3,vec2,tmp); /* hmm, maybe up_vec should be used instead of tmp */
	}
	else {
		cross_v3_v3v3(vec3,tmp,vec2); /* hmm, maybe up_vec should be used instead of tmp */
	}
	normalize_v3(vec3);
	c1 = dot_v3v3(vec3,vec1);
	c2 = dot_v3v3(vec1,vec1);
	if (fabs(c1) < 0.000001f || fabs(c2) < 0.000001f) {
		factor = 0.0f;
	}
	else {
		factor = c2/c1;
	}

	return factor;
}

static float cotan_weight(float *v1, float *v2, float *v3)
{
	float a[3], b[3], c[3], clen;

	sub_v3_v3v3(a, v2, v1);
	sub_v3_v3v3(b, v3, v1);
	cross_v3_v3v3(c, a, b);

	clen = len_v3(c);

	if (clen == 0.0f)
		return 0.0f;
	
	return dot_v3v3(a, b)/clen;
}