First recall that
$$
f(X|X \in [b-\epsilon,b+\epsilon])= \frac{1}{b+\epsilon - (b-\epsilon) } =\frac{1}{2\epsilon}
$$
and
$$
f(X) = f(X|X \in [b-\epsilon,b+\epsilon]) \times Pr(X \in [b-\epsilon,b+\epsilon])+ f(X|X \not\in [b-\epsilon,b+\epsilon]) \times Pr(X \not\in [b-\epsilon,b+\epsilon])
$$
but since $f(X|X \not\in [b-\epsilon,b+\epsilon]) =0$ the above simplifies to
$$
f(X) = f(X|X \in [b-\epsilon,b+\epsilon]) \times Pr(X \in [b-\epsilon,b+\epsilon])
$$
Note $ X \in [b-\epsilon,b+\epsilon] \equiv b \in [X-\epsilon,X+\epsilon]$ and that
$$
Pr(b \in [X-\epsilon,X+\epsilon]) = \int_{\max(b_l,X-\epsilon)}^{\min(b_h,X+\epsilon)}\frac{1}{b_h-b_l}db=\frac{\min(b_h,X+\epsilon)-\max(b_l,X-\epsilon)}{b_h-b_l}
$$
Thus;
$$
f(X)=\frac{\min(b_h,X+\epsilon)-\max(b_l,X-\epsilon)}{2\epsilon(b_h-b_l)}
$$
The above can be broken up into a piece-wise function differently depending on the values of $\epsilon,b_l,$ and $b_h$.
Simulation in R
To validate this answer I also ran a simulation where $\epsilon=2,b_l=0,$ and $b_h=1$. This results in the following pdf
$$
f(X) = \begin{cases}
\frac{X+\epsilon-b_l}{2\epsilon(b_h-b_l)} = \frac{X+2}{4}\;\;\mathrm{if}\;\;X \in [-2,-1)\\
\frac{b_h-b_l}{2\epsilon(b_h-b_l)} = \frac{1}{4}\;\;\mathrm{if}\;\;X \in [-1,2]\\
\frac{b_h+\epsilon-X}{2\epsilon(b_h-b_l)} = \frac{3-X}{4}\;\;\mathrm{if}\;\;X \in (2,3]\\
0\;\;\;\;\mathrm{otherwise}
\end{cases}
$$
The trapezoidal shape of $f(X)$ is validated with simulation and a kernel density estimator which I plot below
The code in R for the above plot is;
e=2
bl=0;
bh=1;
b=runif(1e5,bl,bh)
x=runif(1e5,min=b-e,max=b+e)
plot(density(x,from=-2,to=3),main="Density Estimate (KDE)",lwd=2)
abline(h=.25,col=2,lty=2,lwd=2)
abline(v=-1,col=4,lty=4,lwd=2)
abline(v=2,col=4,lty=4,lwd=2)