As the technology scales down, two-dimensional (2D) NAND flash memory has reached its bottleneck. Three-dimensional (3D) NAND flash memory was proposed to further increase the storage capacity by vertically stacking multiple layers. However, the new architecture of 3D flash memory leads to new sources of errors, which severely affects the reliability of the system. In this paper, for the first time, we derive the channel probability density function of 3D NAND flash memory by taking major sources of errors. Based on the derived channel probability density function, the mutual information (MI) for 3D flash memory with multiple layers is derived and used as a metric to design the quantization. Specifically, we propose a dynamic programming algorithm to jointly optimize read-voltage thresholds for all layers by maximizing the MI (MMI). To further reduce the complexity, we develop an MI derivative (MID)-based method to obtain read-voltage thresholds for hard-decision decoding (HDD) of error correction codes (ECCs). Simulation results show that the performance with jointly optimized read-voltage thresholds can closely approach that with read-voltage thresholds optimized for each layer, with much less read latency. Moreover, the MID-based MMI quantizer almost achieves the optimal performance for HDD of ECCs.