How to work with PCA on BW Image?

Hello,

I've been trying to use PCA on my BW image with no success at all. What I'm trying to achieve is for each blob in the image, retrieve is principal component axis (which is different from the major axis length given by regionprops), draw it on top of the current blob and then calculate the principal axis angle in relation to an horizontal line. The principal component axis would give me a major lead on calculating the line that better represents my blob and with that, calculating its angle should be rather simple, which is why I thought of using PCA.

Here's my code so far:

imshow(edges);
hold on
[labeledImage, numberOfBlobs] = bwlabel(edges);
for blob = 1 : numberOfBlobs
    thisBlob = ismember(labeledImage, blob);
    [rows, cols] = find( thisBlob );
    
    A = [cols, rows];
    uCols = unique(cols);
    uRows = unique(rows);
    
    colsLength = uCols(end)-uCols(1);
    rowsLength = uRows(end)-uRows(1); 
    
    C=cov(A);
    vtxmean=mean(A);

    [eVect,~]=eig(C);
    [~,ind] = sort(diag(eVect));
    
    %https://stackoverflow.com/questions/27320142/oriented-bounding-box-is-misshapen-and-the-wrong-size-in-opengl
    
    [xa, ya] = deal(eVect(1,ind(end)), eVect(2,ind(end)));
    angle = cart2pol(xa, ya)*180/pi; %from https://www.mathworks.com/matlabcentral/fileexchange/7432-pcaangle

    currLine.point1(1) = vtxmean(1); %x
    currLine.point1(2) = vtxmean(2); %y
    
    currLine.point2(1) = vtxmean(1)+eVect(2,1)*colsLength/2;   %x 
    currLine.point2(2) = vtxmean(2)+eVect(2,2)*rowsLength/2;   %y 

    plot([currLine.point1(1), currLine.point2(1)], ...
        [currLine.point1(2), currLine.point2(2)], 'b','LineWidth',2);
end
hold off

Although this seems close to the expected result, the line draws are always off. Their center seems correct, but the final point is just wrong, as seen below. Also, the angle I'm getting is 77.2428 which also doesn't make that much sense since this case clearly has more than 90 degrees in respect to the horizontal line. The angle always have to be calculated counter clockwise.

Some more thoughts I had before this approach:

1. Regionprops Orientation doesn't work because the ellipse calculated sometimes won't be aligned with the line segment itself.
2. Feret Diamater won't work aswell since that would give me the biggest line possible inside the line segment, which is different from having a line that better represents the segment.
3. The following code got close to work, but I figured that cases where the angle was 90 degrees would make it fail.

    [p, S, mu] = polyfit(rows, cols, 2);
    cols2 = polyval(p, rows, S, mu);
    %plot(cols2, rows, 'r');
    currLine.point1(1) = min(cols2);
    currLine.point1(2) = min(rows);
    currLine.point2(1) = max(cols2);
    currLine.point2(2) = max(rows);

Angle was calculated with

    function angle = CalcAngle(v1, v2)
        lengthH = sqrt(power(v2(1)-v1(1),2) + power(v2(2)-v1(2),2));
        lengthX = v2(1)-v1(1);
        angle = asind(lengthX/lengthH)+90;
    end

Or even

    function angle = CalcAngle(v1, v2)
        lengthX = v2(1)-v1(1);
        lengthY = v2(2)-v1(2);
        angle = atan2d(lengthY, lengthX);
        if( angle < 0 )
            angle = angle + 180;
        end
    end

Can someone help me out?

Answer 1

[Image Analyst]

The angle is given by the 'Orientation' parameter, so if you want the angle, just ask for the orientation.

If you want a line through the blob then you can determine the orientation of the blob (mostly vertical or mostly horizontal) and then just feed all the points in it (from PixelList property) into polyfit().

[Thiago Motta]

You are absolutely right. The values given by RegionProps.Orientation are indeed the ones that seems most accurate so far. I dont know why it didn't work before.

Polyfit is also giving an accurate result.

Thank you for your comments.