clear();
# deklaracja zmiennych
raster R, G, B, Ir,Pan;
raster Rc, Gc, Bc, Irc ;
raster Rres, Gres, Bres, Ires;
raster Rpan, Gpan, Bpan, Pancor,W;
# otwarcie rastrów wejściowych
GetInputRasters(R,G,B,Ir,Pan);
# odczytanie ilości wierszy, kolumn i typu rastru panchromatycznego
numeric lin=NumLins(Pan);
numeric col=NumCols(Pan);
string type$=RastType(Pan);
# utworzenie rasstrów tymczasowych - spektralnych  dopasowanych do rastra Panchromatycznego
CreateTempRaster(Rres,lin,col,type$);
CreateTempRaster(Gres,lin,col,type$);
CreateTempRaster(Bres,lin,col,type$);
CreateTempRaster(Ires,lin,col,type$);
ResampleRasterToMatch(R,Pan,Rres, "affine","cubic");
ResampleRasterToMatch(G,Pan,Gres, "affine","cubic");
ResampleRasterToMatch(B,Pan,Bres, "affine", "cubic");
ResampleRasterToMatch(Ir,Pan,Ires, "affine", "cubic");
# filtr usredniajacy majacy na celu eliminacje wartosci skrajnych, poprawe kształtu histogramu
CreateTempRaster(Rresf,lin,col,type$);
CreateTempRaster(Gresf,lin,col,type$);
CreateTempRaster(Bresf,lin,col,type$);
CreateTempRaster(Iresf,lin,col,type$);
CreateTempRaster(Panf,lin,col,type$);
Rresf=FocalMedian(Rres,3,3);
Gresf=FocalMedian(Gres,3,3);
Bresf=FocalMedian(Bres,3,3);
Iresf=FocalMedian(Ires,3,3);
Panf=FocalMedian(Pan,3,3);
## Obliczenia związane z Pansharp ESRI
#numeric r, g, b, i; wskazniki wagowe, korekcja Pan przez Ir 
#r=0.166; b= 0.167; g=0.167; i=0.500;
GetOutputRaster(Rpan,lin,col,type$);
GetOutputRaster(Gpan,lin,col,type$);
GetOutputRaster(Bpan,lin,col,type$);
GetOutputRaster(Pancor,lin,col,"32-bit float");
# korekcja kanału panchromatycznego
Pancor=(Panf)/(Rresf+Gresf+Bresf);
# korekcja kanałów spektralnych
#!!!!!!!!!!!!!problem obliczeniowy z RPAN, dostepem obliczeyc piramide; jest skuteczna
Rpan=Rresf*Pancor; CreatePyramid(Rpan);
Gpan=Gresf*Pancor; CreatePyramid(Gpan);
Bpan=Bresf*Pancor; CreatePyramid(Bpan);
#przeskalowanie do skali 8-bitowej;
raster Rs, Gs, Bs;
CreateTempRaster(Rs,lin,col,"8-bit unsigned", 0);
CreateTempRaster(Gs,lin,col,"8-bit unsigned", 0);
CreateTempRaster(Bs,lin,col,"8-bit unsigned", 0);
numeric rmax=GlobalMax(Rpan); 
	numeric rmin=GlobalMin(Rpan);
		numeric rrange=rmax-rmin;print(rmin,rmax,rrange);
			Rs=(((Rpan-rmin)/rrange)*254)^1.04;
numeric gmax=GlobalMax(Gpan);
	numeric gmin=GlobalMin(Gpan);
		numeric grange=gmax-gmin;
			Gs=(((Gpan-gmin)/grange)*255)^1.07;
numeric bmax=GlobalMax(Bpan);
	numeric bmin=GlobalMin(Bpan);
		numeric brange=bmax-bmin;
			Bs=(((Bpan-bmin)/brange)*255)^0.96;
# Złozenie i wyświetlenie rastra złożonego;
GetOutputRaster(PanRGB, lin,col, "24-bit color RGB");
file$=GetObjectFileName(PanRGB);
onumber=GetObjectNumber(PanRGB);
object$=GetObjectName(file$,onumber);
print(file$);
RasterRGBToComposite(Rs,Gs,Bs,PanRGB,file$,object$, "Composite 24-bit raster");
class DISP disp1;
disp1=DispOpen(450,450,200,200);
#DispSetTitle(disp1, "Display 1"); # set the titles
DispQuickAddRasterVar(disp1,PanRGB);
DispRedraw(disp1);