SVD line fitting or ransac line fitting in multidimensionl image

i have a multidimensional image of size 1024*512*128. For each slice(1024*512), I have single point from the mid slice of an image say from slice 40 to 128. So, i have 89 points in my multidimensional (volumetric) image. how can i fit the straight line using svd/ ransac or any other method and then extrapolate that line to the end of the volumetric image?

Answer 1

[Matt J]

This will give you the equations Aeq*x=beq for the line,

       Q=[X(:),Y(:),Z(:)];
       m=mean(Q,1).';
       Q=Q-m.';

       [~,~,Z]=svd(Q,0);

   Aeq=Z(:,2:3).';
   beq=Aeq*m;

Because you know the vertices of the volume's bounding box (put them in the rows of a matrix V), you can then use intersectionHull (download here) to find the intersection points of the line with the box,

   I=intersectionHull('vert',V,'lcon',[],[],Aeq,beq);

[prashant singh]

i have tried this approach but i am getting this error for svd function that Input must be 2-D. My image size is 1024 *512 *128. i can't use the transpose directly, so i have used permute function.

Q = array3d;
m=mean(permute(Q,[2 1 3]),1);
Q=Q-permute(m,[2 1 3]);
[~,~,Z]=svd(Q,0);

Error using svd
Input must be 2-D.

[Matt J]

The error is in the line

Q=array3d;

It's not what I directed you to input.

[prashant singh]

my array3d is my multidimensional data. I am generating it from another code that's why i am assigning it to another variable to work on so, I don't have to generate it again. Now my Q becomes my volumetric image of size 1024*512*128.

[Matt J]

Q should be an 89x3 matrix whos rows are the coordinates of your 89 points.

[prashant singh]

yep that was dumb on my part. Now i have put the correct matrix , i am getting the error on finding the intersection hull part.

Error using lcon2vert (line 121)
A and Aeq must have the same number of columns if both non-empty

 Error in intersectionHull (line 132)
   [V,nr,nre]=lcon2vert(A,b,Aeq,beq,TOL);

Uploading my mat file which I use as an input.

[Matt J]

I don't have that problem. What did you provide for V? I used the following

V =

           1           1           1
        1024           1           1
           1         512           1
        1024         512           1
           1           1         128
        1024           1         128
           1         512         128
        1024         512         128

and I obtain as follows

I=intersectionHull('vert',V,'lcon',[],[],Aeq,beq)

    >> I.vert

    ans =

      242.7437  512.0000   61.3366
      211.6177    1.0000  109.9629

[prashant singh]

Thanks, Matt. Is there function to plot this Hull on my volumetric image to visualise the result?

[Matt J]

You're welcome, but please click "Accept" if this has solved your issue.

[prashant singh]

Done but please let me know if I can visualize the convex hull on my volumetric image to validate the result.

[Matt J]

I'm not sure how you would do that. Possibly you could use a slice viewer from the File Exchange, e.g.,

https://www.mathworks.com/matlabcentral/fileexchange/27983-slicer

and then add lines to the 3D axes using the line() command (with HOLD ON active).