I've run up this alley a couple of times. Let's look at what you've actually accomplished.
1) You picked an interval/DFT frame
2) You built basis vectors with sinusoids the length of the frame, including only the time points where you have samples.
As you observed, this is no longer an orthonormal basis. This means you are going to have to solve a matrix inverse to get your actual coordinates. With a regular DFT, the matrix is the identity matrix so solving the inverse is implied. The coordinates you solve for are very much like bin values at this point.
3) Your basis vectors are actually sampled points from continuous sinusoidal functions. Therefore you can construct a continuous function which includes all your original sample points.
What have you actually accomplished? Well, a very good interpolation of your data.
If it so happens that your signal is composed of tones that have a whole number of cycles, your coordinates/bin values will be spot on. However, if the tone is not a whole number of cycles, the leakage of your uneven sampling will not match the leakage of evenly sampled points, just as leakage doesn't match for different sized DFTs.
So, at this point, if you want to do any serious spectral analysis of your signal using DFT techniques, I agree with the recommendations you have gotten, resample your reconstructed signal at evenly spaced points and work with that.
You will get a different interpolation function for every different interval length. You need to be aware of the "wraparound" nature of the DFT. This means the interpolation function will repeat every interval length. Therefore you will want to make your interval wider than just the width of your samples. When you resample you will likely want to narrow the interval to just around your points. The extrapolations are useless.
Of course, there are much simpler interpolation techniques.