How can I crop and draw a point in a set of dicom images.

I am trying to crop 200 of a serie of 1200 dicom images and draw a dot, square, etc in 20-30 slices of it, which should be seen in a 3D model in the Slicer Software.

The problem is that first, cropped images are not treated as equal as the original in Slicer, because the software thinks that each one of them comes from a different patient. This issue doesn't allow us to build the cropped 3D model with the dot on it.

Also, I want to avoid drawing the point in each of the slices which I want to mark.

The script is the following:

    for p=520:725
        I = dicomread(['IM-0001-0' num2str(p) '.dcm']);
        IC= I(10:100,105:270);
        if (p == 644)
            info = dicominfo(IC);
            dl = info.PixelSpacing; %dl = [dx, dy];
            figure;
            imshow(IC);
            x = ginput(1);
            IC(x(2)-5:x(2)+5,x(1)-5:x(1)+5) = 1; 
            imshow(IC)
        end
        if any(p == 645:665)
            info = dicominfo(IC);
            dl = info.PixelSpacing; 
            IC(x(2)-5:x(2)+5,x(1)-5:x(1)+5) = 1;  
        end
        dicomwrite(IC,['IMcr-0001-0' num2str(p) '.dcm']);
    end

Thanks.

Answer 1

[Walter Roberson]

dicominfo() applied to a numeric array makes no sense.

You can transfer the information about patient and so on by including the info structure in the dicomwrite

    for p=520:725
        dfile_in = sprintf('IM-0001-%04d.dcm', p);
        dfile_out = sprintf('IMcr-0001-%04d.dcm', p);

        info = dicominfo(dfile_in);
        % dl = info.PixelSpacing; %dl = [dx, dy];  %you do not use this
        I = dicomread(info);

        IC= I(10:100,105:270);
        if (p == 644)
            figure;
            imshow(IC);
            x = ginput(1);
            IC(x(2)-5:x(2)+5,x(1)-5:x(1)+5) = 1; 
            imshow(IC)
        end
        if ismember(p, 645:665)
            dl = info.PixelSpacing; 
            IC(x(2)-5:x(2)+5,x(1)-5:x(1)+5) = 1;  
        end
        dicomwrite(IC, dfile_out, info);
    end