function bitmap_fast(filename,nxx,pxx,nyy,pyy)%, knot_vectorx, knot_vectory)
  
tic;

%number of basis functions
compute_nr_basis_functions = @(knot_vector,p) size(knot_vector, 2) - p - 1;

knot_vectorx = simple_knot(nxx,pxx);
knot_vectory = simple_knot(nyy,pyy);

%RGB format
X = imread(filename);

R = X(:,:,1);
G = X(:,:,2);
B = X(:,:,3);

ix = size(X,1);
iy = size(X,2);

px = compute_p(knot_vectorx);
py = compute_p(knot_vectory);

elementsx = number_of_elements(knot_vectorx,px);
elementsy = number_of_elements(knot_vectory,py);

nx = number_of_dofs(knot_vectorx,px);
ny = number_of_dofs(knot_vectory,py);

A = sparse(nx*ny,nx*ny);
FR = zeros(nx*ny,1);
FG = zeros(nx*ny,1);
FB = zeros(nx*ny,1);

init=toc
tic;

splinex = zeros(elementsx,nx,px+1);
spliney = zeros(elementsy,ny,py+1);


for ex = 1:elementsx;
  % suport of basis functions over an element in x directin
  [xl,xh] = dofs_on_element(knot_vectorx,px,ex);
  % elemental span along x axis
  [ex_bound_l,ex_bound_h] = element_boundary(knot_vectorx,px,ex);
  % quadrature points along x axis
  qpx = quad_points(ex_bound_l,ex_bound_h,px+1);
  % quadrature weights along x axis
  qwx = quad_weights(ex_bound_l,ex_bound_h,px+1);
  % loop through non-zero functions
  for bi = xl:xh
    % loop through quadrature points and weights
    for iqx = 1:size(qpx,2)
      splinex(ex,bi,iqx)= compute_spline(knot_vectorx,px,bi,qpx(iqx));
    endfor
  endfor
endfor

for ey = 1:elementsy;
  % suport of basis functions over an element in y directin
  [yl,yh] = dofs_on_element(knot_vectory,py,ey);
  % elemental span along y axis
  [ey_bound_l,ey_bound_h] = element_boundary(knot_vectory,py,ey);
  % quadrature points along y axis
  qpy = quad_points(ey_bound_l,ey_bound_h,py+1);
  % quadrature weights along y axis
  qwy = quad_weights(ey_bound_l,ey_bound_h,py+1);
  % loop through non-zero functions
  for bi = yl:yh
    % loop through quadrature points and weights
    for iqy = 1:size(qpy,2)
      spliney(ey,bi,iqy)= compute_spline(knot_vectory,py,bi,qpy(iqy));
    endfor
  endfor
endfor

init_splines=toc
tic;

%integrals of B^x_i(x) B^y_j(y) B^x_k(x) B^y_l(y)
%(i,k=1,...,Nx; j,l=1,...,Ny)
%loop through elements along x
for ex = 1:elementsx;
  % elemental span along x axis
  [xl,xh] = dofs_on_element(knot_vectorx,px,ex);
  % left and right border of element along x 
  [ex_bound_l,ex_bound_h] = element_boundary(knot_vectorx,px,ex);
  %loop through elements along y
  for ey = 1:elementsy
    % elemental span along y axis
    [yl,yh] = dofs_on_element(knot_vectory,py,ey);
    % left and right border of element along y 
    [ey_bound_l,ey_bound_h] = element_boundary(knot_vectory,py,ey);
    % jacobian = element size
    Jx = ex_bound_h - ex_bound_l;
    Jy = ey_bound_h - ey_bound_l;
    J = Jx*Jy;
    % quadrature points along x axis
    qpx = quad_points(ex_bound_l,ex_bound_h,px+1);
    % quadrature weights along x axis
    qwx = quad_weights(ex_bound_l,ex_bound_h,px+1);
    % quadrature points along y axis
    qpy = quad_points(ey_bound_l,ey_bound_h,py+1);
    % quadrature weights along y axis
    qwy = quad_weights(ey_bound_l,ey_bound_h,py+1);
    % loop through non-zero functions over an element
    for bi = xl:xh
      for bj = yl:yh
        for bk = xl:xh
          for bl = yl:yh
            % loop through quadrature points
            for iqx = 1:size(qpx,2)
              for iqy = 1:size(qpy,2)
                % B^x_k(x)
                funk = splinex(ex,bk,iqx);
                % B^y_l(y)
                funl = spliney(ey,bl,iqy);
                % B^x_i(x)
                funi = splinex(ex,bi,iqx);
                % B^y_j(y)
                funj = spliney(ey,bj,iqy);
                % B^x_i(x) B^y_j(y) B^x_k(x) B^y_l(y)
                fun = funi*funj*funk*funl;
                
                % Integrals of B^x_i(x) B^y_j(y) B^x_k(x) B^y_l(y)
                % (i,k=1,...,Nx; j,l=1,...,Ny)
                int = fun*qwx(iqx)*qwy(iqy)*J;
                if (int!=0)
                  A((bj-1)*nx+bi,(bl-1)*nx+bk) = A((bj-1)*nx+bi,(bl-1)*nx+bk) + int;
                endif
              endfor
            endfor
          endfor
        endfor
      endfor
    endfor
  endfor
