Документ взят из кэша поисковой машины. Адрес оригинального документа : http://lnfm1.sai.msu.ru/~kupry/SOFTWARE/korel.pro Дата изменения: Fri Aug 15 13:05:25 2003 Дата индексирования: Mon Oct 1 22:39:07 2012 Кодировка: Поисковые слова: integral

pro rectangle,xx,yy,sirka,vyska ;nakresli obdelnik
;device,set_graphics=6
px=[xx,xx+sirka,xx+sirka,xx,xx]
py=[yy,yy,yy+vyska,yy+vyska,yy]
plots,px,py,/dev,color=255;lines=0,color=255
;empty
end

pro rectangle1,xx,yy,sirka,vyska
;device,set_graphics=6
px=[xx,xx+sirka,xx+sirka,xx,xx]
py=[yy,yy,yy+vyska,yy+vyska,yy]
plots,px,py,/dev,color=0;lines=0,color=255
;empty
end

function integruj,obr,x,y,w,h ;spocita integral v oblasti x:x+w,y:y+h
sub=float(obr(x:x+w,y:y+h)) & s=size(sub)
iks=findgen(s(1)) & ips=findgen(s(2))
pole=float(findgen(s(2)))
for i=0,s(2)-1 do pole(i)=int_tabulated(iks,sub(*,i))
return,int_tabulated(ips,pole)
end

pro hist2_event,event
COMMON uloz4,histdraw2,corrdraw2,corrbase2
widget_control,event.id,get_uvalue=uvalue
case uvalue of
'SAVE':begin

widget_control,histdraw2,get_value=okno2
wset,okno2 & write_jpeg,'hist1.jpg',tvrd()

if widget_info(corrbase2,/valid_id) eq 1 then begin
widget_control,corrdraw2,get_value=okno1
wset,okno1 & write_jpeg,'corr1.jpg',tvrd()
widget_control,corrbase2,/destroy
endif
widget_control,event.top,/destroy

end
'CANCEL':begin
widget_control,event.top,/destroy
if widget_info(corrbase2,/valid_id) eq 1 then widget_control,corrbase2,/destroy
end
endcase
end

pro cross_event,event
COMMON uloz3,crossdraw2
widget_control,event.id,get_uvalue=uvalue
case uvalue of
'OK':begin
widget_control,crossdraw2,get_value=okno & wset,okno
write_jpeg,'cross.jpg',tvrd()
widget_control,event.top,/destroy
end
'CANCEL':widget_control,event.top,/destroy
endcase
end

pro auto_event,event
COMMON uloz2,autodraw2
widget_control,event.id,get_uvalue=uvalue
case uvalue of
'OK':begin
widget_control,autodraw2,get_value=okno & wset,okno
write_jpeg,'auto.jpg',tvrd()
widget_control,event.top,/destroy
end
'CANCEL':widget_control,event.top,/destroy
endcase
end


pro light_event,event
COMMON uloz6,lightdraw2
widget_control,event.id,get_uvalue=uvalue
case uvalue of
'OK':begin
widget_control,lightdraw2,get_value=okno & wset,okno
write_jpeg,'light.jpg',tvrd()
widget_control,event.top,/destroy
end
'CANCEL':widget_control,event.top,/destroy
endcase
end

pro hist_event,event
COMMON uloz,corrdraw,histdraw
widget_control,event.id,get_uvalue=uvalue
case uvalue of
'OK':begin
widget_control,corrdraw,get_value=okno1
widget_control,histdraw,get_value=okno2
wset,okno1 & write_jpeg,'corr.jpg',tvrd()
wset,okno2 & write_jpeg,'hist.jpg',tvrd()
widget_control,event.top,/destroy
end
'CANCEL':widget_control,event.top,/destroy
endcase
end

pro mnoho_event,event
COMMON uloz5,mnohofield,mnohodraw,mnohodraw2
widget_control,event.id,get_uvalue=uvalue
case uvalue of
'OK':begin
widget_control,mnohodraw,get_value=okno1
widget_control,mnohodraw2,get_value=okno2
wset,okno1 & write_jpeg,'reg.jpg',tvrd()
wset,okno2 & write_jpeg,'count.jpg',tvrd()
widget_control,event.top,/destroy
end
'CANCEL':widget_control,event.top,/destroy
endcase
end


;--------zde se pocitaji integraly pro vsechny oblasti site pro kazdy obrazek
pro integrals_event,event
COMMON integral,poleintegralu,xpocatek,ypocatek,xpocet,ypocet,krok,prumer
COMMON obrazky,soubory,pocet
COMMON oblast,inttext1,inttext2,inttext3
COMMON baze,base,innerbase
widget_control,event.id,get_uvalue=uvalue
case uvalue of

'OK':begin
;-----------nacteni zadani oblasti----------------------
widget_control,inttext1,get_value=text1 & text1=text1(0)
widget_control,inttext2,get_value=text2 & text2=text2(0)
widget_control,inttext3,get_value=text3 & text3=text3(0)
pozice=strpos(text1,',') & delka=strlen(text1)
xpocatek=fix(strmid(text1,0,pozice)) & ypocatek=fix(strmid(text1,pozice+1,delka-pozice))
pozice=strpos(text2,',') & delka=strlen(text2)
xpocet=fix(strmid(text2,0,pozice)) & ypocet=fix(strmid(text2,pozice+1,delka-pozice))
krok=fix(text3)
widget_control,event.top,/destroy

;----------vytvoreni pole integralu--------------------
obrazek=readfits(soubory(pocet/2)) & s=size(obrazek)
widget_control,innerbase,/destroy
innerbase=widget_base(base,/column)
d=widget_draw(innerbase,xsize=s(1),ysize=s(2))
innerlabel=widget_label(innerbase,xsize=900,/align_left)

