K-means via PCA

I want cluster the H&E image to segment cellular, stromal, and luminal region. First I converted RGB image matrix into n*3 matrix. each columns represent RGB. then apply PCA so I got reduce dimensional matrix of 3*3 size. I use this matrix as initialize the cluster. My code works fine but it take long time i.e. 30 to 60 minuets. Any one can modified the code to optimize and speedup the process. The code is as follows:

%k is 3 here
function im = kmeansclusterviapca(filename,k)
  in_image = imread(filename);
  rows  = size(in_image,1);
  cols = size(in_image,2); 
  Data = zeros(rows*cols,3);
  R = in_image(:,:,1);
  G = in_image(:,:,2);
  B = in_image(:,:,3);
  l = 1;
  for i = 1:rows
      for j=1:cols
         Data(l,1) = R(i,j);
         Data(l,2) = G(i,j);
         Data(l,3) = B(i,j);
         l = l + 1;
      end
  end
  clearvars R G B;
  [rows cols] = size(Data);
  m = initializeClusterVectorsUsingPrincipalComponent(Data,k);
  prevM = zeros(k, 3);
  iterNo = 0;
  distance = zeros(size(m,1),1);
  while (iterNo <= 500)
      iterNo = iterNo + 1
      prevM = m;
      BB = zeros(rows, k);
      %----------------------------------------------------------------
      p = size(m,1);
      for i = 1 : rows
          for j = 1 : p
              distance(j) = 0;
              for t = 1 : cols
                  distance(j) = distance(j) + (m(j,t)-Data(i,t))^2;
              end
          end
          [minval minIndex] = min(distance);
          BB(i, minIndex)= 1;
      end
      %---------------------------------------------------------------
      n = zeros(size(m,1),1);
      m=m.* 0;
      for t = 1 :size(BB,1)
          for i = 1 : size(BB,2)
              if (int16(BB(t, i)) == 1) 
                  for j = 1 :size(m,2)
                      m(i, j) = m(i, j) + Data(t, j);
                  end
                  n(i,1)= n(i,1) + 1;
              end
          end
      end
      for i = 1: size(n,1)
          if (n(i, 1) > 0) 
              for j = 1 :size(m,2)
                  m( i, j)= m(i, j) / n(i, 1);
              end
          else       %randomly assign new clustering vector
              t = int16(rnd) *size(Data,1);
              for j = 1 :size( m,2)
                  m(i, j)= Data(t, j);
              end
          end
      end
      %---------------------------------------------------------
      maxDeltaM = maxAbsoluteDifferenceBetweenMatrices(m,prevM);
      if (maxDeltaM <= 0.0001) 
          break
      end
  end
  clearvars BB;
  p = size(m,1);
  distance = zeros(p);
  rows = size(Data,1);
  labels = zeros(rows);

    for i = 1: size(Data,1)
        for j = 1: p
            distance(j) = 0;
            for t = 1:size(Data,2)
                distance(j) = distance(j) + (m(j, t)-Data(i, t))^2;
            end
        end
        [minval minIndex] = min(distance);
        labels(i)= minIndex;
    end
    p = 1;
    rows  = size(in_image,1);
    cols = size(in_image,2); 
    out = zeros(rows,cols);
    for i = 1: rows
        for j = 1: cols
            out(i, j)= labels(p);
            p = p + 1;
        end
    end
    clearvars labels;
    if ((m(1, 1) > m(2, 1)) && (m(1, 1) > m(3, 1))) 
        lm_id = 0;
        if (m(2, 1) > m(3, 1))
            st_id = 1;
            nc_id = 2;
        else
            st_id = 2;
            nc_id = 1;
        end
    else
    if ((m(2, 1) > m(1, 1)) && (m(2, 1) > m(3, 1)))
        lm_id = 1;
        if (m(1, 1) > m(3, 1)) 
            st_id = 0;
            nc_id = 2;
        else
            st_id = 2;
            nc_id = 0;
        end
    else
        lm_id = 2;
        if (m(1, 1) > m(2, 1)) 
            st_id = 0;
            nc_id = 1;
        else
            st_id = 1;
            nc_id = 0;
        end
    end
    end
    for i = 1: rows
        for j = 1: cols
            if (out(i, j) == lm_id) 
                out(i, j)= 2;
            else
                if (out(i, j) == nc_id) 
                out(i, j)= 0;
                else
                    if (out(i, j) == st_id) 
                        out(i, j)= 1;
                    end
                end 
            end
        end
    end
    im = out;
end
-----------------------------------------------------
function m=initializeClusterVectorsUsingPrincipalComponent(Data , k )
    CR = corrcoef(Data);
    [V,D] = eig(CR);
    newdata = Data * V;
    maxentry = max(newdata(:,1));
    minentry = min(newdata(:,1));
    d = (maxentry - minentry)/k;
    c = zeros(k+1);
    c(1) = minentry;
    for i = 2 : k+1
        c(i) = c(i - 1) + d;
    end
    No = zeros(1,k);
    newM = zeros(k, 3);
    [rows cols] = size(newdata);
    for i = 1 : rows
        for j = 1 : k 
            if ((newdata(i, 1) >= c(j)) && (newdata(i, 1) < c(j + 1))) 
                for l = 1 :cols
                    newM(j, l)= newM(j, l) + newdata(i, l);
                end
                No(j) = No(j) + 1;
                break;
            end
        end
    end
    for i = 1 : k
        for j = 1 : cols
            newM(i, j)= newM(i, j) / No(i);
        end
    end
    temp = inv(V);
    m = newM*temp;
end
--------------------------------------------------------------------------
function maxAbsEntry =maxAbsoluteDifferenceBetweenMatrices(A,B)

    Entry = abs(A(1,1) - B(1,1));
    for i = 1 :size(A,1)
        for j = 1 : size(A,2)
            if (Entry < abs(A(i, j) - B(1, 1)))
                Entry = abs(A(i, j) - B(1, 1));
            end
        end
    end
    maxAbsEntry=Entry;
end