endfor

lhs=toc

tic;
            
%Integra;s BITMAPA(x,y) B^x_k(x) B^y_l(y)
%loop through elements along x
for ex = 1:elementsx;
  % elemental span along y axis
  [xl,xh] = dofs_on_element(knot_vectorx,px,ex);
  % left and right border of element along x 
  [ex_bound_l,ex_bound_h] = element_boundary(knot_vectorx,px,ex);
  %loop through elements along y
  for ey = 1:elementsy
    % elemental span along y axis
    [yl,yh] = dofs_on_element(knot_vectory,py,ey);
    % borders of element 
    [ey_bound_l,ey_bound_h] = element_boundary(knot_vectory,py,ey);
    % jacobian = element size
    Jx = ex_bound_h - ex_bound_l;
    Jy = ey_bound_h - ey_bound_l;
    J = Jx * Jy;
    % quadrature points along x axis
    qpx = quad_points(ex_bound_l,ex_bound_h,px+1);
    % quadrature points along y axis
    qpy = quad_points(ey_bound_l,ey_bound_h,px+1);
    % quadrature weights along x axis
    qwx = quad_weights(ex_bound_l,ex_bound_h,py+1);
    % quadrature weights along y axis
    qwy = quad_weights(ey_bound_l,ey_bound_h,py+1);
    %loop through non-zero functions
    for bk = xl:xh
      for bl = yl:yh  
        %loop through quadrature points
        for iqx = 1:size(qpx,2)
          for iqy = 1:size(qpy,2)
            % B^x_k(x)
            funk = splinex(ex,bk,iqx);
            % B^y_l(y)
            funl = spliney(ey,bl,iqy);
            % Integral BITMAPA(x,y) B^x_k(x) B^y_l(y) over RGB colors
            intR = funk*funl*qwx(iqx)*qwy(iqy)*J*bitmp(R,qpx(iqx),qpy(iqy));
            intG = funk*funl*qwx(iqx)*qwy(iqy)*J*bitmp(G,qpx(iqx),qpy(iqy));
            intB = funk*funl*qwx(iqx)*qwy(iqy)*J*bitmp(B,qpx(iqx),qpy(iqy));
            FR((bl-1)*nx+bk) = FR((bl-1)*nx+bk) + intR;
            FG((bl-1)*nx+bk) = FG((bl-1)*nx+bk) + intG;
            FB((bl-1)*nx+bk) = FB((bl-1)*nx+bk) + intB;
          endfor
        endfor
      endfor
    endfor
  endfor
endfor

rhs=toc
tic;

% solution of three systems of equations
RR=A\FR;
GG=A\FG;
BB=A\FB;

%solution reconstruction

%preparing of matrices
R1 = zeros(ix,iy);
G1 = zeros(ix,iy);
B1 = zeros(ix,iy);

factor=toc
tic;

funx_tab = zeros(nx,ix);
funy_tab = zeros(ny,iy);

% loop through functions
for bi = 1:nx
% loop through pixels covering the support of functions
  for i=xx(knot_vectorx(bi)):xx(knot_vectorx(bi+px+1))
% coordinates transformation [1-length] -> [0-1]
    ii = (i-1)/(ix-1);
    % B^x_i(x)
    funx_tab(bi,i) = compute_spline(knot_vectorx,px,bi,ii);
  endfor
endfor


% loop through functions
  for bj = 1:ny
% loop through pixels covering the support of functions
    for j=yy(knot_vectory(bj)):yy(knot_vectory(bj+py+1))  
% coordinates transformation [1-height] -> [0-1]
      jj = (j-1)/(iy-1);
      % B^y_j(y)
      funy_tab(bj,j) = compute_spline(knot_vectory,py,bj,jj);
  endfor
endfor

preprocess=toc
tic;