poleintegralu=fltarr(xpocet,ypocet,pocet)
for o=0,pocet-1 do begin ;pocet-1
obrazek=readfits(soubory(o))
widget_control,innerlabel,set_value=soubory(o)
for i=0,xpocet-1 do begin
for j=0,ypocet-1 do begin
poleintegralu(i,j,o)=integruj(obrazek,xpocatek+i*krok,ypocatek+j*krok,krok,krok)
endfor
endfor
endfor

;pro dany obrazek odecteme od integralu jejich stredni hodnotu
prumer=fltarr(pocet)
for i=0,pocet-1 do begin ;pocet-1
prumer(i)=total(poleintegralu(*,*,i))/float(xpocet*ypocet)
poleintegralu(*,*,i)=poleintegralu(*,*,i)-prumer(i)
endfor

;ulozeni integralu, strednich hodnot a parametru site a vykresleni site
save,prumer,filename='average.dat'
save,poleintegralu,filename='variable.dat'
openw,unit,'region.dat',/get_lun
printf,unit,xpocatek,ypocatek,xpocet,ypocet,krok ;ulozeni parametru oblasti
free_lun,unit

tvscl,obrazek
for i=0,xpocet-1 do begin
for j=0,ypocet-1 do rectangle,xpocatek+i*krok,ypocatek+j*krok,krok,krok
endfor
widget_control,innerlabel,set_value='The integrals have been counted'
end


'CANCEL':widget_control,event.top,/destroy

endcase
end

pro korel_event,event
COMMON uloz4,histdraw2,corrdraw2,corrbase2
COMMON baze,base,innerbase
COMMON obrazky,soubory,pocet
COMMON oblast,inttext1,inttext2,inttext3
COMMON integral,poleintegralu,xpocatek,ypocatek,xpocet,ypocet,krok,prumer
COMMON start,omenu
COMMON uloz,corrdraw,histdraw
COMMON uloz2,autodraw2
COMMON uloz3,crossdraw2
COMMON inner1,autofield
COMMON inner2,crossfield1,crossfield2
COMMON inner3,distfield
COMMON inner4,corectfield1,corectfield2,corectfield3
COMMON inner5,prepinac,prepinac1
COMMON uloz5,mnohofield,mnohodraw,mnohodraw2
COMMON inner6,lightfield
COMMON uloz6,lightdraw2
widget_control,event.id,get_uvalue=uvalue
case uvalue of

'LOAD':begin ;zde se nahrava adresar se serii obrazku
adresar=dialog_pickfile(/directory)
soubory=findfile(adresar+'*.fts',count=pocet)

;pokud jsou spocteny integraly, nahrajem je
r=findfile('variable.dat') ;integraly
if r(0) ne '' then restore,'variable.dat'
reg=findfile('region.dat') ;parametry oblasti
if reg(0) ne '' then begin
openr,unit,'region.dat',/get_lun
readf,unit,xpocatek,ypocatek,xpocet,ypocet,krok
free_lun,unit
endif
pr=findfile('average.dat') ;stredni hodnoty
if pr(0) ne '' then restore,'average.dat'
widget_control,omenu,sensitive=1
end

'QUIT':widget_control,event.top,/destroy

;-----------------------------------------------------------

'INTEGRALS':begin
intbase=widget_base(/column,title='region for integration')
intbase1=widget_base(intbase,/column,/frame)
inttext1=cw_field(intbase1,title='Start point (x,y)')
inttext2=cw_field(intbase1,title='Number of regions (n,m)')
inttext3=cw_field(intbase1,title='Step (k)')
intbase2=widget_base(intbase,/row)
okbut=widget_button(intbase2,value='OK',uvalue='OK')
canbut=widget_button(intbase2,value='CANCEL',uvalue='CANCEL')
widget_control,intbase,/realize
xmanager,'integrals',intbase
end

'KOEFICIENTY':begin ;zde se pocitaji korelacni koeficienty pro vsechny dvojice oblasti
obrazek=readfits(soubory(pocet/2)) & s=size(obrazek)
widget_control,innerbase,/destroy
innerbase=widget_base(base,/column)
d=widget_draw(innerbase,xsize=s(1),ysize=s(2))
tvscl,obrazek


;------------------ukladani korelacnich koeficientu-----------
openw,unit,'koeficienty.dat',/get_lun
for i=0,ypocet-1 do begin ;kombinace oblasti
for j=0,xpocet-1 do begin ;oblast [l,k] a [j,i]
for k=0,ypocet-1 do begin
for l=0,xpocet-1 do begin
if (k gt i) or ((k eq i) and (l gt j)) then begin
x1=xpocatek+l*krok & y1=ypocatek+k*krok
x2=xpocatek+j*krok & y2=ypocatek+i*krok
lag=0 & x=poleintegralu(l,k,*) & y=poleintegralu(j,i,*)
result=c_correlate(x,y,lag)
printf,unit,result(0),x2,y2,x1,y1
rectangle,x2,y2,krok,krok
endif
endfor
endfor
endfor
endfor
free_lun,unit
corlabel=widget_label(innerbase,value='The correlation coefficients has been saved')
end

'HISTOGRAM':begin ;zde se pocita histogram vsech korelacnich koeficientu
widget_control,innerbase,/destroy
innerbase=widget_base(base,/column)
obrazek=readfits(soubory(pocet/2))
s=size(obrazek)
corrdraw=widget_draw(innerbase,xsize=s(1),ysize=s(2))
tvscl,obrazek

