Converting triangular polydata to quadratic


I have segmentation data that I use with MarchingCubes to make an isosurface. Then I use this surface with the probe filter to map intensity values on the mesh. The resulting geometry has triangular elements (see the figures below) but I want to convert it to quadratic elements. I tried vtkQuadricDecimation but it only converts very few elements and if I increase the decimation, the geometry messes up :/. The only remaining option is to convert cell by cell. Yet, the points are not organized in counter clock wise or clock wise necessarily(?). So, I am open to suggestions :))