Mixing indices with subscripts or how can I capture k-means features faster

I'm trying to do something I think is relatively simple. I'm going to perform a k-means clustering (via imsegkmeans) on a grayscale image. For this I'd like the 9 pixels centered on the current pixel (y,x) to be features so I'm taking the 2D image and making it 3D by adding 9 more layers, one for each neighboring pixel (I know, the current pixel is now there twice at :,:,1 and :,:,5 but we'll fix that later). I have two methods that can do this correctly:

    featureSet2 = zeros(size(pImage,1),size(pImage,2),9);
    incr = size(pImage,1);
    for imagePos = incr+1:(incr * (incr -1)) + 1
        
        % skip the outside boarder (mod() == 0 is bottom edge, mod() == 1
        % is top edge). The incrementor takes care of left and right sides
        if mod(imagePos,incr) > 1
            [y,x] = ind2sub(size(featureSet2),imagePos);
            featureSet2(y,x,:) = reshape(pImage(imagePos+offsets),[9,1]);
        end
    end

and

    featureSet = zeros(size(pImage,1),size(pImage,2),9);
    for x = 2:s(2)-1 % s is the size of pImage
        for y = 2:s(1)-1
            coord = x * s(1) + y;
            % offset is a 3x3 matrix of indicies offsets for the surrounding eight pixels
            featureSet(y,x,:) = reshape(pImage(coord+offsets),[9,1]);
        end
    end

The top one seems more Matlab-ish with only one loop and using matrix indexing yet it runs slower than the bottom version (3.196324 seconds versus 1.980120 seconds). I suspect it's the ind2sub being executed each cycle of the loop but I'm not sure of another way to index the storage location, featureSet2.

Thoughts on how to speed this up? I have large number of images to do sequentially and I don't have a GPU so seconds count.

The basis of the code above is loosly based on the gabor filter example here:

https://www.mathworks.com/help/images/ref/imsegkmeans.html

Thanks!

[Walter Roberson]

Since seconds count, move size(featureSet2) out of the for loop.

[John Komp]

Grrr. Can't believe I missed that...

Here's the new code:

    featureSet2 = zeros(size(pImage,1),size(pImage,2),9);
    featSize = size(featureSet2);
    incr = size(pImage,1);
    for imagePos = incr+1:(incr * (incr -1)) + 1
        
        % skip the outside boarder (mod() == 0 is bottom edge, mod() == 1
        % is top edge). The incrementor takes care of left and right sides
        if mod(imagePos,incr) > 1
            [y,x] = ind2sub([featSize(1),featSize(2)],imagePos);
            featureSet2(y,x,:) = reshape(pImage(imagePos+offsets),[9,1]);
        end
    end

Now the run time for this loop is down to 1.7s. So better but not great yet. Thanks for pointing out the unneeded extra processing.

I wish there was a way that I could load the z direction vector [9,1] into a 2D array at each point in that 2D array and then reshape the array to 3D where that vector at each point becomes the values in the z direction.

[Walter Roberson]

Tweaking. First make sure that offsets is a column vector. Then

    featureSet2 = zeros(size(pImage,1),size(pImage,2),9);
    featSize = size(featureSet2);
    featSize12 = featSize(1:2);
    incr = size(pImage,1);
    for imagePos = incr+1:(incr * (incr -1)) + 1    
        % skip the outside boarder (mod() == 0 is bottom edge, mod() == 1
        % is top edge). The incrementor takes care of left and right sides
        if mod(imagePos,incr) > 1
            [y,x] = ind2sub(featSize12, imagePos);
            featureSet2(y,x,:) = pImage(imagePos+offsets);
        end
    end

[John Komp]

Thanks for the tweak. That actually made the double loop a bit faster. I'll probably have to keep that instead of the single loop/ind2sub version as the double loop is still a half second faster (~33%).