;----------nacteni koeficientu ze souboru--------------
n=long(xpocet)*long(ypocet) ;nezapomenout na longint!
pocetkombinaci=n*(n-1)/2
openr,unit,'koeficienty.dat',/get_lun
data=fltarr(5,pocetkombinaci)
readf,unit,data
free_lun,unit

;---------rozdeleni a jeho parametry--------
hist=histogram(data(0,*),min=-1,max=1,bins=0.01)
hodnotyx=findgen(n_elements(hist))*0.01-1
gauss=gaussfit(hodnotyx,hist,koef,nterms=3)
mi=koef(1) & sigma=koef(2)
horni=mi+3.0*sigma & dolni=mi-3.0*sigma


;---------vykresleni oblasti s vysokymi korelacemi-------
indexy=where((data(0,*) lt dolni) or (data(0,*) gt horni),count)
korelace=data(*,indexy)

for i=0,count-1 do begin
rectangle,korelace(1,i),korelace(2,i),krok,krok ;x,y,s,v
rectangle,korelace(3,i),korelace(4,i),krok,krok
endfor

;--------zobrazeni histogramu--------------
histbase=widget_base(/column)
histdraw=widget_draw(histbase,xsize=500,ysize=400)
histbase1=widget_base(histbase,/row)
okbut=widget_button(histbase1,value='Save images',uvalue='OK')
canbut=widget_button(histbase1,value='Cancel',uvalue='CANCEL')
widget_control,histbase,/realize
xmanager,'hist',histbase

widget_control,histdraw,get_value=okno2 & wset,okno2
plot,hodnotyx,hist,xrange=[-1,1],xstyle=1,$
title='DISTRIBUTION OF CORRELATION COEFFICIENTS',$
xtitle='correlation coefficient',ytitle='density'
oplot,hodnotyx,gauss,linestyle=1
xyouts,100,303,'...........',/device
xyouts,150,300,'gauss fit',/device
xyouts,100,280,'mi='+strcompress(string(mi),/remove_all),/device
xyouts,100,265,'sigma='+strcompress(string(sigma),/remove_all),/device

end

'MNOHO':begin
widget_control,innerbase,/destroy
innerbase=widget_base(base,/row)
obrazek=readfits(soubory(pocet/2))
s=size(obrazek)
mnohodraw=widget_draw(innerbase,xsize=s(1),ysize=s(2))
mnohobase=widget_base(innerbase,/column,xsize=250)
mnohofield=cw_field(mnohobase,title='n ')
mnohobutton=widget_button(mnohobase,value='SHOW REGIONS',uvalue='MNOHO2')
end

'MNOHO2':begin ;zde se zjistuje s kolika dalsimi oblastmi oblasti koreluji
widget_control,mnohofield,get_value=kolik & kolik=fix(kolik(0))
widget_control,mnohodraw,get_value=okno & wset,okno
erase

;------nacteni korelacnich koeficientu ze souboru-------
n=long(xpocet)*long(ypocet)
pocetkombinaci=n*(n-1)/2
openr,unit,'koeficienty.dat',/get_lun
data=fltarr(5,pocetkombinaci)
readf,unit,data ;zde jsou vsechny korelacni koeficienty
free_lun,unit

;-----spocteni histogramu a zjisteni mi a sigma---------
hist=histogram(data(0,*),min=-1,max=1,bins=0.01)
hodnotyx=findgen(n_elements(hist))*0.01-1
gauss=gaussfit(hodnotyx,hist,koef,nterms=3)
mi=koef(1) & sigma=koef(2)
horni=mi+3.0*sigma & dolni=mi-3.0*sigma

;------z pole data vytrhneme vysoke korelace-----------
indexy=where((data(0,*) lt dolni) or (data(0,*) gt horni),count)
korelace=data(*,indexy) ;pole s vysokymi korelacemi

;----vytvorime pole oblasti, ve kterem budou pocty oblasti, s kterymi
;----prislusna oblast koreluje---------------------
oblasti=intarr(xpocet,ypocet)
for i=0,count-1 do begin
x1=korelace(1,i) & y1=korelace(2,i)
x2=korelace(3,i) & y2=korelace(4,i)
s1=(x1-xpocatek)/krok & s2=(y1-ypocatek)/krok
oblasti(s1,s2)=oblasti(s1,s2)+1
s1=(x2-xpocatek)/krok & s2=(y2-ypocatek)/krok
oblasti(s1,s2)=oblasti(s1,s2)+1
endfor

;----vykresleni oblasti, ktere koreluji s vice nez n oblastmi
obrazek=readfits(soubory(pocet/2))
tvscl,obrazek
for i=0,xpocet-1 do begin
for j=0,ypocet-1 do rectangle,xpocatek+i*krok,ypocatek+j*krok,krok,krok
endfor
for i=0,xpocet-1 do xyouts,xpocatek+i*krok-5,ypocatek-10,string(i+1),/device,charsize=0.8,alignment=0.5
for i=0,ypocet-1 do xyouts,xpocatek-20,ypocatek+i*krok+5,string(i+1),/device,charsize=0.8,alignment=0.5
indexy1=where(oblasti gt kolik,count1)
if count1 ne 0 then begin
print,'souradnice oblasti:'
for i=0,count1-1 do begin
s1=fix(indexy1(i) mod xpocet)
s2=fix(indexy1(i)/xpocet)
rectangle1,s1*krok+xpocatek,s2*krok+ypocatek,krok,krok
;print,s1+1,s2+1
endfor
endif

