function cvd ( m, n, k )

%% CVD carries out 20 steps of the centroidal Voronoi iteration.
%
%  Parameters:
%
%    Input, integer M, N, the number of pixels in the vertical and horizontal
%    directions.
%
%    Input, integer K, the number of regions to generate.
%
  p = rand ( k, 2 );
  pc(1:k) = uint8 ( 1 : k );
  a = zeros ( m, n );

  for it = 1 : 20

    p_new = zeros(k,2);
    p_count = zeros(k);

    for i = 1 : m
      y = ( m  - i ) / ( m - 1 );
      for j = 1 : n
        x = ( j - 1 ) / ( n - 1 );
        dmin = Inf;
        kkmin = 0;
        for kk = 1 : k
          d = norm ( p(kk,1:2) - [ x, y ] );
          if ( d < dmin )
            dmin = d;
            kkmin = kk;
          end
        end
        p_new(kkmin,1) = p_new(kkmin,1) + x;
        p_new(kkmin,2) = p_new(kkmin,2) + y;
        p_count(kkmin) = p_count(kkmin) + 1;
        a(i,j) = pc(kkmin);
      end
    end

    for kk = 1 : k
      p(kk,1:2) = p_new(kk,1:2) / p_count(kk);
    end

    imagesc ( a );
    pause

  end
  
  return
end
