fft of linear array

[imaggie]

Hi,

I am trying to use IM to produce the frequency spectrum of a linear array.

Code: Select all

w=`wc -l <$DATA`
mid=$(($w/2))

(echo "P2  ${w} 1 65535  " && cat $INTDATA) | \
  convert - -scale ${w}x${w}! -fft -delete 1 -crop ${w}x1+0+$(($w/2)) +repage PGM: - | pnmtopnm -plain

this produces credible output with what seems to be a non critical error.

Code: Select all

convert: improper image header `' @ error/pnm.c/ReadPNMImage/290.
P2
306 1
65535
0 26 24 24 33 30 20 25 17 15 14 6 9 6 14 
8 4 8 7 21 19 16 41 47 41 56 41 29 34 31 
40 19 24 37 8 6 17 28 14 24 40 16 9 13 37 
24 17 19 14 14 36 37 20 12 3 19 21 40 33 15 
19 12 15 6 21 33 30 38 48 39 45 39 22 14 11 
15 10 47 35 18 26 39 15 47 64 19 31 55 103 65 
27 41 20 60 93 24 59 58 66 81 35 39 88 25 85 
119 99 76 69 42 27 36 49 110 125 58 17 118 122 95 
54 43 29 136 119 76 157 226 181 104 373 344 44 201 337 
191 125 663 386 823 994 345 643 1308 1374 65 1329 1701 2774 3476 
1093 1603 4629 36881 4629 1603 1093 3476 2774 1701 1329 65 1374 1308 643 
345 994 823 386 663 125 191 337 201 44 344 373 104 181 226 
157 76 119 136 29 43 54 95 122 118 17 58 125 110 49 
36 27 42 69 76 99 119 85 25 88 39 35 81 66 58 
59 24 93 60 20 41 27 65 103 55 31 19 64 47 15 
39 26 18 35 47 10 15 11 14 22 39 45 39 48 38 
30 33 21 6 15 12 19 15 33 40 21 19 3 12 20 
37 36 14 14 19 17 24 37 13 9 16 40 24 14 28 
17 6 8 37 24 19 40 31 34 29 41 56 41 47 41 
16 19 21 7 8 4 8 14 6 9 6 14 15 17 25 
20 30 33 24 24 0 

I then have to parse this with awk into separate lines since the line breaks here seem somewhat arbitrary. When I plot this (gnuplot) I get the expected central reflected image but I'm not sure how to scale it and get a calibrated frequency scale to add some meaning to the result.

Could I use something like kerneltoimage to plot this directly using IM ?


As a side issue I would like to know why there is the error message, it seems IM is producing something not totally expected by the netpbm tools.

Thanks for any advice.

[fmw42]

I don't know about your NetPBM problem. Anthony would have to reply on that.

But Anthony has a script to plot one-row images, called im_profile. It uses Gnuplot. It can be found at http://www.imagemagick.org/Usage/scripts/ See examples at http://www.imagemagick.org/Usage/transf ... luate_math

Alternately, I have a bash unix IM script, called profile that can do that also, but without Gnuplot. See link below for my scripts.

[imaggie]

Thanks Fred, that script did not do really what I wanted .

Just for reference I did manage to get IM to do my FFT on some time series data by convincing it it was one pixel high grey-scale image. I had to scale the data to 65535 first

Code: Select all

(echo "P2  ${w} 1 65535  " &&  cat $DATA) |   convert - -scale ${w}x${w}! -fft -delete 1 \
  -crop ${w}x1+0+$(($w/2)) +repage PGM: - | pnmtopnm -plain > $PGMFILE

Then parse the output with awk to format it suitably for gnuplot to get a nice clean graph.

it limits me to integer precision but that serves what I'm doing.

[anthony]

I then have to parse this with awk into separate lines since the line breaks here seem somewhat arbitrary.

Not certain. pbmplus probably only does this for single row images.
One thing you can do is simply filter the ASCII form (either from convert -compress none or from pnmtopnm -plain) through... tr -s ' ' '\012'.
That replaces any sequence of whitespace with a single newline, so each number (including header numbers) is on a line by itself! Should make things easier!

When I plot this (gnuplot) I get the expected central reflected image but I'm not sure how to scale it and get a calibrated frequency scale to add some meaning to the result.

I would probably use a HDRI version of Image magick (using doubles for image values rather than 16 bit integers), as you will get a much more a accurate representation.

Note the scale you have is 16bit and as such 65525 represents white or a normalized value of 1.0
when mapping source image to fft image.

Could I use something like kerneltoimage to plot this directly using IM ?

No. the 'kernel' in this case is a array of doubles used for morphology and convolution, not for FFT which only uses images.

Gnuplot is good, and if you strip the 'header' lines from a pfm (pbmplus floating point - no ascii version, only raw versions) gnuplot can then read the values directly from the saved file (as floats), check the PFM format manpage.

As fred mentioned I have a number of gnuplot scripts (like im_profile) that read raw data output form Imagemagick! This is just like compressed pbmplus but without the header info.

As a side issue I would like to know why there is the error message, it seems IM is producing something not totally expected by the netpbm tools.

Not certain unless I can see the contents of $INTDATA
The error is on read of your generated PGM image!

[imaggie]

Thanks for the hints, Anthony.

The PGM error was my bad, sorry.

I'm using HDRI and built everything using ASCII so it was easier to debug. Since absolute accuracy is probably not critical it'll probably stay like that now. But thanks for pointing out PG format , I was unaware of that and thought PGM was limited to integers.

I'm quite pleased at having been able to get IM to do my FFT on time series data. It saved me having to interface to FFTW myself and is ideal to script processing.

I may rewrite for f.p. later as that would avoid a scaling step I'm having to do now.

Thanks again.

[anthony]

Check the manpage for PFM format. But you can also specify raw output too using channel formats like 'red' or just 'r'
to get floating point you will need to tel IM you want floating point and 32 (floats)or 64 (doubles).
See HDRI file formats...
http://www.imagemagick.org/Usage/basics/#hdri_formats

[imaggie]

Hi Anthony,

there is definitely a problem with a spurious error on PGM input. I thought it was me before since I had let lines go beyond 70 chars, which I found is not strictly applying the spec.

However, having checked it all over and put one datum on each line I still get the error about headers.

What is more telling is that I get the same error if I feed IM output back into IM. It gives the same error message.

Code: Select all

 

$ wc <tmp.dat
121 124 855

$ head !$
head tmp.dat
P2 120 1 65535
341380
321841
355531
284366
255716
280377
326210
266110
251791


$ cat tmp.dat | convert - -scale ${w}x${w}! -fft -delete 1   -crop ${w}x1+0+$(($w/2)) +repage PGM: - | pnmtopnm -plain >tmp.txt
convert: improper image header `' @ error/pnm.c/ReadPNMImage/290.