;---------pocty oblasti a jejich contour------------------
mnohobase2=widget_base(/column)
mnohodraw2=widget_draw(mnohobase2,xsize=400,ysize=600)
mnohobase21=widget_base(mnohobase2,/row)
mnohobut1=widget_button(mnohobase21,value='SAVE',uvalue='OK')
mnohobut2=widget_button(mnohobase21,value='CANCEL',uvalue='CANCEL')
widget_control,mnohobase2,/realize
xmanager,'mnoho',mnohobase2

widget_control,mnohodraw2,get_value=okno & wset,okno
obr=obrazek(xpocatek:xpocatek+xpocet*krok,ypocatek:ypocatek+ypocet*krok)
obr2=obr
;-----------------pocty oblasti-------------
mez=50
konst=fix(float(mez*max(oblasti))/float(255-mez))
for m=0,xpocet-1 do begin
for n=0,ypocet-1 do begin
xx=m*krok & yy=n*krok
jas=fix(255*float(oblasti(m,n)+konst)/float(max(oblasti)+konst))
if oblasti(m,n) eq 0 then jas=0
obr2(xx:xx+krok,yy:yy+krok)=jas ;zde oblasti priradime jas podle poctu
endfor ;oblasti s kterymi koreluje
endfor
tv,congrid(obr2,370,270),15,315
for m=0,xpocet-1 do begin
xkrok=370.0/xpocet & ykrok=270.0/ypocet
for n=0,ypocet-1 do rectangle,m*xkrok+15,n*ykrok+315,xkrok,ykrok
endfor
xyouts,200,590,'COUNT OF CORRELATIONS',charsize=1,alignment=0.5,/device
;----------------contour------------------
tvscl,congrid(obr,370,270),15,15
contour,oblasti,charsize=0.7,title='CONTOUR',pos=[15,15,385,285],/noerase,xmargin=[15,15],$
ymargin=[15,15],xrange=[-0.5,xpocet-1.5],yrange=[0.5,ypocet-1.5],xstyle=1,ystyle=1,$
xticks=xpocet,yticks=ypocet,xtickn=replicate(' ',xpocet+1),ytickn=replicate(' ',ypocet+1),/device

end

'AUTO':begin
;----------zobrazeni obrazku---------------
widget_control,innerbase,/destroy
innerbase=widget_base(base,/row)
obrazek=readfits(soubory(pocet/2))
s=size(obrazek)
autodraw=widget_draw(innerbase,xsize=s(1),ysize=s(2))
tvscl,obrazek
;-------------vykresleni site-----------
for i=0,xpocet-1 do begin
for j=0,ypocet-1 do rectangle,xpocatek+i*krok,ypocatek+j*krok,krok,krok
endfor
for i=0,xpocet-1 do xyouts,xpocatek+i*krok-5,ypocatek-10,string(i+1),/device,charsize=0.8,alignment=0.5
for i=0,ypocet-1 do xyouts,xpocatek-20,ypocatek+i*krok+5,string(i+1),/device,charsize=0.8,alignment=0.5


autobase=widget_base(innerbase,/column,xsize=250)
autofield=cw_field(autobase,title='coordinates of region (x,y)')
autobut=widget_button(autobase,value='CORRELATE',uvalue='AUTOCORR')
end
'AUTOCORR':begin ;zde se pocita autokorelace vybrane oblasti a jeji FFT
widget_control,autofield,get_value=sourad & sourad=sourad(0)
pozice=strpos(sourad,',') & delka=strlen(sourad)
xsourad=fix(strmid(sourad,0,pozice)) & ysourad=fix(strmid(sourad,pozice+1,delka-pozice))

;----graf s korelaci--------
autobase2=widget_base(/column)
autodraw2=widget_draw(autobase2,xsize=400,ysize=600)
autobase3=widget_base(autobase2,/row)
autobut2=widget_button(autobase3,value='SAVE',uvalue='OK')
autocan=widget_button(autobase3,value='CANCEL',uvalue='CANCEL')
widget_control,autobase2,/realize
xmanager,'auto',autobase2

widget_control,autodraw2,get_value=okno & wset,okno
rada=poleintegralu(xsourad-1,ysourad-1,*)
lag=findgen(2*pocet-1)-pocet+1
res=a_correlate(rada,lag)
koef=res(pocet-1) & koef=strcompress(string(koef),/remove_all)

stary=!p.multi & !p.multi=[0,1,2]
plot,lag,res,title='AUTO CORRELATION c0='+koef,xtitle='lag',ytitle='correlation coefficient',$
yrange=[-1,1],ystyle=1
plot,/ylog,abs(fft(res,-1)),title='THE LOG OF THE POWER SPECTRUM'
!p.multi=stary
end

'CROSS':begin
widget_control,innerbase,/destroy
innerbase=widget_base(base,/row)
obrazek=readfits(soubory(pocet/2))
s=size(obrazek)
crossdraw=widget_draw(innerbase,xsize=s(1),ysize=s(2))
tvscl,obrazek
;------------zobrazeni site---------------
for i=0,xpocet-1 do begin
for j=0,ypocet-1 do rectangle,xpocatek+i*krok,ypocatek+j*krok,krok,krok
endfor
for i=0,xpocet-1 do xyouts,xpocatek+i*krok-5,ypocatek-10,string(i+1),/device,charsize=0.8,alignment=0.5
for i=0,ypocet-1 do xyouts,xpocatek-20,ypocatek+i*krok+5,string(i+1),/device,charsize=0.8,alignment=0.5


crossbase=widget_base(innerbase,/column,xsize=250)
crossfield1=cw_field(crossbase,title='coordinates of 1st region (x,y)')
crossfield2=cw_field(crossbase,title='coordinates of 2nd region (x,y)')
crossbut=widget_button(crossbase,value='CORRELATE',uvalue='CROSSCORR')
end



