[Bug] Line Duplication Due To Line-Plane Intersection

[AhmedLilah]

When drawing a line and a plane that intersect each other, the line might show up as two lines. especially when there is a high contrast between the colors.

This is probably due to the method used to draw lines (i.e., two triangles to make a quad) and planes being the same.

Somehow, while sorting the triangles by the ZBuffer, we suspect that the sorting algorithm doesn't take into consideration the order in which lines or planes were created, leading to a reordering of the triangle that makes part of the plane be drawn between the line quad.

or,

That is because the line quad is always being drawn perpendicular to the camera, the plane is not at the same angle has parts that cross the line plane.

The part with the difference in Z-axis value from the camera would show in the middle of the line quad.

[AhmedLilah]

A quick fix to this problem that might lead to better 3D visualization of lines would be to represent the line as an octagonal prism.

That would get rid of the line splitting problem, as the sides will hide the edge of the plane.

It will also remove the need for the line quads to be drawn perpendicular to the camera.

[AhmedLilah]

@brenocq
I was trying to make a simple function for the prism plotting.

However, I encountered issues with plotting using PlotQuad.

I ignored the full implementation and tried to generate a stack of quads with different z values to plot them.

using this function

template<typename T> 
void PlotPrismLine(const char* label_id, const T* xs, const T* ys, const T* zs, int count, ImPlot3DLineFlags flags = 0, int offset = 0, int stride = sizeof(T)) {
	auto quad_xs = (float*)malloc(count * 8 * 4 * sizeof(float));
	auto quad_ys = (float*)malloc(count * 8 * 4 * sizeof(float));
	auto quad_zs = (float*)malloc(count * 8 * 4 * sizeof(float));

	for (int i = 0; i < count - 1; ++i) {
		constexpr auto sides = 8;
		constexpr auto radius = 1.0f;

		Vec3 octagonPoints1[sides];
		Vec3 octagonPoints2[sides];

		for (int j = 0; j < sides; ++j) {
			float const angle = 2.0f * 3.1415f * ((float)j / (float)sides);
			float const x = cosf(angle) * radius;
			float const y = sinf(angle) * radius;
			octagonPoints1[j] = Vec3(x, y, 0.0f) + Vec3(xs[i], ys[i], zs[i]);
			octagonPoints2[j] = Vec3(x, y, 0.0f) + Vec3(xs[i + 1], ys[i + 1], zs[i + 1]);
		}

		for (int j = 0; j < sides; ++j) {
			Vec3 vertices[4];
			vertices[0] = Vec3(0.0f + xs[i], 0.0f + ys[i], (float)j);
			vertices[1] = Vec3(1.0f + xs[i + 1], 0.0f + ys[i], (float)j);
			vertices[2] = Vec3(1.0f + xs[i + 1], 1.0f + ys[i + 1], (float)j);
			vertices[3] = Vec3(0.0f + xs[i], 1.0f + ys[i + 1], (float)j);

			for (int k = 0; k < 4; ++k) {
				quad_xs[(i * 8 * 4) + (j * 4) + (k)] = vertices[k].x;
				quad_ys[(i * 8 * 4) + (j * 4) + (k)] = vertices[k].y;
				quad_zs[(i * 8 * 4) + (j * 4) + (k)] = vertices[k].z;
			}

			if (i == 0) {
				printf("Quad %d:\n", j);
				for (int k = 0; k < 4; ++k)
					printf("  Vertex %d: (%.2f, %.2f, %.2f)\n", k, quad_xs[(i * 8 * 4) + (j * 4) + (k)], quad_ys[(i * 8 * 4) + (j * 4) + (k)], quad_zs[(i * 8 * 4) + (j * 4) + (k)]);
			}
		}
	}

	ImPlot3D::PlotQuad(label_id, quad_xs, quad_ys, quad_zs, (count - 1) * 8 * 4, ImPlot3DQuadFlags_None, 0, 1);

	free(quad_xs);
	free(quad_ys);
	free(quad_zs);
}

However, I keep getting these weird plots like this one.

While it should be something like this:

And the vertices are supposed to be done right as this is the output printed by the printf:

Quad 0:
  Vertex 0: (0.00, 0.00, 0.00)
  Vertex 1: (2.00, 0.00, 0.00)
  Vertex 2: (2.00, 2.00, 0.00)
  Vertex 3: (0.00, 2.00, 0.00)
Quad 1:
  Vertex 0: (0.00, 0.00, 1.00)
  Vertex 1: (2.00, 0.00, 1.00)
  Vertex 2: (2.00, 2.00, 1.00)
  Vertex 3: (0.00, 2.00, 1.00)
Quad 2:
  Vertex 0: (0.00, 0.00, 2.00)
  Vertex 1: (2.00, 0.00, 2.00)
  Vertex 2: (2.00, 2.00, 2.00)
  Vertex 3: (0.00, 2.00, 2.00)
