In this paper, we consider the problem of recovering the phase information of the multiple images from the multiple mixed phaseless Short-Time Fourier Transform (STFT) image measurements, which is called the blind multiple input multiple output image phase retrieval (BMIPR) problem. It is an inherently ill-posed problem due to the lack of the phase and mixing information, and the existing phase retrieval algorithms are not explicitly designed for this case. To address the BMIPR phase retrieval problem, an integrated algorithm is presented, which combines a gradient descent (GD) algorithm by minimizing a nonconvex loss function with an independent component analysis (ICA) algorithm and a non-local means (NM) algorithm. Experimental evaluation has been conducted to show that under appropriate conditions the proposed algorithms can explicitly recover the images, the phases of the images and the mixing matrix. In addition, the algorithm is robust to noise. Index Terms-Blind multiple input multiple output image phase retrieval (BMIPR), short-time Fourier transform (STFT), non-convex optimization, independent component analysis (ICA), non-local means (NM).