'CROSSCORR':begin ;zde se pocita korelace dvou vybranych oblasti a jeji FFT
widget_control,crossfield1,get_value=sourad & sourad=sourad(0)
pozice=strpos(sourad,',') & delka=strlen(sourad)
xsourad1=fix(strmid(sourad,0,pozice)) & ysourad1=fix(strmid(sourad,pozice+1,delka-pozice))
widget_control,crossfield2,get_value=sourad & sourad=sourad(0)
pozice=strpos(sourad,',') & delka=strlen(sourad)
xsourad2=fix(strmid(sourad,0,pozice)) & ysourad2=fix(strmid(sourad,pozice+1,delka-pozice))

;----graf s korelaci--------
crossbase2=widget_base(/column)
crossdraw2=widget_draw(crossbase2,xsize=400,ysize=600)
crossbase3=widget_base(crossbase2,/row)
crossbut2=widget_button(crossbase3,value='SAVE',uvalue='OK')
crosscan=widget_button(crossbase3,value='CANCEL',uvalue='CANCEL')
widget_control,crossbase2,/realize
xmanager,'cross',crossbase2

widget_control,crossdraw2,get_value=okno & wset,okno

rada1=poleintegralu(xsourad1-1,ysourad1-1,*)
rada2=poleintegralu(xsourad2-1,ysourad2-1,*)
lag=findgen(2*pocet-1)-pocet+1
res=c_correlate(rada1,rada2,lag)
koef=res(pocet-1) & koef=strcompress(string(koef),/remove_all)

stary=!p.multi & !p.multi=[0,1,2]
plot,lag,res,title='CROSS CORRELATION c0='+koef,xtitle='lag',ytitle='correlation coefficient',$
yrange=[-1,1],ystyle=1
plot,/ylog,abs(fft(res,-1)),title='THE LOG OF THE POWER SPECTRUM'
!p.multi=stary
end

'LIGHT':begin
widget_control,innerbase,/destroy
innerbase=widget_base(base,/row)
obrazek=readfits(soubory(pocet/2))
s=size(obrazek)
lightdraw=widget_draw(innerbase,xsize=s(1),ysize=s(2))
tvscl,obrazek
;-----------zobrazeni site---------------
for i=0,xpocet-1 do begin
for j=0,ypocet-1 do rectangle,xpocatek+i*krok,ypocatek+j*krok,krok,krok
endfor
for i=0,xpocet-1 do xyouts,xpocatek+i*krok-5,ypocatek-10,string(i+1),/device,charsize=0.8,alignment=0.5
for i=0,ypocet-1 do xyouts,xpocatek-20,ypocatek+i*krok+5,string(i+1),/device,charsize=0.8,alignment=0.5


lightbase=widget_base(innerbase,/column,xsize=250)
lightfield=cw_field(lightbase,title='coordinates of region (x,y)')
lightbut=widget_button(lightbase,value='SHOW LIGHT CURVE',uvalue='LIGHT2')
end

'LIGHT2':begin ;zde se zobrazi svetelna krivka vybrane oblasti a jeji FFT
widget_control,lightfield,get_value=sourad & sourad=sourad(0)
pozice=strpos(sourad,',') & delka=strlen(sourad)
xsourad=fix(strmid(sourad,0,pozice)) & ysourad=fix(strmid(sourad,pozice+1,delka-pozice))

;----graf s se svetelnou krivkou a FFT--------
lightbase2=widget_base(/column)
lightdraw2=widget_draw(lightbase2,xsize=700,ysize=600)
lightbase3=widget_base(lightbase2,/row)
lightbut2=widget_button(lightbase3,value='SAVE',uvalue='OK')
lightcan=widget_button(lightbase3,value='CANCEL',uvalue='CANCEL')
widget_control,lightbase2,/realize
xmanager,'light',lightbase2

widget_control,lightdraw2,get_value=okno & wset,okno
rada=poleintegralu(xsourad-1,ysourad-1,*) ;rada po odecteni str. hodnoty
radapuv=rada+prumer ;rada puvodni - pred odectenim str. hodnoty integralu

stary=!p.multi & !p.multi=[0,2,2]
plot,rada,title='LIGHT CURVE',charsize=1,xtitle='after substracting the average value',$
yrange=[min(rada)-500,max(rada)+500],ystyle=1
plot,/ylog,abs(fft(rada,-1)),title='THE LOG OF THE POWER SPECTRUM',charsize=1
plot,radapuv,title='LIGHT CURVE',charsize=1,xtitle='before substracting the average value',$
yrange=[min(radapuv)-500,max(radapuv)+500],ystyle=1
plot,/ylog,abs(fft(radapuv,-1)),title='THE LOG OF THE POWER SPECTRUM',charsize=1

!p.multi=stary
end

'DIST':begin
widget_control,innerbase,/destroy
innerbase=widget_base(base,/row)
obrazek=readfits(soubory(pocet/2))
s=size(obrazek)
crossdraw=widget_draw(innerbase,xsize=s(1),ysize=s(2))
tvscl,obrazek
;-----------zobrazeni site---------------
for i=0,xpocet-1 do begin
for j=0,ypocet-1 do rectangle,xpocatek+i*krok,ypocatek+j*krok,krok,krok
endfor
for i=0,xpocet-1 do xyouts,xpocatek+i*krok-5,ypocatek-10,string(i+1),/device,charsize=0.8,alignment=0.5
for i=0,ypocet-1 do xyouts,xpocatek-20,ypocatek+i*krok+5,string(i+1),/device,charsize=0.8,alignment=0.5