% loop through functions
for bi = 1:nx
  for bj = 1:ny
% loop through pixels covering the support of functions
    for i=xx(knot_vectorx(bi)):xx(knot_vectorx(bi+px+1))
      for j=yy(knot_vectory(bj)):yy(knot_vectory(bj+py+1))  
        % B^x_i(x)
        funi = funx_tab(bi,i);%compute_spline(knot_vectorx,px,bi,ii);
        % B^y_j(y)
        funj = funy_tab(bj,j);%compute_spline(knot_vectory,py,bj,jj);
        ff = funi*funj;
        R1(i,j) = R1(i,j) + floor(ff*RR((bj-1)*nx+bi));
        G1(i,j) = G1(i,j) + floor(ff*GG((bj-1)*nx+bi));
        B1(i,j) = B1(i,j) + floor(ff*BB((bj-1)*nx+bi));
      endfor
    endfor
  endfor
endfor

RGB=X;
RGB(:,:,1) = R1;
RGB(:,:,2) = G1;
RGB(:,:,3) = B1;

rebuild=toc

imshow(RGB);





function resx=xx(x)
  resx = floor((ix-1)*x+1);
endfunction



function resy=yy(y)
  resy = floor((iy-1)*y+1);
endfunction


function val=bitmp(M,x,y)
  val = zeros(size(x));
  for i=1:size(x,1)
    for j=1:size(x,1)
      val(i,j)=M(xx(x(1,i)),yy(y(1,j)));
    endfor
  endfor
endfunction


%the function computes the order of function
function p=compute_p(knot_vector)

%first knot point
  initial = knot_vector(1);
%length of knot_vector
  kvsize = size(knot_vector,2);
  p = 0;

%checking the repetitions of the first knot point
  while (p+2 <= kvsize) && (initial == knot_vector(p+2))
    p = p+1;
  endwhile
  
  return
endfunction
  
  
  
  
  
%function checking the correctness of knot_vector
function t=check_sanity(knot_vector,p)

  initial = knot_vector(1);
  kvsize = size(knot_vector,2);

  t = true;
  counter = 1;

%warning of number of repetitions does not correspond with order 
  for i=1:p+1
    if (initial != knot_vector(i))
%return false
      t = false;
      return
    endif
  endfor

%warning if too many repetitions inside
  for i=p+2:kvsize-p-1
    if (initial == knot_vector(i))
      counter = counter + 1;
      if (counter > p)
%return false
        t = false;
        return
      endif
    else
      initial = knot_vector(i);
      counter = 1;
    endif
  endfor

  initial = knot_vector(kvsize);
  
%warning of number of repetitions does not correspond with order 
  for i=kvsize-p:kvsize
    if (initial != knot_vector(i))
%return false
      t = false;
      return
    endif
  endfor
  
%if the next knot point is different from previous ones
  for i=1:kvsize-1
    if (knot_vector(i)>knot_vector(i+1))
%return false
      t = false;
      return
    endif
  endfor

  return

endfunction



%the functions computes B-splines recursively according to the Cox-de-Boor
function y=compute_spline(knot_vector,p,nr,x)
  
%function (x-x_i)/(x_{i-p}-x_i)
  fC= @(x,a,b) (x-a)/(b-a);
%function (x_{i+p+1}-x)/(x_{i+p+1}-x_{i+1})
  fD= @(x,c,d) (d-x)/(d-c);
  
%x_i  
  a = knot_vector(nr);
%x_{i-p} 
  b = knot_vector(nr+p);
%x_{i+1}
  c = knot_vector(nr+1);
%x_{i+p+1}
  d = knot_vector(nr+p+1);

%linear function for p=0
  if (p==0)
    y = 0 .* (x < a) + 1 .* (a <= x & x <= d) .+ 0 .* (x > d);
    return
  endif

%B_{i,p-1}  
  lp = compute_spline(knot_vector,p-1,nr,x);
%B_{i+1,p-1}
  rp = compute_spline(knot_vector,p-1,nr+1,x);
  
%(x-x_i)/(x_{i-p)-x_i)*B_{i,p-1}  
  if (a==b)
%if knot points are repeated 
    y1 = 0 .* (x < a) + 1 .* (a <= x & x <= b) + 0 .* (x > b);
  else
    y1 = 0 .* (x < a) + fC(x,a,b) .* (a <= x & x <= b) + 0 .* (x > b);
  endif
  
%(x_{i+p+1}-x)/(x_{i+p+1)-x_{i+1})*B_{i+1,p-1}
  if (c==d)
