Geometry problem to prove Euler line of a triangle $I$ is the incenter of $\triangle{ABC}$ and $D$ is the tangent point of the incircle with side $BC$. Points $E$ and $F$ are on side $BC$ so that $\angle{BAE}=\angle{ACB}, \angle{CAF}=\angle{ABC}$. $G$ and $H$ are incenters of $\triangle{ABE}$ and $\triangle{ACF}$ respectively. Prove that $AD$ is the Euler Line of $\triangle{GHI}$.
I tried to find related material to solve this problem but in vain. Any hints or helpful suggestion? Thanks!

With hint to use barycentric coordinates, here is what I have got:
Let $s=\dfrac{a+b+c}{2}, t=\sin A+\sin B+\sin C$, then
$A(1:0:0)\\
B(0:1:0)\\
C(0:0:1)\\
I(\dfrac{a}{2s}:\dfrac{b}{2s}:\dfrac{c}{2s})\\
D(0:\dfrac{s-c}{a}:\dfrac{s-b}{a})\\
E(0:1-\dfrac{c\sin C}{a\sin A}:\dfrac{c\sin C}{a\sin A})\\
F(0:\dfrac{b\sin B}{a\sin A}:1-\dfrac{b\sin B}{a\sin A})\\
G(\dfrac{\sin C}{t}:1-\dfrac{(a+c)\sin C}{at}:\dfrac{c\sin C}{at})\\
H(\dfrac{\sin B}{t}:\dfrac{b\sin B}{at}:1-\dfrac{(a+b)\sin B}{at})
$
 A: $\begin{array}{} ∠(ACB)=∠(BAE)=γ & ∠(ABC)=∠(CAF)=β & ∠(AEK)=β+γ \end{array}$