distbase=widget_base(innerbase,/column,xsize=250)
distfield=cw_field(distbase,title='coordinates of region (x,y)')
togbase1=widget_base(distbase,/column,/exclusive)
prep1=widget_button(togbase1,value='histogram from all the coefficients',uvalue='togg1',/no_release)
prep2=widget_button(togbase1,value='histogram only for chosen region',uvalue='togg2',/no_release)
distbut=widget_button(distbase,value='SHOW CORRELATION',uvalue='DISTCORR')
widget_control,togbase1,/realize
widget_control,prep1,/set_button
prepinac1=1
end

'togg1':prepinac1=1
'togg2':prepinac1=2


'DISTCORR':begin ;zde se zjistuje s kterymi oblastmi zadana oblast koreluje
obrazek=readfits(soubory(pocet/2))
s=size(obrazek)
corrbase2=widget_base(/column)
corrdraw2=widget_draw(corrbase2,xsize=s(1),ysize=s(2))
widget_control,corrbase2,/realize

histbase2=widget_base(/column)
histdraw2=widget_draw(histbase2,xsize=500,ysize=400)
histbase21=widget_base(histbase2,/row)
histbut2=widget_button(histbase21,value='SAVE IMAGES',uvalue='SAVE')
histcan2=widget_button(histbase21,value='CANCEL',uvalue='CANCEL')
widget_control,histbase2,/realize
xmanager,'hist2',histbase2

;----------zadane hodnoty-----------------
widget_control,distfield,get_value=sourad & sourad=sourad(0)
pozice=strpos(sourad,',') & delka=strlen(sourad)
xsourad=fix(strmid(sourad,0,pozice)) & ysourad=fix(strmid(sourad,pozice+1,delka-pozice))
x=(xsourad-1)*krok+xpocatek & y=(ysourad-1)*krok+ypocatek

;----------nacteni korelacnich koeficientu-----------------
n=long(xpocet)*long(ypocet) ;nezapomenout na longint!!
pocetkombinaci=n*(n-1)/2
openr,unit,'koeficienty.dat',/get_lun
data=fltarr(5,pocetkombinaci)
readf,unit,data
free_lun,unit

;---------vytrhnuti koeficientu tykajicich se zadane oblasti----------
indexy=where(((data(1,*) eq x) and (data(2,*) eq y)) or $
((data(3,*) eq x) and (data(4,*) eq y)),count)

if prepinac1 eq 1 then korelaceh=data else korelaceh=data(*,indexy)
korelace=data(*,indexy) ;pokud prepinc1=1, pak histogram ze vsech koeficientu
;jinak pouze z koef. tykajicich se nasi oblasti

;-----------histogram a jeho parametry-----------------------;
hist=histogram(korelaceh(0,*),min=-1,max=1,bins=0.01)
hodnotyx=findgen(n_elements(hist))*0.01-1
gauss=gaussfit(hodnotyx,hist,koef,nterms=3)
mi=koef(1) & sigma=koef(2)
horni=mi+3.0*sigma & dolni=mi-3.0*sigma ;


;----------nakresleni korelaci------------------
widget_control,corrdraw2,get_value=okno1 & wset,okno1
indexy1=where((korelace(0,*) lt dolni) or (korelace(0,*) gt horni),count)
tvscl,obrazek
if count ne 0 then begin
korelace1=korelace(*,indexy1)
for i=0,count-1 do begin
rectangle,korelace1(1,i),korelace1(2,i),krok,krok ;x,y,s,v
rectangle,korelace1(3,i),korelace1(4,i),krok,krok
endfor
endif
rectangle1,x,y,krok,krok

;-------------nakresleni histogramu-----------
widget_control,histdraw2,get_value=okno2 & wset,okno2
plot,hodnotyx,hist,xrange=[-1,1],xstyle=1,$
title='DISTRIBUTION OF CORRELATION COEFFICIENTS',$
xtitle='correlation coefficient',ytitle='density'
oplot,hodnotyx,gauss,color=255
arrow,100,303,140,303,hsize=0,color=255;,/device
;xyouts,100,303,'',/device
xyouts,150,300,'gauss fit',/device
xyouts,100,280,'mi='+strcompress(string(mi),/remove_all),/device
xyouts,100,265,'sigma='+strcompress(string(sigma),/remove_all),/device
end

'toggle1':prepinac=1
'toggle2':prepinac=2

'CORECT':begin
;---------smazani predchozi sekce, nahresleni site a textovych poli-----------
widget_control,innerbase,/destroy
innerbase=widget_base(base,/row)
obrazek=readfits(soubory(pocet/2))
s=size(obrazek)
crossdraw=widget_draw(innerbase,xsize=s(1),ysize=s(2))
tvscl,obrazek

for i=0,xpocet-1 do begin
for j=0,ypocet-1 do rectangle,xpocatek+i*krok,ypocatek+j*krok,krok,krok
endfor
for i=0,xpocet-1 do xyouts,xpocatek+i*krok-5,ypocatek-10,string(i+1),/device,charsize=0.8,alignment=0.5
for i=0,ypocet-1 do xyouts,xpocatek-20,ypocatek+i*krok+5,string(i+1),/device,charsize=0.8,alignment=0.5

