Calculation of spatial frequency

Algorithms
1

Hello ITK,
i want to compute “spatial frequency” for wavelet. In these publication (DOI: 10.17694/bajece.24838,H. Boztoprak ) they define it as:
afbeelding
with M, N the voxel indice and F his value.
So i create a code like this which is suboptimal.

def SpatialFrequency(image):
    SizeMatrix=image.GetSize()
    Square_diff_x=0
    Square_diff_y=0
    Square_diff_z=0
    Nvoxel=0
    for x in range(SizeMatrix[0]-1):
        for y in range(SizeMatrix[1]-1):
            for z in range(SizeMatrix[2]-1):
                    Square_diff_x=Square_diff_x+(image.GetPixel(x+1,y,z)-image.GetPixel(x,y,z))**2
                    Square_diff_y=Square_diff_y+(image.GetPixel(x,y+1,z)-image.GetPixel(x,y,z))**2
                    Square_diff_z=Square_diff_z+(image.GetPixel(x,y,z+1)-image.GetPixel(x,y,z))**2
                    Nvoxel=Nvoxel+1
    SF=(Square_diff_x/(Nvoxel)+Square_diff_y/(Nvoxel)+Square_diff_z/(Nvoxel))**0.5
    return SF

can you help me to improve it?
Thank you,
Cyril Jaudet

2

this new code is really faster but not give the same result. Any clue?

def SpatialFrequencyOptim(image):
    dimension = 3        
    offset =(1,0,0) # offset can be any vector-like data  
    translation_dx = sitk.TranslationTransform(dimension, offset)
    image_dx=reechantillonage(image,translation_dx,0)      
    offset =(0,1,0) # offset can be any vector-like data  
    translation_dy = sitk.TranslationTransform(dimension, offset)
    image_dy=reechantillonage(image,translation_dy,0)
    offset =(0,0,1) # offset can be any vector-like data  
    translation_dz = sitk.TranslationTransform(dimension, offset)
    image_dz=reechantillonage(image,translation_dz,0)
    imageSF=(image-image_dx)**2+(image-image_dy)**2+(image-image_dz)**2
    imageSF=sitk.SumProjection(imageSF,2)
    imageSF=sitk.SumProjection(imageSF,1)
    imageSF=sitk.SumProjection(imageSF,0)
    #SF=imageSF.GetPixel(0,0,0)**0.5 #for testing
    SF=(imageSF.GetPixel(0,0,0)/(image.GetSize()[0]*image.GetSize()[1]*image.GetSize()[2]))**0.5
    return SF
3

@phcerdan might give some advice.

4

Hi! I don’t think I can, not familiar with those sitk functions, or reechantillonage.

@drcjaudet Are you developing those functions for yourself or getting them from elsewhere?
@blowekamp might have a better idea about if you are using correctly those sitk functions.

5

Hello @drcjaudet

Not sure about your code, as the function reechantillonage is not shown. I suspect that it is a call to resample with zero padding. You should be able to do this via numpy (we don’t want to enforce the spatial overlap between images). I believe the code below does what you want for an arbitrary dimensional image (haven’t tested it so be cautious):

import SimpleITK as sitk
import numpy as np


def SpatialFrequencyOptim(image):
    dim = image.GetDimension()
    size = image.GetSize()
    #work in numpy because we don't want to enforce spatial overlap
    npa = sitk.GetArrayViewFromImage(image)
    sq_diff = 0.0
    for i in range(dim): #iterate over all image dimensions
        slc1 = [slice(None)]*dim
        slc1[i] = slice(0,size[i]-1)
        slc2 = [slice(None)]*dim
        slc2[i] = slice(1,size[i])
        sq_diff+= np.sum((npa[tuple(slc2)] - npa[tuple(slc1)])**2)
    return sq_diff/np.prod(size)
        
dim=3
img = sitk.GaussianSource(outputPixelType=sitk.sitkUInt8, size=[128]*dim, sigma=[20]*dim, mean=[60]*dim)
sitk.Show(img)

res = SpatialFrequencyOptim(img)
2 Likes
6

Thank a lot for the tips and code.
There was a problem with my image since npa.shape=(432L, 288L, 288L) was different than size(288,288,432); It was resolved with the following code were i just pass the an np matrix.

def SpatialFrequencyOptim2(matrix):
    sq_diff = 0.0
    size=matrix.shape
    dim=len(size)  
    for i in range(dim): #iterate over all image dimensions
        slc1 = [slice(None)]*dim
        slc1[i] = slice(0,size[i]-1)
        slc2 = [slice(None)]*dim
        slc2[i] = slice(1,size[i])
        sq_diff+= np.sum((matrix[tuple(slc2)]- matrix[tuple(slc1)])**2)
    return sq_diff/np.prod(size)

the different time calculation were:
solution 1, time=60.38s
solution 2, time=1.96 secondes
solution 3, Zivy one, time=1.87 secondes
thank a lot :slight_smile: )

2 Likes