Let us first consider the case $D=I$. An invertible matrix $S$ can be (uniquely) factorized as $S=RQ$, where $R$ is upper triangular with positive diagonal entries and $Q$ is orthogonal, $QQ^\top=I$ (this is QR decomposition). Then $SS^\top = RQQ^\top R^\top = RR^\top$. On the other hand, by Cholesky decomposition, every positive definite symmetric matrix $A$ can be uniquely decomposed as $A=RR^\top$, where $R$ is as above. So in this case the space of such $S$ is $R\cdot O(n)$. This reflects the fact that the space of positive definite matrices is isomorphic (in a suitable sense) to $GL(n,\mathbb{R})/O(n)$ (through the mapping $X\mapsto XX^\top$).
What follows lacks some (many) details and should definitely not be taken for granted.
Now for a general non-singular $D$ we would like to use the same consideration, and for that we need $QDQ^\top=D$. This means that $Q$ is an element of the corresponding indefinite orthogonal group $O(D)$.
The analogue for QR decomposition in this case is called HR decomposition (sometimes also hyperbolic QR). Denote $\mathcal{X}$ the set of all matrices $X$ such that the leading principal minors of $X^\top D X$ have the same signs as the corresponding minors of $D$. Then every matrix $X\in\mathcal{X}$ has a decomposition $X=HR$, where $HDH^\top=D$ and $R$ is upper triangular with positive entries (similarly, one can take $R$ to be lower triangular).
If the leading principal minors of $A$ are non-zero, it posesses LDL decomposition $A=LCL^\top$, where $L$ is lower triangular with $1$ on the diagonal and $C$ is diagonal. This can be rewritten as $A=L'DL'^\top$, where $D$ is as above (up to a different order of the diagonal entries) and $L'$ is lower triangular with positive entries on the diagonal (namely, the square roots of the absolute values of the diagonal entries of $C$).
Again, if $D$ is ordered as the inertia matrix obtained from LDL decomposition, $A = SDS^\top = LDL^\top$, hence $D = L^{-1}SDS^\top L^{-\top}$ and $L^{-1}S \in O(D)$. If $D$ is ordered differently, the answer differs by a multiplication with a permutation matrix.
The relation with the HR decomposition is based on two observations: 1) the entries of $C$ can be chosen as $C = \operatorname{diag}(\Delta_1,\frac{\Delta_2}{\Delta_1},\ldots,\frac{\Delta_n}{\Delta_{n-1}})$, where $\Delta_i$ are the leading principal minors of $A$ (this is Jacobi method for the diagonalization of quadratic forms); 2) the upper-left corner submatrix of $LML^\top$ for any $M$ and lower-triangular $L$ is the product of the corresponding corner submatrices, and so in case the diagonal entries of $L$ are positive, the leading principal minors of $LML^\top$ and $M$ have the same signs.
I am not sure what to do in case $D$ is singular. Apparently, one can devise some sort of combination of HR and thin QR decompositions, but I do not see this in the literature.