Search code examples
matlabcontainersvectorizationtriangulation

Crop a triangulation by indices


I have a large 3D triangulation (~2000 faces, ~1000 vertices), and I want to "crop" it, i.e. create a new tirangulation contains only the edges connecting a given list of vertex indices.

I work in MATLAB, and I tried to export the points and connectivity matrix to variables, and crop them manually:

% tri is the triangulation to be cropped, and inds is the vector of vertex 
% indices to remain in the final triangulation.

%% Create list of final vertices
p = tri.Points(inds, :);

%% Create new connectivity matrix
% start by cropping the connectivity matrix to contain only the final vertices
conn = [];
for j = 1:height(tri.ConnectivityList)
    if all( ismember(tri.ConnectivityList(j,:),inds) )
        conn = [conn; tri.ConnectivityList(j,:)];
    end
end

% then fix the indices to the new ones, since the index of each point changed
for i = 1:height(conn)
    for j = 1:width(conn)
        conn(i,j) = find(inds==conn(i,j),1);
    end
end

%% Create final triangulation
tri = triangulation(conn, p);

This solution works, but it is rather slow. Specifically, the cropping part using the ismember() function is a strong bottleneck.

After searching for a bit I found out somtimes people use containers.Map() to do similar jobs, but I couldn't understand if and how is this relevant to my case.

Any help is appreciated!

EDIT: here is an example conn matrix to work with:

1   2   3
4   2   5
2   4   3
6   5   7
8   9   10
11  12  7
13  10  14
15  16  17
14  16  15
18  19  17
20  21  22
23  24  25
24  26  25
27  21  24
20  28  21
27  24  23
29  13  12
10  23  14
14  23  25
30  29  31
31  11  1
30  31  32
33  34  8
34  35  9
11  7   2
7   36  6
17  19  37
36  17  37
38  39  40
41  22  39
41  42  22
42  20  22
43  44  28
45  20  42
40  39  35
39  22  35
46  33  8
47  48  49
47  50  51
49  52  50
53  54  18
55  53  18
55  56  53
57  58  53
53  58  59
25  56  55
25  26  56
24  60  61
62  63  57
56  64  57
65  66  67
67  68  65
68  69  62
68  62  64
61  64  56
65  64  61
69  70  62
19  18  48
48  54  49
71  72  49
71  73  72
71  58  73
59  71  54
59  58  71
73  74  75
76  77  74
63  78  79
37  47  80
2   7   5
81  82  6
81  6   80
83  81  80
5   6   82
46  84  85
29  84  13
85  84  29
84  46  13
46  8   13
8   10  13
8   34  9
7   12  36
36  12  15
12  14  15
13  14  12
36  15  17
14  25  16
17  16  18
16  25  18
9   86  10
10  86  27
10  27  23
9   22  86
86  21  27
22  21  86
24  21  60
21  28  60
29  12  11
87  31  1
87  32  31
31  29  11
30  85  29
88  39  89
40  35  34
11  2   1
9   35  22
38  89  39
41  39  88
20  45  43
41  88  42
43  28  20
47  37  48
51  80  47
47  49  50
49  72  52
18  25  55
53  59  54
56  26  61
60  28  66
61  60  65
60  66  65
90  57  63
64  62  57
64  65  68
26  24  61
63  91  92
62  70  78
63  62  78
63  92  90
53  56  57
58  57  90
37  19  48
18  54  48
49  54  71
73  58  74
74  58  76
58  90  76
75  74  77
91  63  79
80  51  83
44  66  28
6   36  80
80  36  37

pt is a points array of 92 points in 3D (92 by 3 matrix).


Solution

  • Your first for loop can be replaced by a vectorized call to ismember and use of logical indexing:

    conn = tri.ConnectivityList(all(ismember(tri.ConnectivityList, inds), 2), :);
    

    In this case, all is applied across the columns of the logical output from ismember. This results in a logical vector the same length as tri.ConnectivityList that can be used to extract just the valid rows.

    There are likely other ways to accomplish this. Notably, more general code might check if the number of points being "cropped" is more or less than half the number of original points. In that case, it might be more efficient to use logic to "crop-in" or "crop-out" points.