Discussion :

[laritale15]

Bonjour,

Débutante en python, je souhaiterais avoir un peu d'aide sur le calcul+plot du NPS (Noise power spectra) à partir d'une image Dicom.

En vous remerciant par avance

[lg_53]

Si vous débuté en python, commencez par des problèmes simples.

A moins que vous ne sachiez déjà comment traiter du DICOM et du noise spectra dans un autre langage, et dans ce cas vous avez des choses à présenter pour que l'on puisse commencer à vous aider. Sinon à part vous suggérer une recherche google, on ne peut pas grand chose pour vous.

[VinsS]

Salut,

En général, pour face à ce genre de question on commence par rechercher avec son moteur de recherche préféré en utilisant py suivi de machin.

Et dans ton cas: https://pydicom.github.io/pydicom/st...g_started.html

[laritale15]

Effectivement, j'ai réussi à faire marcher un code sous Matlab et j'aimerai le retranscrire en python :

Code MATLAB :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
clear all
 
fid = fopen('I3233964.005.DCM');
img = fread(fid,[512 512],'uint16');
roi = zeros(64);
avg = zeros(64);
Ux= 0.488281;% Pixel size mm
Uy = 0.488281;% Pixel size mm
 
%% This block of code is subdividing imaage into smaller ROI and averaging purpose
for r = 1:8 
    r_shift = (r-1)*64;
    for c = 1:8 
        c_shift = (c-1)*64;
        for i = 1:64 
            for j = 1:64
                p = img(i+r_shift,j+c_shift);
                roi(i,j) = p;
            end
        end
        avg = avg+roi;
    end
end
avg = avg./64;
 
%% Actual process of NPS calculation
scale = (Ux*Uy)./(64*64);%Scaling fator as per NPS calculation formula
f_x = 1/(2*Ux);%Nyquiest frequency along x direction
f_y = 1/(2*Ux);%Nyquiest frequency along y direction
FFT_2d = (fftshift(fft2(avg))).^2;% Power spectrum calculation
NPS = abs(FFT_2d).*scale; %% 2D NPS 
f = linspace(-f_x,f_y,64);% x-axis  vector for 1D NPS 
X_img = linspace(-f_x,f_x,512);% X axis of NPS image
Y_img = linspace(-f_x,f_x,512);% Y axis of NPS image
 
figure(1)
subplot(2,2,1)
imagesc(X_img,Y_img,img)
colormap gray
xlabel('X [cm]'); ylabel('Y [cm]')
title('noise image')
subplot(2,2,2)
imagesc(f,f,log(NPS))
colormap gray
xlabel('frequency [mm^{-1}]'); ylabel('frequency [mm^{-1}]');
title('2D NPS')
subplot(2,2,3)
plot(f,NPS(:,32))
xlabel('frequency [mm^{-2}]'); ylabel('NPS [mm^{-2}]')
title('1D NPS from central slice')
%xlim([0 1])
 
subplot(2,2,4)
plot(f,mean(NPS,2))
%xlim([0 1])
 
xlabel('frequency [mm^{-2}]'); ylabel('NPS [mm^2]')
title('1D NPS along X direction')

Je fais quelques essais, mais déjà ce petit bout de code ne fonctionne pas. Je n'ai pas de message d'erreur mais juste un plot vide

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
from scipy import fftpack, ndimage
import matplotlib.pyplot as plt
import SimpleITK as sitk
 
image = sitk.ReadImage('./I3233964.013.DCM')
arr_im1 = sitk.GetArrayFromImage(image)
fft2 = fftpack.fft2(arr_im1[:,:,1])
 
 
plt.figure(figsize=(15,8), dpi=80)
subplot2=plt.subplot(1,1,1)
plt.plot(abs(fft2))

[lg_53]

Vous devriiez ouvrir un tuto Python et Matplotlib, car ce n'est pas du tout comme cela que l'on affiche des images avec matplotlib. La commande plot est fait pour tracer une courbe...

Vous devriiez également vous plongez dans la doc de SimpleITK ainsi que dans des exemples d'utilisations, car il y a déjà plein de chose de faites, comme celle là :
http://simpleitk.github.io/SimpleITK...7s_imshow.html