Perl 5.10, 173 132 bytes

@c=[2,$l=2];while($l/=2){@c=grep$$_[0]**2+$$_[1]**2<4?$n++*0:1,map{($x,$y)=@$_;[$x,$y],[$x-$l,$y],[$x-$l,$y-$l],[$x,$y-$l]}@c;say$n}

Try it online!

...outputs 1, 5, 14, 33, 71, 140, 274, 563, 1083, 2159, 4293, 8391 before it times out after one minute.

@c = ( [ 2, 2 ] );                 #starting corner
$l = 2;                            #starting length
while( $l /= 2 ){                  #reduce length by half each round
  @c =                             #new corners for next round
    grep $$_[0]**2 + $$_[1]**2 < 4 #pyt dude
         ? $n++ * 0                #count up if inside circle, dont keep corner
         : 1,                      #keep corner, but dont count up
    map {                          #for each current corner...
      ($x,$y) = @$_;               #...split current square into four:
      (
        [ $x,    $y    ],
        [ $x-$l, $y    ],
        [ $x-$l, $y-$l ],
        [ $x,    $y-$l ]
      )
    }
    @c;                            #current corners
  say $n                           #output count each round
}

APL (Dyalog Extended), 31 30 bytes

A train that takes \$n\$ as input and outputs the \$n\$th term.

(⊢⌿-1⊥3×1↓⌽)2{≢⍸⍵>|∘.⌾⍨⍳⍵}¨⍤*⍳


             {≢⍸⍵>|∘.⌾⍨⍳⍵}     ⍝ Number of points resides in the first quarter

            2             ¨⍤*⍳ ⍝ For radius of r in [1, 2, 4, 8, 16, ..., 2^n]

(⊢⌿-1⊥3×1↓⌽)                   ⍝ The last one minus three times of the sum of the rest

Let \$G_i\$ be the number of integer grid points that reside in the first quarter of a circle with radius of \$2^i\$ unit. (Not on the circle, nor on the axis.) Then the sequence \$ S \$ can be written as \$ S_i = \sum_{j=0}^i G_j-4G_{j-1} = G_i - 3\sum_{j=0}^{i-1}G_j\$. Since in the \$i\$th level, \$ 4G_{i-1} \$ of squares within \$ G_i \$ were already covered by levels \$j<i\$.

Try it online!

Try it online! 31 bytes

APL(Dyalog Unicode), 38 33 bytes ^SBCS

(⊢⌿-1⊥3×1↓⌽)2(≢∘⍸×⍨>⍳∘.+⍨⍤×⍳)¨⍤*⍳

(⊢⌿-1⊥3×1↓⌽)2(≢∘⍸×>∘.+⍨⍤×⍥⍳)⍨¨⍤*⍳

Try it on APLgolf!

Try it on APLgolf! 38 bytes

Javascript (Node.JS), 121 114 69 bytes

Versions:
Original: Try it online!
Version 2 (-7 bytes): (Thanks to emanresu) Try it online!
Version 3, with Recursion sorcery (-45 bytes) Try it online!

g=n=>n&&(F=(m,i=1<<m)=>i&&F(m,i-1)+(4**m-i*i)**.5|0)(n)+g(--n)-4*F(n)

Charcoal, 35 bytes

ＩΣＥＸ²⮌…·¹Ｎ×ΣＥιΣＥι‹⁺Ｘ⊕λ²Ｘ⊕ν²Ｘι²∨¬κ±³

Try it online! Link is to verbose version of code. Explanation:

    ²                               Literal integer `2`
   Ｘ                                Vectorised raised to power
      …·                            Inclusive range from
        ¹                           Literal integer `1` to
         Ｎ                          Input as an integer
     ⮌                              Reversed
  Ｅ                                 Map over values
             ι                      Current value
            Ｅ                       Map over implicit range
                ι                   Current value
               Ｅ                    Map over implicit range
                     λ              Inner value
                    ⊕               Incremented
                   Ｘ                Raised to power
                      ²             Literal integer `2`
                  ⁺                 Plus
                         ν          Innermost value
                        ⊕           Incremented
                       Ｘ            Raised to power
                          ²         Literal integer `2`
                 ‹                  Is less than
                            ι       Current value
                           Ｘ        Raised to power
                             ²      Literal integer `2`
              Σ                     Take the sum
           Σ                        Take the sum
          ×                         Multiplied by
                                κ   Current index
                               ¬    Is zero
                              ∨     Logical Or
                                  ³ Literal integer `3`
                                 ±  Negated
 Σ                                  Take the sum
