ARM64 machine code, 52 bytes

Hex dump:

cb4413ea
d2a03feb
f2ffe02b
93cb416b
8b8be00c
8b2b802d
eb02019f
fa4331a2
54000042
a881348c
b6ffff2b
8b441140
d65f03c0

Try on godbolt

Register arguments match

size_t four_adjacent(uint64_t x,        /* x0 */
                     uint64_t y,        /* x1 */
                     uint64_t x_size,   /* x2 */
                     uint64_t y_size,   /* x3 */
                     uint64_t out[4][2] /* pointer in x4 */ );

Stores the output points as consecutive (x,y) pairs in the out array, and returns the number that were stored.

Disassembly:

   0:   cb4413ea        neg     x10, x4, lsr #4
   4:   d2a03feb        mov     x11, #0x1ff0000  
   8:   f2ffe02b        movk    x11, #0xff01, lsl #48
   c:   93cb416b        ror     x11, x11, #16
  10:   8b8be00c        add     x12, x0, x11, asr #56
  14:   8b2b802d        add     x13, x1, w11, sxtb
  18:   eb02019f        cmp     x12, x2
  1c:   fa4331a2        ccmp    x13, x3, #0x2, cc
  20:   54000042        b.cs    0x28 
  24:   a881348c        stp     x12, x13, [x4], #16
  28:   b6ffff2b        tbz     x11, #63, 0xc
  2c:   8b441140        add     x0, x10, x4, lsr #4
  30:   d65f03c0        ret

Some notes:

• x11 is initialized to the signed byte vector {0, 0, -1, 1, 0, 0, 1, -1} (LSB first), which takes two instructions. We rotate right by 16 bits on each iteration. The high byte is added to x, and the low byte is added to y, which is convenient to do with the shifted and sign-extended forms of add. We stop when we get back to the initial value, which is the only one in which the high bit is set (high byte = -1).

• We can check if the new x value is between 0 and x_max in one instruction by doing an unsigned compare. If the new value is negative, it behaves as a very large unsigned number, and so either out-of-bounds case results in the carry set (unsigned higher or same; recall the carry flag on subtract means "no borrow"). Using the conditional compare instruction ccmp, we can check both x and y against both bounds in two instructions.

• The post-increment addressing mode for stp lets us increment the output pointer for free. To recover the count of times it was incremented (i.e. (final - orig) / 16), we shift and negate the original pointer up front, saving it in x10, so that we can do the rest with a single shifted add at the end. (The negation is needed because shifted sub can only shift the subtrahend, not the minuend.) So it's effectively -(orig / 16) + (final / 16).

Dyalog APL, 19 bytes

Vector of (n, m) passed on the left; vector of (x, y) passed on the right.

{⍸1=+/¨2*⍨(⊂⌽⍺)-⍳⍵}
          (⊂⌽⍺)            Reverse and enclose the coordinate
                ⍳⍵         Get all indices for the width/height pair
               -           Subtract
       2*⍨                 Square
    +/¨                    Add over each
  1=                       Where the distance equals 1 are the four neighbors
 ⍸                         "where" to get the coordinates from the bitmask

Pip -xp, 18 17 bytes

SN_ADc=,2FIFLaCGb

Takes three command-line arguments: m, n, and a list containing x & y. Attempt This Online!

Explanation

Generates a list of all points in the valid rectangle and filters for the ones that are one unit away from the given point.

SN_ADc=,2FIFLaCGb
             aCGb  Generate a coordinate grid of a rows and b columns
           FL      Flatten once, giving a list of coord pairs
         FI        Filter on this function:
  _                  The function's argument
   AD                Absolute difference (element-wise) with
     c               The program's third argument
SN                   Sort the resulting pair in numerical ascending order
      =              Must equal
       ,2            Range(2), i.e. [0;1]

Wolfram Language (Mathematica), 42 37 bytes

