Extracting points from a matrix tangential to a centreline

Hello,

I have a tube with two branches (please see attached point cloud showing this model where the columns are x y z and a variable respectively). I also have a centreline that runs through the main branch of the model. Please also see attached centreline points where the columns are x y z respectively. I have also attached images showing the model and centreline. The red outline drawn on the model geometry highlights where the main branch is located.

At each centreline point I am trying to extract the points located tangentially along the model to that centreline point to a new matrix along with their corresponding variable value. I have attached two images that hopefully clarify better what I want to achieve. I am pretty new to Matlab so in depth explanations would be great.

Thanks in advance.

[darova]

The best idea i have: check each point of centerline (red) with data

ix = (xdata-xc(i)).^2 + (ydata-(yc(i)).^2 < R.^2;        % indices inside circle

[HB]

Thanks!

I tried the following script to check out your suggestion however I get the error "Array indices must be positive integers or logical values".

CL = load('Centerline.txt');
CLx = CL(:,1);
CLy = CL(:,2);
CLz = CL(:,3);

geom_variable=load('geom_variable.txt');
geomx = geom_variable(:,1);
geomy = geom_variable(:,2);
geomz = geom_variable(:,3);



ix = (geomx-CLx(i)).^2 + (geomy-(CLy(i)).^2 < R.^2;

[darova]

What are these?

[HB]

Yes! Sorry. i = 46 and R= 1.5. ix matrix produced has the same number of rows as geomx, geomy and geomz. How do I then extract values that intersect with the radius ?

[darova]

Maybe better be to use pdist2. IT calculates every possible combination of distances

You have two set of data (46 and 113502). So it creates matrix of size 46x11502

Try this

D = pdist2([CLx CLy CLz],[geomx geomy geomz]);      % matrix of size 46x113502
D1 = D < 2;                                         % distances smaller than "2" (matrix "0" and "1")
[~,ix] = find(sum(D1)>0);                           % if columns has one "1" - find it

plot3(CLx,CLy,CLz,'r')                              % centerline
hold on
plot3(geomx,geomy,geomz,'.b')                       % all data
plot3(geomx(ix),geomy(ix),geomz(ix),'.r')           % data close to centerline
hold off

[HB]

Great!

So D1 now contains all the points excluding the side branches and points outside the centreline (see attached). What would be the best way to now get the cross section points at every centreline point as seen in desired_output?

Thanks!

Answer 1

[darova]

Here is what i did

I moved data to origin and rotated at some points of centerline

Once centerline is vertical i extract data -1<z<1 (for example)

I created surface at some points of centerline

countour was used to get crossection

I used Rodrigues rotation matrix/formula to rotate data

Result

See attached script

[HB]

Thanks, that is such a great way of doing it, very helpful!!

What do you think would be the best way of exporting the x,y,z values of the cross-sections to a text file and their corresponding variable value which are located in geom_variable(1:1:end,4)?

[darova]

To write data

DATA = [];
for i = 1...
    % code
    V(:,1) = [];
    DATA = [DATA V'];
end
dlmwrite('myFile.txt',DATA,'delimiter','\t')

• and their corresponding variable value which are located in geom_variable(1:1:end,4)?

Im afraid crossection line doesn't have corresponding values of geom_variable.

Closest points/values (red) can be found. But it's always different number

Change these lines in a loop

ix = abs(V(3,:)) < 0.1;                         % data for crossection creating
length(find(ix))

[darova]

Answer 2

[HB]

Ah I see, what I need is cross sections that consists of points from my geometry and the corresponding variable which is shear stress. The data is actually an artery model from a CFD.

What does the values length(find(ix)) values refer to? And how would I export the points shown in your latest plot?

Thanks

[darova]

This line means the number of points

length(find(ix))

Once the number of points is different (very simple scheme)

cross1 cross2 cross3 ...
    1       2      1
    5       4      2
            3      2
            5

How do you want to have points? Separate txt?

[HB]

Yes a text file would be the best. Thanks

[darova]

this should work

ix = abs(V(3,:)) < 0.5;                         % data for crossection creating
dlmwrite([num2str(i) '.txt'],V(:,ix)','delimiter','\t')

[HB]

It seems to have written one text file 45.txt containing one of the cross sections?

Also this cross section has 14 points, is there an easy way of extracting more points (as many as possible)?

Thanks

[darova]

Put those lines inside for loop

Change height of crossection to grab more points

[HB]

Hello,

I have come back to this and I have put those lines you specified inside the loop:

However when I view the ouput files it seems that the data extracted is not at every centreline point (as seen below)..any suggestions why? Thanks!

[darova]

You should have a few files (separate crossection). Do you have them?

[HB]

Yep, I do have 46 separate txt files (which is the number of centreline points). What is shown in the above image is all those files loaded together..as shown it strangely doesn't seem to be the data at every centreline point.

[darova]

Sorry. I made a mistake

Try this script

[HB]

Great thanks, that works!! Is there any way to export the corresponding variable values which are in geom_variable(:,4)...

[darova]

There are no 'corresponding' values

They can be interpolated

See script

[HB]

Thanks! Yes by corresponding variable, I meant that for every xyz position I have a variable value corresponding to that position.

Apologies for another question: I triedtthe script on another case (see attached files), but for some reason I am recieving the following error

"Error using linspace (line 22)

Inputs must be scalars.

Error in crossection (line 37)

tt = linspace(min(t)+0.02,max(t)-0.02,3000);"

Any suggestions why? Thanks!

[darova]

The data looks strange. Maybe units?

[darova]

I scaled (100x) geometry. Just a few similar points

[HB]

Ah yes, you're correct, it is a unit issue!

Quick question - From the script I don't wholly understand how the the variable values are interpolated for each point? Do you mind clarifying? Thanks!

[darova]

Imagine this is your data (x,y,z and c - some value - color)

Your data is like cylinder. I usually use interp2 or griddata for interpolation. Your data should like this (for each X Y only one Z)

So to make your data more 'flat' i unwided it (don't know if it's correct word). Look - >

[t0,r0] = cart2pol(x0,y0);          % convert data from cartesian to polar

I used griddata to interpolate your data. Then converted back to cartesian

R = griddata(t0,z0,r0,T,Z);     % calculate radius for grid points
C = griddata(t0,z0,c0,T,Z);     % calculate color for grid points
[X,Y] = pol2cart(T,R);          % convert back to cartesian

[HB]

Thank you, that's a great explanation!

I was actually trying the script on a some other cases, and for some reason on a few of them points of the cross section line are missing as seen in image below. I don't know if there is something that can be modified in the script to fix this? I have attached the files FYI. Thanks!

[darova]

Is it acceptable?

[HB]

Yes thanks, now there are no gaps! however for the variable value there are NaN values for some xyz positions..can these values be obtained/interpolated?

[darova]

My bad. griddata has some problem at the boundaries. scatteredInterpolant uses extrapolation

I changed these lines

F     = scatteredInterpolant(t(:),5*gz1(:),gc(ix));
gc1   = F(ct,ct*0);

[HB]

Thanks, thats great!

I actually hope you don't mind me asking you for your opinion on something. I have second approach to achieving the above. Please see attached script. However it doesn't quite work for cases where there is a branch coming of the primary branch such as the cases I've shown you (because it also tries to extract the secondary branch). However it works fine when there is only the primary branch. It would be very helpful if you could maybe see if you can spot anything in the script that could be modified to make it work. If you don't think it can be, what are the issues with it in your opinion? Thanks a lot!

[darova]

Sure, look on this picture

As i understood correctly: the script calculates red distance and if it's small enough the point is taken

You can calculate green distance too and compare it

d2 = pdist2(geomp(igeomp,1:3), CLp(iCLp,1:3));  % green distance
if dtmp<dist && d2 < 5                          % compare red and green distance

Is it that you wanted?

[HB]

Well the script extracts cross sections of the geometry at each CL point. So I think as well as the red distance it also calculated distances to the wall somehow.

If you were to run the attached case (no branches coming off primary) it does so fine. But not for the case I attached in the previous message, which I am hoping I can get to work with the script or figure out why it's not the best approach for cases with branches.

Your suggestions/assistance would be great! Thanks

[darova]

• So I think as well as the red distance it also calculated distances to the wall somehow.

It doesn't that's why second branch is taken. I changed your script a bit

[HB]

Thank you, that is so helpful!

I was wondering in your script how do you actually round off/close the gaps in the cross sections?

[darova]

Don't see any gap

[HB]

Ah sorry! I was referring to the script that you gave from my question on the 4th May. You edited the script to close gaps in the cross sections, I was just wondering how this was done? Also the script 'rounds off' the geometry where there is the branch to create the cross section, how is it done exactly? Is it the same way as the gaps were closed? Thanks!

[darova]

I used griddata first to interpolate values to get crossection line. See this post: LINK

But griddata has some problems at boundaries

So i replaced griddata with scatteredInterpolant, it uses extrapolation: LINK

[HB]

Hello,

Hope you’re well!

I’ve been running cases with the latest script (see attached crossection_latest.m). It seems able to extract some cross sections but not others. For instance if I change line 26 from “i = 1:5:size(CL,1)-1” to “i = 1:size(CL,1)-1” therefore to include all the cross sections I get the error:

Error using linspace (line 22)
Inputs must be scalars.
Error in crossection_latest (line 43)
    zz    = linspace(min(gz1),max(gz1),10);

Any suggestions on how this can be overcome this (I have attached the files for a case)? It happens for a few of my cases, many thanks in advance!

[darova]

Hello!

I think the problem you described is cause by radius you choosed. You didn't found all the points of an arteria

ix = sum(D<2) > 0;                 % if columns has one "1" - find it

You arteria is thicker than you think. Try to increase the radius

ix = sum(D<2.5) > 0;                 % if columns has one "1" - find it