$ cat tmp.txt | convert - -scale ${w}x${w}! -fft -delete 1   -crop ${w}x1+0+$(($w/2)) +repage PGM: - | pnmtopnm -plain
convert: improper image header `' @ error/pnm.c/ReadPNMImage/290.
P2
120 1
65535
0 192 58 116 73 3 33 89 104 24 93 69 32 205 103 
69 55 87 88 184 242 355 467 481 599 413 467 324 214 287 
254 312 560 669 763 915 1082 997 844 918 713 758 708 677 702 
589 567 856 825 772 838 921 665 835 796 928 1115 1436 1662 2271 
3972 2271 1662 1436 1115 928 796 835 665 921 838 772 825 856 567 
589 702 677 708 758 713 918 844 997 1082 915 763 669 560 312 
254 287 214 324 467 413 599 481 467 355 242 184 88 87 55 
69 103 205 32 69 93 24 104 89 33 3 73 116 58 0 

It seems it is never happy !

[anthony]

Your image values are larger than the PBM limit both by spec (16-bit) and by the file MaxValue (65535)

That is the first value 341380 is TOO BIG!

However even on my machine I get NO error!!!

Code: Select all

$ cat tmp.dat
P2 1 1 65535
341380

$ cat tmp.dat | convert - txt:
# ImageMagick pixel enumeration: 1,1,65535,rgb
0,0: (13700,13700,13700)  #358435843584  rgb(20.9049%,20.9049%,20.9049%)

No errors! though the value obvious got truncated, as it can't fit in my Q16 version of IM (which also has a 16-bit limit on image values).

If I change the value to something 'in range', it comes out correct!

[imaggie]

That is the first value 341380 is TOO BIG!
However even on my machine I get NO error!!!

I did not spot that but it's irrelevant to this issue. Spec says to be as flexible as possible on input and IM seems to be complying and making a sensible truncation.


OK , it looks like I messed up trying to extract a test data set. Let's try again :

Code: Select all

P2 512 1 65535
37163
40185
39862
40118
38655
37769
38918
39862
37072
36675
35821
37949
36741
38971
41345
39348
43167
46731
47874
49393
49793
51812
52037
52209
49762
52232
52342
52664
52922
51859
54063
52738
44679
44578
45633
43726
39211
41390
43595
42848
40724
41447
41668
39974
41678
38744
39827
38750
39639
39557
40235
42113
41379
40861
41232
44560
43595
42076
43575
40241
41318
42831
47781
46649
44085
47810
48021
49264
48197
47431
47476
45275
47607
46186
51380
48689
48615
49405
46952
47906
46172
48529
45383
48500
46778
47751
49270
46925
46759
49254
50448
52441
51859
51183
48787
45654
47480
47190
45654
48162
45754
44480
48316
49534
49297
47968
48005
46759
47409
46239
48101
48826
49176
50407
51597
53432
50524
49203
46968
44675
44904
45973
44244
46399
47044
47458
48719
51996
49702
49481
49496
49168

If there's no cut and paste errors that should conform to the following:

Code: Select all

$ wc <IMtest.dat
133 136 807

ie 132 lines of data plus the header.

Code: Select all

$ echo $w
132

$ cat IMtest.dat| convert - -scale ${w}x${w}! -fft -delete 1  -crop ${w}x1+0+$(($w/2)) +repage PGM: - | pnmtopnm -plain
convert: improper image header `' @ error/pnm.c/ReadPNMImage/290.
P2
132 1
65535
0 131 85 126 85 104 39 57 172 174 169 80 142 97 42 
77 157 25 117 48 92 58 99 174 169 66 113 95 233 169 
26 111 170 322 72 72 122 235 71 109 38 144 60 169 247 
127 144 46 129 239 105 189 158 58 366 730 280 203 263 327 
644 617 1264 1799 1010 1171 45787 1171 1010 1799 1264 617 644 327 263 
203 280 730 366 58 158 189 105 239 129 46 144 127 247 169 
60 144 38 109 71 235 122 72 72 322 170 111 26 169 233 
95 113 66 169 174 99 58 92 48 117 25 157 77 42 97 
142 80 169 174 172 57 39 104 85 126 85 0 

