On the evaluation of the integral $\int_{-\frac{b}{a}}^{\frac{1-b}{a}}\log\left(ax+b\right)\exp\left(-\frac{1}{2}x^2\right)\mathrm{d}x$.
First Approach: EDIT: I'm not too fond of this one, the second one seems better to me.I'll write $a = -\frac{\mu}{\sigma}$ and $b = \frac{1-\mu}{\sigma}$ and assume $\mu \neq 0$.
$$I = \frac{1}{\sqrt{2\pi}D}\int_{a}^{b}\log\left(\sigma u+\mu\right)\exp\left(-\frac{1}{2}u^2\right)\mathrm{d}u.$$
If $\mu \neq 0$, we can rewrite $\log(\sigma u + \mu) = \log(\frac{1}{\mu}) + \log(1+\frac{\sigma}{\mu}u)$. I'll denote $c := \frac{\sigma}{\mu}$. We can split this up into two integrals $I = I_1 + I_2$, where
$$ I_1 = \frac{1}{\sqrt{2\pi}D} \int_{a}^{b} -\log(\mu) \exp\left(-\frac{1}{2}u^2\right)\mathrm{d}u = \frac{-\log(\mu)}{D}\text{erf}(u)\Big|_{a}^{b},$$ and $$I_2 = \frac{1}{\sqrt{2\pi}D}\int_{a}^{b}\log\left(1+ cu\right)\exp\left(-\frac{1}{2}u^2\right)\mathrm{d}u.$$ Here I see two possiblites to continue: Plugging in the taylor expansion for $\log(1+cu)$ or the taylor expansion for $\exp(u)$. I'll begin with the former. Exchanging summation and integral (which i'm not yet sure is legal), and integrating by parts we obtain \begin{align} I_2 &= \frac{1}{\sqrt{2\pi}D}\int_{a}^{b}\sum_{n=1}^{\infty} \frac{(-1)^n}{n} (cu)^n\exp\left(-\frac{1}{2}u^2\right)\mathrm{d}u \\ &=\frac{1}{\sqrt{2\pi}D}\sum_{n=1}^{\infty}\int_{a}^{b} \frac{(-c)^n}{n} u^n\exp\left(-\frac{1}{2}u^2\right)\mathrm{d}u\\ &=\frac{1}{D}\sum_{n=1}^{\infty} (-c)^n \left[ \frac{u^n}{n}\text{erf}(u)\Big|_{a}^{b} - \int_{a}^{b}u^{n-1}\text{erf}(u) \mathrm{d}u\right]\\ \end{align}
Looking up the last integral in this integral table (page 5, number 7) , we obtain
\begin{align} I_2 = \frac{1}{D}\sum_{n=1}^{\infty} \frac{(-c)^n}{\sqrt{\pi}(n+1)} \left[ e^{-u^2}\sum_{k=0}^{l-1}\frac{\Gamma(\frac{n}{2}+1)}{\Gamma(\frac{n}{2}-k+1)} u^{n-2k} - (1-j)\Gamma\left(l+\frac{1}{2}\right)\text{erf}(u)\right]_{a}^{b} \end{align}
where $j = 0$ or $j=1$ such that $2l-j = n+1$. Now for the second approach, that yields a much nicer-looking sum.
Second Approach: Starting with
$$I = \frac{1}{\sqrt{2\pi\sigma^2}D}\int_{0}^{1}\log(x)\exp\left(-\frac{1}{2}\left(\frac{x-\mu}{\sigma}\right)^2\right)\mathrm{d}x$$
I'll expand the exponential function: \begin{align} \sqrt{2\pi\sigma}D\cdot I &= \int_{0}^{1} \log(x) \sum_{n=0}^{\infty} \frac{(-1)^n}{2^n \cdot n!}\left(\frac{x-\mu}{\sigma}\right)^{2n} \mathrm{d}x \\ &= \sum_{n=0}^{\infty} \frac{(-1)^n}{\sigma^{2n}\cdot 2^n \cdot n!} \int_{0}^{1} \left(x-\mu\right)^{2n} \log(x) \mathrm{d}x \\ &= \sum_{n=0}^{\infty} \frac{(-1)^n}{\sigma^{2n}\cdot 2^n \cdot n!} \sum_{k=0}^{2n} \binom{2n}{k}(-\mu)^{2n-k} \int_{0}^{1}x^{k} \log(x) \mathrm{d}x \\ \end{align}
Applying the formula $$\int x^{k}\ln x\,dx=\frac{x^{k+1}((k+1)\ln x-1)}{(k+1)^2}$$ we obtain
$$ I = \frac{1}{\sqrt{2\pi\sigma^2}D} \sum_{n=0}^{\infty} \frac{(-1)^{n+1}}{\sigma^{2n}\cdot 2^n \cdot n!} \sum_{k=0}^{2n} \binom{2n}{k}(-\mu)^{2n-k} \frac{1}{(k+1)^2} $$
It's not too hard to see that interchanging integral and sum is allowed using Fubini. We can find yet another representation for $I$: Using the binomial theorem, we see that
$$ \int_{0}^{1}\int_{0}^{1} (xy - \mu)^n \mathrm{d}x\mathrm{d}y = \sum_{k=0}^{2n} \binom{2n}{k}(-\mu)^{2n-k} \frac{1}{(k+1)^2}$$
Thus, we can write \begin{align} I &= \frac{1}{\sqrt{2\pi\sigma^2}D} \sum_{n=0}^{\infty} \frac{(-1)^{n+1}}{\sigma^{2n}\cdot 2^n \cdot n!} \int_{0}^{1}\int_{0}^{1} (xy - \mu)^n \mathrm{d}x\mathrm{d}y \\ &= -\int_{0}^{1} \int_{0}^{1} \frac{1}{\sqrt{2\pi\sigma^2}D} \exp\left(-\frac{1}{2}\left(\frac{xy - \mu}{\sigma}\right)^2\right) \mathrm{d}x\mathrm{d}y \end{align} Now we have a smooth integrand on our domain. Therefore, it should be possible to obtain good numerical results by applying high-degree 2-dimensional tensor product formulas of univariate Gauss-Legendre Quadrature formulas.
Hint:
$$\int_0^1\log(x)\exp\left(-\frac{1}{2}\left(\frac{x-\mu}\sigma\right)^2\right)\mathrm{d}x$$
can be seen as a convolution of $\log(x)$ restricted to the domain $(0,1]$, with a Gaussian, i.e. it is a blurred version of $\log(x)$.
As $\log(x)$ is integrable over $(0,1]$ (which yields $-1$), so is the given integral, for any $\mu,\sigma$.
As a consequence of the central limit theorem, convolution with a Gaussian can be approximated by a succession of convolutions with a uniform filter of appropriate width (the resulting variance is $nw^2/12$). http://oscar6echo.blogspot.be/2012/10/convolve-n-square-pulses-to-gaussian.html.
Also see Wells, W.M.: Efficient synthesis of Gaussian filters by cascaded uniform filters. IEEE Transactions on Pattern Analysis and Machine Intelligence 8(2), 234–239 (Mar 1986).
The uniform convolutions can be performed analytically, as
$$\int x^{k-1}\ln x\,dx=\frac{x^k(k\ln x-1)}{k^2}.$$
Stated differently, you can obtain good estimates of the integral by using the uniform B-spline approximation of the Gaussian, which leads to integrands of the form $\ln xP(x)$.
You find the relevant formula here: http://sepwww.stanford.edu/public/docs/sep105/sergey2/paper_html/node5.html#sergey2_eqn:expl
$$\beta_n(x)=\frac1{n!}\sum_{k=0}^n\binom{n+1}k(-1)^k\left(x+\frac{n+1}2-k\right)_+^n$$ where the subscript $_+$ denotes the restriction to positive arguments.
Notice that the approximation is only useful where the Gaussian is significantly nonzero.