Ｉ                                   Cast to string
                                    Implicitly print

I tried a more "graphical" version but it came out to be 56 54 bytes:

≔Ｘ²Ｎθ⊞υＥ³θＦυ¿›ΣＸ…鲦²Ｘ貫≔÷⌊ι²η¿ηＦ⁴⊞υＥι⎇＆κＸ²μλ⁻λη»Ｍ→Ｉⅈ

Try it online! Link is to verbose version of code. Explanation:

≔Ｘ²Ｎθ

Take 2 to the power of the input.

⊞υＥ³θＦυ

Construct a square of that size and loop over squares are they are generated.

¿›ΣＸ…鲦²Ｘ貫

If this square does not fit inside the circle, then...

≔÷⌊ι²η¿η

Halve the size of the square; if the square still has a size, then...

Ｆ⁴⊞υＥι⎇＆κＸ²μλ⁻λη

... push the four quarters of the original square to be checked.

»Ｍ→Ｉⅈ

Output the total number of squares found.

Vyxal, 16 bytes

ƛEDẊ²Ṡ√>∑;¦÷$4*ε

Try it Online!

Mostly a port of UnrelatedString's Jelly answer. There's got to be a way to save a byte on the last bit -1 thanks to UnrelatedString.

ƛ        ;       # Over 1...n
        ∑        # Count the number of 
  DẊ             # Pairs
 E               # Of 1...2^m
    ²Ṡ√          # Whose squared sums, sqrt'd
       >         # Are less than
 E               # 2^m
          ¦      # Take the cumulative sum
           ÷   ε # Subtract from the last item
            $4*  # The second-to-last item times 4

Python 3, 176 bytes

x=a=0;y=1
while 1:x*=4;z=len([*filter(lambda x:x[0]**2+x[1]**2<4,(lambda i:[(x*i,y*i)for x in range(1,int(2/i+1))for y in range(1,int(2/i+1))])(y))])-x;a+=z;print(a);y/=2;x+=z;

Try it online!

Outputs the sequence infinitely.

x and a are initialized to \$0\$, and y is initialized to \$1\$.

It does this next paragraph forever:

It multiplies x by 4. It then creates a list of the possible top-right corners of the squares with size y, then removes the vertices which do not satisfy the condition \$x^2+y^2<4\$. It sets z to the number of vertices which haven't been removed. minus x from it. z is added to a. It prints the current value of a, and then divides y by \$2\$, and adds z to x.

Jelly, 16 bytes

4*p`ḋ`<Ɗ)§Äṫ-ḅ-4

Try it online!

VERY slow--times out for n=6, brute forcing \$16^6 = 2^{24}\$ calculations.

  p`        For every pair of numbers from
4*          [1 .. 4^n],
    ḋ`      is its dot product with itself (sum of squares)
      <Ɗ    less than 4^n?

See below for further explanation:

Jelly, 18 17 bytes

4*ƲƇạ$)½Ḟ§Äṫ-ḅ-4

Try it online!

-1 thanks to Jonathan Allan

Takes n as input and outputs the nth element, computing the sequence as the cumulative sums of A177144 (number of squares added to areas not covered by squares in the previous iteration of A156790). Just barely beats the first-n-outputting 4*ƲƇạ$)½Ḟ§×4ạḊƲÄṖ and the somewhat faster R2*R²ạ²Ʋ€½Ḟ§Äṫ-ḅ-4.

       )             For each m in [1 .. n],
          §          sum
        ½            the square roots
         Ḟ           floored
4*   ạ               of 4^n minus
  ƲƇ $              every perfect square less than or equal to it.
            ṫ-       Last two
           Ä         cumulative sums
              ḅ-4    converted from base -4.

4*ƲƇạ$)½Ḟ§ḋṬo-3Ɗ and a few variants thereof tie, borrowing a little mathematical manipulation from David Chew's Python solution. Likewise 4*ṖṚƲT)½Ḟ§Äṫ-ḅ-4. A faster 16 sounds almost possible.

...Come to think of it, 2*R²ạ²Ʋ)½Ḟ§Äṫ-ḅ-4 also ties. Not sure how I missed that one, but I still like the suggested version too :P

