  We study the interior regularity of solutions to the Dirichlet problem $Lu=g$
in $\Omega$, $u=0$ in $\R^n\setminus\Omega$, for anisotropic operators of
fractional type $$ Lu(x)= \int_{0}^{+\infty}\,d\rho
\int_{S^{n-1}}\,da(\omega)\,
  \frac{ 2u(x)-u(x+\rho\omega)-u(x-\rho\omega)}{\rho^{1+2s}}.$$ Here, $a$ is
any measure on~$S^{n-1}$ (a prototype example for~$L$ is given by the sum of
one-dimensional fractional Laplacians in fixed, given directions).
  When $a\in C^\infty(S^{n-1})$ and $g$ is $C^\infty(\Omega)$, solutions are
known to be $C^\infty$ inside~$\Omega$ (but not up to the boundary). However,
when $a$ is a general measure, or even when $a$ is $L^\infty(S^{n-1})$,
solutions are only known to be $C^{3s}$ inside $\Omega$.
  We prove here that, for general measures $a$, solutions are
$C^{1+3s-\epsilon}$ inside $\Omega$ for all $\epsilon>0$ whenever $\Omega$ is
convex. When $a\in L^{\infty}(S^{n-1})$, we show that the same holds in all
$C^{1,1}$ domains. In particular, solutions always possess a classical first
derivative.
  The assumptions on the domain are sharp, since if the domain is not convex
and the spectral measure is singular, we construct an explicit counterexample
for which $u$ is \emph{not} $C^{3s+\epsilon}$ for any $\epsilon>0$ -- even if
$g$ and $\Omega$ are $C^\infty$.
