Creating Mosaics from Pluto Imagery

原文链接: https://www.nv5geospatialsoftware.com/Learn/Blogs/Blog-Details/creating-mosaics-from-pluto-imagery

18800 Rate this article:

No rating

Creating Mosaics from Pluto Imagery

Zachary Norman Monday, January 25, 2016

For a good portion of my life I have been a space nerd and I have always enjoyed looking at images of other planets in our solar system, nebulae, and other galaxies. Recently there was a new addition to the group of photographed celestial objects: Pluto. On July 14th, 2015 the New Horizons probe flew within 12,600 km of the surface of Pluto and took some awesome images. With these new images, I thought it would be a neat idea to try and mosaic a few of the images together using ENVI to see if it could be done.

In order to get some images that would be suitable for mosaicking, I needed to make sure that they were taken in quick succession so that they would have a good amount of overlap with each other. This allows ENVI to be able to register the images to one another without too much trouble. To find these images I got data from New Horizon's instrument LORRI at the following link:

https://pluto.jhuapl.edu/soc/Pluto-Encounter/

It was really important when selecting images to mosaic that the images needed to have similar acquisition times. After looking through the available images from LORRI, I found 13 images from the acquisition time from 11:36:00 UTC to 11:36:36 UTC. It is great that the images of Pluto are available, but unfortunately there was no metadata associated with the data to accurately geolocate the images. To work around this missing metadata, I decided that I would just create a dummy spatial reference centered around a reference image. For my base image, I chose to use the picture taken at 11:36:00 UTC and gave this image a basic WGS-84 coordinate system with fake pixel sizes of 1 meter by 1 meter. At this point it is important to note that my fake coordinate system and pixel sizes likely do not correspond at all to the actual coordinate system or pixel sizes, I just needed to create a spatial reference that I could use to register and mosaic the images using ENVI.

With my reference image and fake spatial reference, the next step was to figure out how to get the images registered properly to one another. After doing some investigation with the Image Registration Workflow, I realized that in order to get the best results I needed to have the images placed as close to their actual relative positions as possible. To do this, I looked at my reference image and first image to be registered in the series (reference image had the acquisition time of 2015-07-14 11:36:00 UTC and the first image to be registered was acquired at 2015-07-14 11:36:03 UTC).

After looking closely at the two images and tweaking things by hand a bit, I found that the approximate pixel offset of the second image with respect to the first was about [-580, -18] pixels in the x and y directions. This mean that the first image in series after the base image had an offset of [-580, -18] pixels and the second image had an offset of [-1160, -36] pixels. At this point I decided to use the ENVI API to programmatically apply this offset to each image in the series. A big reason for this choice to use the ENVI API is that I wouldn't have to step through the Image Registration Workflow for every image in my series and, if I changed the offset at all, it would be quick to apply to a different image.

After I re-projected each image based on my pixel offset the images were relatively close at this point. All that was left was to apply color balancing (through a custom histogram matching function), register the images to one another, and mosaic them together. Here is a quick summary of the steps that I performed using the ENVI API:

Apply color matching to reference image -> Use task GenerateTiePointsByCrossCorrelation -> Use task FilterTiePointsByGlobalTransform -> Use task ImageToImageRegistration -> Use task BuildMosaicRaster.

Below is the output mosaic and the IDL code used to create the image.

;+

; :Author: Zach Norman

;

; :Description:

;   Small example program to register and mosaic some Pluto images from the New Horizon's

;   LORRI instrument. Save this code as pluto_mosaicking.pro and to run you can press 'Run' 

;   in the IDL Workbench. See the output section below for where to find the results of the

;   program and details on the results.

;  

;   The code is written to automatically download the images of Pluto to a directory called images

;   that will be created in the same location that this .pro file will reside. If there are any

;   issues downloading the files then you may have to download them from the web by hand. This

;   could be necessary if the URLs change for the images from NASA. You may need to

;   make some additional changes as well if you download the images by hand.

;  

;   Specifically, you may have to change some of the parameters in the procedure pluto_mosaicking

;   so that they are current for your machine. If you download the Pluto images from

;   http://pluto.jhuapl.edu/soc/Pluto-Encounter/index.php then you will need to specify

;   the value of "imagedir" in the pluto_mosaicking procedure.

;  

;   If you want to or have to download the images by hand, there are 13 images which have

;   acquisition times ranging from 2015-07-14 11:36:00 UTC to 2015-07-14 11:36:36 UTC. As

;   is, the code is written so that you place this .pro file next to the directory which

;   contains the images.

;  

;  

; :Output:

;   Three files will be generated in the same directory as this .pro code and they will be

;   named: pluto_mosaic, pluto_mosaic.hdr, pluto_mosaic.tif.

;  

;   The first file is an ENVI format image that can be read into ENVI with the header file (.hdr)

;   and the third file is a TIFF images so that you can view the mosaic with a general photo viewer.

;   Note that the ENVI header file contains fake projection information. I created a dummy spatial

;   reference for each image since there was no metadata provided with the LORRI images.

;  

;  

;-

;small function for the histogram matching between images

;this is useful because some of the images have bad automatic scaling from NASA and

;they appear very dark

function histogram_matching, match, base

compile_opt idl2

maxval = (max(match))>(max(base))

;calculate the histograms and cumulatide densities

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)

;create the transform function z to switch our data space in the match image

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

;do the histogram matching

matched = z[base]

return, matched

end

;method for ENVI to close all data sources and

pro ENVI::reset

compile_opt idl2

e = envi(/current)

if (e eq !NULL) then return

if (e.WIDGET_ID ne 0) then begin