$\begin{array}{} A=\left( c\,cos(β),c\,sin(β)  \right)  & B=(0,0) & C=(a,0) \end{array}$
$\begin{array}{} ΔABC & cos(β)=\frac{a^2-b^2+c^2}{2ac} & cos(γ)=\frac{a^2+b^2-c^2}{2ab} \\ sin(β)=\frac{2S}{ac} & sin(β)=\frac{2S}{ab} & S=Δ(ABC)=\sqrt{s(s-a)(s-b)(s-c)} & s=\frac{a+b+c}{2}\end{array}$
$\begin{array}{} x_{A}=\frac{a^2-b^2+c^2}{2a} & y_{A}=\frac{2S}{a}  \end{array}$
$\begin{array}{} ΔAEK & AK=c\,sin(β) & EK=\frac{AK}{tan(β+γ)} \\ tan(β+γ)=\frac{cos(β)sin(β)+cos(γ)sin(γ)}{cos(β)^2+cos(γ)^2-1}=\frac{4S}{a^2-b^2-c^2} & EK=\frac{a^2-b^2-c^2}{2a} & BE=BK-EK \\ \end{array}$
$\begin{array}{} BE=\frac{c^2}{a} & E=\left( \frac{c^2}{a},0  \right)  & ∠(AEK)=∠(AFK) & AE=AF & AE=\frac{bc}{a} \end{array}$
$\begin{array}{} BF=BK+EK & BF=\frac{a^2-b^2}{a} & F=\left(  \frac{a^2-b^2}{a},0  \right)  \end{array}$
$\begin{array}{} \text{G is the incenter ΔABE } & x_{G}=\frac{x_{A}·BE+x_{B}·AE+x_{E}·c}{BE+AE+c} & x_{G}=\frac{c(a-b+c)}{2a} & y_{G}=\frac{2Sc}{a(a+b+c)} \end{array}$
$\begin{array}{} \text{H is the incenter ΔACF} & FC=a-BF & FC=\frac{b^2}{a} \\ x_{H}=\frac{2a^2-b(a+b-c)}{2a} & y_{H}=\frac{2Sb}{a(a+b+c)} & \, \\ \text{I is the incenter ΔABC} & x_{I}=\frac{a-b+c}{2} & y_{I}=\frac{2S}{a+b+c}\end{array}$
$l.x+m.y+n=0$ is the Euler line of $ΔGHI$. $l$ and $m$ are the coefficients of the Euler line in Cartesian coordinates. With CAS we simplify the matrix and get:
$\begin{array}{} l=\left| \begin{array}{} x_{G} & y_{G}^2+2y_{I}\,y_{H}+3x_{G}^2 & 1 \\ x_{I} & y_{I}^2+2y_{G}\,y_{H}+3x_{I}^2 & 1\\ x_{H} & y_{H}^2+2y_{G}\,y_{I}+3x_{H}^2 & 1 \\ \end{array} \right|  & m=\left| \begin{array}{} y_{G} & x_{G}^2+2x_{I}\,x_{H}+3y_{G}^2 & 1 \\ y_{I} & x_{I}^2+2x_{G}\,x_{H}+3y_{I}^2 & 1\\ y_{H} & x_{H}^2+2x_{G}\,x_{I}+3y_{H}^2 & 1 \\ \end{array} \right|   \end{array}$
$\begin{array}{} l=\frac{(a-c)(a-b)(a^2-b^2-c^2+2bc)(a^2-b^2-c^2+bc)}{a^2(a+b+c)} &  \\ m=\frac{4S(a-b)(a-c)(b-c)(a^2-b^2-c^2+bc)}{a^2(a+b+c)^2} &  \\ \end{array}$
$\begin{array}{} m_{Euler}=\frac{-l}{m} & m_{Euler}=\frac{-1}{4}\frac{(a+b+c)(a-b+c)(a+b-c)}{(b-c)S} \\ (a+b+c)(a-b+c)(a+b-c)=\frac{16S^2}{-a+b+c} & m_{Euler}=\frac{4S}{(a-b-c)(b-c)} \\ \end{array}$
$\begin{array}{} \text{J is the centroid of ΔGIH} & x_{J}=\frac{3a^2-2ab-b^2+2ac+c^2}{6a} & y_{J}=\frac{2S}{3a} \end{array}$
$\begin{array}{} \text{Euler line} & y-y_{J}=m_{Euler}(x-x_{J}) & \frac{y-y_{J}}{x-x_{J}}=m_{Euler} & \text{(x,y) is any point on the Euler line }\end{array}$
$\begin{array}{} \text{checking if the points belong to the line} & \frac{y_{A}-y_{J}}{x_{A}-x_{J}}-m_{Euler} =0& A=\left( \frac{a^2-b^2+c^2}{2a},\frac{2S}{a}  \right)  & \text{true} \\ \frac{y_{D}-y_{J}}{x_{D}-x_{J}}-m_{Euler}=0 & D=\left( \frac{a-b+c}{2},0  \right)  & \text{true}  \\ \end{array}$
A: Oh, this is a nice exercise, i will remember it for didactic purposes. Below, all details will be given, and the computations are done explicitly.
As references i am using
bary-full and bary-short. Notations like $X(2)$ for the centroid are as in ETC.
Fix a (non-degenerated) triangle $\Delta ABC$ with sides $a,b,c$.
Below, $s=(a+b+c)/2$ is the semiperimeter.
Each point in plane is written uniquely as $P=xA+yB+zC$ with $x+y+z=1$.
This makes sense by replacing $P,A,B,C$ by their affixes in $\Bbb C$
or cartesian coordinates in $\Bbb R^2$. Also, it makes sense in vectorial notation after fixing some / any origin $O$, and plugging it in in the form
$OP=xOA+yOB+zOC$ (with vectors not explicitly decorated with arrows).
Then in (normalized) notation $P=(x,y,z)$. Sometimes, we use also $x,y,z$ with non-zero sum to describe a point. Then $P=[x:y:z]$ means that $P=(x/(x+y+z),\ y/(x+y+z),\ z/(x+y+z))$. Abusively, i will also write
$\displaystyle P=\frac 1{x+y+z}(x,y,z)$ in such situations.
Note that $\Delta EBA\sim \Delta ABC$, so
$\displaystyle
\frac{BE}{AB}=\frac{AB}{BC}$, so
$\displaystyle BE=\frac{c^2}a$, and
$\displaystyle EC=a-\frac{c^2}a$, which gives
$\displaystyle
BE:EC=\frac{c^2}{a^2}:\left(1-\frac{c^2}{a^2}\right)=c^2:(a^2-c^2)$.
At this point, we already know the coordinates of $E$. Let us make the full list so far:
$$
\begin{aligned}
A &=(1,0,0)\\
B &=(0,1,0)\\
C &=(0,0,1)\\
I &= [a:b:c]=\frac 1{a+b+c}(a,b,c)\\
D &= [0:s-c:s-b]\\[3mm]
&\qquad\text{ so $AD$ has the equation }
  0 = \begin{vmatrix}x&y&z\\1&0&0\\0&s-c&s-b\end{vmatrix}
