function [totalcoef,out] = jpg2000(file,N)

% im = imread(file);
im = file;
width = size(im,2);
height = size(im,1);
siz1 = max(width,height);
b = round(log2(siz1));
siz2 = 2^b;
im = double(imresize(im,[siz2 siz2],'bicubic'));
original=im;
width = size(im,2);
height = size(im,1);

if length(size(im))==3
    color=1;
    numbits=width*height*8*3;
else
    color=0;
    numbits=width*height*8;
end

% figure(1);cla;
% if color==1
%     imshow(uint8(original));
% else
%     % For grayscale images: scale and give the colormap
%     imshow(im);
% end
% title(sprintf('Originalna slika %s.',file));


if color==1  
    
% im = RGB2YUV(im,width,height);  
    im = round(im*255)-127;
    for k=1:3
        yrange = 1:1*width;          
        xrange = 1:1*width;               
        [block,dwtcoef{1,1,k}]=fdwt(im(xrange,yrange,k),N,width);
        reconstructedtile{1,1,k} = bdwt(block,N,width);
    end
    
    for k=1:3
        yrange = 1:1*width;
        xrange = 1:1*width;     
        reconstructedim(xrange,yrange,k) = reconstructedtile{1,1,k};
        totalcoef(xrange,yrange,k) = dwtcoef{1,1,k}/max(max(abs(dwtcoef{1,1,k})));
    end
    out=(reconstructedim+127)/255;
%     out = YUV2RGB(out,width,height);
%     totalcoef = YUV2RGB(totalcoef,width,height);
    
%     figure(4);cla;
%     imshow(totalcoef);
%     
%     figure(5);cla;
%     imshow(uint8(out));
%     title('Reconstructed image.');
    
else
    im = round(im*255)-127;

    xrange = 1:1*width;
    yrange = 1:1*width;
    [block,dwtcoef{1,1}]=fdwt(im(xrange,yrange),N,width);
    [reconstructedtile{1,1}] = bdwt(block,N,width);

    xrange = 1:1*width;
    yrange = 1:1*width;
            
    reconstructedim(xrange,yrange) = reconstructedtile{1,1};
    totalcoef(xrange,yrange) = dwtcoef{1,1};

    out=reconstructedim+127;
    
%     figure(4);cla; 
%     imshow(totalcoef);     
%     title('Wavelet koeficijenti.');
%     out = out/255
%     figure(5);cla;
%     imshow(out); 
%     title('Reconstructed image.');
      
end 

function [block,dwtcoef]=fdwt(im,N,width)

    [C,S]=wavedec2(im,N,'bior4.4');
    block{N,4}=appcoef2(C,S,'bior4.4',N);
    
    sz = width*(2^-N);
    dwtcoef= wkeep(block{N,4},[sz,sz]);
    for i=N:-1:1
        sz = width*(2^-i);
        [block{i,1},block{i,2},block{i,3}]= detcoef2('all',C,S,i);
        dwtcoef=[dwtcoef, wkeep(block{i,1},[sz,sz]);
        wkeep(block{i,2},[sz,sz]),  wkeep(block{i,3},[sz,sz])];
    end

function im=bdwt(block,N,width)

    C = block{N,4}(:);
    sz = size(block{N,4},1);
    S(1,:)= [sz,sz];
    
    for i=N:-1:1
        sz = size(block{i,1},1);
        S(N-i+2,:)= [sz,sz];
        C = [C;block{i,1}(:);block{i,2}(:);block{i,3}(:)];
    end
    S(N+2,:)=[width,width];
    
    im = waverec2(C,S,'bior4.4');    
    
% %
% function x=YUV2RGB(x,width,height)
%     rgbmat=[1,  0,          1.402; 
%             1,  -0.34414,   -0.71414;
%             1,  1.772,      0];
%        
%     % apply transformation to all RGB tripples
%     for i=1:width
%         for j=1:height
%             x(j,i,:) = rgbmat*(squeeze(x(j,i,:))-[0;.5;.5]);
%         end
%     end
%     
%     %
% %
% function x=RGB2YUV(x,width,height)
%     yuvmat = [  0.299   ,0.587      ,0.114; 
%                 -0.1687 ,-0.3313    ,0.5; 
%                 0.5     ,-0.4187    ,-0.0813];
% 
%     % apply transformation to all RGB tripples
%     for i=1:width
%         for j=1:height
%             x(j,i,:) = yuvmat*squeeze(x(j,i,:))+[0;.5;.5];
%         end
%     end
    