clear();
raster A, B, C, NA, RL,GL;
GetInputRaster(C);
string s$=GetObjectFileName(C);
array objnums[20];
CloseRaster(C);     
numeric i;
numeric max=GetAllObjectNumbers(s$,objnums,"Raster");
data=Date();
print("Przetwarzanie wykonane  - ", data);
print("-----------------------");
numeric f, g, m, n;
numeric dmin, dzakr;
numeric srodlin;numeric srodcol; numeric lins; numeric cols;
lins = NumLins(C);
cols = NumCols(C);
# domyslna rozdzielczość
numeric resdef=lins/(24/2.54);
#domyslny numer wiersza punktu głównego
numeric srodlindef=lins/2;
#domyślny numer kolumny punktu głównego
numeric srodcoldef=cols/2;
# domylsna wartość numeru srodka punktu głownego
f = PopupNum("Wprowadź ogniskową kamery lotniczej w [mm] f: ", 152.5, 50, 300, 3);
g= PopupNum("Wprowadź współczynnik gamma: ", 0.60, 0.25, 4, 2);
n = PopupNum("Wprowadź wykładnik potęgi dla funkcji cosinus n: ", 4.0,1.0, 8.0, 1);
m = PopupNum("Wprowadź rozdzielczość skanowania w punktach na cal [dpi]: ", resdef, 100, 5000, 0);
srodlin = PopupNum("Wprowadź numer linii dla punktu głównego zdjęcia srodlin: ", srodlindef,1,25000,0);
srodcol = PopupNum("Wprowadź numer kolumny dla punktu głównego zdjęcia srodcol: ", srodcoldef, 1, 25000,0);
dzakr = PopupNum("Wprowadź średni zakres gęstości optycznej:", 2.5, 1.0, 4.0, 1);
string path$ = GetDirectory( defaultpath$, "Wskaż katalog obliczeniowy" );
string outputfile$ = GetOutputFileName(path$, "Wprowadź nazwę pliku rvc dla obrazów skorygowanych", ".rvc" );
string desc1$="korekcja winietowania";
print(outputfile$);CreateProjectFile(outputfile$,desc$);
#--------------------------
#szacowanie wykladnika potegi
# gromadzenie danych o jasności na kolejnych liniach profilowych przechodzących przez srodoek zdjęcia w trzech kanalach
# 	- wyznacznie poczatku i konca linii profilewej w oparciu o krok katowy
#		- zapis wyników pomiaru jasności do macierzy (odleglośc, jansośc)
#		- oszacowanie parametrów wielomianu drugiego rzędu dl alinii profilowej
#   - obliczenie janosci w zaleznosci od odleglosci od srodka zdjęcia zyskana funkja oddielnie dla dwóch polprostych
#   - oddzielne oszacowanie wykladnia n dla polprostych równaniem z wykorzystanie logarytmow dziesietnych
#   - usrednienie parametru n dla pojednczego kanalu, nastepnie dla wszystkich trzech kanalow (?)!!!!!!!!


# Tworzenie maski korekcyjnej
SetStatusMessage("Tworzenie maski korekcyjnej");
raster K;
string f2$="C:/maska_kor.rvc";
string k$="mask_kor";
CreateRaster(K, f2$, k$, "maska korekcyjna", lins, cols, "16-bit unsigned");

for each K[lins,cols] begin
  	  	K[lins,cols]=10000*1/pow(cosd(atand(sqrt((pow((lins-srodlin),2)+pow((cols-srodcol),2)))*(25.4/m)/f)), n);
   end
