A systematic way to constructing optimized interpolating operators for two-hadron systems is developed by incorporating inter-hadron spatial wavefunctions. The wavefunctions can be obtained from an iterative process with an appropriate initial guess. To implement these operators, a novel quark smearing technique utilizing Z_3 noise vectors is proposed, which allows for effectively incorporating inter-hadron spatial wavefunctions at the source without using all-to-all quark propagators. Proof-of-principle application to the Ω_cccΩ_ccc system using physical-point lattice configurations with a large size La≃8.1 fm shows that optimized operators enables clear identification of states around 2m_Ω_ccc≃ 9700 MeV with the energy gap as narrow as ∼ 5 MeV. A comparison on correlation functions, effective energies, and HAL QCD potentials between unoptimized operators and optimized operators is given, with a special emphasis on the effects from nearby elastic scattering states. Potential applicability of the optimized operator to various two-hadron systems and its relation to the variational method are also discussed.