Quad 3:
  Vertex 0: (0.00, 0.00, 3.00)
  Vertex 1: (2.00, 0.00, 3.00)
  Vertex 2: (2.00, 2.00, 3.00)
  Vertex 3: (0.00, 2.00, 3.00)
Quad 4:
  Vertex 0: (0.00, 0.00, 4.00)
  Vertex 1: (2.00, 0.00, 4.00)
  Vertex 2: (2.00, 2.00, 4.00)
  Vertex 3: (0.00, 2.00, 4.00)
Quad 5:
  Vertex 0: (0.00, 0.00, 5.00)
  Vertex 1: (2.00, 0.00, 5.00)
  Vertex 2: (2.00, 2.00, 5.00)
  Vertex 3: (0.00, 2.00, 5.00)
Quad 6:
  Vertex 0: (0.00, 0.00, 6.00)
  Vertex 1: (2.00, 0.00, 6.00)
  Vertex 2: (2.00, 2.00, 6.00)
  Vertex 3: (0.00, 2.00, 6.00)
Quad 7:
  Vertex 0: (0.00, 0.00, 7.00)
  Vertex 1: (2.00, 0.00, 7.00)
  Vertex 2: (2.00, 2.00, 7.00)
  Vertex 3: (0.00, 2.00, 7.00)

Could you think of any reason that makes it do that?

[AhmedLilah]

@brenocq

This is a snapshot from the current implementation of the PlotPrismLine.

There is room for improvement and a lot to ensure it's compliant with the API. However, for now, I think this would be a good state to examine to decide whether to proceed further with it.

[brenocq]

Hi @AhmedLilah! First of all, thank you for taking the lead on this issue. Just to make sure I understood this issue, is this happening in (1) and (2) because you are drawing very thick green lines and very thin blue lines exactly on top of it? It would be great if you could provide the code so I could also try to recreate the same problem here.

A quick fix to this problem that might lead to better 3D visualization of lines would be to represent the line as an octagonal prism.

I am a bit concerned about the number of triangles that the octagonal prism will generate, since I believe the octagonal prism may generate around 8x more triangles than the current approach. Ideally, the user should not need to think about the different ways of rendering lines to get the plot without artifacts.

This is a snapshot from the current implementation of the PlotPrismLine.

This looks way better than the first image! Can you show the side-by-side comparison of using octagonal prism lines and the quad lines?

[AhmedLilah]

Just to make sure I understood this issue, is this happening in (1) and (2) because you are drawing very thick green lines and very thin blue lines exactly on top of it? It would be great if you could provide the code so I could also try to recreate the same problem here.

No, In the case you're seeing there, the line in the middle of the green one is the black plane showing through the triangles that makes the line similar to the drawing I showed earlier, as far as I think.

There is no blue line on top of the green one.

I am a bit concerned about the number of triangles that the octagonal prism will generate, since I believe the octagonal prism may generate around 8x more triangles than the current approach. Ideally, the user should not need to think about the different ways of rendering lines to get the plot without artifacts.

I took that into consideration also and made the number of prism side parametrized depending on how much the user wants it to be to avoid the performance hit from a huge number of points.

This looks way better than the first image! Can you show the side-by-side comparison of using octagonal prism lines and the quad lines?

Sure, here you go.
Open the image in a large canvas to be able to see the difference clearly

This is the algorithm used to do the prism:

https://www.desmos.com/3d/az4uf8qcnz

As for the code here it's:

//-----------------------------------------------------------------------------
//[Section] PlotPrismLine
//-----------------------------------------------------------------------------

template<typename T> requires(std::is_floating_point_v<T>) 
struct Vec3 {
	T x, y, z;

	Vec3() = default;
	Vec3(T x, T y, T z) : x(x), y(y), z(z) {}

	Vec3 operator+(const Vec3& b) const { return {x + b.x, y + b.y, z + b.z}; }
	Vec3 operator-(const Vec3& b) const { return {x - b.x, y - b.y, z - b.z}; }
	Vec3 operator*(T s) const { return {x * s, y * s, z * s}; }
	Vec3 operator/(T s) const { return {x / s, y / s, z / s}; }

	T dot(const Vec3& b) const {
		return x * b.x + y * b.y + z * b.z;
	}

	Vec3 cross(const Vec3& b) const {
		return {
			y * b.z - z * b.y,
			z * b.x - x * b.z,
			x * b.y - y * b.x
		};
	}

	static constexpr T EPSILON = 1e-6f;

	T length() const {
		return sqrt(x * x + y * y + z * z); // or sqrtf if you prefer T
	}

