Molecular transport and interaction are of fundamental importance in biology and medicine. The spatiotemporal diffusion map can reflect the regulation of molecular interactions and their intracellular functions. To construct subcellular diffusion maps based on bio-imaging data, we explore a general optimization framework with diffusion equation constraints (OPT-PDE). For the solution of the spatially piecewise constant and anisotropic diffusion tensors, we develop an efficient block solver based on applying Newton's method to the first-order necessary condition for optimality. We characterize the wellposedness of the OPT-PDE model problem and the convergence properties of the solver. We also demonstrate the general utility of the solver in recovering spatially heterogeneous and anisotropic diffusion maps with computer-simulated bio-images. The results indicate that the solver can accurately recover piecewise-constant isotropic and anisotropic diffusion coefficients, while exhibiting efficient convergence and robustness. This work highlights the power of the OPT-PDE model and the solver in recovering the diffusion map from imaging data, and demonstrates that it has significant implications in bio-imaging analysis.