Bag-of-Words for SIFT using vl_feat

Question

I am intending to quantize the SIFT features I have intended to extract from my image set using the BOW. I know how to extract the SIFT features from one image using the vl_sift function, but I don't know how to perform the kmeans operation on all the features extracted from all the images in the big set I have. This is similar to my case except that I don't know what the data he used referred to and I didn't know how to collect all the features extracted from all the images to do a kmeans on them. this is the code I have so far. It extracts the SIFT features from each image but doesn't save them anywhere.

imgds = imageSet('E:\Flower Classification\10\seg','recursive')
for index = 1:imgds.Count
    for i = 1:40
        I = read(imgds(index),i);
        II=single(rgb2gray(I));
        [~,d]= vl_sift(II);
        ???
    end
end

I appreciate any suggestions. I need the following:

1. collect all the features in one list
2. perform kmeans on all these features (there is a function called vl_kmeans but it takes a set of data and I couldn't find in the vl_feat documentation anything about how to form this data)
3. create a vector for each image to represent it. (which is also called the features histogram, so instead of representing the image using a set of features it is represented by the number of features that lie in each kmeans center)