WO2018193517A1 - 画像処理装置 - Google Patents

画像処理装置 Download PDF

Info

Publication number
WO2018193517A1
WO2018193517A1 PCT/JP2017/015585 JP2017015585W WO2018193517A1 WO 2018193517 A1 WO2018193517 A1 WO 2018193517A1 JP 2017015585 W JP2017015585 W JP 2017015585W WO 2018193517 A1 WO2018193517 A1 WO 2018193517A1
Authority
WO
WIPO (PCT)
Prior art keywords
image
correction value
unit
radiance
calculation unit
Prior art date
Application number
PCT/JP2017/015585
Other languages
English (en)
French (fr)
Inventor
真梨子 佐藤
将敬 鈴木
玉川　恭久
Original Assignee
三菱電機株式会社
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by 三菱電機株式会社 filed Critical 三菱電機株式会社
Priority to PCT/JP2017/015585 priority Critical patent/WO2018193517A1/ja
Priority to JP2019512336A priority patent/JP6556409B2/ja
Publication of WO2018193517A1 publication Critical patent/WO2018193517A1/ja

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T1/00—General purpose image data processing
    • H—ELECTRICITY
    • H04—ELECTRIC COMMUNICATION TECHNIQUE
    • H04N—PICTORIAL COMMUNICATION, e.g. TELEVISION
    • H04N1/00—Scanning, transmission or reproduction of documents or the like, e.g. facsimile transmission; Details thereof
    • H04N1/40—Picture signal circuits
    • H04N1/407—Control or modification of tonal gradation or of extreme levels, e.g. background level

