You are now following this question
- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.
Extracting points from a matrix tangential to a centreline
3 views (last 30 days)
Show older comments
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.
6 Comments
darova
on 13 Mar 2020
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
on 14 Mar 2020
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;
HB
on 15 Mar 2020
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
on 15 Mar 2020
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
Answers (2)
darova
on 16 Mar 2020
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
Result
See attached script
3 Comments
HB
on 17 Mar 2020
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
on 17 Mar 2020
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))
HB
on 17 Mar 2020
Edited: HB
on 4 May 2020
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
30 Comments
darova
on 17 Mar 2020
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?
darova
on 18 Mar 2020
this should work
ix = abs(V(3,:)) < 0.5; % data for crossection creating
dlmwrite([num2str(i) '.txt'],V(:,ix)','delimiter','\t')
HB
on 18 Mar 2020
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
on 18 Mar 2020
Put those lines inside for loop
Change height of crossection to grab more points
HB
on 16 Apr 2020
Great thanks, that works!! Is there any way to export the corresponding variable values which are in geom_variable(:,4)...
HB
on 18 Apr 2020
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!
HB
on 28 Apr 2020
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
on 28 Apr 2020
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
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
on 4 May 2020
Edited: darova
on 4 May 2020
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!
HB
on 4 May 2020
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
on 4 May 2020
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
on 7 May 2020
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
on 7 May 2020
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
on 7 May 2020
Edited: HB
on 7 May 2020
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
on 7 May 2020
- 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
on 8 May 2020
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?
HB
on 8 May 2020
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
on 8 May 2020
HB
on 8 Jul 2020
Edited: darova
on 8 Jul 2020
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
on 8 Jul 2020
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
See Also
Categories
Find more on Live Scripts and Functions in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!An Error Occurred
Unable to complete the action because of changes made to the page. Reload the page to see its updated state.
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list
How to Get Best Site Performance
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
Asia Pacific
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)