【IDL】实现ENVI默认显示效果（线性拉伸算法）

发表于 2011-08-25

---

ENVI下无论打开任何一种图像，默认效果都特别的漂亮，而用IDL写代码直接显示的图像，在色彩、亮度方面都不如ENVI好看。而关于ENVI的2%裁剪算法，也没有详细的文档说明，具体原理是什么？应该怎么用IDL实现呢。

据个人字面理解，2%裁剪拉伸是基于直方图分布，对图像DN值分布在2%和98%之间的做线性拉伸。即拉伸时去除小于2%和大于98%的值，这样绝大多数的异常值会在拉伸时舍掉，显示出漂亮直观的效果。

存在一个问题，万一没有正好处于2%和98%的时候怎么办？

找了个影像仔细看下，默认[Scroll] Linear2%的拉伸情况下看直方图，发现下面特点

拉伸起点值为4，直方图统计分布在1.82%，下一个值5所在位置为2.79%；

拉伸终点为34，直方图统计分布在98.14%，前一个值33所在位置为97.63%；

故可认定ENVI的2%线性裁剪拉伸是按照2%和98%为界限，选取最临近分布的值分别作为拉伸数据范围。

编写IDL的测试代码如下：

FUNCTION LINEAR2, inImage,inPer

  COMPILE_OPT idl2

  ;

  IF N_ELEMENTS(inPer) EQ 0 THEN inPer = 0.02

 

  sz = SIZE(inImage)

  ;图像太大的话统计太慢，模仿ENVI，仅对Scroll窗口进行统计，默认大小是256*256

  IF sz[0] EQ 2 THEN BEGIN

   image = CONGRID(inImage, 256, 256)

  ENDIF ELSE BEGIN

   image = CONGRID(inImage, sz[1], 256, 256)

  ENDELSE

  ;图像基本信息

  sz = SIZE(image)

  IF sz[0] EQ 2 THEN BEGIN

   nPlanes = 1

    x = sz[1]

    y = sz[2]

  ENDIF ELSE BEGIN

   nPlanes = 3

    x = sz[2]

    y = sz[3]

  ENDELSE

 

  outImage = inImage

 

  FOR i=0, nPlanes-1 DO BEGIN

   IF nPlanes EQ 3 THEN img = REFORM(image[i,*,*]) ELSE img=image

   ;直方图统计

   array = HISTOGRAM(img,oMax = maxV,oMin = minV)

   arrnumb= N_ELEMENTS(array)

   ;

   percent = TOTAL(array,/CUMULATIVE)/TOTAL(array)

   idx1 = WHERE(percent LE inPer)

   idx2 = WHERE(percent GE inPer)

   number = N_ELEMENTS(idx1)

   ;计算当前inpert对应的数值 （2%）

   ;两个索引，取临近索引

   curIdx = (ABS(percent[idx1[number-1]]-inPer) LE ABS(percent[idx2[0]]-inPer))? idx1[number-1]:idx2[0]

   minvalue = minV +(maxV-minV)*curIdx/(arrnumb-1)

   ;1-inper对应数值 （98%）

   idx1 = WHERE(percent LE (1-inPer))

   idx2 = WHERE(percent GE (1-inPer))

   number = N_ELEMENTS(idx1)

   ;两个索引，取临近索引

   curIdx = (ABS(percent[idx1[number-1]]-1+inPer) LE ABS(percent[idx2[0]]-1+inPer))? idx1[number-1]:idx2[0]

   maxvalue = minV +(maxV-minV)*curIdx/(arrnumb-1)

   ;单波段还是多波段

   IF nPlanes EQ 3 THEN $

     outImage[i,*,*] = BYTSCL(outImage[i,*,*], max=maxvalue, min=minvalue) $

   ELSE outImage = BYTSCL(outImage, max=maxvalue, min=minvalue)

  ENDFOR

 

  IF nPlanes EQ 1 THEN outImage = REFORM(outImage)

 

  RETURN, outImage

END