The killing is (and the bit you failed to notice in the last line of my last messages:

What is more telling is that I get the same error if I feed IM output back into IM. It gives the same error message.

Code: Select all

 $ cat IMtest.dat| convert - -scale ${w}x${w}! -fft -delete 1  -crop ${w}x1+0+$(($w/2)) +repage PGM: - | pnmtopnm -plain >IM_own_output.txt
convert: improper image header `' @ error/pnm.c/ReadPNMImage/290.

$ cat IM_own_output.txt | convert - -scale ${w}x${w}! -fft -delete 1  -crop ${w}x1+0+$(($w/2)) +repage PGM: - | pnmtopnm -plain                   
convert: improper image header `' @ error/pnm.c/ReadPNMImage/290.
P2
132 1
65535
0 327 336 340 324 334 338 338 329 327 328 329 328 334 337 
325 336 344 348 356 354 353 349 343 343 343 347 353 347 361 
360 359 346 335 346 340 323 321 335 333 324 341 347 341 323 
324 308 316 311 313 318 319 330 340 352 350 352 359 360 377 
373 390 409 424 435 460 574 460 435 424 409 390 373 377 360 
359 352 350 352 340 330 319 318 313 311 316 308 324 323 341 
347 341 324 333 335 321 323 340 346 335 346 359 360 361 347 
353 347 343 343 343 349 353 354 356 348 344 336 325 337 334 
328 329 328 327 329 338 338 334 324 340 336 0 

So the issue remains exactly as I posted it last time.


Now if you're seeing none of this maybe we need to look at exactly what build does create it on which platforms.

Code: Select all

Version: ImageMagick 6.7.4-0 2012-01-03 Q16 http://www.imagemagick.org
Copyright: Copyright (C) 1999-2011 ImageMagick Studio LLC
Features:   HDRI 

The above and previous posts were created on this version . I have just now rebuild IM to imagemagick-6.7.4.8

Code: Select all

 $ convert -version
Version: ImageMagick 6.7.4-8 2012-01-31 Q16 http://www.imagemagick.org
Copyright: Copyright (C) 1999-2012 ImageMagick Studio LLC
Features:   HDRI  

$ cat IM_own_output.txt | convert - -scale ${w}x${w}! -fft -delete 1  -crop ${w}x1+0+$(($w/2)) +repage PGM: - | pnmtopnm -plain 
convert: improper image header `' @ error/pnm.c/ReadPNMImage/293.
P2
132 1
65535
0 327 336 340 324 334 338 338 329 327 328 329 328 334 337 
325 336 344 348 356 354 353 349 343 343 343 347 353 347 361 
360 359 346 335 346 340 323 321 335 333 324 341 347 341 323 
324 308 316 311 313 318 319 330 340 352 350 352 359 360 377 
373 390 409 424 435 460 574 460 435 424 409 390 373 377 360 
359 352 350 352 340 330 319 318 313 311 316 308 324 323 341 
347 341 324 333 335 321 323 340 346 335 346 359 360 361 347 
353 347 343 343 343 349 353 354 356 348 344 336 325 337 334 
328 329 328 327 329 338 338 334 324 340 336 0
 

Same result . Seems clear there is a problem.

[magick]

Try this command:

• cat IMtest.dat| convert - -scale ${w}x${w}! -fft -delete 1 -crop ${w}x1+0+$(($w/2)) +repage PGM:- | pnmtopnm -plain >IM_own_output.txt

Did you intend a space between PGM: and -?

[imaggie]

Thanks magick, that was indeed where I was going wrong.

While we're here I should flag that this usage seems to be displaying a significant memory leak.

I have a script doing the above on an "image" of about 1994x1 grey pixels. It seems to be producing the correct results as verified against non FFTW based FFT software. That suggests I'm not doing anything else silly

The script gets run in a loop an equal number of times, so 1994 transforms of 1994x1 image.

After running the loop half a dozen times the system (LInux) becomes very unresponsive other programs start reporting errors and I see a lot missing free memory.

If I close the shell from which the script was run I regain about 250MB (sic) of memory as reported by free command and all else returns to normal .

It maybe worth checking with valgrind or something.

thanks again for the usage tip. BTW where is the doc on that ?