## Copyright (C) 2008 David Barina

## -*- texinfo -*-
## @deftypefn {Function File} {[@var{len}, @var{angle}] =} blur_c(@var{image})
## Rozmaze obrazek (odstiny sedi nebo 3 kanaly) pod zadanou delkou a uhlem.
##
## Pokusi se odhadnout delku a uhel rozmazani. Je urcena
## jen pro obrazky poskozene rovnomernym rozmazanim
## po primce.
## @end deftypefn

## Author: David Barina  <ibarina@fit.vutbr.cz>
## Created: 2008-11-16

function blur_c(fnamein, fnameout='blurred.png', length=10, angle=0)
	% nacist original
	IM = imread(fnamein);
	IM_IN=double(IM);

	% zjistit pocet kanalu
	S=size(IM_IN);
	if(size(S)(2) == 3)
		L=S(3);
	else
		L=1;
	endif

	printf("Pocet vrstev obrazu:\n");
	L

	% horizontalni konvolucni jadro
	h=zeros(length,length);
	h(1,:) = ones(1,length);

	% otocit o dany uhel
	hx = imrotate(h,angle,"bicubic","loose",0);
	% normalizace
	hx ./= sum(hx(:));

	% rozmazani obrazku
	IM_OUT = conv2(IM_IN(:,:,1), hx, 'valid');
	for i = 2:L
		IM_OUT(:,:,i) = conv2(IM_IN(:,:,i), hx, 'valid');
	endfor

	IM_OUT=uint8(IM_OUT);

	% ulozi vysledek
	imwrite(IM_OUT, fnameout);

endfunction

## Copyright (C) 2008 David Barina

## -*- texinfo -*-
## @deftypefn {Function File} {[@var{len}, @var{angle}] =} detect_c(@var{image})
## Detekuje delku a smer rozmazani z obrazku.
##
## Pokusi se odhadnout delku a uhel rozmazani. Je urcena
## jen pro obrazky poskozene rovnomernym rozmazanim
## po primce.
## @end deftypefn

## Author: David Barina  <ibarina@fit.vutbr.cz>
## Created: 2008-11-16

