Category Archives: Matlab

Save Image into PPM ASCII Format in Matlab and Python

1 Reply

Here I wrote a function to save image/matrix into PPM format. I tried to simplify the saving process by using the matlab function ‘dlmwrite’ such that I only need to write the header for this PPM format. However, before streaming the image data, different channels of the input image must be interleaved. I experimented two ways to do this. Both work just fine, and the use of matlab function ‘cat’ (or ‘permute’) make the implementation very simple.

function write_ppm(im, fname, xval)

[height, width, c] = size(im);
assert(c == 3);

fid = fopen( fname, 'w' );

%% write headers
fprintf( fid, 'P3\n' );
fprintf( fid, '%d %d\n', width, height);
fprintf( fid, '%d \n', xval); %maximum values

fclose( fid ); 

%% interleave image channels before streaming
c1 = im(:, :, 1)';
c2 = im(:, :, 2)';
c3 = im(:, :, 3)';
im1 = cat(2, c1(:), c2(:), c3(:));

%% data streaming, could be slow if the image is large
dlmwrite(fname, int32(im1), '-append', 'delimiter', '\n')

The implementation in Python is even simpler. ‘f.write()’ is capable of writing the entire image in one function call with the help of ‘\n’.join(). 

def write_ppm(fname, im, xval):
    
    height, width, nc = im.shape
    assert nc == 3
    
    f = open(fname, 'w')
    
    f.write('P3\n')
    f.write(str(width)+' '+str(height)+'\n')
    f.write(str(xval)+'\n')
    
    # interleave image channels before streaming    
    c1 = np.reshape(im[:, :, 0], (width*height, 1))
    c2 = np.reshape(im[:, :, 1], (width*height, 1))
    c3 = np.reshape(im[:, :, 2], (width*height, 1))
    
    im1 = np.hstack([c1, c2, c3])
    im2 = im1.reshape(width*height*3)

    f.write('\n'.join(im2.astype('str')))
    f.write('\n')

    f.close()

write DICOM images from a 3D array

1 Reply

To write a DICOM image sequence from a 3D data array using MATLAB, you will need to specify ‘PixelSpacing’ and the ‘SliceThickness’, otherwise the 3D reconstruction software would have no idea to determine the length of the volume data. One way is to copy a sample metadata structure from another volume, update it with the new  ‘PixelSpacing’ and ‘SliceThickness’,  and pass it to dicomwrite. dicomwrite can write the correct ‘PixelSpacing’ and ‘SliceThickness’ into the dicom file, however, the reconstruction software (for example, slicer) would still treat ‘SliceThickness’ as 1 even it displays the correct ‘SliceThickness’ value. Another way is to directly specify these matadata values as the following. You just need to find the corresponding SOP Class UID for your dicom data from http://www.apteryx.fr/dicom/dicom_conformance


% T defines the 'PixelSpacing' and the 'SliceThickness'

%% save dicom images
for i = 1 : n_slices

    status = dicomwrite(uint16(data(:,:,i)), ...
    file_name, ...
    'SliceThickness', T(3), ...
    'PixelSpacing', [T(1); T(2)], ...
    'SliceLocation', location_base - (i-1)*T(3), ...
    'ImageOrientationPatient', [1 0 0 0 1 0], ...
    'CreateMode', 'Copy', ...
    'SOPClassUID', '1.2.840.10008.5.1.4.1.1.2');

end

Show the input data for the accumulating function of Matlab accumarray

1 Reply

When using MATLAB accumarray, if the subscripts in 'subs' are not sorted, the accumulating function ‘fun‘ should not depend on the order of the values in its input data. Otherwise you would have weird results. To show the input data for the accumulating function, one can define the function as ‘@(x){x}’. This would save the input data for you. For example,

>> val = 101:105;

>> subs = [1; 2; 4; 2; 4]

>> A = accumarray(subs, val, [], @(x){x})


A =


[       101]

[2x1 double]

[]

[2x1 double]

tsearchn and pointLocation

2 Replies

To use pointLocation() to find the enclosing simplex (triangle/tetrahedron) for each query point location in Matlab, one has to supply this function with a ‘DelaunayTri’ structure as an input. This is a newly introduced class for performing delaunay triangulation. While in my case I need to reuse my saved delaunay triangulation results and assigned it to the Triangulation component of the ‘DelaunayTri’ structure. However I got this error,

??? Error using ==> subsasgn
Cannot assign values to the Triangulation.

I think this is because Triangulation is a private member of the class DelaunayTri, and every time one revises the values of X, this Triangulation  is updated automatically. If you are still allowed, you can use the old function ‘tsearchn’ to find the enclosing simplex of the delaunay triangulation. It is also very easy to write a simple function to compute the barycentric coordinate for each query point location. In case ‘tsearchn’ is not supported by the future Matlab, you can also write your own version of ‘tsearchn’. The brute force way is to check the barycentric coordinate of that point for each triangle. If one of the coordinates is negative, then the point is outside of that triangle.

draw line segments in matlab

3 Replies

After spending some effort searching online, I found one solution:

r = 5;
x1 = linspace(pos.cx – r, pos.cx + r);
y1 = linspace(pos.cy, pos.cy);

x2 = linspace(pos.cx, pos.cx);
y2 = linspace(pos.cy – r, pos.cy + r);

index1 = sub2ind(size(im), round(y1), round(x1));
index2 = sub2ind(size(im), round(y2), round(x2));
im(index1) = 127; % set the valudes to 127 for display
im(index2) = 127; % set the valudes to 127 for display

The above code snippet is used to draw a cross on an image or matrix. Because of over sampling, the line drawn is without aliasing artifacts. For drawing lines on figures, there exists an easier solution, just use function ‘line’. For example:

line(rect(1,:), rect(2,:),’Color’,’r’,’LineWidth’,2);