function Arot=picrotate(A,ang,modus);
%% picrotate - Rotate a Matrix (Picture) by the angle ang [rad]
%%
%% Arot=picrotate(A,ang,[modus])
%%
%% At first the position of non-zero elements are determined.
%% On the positions we perform the rotation. The rotated picture
%% will have other dimensions (depending on the angle), so that
%% a further normalization to new pixel positions will be performed.
%%
%% A new assignment to the pixel with help of sparse functions is done.
%% The filling of possible gaps is done by mean value calculation of
%% surrounding non-zero elements.
%%
%% The optional input argument modus determines, who multiple assigned
%% positions and gaps after rotation are handled:
%%
%% modus=0 (Default)	: The mean value ist taken
%% modus=1			: The maximum value is used
%%
%% A.J. Geissler, February 2010, Seeheim-Jugenheim, Germany

if nargin<3,
	modus=0;
endif

[Z,S]=size(A);									% Determine the size of the input picture
[pz,ps,V]=find(A);								% Determine positions of non-zero elements
PosOld=[pz.';ps.'];								% Represent positions as matrix
R=[cos(ang),-sin(ang);sin(ang),cos(ang)];		% The Transformation Matrix
PosNew=(R*(PosOld-1)).';						% New Positions by multiplication with TR-Matrix
pz2=PosNew(:,1);								% Extract new lines and columns from the result
ps2=PosNew(:,2);

Zneu=(Z-1)*abs(cos(ang)) + (S-1)*abs(sin(ang));		% Calculate the new dimensions of the rotated picture
Sneu=(Z-1)*abs(sin(ang)) + (S-1)*abs(cos(ang));

pz3=floor(Zneu* (pz2 - min(pz2)) ./(max(pz2) - min(pz2))) + 1;	% Normalize the new coordinates to new
ps3=floor(Sneu* (ps2 - min(ps2)) ./(max(ps2) - min(ps2))) + 1;	% dimensions. Take Care: GAPS !
Zneu=Zneu+1;													% Here we adopt the range to starting index 1
Sneu=Sneu+1;													% of lines and columns

if modus!=0,
	Arot=full(sparse(pz3,ps3,V,Zneu,Sneu,"unique"));	% Assign the pixel values to new coordinates
														% with the option "unique", the max value of pixels
														% which are assigned to one position is taken
else
	Nput=full(sparse(pz3,ps3,ones(size(pz3)),Zneu,Sneu));	% Number of Assignments to new positions
	Arot=full(sparse(pz3,ps3,V,Zneu,Sneu));					% Assign pixel values to new coordinates
	Arot(Nput!=0)=Arot(Nput!=0) ./Nput(Nput!=0);			% If several pixels are assigned to one position,
															% we take the mean value of them.
endif

Z=zeros(Zneu-2,Sneu-2);						% Now we determine the number of empty positions
											% by analysis of its 8 neighbours
for p=-1:1:1,
	for s=-1:1:1,
		if !(p==0 & s==0),						% Don't take the origin pixel !
			Aakt=Arot(2+p:1:Zneu-1+p,2+s:1:Sneu-1+s);
			Z=Z + (Aakt==0);				% If a neighbour pixel is zero, we increment
		endif								% the counter
	endfor
endfor
		
[pz,ps]=find(Z);							% These are the positions of gaps
[Zneu,Sneu]=size(Arot);

ArotVek=Arot.';								% The rotated picture is vectorized in order
ArotVek=ArotVek(:);							% to perform further sub-indices-adressings 

pzakt=pz;									% Calculate the subindices for empty pixels
psakt=ps;
pzm0=pzakt((pzakt>1) & (pzakt<(Zneu-1)) & (psakt>1) & (psakt<(Sneu-1)));
psm0=psakt((pzakt>1) & (pzakt<(Zneu-1)) & (psakt>1) & (psakt<(Sneu-1)));
indices0=psm0 + (pzm0 -1) .*Sneu;

Z=[zeros(1,Sneu-2);Z;zeros(1,Sneu-2)];			% Here we fit the dimension of the Gap-Counter
Z=[zeros(Zneu,1),Z,zeros(Zneu,1)];				% matrix to the size of the rotated picture
Zvek=Z.'(:);									% Vectorize -> Subadressing

Sfill=double(ArotVek);
Sfill(indices0)=zeros(length(indices0),1);

for p=-1:1:1,									% Empty pixels will be filled with the
	for s=-1:1:1,								% mean value of non-empty pixels.
		if !(p==0 & s==0),
			pzm=pzm0+p;
			psm=psm0+s;
			indices=psm + (pzm -1) .*Sneu;
			if modus==0,
				Sfill(indices0)=Sfill(indices0) + (ArotVek(indices) ./(8-Zvek(indices0)));
			else
				Sfill(indices0)=max(Sfill(indices0),ArotVek(indices));
			endif
		endif
	endfor
endfor

Arot=reshape(Sfill,Sneu,Zneu).';					% Convert back to picture

endfunction;