矩匹配算法影像去噪

2020-5-2-9:19

常规的矩匹配算法

能起到去除条带噪声的效果，但去噪后影像灰度信息损失过大。
IDL关键代码：

  dim = size(data, /dimensions)
  row = dim[1]
  col = dim[0]

  re_mean = mean(mean(data, dimension=1))
  re_stev = mean(stddev(data, dimension=1))
  
  or_mean = mean(data, dimension=1)
  or_stev = stddev(data, dimension=1)
  
  datar = make_array(col, row, type=4)
  for i=0,col-1 do begin
    for j=0,row-1 do begin
      datar[i,j] = data[i,j]*re_stev / or_stev[j] + re_mean - $
                   or_mean[j]*re_stev / or_stev[j]
    endfor
  endfor

基于均值补偿的矩匹配算法

2020-5-3-16:34
CCD = 40
效果不好，甚至不如常规的矩匹配，可能是因为影像太大了，参考均值和方差不太合适。

  dim = size(data, /dimensions)
  row = dim[1]
  col = dim[0]
  print,row,col
  
  makeup = make_array(1, row, type=4)
  for i=0,row-41 do begin
      if i le 39 then begin
        makeup[0,i] = (total(data[*,i+40])-total(data[*,i])) / col
      endif else begin
        makeup[0,i] = (total(data[*,i+40])-total(data[*,i mod 40])) / col
      endelse
  endfor
  
  print,makeup
  makeup = makeup[*,0:7959]
  help,makeup
  
  re_mean = mean(mean(data, dimension=1))
  re_stev = mean(stddev(data, dimension=1))
  
  ;re_mean = mean(data[*,0])
  ;re_stev = stddev(data[*,0])
  
  or_mean = mean(data, dimension=1)
  or_stev = stddev(data, dimension=1)
  
  datar = make_array(col, row, type=4)
  for i=0,col-1 do begin
    for j=0,row-1 do begin
      case 1 of
       j gt 39 : datar[i,j] = data[i,j]*re_stev / or_stev[j] + re_mean - $
                   or_mean[j]*re_stev / or_stev[j] + makeup[0,j-40]
                   
       else    : datar[i,j] = data[i,j]*re_stev / or_stev[j] + re_mean - $
                   or_mean[j]*re_stev / or_stev[j]                     
      endcase
    endfor
  endfor

2020-5-5-17:25
情况有变，听从我导，暂时停更这个原始数据去噪。

基于影像分割的矩匹配算法

2020-5-12-17:05更新
以40*40的窗口为单位进行去噪，从剖面曲线来看，效果不错。但地物细节模糊严重，窗口间过渡不够平滑，失去原先影像分辨率优势。

  re_raster = e.openraster(re_file)
  or_raster = e.openraster(or_file)
  re_data = re_raster.getdata(bands=0)
  or_data = or_raster.getdata(bands=0)

  dim = size(re_data, /dimensions)
  row = dim[1]
  col = dim[0]

  data = make_array(col, row, type=4)
  subs = block(40, col, row) ;调用自定义分块影像函数,博客中有函数代码
  dim = size(subs, /dimensions)
  count = dim[1] ; 获取分块图像的数量
  
  ;or_mean = mean(or_data, dimension=1)
  ;or_stev = stddev(or_data, dimension=1)
  
  for k = 0,count-1 do begin    
    sub = subs[*,k]  ; 获取第一个分块影像的 左上角 右下角 行 列 坐标
    re_mean = mean(mean(re_data[sub[0]:sub[2], sub[1]:sub[3]], dimension=1))
    re_stev = mean(stddev(re_data[sub[0]:sub[2], sub[1]:sub[3]], dimension=1))  
            
    for i = sub[0],sub[2] do begin
      for j = sub[1],sub[3] do begin 
          or_mean = mean(or_data[sub[0]:sub[2], j]) 
          or_stev = stddev(or_data[sub[0]:sub[2], j])      
          data[i,j] = $
                    or_data[i,j]*re_stev / or_stev[j] + re_mean - $
                    or_mean[j]*re_stev / or_stev[j]                        
      endfor
    endfor    
  endfor

基于参考影像的矩匹配算法

2020-5-14-10:25更新
以没有噪声影像行均值和方差为基准逐行纠正影像，对比较明显的条带噪声效果显著，总体效果一般。

  for i = 0,row-1 do begin
    re_mean = mean(re_data[*,i])
    re_stev = stddev(re_data[*,i])
    or_mean = mean(or_data[*,i])
    or_stev = stddev(or_data[*,i])
    for j = 0,col-1 do begin
       data[j,i] = or_data[j,i]*re_stev / or_stev + re_mean - $
                   or_mean*re_stev / or_stev
    endfor
  endfor