corectbase=widget_base(innerbase,/column,xsize=250)
corectfield1=cw_field(corectbase,title='coordinates of 1st region (x,y)')
corectfield2=cw_field(corectbase,title='coordinates of 2nd region (x,y)')
corectfield3=cw_field(corectbase,title='start image (n)')
togbase=widget_base(corectbase,/column,/frame,/exclusive)
tog1=widget_button(togbase,value='center of brightness',uvalue='toggle1',/no_release)
tog2=widget_button(togbase,value='the brightest points',uvalue='toggle2',/no_release)
corectbut=widget_button(corectbase,value='CORRELATE',uvalue='CROSSCORECT')
widget_control,tog1,/set_button
prepinac=1
end

'CROSSCORECT':begin ;zde se provadi korekce podle teziste nebo nejjasnejsiho bodu
;---------parametry oblasti 1 a 2 - x1,y1,x2,y2,krok------------
widget_control,corectfield3,get_value=poradi & poradi=fix(poradi(0))
widget_control,corectfield1,get_value=sourad & sourad=sourad(0)
pozice=strpos(sourad,',') & delka=strlen(sourad)
xsourad=fix(strmid(sourad,0,pozice)) & ysourad=fix(strmid(sourad,pozice+1,delka-pozice))
x1=(xsourad-1)*krok+xpocatek & y1=(ysourad-1)*krok+ypocatek
widget_control,corectfield2,get_value=sourad & sourad=sourad(0)
pozice=strpos(sourad,',') & delka=strlen(sourad)
xsourad=fix(strmid(sourad,0,pozice)) & ysourad=fix(strmid(sourad,pozice+1,delka-pozice))
x2=(xsourad-1)*krok+xpocatek & y2=(ysourad-1)*krok+ypocatek
w=krok & h=krok
;--------vytvoreni poli pro ukladani integralu a teziste---------
rada1=fltarr(pocet);pro ukladani integralu
rada2=fltarr(pocet)
teziste1=intarr(2,pocet);pro ukladani teziste
teziste2=intarr(2,pocet)

if prepinac eq 1 then begin ;pokud vybrano teziste
iksy=findgen(w+1,h+1) ;pomocne pole pro vypocet teziste
iny=findgen(w+1,h+1)
for i=0,h do iksy(*,i)=findgen(w+1)
for i=0,w do iny(i,*)=findgen(h+1)
endif

obrazek=readfits(soubory(pocet/2))
s=size(obrazek)
animbase=widget_base(/column)
animdraw=widget_draw(animbase,xsize=s(1),ysize=s(2))
widget_control,animbase,/realize
widget_control,animdraw,get_value=okno & wset,okno

;-------------cyklus - pocitani polohy teziste---------------
for i=0,pocet-1 do begin ;4697 ;0,pocet-1

obr=readfits(soubory(i)) ;nacteme obrazek
if i ge poradi then begin ;prepocet souradnic oblasti az od urciteho cisla



;-----poloha teziste nebo nejjasnejsiho bodu pro prvni oblast------
intenzity1=obr(x1:x1+w,y1:y1+h) ;vytrhneme prvni oblast
if prepinac eq 1 then begin ;teziste
xt1=total(intenzity1*iksy)/total(intenzity1)
yt1=total(intenzity1*iny)/total(intenzity1)
xt1=round(xt1) & yt1=round(yt1) ;souradnice teziste v oblasti
endif else begin ;nejjasnejsi body
indexy1=where(intenzity1 eq max(intenzity1))
iksy1=indexy1 mod (w+1)
iny1=indexy1/(w+1)
xt1=total(iksy1)/n_elements(iksy1)
yt1=total(iny1)/n_elements(iny1)
endelse
teziste1(0,i)=x1+xt1 & teziste1(1,i)=y1+yt1 ;ulozeni polohy teziste
x1=x1+xt1-w/2 & y1=y1+yt1-h/2 ;souradnice nove oblasti

;-----poloha teziste nebo nejjasnejsiho bodu pro druhou oblast------
intenzity2=obr(x2:x2+w,y2:y2+h) ;vytrhneme druhou oblast
if prepinac eq 1 then begin ;teziste
xt2=total(intenzity2*iksy)/total(intenzity2)
yt2=total(intenzity2*iny)/total(intenzity2)
xt2=round(xt2) & yt2=round(yt2) ;souradnice teziste v oblasti
endif else begin ;nejjasnejsi bod
indexy2=where(intenzity2 eq max(intenzity2))
iksy2=indexy2 mod (w+1)
iny2=indexy2/(w+1)
xt2=total(iksy2)/n_elements(iksy2)
yt2=total(iny2)/n_elements(iny2)
endelse
teziste2(0,i)=x2+xt2 & teziste2(1,i)=y2+yt2 ;ulozeni polohy teziste
x2=x2+xt2-w/2 & y2=y2+yt2-h/2 ;souradnice nove oblasti

endif else begin
teziste1(0,i)=x1+w/2 & teziste1(1,i)=y1+h/2 ;pevna poloha teziste a oblasti
teziste2(0,i)=x2+w/2 & teziste2(1,i)=y2+h/2
endelse

;---------vypocet integralu v nove oblasti------------
rada1(i)=integruj(obr,x1,y1,w,h)-prumer(i)
rada2(i)=integruj(obr,x2,y2,w,h)-prumer(i) ;ulozeni integralu

;---------animace-------------
tvscl,obr
plots,x1+w/2,y1+h/2,psym=1,color=255,/device
plots,x2+w/2,y2+h/2,psym=1,color=255,/device
rectangle,x1,y1,w,h
rectangle,x2,y2,w,h
print,soubory(i) ;vypis
endfor


