I have a box in 3D space,I can move and rotate it freely and I have access to its homogeneous transformation matrix and form it I know the position and orientation.
The thing I want to do is to rotate the cube so only one of its axis align with the corresponding coordinates axis and with the same direction (let's say the X axis in red )but I want to keep the orientation on the other 2 axis
I tried setting the values in the matrix manually but that didn't work, and I tried rotate it on the Y, followed by Z, followed by Y and didn't work either.
The uploaded images show a simple situation where the first one is the box in an arbitrary position with arbitrary rotation on X, Y, Z and the second image shows the box where its X axis aligned with the X axis of the coordinates system (it was tilted a little in the first image and it's still tilted)
(There is another situation where the axis of the box is not on the same plane as the target axis). Any hints ?

Calling the “before” unit vectors $(\mathbf r,\mathbf g,\mathbf b)$ and the “after” vectors $(\mathbf r',\mathbf g',\mathbf b')$, you’re looking for a rotation that will map $\mathbf r$ to $\mathbf r'$. It’s not clear to me just what other constraints you have in mind for this rotation, so I’ll give you a starting point for constructing it.
The axis of this rotation must lie in the plane that is the angle bisector of $\mathbf r$ and $\mathbf r'$. A normal to this plane is $\mathbf r'-\mathbf r$, so any vector that’s perpendicular to this can serve as the rotation axis. The rotation angle $\theta$ will range from a maximum of $\pi$ when rotating about $\mathbf r+\mathbf r'$ to some minimum obtained by rotating about $\mathbf r\times\mathbf r'$.
Once you’ve chosen a unit vector $\mathbf k$ for the rotation axis, assembling the corresponding rotation matrix is straightforward if you remember that the columns of a transformation matrix are the images of the basis vectors and that the inverse of an orthogonal matrix is its transpose. Thus, if you have two orthonormal frames with respective unit vectors $(\mathbf u, \mathbf v, \mathbf w)$ and $(\mathbf u', \mathbf v', \mathbf w')$, the orthogonal transformation that maps the first onto the second is given by the matrix $$R = \begin{bmatrix}\mathbf u' & \mathbf v' & \mathbf w'\end{bmatrix} \begin{bmatrix} \mathbf u^T \\ \mathbf v^T \\ \mathbf w^T \end{bmatrix}.$$ (Here I’m following the usual mathematical convention of column vectors and left-multiplication by the transformation matrix.) If both of the frames have the same orientation (handedness), then this transformation is a rotation. For your rotation, we want the axis $\mathbf k$ to remain fixed, so set $\mathbf u = \mathbf u' = \mathbf k$. We want $\mathbf r$ to be mapped to $\mathbf r'$. The components of these vectors parallel to $\mathbf k$ remain fixed, so we can use the orthogonal rejections from $\mathbf k$ as the second orthogonal axes for the rotation: set $\mathbf v$ to $\mathbf r - (\mathbf k\cdot\mathbf r)\mathbf k = -\mathbf k\times(\mathbf k\times\mathbf r)$, normalized, and similarly for $\mathbf v'$ and $\mathbf r'$. Finally, for each frame we need a third unit vector orthogonal to the other two. A simple choice is $\mathbf w = \mathbf u\times\mathbf v$ and $\mathbf w'=\mathbf u'\times\mathbf v'$.
The amount by which $\mathbf b$ and $\mathbf g$ would be perturbed by this rotation can be computed via Rodrigues’ rotation formula. As noted above, when a vector $\mathbf v$ is rotated, only the component $\mathbf v_\perp$ orthogonal to the rotation axis $\mathbf k$ is changed. The net change is $$\mathbf v_{\perp\text{rot}}-\mathbf v_\perp = \cos\theta\,\mathbf v_\perp+\sin\theta\,\mathbf k\times\mathbf v - \mathbf v_\perp = (1-\cos\theta)\,\mathbf v_\perp+\sin\theta\,\mathbf k\times\mathbf v.$$ The orthogonal component $\mathbf v_\perp$ is the orthogonal rejection from $\mathbf k$, i.e., $$\mathbf v_\perp = \mathbf v-(\mathbf k\cdot\mathbf v)\mathbf k = -\mathbf k\times(\mathbf k\times\mathbf v).$$ The necessary rotation angle $\theta$ can be expressed as a function of the chosen axis $\mathbf k$, and $\mathbf k$ itself can be parameterized as ${\mathbf r\times\mathbf r' \over \|\mathbf r\times\mathbf r'\|}\cos t+{\mathbf r'-\mathbf r \over \|\mathbf r'-\mathbf r\|}\sin t$, so in principle you could construct a single-variable cost function that represents your constraints on the motion of the other two axes and then optimize based on that, but I suspect that the minimum-angle rotation with $\mathbf k$ parallel to $\mathbf r\times\mathbf r'$ will serve.