PN-AEN indices algorithm

[morando]

I was following this article (http://developer.download.nvidia.com/whitepapers/2010/PN-AEN-Triangles-Whitepaper.pdf), and there is an detail explanation how to build index buffer suitable for this kind of mesh tessellation.
Since i dont get correct results i need someone to tell me if my algorithm is right (as described in article) or i need to debug some things elsewhere.

Implementation details:

1. Create an output IB that is 3 times the size of input IB.
2. For each input Triangle in IB, with indices i0, i1 and i2:
a. Write out an initial output entry of: i0, i1, i2, i0, i1, i1, i2, i2, i0, which sets edges to
initially be neighbors of themselves. This would produce identical results to PN Triangles.
b. Lookup the positions p0, p1, and p2, using i0, i1 and i2 to perform a lookup for
position of the associated vertex in VB.
c. Define 3 Edges, which consist of the two indices and two positions that make up
the corresponding Edge. An Edge should consist of the origin index, the
destination index, the origin position and the destination position.
d. For each edge, store the reverse of that edge in an easily searchable data structure
for the next step. The reference implementation uses an stdext::hash_map<Edge, Edge>
for this purpose. Reverse simply flips the sense of the edge (originating at
the destination position and index and heading to the origin position and index).
3. Walk the output index buffer (OB) constructed in step 2. For each patch of 9 indices:
a. For each Edge in the current Patch, perform a lookup into Edge->Edge mapping
created in step 2d.
b. If found, replace the current indices with the indices found in the map. Note that
two edges should be considered matching if their "from" and "to" indices match,
OR if their "from" and "to" positions match.
c. If not, continue to use the existing indices.

My code:

#define EPSILON 0.001f

// compare 2 floats
bool Equal(float f0, float f1)
{
	return (std::fabs(f0 - f1) <= EPSILON);
}

struct VERTEX
{
    Vec3 pos; // position vector
	Vec3 nrm; // normal vector
	Vec2 tex; // texture coordinate
};

void TestApp::calcAENIndices(const std::vector<WORD>& ind, const std::vector<VERTEX>& verts, std::vector<WORD>& out)
{
	// fill with initial values
	out.resize(ind.size() * 3);
	for (std::size_t i = 0; i < ind.size(); i += 3)
	{
		out[3 * i + 0] = ind[i + 0];
		out[3 * i + 1] = ind[i + 1];
		out[3 * i + 2] = ind[i + 2];

		out[3 * i + 3] = ind[i + 0];
		out[3 * i + 4] = ind[i + 1];
		out[3 * i + 5] = ind[i + 1];

		out[3 * i + 6] = ind[i + 2];
		out[3 * i + 7] = ind[i + 2];
		out[3 * i + 8] = ind[i + 0];
	}

	// define edge
	struct Edge
	{
		Vec3 p[2];
		WORD inx[2];

		bool operator == (const Edge& o) const
		{
			if (inx[0] == o.inx[0] && inx[1] == o.inx[1])
				return true;

			if (Equal(p[0].x, o.p[0].x) && Equal(p[0].y, o.p[0].y) && Equal(p[0].z, o.p[0].z))
			{
				if (Equal(p[1].x, o.p[1].x) && Equal(p[1].y, o.p[1].y) && Equal(p[1].z, o.p[1].z))
				{
					return true;
				}
			}

			return false;
		}
	};

	// fill reversed
	std::vector<Edge> edges(ind.size());
	for (std::size_t i = 0; i < ind.size(); i += 3)
	{
		edges[i + 0].p[1]   = verts[ind[i + 0]].pos;
		edges[i + 0].p[0]   = verts[ind[i + 1]].pos;
		edges[i + 0].inx[1] = ind[i + 0];
		edges[i + 0].inx[0] = ind[i + 1];

		edges[i + 1].p[1]   = verts[ind[i + 1]].pos;
		edges[i + 1].p[0]   = verts[ind[i + 2]].pos;
		edges[i + 1].inx[1] = ind[i + 1];
		edges[i + 1].inx[0] = ind[i + 2];

		edges[i + 2].p[1]   = verts[ind[i + 2]].pos;
		edges[i + 2].p[0]   = verts[ind[i + 0]].pos;
		edges[i + 2].inx[1] = ind[i + 2];
		edges[i + 2].inx[0] = ind[i + 0];
	}

	// compare Edge - Edge
	for (std::size_t i = 0, j = 0; i < out.size(); i += 9, j += 3)
	{
		for (std::size_t k = 0; k < 9; ++k)
		{
			std::size_t first  = k;
			std::size_t second = k + 1;
			if (second == 9)
				second = 0;

			Edge e;

			e.p[0]   = verts[out[i + first]].pos;
			e.p[1]   = verts[out[i + second]].pos;
			e.inx[0] = out[i + first];
			e.inx[1] = out[i + second];

			for (std::size_t m = 0; m < 3; ++m)
			{
				if (e == edges[j + m])
				{
					out[i + first]  = edges[j + m].inx[0];
					out[i + second] = edges[j + m].inx[1];
				}
			}
		}
	}
}

If you need more details, please let me know.
Thanks.