First of all, I would like to state that handling the boundary conditions correctly is challenging, both from mathematical perspective and from an implementation perspective. Additionally, I think most of the papers don't explain the handling of boundary conditions sufficiently.
Furthermore, it is worth noticing that you have to apply boundary conditions to both the velocity field and the pressure field (and they are linked to each other). Additionally, you have to distinguish between is the pressure solve (or projection step) and the diffusion step, which must be treated separately.
In the very beginning you need to decide what kind of physical behavior you want to have at the boundaries. The standard example is a square surrounded by solid walls where no fluid can enter or leave. In this case you constrain the normal component of the velocity field to be zero ( a dirichlet boundary condition for the diffusion step). As we don't want changes for the normal component of the velocity in the projection step, we also need to constrain the pressure. If you check the pressure update routine (which makes the velocity field divergency free) you can see that the derivative of the pressure field in direction of the normal must be zero, I.e. a Neumann boundary condition for the pressure field.
Let's assume you are talking about the boundary conditions in the pressure solve/ projection step. Stam states that for pressure there is a Neumann boundary condition. As written above this condition implies that the velocity won't be modified at the last stage of the projection step (the pressure update). In order to achieve the dirichlet boundaries for the velocity as bridson states it, you require to have homogenous Neumann conditions. So for pressure, both authors have equivalent boundary conditions.
With these considerations you need to set up your system of linear equations. For this it is worth to note the equations for some cells (I.e. one row of the system of linear equations) and check how the matrix and the right hand side must be modified to account for the boundary conditions. In a 2d scenario we have a 5 point Laplace stencil that results in a 4 on the diagonal and -1 for each of the adjacent cells (times some constants). In cells that are adjacent to the boundary there is at least one value that is outside of the simulation domain and care has to be taken what we do with this value. For homogenous dirichlet conditions this value is simply zero and no changes are required for the matrix and the right hand side. For Neumann boundaries the constraint is that there is no change across the boundary, I.e. the value should be identical to the one of the cell. This reduces the value on the diagonal of the matrix by -1, I.e. 3 for a cell with one Neumann boundary. These considerations result in the linear system that you have to solve.
The set bnd method (that doesn't appear in the paper by the way) is a trick to make the computations without the necessity to modify the matrix. This trick only works with gauss seidel or Jacobi iterations (As commented on your other question). The set bnd method is part of the code that stam published. Personally I find it difficult to read and understand because of the variable names and the "b" that changes the behavior significantly and that seems to confuse you as well (see your other question).
Edit:
At first glance it looks like stam and bridson apply different boundary conditions, but they are actually the same. If we assume a solid body that is not moving and the velocity is zero before the pressure correction then the difference between them is zero. The term in parentheses in (4.23) will then be evaluated to -p_{i,j} which is exactly a homogeneous Neumann boundary condition for the pressure. I've got the second edition of bridsons book, where this is explained in equation 5.3 and at the bottom of the page in more detail.
Edit2: up to know we discussed the Neumann boundary conditions for the pressure. No slip/slip boundary conditions are applied to the velocity and for both (with non-moving walls) you need homogenized Neumann boundary conditions for the pressure. With stam's approach of discretization, I.e. storing all velocity components at the cell center it is easier to realize no slip boundary conditions as the boundary for all velocity components are at cell faces. So to realize no slip boundary conditions in the diffusion step (mathematically you can't realize no slip conditions with the Euler equations ( no diffusion)) you need to use homogeneous dirichlet boundary conditions for all velocity components.
Edit3: It is not straight forward to use the set_bnd method on a staggered grid for no slip boundary conditions. The storage for the tangential (to the boundary) components is not ideal.
Edit 4: A Neumann boundary condition for the pressure implies that there are no changes for the (normal) velocity at this boundary. If you subtract the pressure gradient at such a boundary, the pressure gradient is simply zero due to the boundary condition.
Edit 5: Slip / no slip boundary conditions (for the velocity) are not(!) realized in the pressure correction step. For non-moving solid walls there is always a homogeneous Neumann boundary condition for the pressure no matter if you impose a slip or no-slip boundary condition.
Edit 6: I don't "have" facebook ;)
There are more subtleties like singular systems for Neumann- only boundary conditions, central differences in staggered or standard grids, distance to the boundary (at cell border or at cell center outside three domain) etc., that I could address if you detail your question.
