In this answer, we obtain
$$ \bbox[padding:5px; border:1px solid navy; background-color:azure;]{ I(\omega) = -\omega \log \omega + \omega \left( 1 + \log(2\pi) \right)
+ 2\pi \operatorname{Im} \left[ \log\Gamma\left(\frac{i\omega}{2\pi}+\frac{1}{2}\right) \right] } $$
Step 1 - Preliminary. We introduce some key identities. First, we have
$$ \int_{0}^{\infty} \frac{\sin(at)\sin(bt)e^{-\varepsilon t}}{t} \, \mathrm{d}t = \frac{1}{2}\log\left| \frac{a+b + i\varepsilon}{a-b + i\varepsilon} \right|, \tag{1}\label{e:1} $$
where $a, b \in \mathbb{R}$ and $\varepsilon > 0$. Next, for $t \in \mathbb{R}$, we have
$$
\begin{align*}
\int_{0}^{\infty} \frac{\sin(xt)}{e^x + 1} \, \mathrm{d}x
&= \int_{0}^{\infty} \frac{\sin(xt) e^{-x}}{1 + e^{-x}} \, \mathrm{d}x \\
&= \sum_{n=1}^{\infty} (-1)^{n-1} \int_{0}^{\infty} \sin(xt) e^{-nx} \, \mathrm{d}x \\
&= \sum_{n=1}^{\infty} (-1)^{n-1} \frac{t}{x^2 + n^2} \\
&= \frac{1}{2} \left( \frac{1}{t} - \frac{\pi}{\sinh(\pi t)} \right).
\end{align*}
$$
Finally, a similar computation shows that for $\omega \in \mathbb{R}$ and $\varepsilon > 0$,
$$
\begin{align*}
\int_{0}^{\infty} \frac{\pi}{\sinh(\pi t)} \sin(\omega t)e^{-\varepsilon t} \, \mathrm{d}x
&= \operatorname{Im} \left[ \psi\left(\frac{\varepsilon+i\omega}{2\pi}+\frac{1}{2}\right) \right],
\end{align*}
$$
where $\psi$ is the digamma function. When obtaining this equality, the series formula for $\psi$ is particularly useful.
Step 2. By Dominated Convergene Theorem and Fubini's Theorem,
$$
\begin{align*}
I(\omega)
&= \lim_{\varepsilon \to 0^+} \int_{0}^{\infty} \frac{1}{e^x + 1} \log\left|\frac{x+\omega+i\varepsilon}{x-\omega+i\varepsilon}\right| \, \mathrm{d}x \\
&= \lim_{\varepsilon \to 0^+} \int_{0}^{\infty} \frac{1}{e^x + 1} \left( \int_{0}^{\infty} \frac{2\sin(xt)\sin(\omega t)e^{-\varepsilon t}}{t} \, \mathrm{d}t \right) \, \mathrm{d}x \\
&= \lim_{\varepsilon \to 0^+} \int_{0}^{\infty} \frac{\sin(\omega t)e^{-\varepsilon t}}{t} \left( \int_{0}^{\infty} \frac{2\sin(xt)}{e^x + 1} \, \mathrm{d}x \right) \, \mathrm{d}t \\
&= \lim_{\varepsilon \to 0^+} \int_{0}^{\infty} \frac{\sin(\omega t)e^{-\varepsilon t}}{t}\left( \frac{1}{t} - \frac{\pi}{\sinh(\pi t)} \right) \, \mathrm{d}t.
\end{align*}
$$
Now, let $J(\varepsilon)$ denote the integral in the last line. Then
$$
\begin{align*}
J'(\varepsilon)
&= -\int_{0}^{\infty} \sin(\omega t)\left( \frac{1}{t} - \frac{\pi}{\sinh(\pi t)} \right)e^{-\varepsilon t} \, \mathrm{d}t \\
&= -\arctan\left(\frac{\omega}{\varepsilon}\right) + \operatorname{Im} \left[ \psi\left(\frac{\varepsilon+i\omega}{2\pi}+\frac{1}{2}\right) \right]
\end{align*}
$$
Using the boundary condition $J(\infty) = 0$, it follows that
$$
\begin{align*}
I(\omega)
&= \lim_{R\to\infty} \left( -\int_{0}^{R} J(\varepsilon) \, \mathrm{d}\varepsilon \right) \\
&= \lim_{R\to\infty} \biggl( \int_{0}^{R} \arctan\left(\frac{\omega}{\varepsilon}\right) \, \mathrm{d}\varepsilon \\
&\qquad
- 2\pi \operatorname{Im} \left[ \log\Gamma\left(\frac{R+i\omega}{2\pi}+\frac{1}{2}\right) \right]
+ 2\pi \operatorname{Im} \left[ \log\Gamma\left(\frac{i\omega}{2\pi}+\frac{1}{2}\right) \right]
\biggr).
\end{align*}
$$
By the Stirling's approximation, we can check that
$$
- 2\pi \operatorname{Im} \left[ \log\Gamma\left(\frac{R+i\omega}{2\pi}+\frac{1}{2}\right) \right]
= -\omega \log R + \omega \log(2\pi) + \mathcal{O}(R^{-1}) $$
as $R \to \infty$. Plugging this back and substituting $\varepsilon = \omega x$,
$$
\begin{align*}
I(\omega)
&= \omega \lim_{R\to\infty} \biggl(
\int_{0}^{R/\omega} \arctan\left(\frac{1}{x}\right) \, \mathrm{d}x
- \log R \biggr) \\
&\qquad
+ \omega \log(2\pi)
+ 2\pi \operatorname{Im} \left[ \log\Gamma\left(\frac{i\omega}{2\pi}+\frac{1}{2}\right) \right] \\
&= -\omega \log \omega + \omega \biggl(
\int_{0}^{\infty} \left[ \arctan\left(\frac{1}{x}\right) - \frac{\mathbf{1}[x \geq 1]}{x} \right] \, \mathrm{d}x \biggr) \\
&\qquad
+ \omega \log(2\pi)
+ 2\pi \operatorname{Im} \left[ \log\Gamma\left(\frac{i\omega}{2\pi}+\frac{1}{2}\right) \right] \\
&= \bbox[padding:5px; border:1px solid navy; background-color:azure;]{ -\omega \log \omega + \omega \left( 1 + \log(2\pi) \right)
+ 2\pi \operatorname{Im} \left[ \log\Gamma\left(\frac{i\omega}{2\pi}+\frac{1}{2}\right) \right] }
\end{align*}
$$
Using this, we can obtain the series expansion of $I(\omega)$:
$$
\begin{align*}
I(\omega)
&= -\omega \log \omega
+ \omega \left( 1 - \gamma + \log\left(\frac{\pi}{2}\right) \right)
+ \frac{7\zeta(3)}{12\pi^2} \omega^3
- \frac{31\zeta(5)}{80\pi^4} \omega^5
+ \frac{127\zeta(7)}{448\pi^6} \omega^6
- \cdots.
\end{align*}
$$
Addendum. Here we provide more details in some non-trivial steps in the above computation.
Lemma 1. Let $a > 0$ and $b \in \mathbb{R}$. Then, as $ a \to \infty$,
$$ \operatorname{Im}\log\Gamma(a+ib) = b \log a + \mathcal{O}(a^{-1}). $$
Proof. By the Stirling's formula,
$$
\begin{align*}
\operatorname{Im}\log\Gamma(a+ib)
&= \operatorname{Im}\left[ (a + ib - \tfrac{1}{2})\log(a + ib) - (a + ib) + \tfrac{1}{2}\log(2\pi) + \mathcal{O}(a^{-1}) \right] \\
&= (a-\tfrac{1}{2})\arctan(b/a) + b\log\left|a + ib\right| - b + \mathcal{O}(a^{-1}) \\
&= b + b\log a - b + \mathcal{O}(a^{-1}) \\
&= b\log a + \mathcal{O}(a^{-1}).
\end{align*}
$$
Lemma 2. We have
$$ \int_{0}^{\infty} \left[ \arctan\left(\frac{1}{x}\right) - \frac{\mathbf{1}[x \geq 1]}{x} \right] \, \mathrm{d}x = 1. $$
Proof. By substituting $x \mapsto 1/x$, it suffices to prove that
$$ \int_{0}^{\infty} \left[ \frac{\arctan x}{x^2} - \frac{\mathbf{1}[x \leq 1]}{x} \right] \, \mathrm{d}x = 1. $$
Let $0 < \varepsilon < 1 < R$. Then
$$
\begin{align*}
\int_{\varepsilon}^{R} \left[ \frac{\arctan x}{x^2} - \frac{\mathbf{1}[x \leq 1]}{x} \right] \, \mathrm{d}x
&= \left[ -\frac{\arctan x}{x} \right]_{\varepsilon}^{R}
+ \int_{\varepsilon}^{R} \frac{1}{x(x^2 + 1)} \, \mathrm{d}x
- \int_{\varepsilon}^{1} \frac{1}{x} \, \mathrm{d}x \\
&= \frac{\arctan \varepsilon}{\varepsilon} - \frac{\arctan R}{R}
+ \int_{1}^{R} \frac{1}{x} \, \mathrm{d}x
- \int_{\varepsilon}^{R} \frac{x}{x^2 + 1} \, \mathrm{d}x \\
&= \frac{\arctan \varepsilon}{\varepsilon} - \frac{\arctan R}{R}
+ \log \frac{R}{\sqrt{R^2 + 1}} + \log \sqrt{\varepsilon^2 + 1}.
\end{align*}
$$
Now letting $R \to \infty$ and $\varepsilon \to 0^+$ proves the desired equality.
Lemma 3. For any $a, b \in \mathbb{R}$ and $\varepsilon > 0$, we have
$$ \int_{0}^{\infty} \frac{\sin(at)\sin(bt)e^{-\varepsilon t}}{t} \, \mathrm{d}t = \frac{1}{2}\log\left| \frac{a+b + i\varepsilon}{a-b + i\varepsilon} \right|. $$
Proof. This is a Frullani's integral in disguise:
$$
\begin{align*}
\int_{0}^{\infty} \frac{\sin(at)\sin(bt)e^{-\varepsilon t}}{t} \, \mathrm{d}t
&= \frac{1}{2} \int_{0}^{\infty} \frac{[\cos(at-bt) - \cos(at+bt)]e^{-\varepsilon t}}{t} \, \mathrm{d}t \\
&= \frac{1}{2} \operatorname{Re} \int_{0}^{\infty} \frac{e^{-(\varepsilon-ia+ib) t} - e^{-(\varepsilon-ia-ib) t}}{t} \, \mathrm{d}t \\
&= \frac{1}{2} \operatorname{Re} \int_{0}^{\infty} \int_{-b}^{b} -i e^{-(\varepsilon-ia-is) t} \, \mathrm{d}s\mathrm{d}t \\
&= \frac{1}{2} \operatorname{Re} \int_{-b}^{b} \frac{-i}{\varepsilon-ia-is} \, \mathrm{d}s \\
&= \frac{1}{2} \log \left| \frac{\varepsilon-ia-ib}{\varepsilon-ia+ib} \right|.
\end{align*}
$$
Now the desired identity follos by multiplying $i$ to both the numerator and denominator of the fraction inside the logarithm.
