IDL实现PCA图像融合

        PCA原理网上很多，这里直接给出实现代码

pro pca
  Compile_opt idl2
;  workdir=SourceRoot()
  
  img_file_b1 = 'E:\yaoganhomework\ch06-4his变换融合idl程序\ch06-4 his变换融合\data\tb1.bmp'
  img_file_b2 = 'E:\yaoganhomework\ch06-4his变换融合idl程序\ch06-4 his变换融合\data\tb2.bmp'
  img_file_b3 = 'E:\yaoganhomework\ch06-4his变换融合idl程序\ch06-4 his变换融合\data\tb3.bmp'
  img_file_b4 = 'E:\yaoganhomework\ch06-4his变换融合idl程序\ch06-4 his变换融合\data\tb4.bmp'
  img_file_b5 = 'E:\yaoganhomework\ch06-4his变换融合idl程序\ch06-4 his变换融合\data\tb5.bmp'
  img_file_b6 = 'E:\yaoganhomework\ch06-4his变换融合idl程序\ch06-4 his变换融合\data\tb6.bmp'
  img_file_pan = 'E:\yaoganhomework\ch06-4his变换融合idl程序\ch06-4 his变换融合\data\spotsub.BMP'
  
  tm_data_b1 = READ_BMP(img_file_b1)
  tm_data_b2 = READ_BMP(img_file_b2)
  tm_data_b3 = READ_BMP(img_file_b3)
  tm_data_b4 = READ_BMP(img_file_b4)
  tm_data_b5 = READ_BMP(img_file_b5)
  tm_data_b6 = READ_BMP(img_file_b6)
  pan_data = READ_BMP(img_file_pan)
  
  dims = size(tm_data_b1, /DIMENSIONS)
  
  rgb_data = intarr(dims[0],dims[1],6)
  
  rgb_data[*,*,0] = tm_data_b1
  rgb_data[*,*,1] = tm_data_b2
  rgb_data[*,*,2] = tm_data_b3
  rgb_data[*,*,3] = tm_data_b4
  rgb_data[*,*,4] = tm_data_b5
  rgb_data[*,*,5] = tm_data_b6

  rgb_strech = intarr(dims[0],dims[1],6)

  img_strech_b1=nonliner_strech(rgb_data[*,*,0])
  img_strech_b2=nonliner_strech(rgb_data[*,*,1])
  img_strech_b3=nonliner_strech(rgb_data[*,*,2])
  img_strech_b4=nonliner_strech(rgb_data[*,*,3])
  img_strech_b5=nonliner_strech(rgb_data[*,*,4])
  img_strech_b6=nonliner_strech(rgb_data[*,*,5])
  
  pan_data_1=nonliner_strech(pan_data)

  rgb_strech[*,*,0] = img_strech_b1
  rgb_strech[*,*,1] = img_strech_b2
  rgb_strech[*,*,2] = img_strech_b3
  rgb_strech[*,*,3] = img_strech_b4
  rgb_strech[*,*,4] = img_strech_b5
  rgb_strech[*,*,5] = img_strech_b6

  X=intarr(6,dims[0]*dims[1])
  X[0,*]=reform(rgb_strech[*,*,0],1,dims[0]*dims[1])
  X[1,*]=reform(rgb_strech[*,*,1],1,dims[0]*dims[1])
  X[2,*]=reform(rgb_strech[*,*,2],1,dims[0]*dims[1])
  X[3,*]=reform(rgb_strech[*,*,3],1,dims[0]*dims[1])
  X[4,*]=reform(rgb_strech[*,*,4],1,dims[0]*dims[1])
  X[5,*]=reform(rgb_strech[*,*,5],1,dims[0]*dims[1])
  
  image_PCA=PCA(rgb_data)
  pca1=round(image_PCA[*,*,0])
  
  match_pan=hist_match(pan_data_1,pca1)
;  match_pan=match_pan
  
  PCA=fltarr(dims[0],dims[1],6)
  PCA[*,*,0]=match_pan
  PCA[*,*,1]=image_PCA[*,*,1]
  PCA[*,*,2]=image_PCA[*,*,2]
  PCA[*,*,3]=image_PCA[*,*,3]
  PCA[*,*,4]=image_PCA[*,*,4]
  PCA[*,*,5]=image_PCA[*,*,5]
  
  A=A(tm_data_b1,tm_data_b2,tm_data_b3,tm_data_b4,tm_data_b5,tm_data_b6)
  PCA_data1 = fltarr(dims[0],dims[1],6)
  X = intarr(1,6)
  for i=0,dims[0]-1 do begin
    for j=0,dims[1]-1 do begin
      x[0,*]=PCA[i,j,*]
      r=A##X
      ;      help,r
      PCA_data1[i,j,*]=r
    endfor
  endfor
  new_data = intarr(dims[0],dims[1],6)
  new_data[*,*,0]=nonliner_strech(PCA_data1[*,*,0])
  new_data[*,*,1]=nonliner_strech(PCA_data1[*,*,1])
  new_data[*,*,2]=nonliner_strech(PCA_data1[*,*,2])
  new_data[*,*,3]=nonliner_strech(PCA_data1[*,*,3])
  new_data[*,*,4]=nonliner_strech(PCA_data1[*,*,4])
  new_data[*,*,5]=nonliner_strech(PCA_data1[*,*,5])
  new_data1 = intarr(dims[0],dims[1],3)
  new_data1[*,*,0]=new_data[*,*,2]
  new_data1[*,*,1]=new_data[*,*,1]
  new_data1[*,*,2]=new_data[*,*,0]
  TV, new_data1, TRUE=3

