Python 2, 54 bytes

f=lambda r,x=0:r-x and-~((r*r-x*x)**.5%1>0)*4+f(r,x+1)

Try it online!

Less golfed (55 bytes) (TIO)

lambda r:8*r-4*sum((r*r-x*x)**.5%1==0for x in range(r))

This estimates the output as 8*r, then corrects for vertex crossings. The result is 8*r-g(r*r), where g(x) counts the number of ways to write x as a sum of two squares (except g(0)=0).

If the circle never went through any vertices, the number of cells touched would equal the number of edges crossed. The circle passes through 2*r vertical gridlines and 2*r horizontal gridlines, passing each one in both directions, for a total of 8*r.

But, each crossing at a vertex counts as two edge crossings while only entering one new cell. So, we compensate by subtracting the number of vertex crossings. This includes the points on axes like (r,0) as well as Pythagorean triples like (4,3) for r=5.

We count for a single quadrant the points (x,y) with x>=0 and y>0 with x*x+y*y==n, then multiply by 4. We do this by counting the numer of sqrt(r*r-x*x) that are whole number for x in the interval [0,r).

Mathematica, 48 bytes

4Count[Range@#~Tuples~2,l_/;Norm[l-1]<#<Norm@l]&

Looks at the first quadrant and counts the number of grid cells for which the input falls between the norms of the cell's lower left and upper right corners (multiplying the result by 4, of course).

Python 2, 72 bytes

lambda n:sum(0<n*n-x*x-y*y<2*(x-~y)for x in range(n)for y in range(n))*4

Try it online!

Jelly, 21 13 12 11 bytes

R²ạ²Æ²SạḤ×4

Try it online!

How it works

R²ạ²Æ²SạḤ×4  Main link. Argument: r

R            Range; yield [1, 2, ..., r].
 ²           Square; yield [1², 2², ..., r²].
   ²         Square; yield r².
  ạ          Absolute difference; yield [r²-1², r²-2², ..., r²-r²].
    Ʋ       Test if each of the differences is a perfect square.
      S      Sum, counting the number of perfect squares and thus the integer
             solutions of the equation x² + y² = r² with x > 0 and y ≥ 0.
        Ḥ    Un-halve; yield 2r.
       ạ     Subtract the result to the left from the result to the right.
         ×4  Multiply by 4.

Perl 6, 61 bytes

->\r{4*grep {my &n={[+] $_»²};n(1 X+$_)>r²>.&n},(^r X ^r)}

How it works

->\r{                                                    } # Lambda (accepts the radius).
                                                (^r X ^r)  # Pairs from (0,0) to (r-1,r-1),
                                                           #   representing the bottom-left
                                                           #   corners of all squares in
                                                           #   the top-right quadrant.
       grep {                                 }            # Filter the ones matching:
             my &n={[+] $_»²};                             #   Lambda to calculate the norm.
                              n(1 X+$_)>r²                 #   Top-right corner is outside,
                                          >.&n             #   and bottom-left is inside.
     4*                                                    # Return length of list times 4.

AWK, 90 bytes

{z=$1*$1
for(x=$1;x>=0;x--)for(y=0;y<=$1;y++){d=z-x*x-y*y
if(d>0&&d<2*(x+y)+2)c++}$0=4*c}1

Usage:

awk '{z=$1*$1
    for(x=$1;x>=0;x--)for(y=0;y<=$1;y++){d=z-x*x-y*y
    if(d>0&&d<2*(x+y)+2)c++}$0=4*c}1' <<< 5525

Just a simple search through quadrant 1 to find all boxes that will intersect the circle. Symmetry allows for the multiply by 4. Could go from -$1 to $1, but that would take a more bytes and be less efficient. Obviously this is not the most time efficient of algorithms, but it only takes about 16 seconds to run the 5525 case on my machine.

Haskell, 74 bytes

f n=sum[4|x<-[0..n],y<-[0..n],(1+n-x)^2+(1+n-y)^2>n^2,(n-x)^2+(n-y)^2<n^2]

Pretty straightforward, count the number of squares between (0,0) and (n,n) where the bottom left is inside the circle and the top right is outside the circle, then multiply by 4.

Pyth, 29 bytes

Lsm*ddb*4lf}*QQrhyTym+1dT^UQ2

Try it!

Explanation

Lsm*ddb*4lf}*QQrhyTym+1dT^UQ2  # implicit input: Q
Lsm*ddb                        # define norm function
 s                             # sum
  m   b                        #     map each coordinate to
   *dd                         #                            its square
                         ^UQ2  # cartesian square of [0, 1, ..., Q - 1]
                               #     -> list of coordinates of all relevant grid points
          f                    # filter the list of coordinates T where:
           }*QQ                # square of Q is in
               r               #     the range [
                hyT            #         1 + norm(T),
                               #                  ^ coordinate of lower left corner
                   ym+1dT      #         norm(map({add 1}, T))
                               #              ^^^^^^^^^^^^^^^ coordinate of upper right corner
                               #     ) <- half-open range
         l                     # size of the filtered list
                               #     -> number of passed-through squares in the first quadrant
       *4                      # multiply by 4
                               # implicit print

Batch, 147 bytes

@set/an=0,r=%1*%1
@for /l %%i in (0,1,%1)do @for /l %%j in (0,1,%1)do @set/a"i=%%i,j=%%j,a=i*i+j*j-r,i+=1,j+=1,a&=r-i*i-j*j,n-=a>>31<<2
@echo %n%

Somewhat inspired by the AWK and Haskell answers.

Bash + Unix utilities, 127 bytes

c()(d=$[(n/r+$1)**2+(n%r+$1)**2-r*r];((d))&&echo -n $[d<0])
r=$1
bc<<<`for((n=0;n<r*r;n++));{ c 0;c 1;echo;}|egrep -c 01\|10`*4

Try it online!

Just go through all the points in the first quadrant, count them up, and multiply by 4. It can be very slow, but it works.