pNorm[{##}-p]&~Array~#~Position~1&

Try it online!

Input [{x,y}][{m,n}]. 1-indexed, following Mathematica's convention.

Haskell, 62 bytes

Using

f m n a b = [(x,y)|x<-[0..m-1],y<-[0..n-1],(x-a)^2+(y-b)^2==1]

Try it online!

Boring approach: 81 bytes

f m n x y=filter (\(x,y)->x>=0&&y>=0&&x<m&&y<n) [(x-1,y),(x+1,y),(x,y-1),(x,y+1)]

Perl 6, 56 49 bytes

-7 bytes thanks to nwellnhof!

{grep 1>(*.reals Z/@^b).all>=0,($^a X+1,-1,i,-i)}

Try it online!

Removes the out of bounds elements by checking if when divided by the array bounds it is between 0 and 1. Takes input and output via complex numbers where the real part is the x coordinate and the imaginary is the y. You can extract these through the .im and .re functions.

Charcoal, 29 bytes

Ｊθη#ＦＩζＦＩε«Ｊικ¿№ＫＶ#⊞υ⟦ικ⟧»⎚Ｉυ

Try it online! Link is to verbose version of code. Takes inputs in the order x, y, width, height. Explanation:

Ｊθη#

Print a # at the provided position.

ＦＩζＦＩε«

Loop over the given rectangle.

Ｊικ

Jump to the current position.

¿№ＫＶ#⊞υ⟦ικ⟧

If there's an adjacent # then save the position.

»⎚Ｉυ

Output the discovered positions at the end of the loop.

Boring answer:

ＦＩζＦＩε¿⁼¹⁺↔⁻ιＩθ↔⁻κＩηＩ⟦ικ

Try it online! Link is to verbose version of code. Works by finding the adjacent positions mathematically.

Jelly,  13  12 bytes

2ḶṚƬNƬẎ+⁸%ƑƇ

A dyadic Link accepting a list of two (0-indexed) integers on the left, [row, column], and two integers on the right, [height, width], which yields a list of lists of integers, [[adjacent_row_1, adjacent_column_1], ...].

Try it online!

How?

2ḶṚƬNƬẎ+⁸%ƑƇ - Link: [row, column]; [height, width]   e.g. [3,2]; [5,3] (the "L" example)
2            - literal 2                                   2
 Ḷ           - lowered range                               [0,1]
   Ƭ         - collect up while distinct, applying:
  Ṛ          -   reverse                                   [[0,1],[1,0]]
     Ƭ       - collect up while distinct, applying:
    N        -   negate                                    [[[0,1],[1,0]],[[0,-1],[-1,0]]]
      Ẏ      - tighten                                     [[0,1],[1,0],[0,-1],[-1,0]]
        ⁸    - chain's left argument ([row, column])       [3,2]
       +     - add (vectorises)                            [[3,3],[4,2],[3,1],[2,2]]
           Ƈ - filter keep if:
          Ƒ  -   is invariant under:
         %   -     modulo ([height, width]) (vectorises)    [3,0] [4,2] [3,1] [2,2]
             - (...and [3,0] is not equal to [3,3] so ->)  [[4,2],[3,1],[2,2]]

JavaScript (ES6), 74 bytes

Boring approach.

(h,w,x,y)=>[x&&[x-1,y],~x+w&&[x+1,y],y&&[x,y-1],++y-h&&[x,y]].filter(_=>_)

Try it online!

JavaScript (Node.js), 74 bytes

Less boring but just as long. Takes input as ([h,w,x,y]).

a=>a.flatMap((_,d,[h,w,x,y])=>~(x+=--d%2)*~(y+=--d%2)&&x<w&y<h?[[x,y]]:[])

Try it online!

JavaScript (V8), 67 bytes

If all standard output methods were allowed, we could just print the valid coordinates with:

(h,w,x,y)=>{for(;h--;)for(X=w;X--;)(x-X)**2+(y-h)**2^1||print(X,h)}

Try it online!

Octave, 90 bytes

This uses a geometric approach: First we create an matrix of zeros of the desired size, and set a 1 to the desired location. Then we convolve with the kernel

[0, 1, 0]
[1, 0, 1]
[0, 1, 0]

which produces a new matrix of the same size with ones at the 4-neighbours of the original point. Then we find() the indices of the nonzero entries of this new matrix.

function [i,j]=f(s,a,b);z=zeros(s);z(a,b)=1;[i,j]=find(conv2(z,(v=[1;-1;1])*v'<0,'same'));

Try it online!

^{_{convolution is the key to success.}}

J, ^{30 29} 28 bytes

(([+.@#~&,1=|@-)j./)~j./&i./

Try it online!

How:

• Turn the right hand m x n arg into a grid of complex numbers j./&i./
• Same for left arg (our point) j./
• Create a mask showing where the distance between our point and the grid is exactly 1 1=|@-
• Use that to filter the grid, after flattening both #~&,
• Turn the result back into real points +.@

C# (Visual C# Interactive Compiler), 91 bytes

(a,b,x,y)=>new[]{(x+1,y),(x-1,y),(x,y+1),(x,y-1)}.Where(d=>((x,y)=d,p:x>=0&x<a&y>=0&y<b).p)

Try it online!

Alternatively:

(a,b,x,y)=>new[]{(x+1,y),(x-1,y),(x,y+1),(x,y-1)}.Where(d=>((x,y)=d).Item1>=0&x<a&y>=0&y<b)

Try it online!

Python 2, 66 bytes

lambda m,n,x,y:[(x-1,y),(x+1,y)][~x:m-x]+[(x,y-1),(x,y+1)][~y:n-y]

Try it online!

Lists the four neighbors, then uses list slicing to remove those that are out-of-bounds.

Python 2, 71 bytes

lambda m,n,x,y:[(k/n,k%n)for k in range(m*n)if(k/n-x)**2+(k%n-y)**2==1]

Try it online!

Instead of checking which of the four neighbors are in-bounds, we do it the slower way of checking all in-bounds points for those that are neighbors, that is have Euclidian distance exactly 1 from (x,y). We also use the classic div-mod trick to iterate over a grid, saving the need to write two loops like for i in range(m)for j in range(n).

I tried using complex arithmetic to write the distance condition, but it turned out longer to write abs((k/n-x)*1j+k%n-y)==1.

Python 2, 70 bytes

lambda m,n,x,y:[(x+t/3,y+t%3-1)for t in-2,0,2,4if m>x+t/3>=0<y+t%3<=n]

Try it online!