# Calculation of spatial frequency

**URL:** <https://discourse.itk.org/t/calculation-of-spatial-frequency/2874>\
**Category:** Algorithms\
**Tags:** simpleitk\
**Created:** [March 26, 2020, 7:35am UTC](https://discourse.itk.org/t/calculation-of-spatial-frequency/2874 "2020-03-26T07:35:26Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![drcjaudet](https://discourse.itk.org/letter_avatar_proxy/v4/letter/d/a698b9/32.png) [@drcjaudet](https://discourse.itk.org/u/drcjaudet)\
**Post date:** [March 26, 2020, 7:35am UTC](https://discourse.itk.org/t/calculation-of-spatial-frequency/2874/1 "2020-03-26T07:35:26Z")

</div>

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](https://discourse.itk.org/uploads/default/original/2X/0/0e5fe6f46f43d133571c738e069150a021c4b6c8.png)  
with M, N the voxel indice and F his value.  
So i create a code like this which is suboptimal.

```python
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

---

<div class="post-metadata">

**Author:** ![drcjaudet](https://discourse.itk.org/letter_avatar_proxy/v4/letter/d/a698b9/32.png) [@drcjaudet](https://discourse.itk.org/u/drcjaudet)\
**Post date:** [March 26, 2020, 12:24pm UTC](https://discourse.itk.org/t/calculation-of-spatial-frequency/2874/2 "2020-03-26T12:24:06Z")

</div>

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

```python
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

```

---

<div class="post-metadata">

**Author:** ![dzenanz](https://discourse.itk.org/user_avatar/discourse.itk.org/dzenanz/32/1093_2.png) [@dzenanz](https://discourse.itk.org/u/dzenanz)\
**Post date:** [March 26, 2020, 12:56pm UTC](https://discourse.itk.org/t/calculation-of-spatial-frequency/2874/3 "2020-03-26T12:56:01Z")

</div>

@phcerdan might give some advice.

---

<div class="post-metadata">

**Author:** ![phcerdan](https://discourse.itk.org/user_avatar/discourse.itk.org/phcerdan/32/286_2.png) [@phcerdan](https://discourse.itk.org/u/phcerdan)\
**Post date:** [March 26, 2020, 1:27pm UTC](https://discourse.itk.org/t/calculation-of-spatial-frequency/2874/4 "2020-03-26T13:27:35Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![zivy](https://discourse.itk.org/user_avatar/discourse.itk.org/zivy/32/1726_2.png) [@zivy](https://discourse.itk.org/u/zivy)\
**Post date:** [March 26, 2020, 1:40pm UTC](https://discourse.itk.org/t/calculation-of-spatial-frequency/2874/5 "2020-03-26T13:40:43Z")

</div>

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):

```auto
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)

```

---

<div class="post-metadata">

**Author:** ![drcjaudet](https://discourse.itk.org/letter_avatar_proxy/v4/letter/d/a698b9/32.png) [@drcjaudet](https://discourse.itk.org/u/drcjaudet)\
**Post date:** [March 27, 2020, 4:43pm UTC](https://discourse.itk.org/t/calculation-of-spatial-frequency/2874/6 "2020-03-27T16:43:18Z")

</div>

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.

```python
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 🙂 )