%if knot points are repeated 
    y2 = 0 .* (x < c) + 1 .* (c < x & x <= d) + 0 .* (d < x);
  else
    y2 = 0 .* (x < c) + fD(x,c,d) .* (c < x & x <= d) + 0 .* (d < x);
  endif
  
  y = lp .* y1 .+ rp .* y2;
  return
  
endfunction

  
  
% number of elements in knot vector
function n = number_of_elements(knot_vector,p)
  
  initial = knot_vector(1);
  kvsize = size(knot_vector,2);
  n = 0;
  
  for i=1:kvsize-1
    if (knot_vector(i) != initial)
      initial = knot_vector(i);
      n = n+1;
    endif
  endfor
  
endfunction
  
% simple knot vector 
function knot = simple_knot(elems, p)
  pad = ones(1, p);
  knot = [0 * pad, 0:elems, elems * pad];
  knot = knot/elems;
endfunction

% number of basis functions
function n = number_of_dofs(knot,p)
  n = length(knot) - p - 1;
endfunction


function first = first_dof_on_element(knot_vector,p,elem_number)
 [l,h] = element_boundary(knot_vector,p,elem_number);
 first = lookup(knot_vector, l) - p;
endfunction
  
function [low,high] = element_boundary(knot_vector,p,elem_number)
  initial = knot_vector(1);
  kvsize = size(knot_vector,2);
  k = 0;
  low=0;
  high=0;
  
  for i=1:kvsize
    if (knot_vector(i) != initial)
      initial = knot_vector(i);
      k = k+1;
    endif
    if (k == elem_number)
      low = knot_vector(i-1);
      high = knot_vector(i);
      return;
    endif
  endfor
endfunction

% computes functions of non-zero span over an element
function [low,high] = dofs_on_element(knot_vector,p,elem_number)
  low = first_dof_on_element(knot_vector,p,elem_number);
  %since we have exactly p+1 non-zero functions over an element
  high = low + p;
endfunction

% Row vector of points of the k-point Gaussian quadrature on [a, b]
function xs = quad_points(a, b, k)
  % points mapping
  map = @(x) 0.5 * (a * (1 - x) + b * (x + 1));
  switch (k)
    case 1
      xs = [0];
    case 2
      xs = [-1/sqrt(3), ...
             1/sqrt(3)];
    case 3
      xs = [-sqrt(3/5), ...
             0,         ...
             sqrt(3/5)];
    case 4
      xs = [-sqrt((3+2*sqrt(6/5))/7), ...
             sqrt((3-2*sqrt(6/5))/7), ...
             sqrt((3-2*sqrt(6/5))/7), ...
             sqrt((3+2*sqrt(6/5))/7)];
    case 5
      xs = [-1/3*sqrt(5+2*sqrt(10/7)), ...
            -1/3*sqrt(5-2*sqrt(10/7)), ...
             0,                        ...
             1/3*sqrt(5-2*sqrt(10/7)), ...
             1/3*sqrt(5+2*sqrt(10/7))];
    otherwise
      xs = [-1/3*sqrt(5+2*sqrt(10/7)), ...
            -1/3*sqrt(5-2*sqrt(10/7)), ...
             0,                        ...
             1/3*sqrt(5-2*sqrt(10/7)), ...
             1/3*sqrt(5+2*sqrt(10/7))];
  endswitch
  xs = map(xs);
endfunction

% Row vector of weights of the k-point Gaussian quadrature on [a, b]
function ws = quad_weights(a, b, k)
  switch (k)
    case 1
      ws = [2];
    case 2
      ws = [1, 1];
    case 3
      ws = [5/9, ...
            8/9, ...
            5/9];
    case 4
      ws = [(18-sqrt(30))/36, ...
            (18+sqrt(30))/36, ...
            (18+sqrt(30))/36, ...
            (18-sqrt(30))/36];
    case 5
      ws = [(322-13.0*sqrt(70))/900, ...
            (322+13.0*sqrt(70))/900, ...
            128/225,                 ...
            (322+13.0*sqrt(70))/900, ...
            (322-13.0*sqrt(70))/900];
    otherwise
      ws = [(322-13.0*sqrt(70))/900, ...
            (322+13.0*sqrt(70))/900, ...
            128/225,                 ...
            (322+13.0*sqrt(70))/900, ...
            (322-13.0*sqrt(70))/900];
  endswitch
  % Gaussian quadrature is defined on [-1, 1], we use [a, b]
  %ws = ws * (b-a)/2;
endfunction



endfunction