im2patch

Below is a demonstration of the features of the im2patch function

Contents

clear; close all; clc;

Syntax

[F,V,C]=im2patch(M,IND,ptype);

Description

This function generates patch data for 3D images. The patch data is only generated for the voxels specified by the indices (logic or linear indices). The patch data is created according to the patch type desired (e.g. voxel type of slice type).

Examples

Plot settings

cMap=parula(250);
faceAlpha1=1;
faceAlpha2=0.5;
edgeColor1='none';
edgeColor2='none';
fontSize=15;

Example: Introduction to using im2patch for voxel plotting

The voxel for which patch data is to be specified can be defined by supplying linear indices or a logic array.

% Simulating an image
M=rand(3,5,7);
% Example supplying linear indices, here all voxels
indPatch=1:numel(M);

Creating patch data for voxel display

[F,V,C]=im2patch(M,indPatch,'v');

cFigure;
title('patch type: v');
xlabel('J - columns');ylabel('I - rows'); zlabel('K - slices'); hold on;
gpatch(F,V,C,'k',faceAlpha2);
colormap(cMap); caxis([min(M(:)) max(M(:))]); colorbar;
axisGeom(gca,fontSize); camlight headlight;
drawnow;

Study the size of the face and vertex arrays to confirm that the patch type 'v' creates vertices and faces for each voxel. Shared vertices and faces are not removed. This option requires the most memory.

disp(num2str(size(F)));
disp(num2str(size(V)));

630    4
192    3

Creating patch data for voxel display with shared vertices and faces removed

[F,V,C]=im2patch(M,indPatch,'vu');

cFigure;
title('patch type: vu');
xlabel('J - columns');ylabel('I - rows'); zlabel('K - slices'); hold on;
gpatch(F,V,C,'k',faceAlpha2);
colormap(cMap); caxis([min(M(:)) max(M(:))]); colorbar;
axisGeom(gca,fontSize); camlight headlight;
drawnow;

Study the size of the face and vertex arrays to confirm that the patch type 'vu' ensures that shared faces and vertices are not shared. Each vertex and face is therefore unique. This saves memory with respect to using the 'v' patch type for voxel display. The figures look identical except that shared faces may appear less dark when transparency is on since now one one face is used. Color information is shared too (averaged).

disp(num2str(size(F)));
disp(num2str(size(V)));

386    4
192    3

Creating patch data for voxel display with only non-shared faces and vertices

[F,V,C]=im2patch(M,indPatch,'vb');

cFigure;
title('patch type: vb');
xlabel('J - columns');ylabel('I - rows'); zlabel('K - slices'); hold on;
gpatch(F,V,C,'k',faceAlpha2);
colormap(cMap); caxis([min(M(:)) max(M(:))]); colorbar;
axisGeom(gca,fontSize); camlight headlight;
drawnow;

Study the size of the face and vertex arrays to confirm that the patch type 'vb' helps to plot only non-shared vertices and faces. In the case of an enclosed shape filled with voxels thus only the boundary faces are displayed. This path type appears no different than the other is transparency is not on however it is much lighter on memory than the above path types.

disp(num2str(size(F)));
disp(num2str(size(V)));

142    4
144    3

Example: Introduction to using im2patch for slice plotting

cFigure;
title('patch type: si, sj, sk');
xlabel('J - columns');ylabel('I - rows'); zlabel('K - slices'); hold on;

%Setting up indices for I direction slices
S=round(size(M,1)./2); %Selection of middle slice
L_plot=false(size(M)); L_plot(S,:,:)=1;
indPatch=find(L_plot);
[F,V,C]=im2patch(M,indPatch,'si'); %Creating patch data for y mid-voxel slices
gpatch(F,V,C);

%Setting up indices for J direction slices
S=round(size(M,2)./2); %Selection of middle slice
L_plot=false(size(M)); L_plot(:,S,:)=1;
indPatch=find(L_plot);
[F,V,C]=im2patch(M,indPatch,'sj'); %Creating patch data for x mid-voxel slices
gpatch(F,V,C);

%Setting up indices for Z direction slices
S=round(size(M,3)./2); %Selection of middle slice
L_plot=false(size(M)); L_plot(:,:,S)=1;
indPatch=find(L_plot);
[F,V,C]=im2patch(M,indPatch,'sk'); %Creating patch data for z mid-voxel slices
gpatch(F,V,C);

colormap(cMap); caxis([min(M(:)) max(M(:))]); colorbar;
axisGeom(gca,fontSize); camlight headlight;

drawnow;

The path type s*u are simular to s* but use shared vertices

Example: Comparison to standard MATLAB imagesc function and the patch type MATLAB functions slice and pcolor