	Vec3 normalized() const {
		T len = length();
		if (len > EPSILON) {
			return *this / len;
		}
		return Vec3{0.0f, 0.0f, 0.0f};
	}
};

template<typename T> requires(std::is_floating_point_v<T>)
void PlotPrismLine(const char* label_id, const T* xs, const T* ys, const T* zs, int count, int const sides, float weight, ImPlot3DLineFlags flags, int offset, int stride) {
	//const ImPlot3DNextItemData& n = GetItemData();
	const auto radius = weight * 0.0005f;

	auto quad_xs = (T*)malloc(count * sides * 4 * sizeof(T));
        auto quad_ys = (T*)malloc(count * sides * 4 * sizeof(T));
        auto quad_zs = (T*)malloc(count * sides * 4 * sizeof(T));

        auto octagonPoints1 = (Vec3<double> *) malloc(sides * sizeof(Vec3<double>));
        auto octagonPoints2 = (Vec3<double> *) malloc(sides * sizeof(Vec3<double>));

	for (int i = 0; i < count - 1; ++i) {
		Vec3<double> p0(xs[i], ys[i], zs[i]);
		Vec3<double> p1(xs[i + 1], ys[i + 1], zs[i + 1]);

		Vec3<double> axis = (p1 - p0);

		// Pick a helper vector that is not parallel to axis
		Vec3<double> helper;

		if (fabs(axis.x) <= fabs(axis.y) && fabs(axis.x) <= fabs(axis.z)) {
			helper = Vec3<double>(1, 0, 0);
		} else if (fabs(axis.y) <= fabs(axis.z)) {
			helper = Vec3<double>(0, 1, 0);
		} else {
			helper = Vec3<double>(0, 0, 1);
		}

		// Project helper onto the plane perpendicular to axis (Gram-Schmidt step)
		Vec3<double> right = (helper - axis * axis.dot(helper)).normalized();

		// Ensure orthonormality
		Vec3<double> up = axis.cross(right).normalized();

		for (int j = 0; j < sides; ++j) {
			T const angle = 2.0f * IM_PI * j / sides;
			T const x = cosf(angle) * radius;
			T const y = sinf(angle) * radius;

			// Generate points in the plane perpendicular to axis
			Vec3<double> offsetVec = right * x + up * y;

			octagonPoints1[j] = p0 + offsetVec;
			octagonPoints2[j] = p1 + offsetVec;
		}

		for (int j = 0; j < sides; ++j) {
			int j_next = (j + 1) % sides;

			Vec3<double> vertices[4];
			vertices[0] = octagonPoints1[j];
			vertices[1] = octagonPoints2[j];
			vertices[2] = octagonPoints2[j_next];
			vertices[3] = octagonPoints1[j_next];

			for (int k =0; k < 4; ++k) {
				quad_xs[(i * sides * 4) + (j * 4) + (k)] = vertices[k].x;
				quad_ys[(i * sides * 4) + (j * 4) + (k)] = vertices[k].y;
				quad_zs[(i * sides * 4) + (j * 4) + (k)] = vertices[k].z;
			}		
		}
	}
	
	PlotQuad(label_id, quad_xs, quad_ys, quad_zs, (count - 1) * sides * 4);

	free(quad_xs);
	free(quad_ys);
	free(quad_zs);
        free(octagonPoints1);
        free(octagonPoints2);
}

#define INSTANTIATE_MACRO(T)                                                                                                                         \
	template IMPLOT3D_API void PlotPrismLine<T>(const char* label_id, const T* xs, const T* ys, const T* zs, int count, int const sides, float weight, ImPlot3DLineFlags flags, int offset, int stride);
#define CALL_INSTANTIATE_FOR_FLOATING_TYPES \
    INSTANTIATE_MACRO(float);               \
    INSTANTIATE_MACRO(double);              \
    INSTANTIATE_MACRO(long double)

CALL_INSTANTIATE_FOR_FLOATING_TYPES
#undef INSTANTIATE_MACRO
#undef CALL_INSTANTIATE_FOR_FLOATING_TYPES

[AhmedLilah]

@brenocq

I submitted a PR from a fork since I couldn't do it directly because I'm not a contributor on the repo.

Some key points to take into consideration.

• This might not be the final implementation for this algorithm since it could use some optimization (some vectorization on modern CPUs, pre-baked LUTs for basic calculations regarding the small number of sides up to 8 or 16, etc)
• The allocation could be avoided for small numbers of sides also by just making a stack array that holds the data for sides up to 8 on a relatively small number of points.
• I believe the allocation using the malloc function isn't fully compatible with ImGui. We should probably use ImGuiMalloc.
• I'm still concerned bout the performance, as you mentioned, that's why I made the number of side parameterized.

[AhmedLilah]

@brenocq
Hey, Checking for any updates regarding this issue.