;reset our views

views = e.getview(/all)

foreach view,views do view.close

endif

;close all open rasters

opendata = (e.data).Get(/raster,/vector, /series)

foreach raster, opendata do raster.close

end

pro pluto_mosaicking

compile_opt idl2

;start ENVI and make it a new state!

e = envi(/current)

if (e eq !NULL) then e = envi(/headless) else e.reset

;approximate pixel offset for each image

px_off = [-580, -18]

;directories for running the code

thisdir = FILE_DIRNAME( (((scope_traceback(/STRUCTURE))[-1]).FILENAME))

;imagedir is the location of the original images

imagedir = thisdir + path_sep() +  'images'

;copydir is the location of copies of the images from imagedir that have a dummy spatial reference

copydir = thisdir + path_sep() +  'image_copies'

;warpeddir is where the registered images go

warpeddir = thisdir + path_sep() + 'images_warped'

;if imagedir does not exist, then make it and download the files to that directory

;if the URLs change  then these will need to be modified by hand or downloaded by hand from the following URL:

; http://pluto.jhuapl.edu/soc/Pluto-Encounter/index.php

if ~file_test(imagedir) then begin

print, 'Did not find the directory ' + imagedir

print, string(9b) + 'Creating directory and downloading images...'

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']

;use the IDLnetURL object to download each image  

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) + 'Downloaded image ' + strtrim(i+1,2) + ' of ' + strtrim(n_elements(urls),2)

endfor

urlobj.cleanup

print, string (9b) + 'Downloaded all images!'

endif

;make our output directories if they don't exist

if ~file_test(copydir) then FILE_MKDIR, copydir

if ~file_test(warpeddir) then FILE_MKDIR, warpeddir

;go to imagedir and look for the Pluto images

cd, imagedir, CURRENT = first_dir

files = file_search('*.jpg')

;check to see if we already made copies fo the images

;create a fake spatial reference

psize = 1d*(180/!CONST.R_EARTH);pixel size of 1 meter

spatialRef = ENVIStandardRasterSpatialRef( $

COORD_SYS_CODE=4326, /GEOGCS, $

PIXEL_SIZE=[psize,psize], TIE_POINT_PIXEL=[0.0D,0.0D], $

TIE_POINT_MAP=[0D,0D])

;make sure to sort the files so that we are arranged by time

files = files[sort(files)]

;copy all of the files and create ENVI rasters witha  base spatial reference

reference = e.openraster(files[0], SPATIALREF_OVERRIDE=spatialref)

reference_data = reference.GetData()

;export original raster to copy raster

outfile = copydir + path_sep() + file_basename(files[0], '.jpg')

if file_test(outfile) then FILE_DELETE, outfile

reference.export, outfile, 'ENVI'

;perform histogram matching on all images and export them to the 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

;go to the copy directory and register our images to one another

cd, copydir

files = file_basename(files, '.jpg')

reference = e.openraster(files[0])

;copy the reference image to the copy directory and open it in 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)

;perform image registration on all images

print, 'Found ' + strtrim(files.length,2) + ' images to register!'

for i=1, n_elements(files)-1 do begin

warp_this = e.openraster(files[i])

; Get the auto tie point generation task from the catalog of ENVITasks

TieTask = ENVITask('GenerateTiePointsByCrossCorrelation')

TieTask.INPUT_RASTER1 = reference

TieTask.INPUT_RASTER2 = warp_this

TieTask.Execute

TieTask_TiePoints = TieTask.OUTPUT_TIEPOINTS

; Get the tie point filter task from the catalog of ENVITasks

FilterTask = ENVITask('FilterTiePointsByGlobalTransform')

FilterTask.INPUT_TIEPOINTS = TieTask_TiePoints

FilterTask.Execute

FilterTask_TiePoints = FilterTask.OUTPUT_TIEPOINTS

; Get the image-to-image registration task from the catalog of 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

;close the raster that was warped

warp_this.close

;delete tie points files

FILE_DELETE, TieTask.OUTPUT_TIEPOINTS_URI

FILE_DELETE, FilterTask.OUTPUT_TIEPOINTS_URI

print, string(9b) + 'Completed ' + strtrim(i,2) + ' of ' + strtrim(files.length-1,2) + ' image registrations'

endfor

print, 'Registered all images!'

print, ''

;close the output raster

reference.close

;got to the warped directory and open each warped raster for the mosaic

cd, warpeddir

scenes = objarr(n_elements(files))

for i=0, n_elements(files)-1 do scenes[i] = e.openraster(files[i])

print, 'Generating mosaic...'

; Get the task from the catalog of ENVITasks

Task = ENVITask('BuildMosaicRaster')

; Define inputs

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)

; Define outputs

outfile = thisdir + path_sep() + 'pluto_mosaic'

if file_test(outfile) then file_delete, outfile

Task.OUTPUT_RASTER_URI = outfile

; Run the task

Task.Execute

;export the raster as a TIF image

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

;return to first directory

cd, first_dir

print, 'Mosaic complete!'

;close ENVI

e.close

end

To run the code, simply create a text file called "pluto_mosaicking.pro" and copy/paste the code into the file. When you run the program, IDL will automatically download the Pluto images, create directories, and produce the mosaicked output in the same directory as the pluto_mosaicking.pro file. If there is an issue while downloading the files, you may have to manually download them from NASA at https://pluto.jhuapl.edu/soc/Pluto-Encounter/. Additional instructions are included at the top of the code for how to run the program along with comments throughout the example to explain the steps taken.

Introducing the ENVIHydratable Interface Using IDL_IDLBridge to run multiple processes in parallel