Jelly, 37 bytes

2*R²ạ²Ʋ½ḞBUz0s"J’2*Ʋ$Ṛ_ṂḤŒHƲ€Ẏ+ɗ\ẎṂ€S

Try it online!

Utterly hopeless, but too fun an idea to pass up :P

Uses the binary expansion of each row to represent greedily filling it with powers of 2, then merges or explodes said powers of 2 into appropriate squares.

Jelly, 31 bytes

U»\U)
2*R²ạ²Ʋ½Rọ2Ç2*«ZgJ€ƲÇFݲS

Try it online!

Hacked together to try to do it graphically. Miraculously seems to work.

Python 3, 100 93 83 bytes

lambda m:sum((4**n-i*i)**.5//1*(n//m*4-3)for n in range(m+1)for i in range(1,2**n))

(Old: Try it online!)

-7 bytes thanks to emanresuA: Try it online!

-10 bytes thanks to GB and emanresuA: Try it online!

(Old: f calculates A156790 and the first lambda function takes the sum of the first n terms of A177144.)

Same as old, but f is combined with g by using a double for loop and (n//m*4-3), which in this context is 1 only for n==m and -3 otherwise.

Ruby, 93 ... 69 bytes

->n{-3*(1..n).sum{|x|n=(1..2**x).sum{|i|((4**x-i*i)**0.5).to_i}}+n*4}

Try it online!

Using the same approach as @David Chew in his Python solution, avoiding recursion makes it shorter and faster.

Wolfram Language (Mathematica), 125 bytes

use quick floorSqrt from this answer.

Try it online!

s@n_:=Floor[Sqrt[SetPrecision[n,RealExponent[1.001n]/2]]]
f@n_:=Sum[s[(4^n-i^2)],{i,1,2^n-1}]
g@n_:=f@n-3*Sum[f[i],{i,0,n-1}]

Scala 3, 173 169 bytes

Saved 4 bytes thanks to @Unrelated String

Golfed version. Attempt This Online!

val p=BigInt(2).pow
def f(n:Int)=(1to p(n).toInt-1).map(i=>math.floor(math.sqrt((p(2*n)-BigInt(i).pow(2)).toDouble)).toLong).sum
def g(n:Int)=f(n)-3*(0to n-1).map(f).sum

Ungolfed version. Attempt This Online!

object Main {
  def main(args: Array[String]): Unit = {
    for (i <- 0 to 20) {
      val result = g(i)
      println(s"g($i) = $result")
    }
  }

  def f(n: Int): Long = {
    var total = 0L
    val maxValue = BigInt(2).pow(2 * n)  // Equivalent to 4^n
    val limit = BigInt(2).pow(n).toInt   // Upper limit for i in the range
    for (i <- 1 until limit) {
      val diff = maxValue - BigInt(i).pow(2)
      val sqrtValue = math.sqrt(diff.toDouble)
      total += math.floor(sqrtValue).toLong
    }
    total
  }

  def g(n: Int): Long = {
    val sumPreviousF = (0 until n).map(f).sum
    f(n) - 3 * sumPreviousF
  }
}

ref Python code

def f(n):
    """
    Computes the sum over i from 1 to (2^n - 1) of the floor of sqrt(4^n - i^2).

    This function calculates the total integer parts of the square roots of (4^n - i^2)
    for each integer i in the range from 1 to 2^n - 1.
    """
    total = 0
    max_value = 2 ** (2 * n)  # This is 4^n, since 4^n = (2^n)^2
    limit = 2 ** n  # The upper limit of the range for i

    for i in range(1, limit):
        # Calculate the square root of (4^n - i^2)
        sqrt_value = (max_value - i * i) ** 0.5
        # Add the floor of the square root to the total sum
        total += int(sqrt_value // 1)
    return total

def g(n):
    """
    Computes the adjusted value g(n) by subtracting three times the sum of f(i) for i from 0 to n-1 from f(n).

    This function effectively isolates the new contributions at each level of n by removing
    overlapping counts from previous levels.
    """
    sum_previous_f = sum(f(i) for i in range(n))  # Sum of f(i) for all i from 0 to n-1
    return f(n) - 3 * sum_previous_f

# Compute and print g(n) for n from 0 to 20
for i in range(21):
    result = g(i)
    print(f"g({i}) = {result}")