end

Function PCA, rgb_data
  dims = size(rgb_data, /DIMENSIONS) 
  b=dblarr(6,6)
  m=dblarr(6)
  for i=0,5 do begin
    m[i]=total((rgb_data[*,*,i]))/(dims[0]*dims[1])
  endfor
  PRINT,M
  for i=0,5 do begin
    for j=0,5 do begin
      img_var=total((rgb_data[*,*,i]-m[i])*(rgb_data[*,*,j]-m[j]))/(dims[0]*dims[1])
      b[i,j]=img_var
    endfor
  endfor
  eigenvalues = EIGENQL(b, EIGENVECTORS = evecs, $
    RESIDUAL = residual)
  r=sort(eigenvalues)
  r1=reverse(r)
  A=dblarr(6,6)
  A=evecs[[r1],*]
  X = intarr(1,6)
  new_data=intarr(dims[0],dims[1],6)
  for i=0,dims[0]-1 do begin
    for j=0,dims[1]-1 do begin
      x[0,*]=rgb_data[i,j,*]
      r=A##X
      new_data[i,j,*]=r
    endfor
  endfor
  help,new_data
  return,new_data
end


Function nonliner_strech, img_data

  dims = size(img_data, /DIMENSIONS)
  hist_array = HISTOGRAM(img_data,Min=0,Max=255)

  hist_array = hist_array*1.0/(dims[0]*dims[1])

  cumul_hist = fltarr(256)
  cumul_hist[0] = hist_array[0]

  For i=1,255 do begin

    cumul_hist[i] = hist_array[i]+cumul_hist [i-1]

  Endfor

  givemin=where(cumul_hist lt 0.02)
  givemax=where(cumul_hist gt 0.98)

  per_min = max(givemin)
  per_max = min(givemax)

  img_min = min(img_data)
  img_max = max(img_data)

  img_strech=round((img_data-per_min)*(245-10)/(per_max-per_min)+10.0)* $
    ((img_data ge per_min) and (img_data le per_max)) + $
    round((img_data-img_min)*(9-0)/(per_min-img_min)+0.0)* (img_data lt per_min) +$
    round((img_data-per_max)*(255-246)/(img_max-per_max)+246.0)* (img_data gt per_max)

  return,img_strech

  ;end img_var
End

Function A, tm_data_b1,tm_data_b2,tm_data_b3,tm_data_b4,tm_data_b5,tm_data_b6
  dims = size(tm_data_b1, /DIMENSIONS)     
  rgb_data = intarr(dims[0],dims[1],6)     
  rgb_data[*,*,0] = tm_data_b1             
  rgb_data[*,*,1] = tm_data_b2             
  rgb_data[*,*,2] = tm_data_b3
  rgb_data[*,*,3] = tm_data_b4
  rgb_data[*,*,4] = tm_data_b5
  rgb_data[*,*,5] = tm_data_b6
  help,rgb_data
  b=dblarr(6,6)
  ; help,b
  m=dblarr(6)
  for i=0,5 do begin
    m[i]=total((rgb_data[*,*,i]))/(dims[0]*dims[1])
  endfor
  PRINT,M
  for i=0,5 do begin
    for j=0,5 do begin
      img_var=total((rgb_data[*,*,i]-m[i])*(rgb_data[*,*,j]-m[j]))/(dims[0]*dims[1])
      b[i,j]=img_var
    endfor
  endfor
  eigenvalues = EIGENQL(b, EIGENVECTORS = evecs, $
    RESIDUAL = residual)
  r=sort(eigenvalues)
  r1=reverse(r)
  A=dblarr(6,6)
  A=evecs[[r1],*]
  return,invert(A)
end

Function hist_match, Slave_img, Base_img
  
  dimsS = size(Slave_img, /DIMENSIONS)
  dimsB = size(Base_img, /DIMENSIONS)

  ;**1. stastic histogram of two img
  hist_Slave = HISTOGRAM(Slave_img,Min=0,Max=255)
  hist_Base = HISTOGRAM(Base_img,Min=0,Max=255)

  ; convert to percent
  hist_Slave = temporary(hist_Slave)*1.0/(dimsS[0]*dimsS[1])
  hist_Base = temporary(hist_Base)*1.0/(dimsB[0]*dimsB[1])

  ;**2. stastic cumulative histogram of two img
  cumul_hist_slave = fltarr(256)
  cumul_hist_Base = fltarr(256)

  cumul_hist_slave[0] = hist_Slave[0]
  cumul_hist_Base[0] = hist_Base[0]

  For i=1,255 do begin

    cumul_hist_slave[i] = hist_Slave[i]+cumul_hist_slave [i-1]
    cumul_hist_Base[i] = hist_Base[i]+cumul_hist_Base [i-1]

  Endfor

  ;**3. image matching

  SeMa = intarr(256)
  For i=0,255 do begin

    SeDelta = abs(cumul_hist_slave[i]-cumul_hist_Base)
    Semin=min(SeDelta)
    SeMai=where(SeDelta eq Semin)

    ; print,SeMai

    SeMa[i]=SeMai[0]

  Endfor

  SeDelta =!NULL

  match_img=intarr(dimsS[0],dimsS[1])
  For i=0, dimsS[0]-1 do begin
    For j=0, dimsS[1]-1 do begin

      match_img[i,j] = SeMa(Slave_img[i,j])

    Endfor
  Endfor

  return, match_img
  SeMa = !NULL

  ;end img_var
End