Best way to solve for domain on which a function is real

Say I have a function like

1/(-54 + 1458 x + 6561 y + Sqrt[
   4 (-9 + 162 x)^3 + (-54 + 1458 x + 6561 y)^2])^(
 1/3) + (-54 + 1458 x + 6561 y + Sqrt[
   4 (-9 + 162 x)^3 + (-54 + 1458 x + 6561 y)^2])^(1/3)

what is the best way to solve for the domain in $(x,y)in mathbb{R}^2$ where the function is real?

I have been extracting the expressions in fractional powers and checking positivity of those but I am not sure that is the best approach (I guess it can miss solutions where two imaginary parts cancel). I also tried using explicitly solving Im[...]==0 (with or without applying ComplexExpand first) which definitely was not a success.

---

To add some background in case the problem is easier to solve starting from some different stage: Here I was trying to solve (3 q^2)/4 + (9 q^3)/4 + (3 q x)/2 - (3 y)/4 == 0 for real numbers x,y. If I add Reals to solve I get the following output:

{{q -> ConditionalExpression[Root[-y + 2 x #1 + #1^2 + 3 #1^3 &, 1], 
    x > 1/18 || (x < 1/
        18 && -2 + 54 x + 2 Sqrt[-(-1 + 18 x)^3] + 243 y > 0 && 
       2 - 54 x + 2 Sqrt[-(-1 + 18 x)^3] - 243 y > 0) || (x < 1/
        18 && -2 + 54 x + 2 Sqrt[-(-1 + 18 x)^3] + 243 y < 0) || (x < 
        1/18 && 2 - 54 x + 2 Sqrt[-(-1 + 18 x)^3] - 243 y < 
        0)]}, {q -> 
   ConditionalExpression[
    Root[-y + 2 x #1 + #1^2 + 3 #1^3 &, 
     2], -2 + 54 x + 2 Sqrt[-(-1 + 18 x)^3] + 243 y > 0 && 
     2 - 54 x + 2 Sqrt[-(-1 + 18 x)^3] - 243 y > 0 && 
     x < 1/18]}, {q -> 
   ConditionalExpression[
    Root[-y + 2 x #1 + #1^2 + 3 #1^3 &, 
     3], -2 + 54 x + 2 Sqrt[-(-1 + 18 x)^3] + 243 y > 0 && 
     2 - 54 x + 2 Sqrt[-(-1 + 18 x)^3] - 243 y > 0 && x < 1/18]}}

This does seem to have some useful information on the domain on (x,y) where there are real solutions. The root itself tells me little on the answer of course. (It basically defines the solution by the problem.) Are the three roots for this problem such as Root[-y + 2 x #1 + #1^2 + 3 #1^3 &, 1] the same ones as found for the (complex) case where I do not specify Reals in Solve or reduce. (Is the ordering the same?)

Clear["Global`*"]

eqn = (3 q^2)/4 + (9 q^3)/4 + (3 q x)/2 - (3 y)/4 == 0;

sol = Solve[eqn, q, Reals] // ToRadicals // Simplify;

dom = sol[[All, 1, -1, -1]]

{* {(18 x < 1 && (2 + 2 Sqrt[-(-1 + 18 x)^3] < 54 x + 243 y || 
      54 x + 2 Sqrt[-(-1 + 18 x)^3] + 243 y < 
       2 || (2 + 2 Sqrt[-(-1 + 18 x)^3] > 54 x + 243 y && 
        54 x + 2 Sqrt[-(-1 + 18 x)^3] + 243 y > 2))) || 18 x > 1, 
 54 x + 2 Sqrt[-(-1 + 18 x)^3] + 243 y > 2 && 
  2 + 2 Sqrt[-(-1 + 18 x)^3] > 54 x + 243 y && 18 x < 1, 
 54 x + 2 Sqrt[-(-1 + 18 x)^3] + 243 y > 2 && 
  2 + 2 Sqrt[-(-1 + 18 x)^3] > 54 x + 243 y && 18 x < 1} *)

domAll = (Or @@ dom) // FullSimplify

(* (18 x < 1 && (2 + 2 (1 - 18 x)^(3/2) < 54 x + 243 y || 
     2 (1 - 18 x)^(3/2) + 54 x + 243 y < 
      2 || (2 + 2 (1 - 18 x)^(3/2) > 54 x + 243 y && 
       2 (1 - 18 x)^(3/2) + 54 x + 243 y > 2))) || 18 x > 1 *)

Plot3D[Evaluate[q /. sol],
  {x, -2, 2}, {y, -2, 2},
  WorkingPrecision -> 15,
  AxesLabel -> (Style[#, 14, Bold] & /@ {x, y, q}),
  PlotStyle -> Opacity[0.75],
  PlotPoints -> 100,
  MaxRecursion -> 5,
  PlotLegends -> Automatic] // Quiet

f=(1/(-54 + 1458 x + 6561 y +Sqrt[4 (-9 + 162 x)^3 + (-54 + 1458 x + 6561 y)^2])^(1/3) +(-54 + 1458 x + 6561 y +Sqrt[4 (-9 + 162 x)^3 + (-54 + 1458 x + 6561 y)^2])^(1/3))
RegionPlot[Im[f] == 0, {x, -2, 2}, {y, -2, 2}]

Reduce[Im[f]==0,{x,y},Reals]]    
(*(x < 1/18 &&y >= -(2/243) (-1 + 27 x) +2/243 Sqrt[1 - 54 x + 972 x^2 -5832 x^3]) 
|| (x == 1/18 &&y > -(1/243)) || x > 1/18*)

f[x_, y_] := 
 1/(-54 + 1458 x + 6561 y + 
      Sqrt[4 (-9 + 162 x)^3 + (-54 + 1458 x + 6561 y)^2])^(1/
      3) + (-54 + 1458 x + 6561 y + 
     Sqrt[4 (-9 + 162 x)^3 + (-54 + 1458 x + 6561 y)^2])^(1/3)
Reduce[#, {x, y}, Reals]& @ FunctionDomain[f[x, y], {x, y}, Reals]
(*(x<1/18&&y≥-(2/243) (-1+27 x)+2/243 Sqrt[1-54 x+972 
x^2-5832 x^3])||(x[Equal]1/18&&y>-(1/243))||x>1/18*)