The hyperholomorphic Bergman kernel function \({}_\psi{\mathcal B}\) for a domain \(\Omega\subset\mathbf R ^4\) is introduced as a special quaternionic ``derivative'' of the Green function for \(\Omega\). It is shown that \({}_\psi{\mathcal B}\) is hyperholomorphic, Hermitian symmetric and reproduces hyperholomorphic functions. An integral representation of \({}_\psi{\mathcal B} \) is further obtained which expresses \({}_\psi{\mathcal B}\) as a sum of a function continuous on \(\overline\Omega\times\overline\Omega\) and an integral over \(\partial\Omega\) of a differential form on \(\mathbf R ^4\) which does not depend on \(\Omega\). As a consequence, the equality \({}^\psi B =I-\bar S\cdot S + C\) is proved for the Bergman projector \({}^\psi B \), where \(\bar S,S\) are certain singular integral operators and \(C\) is a compact operator. In more detail, let \(\psi=\{\psi^k\}_{k=0}^3\) be an orthonormal basis for the quaternions \(\mathbf H(\simeq\mathbf R ^4)\) over \({\mathbf R}\) and \(g\) the Green function for the domain \(\Omega\). A function \(f:\Omega\to\mathbf H\) is called (left-) \(\psi\)-hyperholomorphic if \({}^\psi D [f]:=\sum_{k=0}^3 \psi^k \frac {\partial f}{\partial x_k}=0\), and the hyperholomorphic Bergman kernel is defined by \({}_\psi{\mathcal B}(x,\xi)={}^{\bar\psi} D_x \cdot{}^\psi D _\xi[g](x,\xi)\). Let further \({}^\psi B [f](x)=\int_\Omega {}_\psi{\mathcal B} (x,\xi)f(\xi) d\xi\) (the Bergman operator), \({}^\psi T [f](x)=-\int_\Omega K_\psi(\tau-x) f(\tau) d\tau\), where \(K_\psi(x)=(2\pi^2|x|^4)^{-1} \sum_{k=0}^3 \bar\psi^k x_k\) is the fundamental solution of the operator \({}^\psi D \), and for \(\phi=\{\phi^k\}_{k=0}^3\) another orthonormal basis of \(\mathbf H \) such that \(\sum_{k=0}^3 \phi^k \cdot \psi^k =0\) introduce the operator \({}^{\phi,\psi}S [f](x)=2\pi^{-2}\int_\Omega (\sum_{k=0}^3 (\tau_k-x_k)\bar\psi^k) (\sum_{k=0}^3 (\tau_k-x_k)\bar\phi^k) |\tau-x|^{-6}f(\tau) d\tau\). Then \({}^\psi B \), \({}^\psi T \) and \({}^{\phi,\psi}S\) define bounded linear operators on the right \(\mathbf H \)-modules \(L_p(\Omega;\mathbf H )\) (\(p>1\)) and \(C^{0,\mu}(\Omega;\mathbf H )\) (\(0<\mu<1\)), \({}^\psi B ^2={}^\psi B \), \({}^\psi D \cdot{}^\psi T =I\), \({}^\phi D\cdot{}^{\bar\psi}T= {}^{\phi,\psi}S\), and \({}^\psi B =I-{}^{\bar\psi,\bar\phi}S \cdot{}^{\phi,\psi}S+{}^{\phi,\psi}C\) where \({}^{\phi,\psi}C\) is an integral operator whose kernel is continuous on \(\overline\Omega\times\overline\Omega\). An interpretation of these results in two-dimensional complex analysis is also given.