The comparison is what motivates the choice of coordinate system for im2patch i.e. that it is meant to aid in the visualization of image data as is expected of image data.

cFigure;
subplot(2,2,1);
title('MATLAB imagesc function');
xlabel('J - columns');ylabel('I - rows'); zlabel('K - slices'); hold on;
imagesc(M(:,:,S));
colormap(cMap); caxis([min(M(:)) max(M(:))]);
axisGeom(gca,fontSize); view(2);

subplot(2,2,2);
title('MATLAB slice function');
xlabel('J - columns');ylabel('I - rows'); zlabel('K - slices'); hold on;
slice(M,[],[],S);
axisGeom(gca,fontSize);
colormap(cMap); caxis([min(M(:)) max(M(:))]);
axisGeom(gca,fontSize); view(2);

subplot(2,2,3);
title('im2patch function');
xlabel('J - columns');ylabel('I - rows'); zlabel('K - slices'); hold on;
gpatch(F,V,C);
colormap(cMap); caxis([min(M(:)) max(M(:))]);
axisGeom(gca,fontSize); view(2);

subplot(2,2,4);
title('MATLAB pcolor function');
xlabel('J - columns');ylabel('I - rows'); zlabel('K - slices'); hold on;
pcolor(M(:,:,S));
colormap(cMap); caxis([min(M(:)) max(M(:))]);
axisGeom(gca,fontSize); view(2);

drawnow;

Note that the image size is wrong for the slice and pcolor commands. This is because intensities appear to be defined on voxels vertices (and are reinterpolated onto faces) for these functions instead of voxels centres as should be the case for image data.

Example: Creating and plotting combined voxel and slice patch data

Simulating 3D image

[X,Y,Z]=meshgrid(linspace(-4.77,4.77,25));
phi=(1+sqrt(5))/2;
M=1/6*(2 - (cos(X + phi*Y) + cos(X - phi*Y) + cos(Y + phi*Z) + cos(Y - phi*Z) + cos(Z - phi*X) + cos(Z + phi*X)));

Creating and plotting patch data. Last example illustrates the use of a specific logic description (i.e. a mask) for the voxels of interest

cFigure;
title('Combined voxel and slice plotting');
xlabel('J - columns');ylabel('I - rows'); zlabel('K - slices'); hold on;

% Setting up indices for I direction slices
S=round(size(M,1)./2); %Selection of middle slice
L_plot=false(size(M)); L_plot(S,:,:)=1;
indPatch=find(L_plot);
[F,V,C]=im2patch(M,indPatch,'si'); %Creating patch data for y mid-voxel slices
gpatch(F,V,C,'none',faceAlpha2);

% Setting up indices for J direction slices
S=round(size(M,2)./2); %Selection of middle slice
L_plot=false(size(M)); L_plot(:,S,:)=1;
indPatch=find(L_plot);
[F,V,C]=im2patch(M,indPatch,'sj'); %Creating patch data for x mid-voxel slices
gpatch(F,V,C,'none',faceAlpha2);

% Setting up indices for K direction slices
S=round(size(M,3)./2); %Selection of middle slice
L_plot=false(size(M)); L_plot(:,:,S)=1;
indPatch=find(L_plot);
[F,V,C]=im2patch(M,indPatch,'sk'); %Creating patch data for z mid-voxel slices
gpatch(F,V,C,'none',faceAlpha2);

% Setting up indices for voxels to plot
L_mask=M>-0.2 & M<0;
[F,V,C]=im2patch(M,L_mask,'vb'); %Creating patch data for selection of high voxels

gpatch(F,V,C);
colormap(cMap); colorbar; caxis([min(M(:)) max(M(:))]);
axisGeom(gca,fontSize); camlight headlight;
drawnow;

Example: Shrinking patch data through combination with scalePatch function

% Using the function |scalePatch| patch data can be shrunk to aid
% visualisation (can be a tool to avoid memory costly transparency for
% instance).

[Ev,Vv,Cv]=im2patch(M,L_mask,'h'); %This creates a hexahedral element for each voxel
[Evs,Vvs]=scalePatch(Ev,Vv,0.5); %Apply voxel element scaling
[Fvs,Cvs]=element2patch(Evs,Cv); %Convert to quad faces for plotting

Plotting the voxels

cFigure;
title('Scaled (shrunk) patch data');
xlabel('J - columns');ylabel('I - rows'); zlabel('K - slices'); hold on;
gpatch(Fvs,Vvs,Cvs);
colormap(cMap); colorbar; caxis([min(Cvs(:)) max(Cvs(:))]);
axisGeom(gca,fontSize); camlight headlight;
drawnow;