\\
&\qquad\text{ which is simplified }
 y(a-b+c) = z(a+b-c)
\\[3mm]
E &= [0:a^2-c^2:c^2]=\left(0\ :\ 1-\frac{c^2}{a^2}\ :\ \frac{c^2}{a^2}\right)\\
F &= [0:b^2:a^2-b^2]=\left(0\ :\ \frac{b^2}{a^2}\ :\ 1-\frac{c^2}{a^2}\right)\\
G &= \frac 1{a+b+c}(aE+bB+cA)
\\
&\qquad\qquad\text{ written parallel to }
I=\frac 1{a+b+c}(aA+bB+cC)
\text{ since }\Delta EBA\sim\Delta ABC
\\
&= \frac 1{a+b+c}\left(a\left(1-\frac{c^2}{a^2}\right)B +a\frac{c^2}{a^2}C+bB+cA\right)\qquad\text{ so}
\\
G &
=\left[c\ :\ a+b-\frac{c^2}a\ :\ \frac{c^2}a\right]
=\frac 1{a+b+c}\left(c\ ,\ a+b-\frac{c^2}a\ ,\ \frac{c^2}a\right)
\\
H &=\left[b\ :\ \frac{b^2}a\ :\ a+c-\frac{b^2}a\right]
=\frac 1{a+b+c}\left(b\ ,\ \frac{b^2}a\ ,\ a+c-\frac{b^2}a\right)
\\[3mm]
X(2)_{\Delta GHI}
&=\frac 13(G+H+I)\\
&=\frac 1{3(a+b+c)}
\left( a+b+c\ ,\ a+2b+\frac 1a(b^2-c^2)\ ,\ a+2c+\frac 1a(c^2-b^2)\right)
\\
&=\Big[a(a+b+c)\ :\ (a+b)^2-c^2\ :\ (a+c)^2-b^2\Big]
\\
&=[a\ :\ a+b-c\ :\ a-b+c]\ .
\end{aligned}
$$
Note that this centroid $X(2)$ of $\Delta GHI$ easily satisfies the equation of $AD$ computed above.
We need explicitly one more point on the Euler line of $\Delta GHI$, then chek again the equation of $AD$ for it, and we are done. I will compute its circumcenter $X(3)_{\Delta GHI}$ by intersecting its side bisectors. In this case, a picture is our friend:

Let $(x,y,z)$ be the coordinates of $X(3)_{\Delta GHI}$.
Then we  build the displacement vectors for $IH$ and for the vector from the mid point $H'$ of $IH$ to $X(3)$. They are
$$
\begin{aligned}
\overrightarrow{IH}
&= H-I \\
&= \frac 1{a+b+c}\left(b-a\ ,\ \frac{b^2}a-b\ ,\ a-\frac{b^2}a\right)\\
&\sim (a(b-a)\ ,\ b(b-a)\ ,\ -(b-a)(b+a))\\
&\sim (a\ ,\ b\ ,\ -a-b)\\
\ ,
\\
\overrightarrow{H'X(3)}
&= X(3)-\frac 12(H+I) \\
&= (x,y,z)-
\frac 1{2(a+b+c)}\left(b+a\ ,\ \frac{b^2}a+b\ ,\ a+2c-\frac{b^2}a\right)
\\
&\sim (4asx,\ 4asy,\ 4asz) -(a^2+ab, \ b^2+ab\ ,\ a^2+2ac-b^2)
\ .
\end{aligned}
$$
We write the equation from EFFT (from bary-short) and get the point:
$$
X(3)
=\Big[ a(b^2+c^2-a^2)\ :\ bc(a+b-c)\ :\ bc(a-b+c)\Big]
\ .
$$
(If we already know the point, we may instead just verify it is indeed the circumcenter of $\Delta GHI$.)
It also easily satisfies the equation of the line $AD$. Job done.
$\square$