function [a b]=detect_c(im)
	% zjisteni poctu kanalu obrazu
	S=size(im);
	if(size(S)(2) == 3)
		L=S(3);
	else
		L=1;
	endif
	printf("Pocet vrstev obrazu:\n");
	L

	% RGB obraz pripadne prevest na intenzitu
	if(L == 3)
		IMG=double(im)(:,:,1);
		IMG(:,:) = 0.299.*double(im)(:,:,1) .+ 0.114.*double(im)(:,:,2) .+ 0.587.*double(im)(:,:,3);
	else
		IMG=double(im)(:,:,1);
	endif

	% najit minima/minimum v cepstru
	MIN=min(min( real(fftshift(ifft2(log10(abs(fft2(IMG)))))) ));

	% na pozici minim dat 1, jinak vsude 0
	IMG2=uint8(real(fftshift(ifft2(log10(abs(fft2(IMG))))))<=MIN);

	% vyska cepstra
	Y=size(IMG2)(1);

	%minumum by melo byt ve spodni polovine
	IMG3=IMG2(floor(Y/2)+1:Y,:);

	% minimum je v horni polovine
	if(max(max(IMG3))==min(min(IMG3)))
		IMG3 = IMG2(1:floor(Y/2)+1,:);
		IMG3 = imrotate(IMG3,180);
	endif

	% neni minumum?
	if(max(max(IMG3))==min(min(IMG3)))
		printf("Chyba: Neni minimum ve spektru?!?! Co to znamena???\n");
	endif

	% zjistit pozici minima
	[mx,ix]=max(max(IMG3));
	[mx,iy]=max(max(IMG3'));
	[sy,sx]=size(IMG3);
	ix = ix-1;
	iy = iy-1;

	% poloha minima vuci stredu
	dx = ceil(ix-(sx/2));
	dy = iy;

	% spocitat vzdalenost od stredu a uhel
	delka = abs(dx+dy*i);
	uhel = (angle(dx+dy*i)/pi)*180-180;
	if(uhel==-180)
		uhel=0;
	endif

	% vratit hodnoty
	a = delka;
	b = uhel;
endfunction

## Copyright (C) 2008 David Barina

## -*- texinfo -*-
## @deftypefn {Function File} deblur_c(@var{infilename}, @var{outfilename}, @var{K}, @var{length}, @var{angle})
## Odstrani rozmazani obrazku (stupne sedi nebo barevny).
##
## Pokusi se odhadnout delku a uhel rozmazani. Je urcena
## jen pro obrazky poskozene rovnomernym rozmazanim
## po primce.
## @end deftypefn

## Author: David Barina  <ibarina@fit.vutbr.cz>
## Created: 2008-11-16

function deblur_c(fnamein, fnameout='deblur.png', K=1e-2, length=-1, angle=-1)
	% nacist obrazek
	IM = imread(fnamein);

	% zjisteny poctu kanalu (vrstev)
	S=size(IM);
	if(size(S)(2) == 3)
		L=S(3);
	else
		L=1;
	endif
	printf("Pocet vrstev obrazu:\n");
	L

	% originalni obrazek v doublech
	IM_BLUR=double( IM );

	% zobrazit (prvni kanal pouze)
	figure(1); imshow(uint8(IM_BLUR(:,:,1)), [0 255]); title('blurred');

	% zjisteni rozmeru obrazku
	if(L == 1)
		[i1 j1] = size(IM_BLUR);
	else
		[i1 j1 k1] = size(IM_BLUR);
	endif

	% pridat okraje kvuli nasobeni ve frekvencni oblasti
	IM_BLUR_9 = [repmat(IM_BLUR(1,1,:),i1,j1) repmat(IM_BLUR(1,:,:),i1,1) repmat(IM_BLUR(1,end,:),i1,j1);
		repmat(IM_BLUR(:,1,:),1,j1) IM_BLUR repmat(IM_BLUR(:,end,:),1,j1);
		repmat(IM_BLUR(end,1,:),i1,j1) repmat(IM_BLUR(end,:,:),i1,1) repmat(IM_BLUR(end,end,:),i1,j1)];

	% detekovat poskozeni
	printf("Detekovano: \n");
	[delka uhel] = detect_c(IM_BLUR)

	% pouzit pripadne vynucene parametry
	if(length == -1)
		length = delka;
	endif

	length = round(length);

	if(angle == -1)
		angle = uhel;
	endif

	printf("Pouziji tyto hodnoty:\n");
	K
	length
	angle

	% vytvorit konvolucni jadro simulujici poskozeni
	hr=zeros(length,length);
	hr(1,:) = ones(1,length);
	hrx = imrotate(hr,-angle,"bicubic","loose",0);
	hrx ./= sum(hrx(:));

	% velikost obrazku i s okrajema
	if(L == 1)
		[i j]=size(IM_BLUR_9);
	else
		[i j k]=size(IM_BLUR_9);
	endif

	% vytvorit filtr simulujici poskozeni
	h3=zeros(i,j);
	h3(1:size(hrx)(1),1:size(hrx)(2)) = hrx;
	H3=fft2(h3);

	% vytvorit inverzni filtr ("parametric version of the Wiener filter")
	IF=conj(H3)./(abs(H3).^2 + K);

	% provest inverzni filtraci (dekonvoluci)
	if(L == 1)
		IMd=abs(ifft2(fft2(IM_BLUR_9).*(IF)));
	else
		IMd=abs(ifft2(fft2(IM_BLUR_9(:,:,1)).*(IF)));
		IMd(:,:,2)=abs(ifft2(fft2(IM_BLUR_9(:,:,2)).*(IF)));
		IMd(:,:,3)=abs(ifft2(fft2(IM_BLUR_9(:,:,3)).*(IF)));
	endif

	% vyriznout jen puvodni obrazek bez okraju
	IMd1 = IMd(1+i1:i1+i1,1+j1:j1+j1,:);

	% zobrazit vysledek (jen prvni kanal)
	figure(6); imshow( uint8( IMd1(:,:,1) ), [0 255]); title('deblurred');

	% ulozi vysledek
	imwrite(uint8(IMd1), fnameout);

endfunction