利用冥王星影像创建镶嵌图
原文链接: https://www.nv5geospatialsoftware.com/Learn/Blogs/Blog-Details/creating-mosaics-from-pluto-imagery
18800 评价这篇文章:
暂无评分
利用冥王星影像创建镶嵌图
Zachary Norman 2016年1月25日，星期一
在我生命的大部分时间里，我一直是个太空迷，总是喜欢观看太阳系中其他行星、星云和其他星系的图像。最近，被拍摄的天体家族迎来了一位新成员：冥王星。2015年7月14日，新视野号探测器飞抵距冥王星表面12，600公里处，拍摄了一些令人惊叹的图像。有了这些新图像，我认为尝试使用 ENVI 将其中几张图像拼接成镶嵌图会是一个很棒的想法，看看是否能实现。
为了获取适合镶嵌的图像，我需要确保它们是连续快速拍摄的，以便彼此之间有足够的重叠区域。这样可以让 ENVI 能够相对轻松地将图像相互配准。为了找到这些图像，我从新视野号探测器的 LORRI 仪器获取了数据，链接如下：
https://pluto.jhuapl.edu/soc/Pluto-Encounter/
在选择要镶嵌的图像时，非常重要的是图像需要有相似的采集时间。在浏览了 LORRI 提供的可用图像后，我找到了13张图像，采集时间从UTC时间11:36:00到UTC时间11:36:36。能够获得冥王星的图像固然很好，但遗憾的是，这些数据没有附带元数据来精确地进行地理定位。为了解决元数据缺失的问题，我决定围绕一个参考图像创建一个虚拟的空间参考。我选择使用UTC时间11:36:00拍摄的图像作为基础图像，并为该图像赋予了一个基本的 WGS-84 坐标系，假设像素大小为1米 x 1米。此时需要指出的是，我虚拟的坐标系和像素大小很可能与实际的情况完全不符，我只是需要创建一个空间参考，以便用 ENVI 来配准和镶嵌图像。
有了参考图像和虚拟空间参考后，下一步是弄清楚如何使图像彼此之间正确配准。在对图像配准工作流进行一番研究后，我意识到为了获得最佳结果，需要让图像尽可能接近它们实际的相对位置。为此，我查看了我的参考图像和系列中待配准的第一张图像（参考图像的采集时间为2015-07-14 11:36:00 UTC，待配准的第一张图像采集于2015-07-14 11:36:03 UTC）。
在仔细查看这两张图像并手动调整后，我发现第二张图像相对于第一张图像的近似像素偏移量在x和y方向上约为[-580, -18]像素。这意味着，基础图像之后系列中的第一张图像的偏移量为[-580, -18]像素，第二张图像的偏移量为[-1160, -36]像素。此时，我决定使用 ENVI API 以编程方式将此偏移量应用到系列中的每张图像。选择使用 ENVI API 的一个重要原因是，我不必为系列中的每张图像逐步执行图像配准工作流，并且，如果我更改了偏移量，也能快速应用到不同的图像上。
根据像素偏移量对每张图像进行重投影后，图像此时已经相对接近了。剩下的就是应用色彩平衡（通过自定义的直方图匹配功能）、将图像相互配准，然后镶嵌在一起。以下是我使用 ENVI API 执行的步骤的简要总结：
对参考图像应用色彩匹配 -> 使用任务 GenerateTiePointsByCrossCorrelation -> 使用任务 FilterTiePointsByGlobalTransform -> 使用任务 ImageToImageRegistration -> 使用任务 BuildMosaicRaster。
下方是输出的镶嵌图以及用于创建该图像的 IDL 代码。
;+
; :Author: Zach Norman
;
; :Description:
; 一个注册和镶嵌来自新视野号 LORRI 仪器部分冥王星图像的小型示例程序。将此代码另存为 pluto\_mosaicking.pro，要运行可以在 IDL Workbench 中按 'Run'。
; 请参阅下面的输出部分以了解查找程序结果的位置以及结果的详细信息。
;
; 代码编写为自动将冥王星图像下载到名为 images 的目录中，该目录将与此 .pro 文件所在的同一位置创建。如果下载文件有任何问题，
; 则可能需要手动从网上下载。如果来自 NASA 的图像的 URL 发生更改，这可能是有必要的。如果手动下载图像，也可能需要做一些额外的更改。
;
; 具体来说，您可能必须更改程序 pluto\_mosaicking 中的一些参数，使其与您的计算机当前状态匹配。如果您从
; http://pluto.jhuapl.edu/soc/Pluto-Encounter/index.php 下载冥王星图像，则需要指定
; pluto\_mosaicking 过程中 "imagedir" 的值。
;
; 如果您想或不得不手动下载图像，共有13张图像，其采集时间范围从2015-07-14 11:36:00 UTC 到 2015-07-14 11:36:36 UTC。目前，
; 代码的编写方式是让您将此 .pro 文件放在包含图像的目录旁边。
;
;
; :Output:
; 将在此 .pro 代码所在的同一目录中生成三个文件，它们将被命名为：pluto\_mosaic, pluto\_mosaic.hdr, pluto\_mosaic.tif。
;
; 第一个文件是 ENVI 格式的图像，可以使用头文件 (.hdr) 读入 ENVI，第三个文件是 TIFF 图像，因此您可以使用通用照片查看器查看镶嵌图。
; 请注意，ENVI 头文件包含虚构的投影信息。由于 LORRI 图像未提供元数据，我为每张图像创建了一个虚拟的空间参考。
;
;-
;用于图像间直方图匹配的小函数
;这很有用，因为一些图像来自 NASA 的自动缩放效果不佳，它们看起来非常暗
**function** **histogram\_matching**, match, base
**compile\_opt** idl2
maxval = (**max**(match))>(**max**(base))
;计算直方图和累积密度
match\_histogram = **Histogram**(match, Min=**0**, Max=maxval, Binsize=**1**)
match\_cdf = **TOTAL**(match\_histogram, /CUMULATIVE) / **N\_ELEMENTS**(match\_histogram)
base\_histogram = **Histogram**(base, Min=**0**, Max=maxval, Binsize=**1**)
base\_cdf = **TOTAL**(base\_histogram, /CUMULATIVE) / **N\_ELEMENTS**(base\_histogram)
;创建变换函数 z 来转换匹配图像中的数据空间
z = **bytscl**(base\_histogram)
**FOR** j=**0**,**n\_elements**(z)-**1** **DO** **BEGIN**
i = **Where**(match\_cdf **LT** base\_cdf[j], count)
**IF** count **GT** **0** **THEN** z[j] = i[-**1**] **ELSE** z[j]=**0**
**ENDFOR**
;执行直方图匹配
matched = z[base]
**return**, matched
**end**
;用于 ENVI 关闭所有数据源的方法
**pro** **ENVI**::reset
**compile\_opt** idl2
e = **envi**(/current)
**if** (e **eq** !NULL) **then** **return**
**if** (e.*WIDGET\_ID* **ne** **0**) **then** **begin**
;重置视图
views = e.**getview**(/all)
**foreach** view,views **do** view.*close*
**endif**
;关闭所有打开的栅格数据
opendata = (e.*data*).**Get**(/raster,/vector, /series)
**foreach** raster, opendata **do** raster.*close*
**end**
**pro** **pluto\_mosaicking**
**compile\_opt** idl2
;启动 ENVI 并使其处于新状态！
e = **envi**(/current)
**if** (e **eq** !NULL) **then** e = **envi**(/headless) **else** e.*reset*
;每张图像的近似像素偏移量
px\_off = [-**580**, -**18**]
;运行代码的目录
thisdir = **FILE\_DIRNAME**( (((**scope\_traceback**(/STRUCTURE))[-**1**]).*FILENAME*))
;imagedir 是原始图像的位置
imagedir = thisdir + **path\_sep**() + 'images'
;copydir 是来自 imagedir 的图像副本的位置，这些副本具有虚拟的空间参考
copydir = thisdir + **path\_sep**() + 'image\_copies'
;warpeddir 是配准后的图像存放位置
warpeddir = thisdir + **path\_sep**() + 'images\_warped'
;如果 imagedir 不存在，则创建它并将文件下载到该目录
;如果 URL 发生变化，则需要手动修改这些 URL 或从以下 URL 手动下载：
; http://pluto.jhuapl.edu/soc/Pluto-Encounter/index.php
**if** ~**file\_test**(imagedir) **then** **begin**
**print**, '没有找到目录 ' + imagedir
**print**, **string**(**9b**) + '正在创建目录并下载图像...'
**FILE\_MKDIR**, imagedir
urls = ['http://pluto.jhuapl.edu/soc/Pluto-Encounter/data/pluto/level2/lor/jpeg/029917/lor\_0299179679\_0x636\_sci\_3.jpg',$
'http://pluto.jhuapl.edu/soc/Pluto-Encounter/data/pluto/level2/lor/jpeg/029917/lor\_0299179682\_0x636\_sci\_3.jpg',$
'http://pluto.jhuapl.edu/soc/Pluto-Encounter/data/pluto/level2/lor/jpeg/029917/lor\_0299179685\_0x636\_sci\_2.jpg',$
'http://pluto.jhuapl.edu/soc/Pluto-Encounter/data/pluto/level2/lor/jpeg/029917/lor\_0299179688\_0x636\_sci\_4.jpg',$
'http://pluto.jhuapl.edu/soc/Pluto-Encounter/data/pluto/level2/lor/jpeg/029917/lor\_0299179691\_0x636\_sci\_4.jpg',$
'http://pluto.jhuapl.edu/soc/Pluto-Encounter/data/pluto/level2/lor/jpeg/029917/lor\_0299179694\_0x636\_sci\_4.jpg',$
'http://pluto.jhuapl.edu/soc/Pluto-Encounter/data/pluto/level2/lor/jpeg/029917/lor\_0299179697\_0x636\_sci\_4.jpg',$
'http://pluto.jhuapl.edu/soc/Pluto-Encounter/data/pluto/level2/lor/jpeg/029917/lor\_0299179700\_0x636\_sci\_4.jpg',$
'http://pluto.jhuapl.edu/soc/Pluto-Encounter/data/pluto/level2/lor/jpeg/029917/lor\_0299179703\_0x636\_sci\_4.jpg',$
'http://pluto.jhuapl.edu/soc/Pluto-Encounter/data/pluto/level2/lor/jpeg/029917/lor\_0299179706\_0x636\_sci\_4.jpg',$
'http://pluto.jhuapl.edu/soc/Pluto-Encounter/data/pluto/level2/lor/jpeg/029917/lor\_0299179709\_0x636\_sci\_3.jpg',$
'http://pluto.jhuapl.edu/soc/Pluto-Encounter/data/pluto/level2/lor/jpeg/029917/lor\_0299179712\_0x636\_sci\_3.jpg',$
'http://pluto.jhuapl.edu/soc/Pluto-Encounter/data/pluto/level2/lor/jpeg/029917/lor\_0299179715\_0x636\_sci\_3.jpg']
;使用 IDLnetURL 对象下载每张图像
urlobj = **IDLnetURL**()
**for** i=**0**, **n\_elements**(urls)-**1** **do** **begin**
outfile = imagedir + **path\_sep**() + **strmid**(urls[i],**strpos**(urls[i],'/',/REVERSE\_SEARCH)+**1**)
void = urlobj.**Get**(FILENAME = outfile,$
URL = urls[i])
**print**, **string**(**9b**) + **string**(**9b**) + '已下载图像 ' + **strtrim**(i+**1**,**2**) + ' of ' + **strtrim**(**n\_elements**(urls),**2**)
**endfor**
urlobj.**cleanup**
**print**, string (**9b**) + '已下载所有图像！'
**endif**
;如果输出目录不存在，则创建它们
**if** ~**file\_test**(copydir) **then** **FILE\_MKDIR**, copydir
**if** ~**file\_test**(warpeddir) **then** **FILE\_MKDIR**, warpeddir
;转到 imagedir 并查找冥王星图像
**cd**, imagedir, CURRENT = first\_dir
files = **file\_search**('\*.jpg')
;检查我们是否已经制作了图像的副本
;创建一个虚拟的空间参考
psize = **1d**\*(**180**/!CONST.*R\_EARTH*);1米的像素大小
spatialRef = **ENVIStandardRasterSpatialRef**( $
COORD\_SYS\_CODE=**4326**, /GEOGCS, $
PIXEL\_SIZE=[psize,psize], TIE\_POINT\_PIXEL=[**0.0D**,**0.0D**], $
TIE\_POINT\_MAP=[**0D**,**0D**])
;确保对文件进行排序，以便按时间排列
files = files[**sort**(files)]
;复制所有文件并创建具有基础空间参考的 ENVI 栅格
reference = e.**openraster**(files[**0**], SPATIALREF\_OVERRIDE=spatialref)
reference\_data = reference.**GetData**()
;将原始栅格导出到副本栅格
outfile = copydir + **path\_sep**() + **file\_basename**(files[**0**], '.jpg')
**if** **file\_test**(outfile) **then** **FILE\_DELETE**, outfile
reference.**export**, outfile, 'ENVI'
;对所有图像执行直方图匹配并将它们导出到 copydir
**for** i=**1**, **n\_elements**(files)-**1** **do** **begin**
spatialRef = **ENVIStandardRasterSpatialRef**( $
COORD\_SYS\_CODE=**4326**, /GEOGCS, $
PIXEL\_SIZE=[psize,psize], TIE\_POINT\_PIXEL=[**0.0D**,**0.0D**], $
TIE\_POINT\_MAP=i\*psize\*px\_off)
warp = e.**openraster**(files[i], SPATIALREF\_OVERRIDE=spatialRef)
warp\_data = warp.**GetData**()
warp\_match = **histogram\_matching**(reference\_data, warp\_data)
outfile = copydir + **path\_sep**() + **file\_basename**(files[i], '.jpg')
**if** **file\_test**(outfile) **then** **FILE\_DELETE**, outfile
warp\_copy = **ENVIRaster**(warp\_data, INHERITS\_FROM=warp, URI = outfile)
warp.**close**
warp\_copy.**save**
warp\_copy.**close**
**endfor**
reference.**close**
;转到复制目录并将图像彼此配准
**cd**, copydir
files = **file\_basename**(files, '.jpg')
reference = e.**openraster**(files[**0**])
;将参考图像复制到副本目录并在 ENVI 中打开它
outfile = warpeddir + **path\_sep**() + files[**0**]
**if** **file\_test**(outfile) **then** **FILE\_DELETE**, outfile
reference.**export**,outfile, 'ENVI'
reference.**close**
reference = e.**openraster**(outfile)
;对所有图像执行图像配准
**print**, '发现 ' + **strtrim**(files.*length*,**2**) + ' 张需要配准的图像！'
**for** i=**1**, **n\_elements**(files)-**1** **do** **begin**
warp\_this = e.**openraster**(files[i])
; 从 ENVITasks 目录中获取自动连接点生成任务
TieTask = **ENVITask**('GenerateTiePointsByCrossCorrelation')
TieTask.*INPUT\_RASTER1* = reference
TieTask.*INPUT\_RASTER2* = warp\_this
TieTask.**Execute**
TieTask\_TiePoints = TieTask.*OUTPUT\_TIEPOINTS*
; 从 ENVITasks 目录中获取连接点过滤任务
FilterTask = **ENVITask**('FilterTiePointsByGlobalTransform')
FilterTask.*INPUT\_TIEPOINTS* = TieTask\_TiePoints
FilterTask.**Execute**
FilterTask\_TiePoints = FilterTask.*OUTPUT\_TIEPOINTS*
; 从 ENVITasks 目录中获取图像到图像配准任务
RegistrationTask = **ENVITask**('ImageToImageRegistration')
RegistrationTask.*INPUT\_TIEPOINTS* = FilterTask\_TiePoints
RegistrationTask.*WARPING* = 'RST'
outfile = warpeddir + **path\_sep**() + files[i]
**if** **file\_test**(outfile) **then** **FILE\_DELETE**, outfile
RegistrationTask.*OUTPUT\_RASTER\_URI* = warpeddir + **path\_sep**() + files[i]
RegistrationTask.**Execute**
Reference.**close**
Reference = RegistrationTask.*OUTPUT\_RASTER*
;关闭已配准的栅格
warp\_this.**close**
;删除连接点文件
**FILE\_DELETE**, TieTask.*OUTPUT\_TIEPOINTS\_URI*
**FILE\_DELETE**, FilterTask.*OUTPUT\_TIEPOINTS\_URI*
**print**, **string**(**9b**) + '已完成 ' + **strtrim**(i,**2**) + ' of ' + **strtrim**(files.*length*-**1**,**2**) + ' 个图像配准'
**endfor**
**print**, '已配准所有图像！'
**print**, ''
;关闭输出栅格
reference.**close**
;转到配准后的目录并为镶嵌打开每张配准后的栅格
**cd**, warpeddir
scenes = **objarr**(**n\_elements**(files))
**for** i=**0**, **n\_elements**(files)-**1** **do** scenes[i] = e.**openraster**(files[i])
**print**, '正在生成镶嵌图...'
; 从 ENVITasks 目录中获取任务
Task = **ENVITask**('BuildMosaicRaster')
; 定义输入
Task.*INPUT\_RASTERS* = scenes
Task.*COLOR\_MATCHING\_METHOD* = 'Histogram Matching'
Task.*COLOR\_MATCHING\_STATISTICS* = 'Entire Scene'
Task.*FEATHERING\_METHOD* = 'Edge'
Task.*FEATHERING\_DISTANCE* = **make\_array**(**n\_elements**(files), /INTEGER, VALUE = **40**)
; 定义输出
outfile = thisdir + **path\_sep**() + 'pluto\_mosaic'
**if** **file\_test**(outfile) **then** **file\_delete**, outfile
Task.*OUTPUT\_RASTER\_URI* = outfile
; 运行任务
Task.**Execute**
;将栅格导出为 TIF 图像
outfile = thisdir + **path\_sep**() + 'pluto\_mosaic.tif'
**if** **file\_test**(outfile) **then** **file\_delete**, outfile
(Task.*OUTPUT\_RASTER*).**export**, outfile, 'TIFF'
(Task.*OUTPUT\_RASTER*).**close**
;返回到第一个目录
**cd**, first\_dir
**print**, '镶嵌完成！'
;关闭 ENVI
e.**close**
**end**
要运行此代码，只需创建一个名为 "pluto_mosaicking.pro" 的文本文件，并将代码复制/粘贴到文件中。运行程序时，IDL 将自动下载冥王星图像、创建目录，并在与 pluto_mosaicking.pro 文件相同的目录中生成镶嵌输出。如果下载文件时出现问题，您可能需要从 NASA 网站手动下载，地址是 https://pluto.jhuapl.edu/soc/Pluto-Encounter/。代码顶部包含了如何运行程序的附加说明，示例中的注释也解释了所采取的步骤。