Example: Medical image data and coordinate manipulation due to voxel size

% Get a 3D image
load mri;
M=squeeze(D); %example image data set
v=2./[1,1,.4]; %example voxel size

The voxels to display can be specified as a list (vector) of voxels numbers (linear indices) or using a mask (logic array).

%Defining row, column and slice indicices for slice patching
sliceIndexI=round(size(M,1)/2); %(close to) middle row
sliceIndexJ=round(size(M,2)/2); %(close to) middle column
sliceIndexK=round(size(M,3)/2); %(close to) middle slice

%Defining "masks" i.e. logic arrays with ones for voxels of interest
logicSliceI=false(size(M));
logicSliceI(sliceIndexI,:,:)=1;
logicSliceI=logicSliceI & M>0;

logicSliceJ=false(size(M));
logicSliceJ(:,sliceIndexJ,:)=1;
logicSliceJ=logicSliceJ & M>0;

logicSliceK=false(size(M));
logicSliceK(:,:,sliceIndexK)=1;
logicSliceK=logicSliceK & M>0;

%Defining voxel indices for voxels of interest
T_low=min(M(:))+((max(M(:))-min(M(:)))/10); %Threshold example
logicVoxels=(M>T_low);
logicVoxels(:,1:sliceIndexJ,:)=0;

Creating patch data The patch data consists of a matrix array defining the faces, a matrix array defining the vertices and a vector for the colour data. The vertices are based on the image coordinates however they are formatted as: [X(:) Y(:) Z(:)]. X relates to columns, Y to rows and Z to slices. Use a function like im2cart , or im2mrcart to convert image to cartesian coordinates.

[Fv,Vv,Cv]=im2patch(M,logicVoxels,'vb');
[Fx,Vx,Cx]=im2patch(M,logicSliceJ,'sj');
[Fy,Vy,Cy]=im2patch(M,logicSliceI,'si');
[Fz,Vz,Cz]=im2patch(M,logicSliceK,'sk');

% Convert image coordinates to cartesian coordinates
[Vv(:,1),Vv(:,2),Vv(:,3)]=im2cart(Vv(:,2),Vv(:,1),Vv(:,3),v);
[Vx(:,1),Vx(:,2),Vx(:,3)]=im2cart(Vx(:,2),Vx(:,1),Vx(:,3),v);
[Vy(:,1),Vy(:,2),Vy(:,3)]=im2cart(Vy(:,2),Vy(:,1),Vy(:,3),v);
[Vz(:,1),Vz(:,2),Vz(:,3)]=im2cart(Vz(:,2),Vz(:,1),Vz(:,3),v);

cFigure;
title('MRI visualisation, slices and voxels in cartesian coordinates with aid of voxel size');
xlabel('X (mm)');ylabel('Y (mm)'); zlabel('Z (mm)'); hold on;
gpatch(Fv,Vv,Cv,edgeColor1,faceAlpha2);
gpatch(Fx,Vx,Cx,edgeColor1,faceAlpha1);
gpatch(Fy,Vy,Cy,edgeColor1,faceAlpha1);
gpatch(Fz,Vz,Cz,edgeColor1,faceAlpha1);
colormap(gray(250)); colorbar;
axisGeom(gca,fontSize); camlight headlight;
drawnow;

Example: Combining colormap and RGB driven patch colours

N.B. The figure renderer might have to be set to OPENGL.

Convert voxels colouring to RGB type, here a simple conversion to a gray scale description is used

Cv=(Cv*ones(1,3))./max(Cv(:));

Plotting the voxels

cFigure;
title('MRI visualisation, slices and voxels, colormap and RGB driven respectively');
xlabel('X (mm)');ylabel('Y (mm)'); zlabel('Z (mm)'); hold on;

%RGB driven
gpatch(Fv,Vv,Cv,edgeColor1,faceAlpha2);

%Colormap driven
gpatch(Fx,Vx,Cx,edgeColor1,faceAlpha1);
gpatch(Fy,Vy,Cy,edgeColor1,faceAlpha1);
gpatch(Fz,Vz,Cz,edgeColor1,faceAlpha1);

colormap(cMap); colorbar;
axisGeom(gca,fontSize); camlight headlight;
drawnow;

GIBBON www.gibboncode.org

Kevin Mattheus Moerman, [email protected]

GIBBON footer text

License: https://github.com/gibbonCode/GIBBON/blob/master/LICENSE

GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for image segmentation, image-based modeling, meshing, and finite element analysis.

Copyright (C) 2006-2020 Kevin Mattheus Moerman

This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.

You should have received a copy of the GNU General Public License along with this program. If not, see http://www.gnu.org/licenses/.

Published with MATLAB® R2020a