diff options
| author | Willem Jan Palenstijn <Willem.Jan.Palenstijn@cwi.nl> | 2018-12-12 16:14:00 +0100 |
|---|---|---|
| committer | Willem Jan Palenstijn <Willem.Jan.Palenstijn@cwi.nl> | 2018-12-12 16:14:00 +0100 |
| commit | c549714416a8b8500820d0ceb2cb1ebc6ff38e55 (patch) | |
| tree | 739273353f33f31fb4210937174218dc36987136 /matlab | |
| parent | c712d102e7a5a67468bc46db14f0510b8f6f4381 (diff) | |
| parent | 3276485c96636cd38248908ff3575282654ff335 (diff) | |
| download | astra-c549714416a8b8500820d0ceb2cb1ebc6ff38e55.tar.gz astra-c549714416a8b8500820d0ceb2cb1ebc6ff38e55.tar.bz2 astra-c549714416a8b8500820d0ceb2cb1ebc6ff38e55.tar.xz astra-c549714416a8b8500820d0ceb2cb1ebc6ff38e55.zip | |
Merge branch 'geom_visualizer'
This adds a matlab geometry visualizer, and a sample showing how to use it.
Diffstat (limited to 'matlab')
35 files changed, 2058 insertions, 0 deletions
diff --git a/matlab/algorithms/plot_geom/+draw/draw_cad_phantom.m b/matlab/algorithms/plot_geom/+draw/draw_cad_phantom.m new file mode 100644 index 0000000..55dedc4 --- /dev/null +++ b/matlab/algorithms/plot_geom/+draw/draw_cad_phantom.m @@ -0,0 +1,35 @@ +function [] = draw_cad_phantom(filename, magn) +%% draw_cad_phantom.m +% brief render an stl model into a 3d axis object +% param vol_geom volume geometry describing the phantom +% param h_ax handle to axis to plot into +% param magn magnification multiplier of the phantom. default = 1 +% +% date 02.07.2018 +% author Alice Presenti +% imec VisionLab +% University of Antwerp +% Modified by Tim Elberfeld +%% + h_ax = gca; + if nargin == 1 + magn = 1; + end + + [v,f,~,~] = stlTools.stlRead(filename); + m = mean(v); % to center the CAD model! + for i=1:3 + v(:,i) = (v(:,i)- m(i)) .* magn; + end + object.vertices = v; + object.faces = f; + patch(h_ax, object,'FaceColor', [0.8 0.8 1.0], ... + 'EdgeColor', 'none', ... + 'FaceLighting', 'gouraud', ... + 'AmbientStrength', 0.15); + % Add a camera light, and tone down the specular highlighting + camlight('headlight'); + material('dull'); + hold off; + +end
\ No newline at end of file diff --git a/matlab/algorithms/plot_geom/+draw/draw_cone_geom.m b/matlab/algorithms/plot_geom/+draw/draw_cone_geom.m new file mode 100644 index 0000000..ae98294 --- /dev/null +++ b/matlab/algorithms/plot_geom/+draw/draw_cone_geom.m @@ -0,0 +1,73 @@ +function [] = draw_cone_geom(h_ax, geom, options) +%% draw_cone_geom.m +% draw an astra cone beam projection geometry +% +% param h_ax handle to axis to draw into +% param geom the geometry to draw +% param options struct containing the options for this function +% SourceMarker marker to use to mark the source location +% SourceMarkerColor Color of the marker for the source +% DetectorLineColor Color of the lines that draw the detector +% DetectorLineWidth Width of the lines that draw the detector +% OpticalAxisColor color of the lines representing the optical axis +% +% date 21.06.2018 +% author Tim Elberfeld +% imec VisionLab +% University of Antwerp +%% + hold on; + % draw origin + h_origin = scatter3(h_ax, 0,0,0, '+k'); + h_origin.SizeData = 120; + h_origin.LineWidth = 2; + + % draw lines between source, origin and detector + line(h_ax, [0, 0], [0, geom.DistanceOriginDetector], [0, 0], 'Color', options.OpticalAxisColor); + line(h_ax, [0, 0], [0, -geom.DistanceOriginSource], [0, 0], 'Color', options.OpticalAxisColor); + + % draw source + scatter3(h_ax, 0, -geom.DistanceOriginSource, 0, options.SourceMarker,... + options.SourceMarkerColor); + + detector = draw_detector(h_ax, geom, options); + + % connect source to detector edges + for idx = 1:4 + line(h_ax,[0, detector.vertices(1, idx)], ... + [-geom.DistanceOriginSource, geom.DistanceOriginDetector],... + [0, detector.vertices(3, idx)],... + 'Color', 'k', 'LineStyle', ':'); + end + + % draw rotation axis + line(h_ax, [0,0],[0,0], 0.6*[-detector.height, detector.height],... + 'LineWidth', options.DetectorLineWidth, 'Color',... + 'k', 'LineStyle', '--'); + + perc = 0.05; + text(h_ax, perc*detector.width, perc*detector.width,... + 0.8*detector.height, 'rotation axis'); + text(h_ax, detector.width*perc, 0, perc*detector.height, 'origin'); + text(h_ax, detector.width*perc, -geom.DistanceOriginSource,... + perc*detector.height, 'x-ray source'); + text(h_ax, 0, geom.DistanceOriginDetector, 0, 'detector'); + hold off; + + function [detector] = draw_detector(h_ax, geom, options) + detector = struct; + detector.height = geom.DetectorRowCount * geom.DetectorSpacingY; + detector.width = geom.DetectorColCount * geom.DetectorSpacingX; + + vertices = zeros(3, 5); + vertices(1, :) = 0.5*[-detector.width, -detector.width,... + detector.width, detector.width, -detector.width]; + vertices(2, :) = repmat([geom.DistanceOriginDetector], 5, 1); + vertices(3, :) = 0.5*[detector.height, -detector.height,... + -detector.height, detector.height, detector.height]; + + detector.vertices = vertices; + plot3(h_ax, detector.vertices(1, :), detector.vertices(2, :),... + detector.vertices(3, :), options.DetectorLineColor); + end +end diff --git a/matlab/algorithms/plot_geom/+draw/draw_cone_vec_geom.m b/matlab/algorithms/plot_geom/+draw/draw_cone_vec_geom.m new file mode 100644 index 0000000..dac7410 --- /dev/null +++ b/matlab/algorithms/plot_geom/+draw/draw_cone_vec_geom.m @@ -0,0 +1,108 @@ +function [] = draw_cone_vec_geom(h_ax, geom, options) +%% draw_cone_vec_geom.m +% draw an astra cone beam vectorized projection geometry +% param h_ax handle to axis to draw into +% param geom the geometry to draw +% param options struct holding options for drawing +% RotationAxis if specified, will change the rotation axis +% from z-axis to provided axis. +% Must be 3-vector. Default value [0, 0, 1] +% VectorIdx the index of the angle to draw +% Color global color for all lines and markers. +% is overwritten by specialized color options +% SourceMarker Marker to use to mark the source +% SourceMarkerColor Color of the source marker +% DetectorMarker Marker to use to mark the source +% DetectorMarkerColor Color of the source marker +% DetectorLineColor Color of the outline of the detector +% OpticalAxisColor color for drawing the optical axis +% +% date 21.06.2018 +% author Tim Elberfeld, Van Nguyen +% imec VisionLab +% University of Antwerp +% +% - last update 07.11.2018 +%% + vectors = geom.Vectors; + + % source + xray_source = vectors(:, 1:3); + % center of detector + detector_center = vectors(:, 4:6); + + % draw the points and connect with lines + hold on; + num_angles = size(vectors, 1); + s_source = scatter3(h_ax, xray_source(:, 1), xray_source(:, 2),... + xray_source(:, 3)); + s_source.Marker = options.SourceMarker; + s_source.MarkerEdgeColor = options.SourceMarkerColor; + + s_det = scatter3(h_ax, detector_center(:, 1),... + detector_center(:, 2), detector_center(:, 3)); + s_det.MarkerEdgeColor = options.DetectorMarkerColor; + s_det.Marker = options.DetectorMarker; + + detector = struct; + detector.u = vectors(options.VectorIdx, 7:9); + detector.v = vectors(options.VectorIdx, 10:12); + detector.height = geom.DetectorColCount; + detector.width = geom.DetectorRowCount; + detector.origin = detector_center(options.VectorIdx, :); + + vertices = draw.draw_detector_vec(h_ax, detector, options); + connect_source_detector(h_ax, vertices, detector_center, ... + xray_source, options); + + % rotation axis will be roughly as long as the source detector distance + distances = eucl_dist3d(detector_center, xray_source); + mean_sdd = mean(distances(:)); % mean source detector distance + draw_rotation_axis(h_ax, mean_sdd, options); + + text(h_ax, xray_source(options.VectorIdx, 1),... + xray_source(options.VectorIdx, 2),... + xray_source(options.VectorIdx, 3), 'x-ray source'); + text(h_ax, detector_center(options.VectorIdx, 1),... + detector_center(options.VectorIdx, 2),... + detector_center(options.VectorIdx, 3), 'detector'); + hold off; + + function [] = connect_source_detector(h_ax, vertices,... + detector_center, xray_source, options) + % connect source to detector origin + idx = options.VectorIdx; + line(h_ax, [detector_center(idx, 1), xray_source(idx, 1)],... + [detector_center(idx, 2), xray_source(idx, 2)],... + [detector_center(idx, 3), xray_source(idx, 3)],... + 'Color', options.OpticalAxisColor, 'LineStyle', '--'); + + % connect source to detector edges + for kk = 1:4 + line(h_ax,[xray_source(idx, 1), vertices(1, kk)], ... + [xray_source(idx, 2), vertices(2, kk)],... + [xray_source(idx, 3), vertices(3, kk)],... + 'Color', 'k', 'LineStyle', ':'); + end + end + + function [] = draw_rotation_axis(h_ax, scaling, options) + % draw rotation axis + rot_axis = options.RotationAxis; + if(~isnan(rot_axis(1))) + rot_axis = options.RotationAxis + options.RotationAxisOffset; + origin = options.RotationAxisOffset; + % origin of the geometry is assumed to be [0, 0, 0] always! + line(h_ax, [origin(1), (scaling/2)*rot_axis(1)],... + [origin(2), (scaling/2)*rot_axis(2)],... + [origin(3), (scaling/2)*rot_axis(3)],... + 'Color', options.OpticalAxisColor,... + 'LineStyle', '-.'); + line(h_ax, [origin(1), -(scaling/2)*rot_axis(1)],... + [origin(2), -(scaling/2)*rot_axis(2)],... + [origin(3), -(scaling/2)*rot_axis(3)],... + 'Color', options.OpticalAxisColor,... + 'LineStyle', '-.'); + end + end +end diff --git a/matlab/algorithms/plot_geom/+draw/draw_detector_vec.m b/matlab/algorithms/plot_geom/+draw/draw_detector_vec.m new file mode 100644 index 0000000..5c293d4 --- /dev/null +++ b/matlab/algorithms/plot_geom/+draw/draw_detector_vec.m @@ -0,0 +1,42 @@ +function [ vertices] = draw_detector_vec( h_ax, detector, options) +%% draw_detector_vec.m +% draw a detector for a vector geometry +% param h_ax handle to axis to draw into +% param detector the struct specifying the detector +% .origin 3d coordinates of the detector origin +% .u vector pointing from detector origin to +% pixel (1,0) +% .v vector pointing from detector origin to +% pixel (0,1) +% .width width of the detector (number of px in u +% direction) +% .height height of the detector (number of px in v +% direction) +% param options struct with options +% .Color Color for the line work +% .DetectorLineColor Color of the detector rectangle outline +% return The vertices of the detector rectangle that +% were plotted +% +% date 09.07.2018 +% author Tim Elberfeld +% imec VisionLab +% University of Antwerp +%% + % draw the detector rectangle + vertices = zeros(3, 5); + vertices(:, 1) = detector.origin - detector.u * detector.width / 2 + ... + detector.v * detector.height / 2; + vertices(:, 2) = detector.origin + detector.u * detector.width / 2 + ... + detector.v * detector.height / 2; + vertices(:, 3) = detector.origin + detector.u * detector.width / 2 - ... + detector.v * detector.height / 2; + vertices(:, 4) = detector.origin - detector.u * detector.width / 2 - ... + detector.v * detector.height / 2; + vertices(:, 5) = vertices(:, 1); + + detector.vertices = vertices; + plot3(h_ax, detector.vertices(1, :), detector.vertices(2, :),... + detector.vertices(3, :), options.DetectorLineColor); + +end diff --git a/matlab/algorithms/plot_geom/+draw/draw_fanflat_geom.m b/matlab/algorithms/plot_geom/+draw/draw_fanflat_geom.m new file mode 100644 index 0000000..a5f28ea --- /dev/null +++ b/matlab/algorithms/plot_geom/+draw/draw_fanflat_geom.m @@ -0,0 +1,25 @@ +function [] = draw_fanflat_geom( h_ax, geom, options) +%% draw_fanflat_geom.m +% draw an astra cone beam projection geometry +% param h_ax handle to axis to draw into +% param geom the geometry to draw +% param options struct holding options for drawing +% SourceMarker Marker to use to mark the source +% SourceMarkerColor Color of the source marker +% DetectorMarker Marker to use to mark the source +% DetectorMarkerColor Color of the source marker +% DetectorLineColor Color of the outline of the detector +% OpticalAxisColor color for drawing the optical axis +% date 28.06.2018 +% author Tim Elberfeld +% imec VisionLab +% University of Antwerp +%% + % convert to faux cone geometry so we don't have to write more code :)! + cone_geom = astra_create_proj_geom('cone', geom.DetectorWidth,... + geom.DetectorWidth, 1, geom.DetectorCount, geom.ProjectionAngles,... + geom.DistanceOriginSource, geom.DistanceOriginDetector); + + draw.draw_cone_geom(h_ax, cone_geom, options); +end + diff --git a/matlab/algorithms/plot_geom/+draw/draw_fanflat_vec_geom.m b/matlab/algorithms/plot_geom/+draw/draw_fanflat_vec_geom.m new file mode 100644 index 0000000..761d96a --- /dev/null +++ b/matlab/algorithms/plot_geom/+draw/draw_fanflat_vec_geom.m @@ -0,0 +1,102 @@ +function [ output_args ] = draw_fanflat_vec_geom( h_ax, geom, options) +%% draw_fanflat_vec_geom.m +% draw an astra cone beam projection geometry +% param h_ax handle to axis to draw into +% param geom the geometry to draw +% param options struct holding options for drawing +% VectorIdx the index of the angle to draw +% SourceMarker Marker to use to mark the source +% SourceMarkerColor Color of the source marker +% DetectorMarker Marker to use to mark the source +% DetectorMarkerColor Color of the source marker +% DetectorLineColor Color of the outline of the detector +% OpticalAxisColor color for drawing the optical axis +% +% date 28.06.2018 +% author Tim Elberfeld +% imec VisionLab +% University of Antwerp +% +% - last update 09.07.2018 +%% + vectors = geom.Vectors; + + % source + num_angles = size(vectors, 1); + xray_source = [vectors(:, 1:2), zeros(num_angles, 1)]; + % center of detector + detector_center = [vectors(:, 3:4), zeros(num_angles, 1)]; + + % draw the points and connect with lines + hold on; + scatter3(h_ax, xray_source(:, 1), xray_source(:, 2),... + xray_source(:, 3), options.SourceMarker,... + options.SourceMarkerColor); + scatter3(h_ax, detector_center(:, 1),... + detector_center(:, 2), detector_center(:, 3),... + options.DetectorMarker, options.DetectorMarkerColor); + + detector = struct; + detector.u = [vectors(options.VectorIdx, 5:6), 0]; + detector.v = fliplr(detector.u); + detector.height = 1; + detector.width = geom.DetectorCount; + detector.origin = detector_center(options.VectorIdx, :); + + vertices = draw.draw_detector_vec(h_ax, detector, options); + connect_source_detector(h_ax, vertices, detector_center,... + xray_source, options); + + % rotation axis will be roughly as long as the source detector distance + distances = eucl_dist3d(detector_center, xray_source); + mean_sdd = mean(distances(:)); % mean source detector distance + draw_rotation_axis(h_ax, mean_sdd, options); + + text(h_ax, xray_source(options.VectorIdx, 1),... + xray_source(options.VectorIdx, 2),... + xray_source(options.VectorIdx, 3), 'x-ray source'); + text(h_ax, detector_center(options.VectorIdx, 1),... + detector_center(options.VectorIdx, 2),... + detector_center(options.VectorIdx, 3), 'detector'); + + hold off; + + function [] = connect_source_detector(h_ax, vertices,... + detector_center, xray_source, options) + idx = options.VectorIdx; + % connect source to detector origin + line(h_ax, [detector_center(idx, 1), xray_source(idx, 1)],... + [detector_center(idx, 2), xray_source(idx, 2)],... + [detector_center(idx, 3), xray_source(idx, 3)],... + 'Color', options.OpticalAxisColor, 'LineStyle', '--'); + + % connect source to detector edges + for kk = 1:4 + line(h_ax,[xray_source(idx, 1), vertices(1, kk)], ... + [xray_source(idx, 2), vertices(2, kk)],... + [xray_source(idx, 3), vertices(3, kk)],... + 'Color', 'k', 'LineStyle', ':'); + end + end + + function [] = draw_rotation_axis(h_ax, scaling, options) + % draw rotation axis + rot_axis = options.RotationAxis; + if(~isnan(rot_axis(1))) + rot_axis = options.RotationAxis + options.RotationAxisOffset; + origin = options.RotationAxisOffset; + % origin of the geometry is assumed to be [0, 0, 0] always! + line(h_ax, [origin(1), (scaling/2)*rot_axis(1)],... + [origin(2), (scaling/2)*rot_axis(2)],... + [origin(3), (scaling/2)*rot_axis(3)],... + 'Color', options.OpticalAxisColor,... + 'LineStyle', '-.'); + line(h_ax, [origin(1), -(scaling/2)*rot_axis(1)],... + [origin(2), -(scaling/2)*rot_axis(2)],... + [origin(3), -(scaling/2)*rot_axis(3)],... + 'Color', options.OpticalAxisColor,... + 'LineStyle', '-.'); + end + end +end + diff --git a/matlab/algorithms/plot_geom/+draw/draw_parallel3d_geom.m b/matlab/algorithms/plot_geom/+draw/draw_parallel3d_geom.m new file mode 100644 index 0000000..45a8d60 --- /dev/null +++ b/matlab/algorithms/plot_geom/+draw/draw_parallel3d_geom.m @@ -0,0 +1,88 @@ +function [ ] = draw_parallel3d_geom( h_ax, geom, options) +%% draw_parallel3d_geom.m +% draw an astra parallel3d projection geometry +% param h_ax handle to axis to draw into +% param geom the geometry to draw +% param options options struct with additional settings +% .SourceDistance distance of source to origin +% .SourceMarker marker for the source locations +% .SourceMarkerColor color specifier for the source marker +% .DetectorMarker marker for the detector locations. +% Default = '.' +% .DetectorMarkerColor color specifier for the detector marker. +% Default = 'k' +% .DetectorLineColor color for the lines drawing the detector +% outline +% .OpticalAxisColor Color of the line representing the optical +% axis. +% date 22.06.2018 +% author Tim Elberfeld +% imec VisionLab +% University of Antwerp +% +% - last update 09.07.2018 +%% + dist_origin_detector = options.SourceDistance; + dist_origin_source = options.SourceDistance; + hold on; + + % draw source + scatter3(h_ax, 0, -dist_origin_source, 0, options.SourceMarker,... + options.SourceMarkerColor); + + % draw detector + detector = draw_detector(h_ax, geom, dist_origin_detector, options); + + % draw origin + h_origin = scatter3(h_ax, 0,0,0, '+k'); + h_origin.SizeData = 120; + h_origin.LineWidth = 2; + + % draw lines between source, origin and detector + line(h_ax, [0, 0], [0, dist_origin_detector], [0, 0],... + 'Color', options.OpticalAxisColor); + line(h_ax, [0, 0], [0, -dist_origin_source], [0, 0],... + 'Color', options.OpticalAxisColor); + + % connect source to detector edges + for idx = 1:4 + line(h_ax,[detector.vertices(1, idx),... + detector.vertices(1, idx)],... + [-dist_origin_source, dist_origin_detector],... + [detector.vertices(3, idx), detector.vertices(3, idx)],... + 'Color', 'k', 'LineStyle', ':'); + end + + % draw rotation axis + line(h_ax, [0,0],[0,0], 0.6*[-detector.height, detector.height],... + 'LineWidth', 2, 'Color', 'k', 'LineStyle', '--'); + + perc = 0.05; + text(h_ax, perc*detector.width, perc*detector.width,... + 0.8*detector.height, 'rotation axis'); + text(h_ax, detector.width*perc, 0, perc*detector.height, 'origin'); + text(h_ax, detector.width*perc, -dist_origin_source,... + perc*detector.height, 'x-ray source'); + text(h_ax, 0, dist_origin_detector, 0, 'detector'); + hold off; + + function [detector] = draw_detector(h_ax, geom,... + dist_origin_detector, options) + detector = struct; + detector.height = geom.DetectorRowCount * geom.DetectorSpacingY; + detector.width = geom.DetectorColCount * geom.DetectorSpacingX; + + vertices = zeros(3, 5); + vertices(1, :) = 0.5*[-detector.width, -detector.width,... + detector.width, detector.width, -detector.width]; + vertices(2, :) = repmat(dist_origin_detector, 5, 1); + vertices(3, :) = 0.5*[detector.height, -detector.height,... + -detector.height, detector.height, detector.height]; + + detector.vertices = vertices; + plot3(h_ax, detector.vertices(1, :), detector.vertices(2, :),... + detector.vertices(3, :), options.DetectorLineColor); + end + +end + diff --git a/matlab/algorithms/plot_geom/+draw/draw_parallel3d_vec_geom.m b/matlab/algorithms/plot_geom/+draw/draw_parallel3d_vec_geom.m new file mode 100644 index 0000000..e59fa91 --- /dev/null +++ b/matlab/algorithms/plot_geom/+draw/draw_parallel3d_vec_geom.m @@ -0,0 +1,100 @@ +function [ ] = draw_parallel3d_vec_geom( h_ax, geom, options) +%% draw_parallel3d_vec_geom.m +% draw an astra parallel3d projection geometry +% param h_ax handle to axis to draw into +% param geom the geometry to draw +% param options struct containing the options for this function +% idx index of the vector to draw in more detail +% SourceDistance distance of the source to the origin +% +% date 22.06.2018 +% author Tim Elberfeld +% imec VisionLab +% University of Antwerp +%% + vectors = geom.Vectors; + + % source + xray_source = vectors(:, 1:3)*options.SourceDistance; + + % center of detector + detector_center = vectors(:, 4:6); + + % draw the points and connect with lines + hold on; + num_angles = size(vectors, 1); + scatter3(h_ax, xray_source(:, 1), xray_source(:, 2),... + xray_source(:, 3), options.SourceMarkerColor,... + options.SourceMarker); + scatter3(h_ax, detector_center(:, 1), detector_center(:, 2),... + detector_center(:, 3), 'k.'); + + detector = struct; + detector.u = vectors(options.VectorIdx, 7:9); + detector.v = vectors(options.VectorIdx, 10:12); + detector.height = geom.DetectorColCount; + detector.width = geom.DetectorRowCount; + detector.origin = detector_center(options.VectorIdx, :); + detector.vertices = draw.draw_detector_vec(h_ax, detector, options); + + connect_source_detector(h_ax, detector, detector_center, xray_source,... + options); + + % rotation axis will be roughly as long as the source detector distance + distances = eucl_dist3d(detector_center, xray_source); + mean_sdd = mean(distances(:)); % mean source detector distance + draw_rotation_axis(h_ax, mean_sdd, options); + + text(h_ax, xray_source(options.VectorIdx, 1),... + xray_source(options.VectorIdx, 2), ... + xray_source(options.VectorIdx, 3), 'x-ray source'); + text(h_ax, detector_center(options.VectorIdx, 1),... + detector_center(options.VectorIdx, 2),... + detector_center(options.VectorIdx, 3), 'detector'); + hold off; + + function [] = connect_source_detector(h_ax, detector,... + detector_center, xray_source, options) + % connect source to detector origin + idx = options.VectorIdx; + line(h_ax, [detector_center(idx, 1), xray_source(idx, 1)],... + [detector_center(idx, 2), xray_source(idx, 2)],... + [detector_center(idx, 3), xray_source(idx, 3)],... + 'Color', options.OpticalAxisColor, 'LineStyle', '--'); + + % compute normal of detector plane + n = null([detector.u; detector.v]); + + % connect source to detector edges + for kk = 1:4 + a = detector.vertices(1, kk) - n(1)*xray_source(idx, 1); + b = detector.vertices(2, kk) - n(2)*xray_source(idx, 2); + c = detector.vertices(3, kk) - n(3)*xray_source(idx, 3); + line(h_ax,[a, detector.vertices(1, kk)], ... + [b, detector.vertices(2, kk)],... + [c, detector.vertices(3, kk)],... + 'Color', 'k',... + 'LineStyle', ':'); + end + end + + function [] = draw_rotation_axis(h_ax, scaling, options) + % draw rotation axis + rot_axis = options.RotationAxis; + if(~isnan(rot_axis(1))) + rot_axis = options.RotationAxis + options.RotationAxisOffset; + origin = options.RotationAxisOffset; + % origin of the geometry is assumed to be [0, 0, 0] always! + line(h_ax, [origin(1), (scaling/2)*rot_axis(1)],... + [origin(2), (scaling/2)*rot_axis(2)],... + [origin(3), (scaling/2)*rot_axis(3)],... + 'Color', options.OpticalAxisColor,... + 'LineStyle', '-.'); + line(h_ax, [origin(1), -(scaling/2)*rot_axis(1)],... + [origin(2), -(scaling/2)*rot_axis(2)],... + [origin(3), -(scaling/2)*rot_axis(3)],... + 'Color', options.OpticalAxisColor,... + 'LineStyle', '-.'); + end + end +end diff --git a/matlab/algorithms/plot_geom/+draw/draw_proj_geom.m b/matlab/algorithms/plot_geom/+draw/draw_proj_geom.m new file mode 100644 index 0000000..9691448 --- /dev/null +++ b/matlab/algorithms/plot_geom/+draw/draw_proj_geom.m @@ -0,0 +1,138 @@ +function [] = draw_proj_geom(geom, varargin) +%% draw_proj_geom.m +% brief rendering function for astra geometries. +% param geom the geometry to plot. If geometry type +% is not supported, throws error +% ------------------------------ +% optional parameters that can be provided as string value pairs: +% +% param RotationAxis if specified, will change the drawn +% rotation axis to provided axis. +% Must be 3-vector. Default value is +% [NaN, NaN, NaN], (meaning do not draw). +% param RotationAxisOffset if specified, will translate the drawn +% rotation axis by the provided vector. +% Default = [0, 0, 0] +% param VectorIdx index of the vector to visualize if geom +% is a vector geometry type. Default = 1 +% param Color Color for all markers and lines if not +% otherwise specified +% param DetectorMarker marker for the detector locations. +% Default = '.' +% param DetectorMarkerColor color specifier for the detector marker. +% Default = 'k' +% param DetectorLineColor color for the lines drawing the detector +% outline +% param DetectorLineWidth line width of detector rectangle +% param SourceMarker marker for the source locations +% param SourceMarkerColor color specifier for the source marker +% param SourceDistance (only for parallel3d and parallel3d_vec) +% distance of source to origin +% param OpticalAxisColor Color for drawing the optical axis +% +% date 20.06.2018 +% author Tim Elberfeld +% imec VisionLab +% University of Antwerp +% +% - last update 07.11.2018 +%% + h_ax = gca; + options = parseoptions(varargin); + + switch geom.type + case 'parallel3d' + disp('type: parallel3d') + disp(['detector spacing: [' num2str(geom.DetectorSpacingX), ', '... + num2str(geom.DetectorSpacingY) ']']); + disp(['detector px: [' num2str(geom.DetectorRowCount), ', ' ... + num2str(geom.DetectorColCount) ']']); + disp(['angle lo: ' num2str(geom.ProjectionAngles(1))]); + disp(['angle hi: ' num2str(geom.ProjectionAngles(end))]); + disp(['# angles: ' num2str(numel(geom.ProjectionAngles))]); + disp('DistanceOriginDetector inf'); + disp('DistanceOriginSource inf'); + + draw.draw_parallel3d_geom(h_ax, geom, options); + case 'parallel3d_vec' + disp('type: parallel3d_vec') + disp(['detector px: [' num2str(geom.DetectorRowCount), ', '... + num2str(geom.DetectorColCount) ']']); + disp(['# angles: ' num2str(size(geom.Vectors, 1))]); + + draw.draw_parallel3d_vec_geom(h_ax, geom, options); + case 'cone' + disp('type: cone'); + disp(['detector spacing: [' num2str(geom.DetectorSpacingX), ', '... + num2str(geom.DetectorSpacingY) ']']); + disp(['detector px: [' num2str(geom.DetectorRowCount), ', ' ... + num2str(geom.DetectorColCount) ']']); + disp(['angle lo: ' num2str(geom.ProjectionAngles(1))]); + disp(['angle hi: ' num2str(geom.ProjectionAngles(end))]); + disp(['# angles: ' num2str(numel(geom.ProjectionAngles))]); + disp(['DistanceOriginDetector ' num2str(geom.DistanceOriginDetector)]); + disp(['DistanceOriginSource ' num2str(geom.DistanceOriginSource)]); + + draw.draw_cone_geom(h_ax, geom, options); + case 'cone_vec' + disp('type: cone_vec'); + disp(['detector px: [' num2str(geom.DetectorRowCount), ', ' ... + num2str(geom.DetectorColCount) ']']); + disp(['# angles: ' num2str(size(geom.Vectors, 1))]); + + draw.draw_cone_vec_geom(h_ax, geom, options); + case 'fanflat' + disp('type: fanflat'); + disp(['detector px: ' num2str(geom.DetectorCount)]); + disp(['angle lo: ' num2str(geom.ProjectionAngles(1))]); + disp(['angle hi: ' num2str(geom.ProjectionAngles(end))]); + disp(['# angles: ' num2str(numel(geom.ProjectionAngles))]); + disp(['DistanceOriginDetector '... + num2str(geom.DistanceOriginDetector)]); + disp(['DistanceOriginSource '... + num2str(geom.DistanceOriginSource)]); + + draw.draw_fanflat_geom(h_ax, geom, options); + case 'fanflat_vec' + disp('type: fanflat_vec'); + disp(['detector px: ' num2str(geom.DetectorCount)]); + disp(['# angles: ' num2str(size(geom.Vectors, 1))]); + + draw.draw_fanflat_vec_geom(h_ax, geom, options); + otherwise + error(['Unknown geometry type ' geom.type]) + end + view(45, 25); % gives nicer default view angle + + function [options] = parseoptions(input_args) + % make an options struct + options = struct; + options.RotationAxis = [NaN, NaN, NaN]; + options.RotationAxisOffset = [0, 0, 0]; + options.VectorIdx = 1; + options.Color = 'k'; + options.DetectorMarker = '.'; + options.DetectorMarkerColor = ''; + options.DetectorLineColor = ''; + options.DetectorLineWidth = 1; + options.SourceMarker = '*'; + options.SourceMarkerColor = ''; + options.SourceDistance = 100; + options.OpticalAxisColor = ''; + options = parseargs.parseargs(options, input_args{:}); + + % if the color is still empty, replace by global color + if strcmpi(options.DetectorMarkerColor , '') + options.DetectorMarkerColor = options.Color; + end + if strcmpi(options.DetectorLineColor , '') + options.DetectorLineColor = options.Color; + end + if strcmpi(options.SourceMarkerColor , '') + options.SourceMarkerColor = options.Color; + end + if strcmpi(options.OpticalAxisColor , '') + options.OpticalAxisColor = options.Color; + end + end +end diff --git a/matlab/algorithms/plot_geom/+draw/draw_vol_geom.m b/matlab/algorithms/plot_geom/+draw/draw_vol_geom.m new file mode 100644 index 0000000..65126c2 --- /dev/null +++ b/matlab/algorithms/plot_geom/+draw/draw_vol_geom.m @@ -0,0 +1,105 @@ +function [] = draw_vol_geom( vol_geom, varargin) +%% draw_vol_geom.m +% brief rendering function for astra volume geometries +% describing a phantom. +% param vol_geom volume geometry describing the phantom +% param vx_size voxel size in unit of preference. must be same unit +% that was used to scale the projection geometry +% ------------------------------ +% optional parameters that can be provided as string value pairs: +% +% param Magnification magnification factor for the phantom. for small +% phantoms it might be necessary to scale the render +% up as otherwise it won't show up in the plot. +% Default = 1 +% param LineWidth line width for the box wireframe. Default = 2 +% param Color color of the wireframe. Default = 'r' +% +% date 20.06.2018 +% author Tim Elberfeld +% imec VisionLab +% University of Antwerp +% +% - last update 16.11.2018 +%% + h_ax = gca; + + if mod(size(varargin), 2) ~= 0 + vx_size = varargin{1}; + varargin = varargin(2:end); % consumed vx_size from arg list + else + vx_size = 1; + end + + + options = struct; + options.Color = 'r'; + options.LineWidth = 2; + options.Magnification = 1; + options = parseargs.parseargs(options, varargin{:}); + + hold on; + phantom_height = vol_geom.GridRowCount * vx_size; + phantom_width = vol_geom.GridColCount * vx_size; + phantom_depth = vol_geom.GridSliceCount * vx_size; + + if isfield(vol_geom, 'option') + minx = vol_geom.option.WindowMinX * vx_size; + maxx = vol_geom.option.WindowMaxX * vx_size; + + miny = vol_geom.option.WindowMinY * vx_size; + maxy = vol_geom.option.WindowMaxY * vx_size; + + minz = vol_geom.option.WindowMinZ * vx_size; + maxz = vol_geom.option.WindowMaxZ * vx_size; + else + minx = phantom_width / 2 * vx_size; + maxx = phantom_width / 2 * vx_size; + + miny = phantom_height / 2 * vx_size; + maxy = phantom_height / 2 * vx_size; + + minz = phantom_depth / 2 * vx_size; + maxz = phantom_depth / 2 * vx_size; + end + + xx_phantom = options.Magnification*[minx, minx, minx, minx, maxx, maxx, maxx, maxx]; + yy_phantom = options.Magnification*[miny, miny, maxy, maxy, miny, miny, maxy, maxy]; + zz_phantom = options.Magnification*[minz, maxz, minz, maxz, minz, maxz, minz, maxz]; + + face1 = [xx_phantom(1:4); yy_phantom(1:4); zz_phantom(1:4)]; + face2 = [[xx_phantom(1:2), xx_phantom(5:6)];... + [yy_phantom(1:2), yy_phantom(5:6)];... + [zz_phantom(1:2), zz_phantom(5:6)]]; + face3 = [[xx_phantom(3:4), xx_phantom(7:8)];... + [yy_phantom(3:4), yy_phantom(7:8)];... + [zz_phantom(3:4), zz_phantom(7:8)]]; + face4 = [[xx_phantom(5:6), xx_phantom(7:8)];... + [yy_phantom(5:6), yy_phantom(7:8)];... + [zz_phantom(5:6), zz_phantom(7:8)]]; + + % as we draw only a wire frame, we need only to draw 4 of the faces + draw_face(h_ax, face1, options); + draw_face(h_ax, face2, options); + draw_face(h_ax, face3, options); + draw_face(h_ax, face4, options); + + hold off; + + function [] = draw_face(h_ax, face_coords, options) + line(h_ax, face_coords(1, 1:2), face_coords(2, 1:2),... + face_coords(3, 1:2), 'LineWidth', options.LineWidth,... + 'Color', options.Color); + line(h_ax, face_coords(1, 3:4), face_coords(2, 3:4),... + face_coords(3, 3:4), 'LineWidth', options.LineWidth,... + 'Color', options.Color); + line(h_ax, [face_coords(1, 4),face_coords(1, 2)],... + [face_coords(2, 4),face_coords(2, 2)],... + [face_coords(3, 4),face_coords(3, 2)],... + 'LineWidth', options.LineWidth, 'Color', options.Color); + line(h_ax, [face_coords(1, 1),face_coords(1, 3)],... + [face_coords(2, 1),face_coords(2, 3)],... + [face_coords(3, 1),face_coords(3, 3)],... + 'LineWidth', options.LineWidth, 'Color', options.Color); + end +end diff --git a/matlab/algorithms/plot_geom/+draw/private/eucl_dist3d.m b/matlab/algorithms/plot_geom/+draw/private/eucl_dist3d.m new file mode 100644 index 0000000..57d6cb4 --- /dev/null +++ b/matlab/algorithms/plot_geom/+draw/private/eucl_dist3d.m @@ -0,0 +1,16 @@ +function [dist] = eucl_dist3d(a, b) +%% eucl_dist3d.m +% 3d euclidean distance for a nx3 matrix holding n 3-vectors +% param a - first vectors +% param b - second vectors +% date 07.11.2018 +% author Tim Elberfeld +% imec VisionLab +% University of Antwerp +% last update 07.11.2018 +%% + dist = sqrt((a(:, 1) - b(:, 1)).^2 + ... + (a(:, 2) - b(:, 2)).^2 + ... + (a(:, 3) - b(:, 3)).^2); +end + diff --git a/matlab/algorithms/plot_geom/+parseargs/license.txt b/matlab/algorithms/plot_geom/+parseargs/license.txt new file mode 100644 index 0000000..cae4cd5 --- /dev/null +++ b/matlab/algorithms/plot_geom/+parseargs/license.txt @@ -0,0 +1,27 @@ +Copyright (c) 2016, The MathWorks, Inc. +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are +met: + + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + * Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in + the documentation and/or other materials provided with the distribution. + * In all cases, the software is, and all modifications and derivatives + of the software shall be, licensed to you solely for use in conjunction + with MathWorks products and service offerings. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE +LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +POSSIBILITY OF SUCH DAMAGE. diff --git a/matlab/algorithms/plot_geom/+parseargs/parseargs.m b/matlab/algorithms/plot_geom/+parseargs/parseargs.m new file mode 100644 index 0000000..05ad613 --- /dev/null +++ b/matlab/algorithms/plot_geom/+parseargs/parseargs.m @@ -0,0 +1,97 @@ +function X = parseargs(X,varargin) +%PARSEARGS - Parses name-value pairs +% +% Behaves like setfield, but accepts multiple name-value pairs and provides +% some additional features: +% 1) If any field of X is an cell-array of strings, it can only be set to +% one of those strings. If no value is specified for that field, the +% first string is selected. +% 2) Where the field is not empty, its data type cannot be changed +% 3) Where the field contains a scalar, its size cannot be changed. +% +% X = parseargs(X,name1,value1,name2,value2,...) +% +% Intended for use as an argument parser for functions which multiple options. +% Example usage: +% +% function my_function(varargin) +% X.StartValue = 0; +% X.StopOnError = false; +% X.SolverType = {'fixedstep','variablestep'}; +% X.OutputFile = 'out.txt'; +% X = parseargs(X,varargin{:}); +% +% Then call (e.g.): +% +% my_function('OutputFile','out2.txt','SolverType','variablestep'); + +% The various #ok comments below are to stop MLint complaining about +% inefficient usage. In all cases, the inefficient usage (of error, getfield, +% setfield and find) is used to ensure compatibility with earlier versions +% of MATLAB. + +% Copyright 2006-2010 The MathWorks, Inc. + +remaining = nargin-1; % number of arguments other than X +count = 1; +fields = fieldnames(X); +modified = zeros(size(fields)); +% Take input arguments two at a time until we run out. +while remaining>=2 + fieldname = varargin{count}; + fieldind = find(strcmp(fieldname,fields)); + if ~isempty(fieldind) + oldvalue = getfield(X,fieldname); %#ok + newvalue = varargin{count+1}; + if iscell(oldvalue) + % Cell arrays must contain strings, and the new value must be + % a string which appears in the list. + if ~iscellstr(oldvalue) + error(sprintf('All allowed values for "%s" must be strings',fieldname)); %#ok + end + if ~ischar(newvalue) + error(sprintf('New value for "%s" must be a string',fieldname)); %#ok + end + if isempty(find(strcmp(oldvalue,newvalue))) %#ok + error(sprintf('"%s" is not allowed for field "%s"',newvalue,fieldname)); %#ok + end + elseif ~isempty(oldvalue) + % The caller isn't allowed to change the data type of a non-empty property, + % and scalars must remain as scalars. + if ~strcmp(class(oldvalue),class(newvalue)) + error(sprintf('Cannot change class of field "%s" from "%s" to "%s"',... + fieldname,class(oldvalue),class(newvalue))); %#ok + elseif numel(oldvalue)==1 & numel(newvalue)~=1 %#ok + error(sprintf('New value for "%s" must be a scalar',fieldname)); %#ok + end + end + X = setfield(X,fieldname,newvalue); %#ok + modified(fieldind) = 1; + else + error(['Not a valid field name: ' fieldname]); + end + remaining = remaining - 2; + count = count + 2; +end +% Check that we had a value for every name. +if remaining~=0 + error('Odd number of arguments supplied. Name-value pairs required'); +end + +% Now find cell arrays which were not modified by the above process, and select +% the first string. +notmodified = find(~modified); +for i=1:length(notmodified) + fieldname = fields{notmodified(i)}; + oldvalue = getfield(X,fieldname); %#ok + if iscell(oldvalue) + if ~iscellstr(oldvalue) + error(sprintf('All allowed values for "%s" must be strings',fieldname)); %#ok + elseif isempty(oldvalue) + error(sprintf('Empty cell array not allowed for field "%s"',fieldname)); %#ok + end + X = setfield(X,fieldname,oldvalue{1}); %#ok + end +end + + diff --git a/matlab/algorithms/plot_geom/+stlTools/license.txt b/matlab/algorithms/plot_geom/+stlTools/license.txt new file mode 100644 index 0000000..b2d00b0 --- /dev/null +++ b/matlab/algorithms/plot_geom/+stlTools/license.txt @@ -0,0 +1,34 @@ +Copyright (c) 2017, Pau Micó +Copyright (c) 2015, Sven Holcombe +Copyright (c) 2011, Eric Johnson +Copyright (c) 2013, Adam H. Aitkenhead +Copyright (c) 2011, Francis Esmonde-White +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are +met: + + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + * Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in + the documentation and/or other materials provided with the distribution + * Neither the name of the The MathWorks, Inc. nor the names + of its contributors may be used to endorse or promote products derived + from this software without specific prior written permission. + * Neither the name of the The Christie NHS Foundation Trust nor the names + of its contributors may be used to endorse or promote products derived + from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE +LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +POSSIBILITY OF SUCH DAMAGE. diff --git a/matlab/algorithms/plot_geom/+stlTools/readme.txt b/matlab/algorithms/plot_geom/+stlTools/readme.txt new file mode 100644 index 0000000..b8037d6 --- /dev/null +++ b/matlab/algorithms/plot_geom/+stlTools/readme.txt @@ -0,0 +1,16 @@ +The 'stlTools' toolbox is a collection of functions, samples and demos to illustrate how to deal with STL files. Some of them are contributions published in Matlab Central and properly referenced here. +This toolbox contains the following files: + +stlGetFormat: identifies the format of the STL file and returns 'binary' or 'ascii'. This file is inspired in the 'READ-stl' file written and published by Adam H. Aitkenhead (http://www.mathworks.com/matlabcentral/fileexchange/27390-mesh-voxelisation). Copyright (c) 2013, Adam H. Aitkenhead. +stlReadAscii: reads an STL file written in ascii format. This file is inspired in the 'READ-stl' file written and published by Adam H. Aitkenhead (http://www.mathworks.com/matlabcentral/fileexchange/27390-mesh-voxelisation). Copyright (c) 2013, Adam H. Aitkenhead +stlReadBinary: reads an STL file written in binary format. This file is inspired in the 'READ-stl' file written and published by Adam H. Aitkenhead (http://www.mathworks.com/matlabcentral/fileexchange/27390-mesh-voxelisation). Copyright (c) 2013, Adam H. Aitkenhead +stlRead: uses 'stlGetFormat', 'stlReadAscii' and 'stlReadBinary' to make STL reading independent of the format of the file +stlWrite: writes an STL file in 'ascii' or 'binary' formats. This is written and published by Sven Holcombe (http://www.mathworks.com/matlabcentral/fileexchange/20922-stlwrite-filename--varargin-). Copyright (c) 2012, Grant Lohsen. Copyright (c) 2015, Sven Holcombe. +stlSlimVerts: finds and removes duplicated vertices. This function is written and published by Francis Esmonde-White as PATCHSLIM (http://www.mathworks.com/matlabcentral/fileexchange/29986-patch-slim--patchslim-m-). Copyright (c) 2011, Francis Esmonde-White. +stlGetVerts: returns a list of vertices that are 'opened' or 'closed' depending on the 'mode' input parameter. An 'open' vertice is the one that defines an open side. An open side is the one that only takes part of one triangle +stlDelVerts: removes a list of vertices from STL files +stlAddVerts: adds the new vertices from a list (and consequently, new faces) to a STL object +stlPlot: is an easy way to plot an STL object +stlDemo: is a collection of examples about how to use stlTools +femur_binary: is an ascii STL sample used in 'stlDemo'. It is published by Eric Johnson (http://www.mathworks.com/matlabcentral/fileexchange/22409-stl-file-reader). Copyright (c) 2011, Eric Johnson. +sphere_ascii: is a binary STL sample
\ No newline at end of file diff --git a/matlab/algorithms/plot_geom/+stlTools/stlAddVerts.m b/matlab/algorithms/plot_geom/+stlTools/stlAddVerts.m new file mode 100644 index 0000000..1894cd7 --- /dev/null +++ b/matlab/algorithms/plot_geom/+stlTools/stlAddVerts.m @@ -0,0 +1,17 @@ +function [vnew, fnew] = stlAddVerts(v, f, list) +%STLADDVERTS adds new vertices (and consequently, new faces) to a STL object +%V is the Nx3 array of vertices +%F is the Mx3 array of faces +%LIST is the list of vertices to be added to the object +%VNEW is the new array of vertices +%FNEW is the new array of faces + +import stlTools.* + +% triangulation just with the slice +faces = delaunay(list(:,1),list(:,2)); % calculate new faces +% update object +nvert = length(v); % number of original vertices +v = [v; list]; % update vertices with the ones in the list +f = [f; faces+nvert]; % update faces with the new ones +[vnew,fnew] = stlSlimVerts(v,f); % clear repeated vertices
\ No newline at end of file diff --git a/matlab/algorithms/plot_geom/+stlTools/stlDelVerts.m b/matlab/algorithms/plot_geom/+stlTools/stlDelVerts.m new file mode 100644 index 0000000..718ac53 --- /dev/null +++ b/matlab/algorithms/plot_geom/+stlTools/stlDelVerts.m @@ -0,0 +1,26 @@ +function [vnew, fnew] = stlDelVerts(v, f, list) +%STLDELVERT removes a list of vertices from STL files +%V is the Nx3 array of vertices +%F is the Mx3 array of faces +%LIST are the vertices (rows) to delete, where length(LIST) < N +%VNEW is the new array of vertices +%FNEW is the new array of faces + +% find (on the global set) the position (rows) of the vertices to be deleted +[~,vdel] = ismember(list,v,'rows'); + +% delete vertices and get new tags +vnew = v; +tags = 1:length(v); +vnew(vdel,:) = []; +tags(vdel) = []; + +% delete faces +fnew = f; +[fdel,~] = find(ismember(f,vdel)); % find the position (rows) of the faces to delete +fnew(fdel,:) = []; + +% rename faces, as some of the vertices have been deleted +flist = reshape(fnew,numel(fnew),1); +[~,ind] = ismember(flist,tags); +fnew = reshape(ind,numel(flist)/3,3);
\ No newline at end of file diff --git a/matlab/algorithms/plot_geom/+stlTools/stlDemo.m b/matlab/algorithms/plot_geom/+stlTools/stlDemo.m new file mode 100644 index 0000000..af46315 --- /dev/null +++ b/matlab/algorithms/plot_geom/+stlTools/stlDemo.m @@ -0,0 +1,41 @@ +%% STLDEMO shows how to use the functions included in the toolbox STLTOOLS + +import stlTools.* + +%% EXAMPLE 1.- How to cut a sphere and close the base to get a semisphere + +% load an ascii STL sample file (STLGETFORMAT and STLREADASCII) +[vertices,faces,normals,name] = stlRead('sphere300faces.stl'); +stlPlot(vertices,faces,name); + +% the sphere is centered in the origin +% now we get a list of vertices to be deleted if (x,y,z<0) +minZ = 0; +[rows, ~] = find(vertices(:,3) < minZ); +list = vertices(rows,:); + +% if we delete the list of vertices with z<0, we get an opened semisphere +% (as the base is not closed) +[newv,newf] = stlDelVerts(vertices,faces,list); +stlPlot(newv,newf,name); + +% the next step is to identify a new list with the faces that are opened +% (that means all the sides that belong only to a unique triangle) +list = stlGetVerts(newv,newf,'opened'); + +% finally we generate all the new faces that are needed just to close the +% base of the semisphere +[vsemi,fsemi] = stlAddVerts(newv,newf,list); +stlPlot(vsemi,fsemi,'closed semisphere'); + +%% EXAMPLE 2.- How to get a section of a femur + +[vertices,faces,normals,name] = stlRead('femur_binary.stl'); +stlPlot(vertices,faces,name); + +minX = 1.2; +[rows, ~] = find(vertices(:,1) < minX); +list = vertices(rows,:); + +[newv,newf] = stlDelVerts(vertices,faces,list); +stlPlot(newv,newf,'section of the femur'); diff --git a/matlab/algorithms/plot_geom/+stlTools/stlGetFormat.m b/matlab/algorithms/plot_geom/+stlTools/stlGetFormat.m new file mode 100644 index 0000000..6901b26 --- /dev/null +++ b/matlab/algorithms/plot_geom/+stlTools/stlGetFormat.m @@ -0,0 +1,38 @@ +function format = stlGetFormat(fileName) +%STLGETFORMAT identifies the format of the STL file and returns 'binary' or +%'ascii' + +fid = fopen(fileName); +% Check the file size first, since binary files MUST have a size of 84+(50*n) +fseek(fid,0,1); % Go to the end of the file +fidSIZE = ftell(fid); % Check the size of the file +if rem(fidSIZE-84,50) > 0 + format = 'ascii'; +else + % Files with a size of 84+(50*n), might be either ascii or binary... + % Read first 80 characters of the file. + % For an ASCII file, the data should begin immediately (give or take a few + % blank lines or spaces) and the first word must be 'solid'. + % For a binary file, the first 80 characters contains the header. + % It is bad practice to begin the header of a binary file with the word + % 'solid', so it can be used to identify whether the file is ASCII or + % binary. + fseek(fid,0,-1); % go to the beginning of the file + header = strtrim(char(fread(fid,80,'uchar')')); % trim leading and trailing spaces + isSolid = strcmp(header(1:min(5,length(header))),'solid'); % take first 5 char + fseek(fid,-80,1); % go to the end of the file minus 80 characters + tail = char(fread(fid,80,'uchar')'); + isEndSolid = findstr(tail,'endsolid'); + + % Double check by reading the last 80 characters of the file. + % For an ASCII file, the data should end (give or take a few + % blank lines or spaces) with 'endsolid <object_name>'. + % If the last 80 characters contains the word 'endsolid' then this + % confirms that the file is indeed ASCII. + if isSolid & isEndSolid + format = 'ascii'; + else + format = 'binary'; + end +end +fclose(fid);
\ No newline at end of file diff --git a/matlab/algorithms/plot_geom/+stlTools/stlGetVerts.m b/matlab/algorithms/plot_geom/+stlTools/stlGetVerts.m new file mode 100644 index 0000000..179eff3 --- /dev/null +++ b/matlab/algorithms/plot_geom/+stlTools/stlGetVerts.m @@ -0,0 +1,31 @@ +function list = stlGetVerts(v, f, mode) +%GETVERTS returns the vertices that are 'opened' or 'closed' depending on +%the 'mode'. An 'open' vertice is the one that defines an open side. An +%open side is the one that only takes part of one triangle +%V is the Nx3 array of vertices +%F is the Mx3 array of faces +%MODE can be 'opened' or 'closed' depending of the kind of vertices to list +%LIST is the list of 'opened' or 'closed' vertices + +sides = sort([[f(:,1) f(:,2)]; ... + [f(:,2) f(:,3)]; ... + [f(:,3) f(:,1)]],2); + +[C,ia,ic] = unique(sides,'rows'); +ind_all = sort(ic); % open and closed sides +ind_rep = find(diff(ind_all) == 0); +ind_cls = ind_all(ind_rep); % closed sides +sides_cls = C(ind_cls,:); +ind_rep = [ind_rep; ind_rep+1]; +ind_opn = ind_all; +ind_opn(ind_rep) = []; % open sides +sides_opn = C(ind_opn,:); + +switch mode, + case'opened', + list = v(unique(sides_opn(:)),:); + case 'closed', + list = v(unique(sides_cls(:)),:); + otherwise, + error('getVerts:InvalidMode','The ''mode'' valid values are ''opened'' or ''closed'''); +end
\ No newline at end of file diff --git a/matlab/algorithms/plot_geom/+stlTools/stlPlot.m b/matlab/algorithms/plot_geom/+stlTools/stlPlot.m new file mode 100644 index 0000000..e9f17c5 --- /dev/null +++ b/matlab/algorithms/plot_geom/+stlTools/stlPlot.m @@ -0,0 +1,23 @@ +function stlPlot(v, f, name) +%STLPLOT is an easy way to plot an STL object +%V is the Nx3 array of vertices +%F is the Mx3 array of faces +%NAME is the name of the object, that will be displayed as a title + +figure; +object.vertices = v; +object.faces = f; +patch(object,'FaceColor', [0.8 0.8 1.0], ... + 'EdgeColor', 'none', ... + 'FaceLighting', 'gouraud', ... + 'AmbientStrength', 0.15); + +% Add a camera light, and tone down the specular highlighting +camlight('headlight'); +material('dull'); + +% Fix the axes scaling, and set a nice view angle +axis('image'); +view([-135 35]); +grid on; +title(name); diff --git a/matlab/algorithms/plot_geom/+stlTools/stlRead.m b/matlab/algorithms/plot_geom/+stlTools/stlRead.m new file mode 100644 index 0000000..5775e50 --- /dev/null +++ b/matlab/algorithms/plot_geom/+stlTools/stlRead.m @@ -0,0 +1,15 @@ +function [v, f, n, name] = stlRead(fileName) +%STLREAD reads any STL file not depending on its format +%V are the vertices +%F are the faces +%N are the normals +%NAME is the name of the STL object (NOT the name of the STL file) + +import stlTools.* + +format = stlGetFormat(fileName); +if strcmp(format,'ascii') + [v,f,n,name] = stlReadAscii(fileName); +elseif strcmp(format,'binary') + [v,f,n,name] = stlReadBinary(fileName); +end
\ No newline at end of file diff --git a/matlab/algorithms/plot_geom/+stlTools/stlReadAscii.m b/matlab/algorithms/plot_geom/+stlTools/stlReadAscii.m new file mode 100644 index 0000000..b035aa8 --- /dev/null +++ b/matlab/algorithms/plot_geom/+stlTools/stlReadAscii.m @@ -0,0 +1,58 @@ +function [v, f, n, name] = stlReadAscii(fileName) +%STLREADASCII reads a STL file written in ASCII format +%V are the vertices +%F are the faces +%N are the normals +%NAME is the name of the STL object (NOT the name of the STL file) + +import stlTools.* + +%====================== +% STL ascii file format +%====================== +% ASCII STL files have the following structure. Technically each facet +% could be any 2D shape, but in practice only triangular facets tend to be +% used. The present code ONLY works for meshes composed of triangular +% facets. +% +% solid object_name +% facet normal x y z +% outer loop +% vertex x y z +% vertex x y z +% vertex x y z +% endloop +% endfacet +% +% <Repeat for all facets...> +% +% endsolid object_name + +fid = fopen(fileName); +cellcontent = textscan(fid,'%s','delimiter','\n'); % read all the file and put content in cells +content = cellcontent{:}(logical(~strcmp(cellcontent{:},''))); % remove all blank lines +fclose(fid); + +% read the STL name +line1 = char(content(1)); +if (size(line1,2) >= 7) + name = line1(7:end); +else + name = 'Unnamed Object'; +end + +% read the vector normals +normals = char(content(logical(strncmp(content,'facet normal',12)))); +n = str2num(normals(:,13:end)); + +% read the vertex coordinates (vertices) +vertices = char(content(logical(strncmp(content,'vertex',6)))); +v = str2num(vertices(:,7:end)); +nvert = size(vertices,1); % number of vertices +nfaces = sum(strcmp(content,'endfacet')); % number of faces +if (nvert == 3*nfaces) + f = reshape(1:nvert,[3 nfaces])'; % create faces +end + +% slim the file (delete duplicated vertices) +[v,f] = stlSlimVerts(v,f);
\ No newline at end of file diff --git a/matlab/algorithms/plot_geom/+stlTools/stlReadBinary.m b/matlab/algorithms/plot_geom/+stlTools/stlReadBinary.m new file mode 100644 index 0000000..3e04f7a --- /dev/null +++ b/matlab/algorithms/plot_geom/+stlTools/stlReadBinary.m @@ -0,0 +1,62 @@ +function [v, f, n, name] = stlReadBinary(fileName) +%STLREADBINARY reads a STL file written in BINARY format +%V are the vertices +%F are the faces +%N are the normals +%NAME is the name of the STL object (NOT the name of the STL file) + +import stlTools.* + +%======================= +% STL binary file format +%======================= +% Binary STL files have an 84 byte header followed by 50-byte records, each +% describing a single facet of the mesh. Technically each facet could be +% any 2D shape, but that would screw up the 50-byte-per-facet structure, so +% in practice only triangular facets are used. The present code ONLY works +% for meshes composed of triangular facets. +% +% HEADER: +% 80 bytes: Header text +% 4 bytes: (int) The number of facets in the STL mesh +% +% DATA: +% 4 bytes: (float) normal x +% 4 bytes: (float) normal y +% 4 bytes: (float) normal z +% 4 bytes: (float) vertex1 x +% 4 bytes: (float) vertex1 y +% 4 bytes: (float) vertex1 z +% 4 bytes: (float) vertex2 x +% 4 bytes: (float) vertex2 y +% 4 bytes: (float) vertex2 z +% 4 bytes: (float) vertex3 x +% 4 bytes: (float) vertex3 y +% 4 bytes: (float) vertex3 z +% 2 bytes: Padding to make the data for each facet 50-bytes in length +% ...and repeat for next facet... + +fid = fopen(fileName); +header = fread(fid,80,'int8'); % reading header's 80 bytes +name = deblank(native2unicode(header,'ascii')'); +if isempty(name) + name = 'Unnamed Object'; % no object name in binary files! +end +nfaces = fread(fid,1,'int32'); % reading the number of facets in the stl file (next 4 byters) +nvert = 3*nfaces; % number of vertices +% reserve memory for vectors (increase the processing speed) +n = zeros(nfaces,3); +v = zeros(nvert,3); +f = zeros(nfaces,3); +for i = 1 : nfaces % read the data for each facet + tmp = fread(fid,3*4,'float'); % read coordinates + n(i,:) = tmp(1:3); % x,y,z components of the facet's normal vector + v(3*i-2,:) = tmp(4:6); % x,y,z coordinates of vertex 1 + v(3*i-1,:) = tmp(7:9); % x,y,z coordinates of vertex 2 + v(3*i,:) = tmp(10:12); % x,y,z coordinates of vertex 3 + f(i,:) = [3*i-2 3*i-1 3*i]; % face + fread(fid,1,'int16'); % Move to the start of the next facet (2 bytes of padding) +end +fclose(fid); +% slim the file (delete duplicated vertices) +[v,f] = stlSlimVerts(v,f);
\ No newline at end of file diff --git a/matlab/algorithms/plot_geom/+stlTools/stlSlimVerts.m b/matlab/algorithms/plot_geom/+stlTools/stlSlimVerts.m new file mode 100644 index 0000000..2781067 --- /dev/null +++ b/matlab/algorithms/plot_geom/+stlTools/stlSlimVerts.m @@ -0,0 +1,29 @@ +function [vnew, fnew]= stlSlimVerts(v, f) +% PATCHSLIM removes duplicate vertices in surface meshes. +% +% This function finds and removes duplicate vertices. +% +% USAGE: [v, f]=patchslim(v, f) +% +% Where v is the vertex list and f is the face list specifying vertex +% connectivity. +% +% v contains the vertices for all triangles [3*n x 3]. +% f contains the vertex lists defining each triangle face [n x 3]. +% +% This will reduce the size of typical v matrix by about a factor of 6. +% +% For more information see: +% http://www.esmonde-white.com/home/diversions/matlab-program-for-loading-stl-files +% +% Francis Esmonde-White, May 2010 + +if ~exist('v','var') + error('The vertex list (v) must be specified.'); +end +if ~exist('f','var') + error('The vertex connectivity of the triangle faces (f) must be specified.'); +end + +[vnew, indexm, indexn] = unique(v, 'rows'); +fnew = indexn(f); diff --git a/matlab/algorithms/plot_geom/+stlTools/stlWrite.m b/matlab/algorithms/plot_geom/+stlTools/stlWrite.m new file mode 100644 index 0000000..0ce5d46 --- /dev/null +++ b/matlab/algorithms/plot_geom/+stlTools/stlWrite.m @@ -0,0 +1,252 @@ +function stlWrite(filename, varargin) +%STLWRITE Write STL file from patch or surface data. +% +% STLWRITE(FILE, FV) writes a stereolithography (STL) file to FILE for a +% triangulated patch defined by FV (a structure with fields 'vertices' +% and 'faces'). +% +% STLWRITE(FILE, FACES, VERTICES) takes faces and vertices separately, +% rather than in an FV struct +% +% STLWRITE(FILE, X, Y, Z) creates an STL file from surface data in X, Y, +% and Z. STLWRITE triangulates this gridded data into a triangulated +% surface using triangulation options specified below. X, Y and Z can be +% two-dimensional arrays with the same size. If X and Y are vectors with +% length equal to SIZE(Z,2) and SIZE(Z,1), respectively, they are passed +% through MESHGRID to create gridded data. If X or Y are scalar values, +% they are used to specify the X and Y spacing between grid points. +% +% STLWRITE(...,'PropertyName',VALUE,'PropertyName',VALUE,...) writes an +% STL file using the following property values: +% +% MODE - File is written using 'binary' (default) or 'ascii'. +% +% TITLE - Header text (max 80 chars) written to the STL file. +% +% TRIANGULATION - When used with gridded data, TRIANGULATION is either: +% 'delaunay' - (default) Delaunay triangulation of X, Y +% 'f' - Forward slash division of grid quads +% 'b' - Back slash division of quadrilaterals +% 'x' - Cross division of quadrilaterals +% Note that 'f', 'b', or 't' triangulations now use an +% inbuilt version of FEX entry 28327, "mesh2tri". +% +% FACECOLOR - Single colour (1-by-3) or one-colour-per-face (N-by-3) +% vector of RGB colours, for face/vertex input. RGB range +% is 5 bits (0:31), stored in VisCAM/SolidView format +% (http://en.wikipedia.org/wiki/STL_(file_format)#Color_in_binary_STL) +% +% Example 1: +% % Write binary STL from face/vertex data +% tmpvol = false(20,20,20); % Empty voxel volume +% tmpvol(8:12,8:12,5:15) = 1; % Turn some voxels on +% fv = isosurface(~tmpvol, 0.5); % Make patch w. faces "out" +% stlwrite('test.stl',fv) % Save to binary .stl +% +% Example 2: +% % Write ascii STL from gridded data +% [X,Y] = deal(1:40); % Create grid reference +% Z = peaks(40); % Create grid height +% stlwrite('test.stl',X,Y,Z,'mode','ascii') +% +% Example 3: +% % Write binary STL with coloured faces +% cVals = fv.vertices(fv.faces(:,1),3); % Colour by Z height. +% cLims = [min(cVals) max(cVals)]; % Transform height values +% nCols = 255; cMap = jet(nCols); % onto an 8-bit colour map +% fColsDbl = interp1(linspace(cLims(1),cLims(2),nCols),cMap,cVals); +% fCols8bit = fColsDbl*255; % Pass cols in 8bit (0-255) RGB triplets +% stlwrite('testCol.stl',fv,'FaceColor',fCols8bit) + +% Original idea adapted from surf2stl by Bill McDonald. Huge speed +% improvements implemented by Oliver Woodford. Non-Delaunay triangulation +% of quadrilateral surface courtesy of Kevin Moerman. FaceColor +% implementation by Grant Lohsen. +% +% Author: Sven Holcombe, 11-24-11 + + +% Check valid filename path +path = fileparts(filename); +if ~isempty(path) && ~exist(path,'dir') + error('Directory "%s" does not exist.',path); +end + +% Get faces, vertices, and user-defined options for writing +[faces, vertices, options] = parseInputs(varargin{:}); +asciiMode = strcmp( options.mode ,'ascii'); + +% Create the facets +facets = single(vertices'); +facets = reshape(facets(:,faces'), 3, 3, []); + +% Compute their normals +V1 = squeeze(facets(:,2,:) - facets(:,1,:)); +V2 = squeeze(facets(:,3,:) - facets(:,1,:)); +normals = V1([2 3 1],:) .* V2([3 1 2],:) - V2([2 3 1],:) .* V1([3 1 2],:); +clear V1 V2 +normals = bsxfun(@times, normals, 1 ./ sqrt(sum(normals .* normals, 1))); +facets = cat(2, reshape(normals, 3, 1, []), facets); +clear normals + +% Open the file for writing +permissions = {'w','wb+'}; +fid = fopen(filename, permissions{asciiMode+1}); +if (fid == -1) + error('stlwrite:cannotWriteFile', 'Unable to write to %s', filename); +end + +% Write the file contents +if asciiMode + % Write HEADER + fprintf(fid,'solid %s\r\n',options.title); + % Write DATA + fprintf(fid,[... + 'facet normal %.7E %.7E %.7E\r\n' ... + 'outer loop\r\n' ... + 'vertex %.7E %.7E %.7E\r\n' ... + 'vertex %.7E %.7E %.7E\r\n' ... + 'vertex %.7E %.7E %.7E\r\n' ... + 'endloop\r\n' ... + 'endfacet\r\n'], facets); + % Write FOOTER + fprintf(fid,'endsolid %s\r\n',options.title); + +else % BINARY + % Write HEADER + fprintf(fid, '%-80s', options.title); % Title + fwrite(fid, size(facets, 3), 'uint32'); % Number of facets + % Write DATA + % Add one uint16(0) to the end of each facet using a typecasting trick + facets = reshape(typecast(facets(:), 'uint16'), 12*2, []); + % Set the last bit to 0 (default) or supplied RGB + facets(end+1,:) = options.facecolor; + fwrite(fid, facets, 'uint16'); +end + +% Close the file +fclose(fid); +%fprintf('Wrote %d facets\n',size(facets, 2)); + + +%% Input handling subfunctions +function [faces, vertices, options] = parseInputs(varargin) +% Determine input type +if isstruct(varargin{1}) % stlwrite('file', FVstruct, ...) + if ~all(isfield(varargin{1},{'vertices','faces'})) + error( 'Variable p must be a faces/vertices structure' ); + end + faces = varargin{1}.faces; + vertices = varargin{1}.vertices; + options = parseOptions(varargin{2:end}); + +elseif isnumeric(varargin{1}) + firstNumInput = cellfun(@isnumeric,varargin); + firstNumInput(find(~firstNumInput,1):end) = 0; % Only consider numerical input PRIOR to the first non-numeric + numericInputCnt = nnz(firstNumInput); + + options = parseOptions(varargin{numericInputCnt+1:end}); + switch numericInputCnt + case 3 % stlwrite('file', X, Y, Z, ...) + % Extract the matrix Z + Z = varargin{3}; + + % Convert scalar XY to vectors + ZsizeXY = fliplr(size(Z)); + for i = 1:2 + if isscalar(varargin{i}) + varargin{i} = (0:ZsizeXY(i)-1) * varargin{i}; + end + end + + % Extract X and Y + if isequal(size(Z), size(varargin{1}), size(varargin{2})) + % X,Y,Z were all provided as matrices + [X,Y] = varargin{1:2}; + elseif numel(varargin{1})==ZsizeXY(1) && numel(varargin{2})==ZsizeXY(2) + % Convert vector XY to meshgrid + [X,Y] = meshgrid(varargin{1}, varargin{2}); + else + error('stlwrite:badinput', 'Unable to resolve X and Y variables'); + end + + % Convert to faces/vertices + if strcmp(options.triangulation,'delaunay') + faces = delaunay(X,Y); + vertices = [X(:) Y(:) Z(:)]; + else + if ~exist('mesh2tri','file') + error('stlwrite:missing', '"mesh2tri" is required to convert X,Y,Z matrices to STL. It can be downloaded from:\n%s\n',... + 'http://www.mathworks.com/matlabcentral/fileexchange/28327') + end + [faces, vertices] = mesh2tri(X, Y, Z, options.triangulation); + end + + case 2 % stlwrite('file', FACES, VERTICES, ...) + faces = varargin{1}; + vertices = varargin{2}; + + otherwise + error('stlwrite:badinput', 'Unable to resolve input types.'); + end +end + +if ~isempty(options.facecolor) % Handle colour preparation + facecolor = uint16(options.facecolor); + %Set the Valid Color bit (bit 15) + c0 = bitshift(ones(size(faces,1),1,'uint16'),15); + %Red color (10:15), Blue color (5:9), Green color (0:4) + c0 = bitor(bitshift(bitand(2^6-1, facecolor(:,1)),10),c0); + c0 = bitor(bitshift(bitand(2^11-1, facecolor(:,2)),5),c0); + c0 = bitor(bitand(2^6-1, facecolor(:,3)),c0); + options.facecolor = c0; +else + options.facecolor = 0; +end + +function options = parseOptions(varargin) +IP = inputParser; +IP.addParamValue('mode', 'binary', @ischar) +IP.addParamValue('title', sprintf('Created by stlwrite.m %s',datestr(now)), @ischar); +IP.addParamValue('triangulation', 'delaunay', @ischar); +IP.addParamValue('facecolor',[], @isnumeric) +IP.addParamValue('facecolour',[], @isnumeric) +IP.parse(varargin{:}); +options = IP.Results; +if ~isempty(options.facecolour) + options.facecolor = options.facecolour; +end + +function [F,V]=mesh2tri(X,Y,Z,tri_type) +% function [F,V]=mesh2tri(X,Y,Z,tri_type) +% +% Available from http://www.mathworks.com/matlabcentral/fileexchange/28327 +% Included here for convenience. Many thanks to Kevin Mattheus Moerman +% kevinmoerman@hotmail.com +% 15/07/2010 +%------------------------------------------------------------------------ + +[J,I]=meshgrid(1:1:size(X,2)-1,1:1:size(X,1)-1); + +switch tri_type + case 'f'%Forward slash + TRI_I=[I(:),I(:)+1,I(:)+1; I(:),I(:),I(:)+1]; + TRI_J=[J(:),J(:)+1,J(:); J(:),J(:)+1,J(:)+1]; + F = sub2ind(size(X),TRI_I,TRI_J); + case 'b'%Back slash + TRI_I=[I(:),I(:)+1,I(:); I(:)+1,I(:)+1,I(:)]; + TRI_J=[J(:)+1,J(:),J(:); J(:)+1,J(:),J(:)+1]; + F = sub2ind(size(X),TRI_I,TRI_J); + case 'x'%Cross + TRI_I=[I(:)+1,I(:); I(:)+1,I(:)+1; I(:),I(:)+1; I(:),I(:)]; + TRI_J=[J(:),J(:); J(:)+1,J(:); J(:)+1,J(:)+1; J(:),J(:)+1]; + IND=((numel(X)+1):numel(X)+prod(size(X)-1))'; + F = sub2ind(size(X),TRI_I,TRI_J); + F(:,3)=repmat(IND,[4,1]); + Fe_I=[I(:),I(:)+1,I(:)+1,I(:)]; Fe_J=[J(:),J(:),J(:)+1,J(:)+1]; + Fe = sub2ind(size(X),Fe_I,Fe_J); + Xe=mean(X(Fe),2); Ye=mean(Y(Fe),2); Ze=mean(Z(Fe),2); + X=[X(:);Xe(:)]; Y=[Y(:);Ye(:)]; Z=[Z(:);Ze(:)]; +end + +V=[X(:),Y(:),Z(:)];
\ No newline at end of file diff --git a/matlab/algorithms/plot_geom/astra_create_example_cone.m b/matlab/algorithms/plot_geom/astra_create_example_cone.m new file mode 100644 index 0000000..598899f --- /dev/null +++ b/matlab/algorithms/plot_geom/astra_create_example_cone.m @@ -0,0 +1,165 @@ +function [ proj_geom ] = create_example_cone(type) +%% create_example_cone.m +% create an example standard cone beam geometry +% param type type of geometry to create. provided as one of the +% following strings (not case sensitive): +% 'normal' - standard (non vector) geometry +% 'vec' - example vector geometry +% 'helix' - helical trajectory vector geometry +% 'deform_vec' - deformed vector geometry example, +% courtesy of Van Nguyen +% return proj_geom - the geometry that was created +% +% date 21.06.2018 +% author Tim Elberfeld, Alice Presenti, Van Nguyen +% imec VisionLab +% University of Antwerp +% last mod 07.11.2018 +%% + if(nargin < 1) + type = 'normal'; + end + if strcmpi(type, 'normal') + % first, give measurements in mm + det_spacing = [0.035, 0.035]; + detector_px = [1200, 1200]; + angles = linspace2(0, 2*pi, 100); + source_origin = 30; + origin_det = 200; + phantom_size = 5; + + phantom_px = 150; % voxels for the phantom + vx_size = phantom_size / phantom_px; % voxel size + + % now express all measurements in terms of the voxel size + det_spacing = det_spacing ./ vx_size; + origin_det = origin_det ./ vx_size; + source_origin = source_origin ./ vx_size; + + proj_geom = astra_create_proj_geom('cone', det_spacing(1), ... + det_spacing(2), detector_px(1), detector_px(2), angles,... + source_origin, origin_det); + + elseif strcmpi(type, 'vec') + proj_geom = create_standard_vec_geom(); + + elseif strcmpi(type, 'deform_vec') + proj_geom = create_deform_vec_geom(); + + elseif strcmpi(type, 'helix') + detRows = 220; + detCols = 220; + detPitchX = 10^-3 * 11; % microns to mm + detPitchY = 10^-3 * 11; % microns to mm + objSrcDist = 500; % mm + rotStep = 3.6; % deg + numAngles = 200; + zTranslation = 0.5; % mm + zDist = zTranslation * numAngles; + vectors = zeros(numAngles, 12); + translation = -zDist + (zTranslation * (0:(numAngles-1))); + + minAngle = 0; % just assume we start at 0 + maxAngle = (numAngles * rotStep)/180 * pi; % convert deg to rad + angles = linspace(minAngle, maxAngle, 200); + + % source position per angle + vectors(:, 1) = sin(angles) * objSrcDist; + vectors(:, 2) = -cos(angles) * objSrcDist; + vectors(:, 3) = translation; + + % detector position per angle + vectors(:, 4) = 0; + vectors(:, 5) = 0; + vectors(:, 6) = translation; + + % vector from detector pixel (0,0) to (0,1) + vectors(:, 7) = cos(angles); + vectors(:, 8) = sin(angles); + vectors(:, 9) = zeros(numAngles, 1); + + % vector from detector pixel (0,0) to (1, 0) + vectors(:, 10) = zeros(numAngles, 1); + vectors(:, 11) = zeros(numAngles, 1); + vectors(:, 12) = ones(numAngles, 1); + + proj_geom = astra_create_proj_geom('cone_vec', detCols, detRows, vectors); + + else + proj_geom = create_standard_vec_geom(); + end + + function [proj_geom, z_axis] = create_standard_vec_geom() + % geometry settings taken from code from Alice Presenti + settings = struct; + settings.detectorPixelSize = 0.0748; + % detector size and number of projections + settings.projectionSize = [1536 1944 21]; + settings.SOD = 679.238020; %[mm] + settings.SDD = 791.365618; %[mm] + settings.voxelSize = 1; %[mm] + settings.gamma = linspace(0,300,21)*pi/180; + + S0 = zeros(settings.projectionSize(3), 12); + Sorig = [-settings.SOD, 0, 0,... % the ray origin vector (source) + (settings.SDD-settings.SOD), 0, 0,... % detector center + 0, -settings.detectorPixelSize, 0,... % detector u axis + 0, 0, settings.detectorPixelSize]; % detector v axis + + z_axis = [0, 0, 1]; + for i = 1:settings.projectionSize(3) + S0(i,:) = Sorig(:); + S0(i,1:3) = rotate_around3d(S0(i,1:3),... + z_axis, settings.gamma(i)); + S0(i,4:6) = rotate_around3d(S0(i,4:6), z_axis,... + settings.gamma(i)); + S0(i,7:9) = rotate_around3d(S0(i,7:9), z_axis,... + settings.gamma(i)); + S0(i,10:12) = rotate_around3d(S0(i,10:12),... + z_axis, settings.gamma(i)); + end + proj_geom = astra_create_proj_geom('cone_vec', ... + settings.projectionSize(2), settings.projectionSize(1), S0); + end + + function [proj_geom, z_axis] = create_deform_vec_geom() + settings = struct; + settings.detectorPixelSize = 0.0748; + settings.projectionSize = [1536 1944 21]; + settings.SOD = 679.238020; %[mm] + settings.SDD = 791.365618; %[mm] + settings.voxelSize = 1; %[mm] + settings.gamma = linspace(0,300,21)*pi/180; + + S0 = zeros(settings.projectionSize(3), 12); + Sorig = [-settings.SOD, 0, 0,... % the ray origin vector (source) + (settings.SDD-settings.SOD), 0, 0,... % detector center + 0, -settings.detectorPixelSize, 0,... % detector u axis + 0, 0, settings.detectorPixelSize]; % detector v axis + + z_axis = [0, 0, 1]; + for i = 1:settings.projectionSize(3) + S0(i,:) = Sorig(:); + end + + S0(:, 1:3) = translate_3d(S0(:, 1:3), [100, 150, 0]); + S0(:, 4:6) = translate_3d(S0(:, 4:6), [100, 150, 0]); + + S0 = rotate_detector(S0, [0.48,0.32,0]); + S0 = magnify_proj(S0, 100); + for i = 1:settings.projectionSize(3) + S0(i,1:3) = rotate_around3d(S0(i,1:3), z_axis,... + settings.gamma(i)); + S0(i,4:6) = rotate_around3d(S0(i,4:6), z_axis,... + settings.gamma(i)); + S0(i,7:9) = rotate_around3d(S0(i,7:9), z_axis,... + settings.gamma(i)); + S0(i,10:12) = rotate_around3d(S0(i,10:12), z_axis,... + settings.gamma(i)); + end + + proj_geom = astra_create_proj_geom('cone_vec',... + settings.projectionSize(2), settings.projectionSize(1), S0); + end + +end diff --git a/matlab/algorithms/plot_geom/astra_create_example_fanflat.m b/matlab/algorithms/plot_geom/astra_create_example_fanflat.m new file mode 100644 index 0000000..251cb47 --- /dev/null +++ b/matlab/algorithms/plot_geom/astra_create_example_fanflat.m @@ -0,0 +1,50 @@ +function [ geom ] = create_example_fanflat( type ) +%% create_example_fanflat.m +% brief create an example geometry of type 'fanflat' +% param type type of geometry to create. provided as one of the +% following strings (not case sensitive): +% 'normal' - standard fanflat geometry +% 'vec' - example vector geometry +% return proj_geom - the geometry that was created +% date 22.06.2018 +% author Tim Elberfeld +% imec VisionLab +% University of Antwerp +% last mod 07.11.2018 +%% + if nargin < 1 + type = 'normal'; + end + + if strcmpi(type, 'normal') + geom = make_normal_geometry(); + elseif strcmpi(type, 'vec') + geom = astra_create_example_fanflat('normal'); + geom = astra_geom_2vec(geom); + + else + geom = make_normal_geometry(); + end + + function [geom] = make_normal_geometry() + % first, give measurements in mm + det_spacing = 0.035; + detector_px = 1200; + angles = linspace2(0, 2*pi, 100); + source_origin = 30; + origin_det = 200; + phantom_size = 5; + + phantom_px = 150; % voxels for the phantom + vx_size = phantom_size / phantom_px; % voxel size + + % now express all measurements in terms of the voxel size + det_spacing = det_spacing ./ vx_size; + origin_det = origin_det ./ vx_size; + source_origin = source_origin ./ vx_size; + + geom = astra_create_proj_geom('fanflat', det_spacing, ... + detector_px, angles, source_origin, origin_det); + end +end + diff --git a/matlab/algorithms/plot_geom/astra_create_example_parallel3d.m b/matlab/algorithms/plot_geom/astra_create_example_parallel3d.m new file mode 100644 index 0000000..e6fdf8e --- /dev/null +++ b/matlab/algorithms/plot_geom/astra_create_example_parallel3d.m @@ -0,0 +1,30 @@ +function [ proj_geom ] = create_example_parallel3d( type ) +%% create_example_parallel3d.m +% brief create an example geometry of type 'parallel3d' +% param type type of geometry to create. provided as one of the +% following strings (not case sensitive). +% no arguments will create a standard fanflat geometry. +% 'vec' - example vector geometry + +% return proj_geom - the geometry that was created +% date 22.06.2018 +% author Tim Elberfeld +% imec VisionLab +% University of Antwerp +% last mod 07.11.2018 +%% + if nargin < 1 + type = 'nothing'; + end + + det_spacing = [0.035, 0.035]; + detector_px = [1000, 1000]; + angles = linspace(0, 2*pi, 100); + + proj_geom = astra_create_proj_geom('parallel3d', det_spacing(1),... + det_spacing(2), detector_px(1), detector_px(2), angles); + + if strcmp(type, 'vec') + proj_geom = astra_geom_2vec(proj_geom); + end +end diff --git a/matlab/algorithms/plot_geom/private/magnify_proj.m b/matlab/algorithms/plot_geom/private/magnify_proj.m new file mode 100644 index 0000000..73e892a --- /dev/null +++ b/matlab/algorithms/plot_geom/private/magnify_proj.m @@ -0,0 +1,23 @@ + +function [ magnified_vec_geom ] = magnify_proj( vec_geom, dsdd ) +%% generate magnified vector geometry +% param type vec_geom - example vector geometry +% dsdd - deviation of SDD +% return magnified_vec_geom - the geometry that was created +% +% date 09.07.2018 +% author Van Nguyen +% imec VisionLab +% University of Antwerp +%% + magnified_vec_geom = vec_geom; + vec_sd_direction = vec_geom(:,1:3) - vec_geom(:,4:6); + norm_sd = sqrt(sum(vec_sd_direction.^2,2)); + vec_norm_sd(:,1) = norm_sd; + vec_norm_sd(:,2) = norm_sd; + vec_norm_sd(:,3) = norm_sd; + vec_sd_direction = vec_sd_direction ./ vec_norm_sd; + magnified_vec_geom(:,4:6) = vec_geom(:,4:6) - dsdd * vec_sd_direction; + clearvars -except magnified_vec_geom +end + diff --git a/matlab/algorithms/plot_geom/private/rotate_around3d.m b/matlab/algorithms/plot_geom/private/rotate_around3d.m new file mode 100644 index 0000000..e9f152b --- /dev/null +++ b/matlab/algorithms/plot_geom/private/rotate_around3d.m @@ -0,0 +1,27 @@ +function [rot_vec] = rotate_around3d(vec, ax, angle) +%% rotate_around.m +% rotate a 3d vector around an axis by an angle +% param vec: 3 x 1 vector to be rotated +% param ax: 3 x 1 vector specifying the axis +% param angle: scalar, angle to rotate by +% return: rotated vector +% +% date: 21.06.2018 +% author: someone at VisionLab, modified by Tim Elberfeld +% imec VisionLab +% University of Antwerp +%% + rot_vec = zeros(3, 1); + + rot_vec(1) = (ax(1) * ax(1) * (1-cos(angle)) + cos(angle)) * vec(1) +... + (ax(1) * ax(2) * (1-cos(angle)) - ax(3) * sin(angle)) * vec(2) +... + (ax(1) * ax(3) * (1-cos(angle)) + ax(2) * sin(angle)) * vec(3); + + rot_vec(2) = (ax(1) * ax(2) * (1-cos(angle)) + ax(3) * sin(angle)) * vec(1) +... + (ax(2) * ax(2) * (1-cos(angle)) + cos(angle)) * vec(2) +... + (ax(2) * ax(3) * (1-cos(angle)) - ax(1) * sin(angle)) * vec(3); + + rot_vec(3) = (ax(1) * ax(3) * (1-cos(angle)) - ax(2) * sin(angle)) * vec(1) +... + (ax(2) * ax(3) * (1-cos(angle)) + ax(1) * sin(angle)) * vec(2) +... + (ax(3) * ax(3) * (1-cos(angle)) + cos(angle)) * vec(3); +end diff --git a/matlab/algorithms/plot_geom/private/rotate_detector.m b/matlab/algorithms/plot_geom/private/rotate_detector.m new file mode 100644 index 0000000..019f92a --- /dev/null +++ b/matlab/algorithms/plot_geom/private/rotate_detector.m @@ -0,0 +1,19 @@ +function [ vectors_rot ] = rotate_detector(vectors, rot_angles) +%% rotate the detector part of the vectors of a vector geometry +% param vectors - vectors to transform +% param rot_angles - rotation euler angles +% +% return vectors_rot - copy of input vectors, but rotated +% +% date 09.07.2018 +% author Van Nguyen, Tim Elberfeld +% imec VisionLab +% University of Antwerp +% last mod 07.11.2018 +%% + vectors_rot = vectors; + vectors_rot(:, 1: 3) = rotate_euler3d(vectors(:, 1: 3), rot_angles); + vectors_rot(:, 4: 6) = rotate_euler3d(vectors(:, 4: 6), rot_angles); + vectors_rot(:, 7: 9) = rotate_euler3d(vectors(:, 7: 9), rot_angles); + vectors_rot(:, 10:12) = rotate_euler3d(vectors(:, 10:12), rot_angles); +end diff --git a/matlab/algorithms/plot_geom/private/rotate_euler3d.m b/matlab/algorithms/plot_geom/private/rotate_euler3d.m new file mode 100644 index 0000000..6fddb57 --- /dev/null +++ b/matlab/algorithms/plot_geom/private/rotate_euler3d.m @@ -0,0 +1,35 @@ +function [ vectors_rot ] = rotate_euler3d(vectors, rot_angles) +%% rotate some vectors by euler angles around an axis +% param vectors - vectors to transform +% param rot_angles - rotation euler angles +% +% return vectors_rot - copy of input vectors, but rotated +% +% date 09.07.2018 +% author Van Nguyen, Tim Elberfeld +% imec VisionLab +% University of Antwerp +% last mod 07.11.2018 +%% + roll = rot_angles(1); + yaw = rot_angles(2); + pitch = rot_angles(3); + roll_mat = [1 0 0; ... + 0 cos(roll) -sin(roll); ... + 0 sin(roll) cos(roll)]; + + yaw_mat = [cos(yaw) 0 -sin(yaw); ... + 0 1 0; ... + sin(yaw) 0 cos(yaw)]; + + pitch_mat = [cos(pitch) -sin(pitch) 0; ... + sin(pitch) cos(pitch) 0; ... + 0 0 1]; + % simulate rotation and translation of the DETECTOR + rot_mat = roll_mat * yaw_mat * pitch_mat; + vectors_rot = vectors; + + for i = 1:size(vectors, 1) + vectors_rot(i, :) = (rot_mat * vectors(i, :)')'; + end +end diff --git a/matlab/algorithms/plot_geom/private/translate_3d.m b/matlab/algorithms/plot_geom/private/translate_3d.m new file mode 100644 index 0000000..bd51485 --- /dev/null +++ b/matlab/algorithms/plot_geom/private/translate_3d.m @@ -0,0 +1,18 @@ +function [ transl_vectors ] = translate_3d( vectors, transl_vec) +%% translate some vectors by a translation vector +% param vectors - the vectors to translate +% param transl_vec - vector geometry translation +% return translated_vec_geom - copy of input geometry with translated +% vectors +% +% date 09.07.2018 +% author Van Nguyen, Tim Elberfeld +% imec VisionLab +% University of Antwerp +% last mod 07.11.2018 +%% + transl_vectors = vectors; % copy input + transl = repmat(transl_vec, size(vectors, 1), 1); + transl_vectors = transl_vectors + transl; +end + diff --git a/matlab/tools/astra_plot_geom.m b/matlab/tools/astra_plot_geom.m new file mode 100644 index 0000000..62eed79 --- /dev/null +++ b/matlab/tools/astra_plot_geom.m @@ -0,0 +1,93 @@ +function [] = astra_plot_geom(geometry, varargin) +%-------------------------------------------------------------------------- +% [] = astra_plot_geometry(geometry, varargin) +% +% plot an astra geometry +% +% geometry: any astra geometry, either volume geometry, projection +% geometry or an *.stl file (powered by stlRead). +% varargin: supports a variable number of (ordered and unordered) +% arguments. +% +% the remaining arguments depend on the input: +% if 'geometry' is +% - a volume geometry +% vx_size voxel size in unit of preference. Must +% be same unit that was used to scale the +% projection geometry. +% and as unorderd string-value-pairs +% Magnification magnification factor for the phantom. +% For small phantoms it might be +% necessary to scale the render up as +% otherwise it won't show up in the plot. +% Default = 1 +% LineWidth line width for the box wireframe. +% Default = 2 +% Color color of the wireframe. Default = 'r' +% +% - a projection geometry (as unordered string-value-pairs) +% RotationAxis if specified, will change the drawn +% rotation axis to provided axis. +% Must be 3-vector. Default value is +% [NaN, NaN, NaN], (meaning do not draw). +% RotationAxisOffset if specified, will translate the drawn +% rotation axis by the provided vector. +% Default = [0, 0, 0] +% VectorIdx index of the vector to visualize if +% vector geometry type. Default = 1 +% Color Color for all markers and lines if not +% otherwise specified +% DetectorMarker marker for the detector locations. +% Default = '.' +% DetectorMarkerColor color specifier for the detector marker. +% Default = 'k' +% DetectorLineColor color for the lines drawing the +% detector outline +% DetectorLineWidth line width of detector rectangle +% SourceMarker marker for the source locations +% SourceMarkerColor color specifier for the source marker +% SourceDistance (only for parallel3d and parallel3d_vec) +% distance of source to origin +% OpticalAxisColor Color for drawing the optical axis +% +% - a path to an *.stl file +% magn - magnification factor for vertices in +% CAD file. Default value = 1 +% +%-------------------------------------------------------------------------- +%-------------------------------------------------------------------------- +% This file is part of the ASTRA Toolbox +% +% Copyright: 2010-2018, imec Vision Lab, University of Antwerp +% 2014-2018, CWI, Amsterdam +% License: Open Source under GPLv3 +% Contact: astra@astra-toolbox.com +% Website: http://www.astra-toolbox.com/ +%-------------------------------------------------------------------------- + if exist('astra_create_example_cone') ~= 2 + error('Please add astra/algorithms/plot_geom to your path to use this function') + end + + if is_vol_geom(geometry) + draw.draw_vol_geom(geometry, varargin{:}); + elseif is_proj_geom(geometry) + draw.draw_proj_geom(geometry, varargin{:}); + elseif ischar(geometry) % assume 'geometry' is a path to a CAD file + draw.draw_cad_phantom(geometry, varargin{:}); + end + + % ---- helper functions ---- + function [ res ] = is_vol_geom(geom) + res = false; + if sum(isfield(geom, {'GridRowCount', 'GridColCount'})) == 2 + res = true; + end + end + + function [ res ] = is_proj_geom(geom) + res = false; + if isfield(geom, 'type') + res = true; + end + end +end |