The quadtree structure is generated from the input heightmap. It is of constant depth, predetermined by memory and granularity requirements.
Does this mean I calculate every node and child nodes to the finest level at the setup phase?
Yes, for the "quadtree data" that contains only min/max height information.
I do not split and unsplit at runtime?
At runtime, you may load/generate heightmap "terrain data" for the chunks that need to be visible at a particular level of detail, and recycle chunks of heightmap that are no longer in view / no longer of the needed level of detail. So you do not need to have all of the terrain loaded at all times, only a minimal bit of metadata about its bounds.
What's the roughest level I start with?
You start at the root of the tree - the largest square that encompasses your terrain. For the case of a planet, you might have 6 roots: one for each face of a cube that you project into the spherical domain.
You continue down to the finest level of granularity your game will support. On planetary scales with procedural geometry, this might span a truly stupendous range, so you will likely want to handle this as an offline pre-process if you can, and set it up so you keep only a fixed number of levels "in core" at a time, writing out the results to the metadata structure as you go and discarding any bit of heightmap you're not actively investigating.
If you need this for a planet that is not known at build time, but generated at runtime, then you will need to use a different algorithm, or adapt it to skip this complete top-to-bottom scan, possibly using looser bounds as discussed below.
How do I get the minimum and maximum height from my Perlin noise height map for a node area?
Since the description of the algorithm asks you to evaluate all the way down to the finest granularity leaf node as part of the process, the trivial answer is: generate the height for every point in the terrain mesh in that leaf, and store the min and max that you generated.
Then you can throw out that leaf, and proceed to the next one. Once you have four leaves, you can take the minimum of the four leaf minima, and use that as the parent node's minimum. Then take the maximum of the four leaf maxima, and use that as the parent node's maximum. Then you can compress the leaf minima & maxima into the normalized range of the parent's min & max for compression purposes, as described in the paper.
The process then repeats with four parent nodes' minima & maxima informing the grandparent node's min & max, recursively all the way up the tree.
You may be able to approximate this without visiting every tilemap point:
Most noise-based terrain algorithms use a concept of "octaves", where you add together multiple noise samples over increasingly fine spatial frequencies and diminishing amplitude. In these cases, the first octave dominates the resulting height, and changes the least often over a range.
So, given a node in your quadtree, you can evaluate how many grid cells of the first octave it covers.
If it covers more than a threshold number of cells, you can conclude "the terrain probably hits both the global min and global max over this range" and just use your most conservative limits (ie. assuming every octave contributes its full amplitude at one point, and every octave hits its lowest possible value at another point).
If it covers a smaller number of grid cells, then you can walk through those cells and determine for each what minimum or maximum this first octave hits in each cell. Multiply that by the octave's amplitude, then make it a conservative estimate by adding the maximum possible impact of all the other octaves added together.
This gives you a ballpark for the min/max the noise function will take over that area, with a bounded amount of computation, even for planetary hexant scale quadtree nodes of indeterminate maximum detail.
How would CDLOD work, when I zoom into/out of a big planet?
Your computed distance from each node changes, leading you to either morph it toward and ultimately replace it with a coarser mesh, or subdivide it into finer meshes.
The main difference from the flat plane case is that your distance calculation will need to take into account the curvature of the planet when you're a long distance away, so that chunks that curve away from you are adequately down-rezzed.
You may also need to adjust the formula for morphing between detail levels in the vertex shader to apply the spherical projection after the morphing. I expect trying to apply the morphing in the curved form could get more complex. Doing it this way, the algorithm can basically pretend it's working on a flat plane, then you can just bend the result to fit your real application as a kind of post-process.