CloseRaster(K);
#stworzenie macierzy przeliczeń jasności na wartości względnego naświetlenia
numeric ar;
ar = dzakr/255;
array dn[255]; hmax=255;
for i=0 to hmax
{
 dn[i]=i;print(dn[i]);
}
numeric h;
array h[255]
for i=0 to hmax
{
h[i]=pow(10,((dn[i]*ar)/g));print(h[i]);
}
for i=1 to max
{
numeric a=objnums[i];
string s3$=GetObjectName(s$,a);
print("Zdjęcie - ", s3$);
OpenRaster(A,s$,s3$);
string f1$="C:/temp-rgb.rvc";
string desc$="Składowe RGB";
CreateProjectFile(f1$,desc$);
numeric lins = NumLins(A);
numeric cols = NumCols(A); 
string datatype$="8-bit unsigned";
raster R,G,B;
RasterCompositeToRGB(A, R, G, B,f1$,desc$);
CloseRaster(A);
 
# Zamiana jasnosci pikseli na wartości gęstości optycznej standardowej
SetStatusMessage("Korekcja składowych RGB - obliczenia parametrów");
raster R1, G1, B1, K1;
OpenRaster(R1, f1$, "Red");
OpenRaster(G1,f1$, "Green");
OpenRaster(B1, f1$, "Blue");
OpenRaster(K1,f2$,k$);

raster R2, G2, B2;
string f3$= "C:/korekcja.rvc"
CreateRaster(R2, f3$, "R2", "czerwony_gestosc_kor", lins, cols, "8-bit unsigned");  
CreateRaster(G2, f3$, "G2", "zielony_gestosc_kor", lins, cols, "8-bit unsigned");  
CreateRaster(B2, f3$, "B2", "niebieski_gestosc_kor", lins, cols, "8-bit unsigned");
SetStatusMessage("Korekcja składowej czerwonej");
for each R2[lins,cols] begin
	br=R1[lins,cols];
	R2[lins,cols] = (1/ar)*(g*log10(h[br]*K1/10000));
	end
CloseRaster(R2);CloseRaster(R1);
SetStatusMessage("Korekcja składowej zielonej");
for each G2[lins,cols] begin
	bg=G1[lins,cols];
	G2[lins,cols] = (1/ar)*(g*log10(h[bg]*K1/10000));
	end
CloseRaster(G2);CloseRaster(G1);
SetStatusMessage("Korekcja składowej niebieskiej");
for each B2[lins,cols] begin
	bb=B1[lins,cols];
	B2[lins,cols] = (1/ar)*(g*log10(h[bb]*K1/10000));
	end
CreateHistogram(B2);
CreatePyramid(B2);
CloseRaster(B2);CloseRaster(B1);
CloseRaster(K1);
# wykonanie kompozycji do rastra złożonego w pliku o nazwie podobnej - szer1-korwiniet
SetStatusMessage("Konwersja do postacji rastra złożonego 24-bit RGB");
string s3dod$="-kor";
string s4$=s3$+s3dod$;
string desc5$="skorygowane winietowanie";
raster R3, G3,B3, COMP;
OpenRaster(R3, f3$, "R2");  
OpenRaster(G3, f3$, "G2");  
OpenRaster(B3, f3$, "B2");
RasterRGBToComposite(R3,G3,B3,COMP,outputfile$,s4$,desc5$);
CreateHistogram(COMP);
CreatePyramid(COMP);
CloseRaster(R3);CloseRaster(G3);CloseRaster(B3);
raster A1;
OpenRaster(A1,s$,s3$);
CopySubobjects(A1,COMP, "georef");
CloseRaster(A1);
# usuwanie obiektow tymczasowych
DeleteFile(f1$);
DeleteFile(f3$);
print("-----------------");
# kompresja rastra wynikowego
SetStatusMessage("Kompresja rastra złożonego 24-bit RGB");
string s4dod$="comp";
string s5$=s4$+s4dod$;
numeric t=strlen(outputfile$);
numeric t1=t-7;
string r$=left$(outputfile$,t1);
string s5dod$=".rvc";
string outcomp$=r$+s4dod$+s5dod$;
raster COMP1, COMP2;
OpenRaster(COMP1,outputfile$,s4$);
RasterCompress(COMP1,COMP2, outcomp$, s5$, desc5$, "JP2", 10);
CloseRaster(COMP1);
}
DeleteFile(f2$);