;------------vykresleni grafu----------------
;----------poloha teziste------------
stary=!p.multi & !p.multi=[0,2,2,0,1]
plot,teziste1(0,*),teziste1(1,*),psym=2,title='1st region', $
xtitle='x position',ytitle='y position',$
xrange=[min(teziste1(0,*))-5,max(teziste1(0,*))+5],$
yrange=[min(teziste1(1,*))-5,max(teziste1(1,*))+5],xstyle=1,ystyle=1
plot,teziste2(0,*),teziste2(1,*),psym=2,title='2nd region', $
xtitle='x position',ytitle='y position', $
xrange=[min(teziste2(0,*))-5,max(teziste2(0,*))+5],$
yrange=[min(teziste2(1,*))-5,max(teziste2(1,*))+5],xstyle=1,ystyle=1
!p.multi=[4,2,4,0,1]
plot,teziste1(0,*),title='1st region - x position',xrange=[0,pocet],charsize=1.5,$
yrange=[min(teziste1(0,*))-5,max(teziste1(0,*))+5],xstyle=1,ystyle=1
plot,teziste1(1,*),title='1st region - y position',xrange=[0,pocet],charsize=1.5,$
yrange=[min(teziste1(1,*))-5,max(teziste1(1,*))+5],xstyle=1,ystyle=1
plot,teziste2(0,*),title='2nd region - x position',xrange=[0,pocet],charsize=1.5,$
yrange=[min(teziste2(0,*))-5,max(teziste2(0,*))+5],xstyle=1,ystyle=1
plot,teziste2(1,*),title='2nd region - y position',xrange=[0,pocet],charsize=1.5,$
yrange=[min(teziste2(1,*))-5,max(teziste2(1,*))+5],xstyle=1,ystyle=1
!p.multi=stary
write_jpeg,'correct2.jpg',tvrd()

;-------------korelace poloh teziste--------------------
tezbase=widget_base(/column)
tezdraw=widget_draw(tezbase,xsize=400,ysize=600)
widget_control,tezbase,/realize
widget_control,tezdraw,get_value=okno & wset,okno
lag=findgen(2*pocet-1)-pocet+1
resx=c_correlate(teziste1(0,*),teziste2(0,*),lag)
resy=c_correlate(teziste1(1,*),teziste2(1,*),lag)
koefx=resx(pocet-1) & koefx=strcompress(string(koefx),/remove_all)
koefy=resy(pocet-1) & koefy=strcompress(string(koefy),/remove_all)
stary=!p.multi & !p.multi=[0,1,2]
plot,lag,resx,title='CORRELATION OF X COORDINATE c0='+koefx,xtitle='lag',ytitle='correlation coefficient',$
yrange=[-1,1],ystyle=1
plot,lag,resy,title='CORRELATION OF y COORDINATE c0='+koefy,xtitle='lag',ytitle='correlation coefficient',$
yrange=[-1,1],ystyle=1
!p.multi=stary
write_jpeg,'correct3.jpg',tvrd()

;----------korelace oblasti po korekci
korelbase=widget_base(/column)
koreldraw=widget_draw(korelbase,xsize=400,ysize=300)
widget_control,korelbase,/realize
widget_control,koreldraw,get_value=okno & wset,okno
lag=findgen(2*pocet-1)-pocet+1
res=c_correlate(rada1,rada2,lag)
koef=res(pocet-1) & koef=strcompress(string(koef),/remove_all)
plot,lag,res,title='CROSS CORRELATION c0='+koef,xtitle='lag',ytitle='correlation coefficient',$
yrange=[-1,1],ystyle=1
write_jpeg,'correct1.jpg',tvrd()

end
endcase
end


pro korel
COMMON start,omenu
COMMON baze,base,innerbase
base=widget_base(xsize=1000,ysize=600,mbar=bar)
fmenu=widget_button(bar,value='FILE',/menu)
f1=widget_button(fmenu,value='LOAD WORKING DIRECTORY',uvalue='LOAD')
f2=widget_button(fmenu,value='QUIT',uvalue='QUIT')
omenu=widget_button(bar,value='OPTIONS',/menu)
o1=widget_button(omenu,value='whole image',/menu)
o11=widget_button(o1,value='count integrals',uvalue='INTEGRALS')
o12=widget_button(o1,value='count correlation coefficients',uvalue='KOEFICIENTY')
o13=widget_button(o1,value='show regions with high correlation',uvalue='HISTOGRAM')
o14=widget_button(o1,value='show regions which correlate with too many regions',uvalue='MNOHO')
o2=widget_button(omenu,value='chosen regions',/menu)
o21=widget_button(o2,value='auto correlation',uvalue='AUTO')
o22=widget_button(o2,value='cross correlation',uvalue='CROSS')
o23=widget_button(o2,value='light curve',uvalue='LIGHT')
o24=widget_button(o2,value='distribution',uvalue='DIST')
o3=widget_button(omenu,value='correction',/menu)
o31=widget_button(o3,value='correlation after correction',uvalue='CORECT')
innerbase=widget_base(base)

widget_control,base,/realize
widget_control,omenu,sensitive=0

d=widget_draw(innerbase,xsize=1000,ysize=600)
for i=0,999 do arrow,i,0,i,600,hsize=0,color=25100+(i/4.1)
tl=3 & vel=15 & ypos=250
xyouts,200,ypos,'K',alignment=0.5,/device,charsize=vel,charthick=tl,color=25340
xyouts,350,ypos,'O',alignment=0.5,/device,charsize=vel,charthick=tl,color=25280
xyouts,500,ypos,'R',alignment=0.5,/device,charsize=vel,charthick=tl,color=25250
xyouts,650,ypos,'E',alignment=0.5,/device,charsize=vel,charthick=tl,color=25180
xyouts,800,ypos,'L',alignment=0.5,/device,charsize=vel,charthick=tl,color=25120

xmanager,'korel',base
end