#!!!!!!!!!!!!!!!!!!!!!!!!!
#POLECENIE:
#=wykonać  opis działania skryptu w miejscach, gdzie zaznaczono #...
#prosze zrocić uwagę na obliczenia wykonywane na macierzach, bezpośrednio w RAM-ie, 
#co przyspiesza działanie skryptu wielokrotnie (100-1000 razy);
clear();
raster Bin, Oper, Rcor;#...
GetInputRasters(Bin, Oper);
numeric lbin, cbin;
lbin=NumLins(Bin);#...
cbin=NumCols(Bin);#...
numeric loper, coper;#...
loper=NumLins(Oper);#...
coper=NumCols(Oper);#...
string typebin$, typeoper$;#...
typebin$=RastType(Bin);#...
typeoper$=RastType(Oper);
if (typebin$==typeoper$) then goto label1;
else goto labelend;
label1:
numeric lmax, cmax;
lmax=lbin-loper;print(lmax);
cmax=cbin-coper;print(cmax);
GetOutputRaster(Rcor, lbin, cbin, "32-bit float");
Rcor=0;
numeric korwar;
korwar=PopupNum("Wpisz wartość graniczna dla współczynnika korelacji",0.8,0,1,2);
#...
string fb$=GetObjectFileName(Bin);#...
numeric ob=GetObjectNumber(Bin);#...
string nameb$=GetObjectName(fb$,ob);#...
string pathb$=FileNameGetPath(fb$);##...
#nazwy obiektów i plików Oper
string fo$=GetObjectFileName(Oper);#...
numeric oo=GetObjectNumber(Oper);#...
string nameo$=GetObjectName(fo$,oo);#...
string patho$=FileNameGetPath(fo$);#...
#nazwy obiektów i plików Oper
string fc$=GetObjectFileName(Rcor);
numeric oc=GetObjectNumber(Rcor);
string namec$=GetObjectName(fc$,oc);
string pathc$=FileNameGetPath(fc$);

#Obliczenia dla wskazanego rastra elementu strukturalnego
#...
raster Bin1, Oper1, Rcor1;
OpenRaster(Bin1, fb$, nameb$);
OpenRaster(Oper1, fo$, nameo$);
OpenRaster(Rcor1, fc$, namec$);
#...
array numeric  oper[loper,coper]; 
numeric l1, c1;
numeric sumoper1=0;
numeric sumoper2=0;
numeric pix=loper*coper
for l1=1 to loper
{
	 for c1=1 to coper
	{
		oper[l1,c1]=Oper[l1,c1];#...
		sumoper1=oper[l1,c1];#...
		sumoper2 += sumoper1;#...
	}
}
numeric sredOperm=sumoper2/pix;print(sredOperm);
#macierz różnic wartosci pikseli od sredniej dla rastra Oper 
array numeric operdef[loper,coper];
numeric sumdef1=0;
numeric sumdef2=0;
numeric l3, c3;
for l3=1 to loper
	{
		 for c3=1 to coper
		{
			operdef[l3,c3]=oper[l3,c3]-sredOperm;
			sumdef1=operdef[l3,c3]*operdef[l3,c3];
			sumdef2 += sumdef1;
		}
	}

numeric stddevoper=sqrt(sumdef2/(pix-1)); print(stddevoper);
#-----------------------------------
# główna pętla obliczeniowa
#----------------------------------
numeric l,c;
for l=1 to lmax
{
for c=1 to cmax
{ 
#wycinek rastra Bin porównywany z rastrem Oper, obliczenia sredniej wartości pikseli analizowanego wycinka
array numeric fragment[loper, coper];
numeric l5, c5;numeric sumfr1=0;numeric sumfr2=0;
for l5=1 to loper
{
	 for c5=1 to coper
	{
			numeric w,r;	w=l+l5-1; r=c+c5-1;
			fragment[l5,c5]=Bin1[w,r];
			sumfr1=fragment[l5,c5];			sumfr2 += sumfr1;
	}
}
numeric sredfrm=sumfr2/pix;#print(sredfrm);
#macierz roznic od sredniej fragmentu
numeric l7, c7;
array numeric fragdef[loper,coper];
numeric sumfrdef1=0;numeric sumfrdef2=0;
for l7=1 to loper
	{
		 for c7=1 to coper
		{
			fragdef[l7,c7]=fragment[l7,c7]-sredfrm;
			sumfrdef1=fragdef[l7,c7]*fragdef[l7,c7];			sumfrdef2 += sumfrdef1;
		}
	}
numeric stddevfr=sqrt(sumfrdef2/(pix-1)); #print(stddevfr);
#...
numeric l9, c9;
array numeric iloczyn[loper,coper];
numeric sumil1=0;numeric sumil2=0;
for l9=1 to loper
	{
		 for c9=1 to coper
		{
			iloczyn[l9,c9]=fragdef[l9,c9]*operdef[l9,c9];
			sumil1=iloczyn[l9,c9];			sumil2 += sumil1;
		}
	}
#suma iloczynow roznic
if (stddevfr==0) then 
{numeric kor=0; goto label2;}
else numeric stddevfr1=stddevfr;#...
kor=(sumil2/(pix-1))/(stddevoper*stddevfr1);#...
label2:
Rcor1[l+loper/2,c+coper/2]=kor;
string l$=NumToStr(l); string c$=NumToStr(c);
string m$=l$+"-"+c$;
SetStatusMessage(m$); 
}
}
print("Koniec obliczania korelacji obrazowej");
CloseRaster(Bin1); CloseRaster(Oper1); CloseRaster(Rcor);
#Koniec głównej Pętli obliczeniowej
#----------------------------------------
#czyszczenie rastra bin
#Otwarcie rastrów Bin i Oper i Rcor
raster Bin2, Oper2, Rcor2, OperC, Filt, Filt1;
OpenRaster(Bin2, fb$, nameb$);
OpenRaster(Oper2, fo$, nameo$);
OpenRaster(Rcor2, fc$, namec$);
#wykonanie kopii rastrów Bin2 do czyszczenia
string namefilt$=pathb$+"/filtrowany.rvc";
CreateProjectFile(namefilt$,"Filtrowany");
CreateRaster(Filt, namefilt$, nameb$,"Filtrowany", lbin,cbin,"binary");
OpenRaster(Filt1,namefilt$,nameb$);
Filt1=Bin2;
CloseRaster(Filt1);
CreateTempRaster(OperC,loper, coper, "binary");
OperC=1;
raster Filt2;
OpenRaster(Filt2,namefilt$,nameb$);
IgnoreNull(Filt2);
numeric l11, c11;
for l11=1 to lbin-loper/2#...
	{
		 for c11=1 to cbin-coper/2#...
		{
			if (Rcor2[l11, c11]>korwar) then #...
							{
									print(l11,c11);#...
									for l12=1 to loper#...
											{
		 										for c12=1 to coper
												{
												Filt2[l11-loper/2+l12,c11-coper/2+c12]=1#...
												}
											}
							}
				string l$=NumToStr(l11); string c$=NumToStr(c11);
				string m$=l$+"-"+c$;
				SetStatusMessage(m$); #...
		}
	}
CloseRaster(Bin2);CloseRaster(Oper2); CloseRaster(Rcor2); 
CloseRaster(Filt); CloseRaster(Filt1); CloseRaster(Filt2); CloseRaster(OperC);
goto koniec;
labelend:
print("Raster Bin i oper mają różne kodowanie, koniec skryptu");
 
koniec: