In computer aided diagnosis, a high-accuracy image indicating a temporal difference of a diseased part of a patient is produced from a current image and a past image as follows. First, a plurality of matching ROIs (Regions of Interest) are defined on the current image and the past image. A matching processing is performed on corresponding ROIs of the two images to determine points on the past image corresponding to the respective points on the current image. Obtained shift vectors are fitted to a two-dimensional polynomial. In the fitting, when a matching ROI of the current image is included in or overlapped with a previously-defined region of interest, shift vectors are weighted by a small factor. The current image is then modified according to obtained fitting coefficients. Finally, an image indicating the difference between the modified current image and the past image is produced.