ENVI二次开发的内部（猜测）示例代码分析
之前从RSI网站（现为ITTvis）下载的一部分源码，根据某些信息，看着应该是内部人士写的，从格式、写法上仔细拜读了一下，果然，写的非常严谨。
源码：
;--------------------------------------------------------+
; NAME
;
; GEOREF_RSAT
;
; TYPE
;
; ENVI User Function
;
; CALLING SEQUENCE
;
; GEOREF_RSAT, event
;
; PURPOSE
;
; An ENVI User Function to automatically georeference
; (register) radarsat format image data using the embedded
; ground control points.
;
; MODIFICATION HISTORY
;
; April 1999 -- written under ENVI 3.0, IDL 5.03
; September 2000 -- added About Button and a check for missing GCPs
;
; DISTRIBUTION INFO
;
; Distributed from the RSI Web Site, see Tech Tip
; ENVI249 "Automatically Registering Radarsat
; Imagery from Embedded GCPs" at:
; http://www.rsinc.com/services/search.cfm
;
; AUTHOR INFO
;
; David Gorodetzky
; RSI Profession Services Group
; goro@rsinc.com
;
; PROGRAM NOTES
;
; According to CCRS documentation available at the time this
; routine was written, the lat/lon of the first, middle and
; last pixel in each line of the radarsat image is embedded
; into the line prefix in bytes 132 to 156. The geocoded
; data are 4 byte integers arranged as follows:
;
; 1st_lat, middle_lat, last_lat, 1st_lon, middle_lon, last_lon
;
; The units for the lat/lon data are millionths of a
; degree. Note that the accuracy of the registration
; is limited by the accuracy of the embedded GCPs, which
; are estimated to be about 400 m. Also, the datum of
; the lat/lon coordinates is unknown and assumed to be
; WGS84. The use of the embedded geocoded data to georerference
; the image appears to be a good first step in getting the geometry
; correct, although the absolute positional accuracy of the
; registered result is often in error by as much as 1500 m.
; Translation errors like this can be corrected if an additional
; source od ground control is available, such as an accurate
; vector layer.
;
; This file contains the following routines in this order:
;
; GEOREF_RSAT_EVENT
;
; The event handler for the first GUI. It stores all user
; defined data in a pointer called info.
;
; GEOREF_RSAT
;
; The processing routine. It also makes the second GUI
; which is auto-managed by ENVI.
;
;---------------------------------------------------------------+
FORWARD_FUNCTION ENVI_FILE_TYPE, WIDGET_AUTO_BASE, WIDGET_PMENU, $
AUTO_WID_MNG, WIDGET_PARAM, WIDGET_PROJ, WIDGET_OUTF, WIDGET_MENU, $
WIDGET_OUTFM, ENVI_TRANSLATE_PROJECTION_UNITS, MAGIC_MEM_CHECK, $
ENVI_REAL_MAIN_BASE
;----------------------------------------------------------------
PRO GEOREF_RSAT_EVENT, event
;
; Event handler for the first GUI.
;
;---------------------------------------------------------------+
WIDGET_CONTROL, event.top, get_uvalue=info
WIDGET_CONTROL, event.id, get_uvalue=what
CASE what OF
'OK': BEGIN
;
; transfer user inupted data into info pointer
;
WIDGET_CONTROL, (*info).num_gcp, get_value=skip
(*info).skip = skip
WIDGET_CONTROL, (*info).proj, get_value=sel_proj
(*info).sel_proj = sel_proj
WIDGET_CONTROL, (*info).outf, get_value=outfile
(*info).outfile = outfile
WIDGET_CONTROL, (*info).pts_only, get_value=get_pts_only
(*info).get_pts_only = get_pts_only
;
; destroy the GUI and return
;
WIDGET_CONTROL, event.top, /destroy
RETURN
END
'CANCEL': BEGIN
(*info).cancel = 1
WIDGET_CONTROL, event.top, /destroy
RETURN
END
'ABOUT': BEGIN
;
; display about message text
;
msg = STRARR(28)
msg[0] = "GEOREF_RSAT User Function"
msg[1] = "----------------------------------"
msg[2] = ""
msg[3] = "Written by David Gorodetzky, goro@rsinc.com"
msg[4] = 'Research Systems, Inc. (RSI)'
msg[5] = "Professional Services Group (PSG)"
msg[6] = "http://www.rsinc.com"
msg[7] = ""
msg[8] = "For additional information on this function please"
msg[9] = "see the Tech Tip article 'Automatically Registering"
msg[10] ="Radarsat Imagery from Embedded GCPs' on the RSI web"
msg[11] ="site at http://www.rsinc.com/services/search.cfm"
msg[12] =""
msg[13] ="For more information about consulting services"
msg[14] ="related to ENVI, IDL or remote sensing contact"
msg[15] ="RSI at (303) 786-9900 and ask for the"
msg[16] ="Professional Services Group, or send an e-mail"
msg[17] ="to << consult@rsinc.com >>"
msg[18] = ""
msg[19] ='DISCLAIMER'
msg[20] = ""
msg[21] =' This additional functionality is provided free'
msg[22] =' of charge as a service for our ENVI users.'
msg[23] =' However, this procedure is not supported by'
msg[24] =' Research Systems Technical Support. The'
msg[25] =' procedure GEOREF_RSAT has been tested'
msg[26] =' and we believe that it works correctly, but it'
msg[27] =' is in no way guaranteed.'
ENVI_INFO_WID, msg, title='About RSI User Functions'
END
ELSE:
ENDCASE
END
;----------------------------------------------------------------
PRO GEOREF_RSAT, event
;
; The processing routine. This also makes the second GUI
; which is auto-managed by ENVI.
;
;---------------------------------------------------------------+
; use CATCH to prevent unexpected errors from crashing ENVI
;
!error_state.code = 0
CATCH, error
IF (error NE 0) THEN BEGIN
IF (N_ELEMENTS(rep_base) NE 0) THEN ENVI_REPORT_INIT, base=rep_base, /finish
ok = DIALOG_MESSAGE(!error_state.msg, /error, /cancel)
IF (STRUPCASE(ok) EQ 'CANCEL') THEN BEGIN
IF (N_ELEMENTS(lun) NE 0) THEN FREE_LUN, lun
RETURN
ENDIF
ENDIF
; select the input image, disable spatial and spectral
; subsetting and force the user to select by band -- also
; require that the input image is a RADARSAT file type
;
ft = ENVI_FILE_TYPE("RADARSAT")
ENVI_SELECT, title='Select RADARSAT Input File', $
fid=fid, /no_dims, /no_spec, /band_only,$; file_type=ft, $
dims=dims, pos=pos
IF (fid EQ -1) THEN RETURN
; get some basic info about the file
;
ENVI_FILE_QUERY, fid, ns=ns, nl=nl, offset=offset, $
fname=fname, data_type=dt, sname=sname
OPENR, lun, fname, /Get_LUN
fsize = FSTAT(lun)
; determine the number of bytes per image pixel (bpip)
; based on the difference between the image size and
; the file size, accounting for the line prefix of 192 bytes
;
CASE 1 OF
fsize.size EQ ( (ns+192)*(nl) + offset): bpip = 1
fsize.size EQ ( ((ns*2)+192)*nl + offset): bpip = 2
ELSE: BEGIN
extra_tlb = WIDGET_AUTO_BASE(title='specify data type')
txt1 = 'Unable to compute the data type of the imagery.'
txt2 = 'Select the number of bytes per image pixel. If'
txt3 = 'the data file was scaled to byte when it was'
txt4 = 'opened remember to choose the original data type.'
sbase0 = WIDGET_BASE(extra_tlb, /col, xsize=255, /frame)
msg = WIDGET_LABEL(sbase0, value=txt1, /align_left)
msg = WIDGET_LABEL(sbase0, value=txt2, /align_left)
msg = WIDGET_LABEL(sbase0, value=txt3, /align_left)
msg = WIDGET_LABEL(sbase0, value=txt4, /align_left)
sbase1 = WIDGET_BASE(sbase0, /row)
choose = WIDGET_PMENU(sbase1, /auto, uvalue='bpip', default=0, $
list = ['1 (byte)', '2 (integer)', '4 (long integer/float)'], $
prompt='bytes per image pixel')
lookup = [1,2,4]
extra_result = AUTO_WID_MNG(extra_tlb)
IF (extra_result.accept EQ 0) THEN RETURN
bpip = lookup(extra_result.bpip)
END
ENDCASE
; create the GUI
;
ENVI_CENTER, x, y
TLB = WIDGET_BASE(/col, xoffset=x, yoffset=y-100, /modal, $
group_leader=envi_real_main_base(), title='Extract RADARSAT GCP Parameters')
sbase1 = WIDGET_BASE(TLB, /frame, /col)
line_skip_txt1 = 'Each image line yields 3 GCPs (first, middle, and last pixel).'
line_skip_txt2 = 'Set the number of lines to skip between extracting GCPs.'
label2 = WIDGET_LABEL(sbase1, value=line_skip_txt1, /align_left)
label3 = WIDGET_LABEL(sbase1, value=line_skip_txt2, /align_left)
rbase1 = WIDGET_BASE(sbase1, /row)
num_gcp = WIDGET_PARAM(rbase1, prompt='line skip factor ', $
uvalue = 'skip', xs=3, dt = 2, default=10, floor=0, $
ceil=(nl/2), /all_events)
proj = WIDGET_PROJ(TLB, proj=proj, uvalue='proj', $
prompt='select output projection for registeration', /frame)
sbase2 = WIDGET_BASE(TLB, /frame, /col)
def_fname = STRMID(sname, 0, RSTRPOS(sname, '.'))+'.pts'
outfile = WIDGET_OUTF(sbase2, prompt='output .pts file name', $
uvalue='outfile', default=def_fname)
sbase3 = WIDGET_BASE(sbase2, /row)
pts_only = WIDGET_MENU(sbase2, uvalue='pts_only', $
list=["don't georeference, just write points file"])
row_base = WIDGET_BASE(TLB, /row)
ok = WIDGET_BUTTON(row_base, value=' OK ', uvalue='OK')
can = WIDGET_BUTTON(row_base, value='Cancel', uvalue='CANCEL')
about = WIDGET_BUTTON(row_base, value='About This Function', uvalue='ABOUT')
WIDGET_CONTROL, TLB, /realize
info_data = { $
num_gcp:num_gcp, $ ; widget ID for line skip factor widget_param
skip:10, $ ; default line skip factor
proj:proj, $ ; widget ID for widget_proj
sel_proj:{envi_proj_struct}, $; projection selection
outf:outfile, $ ; widget ID for the OUTF widget for pts file
outfile:def_fname, $ ; default .pts filename
pts_only:pts_only, $ ; widget ID for the checkbox for doing pts only
get_pts_only:0, $ ; 0 = georef, 1=only get the .pts
CANCEL:0 $ ; cancel instructions
}
info = PTR_NEW(info_data, /no_copy)
WIDGET_CONTROL, TLB, set_uvalue=info
XMANAGER, 'GEOREF_RSAT', TLB, /no_block
; return if the user cancels
;
IF (*info).cancel EQ 1 THEN RETURN
IF ( (*info).get_pts_only NE 1) THEN BEGIN ; if doing the full registration
TLB = WIDGET_AUTO_BASE(title='Registration Parameters', /ybig)
sbase0 = WIDGET_BASE(TLB, /col, /frame)
sbase1 = WIDGET_BASE(sbase0, /row)
model = WIDGET_PMENU(sbase1, /auto, default=2, prompt='warp model', $
list=['RST', 'Polynomial', 'Triangulation'], uvalue='model')
sbase2 = WIDGET_BASE(sbase0, /row)
resample = WIDGET_PMENU(sbase2, /auto, default=0, prompt='resampling method', $
list=['Nearest Neighbor', 'Bilinear Interpolation', 'Cubic Convolution'], $
uvalue='resample')
sbase4 = WIDGET_BASE(sbase0, /row)
back = WIDGET_PARAM(sbase4, prompt='background value', $
uvalue = 'back', xs=5, dt = 4, field=1, default=0, /auto)
sbase5 = WIDGET_BASE(sbase0, /row)
pixX = WIDGET_PARAM(sbase5, /auto, prompt='Output Pixel Sizes: X ', $
uvalue='pixX', xs=5, dt=4, field=1, default=100)
pixY = WIDGET_PARAM(sbase5, /auto, prompt=' Y ', $
uvalue='pixY', xs=5, dt=4, field=1, default=100)
warpfile = WIDGET_OUTFM(TLB, /auto, uvalue='warpfile', $
/frame, prompt='Registered Image Filename')
result2 = AUTO_WID_MNG(TLB)
IF result2.accept EQ 0 THEN RETURN
ENDIF
; convert the widget info back into the correct data type
;
(*info).skip = LONG((*info).skip)
; set up varaibles and arrays needed to extract the GCPs
;
num_gcp_lines = LONG( (nl-1) / ((*info).skip+1) + 1 )
full_line = (ns*bpip)+192L
temp = LONARR(6)
lat = LONARR(3,num_gcp_lines)
lon = LONARR(3,num_gcp_lines)
pix_y = LONARR(3,num_gcp_lines)
; set up the status reporting widget
;
rstr = ['Input File: '+fname, 'Output File: '+(*info).outfile]
ENVI_REPORT_INIT, rstr, base=rep_base, /interupt, $
title='extracting GCPs from RADARSAT'
ENVI_REPORT_INC, rep_base, num_gcp_lines+3
ENVI_REPORT_STAT, rep_base, 1, num_gcp_lines+3
; read in the first image line's GCPs
;
POINT_LUN, lun, offset+132
READU, lun, temp
IF ( BYTE(255,1) EQ 0 ) THEN temp = SWAP_ENDIAN(temp)
lat(*,0) = temp(0:2)
lon(*,0) = temp(3:5)
pix_y(*,0) = 1
; read in the rest of the GCPs, except for the last line
;
line = (*info).skip + 2
file_pointer = ( (line-1)*full_line ) + 132 + offset
FOR array_row=1, num_gcp_lines-2 DO BEGIN
;
; update the status reporting widget and
; check for the function being cancelled
;
ENVI_REPORT_STAT, rep_base, array_row+1, num_gcp_lines+3, cancel=cancel
IF (cancel) THEN BEGIN
FREE_LUN, lun
ENVI_REPORT_INIT, base=rep_base, /finish
RETURN
ENDIF
;
; extract the rest of the GCPs except for the last line
;
POINT_LUN, lun, file_pointer
READU, lun, temp
IF ( BYTE(255,1) EQ 0 ) THEN temp = SWAP_ENDIAN(temp)
lat(*,array_row) = temp(0:2)
lon(*,array_row) = temp(3:5)
pix_y(*,array_row) = line
line = line + (*info).skip + 1
file_pointer = file_pointer + ((*info).skip+1)*full_line
ENDFOR
; read in the last image line's GCPs
;
; first, update the status report widget
;
ENVI_REPORT_STAT, rep_base, array_row+1, num_gcp_lines+3, cancel=cancel
IF (cancel) THEN BEGIN
FREE_LUN, lun
ENVI_REPORT_INIT, base=rep_base, /finish
RETURN
ENDIF
;
lastline = offset + (nl-1)*full_line + 132
POINT_LUN, lun, lastline
READU, lun, temp
IF ( BYTE(255,1) EQ 0 ) THEN temp = SWAP_ENDIAN(temp)
lat(*,array_row) = temp(0:2)
lon(*,array_row) = temp(3:5)
pix_y(*,array_row) = nl
; don't forget to free the file pointer
;
FREE_LUN, lun
; make the array that will become the .pts file
;
ENVI_REPORT_STAT, rep_base, array_row+1, num_gcp_lines+3, cancel=cancel
IF (cancel) THEN BEGIN
ENVI_REPORT_INIT, base=rep_base, /finish
RETURN
ENDIF
GCPs = DBLARR(4, num_gcp_lines*3)
;
; define the x-pixel positions (1st, middle, last)
;
pix_x = LONARR(3,num_gcp_lines)
pix_x(0,0) = 1L
pix_x(1,0) = LONG(ns/2)
pix_x(2,0) = ns
pix_x(0,*) = pix_x(0,0)
pix_x(1,*) = pix_x(1,0)
pix_x(2,*) = pix_x(2,0)
;
; convert the lat/lons into the selected output projection,
; remember to convert them from millionths to degrees first
;
ENVI_PROJ_CONVERT, map_x, map_y, REFORM( (lat/1e6), num_gcp_lines*3), $
REFORM( (lon/1e6), num_gcp_lines*3), projection=(*info).sel_proj, /to_map
;
; if the file is missing the embedded GCPs then they usually
; end up being zeros, so do a quick check to see if the data
; are OK -- if more than 50% of the GCPs are zero then
; warn the user
;
w1 = WHERE(map_x EQ 0, c1)
w2 = WHERE(map_y EQ 0, c2)
IF ( c1 GT (0.5*3*num_gcp_lines) OR c2 GT (0.5*3*num_gcp_lines) ) THEN BEGIN
warn = DIALOG_MESSAGE("Most of the GCPs are zeros! You're Radarsat image "+$
'is probably missing the GCPs. Do you want to Continue?', /question)
IF STRUPCASE(warn) EQ 'NO' THEN BEGIN
FREE_LUN, lun
ENVI_REPORT_INIT, base=rep_base, /finish
RETURN
ENDIF
ENDIF
;
; put the GCP data into the .pts array
;
GCPs(0,*) = map_x
GCPs(1,*) = map_y
GCPs(2,*) = REFORM(pix_x, num_gcp_lines*3)
GCPs(3,*) = REFORM(pix_y, num_gcp_lines*3)
; write the .pts file to disk
;
OPENW, lun, (*info).outfile, /Get_LUN
ENVI_REPORT_STAT, rep_base, array_row+1, num_gcp_lines+3, cancel=cancel
IF (cancel) THEN BEGIN
FREE_LUN, lun
ENVI_REPORT_INIT, base=rep_base, /finish
RETURN
ENDIF
PRINTF, lun, 'ENVI Registration GCP File'
PRINTF, lun, GCPs, format="(32000(2(F20.2), 2(I20)/))"
FREE_LUN, lun
ENVI_REPORT_INIT, base=rep_base, /finish
; as long as the user doesn't choose to write
; the points file only, do the warp
;
IF ((*info).get_pts_only NE 1) THEN BEGIN
; x0 and y0 are the map coordinates of the upper left
; pixel of the output warped image -- these should be
; set to the Westernmost (X) and Northernmost (Y) GCPs
;
x0 = MIN( GCPs(0,*) )
y0 = MAX( GCPs(1,*) )
; xsize and ysize are the size of the ouput image in map coordinates
;
xsize = MAX(GCPs(0,*)) - MIN(GCPs(0,*))
ysize = MAX(GCPs(1,*)) - MIN(GCPs(1,*))
; define the map information structure
;
map_info = {envi_map_struct}
map_info.proj = (*info).sel_proj
map_info.ps = [result2.pixX, result2.pixY]
map_info.mc = [0, 0, x0, y0]
CASE map_info.proj.type OF
1: map_info.proj.units = ENVI_TRANSLATE_PROJECTION_UNITS('degrees') ; Geographic Lat/Lon gets degrees
8: map_info.proj.units = ENVI_TRANSLATE_PROJECTION_UNITS('feet') ; State Plane gets feet
ELSE: map_info.proj.units = ENVI_TRANSLATE_PROJECTION_UNITS('meters') ; everything else gets meters
ENDCASE
; if the user wants to output the result
; to memory, make sure it will not exceed
; ENVI's cache limit
;
mem_dims = [-1L, 0, LONG(xsize/result2.pixX), $
0, LONG(ysize/result2.pixY) ]
mem_result = MAGIC_MEM_CHECK(dims=dims, fid=fid, nb=N_ELEMENTS(pos), $
in_memory=result2.warpfile.in_memory, out_dt=dt, $
out_name=result2.warpfile.name)
IF (mem_result.cancel EQ 1) THEN RETURN
; set the warp method
;
CASE result2.model OF
0: BEGIN
CASE result2.resample OF
0: method=0 ; RST w/nearest neighbor
1: method=1 ; RST w/bilinear
2: method=2 ; RST w/cubic convolution
ELSE:
ENDCASE
END
1: BEGIN
CASE result2.resample OF
0: method=3 ; Polynomial w/nearest neighbor
1: method=4 ; Polynomial w/bilinear
2: method=5 ; Polynomial w/cubic convolution
ELSE:
ENDCASE
END
2: BEGIN
CASE result2.resample OF
0: method=6 ; Triangulation w/nearest neighbor
1: method=7 ; Triangulation w/bilinear
2: method=8 ; Triangulation w/cubic convolution
ELSE:
ENDCASE
END
ELSE:
ENDCASE
; make sure the correct save file is restored
;
ENVI_CHECK_SAVE, /registration
; do the registration warp, note that if the user chooses
; to do a Polynomial warp, they always get 1st order
;
REG_WARP_DOIT, fid=fid, pos=pos, dims=dims, map_info=map_info, $
out_name=mem_result.out_name, method=method, degree=1, $
in_memory=mem_result.in_memory, pts=GCPs, background=result2.back, $
to_image=0, xsize=xsize, ysize=ysize, xstart=0, ystart=0, $
x0=x0, y0=y0, pixel_size=map_info.ps
ENDIF
END