Definitions

  • the present invention relates to an image processing apparatus that corrects spectral radiance in a panchromatic image.
  • the contrast of a part of the image acquired from the sensor may be lowered due to the influence of soot or fog contained in the atmosphere. Since a decrease in contrast causes a decrease in the visibility of an image, an image processing apparatus that corrects a decrease in contrast due to the influence of haze or fog has been developed.
  • Patent Document 1 a difference between a signal level of an arbitrary pixel among a plurality of pixels constituting an image and a signal level of a pixel existing in an area around the pixel is measured.
  • An image processing apparatus that raises the signal level of the pixel is disclosed by assuming that the contrast is lowered if is large.
  • any pixel corresponds to the sea area and is present in the area around the pixel May correspond to the coastal area.
  • the reflectance of the sea and the reflectance of the coast there is a big difference between the reflectance of the sea and the reflectance of the coast, so even if it is not affected by the haze or fog contained in the atmosphere, the signal level of any pixel and the surroundings of the pixel
  • the difference from the signal level of the pixels existing in the region may increase. In such a case, assuming that the contrast is lowered and increasing the signal level of the pixel, the visibility of the image is lowered.
  • the image acquired from the sensor is an image obtained by photographing the coastline
  • the image includes a plurality of regions with greatly different reflectances, such as a boundary between a mountainous area and an urban area
  • the present invention has been made to solve the above-described problems. Even when the panchromatic image acquired from the sensor is an image including a plurality of regions having greatly different reflectances, it is possible to reduce the fog or fog.
  • An object of the present invention is to obtain an image processing apparatus capable of correcting a contrast that has been affected and decreased.
  • the image processing apparatus includes an image dividing unit that divides a multispectral image including images of a plurality of wavelength bands into a plurality of regions, and a plurality of wavelength bands of the plurality of images included in the multispectral image.
  • a correction value calculation unit that calculates a correction value of spectral radiance in a panchromatic image that is an image in a single wavelength band including a plurality of wavelength bands, and the correction unit calculates the spectral radiation calculated by the correction value calculation unit
  • the spectral radiance in the panchromatic image is corrected using the luminance correction value.
  • the spectral radiance in a panchromatic image which is an image in a single wavelength band including a plurality of wavelength bands, is calculated from the atmospheric scattered radiance in each region calculated by the scattering luminance calculation unit for each wavelength band. Since the correction value calculation unit for calculating the correction value is provided, and the correction unit is configured to correct the spectral radiance in the panchromatic image using the correction value of the spectral radiance calculated by the correction value calculation unit, Even when the panchromatic image acquired from the sensor is an image including a plurality of regions having greatly different reflectances, there is an effect that the contrast that is lowered due to the influence of haze or fog can be corrected. .
  • FIG. 1 is a configuration diagram illustrating an image processing apparatus according to Embodiment 1 of the present invention; It is a hardware block diagram which shows the image processing apparatus by Embodiment 1 of this invention.
  • FIG. 2 is a hardware configuration diagram of a computer when the image processing apparatus is realized by software or firmware. It is a flowchart which shows the process sequence in case an image processing apparatus is implement
  • FIG. It is explanatory drawing which shows the atmospheric scattering radiance Lscat (b, i, j) modeled by the scattered luminance calculation process part 5.
  • FIG. It is explanatory drawing which shows the atmospheric transmittance
  • FIG. 2 It is a hardware block diagram which shows the image processing apparatus by Embodiment 2 of this invention. It is explanatory drawing which shows the spectral radiance L sensor (b, i, j) modeled by the reflected luminance calculation part 41.
  • FIG. It is a block diagram which shows the image processing apparatus by Embodiment 3 of this invention. It is a hardware block diagram which shows the image processing apparatus by Embodiment 3 of this invention.
  • FIG. FIG. 1 is a block diagram showing an image processing apparatus according to Embodiment 1 of the present invention
  • FIG. 2 is a hardware block diagram showing the image processing apparatus according to Embodiment 1 of the present invention
  • the image dividing unit 1 is realized by, for example, an image dividing circuit 21 shown in FIG.
  • the image dividing unit 1 acquires a multispectral image including images of a plurality of wavelength bands from the multispectral sensor, and performs a process of dividing the multispectral image into a plurality of small regions (regions).
  • the parameter storage unit 2 is realized by, for example, the parameter storage circuit 22 shown in FIG.
  • the parameter storage unit 2 stores scattering characteristic data indicating scattering characteristics in each wavelength band, scattered transmission characteristic data that is a function indicating atmospheric transmittance in a small area, and wavelength characteristic data indicating a wavelength band that can be observed by the sensor. ing.
  • the scattered luminance calculation unit 3 is realized by, for example, the scattered luminance calculation circuit 23 illustrated in FIG. 2, and includes a dark portion luminance calculation unit 4 and a scattered luminance calculation processing unit 5.
  • the scattered luminance calculation unit 3 uses the scattering characteristic data stored in the parameter storage unit 2, and uses the scattering characteristic data stored in the parameter storage unit 2 for each wavelength band of a plurality of images included in the multispectral image. Processing for calculating atmospheric scattering radiance indicating atmospheric scattering in the region is performed.
  • the dark part luminance calculating unit 4 of the scattered luminance calculating unit 3 calculates a histogram of each small region divided by the image dividing unit 1 for each wavelength band of a plurality of images included in the multispectral image, and the histogram Thus, the process of calculating the spectral radiance of the dark area where the spectral radiance is relatively small in each of the small areas divided by the image dividing unit 1 is performed.
  • the scattered luminance calculation processing unit 5 of the scattered luminance calculation unit 3 uses the scattering characteristic data stored in the parameter storage unit 2 and the spectral radiance of the dark area calculated by the dark portion luminance calculation unit 4 to generate a multispectral image. For each wavelength band of a plurality of included images, a process of calculating atmospheric scattering radiance in each small region divided by the image dividing unit 1 is performed.
  • the atmospheric transmittance calculating unit 6 is realized by, for example, the transmittance calculating circuit 24 shown in FIG.
  • the atmospheric transmittance calculation unit 6 uses the scattered transmission characteristic data stored in the parameter storage unit 2 to calculate from the atmospheric scattering radiance in each small region calculated for each wavelength band by the scattering luminance calculation processing unit 5. For each wavelength band, processing for calculating the atmospheric transmittance of atmospheric scattering in each small region divided by the image dividing unit 1 is performed.
  • the correction value calculation unit 7 is realized by, for example, the correction value calculation circuit 25 illustrated in FIG. 2, and includes a luminance correction value calculation unit 8, a transmittance correction value calculation unit 9, an image recombination unit 10, and a correction value interpolation unit. 11 is provided.
  • the correction value calculation unit 7 uses the atmospheric scattering radiance in each small region calculated for each wavelength band by the scattering luminance calculation processing unit 5 in a panchromatic image that is an image in a single wavelength band including a plurality of wavelength bands. Processing for calculating a correction value of the spectral radiance is performed.
  • the correction value calculation unit 7 calculates the correction value of the spectral radiance in the panchromatic image from the atmospheric transmittance of atmospheric scattering in each small region calculated for each wavelength band by the atmospheric transmittance calculation unit 6.
  • the panchromatic image is an image in a single wavelength band including a plurality of image wavelength bands included in the multispectral image.
  • the luminance correction value calculation unit 8 of the correction value calculation unit 7 uses the wavelength characteristic data stored in the parameter storage unit 2 and the atmospheric air in each small region calculated for each wavelength band by the scattering luminance calculation processing unit 5. A process of calculating the first correction value of the spectral radiance in each small region from the scattered radiance is performed.
  • the transmittance correction value calculation unit 9 of the correction value calculation unit 7 uses the wavelength characteristic data stored in the parameter storage unit 2 and uses the wavelength characteristic data stored in the parameter storage unit 2 in each small region calculated for each wavelength band.
  • a process of calculating a second correction value of the spectral radiance in each small region from the atmospheric transmittance of atmospheric scattering is performed.
  • the image recombination unit 10 of the correction value calculation unit 7 combines the first correction values of the spectral radiances in the respective small areas calculated by the luminance correction value calculation unit 8 to obtain a two-dimensional first correction value.
  • a process for generating the array A is performed.
  • the image recombination unit 10 combines the second correction value of the spectral radiance in each small region calculated by the transmittance correction value calculation unit 9 to generate a two-dimensional array B of the second correction values. Perform the process to generate.
  • Each element in the two-dimensional arrays A and B corresponds to each small region divided by the image dividing unit 1.
  • Each element in the two-dimensional array A has a first correction value of the spectral radiance in each small area divided by the image dividing unit 1, and each element in the two-dimensional array B is an image division
  • the second correction value of the spectral radiance in each of the small areas divided by the unit 1 is provided.
  • the correction value interpolation unit 11 of the correction value calculation unit 7 is generated by the image recombination unit 10 so that the resolution of the two-dimensional arrays A and B generated by the image recombination unit 10 matches the resolution of the panchromatic image.
  • a process for interpolating the two-dimensional arrays A and B is performed.
  • the correction unit 12 is realized by the correction circuit 26 shown in FIG. 2 and includes a radiance correction unit 13 and a transmittance correction unit 14.
  • the radiance correction unit 13 of the correction unit 12 performs a process of acquiring a panchromatic image from the panchromatic sensor.
  • the radiance correction unit 13 performs a process of correcting the spectral radiance in the panchromatic image using the first correction value of each element in the two-dimensional array A interpolated by the correction value interpolation unit 11.
  • the transmittance correction unit 14 performs a process of correcting the spectral radiance in the panchromatic image by using the second correction value of each element in the two-dimensional array B interpolated by the correction value interpolation unit 11.
  • This is assumed to be realized by dedicated hardware as shown in FIG. That is, it is assumed to be realized by the image dividing circuit 21, the parameter storage circuit 22, the scattering luminance calculation circuit 23, the transmittance calculation circuit 24, the correction value calculation circuit 25, and the correction circuit 26.
  • the parameter storage circuit 22 includes, for example, a RAM (Random Access Memory), a ROM (Read Only Memory), a flash memory, an EPROM (Erasable Programmable Read Only Memory Memory), an EEPROM (Electrically Erasable Memory), and the like.
  • a volatile semiconductor memory a magnetic disk, a flexible disk, an optical disk, a compact disk, a mini disk, a DVD (Digital Versatile Disc), or the like is applicable.
  • the image dividing circuit 21, the scattered luminance calculation circuit 23, the transmittance calculation circuit 24, the correction value calculation circuit 25, and the correction circuit 26 are, for example, a single circuit, a composite circuit, a programmed processor, a processor programmed in parallel, An application specific integrated circuit (ASIC), a field-programmable gate array (FPGA), or a combination thereof corresponds to the ASIC (Application Specific Integrated Circuit).
  • ASIC Application Specific integrated circuit
  • FPGA field-programmable gate array
  • the components of the image processing apparatus are not limited to those realized by dedicated hardware, and the image processing apparatus may be realized by software, firmware, or a combination of software and firmware.
  • Software or firmware is stored as a program in the memory of a computer.
  • the computer means hardware that executes a program, and includes, for example, a CPU (Central Processing Unit), a central processing unit, a processing unit, an arithmetic unit, a microprocessor, a microcomputer, a processor, a DSP (Digital Signal Processor), and the like. .
  • FIG. 3 is a hardware configuration diagram of a computer when the image processing apparatus is realized by software or firmware.
  • the parameter storage unit 2 is configured on the memory 31 of the computer, the image dividing unit 1, the scattered luminance calculation unit 3, the atmospheric transmittance calculation unit 6, and the correction value calculation.
  • a program for causing the computer to execute the processing procedures of the unit 7 and the correction unit 12 may be stored in the memory 31, and the processor 32 of the computer may execute the program stored in the memory 31.
  • FIG. 4 is a flowchart showing a processing procedure when the image processing apparatus is realized by software or firmware.
  • FIG. 2 shows an example in which each component of the image processing apparatus is realized by dedicated hardware
  • FIG. 3 shows an example in which the image processing apparatus is realized by software, firmware, etc. Some components in the image processing apparatus may be realized by dedicated hardware, and the remaining components may be realized by software, firmware, or the like.
  • the image dividing unit 1 acquires a multispectral image including images of a plurality of wavelength bands from the multispectral sensor, and the radiance correction unit 13 acquires a panchromatic image from the panchromatic sensor.
  • each of the wavelength bands of the plurality of images included in the multispectral image has a narrow band
  • the wavelength band of the panchromatic image includes the wavelength bands of the plurality of images included in the multispectral image.
  • the wavelength band of a plurality of images included in the multispectral image is referred to as a band b.
  • the multispectral image is an image of the band b1, an image of the band b2, and an image of the band b3. It is assumed that the image and the image of band b4 are included.
  • the image dividing unit 1 acquires a multispectral image from the multispectral sensor (step ST1 in FIG. 4), and divides the multispectral image into N small regions S (i, j) as shown in FIG. Step ST2 in FIG.
  • FIG. 5 is an explanatory diagram showing a multispectral image divided by the image dividing unit 1.
  • the division number N of the multispectral image can be arbitrarily set, and is determined from, for example, the size of the subject included in the multispectral image.
  • the number of divisions N of the multispectral image is set so that the subject is included in a small area that includes a subject included in the multispectral image and a dark area such as a shadow of the subject. To do.
  • the multispectral image is divided based on the size of the largest subject among the multiple subjects included in the multispectral image.
  • the number N may be set to reduce the number of small areas S (i, j) that do not include the subject.
  • FIG. 6 is an explanatory diagram showing a multispectral image divided by the image dividing unit 1.
  • the dark region in the dark portion luminance calculation unit 4 is more than in the case where the multispectral image is divided as shown in FIG.
  • the number of calculated spectral radiances can be reduced.
  • the small region S (i, j) indicates the location where the change rate of the density changes.
  • the multispectral image may be divided so as to include it.
  • the dark part luminance calculating unit 4 of the scattered luminance calculating unit 3 displays a histogram of each small region S (i, j) divided by the image dividing unit 1 for each band b of a plurality of images included in the multispectral image. Each is calculated (step ST3 in FIG. 4).
  • FIG. 7 is an explanatory diagram showing a histogram of the small area S (i, j) in the band b calculated by the dark part luminance calculation unit 4.
  • the dark portion luminance calculation unit 4 relatively compares the spectral radiance in each of the small regions S (i, j) divided by the image dividing unit 1 for each of the bands b of the plurality of images included in the multispectral image.
  • FIG. 7 shows an example in which the minimum value of the spectral radiance is the spectral radiance L D (b, i, j) of the dark region in the histogram of the small region S (i, j).
  • the dark portion luminance calculation unit 4 calculates a histogram of the small region S (i, j).
  • the dark portion luminance calculation unit 4 has a pixel value among a plurality of pixels constituting the multispectral image. There may be pixels that are not included.
  • the dark part luminance calculation unit 4 calculates a histogram of the small region S (i, j) by excluding pixels having no pixel value from a plurality of pixels constituting the multispectral image. You may do it.
  • the dark part luminance calculation unit 4 configures the multispectral image when a plurality of pixels constituting the multispectral image include defective pixels whose pixel values are always near zero.
  • the histogram of the small region S (i, j) may be calculated by excluding defective pixels from a plurality of pixels.
  • the scattered luminance calculation processing unit 5 of the scattered luminance calculation unit 3 acquires the scattering characteristic data ⁇ (b) stored in the parameter storage unit 2 (step ST5 in FIG. 4).
  • the scattering characteristic data ⁇ (b) stored in the parameter storage unit 2 is data indicating the scattering characteristic of the band b and is a value unique to the band.
  • the scattering characteristic data ⁇ (b) can be calculated from the radiance ratio of atmospheric scattering at each wavelength using, for example, a light-wave atmospheric propagation calculation simulator such as MODTRAN (Moderate resolution atomic TRANmission).
  • the scattered luminance calculation processing unit 5 performs atmospheric scattering radiance L scat (in each small region S (i, j) divided by the image dividing unit 1 for each band b of the plurality of images included in the multispectral image.
  • b, i, j) is a model of the product of the scattering characteristic data ⁇ (b) and the eigenvalue L scat0 (i, j) of the small region S (i, j) as shown in the following equation (1). Turn into.
  • L scat (b, i, j) ⁇ (b) ⁇ L scat0 (i, j) (1)
  • the atmospheric scattering radiance L scat (b, i, j) is a spectral radiance indicating atmospheric scattering in band b.
  • the eigenvalue L scat0 (i, j) of the small area S (i, j) has no wavelength dependency and is an eigenvalue of the small area S (i, j) depending on the imaging environment.
  • FIG. 8 is an explanatory diagram showing the atmospheric scattered radiance L scat (b, i, j) modeled by the scattered luminance calculation processing unit 5.
  • the solid line is a model of atmospheric scattering radiance L scat (b, i, j), and ⁇ is the spectral radiance of the dark area in the small area S (i, j).
  • the scattered luminance calculation processing unit 5 calculates the spectral radiance L D (b, i, j) of the dark region and the right side of the expression (1) for each of the bands b of the plurality of images included in the multispectral image.
  • the eigenvalue L scat0 (i, j) that is the closest to ⁇ (b) ⁇ L scat0 (i, j) is estimated by regression analysis.
  • the number of measurement points of the spectral radiance L D (b, i, j) in the dark region is the same as the number of bands that the multispectral sensor has.
  • a typical example of a parameter determination method based on regression analysis is a least square method.
  • the scattering intensity calculation processing unit 5 the estimated eigenvalues L scat0 (i, j), the estimated eigenvalues L Scat0 the (i, j) and scattering characteristic data alpha (b) into equation (1), band b Every time, atmospheric scattering radiance L scat (b, i, j) in each small region S (i, j) is calculated (step ST6 in FIG. 4).
  • the scattered luminance calculation processing unit 5 outputs the calculated atmospheric scattered radiance L scat (b, i, j) in each small region S (i, j) to the atmospheric transmittance calculating unit 6 and the luminance correction value calculating unit 8. To do.
  • the scattered luminance calculation processing unit 5 and the measurement point of the spectral radiance L D (b, i, j) in the dark area in a certain small area S (i, j) and the atmospheric scattered radiance L scat (b, i , J) is low in correlation with the model, the atmospheric scattering radiance L scat (b, i, j) in the small region S (i, j) is regarded as an abnormal value, for example, a negative value is output.
  • the atmospheric scattering radiance L scat (b, i) in the small region S (i, j) including the dark region where the measurement point of the spectral radiance L D (b, i, j) is not similar to the model. , J) can be removed.
  • the correlation between the measurement point and the model for example, a case where an R square value, which is a determination coefficient of regression analysis in the least square method, is equal to or less than a threshold value Rth can be cited.
  • the atmospheric transmittance calculation unit 6 acquires the scattered transmission characteristic data f (a) stored in the parameter storage unit 2 (step ST7 in FIG. 4).
  • a (b, L scat (b, i, j)).
  • Scattered transmitted are stored by the parameter storage unit 2
  • characteristic data f (a) is a small area S (i, j) Air scattered radiance L scat showing the atmospheric scattering in (b, i, j) air permeability of ⁇ It is a function for each band b indicating (b, i, j).
  • the scattered transmission characteristic data f (a) can be calculated by a light wave atmospheric propagation calculation simulator such as MODTRAN, for example.
  • FIG. 9 is an explanatory diagram showing the atmospheric transmittance ⁇ (b, i, j) corresponding to the atmospheric scattered radiance L scat (b, i, j) of the small region S (i, j) in the band b.
  • the atmospheric transmittance calculation unit 6 uses the band b and the small region S (output from the scattering luminance calculation processing unit 5) in the variable a of the scattered transmission characteristic data f (a).
  • the atmospheric transmittance ⁇ (b, i) for atmospheric scattering in the small region S (i, j) for each band b. , J) is calculated (step ST8 in FIG. 4).
  • the atmospheric transmittance calculating unit 6 outputs the calculated atmospheric transmittance ⁇ (b, i, j) of atmospheric scattering in each small region S (i, j) to the transmittance correction value calculating unit 9. However, if the atmospheric transmittance ⁇ (b, i, j) in a certain small region S (i, j) is separated from the assumed value by several percent or more, the atmospheric transmittance calculating unit 6 For example, a negative value is output as an abnormal value of the atmospheric transmittance ⁇ (b, i, j) in S (i, j).
  • the luminance correction value calculation unit 8 of the correction value calculation unit 7 acquires the wavelength characteristic data stored in the parameter storage unit 2 (step ST9 in FIG. 4).
  • the wavelength characteristic data stored in the parameter storage unit 2 is data indicating each of the observable wavelength band of the multispectral sensor and the observable wavelength band of the panchromatic sensor.
  • the brightness correction value calculation unit 8 uses the acquired wavelength characteristic data, and the atmospheric scattering radiance L scat (b in each small region S (i, j) for each band b output from the scattering brightness calculation processing unit 5. , I, j) is calculated. Details of the processing for calculating the weighting coefficient ⁇ (b) will be described later.
  • the luminance correction value calculation unit 8 sets the first correction value L scat (Pa ′, i, j) of the spectral radiance in each small region S (i, j).
  • the sum of the atmospheric scattering radiance L scat (b, i, j) in each small region S (i, j) for each band b is calculated using the weight coefficient ⁇ (b) (step ST10 in FIG. 4). ).
  • b b1, b2, b3, b4 Pa ′ is a single wavelength band of an image corresponding to the panchromatic image observed by the panchromatic sensor.
  • the luminance correction value calculation unit 8 outputs the first correction value L scat (Pa ′, i, j) of the spectral radiance in each small region S (i, j) to the image recombination unit 10.
  • the luminance correction value calculation unit 8 is abnormal when the atmospheric scattering radiance L scat (b, i, j) in the small region S (i, j) output from the scattering luminance calculation processing unit 5 is an abnormal value.
  • the first correction value L scat (Pa ′, i, j) of the spectral radiance in the small region S (i, j) is used without using the atmospheric scattering radiance L scat (b, i, j) that is a value. Try to calculate.
  • the brightness correction value calculation unit 8 uses the atmospheric scattering radiance in the small area around the small area S (i, j) where the atmospheric scattered radiance L scat (b, i, j) is an abnormal value, The atmospheric scattering radiance L scat (b, i, j) in the small region S (i, j) is interpolated. Then, the luminance correction value calculation unit 8 calculates the first correction value L scat (Pa ′, i, j) using the interpolated atmospheric scattering radiance L scat (b, i, j). .
  • FIG. 10 is an explanatory diagram showing the observable wavelength band of the multispectral sensor and the observable wavelength band of the panchromatic sensor.
  • the observable wavelength band of the panchromatic sensor is ⁇ 0 (Pa) to ⁇ 1 (Pa).
  • the observable wavelength bands of the multispectral sensor are ⁇ 0 (b1) to ⁇ 1 (b1), ⁇ 0 (b2) to ⁇ 1 (b2), ⁇ 0 (b3) to ⁇ 1 (b3), and ⁇ 0 ( b4) to ⁇ 1 (b4).
  • a wavelength band ⁇ '0 (b) ⁇ ⁇ ' 1 (b) Is performed so as to correspond to each of the bands b1, b2, b3, and b4 of the image included in the multispectral image. That is, the luminance correction value calculation unit 8 assigns the wavelength band ⁇ ′ 0 (b1) to ⁇ ′ 1 (b1) to the band b1, and the wavelength band ⁇ ′ 0 (b2) to ⁇ ′ 1 to the band b2.
  • wavelength bands ⁇ ′ 0 (b3) to ⁇ ′ 1 (b3) are assigned to the band b3.
  • the luminance correction value calculation unit 8 assigns wavelength bands ⁇ ′ 0 (b4) to ⁇ ′ 1 (b4) to the band b4.
  • the wavelength bands ⁇ ′ 0 (b1) and ⁇ ′ 1 (b1) are expressed as in the following equations (5) and (6).
  • the numerator in the weighting coefficient ⁇ (b) simulates wavelength integration with respect to the spectral radiance in the assigned wavelength band ⁇ ′ 0 (b) to ⁇ ′ 1 (b). Further, the denominator in the weighting coefficient ⁇ (b) means that the radiance of the panchromatic sensor is converted into the spectral radiance.
  • the above equation (3) assumes that the spectral radiance in the assigned wavelength bands ⁇ ′ 0 (b) to ⁇ ′ 1 (b) is constant, and is assigned to the wavelength band of the panchromatic sensor. As the number of bands increases, the first correction value L scat (Pa ′, i, j) of the spectral radiance can be calculated with higher accuracy.
  • the luminance correction value calculation unit 8 calculates the weighting coefficient ⁇ (b) by Expression (4).
  • each wavelength in the panchromatic sensor is calculated.
  • the weighting coefficient ⁇ (b) may be calculated based on the spectral sensitivity R ( ⁇ ) of ⁇ .
  • the first correction value L scat (Pa ′, i, j) of the spectral radiance is increased. It can be calculated with accuracy.
  • the luminance correction value calculation unit 8 may calculate the weighting coefficient ⁇ (b) using a panchromatic image of a real image. For example, the luminance correction value calculation unit 8 acquires the spectral radiance L (Pa, x, y) of the panchromatic image of the actual image. (X, y) are coordinates on the panchromatic image, and are different from (i, j) indicating the arrangement of the small regions in the multispectral image. Then, the luminance correction value calculation unit 8 simulates the spectral radiance first correction value L scat (Pa ′, i, j) obtained from Equation (3) from the multispectral image.
  • the transmittance correction value calculation unit 9 of the correction value calculation unit 7 acquires the wavelength characteristic data stored in the parameter storage unit 2 (step ST11 in FIG. 4).
  • the transmittance correction value calculation unit 9 uses the acquired wavelength characteristic data, and the atmospheric transmittance ⁇ of atmospheric scattering in each small region S (i, j) for each band b output from the atmospheric transmittance calculation unit 6.
  • a weight coefficient ⁇ ′ (b) of (b, i, j) is calculated. Details of the processing for calculating the weighting coefficient ⁇ ′ (b) will be described later.
  • the transmittance correction value calculation unit 9 has an abnormal value for the atmospheric transmittance ⁇ (b, i, j) of atmospheric scattering in the small region S (i, j) output from the atmospheric transmittance calculation unit 6.
  • the second correction value ⁇ (Pa ′, i, j) of the spectral radiance in the small region S (i, j) is calculated without using the atmospheric transmittance ⁇ (b, i, j) that is an abnormal value.
  • the transmittance correction value calculation unit 9 uses the atmospheric transmittance in the small region around the small region S (i, j) where the atmospheric transmittance ⁇ (b, i, j) is an abnormal value, to use the small region.
  • the transmittance correction value calculation unit 9 calculates the second correction value ⁇ (Pa ′, i, j) using the interpolated atmospheric transmittance ⁇ (b, i, j).
  • the weighting coefficient ⁇ ′ (b) of the atmospheric transmittance ⁇ (b, i, j) can be calculated by the equation (4) similarly to the luminance correction value calculation unit 8.
  • the transmittance correction value calculation unit 9 uses the spectral sensitivity R ( ⁇ ) of each wavelength ⁇ and the illuminance E ( ⁇ ) of the light source in the panchromatic sensor, as shown in the following formula (9), to give a weighting coefficient ⁇ ′.
  • (B) may be calculated.
  • the illuminance of the sun which is general natural light, can be used as the illuminance E ( ⁇ ) of the light source.
  • Expression (9) simulates the energy ratio of the assigned wavelength to the energy of light from the light source incident on the panchromatic sensor.
  • the transmittance correction value calculation unit 9 may calculate the weighting coefficient ⁇ ′ (b) using the panchromatic image of the actual image. For example, the transmittance correction value calculation unit 9 acquires a panchromatic image that is a real image in which a subject showing a known spectral radiance L sample and a dark area are shown. Then, the transmittance correction value calculation unit 9 calculates a difference dL between the known spectral radiance L sample and the spectral radiance in the dark part region, and the difference dL and the known spectral value are calculated as shown in the following equation (10). The transmittance ⁇ (Pa, x, y) is calculated from the radiance L sample .
  • the transmittance correction value calculation unit 9 uses the second correction value ⁇ (Pa ′, i, j) obtained from Equation (8) as the transmittance ⁇ (Pa ′, x, y) that simulates the multispectral image. Assuming that in equation (8), the weighting coefficient ⁇ ′ (b) is calculated so that the difference between the transmittance ⁇ (Pa ′, x, y) and the transmittance ⁇ (Pa, x, y) is minimized.
  • the image recombination unit 10 of the correction value calculation unit 7 includes the first correction value L scat (Pa ′, i) of the spectral radiance in each small region S (i, j) output from the luminance correction value calculation unit 8. , J) are combined to generate a two-dimensional array A of first correction values (step ST13 in FIG. 4). Further, the image recombination unit 10 calculates the second correction value ⁇ (Pa ′, i, j) of the spectral radiance in each small region S (i, j) output from the transmittance correction value calculation unit 9. By combining, a two-dimensional array B of second correction values is generated (step ST14 in FIG. 4).
  • Each element (i, j) in the two-dimensional arrays A and B corresponds to each small area S (i, j) divided by the image dividing unit 1. If the number of divisions of the multispectral image is N and the number of small regions S (i, j) is N, the number of elements (i, j) in the two-dimensional array is also N.
  • Each element (i, j) in the two-dimensional array A has a first correction value L scat (Pa ′, i, j) of the spectral radiance in each small region S (i, j).
  • Each element (i, j) in the two-dimensional array B has a second correction value ⁇ (Pa ′, i, j) of the spectral radiance in each small region S (i, j). .
  • the image recombination unit 10 outputs the generated two-dimensional arrays A and B to the correction value interpolation unit 11.
  • the correction value interpolation unit 11 of the correction value calculation unit 7 includes an image recombination unit 10 so that each of the resolutions of the two-dimensional arrays A and B output from the image recombination unit 10 matches the resolution of the panchromatic image. 2 are interpolated (step ST15 in FIG. 4). For example, if the number of pixels in the x direction of the panchromatic image is Mx times the number of elements (i, j) in the x direction in the two-dimensional arrays A and B, the x direction of the two-dimensional arrays A and B is multiplied by Mx.
  • an interpolation process is preferable in which correction values between the elements of the two-dimensional arrays A and B are continuous.
  • an interpolation process such as bilinear or cubic conversion is performed. Can be considered.
  • interpolation processing such as a Bezier curve or a Spline curve may be performed so that correction values between elements of the two-dimensional arrays A and B are continuous and draw a curve.
  • the first correction value L scat (Pa ′, i, j) is converted into the first correction value L scat (Pa ′, x, y), and the first correction value L scat ( Pa ′, x, y) is output from the correction value interpolation unit 11 to the radiance correction unit 13.
  • the second correction value ⁇ (Pa ′, i, j) is converted into the second correction value ⁇ (Pa ′, x, y), and the second correction value ⁇ (Pa ′, x) after conversion.
  • Y) is output from the correction value interpolation unit 11 to the transmittance correction unit 14.
  • the correction value interpolation unit 11 exists around the element. Based on the correction value possessed by the element, the correction value possessed by the element is interpolated, and then the resolution of the two-dimensional array A is converted. Further, when the second correction value ⁇ (Pa ′, i, j) of a certain element of the two-dimensional array B is an abnormal value, the correction value interpolation unit 11 exists around the element. Based on the correction value of the element, the correction value of the element is interpolated, and then the resolution of the two-dimensional array B is converted.
  • the radiance correction unit 13 of the correction unit 12 acquires a panchromatic image from the panchromatic sensor (step ST16 in FIG. 4).
  • the radiance correction unit 13 performs the first correction output from the correction value interpolation unit 11 from the spectral radiance L (Pa, x, y) of the acquired panchromatic image, as shown in the following equation (11).
  • the spectral radiance in the panchromatic image is corrected by subtracting the value L scat (Pa ′, x, y) (step ST17 in FIG. 4).
  • L1 (x, y) L (Pa, x, y) -Lscat (Pa ', x, y) (11)
  • L1 (x, y) is the spectral radiance after correction.
  • the radiance correction unit 13 outputs a panchromatic image whose spectral radiance is corrected to the transmittance correction unit 14.
  • the transmittance correction unit 14 of the correction unit 12 corrects the spectral radiance L1 (x, y) after correction of the panchromatic image output from the radiance correction unit 13 as shown in the following equation (12).
  • the spectral radiance in the panchromatic image is further corrected (step ST18 in FIG. 4).
  • L2 (x, y) is the spectral radiance after correction.
  • the transmittance correction unit 14 outputs a panchromatic image with the spectral radiance corrected.
  • a correction value calculation unit 7 that calculates a correction value of spectral radiance in a panchromatic image that is an image is provided, and the correction unit 12 uses the correction value of spectral radiance calculated by the correction value calculation unit 7 to use panchromatic. Since it is configured to correct the spectral radiance in the image, even if the panchromatic image acquired from the panchromatic sensor is an image including a plurality of regions having greatly different reflectances, the influence of haze or fog is affected. There is an effect that it is possible to correct the lowered contrast.
  • the correction value calculation unit 7 calculates the spectral radiance of the panchromatic image from the atmospheric transmittance in each region calculated by the atmospheric transmittance calculation unit 6 for each wavelength band. 2 is calculated, and the correction unit 12 is configured to correct the spectral radiance in the panchromatic image using the second correction value. Therefore, only the first correction value is used to correct the panchromatic image. Compared to correcting spectral radiance in an image, the contrast correction accuracy can be improved.
  • Embodiment 2 uses the first correction value L scat (Pa ′, x, y) and the second correction value ⁇ (Pa ′, x, y) in the panchromatic image.
  • amends spectral radiance is shown.
  • the third correction value is calculated to remove the offset component included in the panchromatic image, and the spectral radiance in the panchromatic image is further calculated using the third correction value. An example of correction will be described.
  • FIG. 11 is a block diagram showing an image processing apparatus according to Embodiment 2 of the present invention
  • FIG. 12 is a hardware block diagram showing an image processing apparatus according to Embodiment 2 of the present invention.
  • the reflection luminance calculation unit 41 is realized by, for example, a reflection luminance calculation circuit 51 shown in FIG.
  • the reflected luminance calculation unit 41 uses the spectral radiance L sensor (b, i, j) in each small region S (i, j) for each band b of the multispectral image and the scattered luminance calculation processing unit 5 for each band b.
  • the atmospheric reflection radiance indicating the atmospheric reflection in each small region S (i, j) from the atmospheric scattering radiance L scat (b, i, j) in each small region S (i, j) calculated in A process of calculating Lr (b, i, j) is performed.
  • the correction value calculation unit 42 is realized by, for example, the correction value calculation circuit 52 illustrated in FIG. 2, and includes a luminance correction value calculation unit 8, a transmittance correction value calculation unit 9, a reflection correction value calculation unit 43, and an image recombination. Unit 44 and correction value interpolation unit 45.
  • the reflection correction value calculation unit 43 of the correction value calculation unit 42 is the atmospheric reflection radiance Lr (b, i, j) in each small region S (i, j) calculated for each band b by the reflection luminance calculation unit 41. From this, a process of calculating the third correction value Lr (Pa ′, i, j) of the spectral radiance in each small region is performed.
  • the image recombination unit 44 of the correction value calculation unit 42 performs the process of generating the two-dimensional array A and the two-dimensional array B, similarly to the image recombination unit 10 of the first embodiment. Further, the image recombination unit 44 combines the third correction values Lr (Pa ′, i, j) of the spectral radiance in the respective small regions S (i, j) calculated by the reflection correction value calculation unit 43. Then, a process of generating a two-dimensional array C of third correction values is performed. Each element in the two-dimensional array C corresponds to each small region S (i, j) divided by the image dividing unit 1.
  • Each element in the two-dimensional array C has a third correction value of the spectral radiance in each small area S (i, j) divided by the image dividing unit 1.
  • the correction value interpolation unit 45 of the correction value calculation unit 42 includes an image recombination unit 44 so that the resolution of the two-dimensional arrays A, B, and C generated by the image recombination unit 44 matches the resolution of the panchromatic image. The process of interpolating the two-dimensional arrays A, B, and C generated by the above is performed.
  • the correction unit 46 is realized by the correction circuit 53 illustrated in FIG. 12, and includes a radiance correction unit 13, a transmittance correction unit 14, and a reflection correction unit 47.
  • the reflection correction unit 47 of the correction unit 46 performs a process of correcting the spectral radiance in the panchromatic image using the third correction value of each element in the two-dimensional array C interpolated by the correction value interpolation unit 45. carry out.
  • the image dividing unit 1, the parameter storage unit 2, the scattered luminance calculation unit 3, the atmospheric transmittance calculation unit 6, the reflection luminance calculation unit 41, the correction value calculation unit 42, and the correction unit 46 that are components of the image processing apparatus.
  • the image dividing circuit 21, the parameter storage circuit 22, the scattering luminance calculation circuit 23, the transmittance calculation circuit 24, the reflection luminance calculation circuit 51, the correction value calculation circuit 52, and the correction circuit 53 are realized.
  • the image dividing circuit 21, the scattered luminance calculating circuit 23, the transmittance calculating circuit 24, the reflected luminance calculating circuit 51, the correction value calculating circuit 52, and the correcting circuit 53 are, for example, a single circuit, a composite circuit, a programmed processor, a parallel program An integrated processor, ASIC, FPGA, or a combination thereof is applicable.
  • the components of the image processing apparatus are not limited to those realized by dedicated hardware, and the image processing apparatus may be realized by software, firmware, or a combination of software and firmware.
  • the parameter storage unit 2 is configured on the memory 31 of the computer shown in FIG. 3, and the image division unit 1, the scattered luminance calculation unit 3, and the atmospheric transmittance calculation unit 6.
  • a program for causing the computer to execute the processing procedure of the reflection luminance calculation unit 41, the correction value calculation unit 42, and the correction unit 46 is stored in the memory 31, and the processor 32 of the computer executes the program stored in the memory 31.
  • What should I do? 12 shows an example in which each component of the image processing apparatus is realized by dedicated hardware
  • FIG. 3 shows an example in which the image processing apparatus is realized by software, firmware, or the like. Some components in the image processing apparatus may be realized by dedicated hardware, and the remaining components may be realized by software, firmware, or the like.
  • the concentration of haze or fog is high, in addition to atmospheric scattering in which light is scattered to atmospheric particles that exist between the subject and the sensor, atmospheric reflection in which light is reflected to atmospheric particles occurs.
  • atmospheric reflection occurs, a reflected light component from the atmosphere is superimposed on the multispectral image and panchromatic image as an offset component.
  • the reflection correction unit 47 of the correction unit 46 performs a process of removing the offset component superimposed on the panchromatic image.
  • the reflection luminance calculation unit 41 acquires the spectral radiance L sensor (b, i, j) in each small region S (i, j) for each band b of the multispectral image divided by the image dividing unit 1.
  • the reflection luminance calculation unit 41 calculates the spectral radiance L sensor (b, i, j) in each small region S (i, j) by the scattered luminance calculation processing unit 5 as shown in the following equation (13).
  • FIG. 13 is an explanatory diagram showing the spectral radiance L sensor (b, i, j) modeled by the reflected luminance calculation unit 41.
  • the solid line is a model of spectral radiance L sensor (b, i, j)
  • ⁇ is the spectral radiance of the dark area in the small area S (i, j).
  • Lr (b, i, j) ⁇ (b) L r0 (i, j) (14)
  • ⁇ (b) is a spectral characteristic component of reflected light
  • L r0 (i, j) is an eigenvalue of the reflected light component.
  • the spectral characteristic component ⁇ (b) of the reflected light can be calculated from the spectral radiance of the white subject, as shown in the following equation (15).
  • ⁇ ( ⁇ ) is the reflectance of the white subject
  • E ( ⁇ ) is the illuminance of the sun as the light source. Since the atmospheric scattering radiance L scat (b, i, j) in Expression (13) is expressed as Expression (1), Expression (13) is expressed as Expression (16) below.
  • L sensor (b, i, j) ⁇ (b) L scat0 (i, j) + ⁇ (b) L r0 (i, j) (16)
  • the reflected luminance calculation unit 41 estimates the eigenvalue L r0 (i, j) of the reflected light component by regression analysis in the model of Expression (16). Reflection luminance calculation unit 41, the eigenvalues L r0 (i, j) of the reflected light component when estimating the eigenvalues L r0 estimated (i, j) and spectral characteristics ingredient beta (b) and the expression for the reflected light (14 ) To calculate the atmospheric reflection radiance Lr (b, i, j) in each small region S (i, j) of the multispectral image. The reflection luminance calculation unit 41 outputs the calculated atmospheric reflection radiance Lr (b, i, j) in each small region S (i, j) to the reflection correction value calculation unit 43.
  • the reflection correction value calculation unit 43 of the correction value calculation unit 42 acquires the wavelength characteristic data stored in the parameter storage unit 2.
  • the reflection correction value calculation unit 43 uses the acquired wavelength characteristic data to reflect the atmospheric reflection radiance Lr (b, i) in each small region S (i, j) for each band b output from the reflection luminance calculation unit 41. , J) is calculated.
  • the weighting coefficient ⁇ ′′ (b) can be calculated in the same manner as the weighting coefficient ⁇ (b) or the weighting coefficient ⁇ ′ (b), for example.
  • the image recombination unit 44 of the correction value calculation unit 42 includes the two-dimensional array A of the first correction values and the two-dimensional array B of the second correction values. Generate each. Further, the image recombination unit 44 combines the third correction values Lr (Pa ′, i, j) of the spectral radiances output from the reflection correction value calculation unit 43 in the respective small regions S (i, j). Then, a two-dimensional array C of third correction values is generated. Each element (i, j) in the two-dimensional array C corresponds to each small region S (i, j) divided by the image dividing unit 1.
  • the number of elements (i, j) of the two-dimensional array C is also N.
  • Each element (i, j) in the two-dimensional array C has a third correction value Lr (Pa ′, i, j) of the spectral radiance in each small region S (i, j).
  • the image recombination unit 44 outputs the generated two-dimensional arrays A, B, and C to the correction value interpolation unit 45.
  • the correction value interpolation unit 45 of the correction value calculation unit 42 interpolates each of the two-dimensional arrays A and B output from the image recombination unit 44, similarly to the correction value interpolation unit 11 of the first embodiment.
  • the correction value interpolation unit 45 interpolates the two-dimensional array C output from the image recombination unit 44. For example, if the number of pixels in the x direction of the panchromatic image is Mx times the number of elements (i, j) in the x direction in the two-dimensional array C, the x direction of the two-dimensional array C is upsampled to Mx times. .
  • the third correction value Lr (Pa ′, i, j) is converted into the third correction value Lr (Pa ′, x, y), and the third correction value Lr (Pa ′, after conversion) is converted.
  • x, y) is output from the correction value interpolation unit 45 to the reflection correction unit 47.
  • the reflection correction unit 47 of the correction unit 46 calculates a correction value from the corrected spectral radiance L1 (x, y) of the panchromatic image output from the radiance correction unit 13 as shown in the following equation (18).
  • L3 (x, y) L1 (x, y) ⁇ Lr (Pa ′, x, y) (18)
  • L3 (x, y) is the spectral radiance after correction.
  • the reflection correction unit 47 outputs a panchromatic image whose spectral radiance is corrected to the transmittance correction unit 14.
  • the reflection correction unit 47 is provided at the subsequent stage of the radiance correction unit 13, but the reflection correction unit 47 may be provided at the previous stage of the radiance correction unit 13. .
  • the correction value calculation unit 42 uses the atmospheric reflection radiance in each region calculated by the reflection luminance calculation unit 41 for each wavelength band, in the panchromatic image.
  • the correction unit 46 corrects the spectral radiance in the panchromatic image using the first and second correction values calculated by the correction value calculation unit 42.
  • the third correction value calculated by the correction value calculation unit 42 is used to correct the spectral radiance in the panchromatic image, the same effects as in the first and second embodiments can be obtained.
  • Embodiment 3 FIG. In the third embodiment, an example is described in which spectral radiance in a multispectral image is also corrected, and color synthesis is performed between a panchromatic image with corrected spectral radiance and a multispectral image with corrected spectral radiance. To do.
  • FIG. 14 is a block diagram showing an image processing apparatus according to Embodiment 3 of the present invention
  • FIG. 15 is a hardware block diagram showing an image processing apparatus according to Embodiment 3 of the present invention.
  • the same reference numerals as those in FIGS. 11 and 12 indicate the same or corresponding parts, and thus description thereof is omitted.
  • the correction value calculation unit 42 calculates the correction value of the spectral radiance in the panchromatic image as well as the correction value of the spectral radiance in the multispectral image, as in the second embodiment. Perform the calculation process.
  • the correction unit 46 performs a process of correcting the spectral radiance in the multispectral image in addition to correcting the spectral radiance in the panchromatic image, as in the second embodiment.
  • the color composition processing unit 48 is realized by the color composition processing circuit 54 shown in FIG.
  • the color composition processing unit 48 performs color composition processing of the panchromatic image whose spectral radiance is corrected by the correction unit 46 and the multispectral image whose spectral radiance is corrected by the correction unit 46.
  • FIG. 14 shows an example in which the color composition processing unit 48 is applied to the image processing apparatus in FIG. 11 in the second embodiment. However, the color composition processing unit 48 in FIG. It may be applied to the image processing apparatus.
  • each of the color composition processing unit 48 is realized by dedicated hardware as shown in FIG. That is, what is realized by the image dividing circuit 21, the parameter storage circuit 22, the scattering luminance calculation circuit 23, the transmittance calculation circuit 24, the reflection luminance calculation circuit 51, the correction value calculation circuit 52, the correction circuit 53, and the color composition processing circuit 54. Is assumed.
  • the image dividing circuit 21, the scattered luminance calculation circuit 23, the transmittance calculation circuit 24, the reflection luminance calculation circuit 51, the correction value calculation circuit 52, the correction circuit 53, and the color composition processing circuit 54 are, for example, a single circuit, a composite circuit, a program An integrated processor, a parallel-programmed processor, an ASIC, an FPGA, or a combination thereof is applicable.
  • the components of the image processing apparatus are not limited to those realized by dedicated hardware, and the image processing apparatus may be realized by software, firmware, or a combination of software and firmware.
  • the parameter storage unit 2 is configured on the memory 31 of the computer shown in FIG. 3, and the image division unit 1, the scattered luminance calculation unit 3, and the atmospheric transmittance calculation unit 6.
  • a program for causing the computer to execute the processing procedures of the reflection luminance calculation unit 41, the correction value calculation unit 42, the correction unit 46, and the color composition processing unit 48 is stored in the memory 31, and the processor 32 of the computer is stored in the memory 31.
  • the program that is running should be executed.
  • the brightness correction value calculation unit 8 uses the spectral radiance first correction value L scat (Pa ′, i, j) in each small region S (i, j) as an image recombination unit.
  • the atmospheric scattering radiance L scat (b, i, j) in each small region S (i, j) is output to the image recombining unit 44 as a fourth correction value.
  • the transmittance correction value calculation unit 9 uses the second correction value ⁇ (Pa ′, i, j) of the spectral radiance in each small region S (i, j) as an image recombination unit.
  • the atmospheric transmittance ⁇ (b, i, j) in each small region S (i, j) is output to the image recombining unit 44 as a fifth correction value.
  • the reflection correction value calculation unit 43 uses the image recombination unit 44 as the third correction value Lr (Pa ′, i, j) of the spectral radiance in each small region S (i, j).
  • the atmospheric reflection radiance Lr (b, i, j) in each small region S (i, j) is output to the image recombining unit 44 as a sixth correction value.
  • the image recombining unit 44 of the correction value calculating unit 42 generates the first two-dimensional array A of correction values, and outputs each of the output from the luminance correction value calculating unit 8.
  • L scat (b, i, j) which is the fourth correction value of spectral radiance in the small region S (i, j), is combined to generate a two-dimensional array D of fourth correction values.
  • the image recombination unit 44 generates the second correction value two-dimensional array B, and in addition to the spectral radiance in each small region S (i, j) output from the transmittance correction value calculation unit 9.
  • the image recombining unit 44 generates the two-dimensional array C of the third correction values, and the spectral radiance in each small region S (i, j) output from the reflection correction value calculating unit 43.
  • the six-correction atmospheric reflection radiance Lr (b, i, j) is combined to generate a two-dimensional array F of sixth correction values.
  • the correction value interpolation unit 45 of the correction value calculation unit 42 interpolates each of the two-dimensional arrays A, B, and C output from the image recombination unit 44 as in the second embodiment. Further, the correction value interpolation unit 45 interpolates each of the two-dimensional arrays D, E, and F output from the image recombination unit 44. As a result, the fourth correction value L scat (b, i, j) is converted to L scat (b, x, y), and the converted fourth correction value L scat (b, x) , Y) is output from the correction value interpolation unit 45 to the radiance correction unit 13.
  • the atmospheric transmittance ⁇ (b, i, j) that is the fifth correction value is converted into ⁇ (b, x, y), and ⁇ (b, x, y) that is the fifth correction value after conversion.
  • the atmospheric reflection radiance Lr (b, i, j) that is the sixth correction value is converted into Lr (b, x, y), and Lr (b, x, y) that is the sixth correction value after conversion. ) Is output from the correction value interpolation unit 45 to the reflection correction unit 47.
  • the radiance correction unit 13 corrects the spectral radiance in the panchromatic image as in the first and second embodiments. Further, the radiance correction unit 13 is output from the correction value interpolation unit 45 from the spectral radiance L sensor (b, i, j) in each small region S (i, j) for each band b of the multispectral image. The spectral radiance in the multispectral image is corrected by subtracting L scat (b, x, y), which is the fourth correction value.
  • the reflection correction unit 47 corrects the spectral radiance in the panchromatic image as in the second embodiment.
  • the reflection correction unit 47 corrects the spectral radiance L sensor (b, i, j) in each small region S (i, j) after correction of the multispectral image output from the radiance correction unit 13.
  • the spectral radiance in the multispectral image is corrected by subtracting Lr (b, x, y), which is the sixth correction value output from the value interpolation unit 45.
  • the transmittance correction unit 14 corrects the spectral radiance in the panchromatic image as in the first and second embodiments. Further, the transmittance correction unit 14 corrects the spectral radiance L sensor (b, i, j) in each small region S (i, j) after correction of the multispectral image output from the reflection correction unit 47. By dividing by the fifth correction value ⁇ (b, x, y) output from the value interpolation unit 45, the spectral radiance in the multispectral image is corrected.
  • the color composition processing unit 48 performs color composition processing of the panchromatic image whose spectral radiance has been corrected by the correction unit 46 and the multispectral image whose spectral radiance has been corrected by the correction unit 46, and after color synthesis processing Output the image.
  • a color composition process called IHS (Intensity Hue Saturation) conversion can be used as the color composition process.
  • IHS conversion luminance, hue, and saturation, which are color spaces, are calculated from a multispectral image, and the calculated luminance is replaced with a panchromatic image, thereby obtaining luminance information having the resolution of the panchromatic image. , Adding hue and saturation information.
  • brevey conversion As the color composition processing, brevey conversion, Gram-Schmidt conversion, conversion processing using principal component analysis, or the like can also be used.
  • An image that combines a high-resolution panchromatic image and a low-resolution multispectral image and has high resolution and color information may be referred to as a pan-sharpened image.
  • the present invention is suitable for an image processing apparatus that corrects spectral radiance in a panchromatic image.

