Navigation

Inscrivez-vous gratuitement
pour pouvoir participer, suivre les réponses en temps réel, voter pour les messages, poser vos propres questions et recevoir la newsletter

Calcul scientifique Python Discussion :

Calcul NPS (Noise power spectra)


Sujet :

Calcul scientifique Python

  1. #1
    Nouveau Candidat au Club
    Calcul NPS (Noise power spectra)
    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

  2. #2
    Membre chevronné
    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.

  3. #3
    Expert éminent
    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

  4. #4
    Nouveau Candidat au Club
    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))

  5. #5
    Membre chevronné
    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-Notebooks/10_matplotlib%27s_imshow.html