This paper describes a computational method using tensor math for higher order singular value decomposition (HOSVD) of registered images. Tensor decomposition is a rigorous way to expose structure embedded in multidimensional datasets. Given a dataset of registered 2-D images, the dataset is represented in tensor format and HOSVD of the tensor is computed to obtain a set of 2-D basis images. The basis images constitute a linear decomposition of the original dataset. HOSVD is data-driven and does not require the user to select parameters or assign thresholds. A specific application uses the basis images for pixel-level fusion of registered images into a single image for visualization. The fusion is optimized with respect to a measure of mean squared error. HOSVD and image fusion are illustrated empirically with four real datasets: (1) visible and infrared data of a natural scene, (2) MRI and x ray CT brain images, and in nondestructive testing (3) x ray, ultrasound, and eddy current images, and (4) x ray, ultrasound, and shearography images.