PRO correlate_images, all, result
  d = size(all,/dim)
  VX = d[1]
  VY = d[2]
  ; normalize the background
  for i = 0, d[0] - 1 do begin
    all[i,*,*] /= mean(all[i,*,*])
  endfor
  ; Rho correlation
  cor_pic = make_array(VX,VY,/double,value=0.0)
  f_pic = make_array(VX,VY,value=0.0)
  for x = 0, VX - 1 do begin
    for y = 0, VY - 1 do begin
      r = r_correlate(reform(all[*,x,y]),result)
      cor_pic[x,y]=r[0]
      f_pic[x,y]=1.0-r[1]
    endfor
  endfor
  ; we are interested in correlations with high variance on gel
  var_pic = make_array(VX,VY,/double,value=0.0)
  for x = 0, VX - 1 do begin
    for y = 0, VY - 1 do begin
      var_pic[x,y]=stddev(all[*,x,y])
    endfor
  endfor
  var_pic <= 1.0
  f_pic *= var_pic
  cor_pic <= 1.0
  cor_pic >= -1.0
  show_correlation, cor_pic, f_pic
end

function gauss2d, sx, sy, cx, cy, wx, wy, a
  im = float(make_array(sx,sy,value=0.0))
  for x = 0, sx - 1 do begin
    for y = 0, sy - 1 do begin
      im[x,y]=float(((cx-x)/wx)^2 + ((cy-y)/wy)^2)
    endfor
  endfor
  im = -im/2
  im = exp(im)
  im *= a
  return, im
end

function create_set, nr, sx, sy, wx1, wx2, wy1, wy2, a1, a2
  all = make_array(nr,sx,sy,value=0.0)
  for i = 0, nr - 1 do begin
    wx = wx1 + i*(wx2-wx1)/nr
    wy = wy1 + i*(wy2-wy1)/nr
    a = a1 + i*(a2-a1)/nr
    all[i,*,*] = gauss2d(sx, sy, sx/2, sy/2, wx, wy, a)
  endfor
  return, all
end

set1 = create_set(15, 600.0, 400.0, 10.0, 100.0, 10.0, 100.0, 5.0, 1.0)
set2 = create_set(15, 300.0, 300.0, 10.0, 40.0, 40.0, 10.0, 5.0, 5.0)
set3 = create_set(15, 300.0, 300.0, 20.0, 20.0, 20.0, 20.0, 5.0, 5.0)
set4 = create_set(15, 300.0, 300.0, 20.0, 20.0, 20.0, 20.0, 1.0, 1.0)

set1[*,0:299,0:299] += set2[*,*,*]
set1[*,300:599,0:299] += set3[*,*,*]
set1[*,300:599,100:399] += set4[*,*,*]
for i = 0, 14 do begin
  set1[i,0+i*10:299+i*10,200:399]+=set3[i,*,50:249]
endfor
for i = 0, 14 do begin
  set1[i,*,*]/=double(i)
endfor
set1 = relative(set1)
result1 = findgen(15)
correlate_images, set1, result1
end

PRO show_correlation, cp, t_pic
  cor_pic = cp
  DDD = size(cp,/dim)
  VX = ddd[0]
  VY = ddd[1]
  ; the normal one
  shown = make_array(3,VX,VY,/double,value=255.0)
  multi = 1.0 / max(abs(cor_pic))
  multi = 1.0 / max(abs(cor_pic))
  shown[0,*,*] += (cor_pic[*,*] < 0) * multi * 55
  shown[1,*,*] += (cor_pic[*,*] < 0) * multi * 155
  shown[2,*,*] += (cor_pic[*,*] < 0) * multi * 255
  shown[0,*,*] -= (cor_pic[*,*] > 0) * multi * 255
  shown[1,*,*] -= (cor_pic[*,*] > 0) * multi * 55
  shown[2,*,*] -= (cor_pic[*,*] > 0) * multi * 205
  window, 1, title='Correlation', ret=2, xsize=vx, ysize=vy
  shown >= 0
  shown <= 255
  tvscl, shown, /true
  ; only the significant correlation
  wcor_pic = double(cor_pic) * double(t_pic)
  shown = make_array(3,VX,VY,/double,value=255.0)
  multi = 255.0 / max(abs(wcor_pic))
  shown[0,*,*] += (wcor_pic[*,*] < 0) * multi * 55
  shown[1,*,*] += (wcor_pic[*,*] < 0) * multi * 155
  shown[2,*,*] += (wcor_pic[*,*] < 0) * multi * 255
  shown[0,*,*] -= (wcor_pic[*,*] > 0) * multi * 255
  shown[1,*,*] -= (wcor_pic[*,*] > 0) * multi * 55
  shown[2,*,*] -= (wcor_pic[*,*] > 0) * multi * 205
  window, 3, title='Significant Correlations', ret=2, xsize=vx, ysize=vy
  tvscl, shown, /true
  shown = bytscl(shown)
  profiles, cor_pic
  ; significance
  window, 2, title='Significance', ret=2, xsize=vx, ysize=vy
  tvscl, t_pic
  profiles, t_pic
end