I have two images showing the same scenery from different angles. What I want to find out is the pixel-offset of every(!) single pixel using correlation, but correlate and c_correlate only correlates two arrays and not a single value and an array. Is there something I can do?
but somehow it didn't really help. I figured I have to use a template consisting of my pixel and the neighboring pixelsto find the offset. But something's wrong.
I tried:
CONV= where(max(CONCOVAR(image1, image2, /Correl)))
where concovar is a function returning the correlationmatrix of the images, but somehow all my offsets are 0. That can't be (I checked on the images and there must be an offset around 1 or 2 pixels).
On the other hand I tried using correl_optimize (I found it on the internet), but those results are even weirder. My template is an array of 5x3 elements and the other image is an array of 11x3 elements. correl_optimize returns a y-offset between -2,5 and 2,5 pixels. But how can the y-offset possibly be in a range of 5 pixels when both images only consist of three pixels in y direction? Same with x direction.
I don't have a clue as to what I do wrong. Has someone an explanation for this?
This sounds like an Optical Flow problem. Do a google scholar search. I know there are a bunch of freeware matlab code packages on the web. I haven't seen OF packages in IDL. =(
On Friday, November 9, 2012 2:53:26 AM UTC-5, haik...@gmail.com wrote:
> Thanks Klemen,
> but somehow it didn't really help. I figured I have to use a template consisting of my pixel and the neighboring pixelsto find the offset. But something's wrong.
> where concovar is a function returning the correlationmatrix of the images, but somehow all my offsets are 0. That can't be (I checked on the images and there must be an offset around 1 or 2 pixels).
> On the other hand I tried using correl_optimize (I found it on the internet), but those results are even weirder. My template is an array of 5x3 elements and the other image is an array of 11x3 elements. correl_optimize returns a y-offset between -2,5 and 2,5 pixels. But how can the y-offset possibly be in a range of 5 pixels when both images only consist of three pixels in y direction? Same with x direction.
> I don't have a clue as to what I do wrong. Has someone an explanation for this?
Hi Max, my code compares one image with two other images and it is not really well documented :( But if you want some exaple here is a part of my old code that is more simple, without pyramids. I am not sure if it functions as it is, as I just tried to simplifiy it...
Cheers, Klemen
PRO cloud_height_correlation, file, modis_path, seviri_path
;Cross-corealtaion settings
d_search_area = 3L ;width of the half of the sqaure defining the search area
d_window = 3 ;width of the half of the search window (sqaure)
d_mask = d_search_area - d_window
max_pixels = 10L*(10L)^5 ;maximal number of pixels to be processed at the same time - avoid to run out of memory
;Seach max correlation for all OK pixels
shift_s1_col = long(mask*0)
shift_s1_lin = shift_s1_col
cross_cor1 = make_array(in_size[1], in_size[2], value=-1.)
;prepare the indexes for the moving window
tmp = lindgen(in_size[1], in_size[2])
tmp = tmp - tmp[d_window,d_window]
indx_mw = tmp[0:2*d_window, 0:2*d_window]
tmp = !null
indx_mw = rebin(reform(indx_mw, 1, (2*d_window+1)^2), max_pixels, (2*d_window+1)^2)
;because of the memory issues do not run everything at once
for i=0L,count_mask-1,max_pixels do begin
print, i
if (i+max_pixels) gt (count_mask-1) then begin
max_pixels = count_mask - i
indx_mw = indx_mw[0:max_pixels-1, *]
endif
;initialize partial results
tmp_cor1 = make_array(max_pixels, value=-1.) ;partial output cross-correlation
tmp_cor2 = tmp_cor1
tmp_s1_col = make_array(max_pixels) ;partial output shift for SEVIRI 1 columns
tmp_s1_lin = tmp_s1_col
;select indexes
indx_i = indx_mask[i:i+max_pixels-1] ;indexes of pixels to be processed
indx_all = rebin(indx_i, max_pixels, (2*d_window+1)^2) + indx_mw ;indexes of these pixels and their vicinities
;prepare cross-correlation
y = im_mo[indx_all] ;MODIS with its vicinity
ym = total(y, 2) / (2*d_window+1)^2 ;average of MODIS within the pixel vicinity
ym = rebin(ym, max_pixels, (2*d_window+1)^2)
sum_y2 = total((y - ym)^2, 2)
sign_col = 1
for sh_lin=0,2 do BEGIN
for sh_col_sign=0,5 do BEGIN sign_col = sign_col * (-1)
sh_col = sh_col_sign / 2 * sign_col
print, sh_lin, sh_col
;SEVIRI 1
T = systime(1)
im_shift = shift(im_s1, sh_col, sh_lin))
x = im_shift[indx_all]
xm = total(x, 2) / (2*d_window+1)^2 ;average of SEVIRI within the pixel vicinity
xm = rebin(xm, max_pixels, (2*d_window+1)^2)
T = systime(1)
sum_x2 = total((x - xm)^2, 2)
T = systime(1)
sum_xy = total((x - xm)*(y-ym), 2)
tmp = sum_xy / sqrt(sum_x2 * sum_y2) ;cross correlation
indx_max = where(tmp gt tmp_cor1, count_max)
if count_max gt 0 then begin
tmp_s1_col[indx_max] = sh_col
tmp_s1_lin[indx_max] = sh_lin
tmp_cor1[indx_max] = tmp[indx_max]
endif
Well, I'm quite new to IDL and that code still looks very complex to me (in fact I didn't expect it to be so complicated^^). I guess in your code x and y are the pixel offsets, is that correct?
Computing correlation using IDL function CORRELATE is ok, if you have to do it once. But here you have to do it within each area of interest a couple of times (considering set of possible offsets in x and y direction) so you have to use 4 FOR loops which will run slow!
Thus you should ask yourself first, how large is the data set you want to process. If you are sure that your offset are not larger than 2 pixels, then you might go for the slow version. That would be something like:
results_corr = make_array(...
results_offx = make_array(...
results_offy = make_array(...
FOR xall=xstart,xend do begin
FOR yall=ystart,yend do begin
search_area = data[xall-search_size:yall+search_size]
template = data[xall-template_size:yall+template_size]
FOR xoff=-search_size:search_size do begin
FOR yoff=-search_size:search_size do begin
data = (shift(search_area, xoff, yoff))[right indeice]
corr = correlate(data[*], template[*])
if corr gt results_corr[xall, yall] then begin
results_corr[xall, yall] = corr
results_offx[xall, yall] = xoff
results_offy[xall, yall] = yoff
endif
endfor
endfor
endfor
endfor
