I am using Marching Squares to generate line segments (isosurface). As expected, there's a lot of triangles. What is the best way to accelerate ray-triangle intersection? I am looking into decimation, as well as a quadtree, but I have neither implemented yet. What do you recommend?

Decimation would speed up the build of the acceleration structure, but not so much the tracing. Also the decimation step is probably more expensive than the savings. So i would focus on making just an acceleration structure first.
Must popular for tracing is BVH or kd-tree, but octree works too. Maybe a simple regular grid would be enough as well. (Quadtree makes sense only if your triangles mostly form a plane.)
I recommend the BVH, because it has many applications, is flexible and can support dynamic scenes.

Thanks for the info!

I should have posted pictures sooner, sorry.

The first pic is linear separation. The next two pics use curved separation. This has to do with statistics.

P.S. There's also this pic, which shows how complicated the boundary can be.

Oh, this looks like a simple texture sample to the initial data would be already enough to answer ‘green or purple’? No ray tracing needed?

Otherwise i'd use a 2D regular grids at the same given resolution of the tessellation. Seems the max triangle count par cell is 4, so you could avoid a compaction step by allowing 4 triangle indices per cell.
Then do a 2D point in triangle test for those 4 trinagles to find the triangle containing the sample. Super fast and easy.

Please allow me to think upon what you've described. I really appreciate you sharing your perspective.

I can share code for a bicubic texture sampler on CPU (On GPU this can be emulated with 4 bilinear samples).
This gives smoother output than bilinear because it combines 3x3 cells. The separations within a single texel would be curved, not straight.

Though, bilinear samples would probably match your triangles exactly. (And i have code for that too.)

I added code to the ray-triangle intersection that I found on stackoverflow. Simply checking against the bounding square takes 1/5th of the time all around. Thank you for the insight. Works good.

float smallest_x = FLT_MAX;
float smallest_y = FLT_MAX;
float greatest_x = -FLT_MAX;
float greatest_y = -FLT_MAX;

if (v0.x > greatest_x) greatest_x = v0.x;
if (v1.x > greatest_x) greatest_x = v1.x;
if (v2.x > greatest_x) greatest_x = v2.x;
if (v0.x < smallest_x) smallest_x = v0.x;
if (v1.x < smallest_x) smallest_x = v1.x;
if (v2.x < smallest_x) smallest_x = v2.x;
if (v0.y > greatest_y) greatest_y = v0.y;
if (v1.y > greatest_y) greatest_y = v1.y;
if (v2.y > greatest_y) greatest_y = v2.y;
if (v0.y < smallest_y) smallest_y = v0.y;
if (v1.y < smallest_y) smallest_y = v1.y;
if (v2.y < smallest_y) smallest_y = v2.y;

// Point lies outside of the triangle bounding square. Abort early
if(orig.x < smallest_x || orig.x > greatest_x || 
   orig.y < smallest_y || orig.y > greatest_y)
{
		return false;
}

I have to resort to a hybrid approach when there are more than 2 class types. There's a lot of black space, which denotes disputed space. So, when on black space, one classifies based on nearest neighbour.

Here there are 3 class types.

Hmm. OK, so I found a minimal code for the quadtree. Except that the quadtree stores points. How do I use this code to store triangles? My googling is not really explaining this well.