2D matrix resize or interpolation

Hello,

I'm new to matlab and I'm not sure what would be the best way to do the following.

I have an image/data matrix, I, of size (rows,cols) = (y,x) = [384,512].

The spatial resolution in the x direction is 0.078021 cm (total span is ~40cm)

The spatial resolution in the y direction is 0.078021 (total span is ~30cm)

Q1.

If I would like to have a finer sampling of the data by doubling the samples (eg Inew = [768, 1024], which function would you suggest?

Q2.

If I had another separate image but with offset x,y values, how would I interpolate this data onto the same grid as the other matrix?

Thank you for your help,

RR

ps. I found some info on the interp2 function? not sure if this is the best function to use?

ps. I have attached a zip file with

• dicom image (of a radiotherapy dose on a single plane - my image)
• excel workbook with the data extracted
• an m script with what I tried to do but I got an error

Error using griddedInterpolant

The grid must be created from grid vectors which are strictly monotonically increasing.

Error in interp2>makegriddedinterp (line 228)

F = griddedInterpolant(varargin{:});

Error in interp2 (line 136)

F = makegriddedinterp(X, Y, V, method,extrap);

Error in LoadDosePlaneBailey2D (line 70)

Id_cgy_q = interp2(X,Y,Id_cgy,Xq,Yq);

Answer 1

[Cris LaPierre]

Q1: B = imresize(A,[numrows numcols])

Q2: What type of images are you working with? Please share an example. You can attach images to your post using the paperclip icon.

[Russell]

Thank you for your help.

I have added the files - not sure if you can see them now?

There will be 2 images to compare

• radiotherapy dose calculated on a plane
• eventually I will need to compare to an x-ray image (portal image) which results in my Q2 (the imager will not exactly be centered - therefore there will be an offset from center in each direction x,y

[Cris LaPierre]

It looks like the pixel info and spreadsheet data are the same.

unzip('SampleFiles.zip')

% Use medical imaging toolbox
MI = medicalImage('ExampleDicomFileDosePlane.dcm')

MI = 
  medicalImage with properties:

          Pixels: [384x512 uint32]
        Colormap: []
    SpatialUnits: "mm"
       FrameTime: []
       NumFrames: 1
    PixelSpacing: [0.7802 0.7793]
        Modality: 'RTDOSE'
    WindowCenter: []
     WindowWidth: []

img = extractFrame(MI,1);
imshow(img,[])

% Resize
B1 = imresize(img,[798,1024]);
figure
imshow(B1,[])

Neither of these approaches take into consideration the actual pixel spacing. Instead, row/col index number is used. To display an image where the X and Y spacing is different, use the XData and YData inputs to imshow.

figure
ext = (size(img)-1).*MI.PixelSpacing

ext = 1×2

298.8198 398.2207

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

imshow(B1,[],"XData",[0 ext(2)],"YData",[0 ext(1)])

It's subtle, but you can compare the difference using imshowpair

RA = imref2d(size(img));
RB = imref2d(size(img),MI.PixelSpacing(2),MI.PixelSpacing(1));
imshowpair(img,RA,img,RB,'montage')

You could just as easily use the spreadsheet values. However, you lose the metadata information, which captures the pixel spacing.

dose = readmatrix('ExampleTxtDataOfDosePlane.xlsx');
imshow(dose,[])

B2 = imresize(dose,[768, 1024]);
figure
imshow(B2,[],"XData",[0 ext(2)],"YData",[0 ext(1)])

[Russell]

Thanks for the help, I will try to look through it to see if i can understand it.

[Russell]

how do i use the function 'medicalimage'? I was using a function called dicomread but when I tried medicalImage I got an undefined function

[Cris LaPierre]

medicalImage is part of the Medical Imaging Toolbox. It sounds like you haven't installed it. If it is included on your license, you can install it using the Add-Ons Explorer.

[Russell]

Hi Cris,

I guess I'm still having problems understanding how to do this. I added another pdf showing what I'm trying to do. I'm unable to upload the calculated dicom planar dose image. I was wondering if what I'm asking is possible with what you showed me earlier? (this will be compared to the actual dose image measured that I showed earlier)

Thanks for your help.

[Cris LaPierre]

Try zipping your second image and attaching it.

[Russell]

Hi again Cris,

I have some sample data.

Image1 is a 43cmx43cm image

I would like to extract a 30cmx40cm (inner portion)

I have positions (x2,y2) that I would like to sample Image1 to get my 30cmx40cm image2.

But I don't seem to get a good result.

Would you be able to help figure out what I'm doing wrong?

I zipped and attached my files.

Thanks so much for any help you can provide

[Cris LaPierre]

You talk about images, but are sharing text files. Do you have actual images?

[Russell]

The text files contain the data of the image - the imager gives the results in a txt file I just kept the data and removed all the header txt

• the 2d array are the values at the corresponding position in the arrays of x,y
• with the original array spanning 43cm in x and y

But I need a subset of x,y only spanning the inner space spanning 40cmx30cm and I wanted to sample it in 1024pixels by 768pixels

[Cris LaPierre]

If your images are dicoms, there is metadata that might make this process easier.

[Russell]

Hi Cris, thanks. The images are from a detector for MV xrays. I couldn't include the files because they were too big. Here is a link to a dropbox folder with the files. Thanks for your help.

There are 2 prediction images and 1 meausured image. Eventually I will have another prediciton image in dicom format that will need to be converted to the company's format so I can review it in their software.

https://www.dropbox.com/scl/fo/xt3b1srt8mppbu3qjm8pz/AEMLzSK1DvBDLAflhzLr0GI?rlkey=o7ct8wahhpurln1jw3d7evc6f&st=e3eay333&dl=0

[Cris LaPierre]

Just confirming that you are working with dxf files, not dicoms, correct?

I'm confused because the most common dxf file is an AutoCAD drawing file.

[Russell]

Yes, sorry I forgot to mention that, yes they use the same extension as Autocad. But these are image files from the MV x-ray imager.

[Cris LaPierre]

I'm defnitely not an expert in this space, so I'm still confused. Based on the metadata in the dsf files

• Image_measured_30x40.dxf is actually 21.5x28.7
• Image_prediction1_43x43.dxf is actually 31.1x31.1
• Image_prediction2_30x40.dxf is actually 21.4x28.5

It looks like you have a 2-step process - first to extract the desired region using indexing, then resizing. Given that, I think the best approachis to use imresize on the extracted iage. Pick the approprite method for your task.

[Cris LaPierre]

Heres a function that will read in your dxf files. Not necessarily the best code, but does the job.

unzip Image_prediction1_43x43.zip

params = readDXF('Image_prediction1_43x43.dxf');
figure
imagesc(params.Pixels)

Xsz = params.Size1*params.Res1 

Xsz = 311.4286

Ysz = params.Size2*params.Res2 

Ysz = 311.4286

function data = readDXF(fname)
data = struct;
fid = fopen(fname);
while ~feof(fid)
    ln = fgetl(fid);
    [param,val] = strtok(ln,"=");
    switch param
        case {'[General]','[Geometry]','[Interpretation]','[Patient]','[Field]','[PortalDose]'}
            continue
        case '[Data]'
            fseek(fid,0,'eof');
            fgetl(fid);
        otherwise
            val = strtok(val,'=');
            pat = asManyOfPattern(characterListPattern("0123456789-+."));
            nums = extract(val,pat);
            if length(val)==length(nums{1})
                data.(param) = str2double(val);
            else
                data.(param) = val;
            end
    end
end
fclose(fid);
val = readmatrix(fname,'FileType','text','NumHeaderLines',48);
data.Pixels = val;
end

Answer 2

[Les Beckham]

For Q1, I would suggest imresize. Can you clarify what you mean by "offset x,y values" in Q2?

[Russell]

Hello, I added some information as requested by Cris. Not sure if this helped?

The first image will be a dose plane

X size = ~40cm with 512pix samples at a res of ~0.078cm

Y size = ~30cm with 384pix samples at res of ~0.078cm

The second image will be a x-ray image but the detector may not exactly be centered (this is where the offset comes from)

xoffset, yoffset

I would like to resample either the dose plane calc on the same positions as the imager or vice versa.

Does this help clarify?

Answer 3

[Walter Roberson]

If I would like to have a finer sampling of the data by doubling the samples (eg Inew = [768, 1024], which function would you suggest?

If you were willing to have Inew = [767, 1023] instead of [768, 1024] then you can average adjacent rows and adjacent columns and manually insert the averaged values into the proper place. To be honest, though, calling interp2() is a heck of a lot easier.

The reason the output would be Inew = [767, 1023] instead of [768, 1024] is that you would be calculating new points between each existing point. If you had 3 points across then you have 2 intermediate points, for a total of 5. Calculating intermediate points gets you output of size 2*N-1 not output of size 2*N

Getting output exactly twice the original size requires interpolation at "just less than half" apart, like positions 1, 1.49, 1.98, 2.47,...

[Russell]

Would i be able to do something like this?

[Russell]

Using interp2 on the finer resolution grid?

[Image Analyst]

If you want the same pixel values in the two pixels, then use the 'nearest' option of imresize. If you don't specify then I believe it uses bicubic interpolation.

img = [1,2,3,4;5,6,7,8]

img = 2×4

1 2 3 4 5 6 7 8

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

img2 = imresize(img, 2, 'nearest')

img2 = 4×8

1 1 2 2 3 3 4 4 1 1 2 2 3 3 4 4 5 5 6 6 7 7 8 8 5 5 6 6 7 7 8 8

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

img3 = imresize(img, 2)

img3 = 4×8

0.5312 0.8047 1.3516 1.8750 2.3750 2.8984 3.4453 3.7188 1.7188 1.9922 2.5391 3.0625 3.5625 4.0859 4.6328 4.9062 4.0938 4.3672 4.9141 5.4375 5.9375 6.4609 7.0078 7.2812 5.2812 5.5547 6.1016 6.6250 7.1250 7.6484 8.1953 8.4688

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

●

Note how it's handling the edge effects -- essentially the pixel value in your top row (1 in my example) is considered in the middle of the upsampled version (your row 2) so that's why in the bicubic version the leftmost is not exactly 1.

[Russell]

Thank you for the help, I'll take a look through this.

[Walter Roberson]

img = uint8([1,2,3,4;5,6,7,8])

img = 2×4

1 2 3 4 5 6 7 8

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

img3 = imresize(img, 2)

img3 = 4×8

1 1 2 2 3 3 4 4 2 2 3 3 4 4 5 5 4 4 5 5 6 6 7 7 5 5 6 6 7 7 8 8

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

img4 = cast(imresize(double(img),2),'uint8')

img4 = 4×8

1 1 1 2 2 3 3 4 2 2 3 3 4 4 5 5 4 4 5 5 6 6 7 7 5 6 6 7 7 8 8 8

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[Russell]

Hello again,

Thanks, I think so far this helps

I have now my matrix with the right size (768,1024) - doubled in each direction (or halved in spatial resolution)

If I know the corresponding grid positions x[i:1-768] and y[j=1:1024] in mm.

Should I be able to use interp2 to reinterpolate the matrix/image on the grid positions i need? (xq, yq)?

Img_interp = interp2(x,y,Img,xq,yq)?

[Walter Roberson]

Yes. However, to double the resolution, xq would need to be

linspace(x(1), x(end), 2*(x(end)-x(1)+1))

which would be

linspace(1, 768, 2*768)

This will not be [1, 1.5, 2, 2.5, 3, 3.5, ...] -- instead it will be [[1, 2302/1535, 3069/1535, 3836/1535, 4603/1535 ...] or roughly [1, 1.49967426710098, 1.99934853420195, 2.49902280130293, 2.99869706840391, ...]

[Russell]

Thanks. I guess I'm a little confusd over the spatial vector? I have something like this now

so I'm not sure how to get the interpolated image Iq from I, x[i], y[j] if I have another set of positions xq[i], yq[j] that I want to sample?

[Russell]

Hi Walter,

I guess I'm still having problems understanding how to do this. I added another pdf showing what I'm trying to do. I'm unable to upload the calculated dicom planar dose image. I was wondering if what I'm asking is possible with what you showed me earlier? (this will be compared to the actual dose image measured that I showed earlier)

Thanks for your help.