Landscapes

  • Engineering & Computer Science (AREA)
  • Multimedia (AREA)
  • Signal Processing (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Image Processing (AREA)
  • Facsimile Image Signal Circuits (AREA)

Abstract

散乱輝度算出部（３）により波長帯毎に算出された各々の領域における大気散乱放射輝度から、複数の波長帯を含む単一波長帯の画像であるパンクロマチック画像における分光放射輝度の補正値を算出する補正値算出部（７）を設け、補正部（１２）が、補正値算出部（７）により算出された分光放射輝度の補正値を用いて、パンクロマチック画像における分光放射輝度を補正する。

Description

画像処理装置
　この発明は、パンクロマチック画像における分光放射輝度を補正する画像処理装置に関するものである。
　センサから取得した画像の一部分のコントラストが、大気中に含まれる靄又は霧などの影響を受けて低下していることがある。
　コントラストの低下は、画像の視認性の低下を招くため、靄又は霧などの影響で低下しているコントラストを補正する画像処理装置が開発されている。
　以下の特許文献１には、画像を構成する複数の画素のうち、任意の画素の信号レベルと、当該画素の周囲の領域に存在している画素の信号レベルとの差を計測し、その差が大きければ、コントラストが低下しているとみなして、当該画素の信号レベルを高める画像処理装置が開示されている。
特開２０１２－０２７５４７号公報
　例えば、センサから取得した画像が海岸線を撮影した画像である場合、画像を構成する複数の画素のうち、任意の画素が海の領域に対応し、当該画素の周囲の領域に存在している画素が海岸の領域に対応していることがある。
　このとき、海の反射率と海岸の反射率とは大きな差異があるため、大気中に含まれる靄又は霧などの影響を受けていない場合でも、任意の画素の信号レベルと、当該画素の周囲の領域に存在している画素の信号レベルとの差が大きくなることがある。
　このような場合に、コントラストが低下しているとみなして、当該画素の信号レベルを高めてしまうと、却って、画像の視認性が低下する。
　このため、センサから取得した画像が海岸線を撮影した画像である場合のほか、山間部と市街地の境界等のように、反射率が大きく異なる複数の領域を含んでいる画像である場合、靄又は霧などの影響を受けて低下しているコントラストを補正して、画像の視認性を高めることが困難であるという課題があった。
　この発明は上記のような課題を解決するためになされたもので、センサから取得したパンクロマチック画像が、反射率が大きく異なる複数の領域を含んでいる画像である場合でも、靄又は霧などの影響を受けて低下しているコントラストを補正することができる画像処理装置を得ることを目的とする。
　この発明に係る画像処理装置は、複数の波長帯の画像を含んでいるマルチスペクトル画像を複数の領域に分割する画像分割部と、マルチスペクトル画像が含んでいる複数の画像の波長帯毎に、画像分割部により分割された各々の領域における大気散乱を示す大気散乱放射輝度をそれぞれ算出する散乱輝度算出部と、散乱輝度算出部により波長帯毎に算出された各々の領域における大気散乱放射輝度から、複数の波長帯を含む単一波長帯の画像であるパンクロマチック画像における分光放射輝度の補正値を算出する補正値算出部とを設け、補正部が、補正値算出部により算出された分光放射輝度の補正値を用いて、パンクロマチック画像における分光放射輝度を補正するようにしたものである。
　この発明によれば、散乱輝度算出部により波長帯毎に算出された各々の領域における大気散乱放射輝度から、複数の波長帯を含む単一波長帯の画像であるパンクロマチック画像における分光放射輝度の補正値を算出する補正値算出部を設け、補正部が、補正値算出部により算出された分光放射輝度の補正値を用いて、パンクロマチック画像における分光放射輝度を補正するように構成したので、センサから取得したパンクロマチック画像が、反射率が大きく異なる複数の領域を含んでいる画像である場合でも、靄又は霧などの影響を受けて低下しているコントラストを補正することができる効果がある。
この発明の実施の形態１による画像処理装置を示す構成図である。 この発明の実施の形態１による画像処理装置を示すハードウェア構成図である。 画像処理装置がソフトウェア又はファームウェアなどで実現される場合のコンピュータのハードウェア構成図である。 画像処理装置がソフトウェア又はファームウェアなどで実現される場合の処理手順を示すフローチャートである。 画像分割部１により分割されたマルチスペクトル画像を示す説明図である。 画像分割部１により分割されたマルチスペクトル画像を示す説明図である。 暗部輝度算出部４により算出されたバンドｂにおける小領域Ｓ（ｉ，ｊ）のヒストグラムを示す説明図である。 散乱輝度算出処理部５によりモデル化された大気散乱放射輝度Ｌｓｃａｔ（ｂ，ｉ，ｊ）を示す説明図である。 バンドｂにおける小領域Ｓ（ｉ，ｊ）の大気散乱放射輝度Ｌｓｃａｔ（ｂ，ｉ，ｊ）に対応する大気透過率τ（ｂ，ｉ，ｊ）を示す説明図である。 マルチスペクトルセンサの観測可能な波長帯及びパンクロマチックセンサの観測可能な波長帯を示す説明図である。 この発明の実施の形態２による画像処理装置を示す構成図である。 この発明の実施の形態２による画像処理装置を示すハードウェア構成図である。 反射輝度算出部４１によりモデル化された分光放射輝度Ｌｓｅｎｓｏｒ（ｂ，ｉ，ｊ）を示す説明図である。 この発明の実施の形態３による画像処理装置を示す構成図である。 この発明の実施の形態３による画像処理装置を示すハードウェア構成図である。
　以下、この発明をより詳細に説明するために、この発明を実施するための形態について、添付の図面に従って説明する。
実施の形態１．
　図１は、この発明の実施の形態１による画像処理装置を示す構成図であり、図２は、この発明の実施の形態１による画像処理装置を示すハードウェア構成図である。
　図１及び図２において、画像分割部１は、例えば図２に示す画像分割回路２１で実現される。
　画像分割部１は、マルチスペクトルセンサから複数の波長帯の画像を含んでいるマルチスペクトル画像を取得し、マルチスペクトル画像を複数の小領域（領域）に分割する処理を実施する。
　パラメータ記憶部２は、例えば図２に示すパラメータ記憶回路２２で実現される。
　パラメータ記憶部２は、各々の波長帯の散乱特性を示す散乱特性データ、小領域の大気透過率を示す関数である散乱透過特性データ及びセンサの観測可能な波長帯を示す波長特性データを記憶している。
　散乱輝度算出部３は、例えば図２に示す散乱輝度算出回路２３で実現され、暗部輝度算出部４及び散乱輝度算出処理部５を備えている。
　散乱輝度算出部３は、パラメータ記憶部２により記憶されている散乱特性データを用いて、マルチスペクトル画像が含んでいる複数の画像の波長帯毎に、画像分割部１により分割された各々の小領域における大気散乱を示す大気散乱放射輝度をそれぞれ算出する処理を実施する。
　散乱輝度算出部３の暗部輝度算出部４は、マルチスペクトル画像が含んでいる複数の画像の波長帯毎に、画像分割部１により分割された各々の小領域のヒストグラムをそれぞれ算出し、そのヒストグラムから、画像分割部１により分割された各々の小領域の中で相対的に分光放射輝度が小さい暗部領域の分光放射輝度をそれぞれ算出する処理を実施する。
　散乱輝度算出部３の散乱輝度算出処理部５は、パラメータ記憶部２により記憶されている散乱特性データ及び暗部輝度算出部４により算出された暗部領域の分光放射輝度を用いて、マルチスペクトル画像が含んでいる複数の画像の波長帯毎に、画像分割部１により分割された各々の小領域における大気散乱放射輝度をそれぞれ算出する処理を実施する。
　大気透過率算出部６は、例えば図２に示す透過率算出回路２４で実現される。
　大気透過率算出部６は、パラメータ記憶部２により記憶されている散乱透過特性データを用いて、散乱輝度算出処理部５により波長帯毎に算出された各々の小領域における大気散乱放射輝度から、波長帯毎に、画像分割部１により分割された各々の小領域における大気散乱の大気透過率をそれぞれ算出する処理を実施する。
　補正値算出部７は、例えば図２に示す補正値算出回路２５で実現されるものであり、輝度補正値算出部８、透過率補正値算出部９、画像再結合部１０及び補正値補間部１１を備えている。
　補正値算出部７は、散乱輝度算出処理部５により波長帯毎に算出された各々の小領域における大気散乱放射輝度から、複数の波長帯を含む単一波長帯の画像であるパンクロマチック画像における分光放射輝度の補正値を算出する処理を実施する。
　また、補正値算出部７は、大気透過率算出部６により波長帯毎に算出された各々の小領域における大気散乱の大気透過率から、パンクロマチック画像における分光放射輝度の補正値を算出する処理を実施する。
　パンクロマチック画像は、マルチスペクトル画像が含んでいる複数の画像の波長帯を含む単一波長帯の画像である。
　補正値算出部７の輝度補正値算出部８は、パラメータ記憶部２により記憶されている波長特性データを用いて、散乱輝度算出処理部５により波長帯毎に算出された各々の小領域における大気散乱放射輝度から、各々の小領域における分光放射輝度の第１の補正値をそれぞれ算出する処理を実施する。
　補正値算出部７の透過率補正値算出部９は、パラメータ記憶部２により記憶されている波長特性データを用いて、大気透過率算出部６により波長帯毎に算出された各々の小領域における大気散乱の大気透過率から、各々の小領域における分光放射輝度の第２の補正値をそれぞれ算出する処理を実施する。
　補正値算出部７の画像再結合部１０は、輝度補正値算出部８により算出された各々の小領域における分光放射輝度の第１の補正値を結合して、第１の補正値の２次元配列Ａを生成する処理を実施する。
　また、画像再結合部１０は、透過率補正値算出部９により算出された各々の小領域における分光放射輝度の第２の補正値を結合して、第２の補正値の２次元配列Ｂを生成する処理を実施する。
　２次元配列Ａ，Ｂにおける各々の要素は、画像分割部１により分割された各々の小領域に対応している。２次元配列Ａにおける各々の要素は、画像分割部１により分割された各々の小領域における分光放射輝度の第１の補正値を有しており、２次元配列Ｂにおける各々の要素は、画像分割部１により分割された各々の小領域における分光放射輝度の第２の補正値を有している。
　補正値算出部７の補正値補間部１１は、画像再結合部１０により生成された２次元配列Ａ，Ｂの分解能が、パンクロマチック画像の分解能と一致するように、画像再結合部１０により生成された２次元配列Ａ，Ｂを補間する処理を実施する。
　補正部１２は、図２に示す補正回路２６で実現されるものであり、放射輝度補正部１３及び透過率補正部１４を備えている。
　補正部１２の放射輝度補正部１３は、パンクロマチックセンサからパンクロマチック画像を取得する処理を実施する。
　放射輝度補正部１３は、補正値補間部１１により補間された２次元配列Ａにおける各々の要素が有する第１の補正値を用いて、パンクロマチック画像における分光放射輝度を補正する処理を実施する。
　透過率補正部１４は、補正値補間部１１により補間された２次元配列Ｂにおける各々の要素が有する第２の補正値を用いて、パンクロマチック画像における分光放射輝度を補正する処理を実施する。
　図１では、画像処理装置の構成要素である画像分割部１、パラメータ記憶部２、散乱輝度算出部３、大気透過率算出部６、補正値算出部７及び補正部１２のそれぞれが、図２に示すような専用のハードウェアで実現されるものを想定している。即ち、画像分割回路２１、パラメータ記憶回路２２、散乱輝度算出回路２３、透過率算出回路２４、補正値算出回路２５及び補正回路２６で実現されるものを想定している。
　ここで、パラメータ記憶回路２２は、例えば、ＲＡＭ（Ｒａｎｄｏｍ　Ａｃｃｅｓｓ　Ｍｅｍｏｒｙ）、ＲＯＭ（Ｒｅａｄ　Ｏｎｌｙ　Ｍｅｍｏｒｙ）、フラッシュメモリ、ＥＰＲＯＭ（Ｅｒａｓａｂｌｅ　Ｐｒｏｇｒａｍｍａｂｌｅ　Ｒｅａｄ　Ｏｎｌｙ　Ｍｅｍｏｒｙ）、ＥＥＰＲＯＭ（Ｅｌｅｃｔｒｉｃａｌｌｙ　Ｅｒａｓａｂｌｅ　Ｐｒｏｇｒａｍｍａｂｌｅ　Ｒｅａｄ　Ｏｎｌｙ　Ｍｅｍｏｒｙ）などの不揮発性又は揮発性の半導体メモリや、磁気ディスク、フレキシブルディスク、光ディスク、コンパクトディスク、ミニディスク、ＤＶＤ（Ｄｉｇｉｔａｌ　Ｖｅｒｓａｔｉｌｅ　Ｄｉｓｃ）などが該当する。
　また、画像分割回路２１、散乱輝度算出回路２３、透過率算出回路２４、補正値算出回路２５及び補正回路２６は、例えば、単一回路、複合回路、プログラム化したプロセッサ、並列プログラム化したプロセッサ、ＡＳＩＣ（Ａｐｐｌｉｃａｔｉｏｎ　Ｓｐｅｃｉｆｉｃ　Ｉｎｔｅｇｒａｔｅｄ　Ｃｉｒｃｕｉｔ）、ＦＰＧＡ（Ｆｉｅｌｄ－Ｐｒｏｇｒａｍｍａｂｌｅ　Ｇａｔｅ　Ａｒｒａｙ）、または、これらを組み合わせたものが該当する。
　画像処理装置の構成要素は、専用のハードウェアで実現されるものに限るものではなく、画像処理装置がソフトウェア、ファームウェア、または、ソフトウェアとファームウェアとの組み合わせで実現されるものであってもよい。
　ソフトウェア又はファームウェアはプログラムとして、コンピュータのメモリに格納される。コンピュータは、プログラムを実行するハードウェアを意味し、例えば、ＣＰＵ（Ｃｅｎｔｒａｌ　Ｐｒｏｃｅｓｓｉｎｇ　Ｕｎｉｔ）、中央処理装置、処理装置、演算装置、マイクロプロセッサ、マイクロコンピュータ、プロセッサ、ＤＳＰ（Ｄｉｇｉｔａｌ　Ｓｉｇｎａｌ　Ｐｒｏｃｅｓｓｏｒ）などが該当する。
　図３は、画像処理装置がソフトウェア又はファームウェアなどで実現される場合のコンピュータのハードウェア構成図である。
　画像処理装置がソフトウェア又はファームウェアなどで実現される場合、パラメータ記憶部２をコンピュータのメモリ３１上に構成するとともに、画像分割部１、散乱輝度算出部３、大気透過率算出部６、補正値算出部７及び補正部１２の処理手順をコンピュータに実行させるためのプログラムをメモリ３１に格納し、コンピュータのプロセッサ３２がメモリ３１に格納されているプログラムを実行するようにすればよい。
　図４は、画像処理装置がソフトウェア又はファームウェアなどで実現される場合の処理手順を示すフローチャートである。
　また、図２では、画像処理装置の構成要素のそれぞれが専用のハードウェアで実現される例を示し、図３では、画像処理装置がソフトウェアやファームウェアなどで実現される例を示しているが、画像処理装置における一部の構成要素が専用のハードウェアで実現され、残りの構成要素がソフトウェアやファームウェアなどで実現されるものであってもよい。
　次に動作について説明する。
　この実施の形態１では、画像分割部１が、マルチスペクトルセンサから複数の波長帯の画像を含んでいるマルチスペクトル画像を取得し、放射輝度補正部１３が、パンクロマチックセンサからパンクロマチック画像を取得する。
　このとき、マルチスペクトル画像が含んでいる複数の画像の波長帯のそれぞれは、帯域が狭く、パンクロマチック画像の波長帯は、マルチスペクトル画像が含んでいる複数の画像の波長帯を含んでいるものとする。
　以下、マルチスペクトル画像が含んでいる複数の画像の波長帯をバンドｂと称し、この実施の形態１では、説明の便宜上、マルチスペクトル画像が、バンドｂ１の画像、バンドｂ２の画像、バンドｂ３の画像及びバンドｂ４の画像を含んでいるものとする。
　画像分割部１は、マルチスペクトルセンサからマルチスペクトル画像を取得し（図４のステップＳＴ１）、図５に示すように、マルチスペクトル画像をＮ個の小領域Ｓ（ｉ，ｊ）に分割する（図４のステップＳＴ２）。
　図５は、画像分割部１により分割されたマルチスペクトル画像を示す説明図である。
　図５は、Ｎ＝１６の例を示しており、マルチスペクトル画像が１６個に分割されている。ｉ＝１，２，３，４、ｊ＝１，２，３，４である。
　ただし、マルチスペクトル画像の分割数Ｎは、任意に設定できるものとし、例えば、マルチスペクトル画像に含まれている被写体の大きさなどから決定される。
　具体的には、マルチスペクトル画像に含まれている被写体と、被写体の陰影などの暗部領域とが含まれる程度の小領域の大きさに分割されるように、マルチスペクトル画像の分割数Ｎを設定するようにする。
　マルチスペクトル画像の中に、様々な大きさの被写体が含まれている場合には、マルチスペクトル画像に含まれている複数の被写体の中で、最も大きな被写体の大きさなどからマルチスペクトル画像の分割数Ｎを設定するようにして、被写体が含まれていない小領域Ｓ（ｉ，ｊ）の数を減らすようにしてもよい。
　また、マルチスペクトル画像に靄又は霧などの濃度の変化を示す濃度分布が描画されている場合、図６に示すように、濃度分布において、濃度の極大値又は極小値を含む領域毎に、マルチスペクトル画像を分割するようにしてもよい。
　図６は、画像分割部１により分割されたマルチスペクトル画像を示す説明図である。
　図６では、図５と同様に、Ｎ＝１６の例を示しているが、Ｎ＜１６であって、濃度の極大値又は極小値を含む領域毎に分割された小領域Ｓ（ｉ，ｊ）が、図５のように分割された小領域Ｓ（ｉ，ｊ）よりも大きな領域であれば、図５のようにマルチスペクトル画像を分割する場合よりも、暗部輝度算出部４における暗部領域の分光放射輝度の算出数を低減することができる。
　また、マルチスペクトル画像において、濃度の変化率が変化する箇所を把握できるように、マルチスペクトル画像が描画されている場合、小領域Ｓ（ｉ，ｊ）が、濃度の変化率が変化する箇所を含むようにマルチスペクトル画像を分割するようにしてもよい。
　散乱輝度算出部３の暗部輝度算出部４は、マルチスペクトル画像が含んでいる複数の画像のバンドｂ毎に、画像分割部１により分割された各々の小領域Ｓ（ｉ，ｊ）のヒストグラムをそれぞれ算出する（図４のステップＳＴ３）。
　図７は、暗部輝度算出部４により算出されたバンドｂにおける小領域Ｓ（ｉ，ｊ）のヒストグラムを示す説明図である。
　暗部輝度算出部４は、マルチスペクトル画像が含んでいる複数の画像のバンドｂ毎に、画像分割部１により分割された各々の小領域Ｓ（ｉ，ｊ）の中で相対的に分光放射輝度が小さい暗部領域の分光放射輝度ＬＤ（ｂ，ｉ，ｊ）をそれぞれ算出する（図４のステップＳＴ４）。
　図７では、小領域Ｓ（ｉ，ｊ）のヒストグラムの中で、分光放射輝度の最小値を、暗部領域の分光放射輝度ＬＤ（ｂ，ｉ，ｊ）としている例を示している。
　ここでは、暗部輝度算出部４が、小領域Ｓ（ｉ，ｊ）のヒストグラムを算出する例を示しているが、マルチスペクトル画像を構成している複数の画素の中に、画素値を有していない画素が含まれている場合がある。このような場合、暗部輝度算出部４は、マルチスペクトル画像を構成している複数の画素から、画素値を有していない画素を除いて、小領域Ｓ（ｉ，ｊ）のヒストグラムを算出するようにしてもよい。
　また、暗部輝度算出部４は、マルチスペクトル画像を構成している複数の画素の中に、画素値が常にゼロ付近の値になる欠陥画素が含まれている場合、マルチスペクトル画像を構成している複数の画素から、欠陥画素を除いて、小領域Ｓ（ｉ，ｊ）のヒストグラムを算出するようにしてもよい。
　散乱輝度算出部３の散乱輝度算出処理部５は、パラメータ記憶部２により記憶されている散乱特性データα（ｂ）を取得する（図４のステップＳＴ５）。
　パラメータ記憶部２により記憶されている散乱特性データα（ｂ）は、バンドｂの散乱特性を示すデータであり、バンド固有の値である。
　散乱特性データα（ｂ）は、例えば、ＭＯＤＴＲＡＮ（ＭＯＤｅｒａｔｅ　ｒｅｓｏｌｕｔｉｏｎ　ａｔｍｏｓｐｈｅｒｉｃ　ＴＲＡＮｓｍｉｓｓｉｏｎ）などの光波大気伝搬計算シミュレータを用いて、各々の波長における大気散乱の放射輝度比から算出することができる。
　散乱輝度算出処理部５は、マルチスペクトル画像が含んでいる複数の画像のバンドｂ毎に、画像分割部１により分割された各々の小領域Ｓ（ｉ，ｊ）における大気散乱放射輝度Ｌｓｃａｔ（ｂ，ｉ，ｊ）を、以下の式（１）に示すように、散乱特性データα（ｂ）と、小領域Ｓ（ｉ，ｊ）の固有値Ｌｓｃａｔ０（ｉ，ｊ）との積でモデル化する。
Ｌｓｃａｔ（ｂ，ｉ，ｊ）＝α（ｂ）×Ｌｓｃａｔ０（ｉ，ｊ）　　　（１）
　大気散乱放射輝度Ｌｓｃａｔ（ｂ，ｉ，ｊ）は、バンドｂの大気散乱を示す分光放射輝度である。
　小領域Ｓ（ｉ，ｊ）の固有値Ｌｓｃａｔ０（ｉ，ｊ）は、波長依存性が無く、撮影環境に依存する小領域Ｓ（ｉ，ｊ）の固有値である。
　図８は、散乱輝度算出処理部５によりモデル化された大気散乱放射輝度Ｌｓｃａｔ（ｂ，ｉ，ｊ）を示す説明図である。
　図８において、実線は、大気散乱放射輝度Ｌｓｃａｔ（ｂ，ｉ，ｊ）のモデルであり、〇は、小領域Ｓ（ｉ，ｊ）内の暗部領域の分光放射輝度である。
　次に、散乱輝度算出処理部５は、マルチスペクトル画像が含んでいる複数の画像のバンドｂ毎に、暗部領域の分光放射輝度ＬＤ（ｂ，ｉ，ｊ）と、式（１）の右辺であるα（ｂ）×Ｌｓｃａｔ０（ｉ，ｊ）とが最も近くなるような固有値Ｌｓｃａｔ０（ｉ，ｊ）を回帰分析によって推定する。
　このとき、暗部領域の分光放射輝度ＬＤ（ｂ，ｉ，ｊ）の計測点は、マルチスペクトルセンサが有しているバンド数と同数である。また、回帰分析によるパラメータ決定手法の代表的な例として、最小二乗法などが挙げられる。
　散乱輝度算出処理部５は、固有値Ｌｓｃａｔ０（ｉ，ｊ）を推定すると、推定した固有値Ｌｓｃａｔ０（ｉ，ｊ）と散乱特性データα（ｂ）を式（１）に代入して、バンドｂ毎に、各々の小領域Ｓ（ｉ，ｊ）における大気散乱放射輝度Ｌｓｃａｔ（ｂ，ｉ，ｊ）をそれぞれ算出する（図４のステップＳＴ６）。
　散乱輝度算出処理部５は、算出した各々の小領域Ｓ（ｉ，ｊ）における大気散乱放射輝度Ｌｓｃａｔ（ｂ，ｉ，ｊ）を大気透過率算出部６及び輝度補正値算出部８に出力する。
　ただし、散乱輝度算出処理部５は、ある小領域Ｓ（ｉ，ｊ）における暗部領域の分光放射輝度ＬＤ（ｂ，ｉ，ｊ）の計測点と、大気散乱放射輝度Ｌｓｃａｔ（ｂ，ｉ，ｊ）のモデルとの相関が低い場合、当該小領域Ｓ（ｉ，ｊ）における大気散乱放射輝度Ｌｓｃａｔ（ｂ，ｉ，ｊ）を異常値として、例えば、負の値などを出力するようにする。
　これにより、分光放射輝度ＬＤ（ｂ，ｉ，ｊ）の計測点がモデルと類似していない暗部領域を含んでいる小領域Ｓ（ｉ，ｊ）における大気散乱放射輝度Ｌｓｃａｔ（ｂ，ｉ，ｊ）を除去することができる。
　なお、計測点とモデルの相関が低い場合としては、例えば、最小二乗法における回帰分析の決定係数であるＲ二乗値が閾値Ｒｔｈ以下の場合などが挙げられる。
　大気透過率算出部６は、パラメータ記憶部２により記憶されている散乱透過特性データｆ（ａ）を取得する（図４のステップＳＴ７）。ａ＝（ｂ，Ｌｓｃａｔ（ｂ，ｉ，ｊ））である。
　パラメータ記憶部２により記憶されている散乱透過特性データｆ（ａ）は、小領域Ｓ（ｉ，ｊ）における大気散乱を示す大気散乱放射輝度Ｌｓｃａｔ（ｂ，ｉ，ｊ）の大気透過率τ（ｂ，ｉ，ｊ）を示すバンドｂ毎の関数である。
　散乱透過特性データｆ（ａ）は、例えば、ＭＯＤＴＲＡＮなどの光波大気伝搬計算シミュレータによって算出することができる。
　図９は、バンドｂにおける小領域Ｓ（ｉ，ｊ）の大気散乱放射輝度Ｌｓｃａｔ（ｂ，ｉ，ｊ）に対応する大気透過率τ（ｂ，ｉ，ｊ）を示す説明図である。
　大気透過率算出部６は、以下の式（２）に示すように、散乱透過特性データｆ（ａ）の変数ａに、バンドｂと、散乱輝度算出処理部５から出力された小領域Ｓ（ｉ，ｊ）における大気散乱放射輝度Ｌｓｃａｔ（ｂ，ｉ，ｊ）とを代入することで、バンドｂ毎に、小領域Ｓ（ｉ，ｊ）における大気散乱の大気透過率τ（ｂ，ｉ，ｊ）を算出する（図４のステップＳＴ８）。
τ（ｂ，ｉ，ｊ）＝ｆ（ｂ，Ｌｓｃａｔ（ｂ，ｉ，ｊ））　　　（２）
　大気透過率算出部６は、算出した各々の小領域Ｓ（ｉ，ｊ）における大気散乱の大気透過率τ（ｂ，ｉ，ｊ）を透過率補正値算出部９に出力する。
　ただし、大気透過率算出部６は、ある小領域Ｓ（ｉ，ｊ）における大気透過率τ（ｂ，ｉ，ｊ）が想定値と数パーセント以上離れているような場合には、当該小領域Ｓ（ｉ，ｊ）における大気透過率τ（ｂ，ｉ，ｊ）を異常値として、例えば、負の値などを出力するようにする。
　補正値算出部７の輝度補正値算出部８は、パラメータ記憶部２により記憶されている波長特性データを取得する（図４のステップＳＴ９）。
　パラメータ記憶部２により記憶されている波長特性データは、マルチスペクトルセンサの観測可能な波長帯及びパンクロマチックセンサの観測可能な波長帯のそれぞれを示すデータである。
　輝度補正値算出部８は、取得した波長特性データを用いて、散乱輝度算出処理部５から出力されたバンドｂ毎の各々の小領域Ｓ（ｉ，ｊ）における大気散乱放射輝度Ｌｓｃａｔ（ｂ，ｉ，ｊ）の重み係数ε（ｂ）を算出する。重み係数ε（ｂ）の算出処理の詳細は後述する。
　輝度補正値算出部８は、各々の小領域Ｓ（ｉ，ｊ）における分光放射輝度の第１の補正値Ｌｓｃａｔ（Ｐａ’，ｉ，ｊ）として、以下の式（３）に示すように、重み係数ε（ｂ）を用いて、バンドｂ毎の各々の小領域Ｓ（ｉ，ｊ）における大気散乱放射輝度Ｌｓｃａｔ（ｂ，ｉ，ｊ）の総和を算出する（図４のステップＳＴ１０）。
Figure JPOXMLDOC01-appb-I000001
ｂ＝ｂ１，ｂ２，ｂ３，ｂ４
　Ｐａ’は、パンクロマチックセンサにより観測されるパンクロマチック画像に相当する画像の単一波長帯である。
　輝度補正値算出部８は、各々の小領域Ｓ（ｉ，ｊ）における分光放射輝度の第１の補正値Ｌｓｃａｔ（Ｐａ’，ｉ，ｊ）を画像再結合部１０に出力する。
　なお、輝度補正値算出部８は、散乱輝度算出処理部５から出力された小領域Ｓ（ｉ，ｊ）における大気散乱放射輝度Ｌｓｃａｔ（ｂ，ｉ，ｊ）が異常値である場合、異常値である大気散乱放射輝度Ｌｓｃａｔ（ｂ，ｉ，ｊ）を用いずに、小領域Ｓ（ｉ，ｊ）における分光放射輝度の第１の補正値Ｌｓｃａｔ（Ｐａ’，ｉ，ｊ）を算出するようにする。
　例えば、輝度補正値算出部８は、大気散乱放射輝度Ｌｓｃａｔ（ｂ，ｉ，ｊ）が異常値である小領域Ｓ（ｉ，ｊ）の周辺の小領域における大気散乱放射輝度を用いて、小領域Ｓ（ｉ，ｊ）における大気散乱放射輝度Ｌｓｃａｔ（ｂ，ｉ，ｊ）を補間する。
　そして、輝度補正値算出部８は、補間した大気散乱放射輝度Ｌｓｃａｔ（ｂ，ｉ，ｊ）を用いて、第１の補正値Ｌｓｃａｔ（Ｐａ’，ｉ，ｊ）を算出するようにする。
　ここで、重み係数ε（ｂ）の算出処理について説明する。
　図１０は、マルチスペクトルセンサの観測可能な波長帯及びパンクロマチックセンサの観測可能な波長帯を示す説明図である。
　図１０の例では、パンクロマチックセンサの観測可能な波長帯は、λ０（Ｐａ）～λ１（Ｐａ）である。
　マルチスペクトルセンサの観測可能な波長帯は、λ０（ｂ１）～λ１（ｂ１）、λ０（ｂ２）～λ１（ｂ２）、λ０（ｂ３）～λ１（ｂ３）及びλ０（ｂ４）～λ１（ｂ４）である。
　輝度補正値算出部８は、パンクロマチックセンサの観測可能な波長帯λ０（Ｐａ）～λ１（Ｐａ）に対して、マルチスペクトルセンサの観測可能な波長帯λ’０（ｂ）～λ’１（ｂ）を割り付けることで、以下の式（４）に示すように、重み係数ε（ｂ）を算出する。
Figure JPOXMLDOC01-appb-I000002
ｂ＝ｂ１，ｂ２，ｂ３，ｂ４
　ここで、輝度補正値算出部８におけるマルチスペクトルセンサの波長帯λ’０（ｂ）～λ’１（ｂ）の割り付けは、例えば、波長帯λ’０（ｂ）～λ’１（ｂ）が、マルチスペクトル画像が含んでいる画像のバンドｂ１，ｂ２，ｂ３，ｂ４のそれぞれと対応するように行う。
　即ち、輝度補正値算出部８は、バンドｂ１に対して波長帯λ’０（ｂ１）～λ’１（ｂ１）を割り付け、バンドｂ２に対して波長帯λ’０（ｂ２）～λ’１（ｂ２）を割り付け、バンドｂ３に対して波長帯λ’０（ｂ３）～λ’１（ｂ３）を割り付ける。
　また、輝度補正値算出部８は、バンドｂ４に対して波長帯λ’０（ｂ４）～λ’１（ｂ４）を割り付ける。
　このような割り付けの場合、例えば、波長帯λ’０（ｂ１），λ’１（ｂ１）は、以下の式（５）及び式（６）のように表される。
Figure JPOXMLDOC01-appb-I000003

Figure JPOXMLDOC01-appb-I000004
　重み係数ε（ｂ）における分子は、割り付けられた波長帯λ’０（ｂ）～λ’１（ｂ）における分光放射輝度に対する波長積分を模擬するものである。
　また、重み係数ε（ｂ）における分母は、パンクロマチックセンサの放射輝度から、分光放射輝度に変換すること意味している。
　上記の式（３）は、割り付けられた波長帯λ’０（ｂ）～λ’１（ｂ）における分光放射輝度を一定とみなしたものであり、パンクロマチックセンサの波長帯に対して、割り付けるバンドの数が多いほど、分光放射輝度の第１の補正値Ｌｓｃａｔ（Ｐａ’，ｉ，ｊ）を高精度に算出することができる。
　ここでは、輝度補正値算出部８が、式（４）によって、重み係数ε（ｂ）を算出する例を示しているが、以下の式（７）に示すように、パンクロマチックセンサにおける各波長λの分光感度Ｒ（λ）に基づいて重み係数ε（ｂ）を算出するようにしてもよい。この場合、パンクロマチックセンサが、各々のバンドｂの波長λで感度が著しく変化するようなセンサであっても、分光放射輝度の第１の補正値Ｌｓｃａｔ（Ｐａ’，ｉ，ｊ）を高精度に算出することができる。
Figure JPOXMLDOC01-appb-I000005
　また、輝度補正値算出部８は、実画像のパンクロマチック画像を用いて、重み係数ε（ｂ）を算出するようにしてもよい。
　例えば、輝度補正値算出部８は、実画像のパンクロマチック画像の分光放射輝度Ｌ（Ｐａ，ｘ，ｙ）を取得する。
　（ｘ，ｙ）は、パンクロマチック画像上の座標であり、マルチスペクトル画像における小領域の配列を示す（ｉ，ｊ）と異なる。
　そして、輝度補正値算出部８は、式（３）から得られる分光放射輝度の第１の補正値Ｌｓｃａｔ（Ｐａ’，ｉ，ｊ）を、マルチスペクトル画像から模擬したパンクロマチックセンサの分光放射輝度Ｌ（Ｐａ’，ｘ，ｙ）とみなし、式（３）において、分光放射輝度Ｌ（Ｐａ’，ｘ，ｙ）と分光放射輝度Ｌ（Ｐａ，ｘ，ｙ）との差が最小になるような重み係数ε（ｂ）を算出する。
　補正値算出部７の透過率補正値算出部９は、パラメータ記憶部２により記憶されている波長特性データを取得する（図４のステップＳＴ１１）。
　透過率補正値算出部９は、取得した波長特性データを用いて、大気透過率算出部６から出力されたバンドｂ毎の各々の小領域Ｓ（ｉ，ｊ）における大気散乱の大気透過率τ（ｂ，ｉ，ｊ）の重み係数ε’（ｂ）を算出する。重み係数ε’（ｂ）の算出処理の詳細は後述する。
　透過率補正値算出部９は、各々の小領域Ｓ（ｉ，ｊ）における分光放射輝度の第２の補正値τ（Ｐａ’，ｉ，ｊ）として、以下の式（８）に示すように、重み係数ε’（ｂ）を用いて、バンドｂ毎の各々の小領域Ｓ（ｉ，ｊ）における大気散乱の大気透過率τ（ｂ，ｉ，ｊ）の総和を算出する（図４のステップＳＴ１２）。
Figure JPOXMLDOC01-appb-I000006
ｂ＝ｂ１，ｂ２，ｂ３，ｂ４
　透過率補正値算出部９は、各々の小領域Ｓ（ｉ，ｊ）における分光放射輝度の第２の補正値τ（Ｐａ’，ｉ，ｊ）を画像再結合部１０に出力する。
　なお、透過率補正値算出部９は、大気透過率算出部６から出力された小領域Ｓ（ｉ，ｊ）における大気散乱の大気透過率τ（ｂ，ｉ，ｊ）が異常値である場合、異常値である大気透過率τ（ｂ，ｉ，ｊ）を用いずに、小領域Ｓ（ｉ，ｊ）における分光放射輝度の第２の補正値τ（Ｐａ’，ｉ，ｊ）を算出するようにする。
　例えば、透過率補正値算出部９は、大気透過率τ（ｂ，ｉ，ｊ）が異常値である小領域Ｓ（ｉ，ｊ）の周辺の小領域における大気透過率を用いて、小領域Ｓ（ｉ，ｊ）における大気散乱の大気透過率τ（ｂ，ｉ，ｊ）を補間する。
　そして、透過率補正値算出部９は、補間した大気透過率τ（ｂ，ｉ，ｊ）を用いて、第２の補正値τ（Ｐａ’，ｉ，ｊ）を算出するようにする。
　ここで、重み係数ε’（ｂ）の算出処理について説明する。
　大気透過率τ（ｂ，ｉ，ｊ）の重み係数ε’（ｂ）は、輝度補正値算出部８と同様に、式（４）によって算出することができる。
　透過率補正値算出部９は、パンクロマチックセンサにおける各波長λの分光感度Ｒ（λ）及び光源の照度Ｅ（λ）を用いて、以下の式（９）に示すように、重み係数ε’（ｂ）を算出するようにしてもよい。
Figure JPOXMLDOC01-appb-I000007
　光源の照度Ｅ（λ）として、一般的な自然光である太陽の照度を用いることができる。
　式（９）は、パンクロマチックセンサに入射される光源からの光のエネルギーに対する割り付け波長のエネルギー比を模擬している。
　また、透過率補正値算出部９は、実画像のパンクロマチック画像を用いて、重み係数ε’（ｂ）を算出するようにしてもよい。
　例えば、透過率補正値算出部９は、既知の分光放射輝度Ｌｓａｍｐｌｅを示す被写体と、暗部領域とが映っている実画像であるパンクロマチック画像を取得する。
　そして、透過率補正値算出部９は、既知の分光放射輝度Ｌｓａｍｐｌｅと暗部領域の分光放射輝度との差分ｄＬを算出し、以下の式（１０）に示すように、差分ｄＬと既知の分光放射輝度Ｌｓａｍｐｌｅとから透過率τ（Ｐａ，ｘ，ｙ）を算出する。
Figure JPOXMLDOC01-appb-I000008
　透過率補正値算出部９は、式（８）から得られる第２の補正値τ（Ｐａ’，ｉ，ｊ）を、マルチスペクトル画像から模擬した透過率τ（Ｐａ’，ｘ，ｙ）とみなし、式（８）において、透過率τ（Ｐａ’，ｘ，ｙ）と透過率τ（Ｐａ，ｘ，ｙ）との差が最小になるような重み係数ε’（ｂ）を算出する。
　補正値算出部７の画像再結合部１０は、輝度補正値算出部８から出力された各々の小領域Ｓ（ｉ，ｊ）における分光放射輝度の第１の補正値Ｌｓｃａｔ（Ｐａ’，ｉ，ｊ）を結合して、第１の補正値の２次元配列Ａを生成する（図４のステップＳＴ１３）。
　また、画像再結合部１０は、透過率補正値算出部９から出力された各々の小領域Ｓ（ｉ，ｊ）における分光放射輝度の第２の補正値τ（Ｐａ’，ｉ，ｊ）を結合して、第２の補正値の２次元配列Ｂを生成する（図４のステップＳＴ１４）。
　２次元配列Ａ，Ｂにおける各々の要素（ｉ，ｊ）は、画像分割部１により分割された各々の小領域Ｓ（ｉ，ｊ）に対応している。マルチスペクトル画像の分割数がＮで、小領域Ｓ（ｉ，ｊ）の個数がＮであれば、２次元配列の要素（ｉ，ｊ）の数もＮである。
　２次元配列Ａにおける各々の要素（ｉ，ｊ）は、各々の小領域Ｓ（ｉ，ｊ）における分光放射輝度の第１の補正値Ｌｓｃａｔ（Ｐａ’，ｉ，ｊ）を有しており、２次元配列Ｂにおける各々の要素（ｉ，ｊ）は、各々の小領域Ｓ（ｉ，ｊ）における分光放射輝度の第２の補正値τ（Ｐａ’，ｉ，ｊ）を有している。
　画像再結合部１０は、生成した２次元配列Ａ，Ｂを補正値補間部１１に出力する。
　補正値算出部７の補正値補間部１１は、画像再結合部１０から出力された２次元配列Ａ，Ｂの分解能のそれぞれが、パンクロマチック画像の分解能と一致するように、画像再結合部１０から出力された２次元配列Ａ，Ｂのそれぞれを補間する（図４のステップＳＴ１５）。
　例えば、パンクロマチック画像のｘ方向の画素数が、２次元配列Ａ，Ｂにおけるｘ方向の要素（ｉ，ｊ）の数のＭｘ倍であれば、２次元配列Ａ，Ｂのｘ方向をＭｘ倍にアップサンプリングする。
　また、パンクロマチック画像のｙ方向の画素数が、２次元配列Ａ，Ｂにおけるｙ方向の要素（ｉ，ｊ）の数のＭｙ倍である場合、２次元配列Ａ，Ｂのｙ方向をＭｙ倍にアップサンプリングする。
　これにより、２次元配列Ａ，Ｂの分解能が、パンクロマチック画像の分解能と同等になる。
　２次元配列Ａ，Ｂの補間処理としては、２次元配列Ａ，Ｂの要素間の補正値が連続的となるような補間処理が望ましく、例えば、バイリニア又はＣｕｂｉｃ　Ｃｏｎｖｏｌｕｔｉｏｎなどの補間処理を実施することが考えられる。また、２次元配列Ａ，Ｂの要素間の補正値が連続的かつ曲線を描くように、ベジエ曲線又はＳｐｌｉｎｅ曲線などの補間処理を実施するようにしてもよい。
　これにより、第１の補正値Ｌｓｃａｔ（Ｐａ’，ｉ，ｊ）は、第１の補正値Ｌｓｃａｔ（Ｐａ’，ｘ，ｙ）に変換され、変換後の第１の補正値Ｌｓｃａｔ（Ｐａ’，ｘ，ｙ）が補正値補間部１１から放射輝度補正部１３に出力される。また、第２の補正値τ（Ｐａ’，ｉ，ｊ）は、第２の補正値τ（Ｐａ’，ｘ，ｙ）に変換され、変換後の第２の補正値τ（Ｐａ’，ｘ，ｙ）が補正値補間部１１から透過率補正部１４に出力される。
　なお、２次元配列Ａの或る要素が有している第１の補正値Ｌｓｃａｔ（Ｐａ’，ｉ，ｊ）が異常値である場合、補正値補間部１１は、当該要素の周囲に存在している要素が有している補正値に基づいて、当該要素が有している補正値を補間してから、２次元配列Ａの分解能を変換する。
　また、２次元配列Ｂの或る要素が有している第２の補正値τ（Ｐａ’，ｉ，ｊ）が異常値である場合、補正値補間部１１は、当該要素の周囲に存在している要素が有している補正値に基づいて、当該要素が有している補正値を補間してから、２次元配列Ｂの分解能を変換する。
　補正部１２の放射輝度補正部１３は、パンクロマチックセンサからパンクロマチック画像を取得する（図４のステップＳＴ１６）。
　放射輝度補正部１３は、以下の式（１１）に示すように、取得したパンクロマチック画像の分光放射輝度Ｌ（Ｐａ，ｘ，ｙ）から、補正値補間部１１から出力された第１の補正値Ｌｓｃａｔ（Ｐａ’，ｘ，ｙ）を減算することで、パンクロマチック画像における分光放射輝度を補正する（図４のステップＳＴ１７）。
Ｌ１（ｘ，ｙ）
＝Ｌ（Ｐａ，ｘ，ｙ）－Ｌｓｃａｔ（Ｐａ’，ｘ，ｙ）　　　（１１）
　式（１１）において、Ｌ１（ｘ，ｙ）は、補正後の分光放射輝度である。
　放射輝度補正部１３は、分光放射輝度を補正したパンクロマチック画像を透過率補正部１４に出力する。
　補正部１２の透過率補正部１４は、以下の式（１２）に示すように、放射輝度補正部１３から出力されたパンクロマチック画像の補正後の分光放射輝度Ｌ１（ｘ，ｙ）を、補正値補間部１１から出力された第２の補正値τ（Ｐａ’，ｘ，ｙ）で除算することで、さらに、パンクロマチック画像における分光放射輝度を補正する（図４のステップＳＴ１８）。
Figure JPOXMLDOC01-appb-I000009
　式（１２）において、Ｌ２（ｘ，ｙ）は、補正後の分光放射輝度である。
　透過率補正部１４は、分光放射輝度を補正したパンクロマチック画像を出力する。
　以上で明らかなように、この実施の形態１によれば、散乱輝度算出部３により波長帯毎に算出された各々の領域における大気散乱放射輝度から、複数の波長帯を含む単一波長帯の画像であるパンクロマチック画像における分光放射輝度の補正値を算出する補正値算出部７を設け、補正部１２が、補正値算出部７により算出された分光放射輝度の補正値を用いて、パンクロマチック画像における分光放射輝度を補正するように構成したので、パンクロマチックセンサから取得したパンクロマチック画像が、反射率が大きく異なる複数の領域を含んでいる画像である場合でも、靄又は霧などの影響を受けて低下しているコントラストを補正することができる効果を奏する。
　また、この実施の形態１によれば、補正値算出部７が、大気透過率算出部６により波長帯毎に算出された各々の領域における大気透過率から、パンクロマチック画像における分光放射輝度の第２の補正値を算出し、補正部１２が、第２の補正値を用いて、パンクロマチック画像における分光放射輝度を補正するように構成したので、第１の補正値だけを用いて、パンクロマチック画像における分光放射輝度を補正する場合よりも、コントラストの補正精度を高めることができる効果を奏する。
実施の形態２．
　上記実施の形態１では、補正部１２が、第１の補正値Ｌｓｃａｔ（Ｐａ’，ｘ，ｙ）及び第２の補正値τ（Ｐａ’，ｘ，ｙ）を用いて、パンクロマチック画像における分光放射輝度を補正する例を示している。
　この実施の形態２では、パンクロマチック画像に含まれているオフセット成分を除去するために第３の補正値を算出し、第３の補正値を用いて、さらに、パンクロマチック画像における分光放射輝度を補正する例を説明する。
　図１１は、この発明の実施の形態２による画像処理装置を示す構成図であり、図１２は、この発明の実施の形態２による画像処理装置を示すハードウェア構成図である。
　図１１及び図１２において、図１及び図２と同一符号は同一または相当部分を示すので説明を省略する。
　反射輝度算出部４１は、例えば図１２に示す反射輝度算出回路５１で実現される。
　反射輝度算出部４１は、マルチスペクトル画像のバンドｂ毎の各々の小領域Ｓ（ｉ，ｊ）における分光放射輝度Ｌｓｅｎｓｏｒ（ｂ，ｉ，ｊ）と、散乱輝度算出処理部５によりバンドｂ毎に算出された各々の小領域Ｓ（ｉ，ｊ）における大気散乱放射輝度Ｌｓｃａｔ（ｂ，ｉ，ｊ）とから、各々の小領域Ｓ（ｉ，ｊ）における大気反射を示す大気反射放射輝度Ｌｒ（ｂ，ｉ，ｊ）をそれぞれ算出する処理を実施する。
　補正値算出部４２は、例えば図２に示す補正値算出回路５２で実現されるものであり、輝度補正値算出部８、透過率補正値算出部９、反射補正値算出部４３、画像再結合部４４及び補正値補間部４５を備えている。
　補正値算出部４２の反射補正値算出部４３は、反射輝度算出部４１によりバンドｂ毎に算出された各々の小領域Ｓ（ｉ，ｊ）における大気反射放射輝度Ｌｒ（ｂ，ｉ，ｊ）から、各々の小領域における分光放射輝度の第３の補正値Ｌｒ（Ｐａ’，ｉ，ｊ）をそれぞれ算出する処理を実施する。
　補正値算出部４２の画像再結合部４４は、上記実施の形態１の画像再結合部１０と同様に、２次元配列Ａ及び２次元配列Ｂを生成する処理を実施する。
　また、画像再結合部４４は、反射補正値算出部４３により算出された各々の小領域Ｓ（ｉ，ｊ）における分光放射輝度の第３の補正値Ｌｒ（Ｐａ’，ｉ，ｊ）を結合して、第３の補正値の２次元配列Ｃを生成する処理を実施する。
　２次元配列Ｃにおける各々の要素は、画像分割部１により分割された各々の小領域Ｓ（ｉ，ｊ）に対応している。２次元配列Ｃにおける各々の要素は、画像分割部１により分割された各々の小領域Ｓ（ｉ，ｊ）における分光放射輝度の第３の補正値を有している。
　補正値算出部４２の補正値補間部４５は、画像再結合部４４により生成された２次元配列Ａ，Ｂ，Ｃの分解能が、パンクロマチック画像の分解能と一致するように、画像再結合部４４により生成された２次元配列Ａ，Ｂ，Ｃを補間する処理を実施する。
　補正部４６は、図１２に示す補正回路５３で実現されるものであり、放射輝度補正部１３、透過率補正部１４及び反射補正部４７を備えている。
　補正部４６の反射補正部４７は、補正値補間部４５により補間された２次元配列Ｃにおける各々の要素が有する第３の補正値を用いて、パンクロマチック画像における分光放射輝度を補正する処理を実施する。
　図１１では、画像処理装置の構成要素である画像分割部１、パラメータ記憶部２、散乱輝度算出部３、大気透過率算出部６、反射輝度算出部４１、補正値算出部４２及び補正部４６のそれぞれが、図１２に示すような専用のハードウェアで実現されるものを想定している。即ち、画像分割回路２１、パラメータ記憶回路２２、散乱輝度算出回路２３、透過率算出回路２４、反射輝度算出回路５１、補正値算出回路５２及び補正回路５３で実現されるものを想定している。
　画像分割回路２１、散乱輝度算出回路２３、透過率算出回路２４、反射輝度算出回路５１、補正値算出回路５２及び補正回路５３は、例えば、単一回路、複合回路、プログラム化したプロセッサ、並列プログラム化したプロセッサ、ＡＳＩＣ、ＦＰＧＡ、または、これらを組み合わせたものが該当する。
　画像処理装置の構成要素は、専用のハードウェアで実現されるものに限るものではなく、画像処理装置がソフトウェア、ファームウェア、または、ソフトウェアとファームウェアとの組み合わせで実現されるものであってもよい。
　画像処理装置がソフトウェア又はファームウェアなどで実現される場合、パラメータ記憶部２を図３に示すコンピュータのメモリ３１上に構成するとともに、画像分割部１、散乱輝度算出部３、大気透過率算出部６、反射輝度算出部４１、補正値算出部４２及び補正部４６の処理手順をコンピュータに実行させるためのプログラムをメモリ３１に格納し、コンピュータのプロセッサ３２がメモリ３１に格納されているプログラムを実行するようにすればよい。
　また、図１２では、画像処理装置の構成要素のそれぞれが専用のハードウェアで実現される例を示し、図３では、画像処理装置がソフトウェアやファームウェアなどで実現される例を示しているが、画像処理装置における一部の構成要素が専用のハードウェアで実現され、残りの構成要素がソフトウェアやファームウェアなどで実現されるものであってもよい。
　次に動作について説明する。
　靄又は霧などの濃度が高い場合、被写体とセンサの間に存在している大気中の粒子に光が散乱される大気散乱のほかに、大気中の粒子に光が反射される大気反射が生じることがある。
　大気反射が生じることで、大気からの反射光成分がオフセット成分としてマルチスペクトル画像及びパンクロマチック画像に重畳される。
　この実施の形態２では、補正部４６の反射補正部４７が、パンクロマチック画像に重畳されているオフセット成分を除去する処理を実施する。
　以下、上記実施の形態１と相違する箇所を説明する。
　反射輝度算出部４１は、画像分割部１により分割されたマルチスペクトル画像のバンドｂ毎の各々の小領域Ｓ（ｉ，ｊ）における分光放射輝度Ｌｓｅｎｓｏｒ（ｂ，ｉ，ｊ）を取得する。
　反射輝度算出部４１は、各々の小領域Ｓ（ｉ，ｊ）における分光放射輝度Ｌｓｅｎｓｏｒ（ｂ，ｉ，ｊ）を、以下の式（１３）に示すように、散乱輝度算出処理部５により算出された各々の小領域Ｓ（ｉ，ｊ）における大気散乱放射輝度Ｌｓｃａｔ（ｂ，ｉ，ｊ）と、各々の小領域Ｓ（ｉ，ｊ）における大気反射放射輝度Ｌｒ（ｂ，ｉ，ｊ）との和でモデル化する。
Ｌｓｅｎｓｏｒ（ｂ，ｉ，ｊ）
＝Ｌｓｃａｔ（ｂ，ｉ，ｊ）＋Ｌｒ（ｂ，ｉ，ｊ）　　　　（１３）
　図１３は、反射輝度算出部４１によりモデル化された分光放射輝度Ｌｓｅｎｓｏｒ（ｂ，ｉ，ｊ）を示す説明図である。
　図１３において、実線は、分光放射輝度Ｌｓｅｎｓｏｒ（ｂ，ｉ，ｊ）のモデルであり、〇は、小領域Ｓ（ｉ，ｊ）内の暗部領域の分光放射輝度である。
　ここで、大気からの反射光成分の分光特性が変化せずに、反射光成分の絶対値のみが変化するとすれば、マルチスペクトル画像の各々の小領域Ｓ（ｉ，ｊ）における大気反射放射輝度Ｌｒ（ｂ，ｉ，ｊ）は、以下の式（１４）のように表すことができる。
Ｌｒ（ｂ，ｉ，ｊ）＝β（ｂ）Ｌｒ０（ｉ，ｊ）　　　　　（１４）
　式（１４）において、β（ｂ）は、反射光の分光特性成分、Ｌｒ０（ｉ，ｊ）は、反射光成分の固有値である。
　反射光の分光特性成分β（ｂ）は、以下の式（１５）に示すように、白色の被写体が持つ分光放射輝度より算出することができる。
Figure JPOXMLDOC01-appb-I000010
　式（１５）において、ρ（λ）は、白色の被写体の反射率、Ｅ（λ）は、光源である太陽の照度である。
　式（１３）における大気散乱放射輝度Ｌｓｃａｔ（ｂ，ｉ，ｊ）は、式（１）のように表されるので、式（１３）は、以下の式（１６）のように表される。
Ｌｓｅｎｓｏｒ（ｂ，ｉ，ｊ）
＝α（ｂ）Ｌｓｃａｔ０（ｉ，ｊ）＋β（ｂ）Ｌｒ０（ｉ，ｊ）　　　（１６）
　反射輝度算出部４１は、式（１６）のモデルにおいて、反射光成分の固有値Ｌｒ０（ｉ，ｊ）を回帰分析によって推定する。
　反射輝度算出部４１は、反射光成分の固有値Ｌｒ０（ｉ，ｊ）を推定すると、推定した固有値Ｌｒ０（ｉ，ｊ）と、反射光の分光特性成分β（ｂ）とを式（１４）に代入して、マルチスペクトル画像の各々の小領域Ｓ（ｉ，ｊ）における大気反射放射輝度Ｌｒ（ｂ，ｉ，ｊ）を算出する。
　反射輝度算出部４１は、算出した各々の小領域Ｓ（ｉ，ｊ）における大気反射放射輝度Ｌｒ（ｂ，ｉ，ｊ）を反射補正値算出部４３に出力する。
　補正値算出部４２の反射補正値算出部４３は、パラメータ記憶部２により記憶されている波長特性データを取得する。
　反射補正値算出部４３は、取得した波長特性データを用いて、反射輝度算出部４１から出力されたバンドｂ毎の各々の小領域Ｓ（ｉ，ｊ）における大気反射放射輝度Ｌｒ（ｂ，ｉ，ｊ）の重み係数ε”（ｂ）を算出する。重み係数ε”（ｂ）は、例えば、重み係数ε（ｂ）又は重み係数ε’（ｂ）と同様に算出することができる。
　反射補正値算出部４３は、各々の小領域Ｓ（ｉ，ｊ）における分光放射輝度の第３の補正値Ｌｒ（Ｐａ’，ｉ，ｊ）として、以下の式（１７）に示すように、重み係数ε”（ｂ）を用いて、バンドｂ毎の各々の小領域Ｓ（ｉ，ｊ）における大気反射放射輝度Ｌｒ（ｂ，ｉ，ｊ）の総和を算出する。
Figure JPOXMLDOC01-appb-I000011
ｂ＝ｂ１，ｂ２，ｂ３，ｂ４
　反射補正値算出部４３は、各々の小領域Ｓ（ｉ，ｊ）における分光放射輝度の第３の補正値Ｌｒ（Ｐａ’，ｉ，ｊ）を画像再結合部４４に出力する。
　補正値算出部４２の画像再結合部４４は、上記実施の形態１の画像再結合部１０と同様に、第１の補正値の２次元配列Ａ及び第２の補正値の２次元配列Ｂのそれぞれを生成する。
　また、画像再結合部４４は、反射補正値算出部４３から出力された各々の小領域Ｓ（ｉ，ｊ）における分光放射輝度の第３の補正値Ｌｒ（Ｐａ’，ｉ，ｊ）を結合して、第３の補正値の２次元配列Ｃを生成する。
　２次元配列Ｃにおける各々の要素（ｉ，ｊ）は、画像分割部１により分割された各々の小領域Ｓ（ｉ，ｊ）に対応している。マルチスペクトル画像の分割数がＮで、小領域Ｓ（ｉ，ｊ）の個数がＮであれば、２次元配列Ｃの要素（ｉ，ｊ）の数もＮである。
　２次元配列Ｃにおける各々の要素（ｉ，ｊ）は、各々の小領域Ｓ（ｉ，ｊ）における分光放射輝度の第３の補正値Ｌｒ（Ｐａ’，ｉ，ｊ）を有している。
　画像再結合部４４は、生成した２次元配列Ａ，Ｂ，Ｃを補正値補間部４５に出力する。
　補正値算出部４２の補正値補間部４５は、上記実施の形態１の補正値補間部１１と同様に、画像再結合部４４から出力された２次元配列Ａ，Ｂのそれぞれを補間する。
　また、補正値補間部４５は、画像再結合部４４から出力された２次元配列Ｃを補間する。
　例えば、パンクロマチック画像のｘ方向の画素数が、２次元配列Ｃにおけるｘ方向の要素（ｉ，ｊ）の数のＭｘ倍であれば、２次元配列Ｃのｘ方向をＭｘ倍にアップサンプリングする。
　これにより、第３の補正値Ｌｒ（Ｐａ’，ｉ，ｊ）は、第３の補正値Ｌｒ（Ｐａ’，ｘ，ｙ）に変換され、変換後の第３の補正値Ｌｒ（Ｐａ’，ｘ，ｙ）が補正値補間部４５から反射補正部４７に出力される。
　補正部４６の反射補正部４７は、以下の式（１８）に示すように、放射輝度補正部１３から出力されたパンクロマチック画像の補正後の分光放射輝度Ｌ１（ｘ，ｙ）から、補正値補間部４５から出力された第３の補正値Ｌｒ（Ｐａ’，ｘ，ｙ）を減算することで、パンクロマチック画像における分光放射輝度を補正する。
Ｌ３（ｘ，ｙ）＝Ｌ１（ｘ，ｙ）－Ｌｒ（Ｐａ’，ｘ，ｙ）　　　（１８）
　式（１８）において、Ｌ３（ｘ，ｙ）は、補正後の分光放射輝度である。
　反射補正部４７は、分光放射輝度を補正したパンクロマチック画像を透過率補正部１４に出力する。
　この実施の形態２では、反射補正部４７が放射輝度補正部１３の後段に設けられている例を示しているが、反射補正部４７が放射輝度補正部１３の前段に設けられていてもよい。
　以上で明らかなように、この実施の形態２によれば、補正値算出部４２が、反射輝度算出部４１により波長帯毎に算出された各々の領域における大気反射放射輝度から、パンクロマチック画像における分光放射輝度の第３の補正値を算出し、補正部４６が、補正値算出部４２により算出された第１及び第２の補正値を用いて、パンクロマチック画像における分光放射輝度を補正するほかに、補正値算出部４２により算出された第３の補正値を用いて、パンクロマチック画像における分光放射輝度を補正するように構成したので、上記実施の形態１，２と同様の効果を奏するほかに、マルチスペクトル画像及びパンクロマチック画像に重畳されている大気からの反射光成分を除去することができる効果を奏する。
実施の形態３．
　この実施の形態３では、マルチスペクトル画像における分光放射輝度についても補正し、分光放射輝度が補正されたパンクロマチック画像と、分光放射輝度が補正されたマルチスペクトル画像との色合成を行う例を説明する。
　図１４は、この発明の実施の形態３による画像処理装置を示す構成図であり、図１５は、この発明の実施の形態３による画像処理装置を示すハードウェア構成図である。
　図１４及び図１５において、図１１及び図１２と同一符号は同一または相当部分を示すので説明を省略する。
　この実施の形態３では、補正値算出部４２は、上記実施の形態２と同様に、パンクロマチック画像における分光放射輝度の補正値を算出するほかに、マルチスペクトル画像における分光放射輝度の補正値を算出する処理を実施する。
　この実施の形態３では、補正部４６は、上記実施の形態２と同様に、パンクロマチック画像における分光放射輝度を補正するほかに、マルチスペクトル画像における分光放射輝度を補正する処理を実施する。
　色合成処理部４８は、図１５に示す色合成処理回路５４で実現される。
　色合成処理部４８は、補正部４６により分光放射輝度が補正されたパンクロマチック画像と補正部４６により分光放射輝度が補正されたマルチスペクトル画像との色合成処理を実施する。
　図１４では、色合成処理部４８が、上記実施の形態２における図１１の画像処理装置に適用されている例を示しているが、色合成処理部４８が、上記実施の形態１における図１の画像処理装置に適用されているものであってもよい。
　図１４では、画像処理装置の構成要素である画像分割部１、パラメータ記憶部２、散乱輝度算出部３、大気透過率算出部６、反射輝度算出部４１、補正値算出部４２、補正部４６及び色合成処理部４８のそれぞれが、図１５に示すような専用のハードウェアで実現されるものを想定している。即ち、画像分割回路２１、パラメータ記憶回路２２、散乱輝度算出回路２３、透過率算出回路２４、反射輝度算出回路５１、補正値算出回路５２、補正回路５３及び色合成処理回路５４で実現されるものを想定している。
　画像分割回路２１、散乱輝度算出回路２３、透過率算出回路２４、反射輝度算出回路５１、補正値算出回路５２、補正回路５３及び色合成処理回路５４は、例えば、単一回路、複合回路、プログラム化したプロセッサ、並列プログラム化したプロセッサ、ＡＳＩＣ、ＦＰＧＡ、または、これらを組み合わせたものが該当する。
　画像処理装置の構成要素は、専用のハードウェアで実現されるものに限るものではなく、画像処理装置がソフトウェア、ファームウェア、または、ソフトウェアとファームウェアとの組み合わせで実現されるものであってもよい。
　画像処理装置がソフトウェア又はファームウェアなどで実現される場合、パラメータ記憶部２を図３に示すコンピュータのメモリ３１上に構成するとともに、画像分割部１、散乱輝度算出部３、大気透過率算出部６、反射輝度算出部４１、補正値算出部４２、補正部４６及び色合成処理部４８の処理手順をコンピュータに実行させるためのプログラムをメモリ３１に格納し、コンピュータのプロセッサ３２がメモリ３１に格納されているプログラムを実行するようにすればよい。
　次に動作について説明する。
　以下、上記実施の形態１，２と相違する箇所を説明する。
　この実施の形態３では、輝度補正値算出部８は、各々の小領域Ｓ（ｉ，ｊ）における分光放射輝度の第１の補正値Ｌｓｃａｔ（Ｐａ’，ｉ，ｊ）を画像再結合部４４に出力するほか、各々の小領域Ｓ（ｉ，ｊ）における大気散乱放射輝度Ｌｓｃａｔ（ｂ，ｉ，ｊ）を第４の補正値として画像再結合部４４に出力する。
　この実施の形態３では、透過率補正値算出部９は、各々の小領域Ｓ（ｉ，ｊ）における分光放射輝度の第２の補正値τ（Ｐａ’，ｉ，ｊ）を画像再結合部４４に出力するほか、各々の小領域Ｓ（ｉ，ｊ）における大気透過率τ（ｂ，ｉ，ｊ）を第５の補正値として画像再結合部４４に出力する。
　この実施の形態３では、反射補正値算出部４３は、各々の小領域Ｓ（ｉ，ｊ）における分光放射輝度の第３の補正値Ｌｒ（Ｐａ’，ｉ，ｊ）として画像再結合部４４に出力するほか、各々の小領域Ｓ（ｉ，ｊ）における大気反射放射輝度Ｌｒ（ｂ，ｉ，ｊ）を第６の補正値として画像再結合部４４に出力する。
　補正値算出部４２の画像再結合部４４は、上記実施の形態２と同様に、第１の補正値の２次元配列Ａを生成するほかに、輝度補正値算出部８から出力された各々の小領域Ｓ（ｉ，ｊ）における分光放射輝度の第４の補正値であるＬｓｃａｔ（ｂ，ｉ，ｊ）を結合して、第４の補正値の２次元配列Ｄを生成する。
　また、画像再結合部４４は、第２の補正値の２次元配列Ｂを生成するほかに、透過率補正値算出部９から出力された各々の小領域Ｓ（ｉ，ｊ）における分光放射輝度の第５の補正値である大気透過率τ（ｂ，ｉ，ｊ）を結合して、第５の補正値の２次元配列Ｅを生成する。
　また、画像再結合部４４は、第３の補正値の２次元配列Ｃを生成するほかに、反射補正値算出部４３から出力された各々の小領域Ｓ（ｉ，ｊ）における分光放射輝度の第６の補正値である大気反射放射輝度Ｌｒ（ｂ，ｉ，ｊ）を結合して、第６の補正値の２次元配列Ｆを生成する。
　補正値算出部４２の補正値補間部４５は、上記実施の形態２と同様に、画像再結合部４４から出力された２次元配列Ａ，Ｂ，Ｃのそれぞれを補間する。
　また、補正値補間部４５は、画像再結合部４４から出力された２次元配列Ｄ，Ｅ，Ｆのそれぞれを補間する。
　これにより、第４の補正値であるＬｓｃａｔ（ｂ，ｉ，ｊ）は、Ｌｓｃａｔ（ｂ，ｘ，ｙ）に変換され、変換後の第４の補正値であるＬｓｃａｔ（ｂ，ｘ，ｙ）が補正値補間部４５から放射輝度補正部１３に出力される。
　第５の補正値である大気透過率τ（ｂ，ｉ，ｊ）は、τ（ｂ，ｘ，ｙ）に変換され、変換後の第５の補正値であるτ（ｂ，ｘ，ｙ）が補正値補間部４５から透過率補正部１４に出力される。
　第６の補正値である大気反射放射輝度Ｌｒ（ｂ，ｉ，ｊ）は、Ｌｒ（ｂ，ｘ，ｙ）に変換され、変換後の第６の補正値であるＬｒ（ｂ，ｘ，ｙ）が補正値補間部４５から反射補正部４７に出力される。
　放射輝度補正部１３は、上記実施の形態１，２と同様に、パンクロマチック画像における分光放射輝度を補正する。
　また、放射輝度補正部１３は、マルチスペクトル画像のバンドｂ毎の各々の小領域Ｓ（ｉ，ｊ）における分光放射輝度Ｌｓｅｎｓｏｒ（ｂ，ｉ，ｊ）から、補正値補間部４５から出力された第４の補正値であるＬｓｃａｔ（ｂ，ｘ，ｙ）を減算することで、マルチスペクトル画像における分光放射輝度を補正する。
　反射補正部４７は、上記実施の形態２と同様に、パンクロマチック画像における分光放射輝度を補正する。
　また、反射補正部４７は、放射輝度補正部１３から出力されたマルチスペクトル画像の補正後の各々の小領域Ｓ（ｉ，ｊ）における分光放射輝度Ｌｓｅｎｓｏｒ（ｂ，ｉ，ｊ）から、補正値補間部４５から出力された第６の補正値であるＬｒ（ｂ，ｘ，ｙ）を減算することで、マルチスペクトル画像における分光放射輝度を補正する。
　透過率補正部１４は、上記実施の形態１，２と同様に、パンクロマチック画像における分光放射輝度を補正する。
　また、透過率補正部１４は、反射補正部４７から出力されたマルチスペクトル画像の補正後の各々の小領域Ｓ（ｉ，ｊ）における分光放射輝度Ｌｓｅｎｓｏｒ（ｂ，ｉ，ｊ）を、補正値補間部４５から出力された第５の補正値であるτ（ｂ，ｘ，ｙ）で除算することで、マルチスペクトル画像における分光放射輝度を補正する。
　色合成処理部４８は、補正部４６により分光放射輝度が補正されたパンクロマチック画像と、補正部４６により分光放射輝度が補正されたマルチスペクトル画像との色合成処理を実施し、色合成処理後の画像を出力する。
　色合成処理として、例えば、ＩＨＳ（Ｉｎｔｅｎｓｉｔｙ　Ｈｕｅ　Ｓａｔｕｒａｔｉｏｎ）変換と呼ばれる色合成処理を用いることができる。ＩＨＳ変換と呼ばれる色合成処理では、マルチスペクトル画像から色空間である輝度、色相及び彩度を算出し、この算出した輝度をパンクロマチック画像に置き換えることで、パンクロマチック画像の分解能を持つ輝度情報に、色相及び彩度の情報を加えるものである。
　色合成処理としては、ｂｒｏｖｅｙ変換、Ｇｒａｍ－Ｓｃｈｍｉｄｔ変換、あるいは、主成分分析を用いた変換処理なども使用することができる。
　高解像度のパンクロマチック画像と低解像度のマルチスペクトル画像とを合成し、高解像度かつ色情報を保有する画像は、パンシャープン画像と呼ばれることがある。
　なお、本願発明はその発明の範囲内において、各実施の形態の自由な組み合わせ、あるいは各実施の形態の任意の構成要素の変形、もしくは各実施の形態において任意の構成要素の省略が可能である。
　この発明は、パンクロマチック画像における分光放射輝度を補正する画像処理装置に適している。
　１　画像分割部、２　パラメータ記憶部、３　散乱輝度算出部、４　暗部輝度算出部、５　散乱輝度算出処理部、６　大気透過率算出部、７　補正値算出部、８　輝度補正値算出部、９　透過率補正値算出部、１０　画像再結合部、１１　補正値補間部、１２　補正部、１３　放射輝度補正部、１４　透過率補正部、２１　画像分割回路、２２　パラメータ記憶回路、２３　散乱輝度算出回路、２４　透過率算出回路、２５　補正値算出回路、２６　補正回路、３１　メモリ、３２　プロセッサ、４１　反射輝度算出部、４２　補正値算出部、４３　反射補正値算出部、４４　画像再結合部、４５　補正値補間部、４６　補正部、４７　反射補正部、４８　色合成処理部、５１　反射輝度算出回路、５２　補正値算出回路、５３　補正回路、５４　色合成処理回路。

Claims (4)

  1. 　複数の波長帯の画像を含んでいるマルチスペクトル画像を複数の領域に分割する画像分割部と、
    　前記マルチスペクトル画像が含んでいる複数の画像の波長帯毎に、前記画像分割部により分割された各々の領域における大気散乱を示す大気散乱放射輝度をそれぞれ算出する散乱輝度算出部と、
    　前記散乱輝度算出部により波長帯毎に算出された各々の領域における大気散乱放射輝度から、前記複数の波長帯を含む単一波長帯の画像であるパンクロマチック画像における分光放射輝度の補正値を算出する補正値算出部と、
    　前記補正値算出部により算出された分光放射輝度の補正値を用いて、前記パンクロマチック画像における分光放射輝度を補正する補正部と
    　を備えた画像処理装置。
  2. 　前記散乱輝度算出部により波長帯毎に算出された各々の領域における大気散乱放射輝度から、前記波長帯毎に、前記画像分割部により分割された各々の領域における大気散乱の大気透過率をそれぞれ算出する大気透過率算出部を備え、
    　前記補正値算出部は、前記パンクロマチック画像における分光放射輝度の補正値として、前記散乱輝度算出部により波長帯毎に算出された各々の領域における大気散乱放射輝度から、前記パンクロマチック画像における分光放射輝度の第１の補正値を算出するほかに、前記大気透過率算出部により波長帯毎に算出された各々の領域における大気散乱の大気透過率から、前記パンクロマチック画像における分光放射輝度の第２の補正値を算出し、
    　前記補正部は、前記補正値算出部により算出された第１の補正値を用いて、前記パンクロマチック画像における分光放射輝度を補正するほかに、前記補正値算出部により算出された第２の補正値を用いて、前記パンクロマチック画像における分光放射輝度を補正することを特徴とする請求項１記載の画像処理装置。
  3. 　前記マルチスペクトル画像の波長帯毎の各々の領域における分光放射輝度と、前記散乱輝度算出部により波長帯毎に算出された各々の領域における大気散乱放射輝度とから、前記波長帯毎に、前記画像分割部により分割された各々の領域における大気反射を示す大気反射放射輝度をそれぞれ算出する反射輝度算出部を備え、
    　前記補正値算出部は、前記第１及び第２の補正値を算出するほかに、前記反射輝度算出部により波長帯毎に算出された各々の領域における大気反射放射輝度から、前記パンクロマチック画像における分光放射輝度の第３の補正値を算出し、
    　前記補正部は、前記補正値算出部により算出された第１及び第２の補正値を用いて、前記パンクロマチック画像における分光放射輝度を補正するほかに、前記補正値算出部により算出された第３の補正値を用いて、前記パンクロマチック画像における分光放射輝度を補正することを特徴とする請求項２記載の画像処理装置。
  4. 　前記補正値算出部は、前記散乱輝度算出部により波長帯毎に算出された各々の領域における大気散乱放射輝度から、前記波長帯毎に、前記マルチスペクトル画像における分光放射輝度の補正値を算出し、
    　前記補正部は、前記補正値算出部により算出されたマルチスペクトル画像における分光放射輝度の補正値を用いて、前記マルチスペクトル画像における分光放射輝度を補正し、
    　前記補正部により分光放射輝度が補正されたパンクロマチック画像と前記補正部により分光放射輝度が補正されたマルチスペクトル画像との色合成を行う色合成処理部を備えたことを特徴とする請求項１記載の画像処理装置。
PCT/JP2017/015585 2017-04-18 2017-04-18 画像処理装置 WO2018193517A1 (ja)

Priority Applications (2)

Application Number Priority Date Filing Date Title
PCT/JP2017/015585 WO2018193517A1 (ja) 2017-04-18 2017-04-18 画像処理装置
JP2019512336A JP6556409B2 (ja) 2017-04-18 2017-04-18 画像処理装置

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
PCT/JP2017/015585 WO2018193517A1 (ja) 2017-04-18 2017-04-18 画像処理装置

Publications (1)

Publication Number Publication Date
WO2018193517A1 true WO2018193517A1 (ja) 2018-10-25

Family

ID=63855659

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2017/015585 WO2018193517A1 (ja) 2017-04-18 2017-04-18 画像処理装置

Country Status (2)

Country Link
JP (1) JP6556409B2 (ja)
WO (1) WO2018193517A1 (ja)

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPH06300845A (ja) * 1993-04-16 1994-10-28 N T T Data Tsushin Kk 画像情報処理装置
US20140270332A1 (en) * 2013-03-15 2014-09-18 Digitalglobe, Inc. Atmospheric compensation in satellite imagery
JP2015215783A (ja) * 2014-05-12 2015-12-03 株式会社日立製作所 画像処理装置、画像処理方法及びプログラムを記録した記録媒体
JP2016126566A (ja) * 2015-01-05 2016-07-11 三菱電機株式会社 画像処理装置及び画像処理方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPH06300845A (ja) * 1993-04-16 1994-10-28 N T T Data Tsushin Kk 画像情報処理装置
US20140270332A1 (en) * 2013-03-15 2014-09-18 Digitalglobe, Inc. Atmospheric compensation in satellite imagery
JP2015215783A (ja) * 2014-05-12 2015-12-03 株式会社日立製作所 画像処理装置、画像処理方法及びプログラムを記録した記録媒体
JP2016126566A (ja) * 2015-01-05 2016-07-11 三菱電機株式会社 画像処理装置及び画像処理方法

Also Published As

Publication number Publication date
JPWO2018193517A1 (ja) 2019-06-27
JP6556409B2 (ja) 2019-08-07

Similar Documents

Publication Publication Date Title
US7826685B2 (en) Spatial and spectral calibration of a panchromatic, multispectral image pair
US8094960B2 (en) Spectral calibration of image pairs using atmospheric characterization
US9342760B1 (en) System and method for combining color information with spatial information in multispectral images
US20100008595A1 (en) Automated atmospheric characterization of remotely sensed multi-spectral imagery
JP6964834B2 (ja) 画像処理装置および画像処理方法
EP2095330A1 (en) Panchromatic modulation of multispectral imagery
CN102341827A (zh) 利用大气和太阳照射建模来融合图像对的系统和方法
US11039076B2 (en) Information processing apparatus, information processing method, and storage medium
CN109643440B (zh) 图像处理设备、图像处理方法和计算机可读记录介质
CN109074637B (zh) 用于从多个相应的输入图像通道生成输出图像的方法和系统
JP6869652B2 (ja) 画像処理装置、撮像装置、画像処理方法、画像処理プログラム、および、記憶媒体
JP5921311B2 (ja) 画像処理装置及び画像処理方法
JP6556409B2 (ja) 画像処理装置
JP2016126566A (ja) 画像処理装置及び画像処理方法
JP6463244B2 (ja) 画像処理装置及び画像処理方法
JP2017198464A (ja) 画像処理装置及び画像処理方法
KR101630856B1 (ko) 다분광 포토메트릭 스테레오 시스템 그리고 이의 동작 방법
KR101952394B1 (ko) 레티넥스 모델 기반 엘이디영상 색 보정 방법
Kaur et al. Enriched image demosaicing using illuminate normalization and content based color filter array
Anshu et al. Evaluation of Fusion Techniques for High Resolution Data-A Worldview-2 Imagery
Jacobsen Effective resolution of digital frame images
Gao et al. Spectral Reconstruction Based on Svm for Cross Calibration
JP2009188806A (ja) 撮像方法および撮像システム
Borel et al. Improved panchromatic sharpening of multi-spectral image data

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 17906436

Country of ref document: EP

Kind code of ref document: A1

ENP Entry into the national phase

Ref document number: 2019512336

Country of ref document: JP

Kind code of ref document: A

NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 17906436

Country of ref document: EP

Kind code of ref document: A1