N images in a respiratory cycle which are radiographed in year P are input, and binary lung field images are extracted from the respective front chest images. Lung field areas S and lung field heights are then calculated. In forming area and height variation waveforms, regions of the N input images are obtained and plotted. Each image is determined as an image belonging to the inspiration mode or expiration mode. The respective images are sorted and stored. Similar processing is performed for N images in a respiratory cycle which are radiographed in year P+1, and the resultant images are stored. Difference images are obtained from the basic images radiographed in year P+1 and the reference images radiographed in year P for each mode by image analysis, thereby extracting changes over time.