home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 18 REXX / 18-REXX.zip / mathapps.zip / Convolution.cmd < prev    next >
OS/2 REXX Batch file  |  1998-08-28  |  9KB  |  305 lines
  1. /* Convolve a square kernel with an array of data.                           */
  2. /* Kernels must be square.  If the sum of the kernel weights does not equal  */
  3. /* 0 or 1, then the weights are adjusted such that the sum of weights=1.     */
  4. /* Doug Rickman  March 10, 1998; mod. Aug. 5, 1998                           */
  5. signal on Halt
  6. signal on NotReady
  7.  
  8. if rxfuncquery('sysloadfuncs') then do           /* this will start rexxutil */
  9.     CALL RxFuncAdd 'SysLoadFuncs', 'RexxUtil', 'SysLoadFuncs' 
  10.     CALL SysLoadFuncs
  11.     end
  12.  
  13. arg inData inKernel out
  14. inData=strip(inData)
  15. out=strip(out)
  16.  
  17. /* Check input arguments. */
  18. if inData='' | inData='?' | inData='-?' | inData='/?' then call Help
  19. if out   ='' then call Help
  20.  
  21. if dosisfile(inData)<>1 then do
  22.    say 'The input file: ' inData' is not a valid file.'
  23.     exit
  24.     end /* do */
  25.  
  26. if dosisfile(inKernel)<>1 then do
  27.    say 'The input file: ' inKernel' is not a valid file.'
  28.     exit
  29.     end /* do */
  30.  
  31. if dosisfile(out) then do
  32.    say 'The file ' out' already exists, do you want it overwritten?'
  33.     say 'Enter a "y" for yes, "h" will give help,'
  34.     say 'any other response will abort processing.'
  35.     key=translate(sysgetkey())
  36.     say
  37.     select
  38.         when key='Y' then do
  39.             rc=dosdel(out)
  40.             if rc=1 then nop
  41.                 else do
  42.                say 'The file 'out' could not be deleted.  Goodbye, Sweet Prince.'
  43.                 exit
  44.                 end /* else do */
  45.             end /* when key='Y' then ... */
  46.  
  47.         when key='H' then call Help
  48.         otherwise exit
  49.         end /* select */
  50.  
  51.     end /* if dosisfile(out) */
  52.  
  53. /* --------------------------------------------------------------------------*/
  54. /* --- begin MAIN                                               -------------*/
  55.  
  56. Call ReadKernel
  57. Call ReadData
  58. Call PadData
  59. Call Convolve
  60. exit
  61. /* --- end MAIN                                                 -------------*/
  62. /* --------------------------------------------------------------------------*/
  63.  
  64.  
  65. /* --------------------------------------------------------------------------*/
  66. /* --- begin subroutine - Kernel:                               -------------*/
  67. ReadKernel:
  68. /* Read the kernel. */
  69. i=0
  70. j=0
  71. NWeights=0
  72. SWeights=0
  73.  
  74. /* Determine the number of weights in a line of the filter.  Filters must    */
  75. /* have the same number of lines as elements.  Minimum kernel size is 3x3.   */
  76.  
  77. data=linein(inKernel)
  78. KernelSize=words(data)
  79. rc=lineout(inKernel)
  80. Pad=(KernelSize-1)/2
  81.  
  82. /* Read in the weights and store in an array centered at 0,0 */
  83. do j=0 to KernelSize-1
  84.    y=j-Pad
  85.    data=linein(inKernel)
  86.    do i=0 to KernelSize-1
  87.       x=i-Pad
  88.       parse var data k.x.y data
  89.       if k.x.y<>0 then do
  90.          NWeights=NWeights+1
  91.          SWeights=SWeights+k.x.y
  92.          end
  93.       end i
  94.    end j
  95. rc=lineout(inKernel)
  96.  
  97. /* echo the kernel */
  98. say 'Kernel weights:'
  99. do j=0 to KernelSize-1
  100.    y=j-Pad
  101.    dataout=''
  102.    do i=0 to KernelSize-1
  103.       x=i-Pad
  104.       dataout=dataout right(k.x.y,3)', '
  105.       end i
  106.    say dataout
  107.    end j
  108. say 'Number of weights: 'NWeights '  Sum of weights: 'SWeights
  109.  
  110. /* Condition the kernel weights */
  111. select
  112.    when SWeights=0 then nop
  113.    when SWeights=1 then nop
  114.    otherwise  
  115.       do j=0 to KernelSize-1
  116.          y=j-Pad
  117.          do i=0 to KernelSize-1
  118.             x=i-Pad
  119.             k.x.y=k.x.y/SWeights
  120.             end i
  121.          end j
  122.    end /* select */
  123. return
  124. /* --- end subroutine   - ReadKernel:                           -------------*/
  125. /* --------------------------------------------------------------------------*/
  126.  
  127. /* --------------------------------------------------------------------------*/
  128. /* --- begin subroutine - ReadData:                             -------------*/
  129. ReadData:
  130. /* Read in the data */
  131.  
  132. j=0
  133. do while lines(inData)>0
  134.    j=j+1
  135.    data=linein(inData)
  136.    if data='' then do /* A blank line was left at the end of the data. */
  137.       j=j-1
  138.       leave
  139.       end
  140.    i=0
  141.    do while data<>''
  142.       i=i+1
  143.       parse var data d.i.j data
  144.       end /* do while data ... */
  145.    end /* do while lines(inData) ... */
  146. rc=lineout(inData)
  147. DataSizei=i
  148. DataSizej=j
  149.  
  150. /* echo the data */
  151. say 
  152. say 'Data values:'
  153. do j=1 to DataSizej
  154.    dataout=''
  155.    do i=1 to DataSizei
  156.       dataout=dataout format(d.i.j,3,3)', '
  157.       end i
  158.    say dataout
  159.    end j
  160. return
  161. /* --- end subroutine  - ReadData:                              -------------*/
  162. /* --------------------------------------------------------------------------*/
  163.  
  164. /* --------------------------------------------------------------------------*/
  165. /* --- begin subroutine - PadData:                              -------------*/
  166. /* Pad the data to provide values to kernel along edges of data. */
  167. PadData:
  168.  
  169. /* Replicate the first line as needed. */
  170. do j=1-Pad to 0
  171.    do i=1 to DataSizei
  172.       d.i.j=d.i.1
  173.       end i
  174.    end j
  175.  
  176. /* Replicate the last line as needed. */
  177. do j=DataSizej+1 to DataSizeJ+Pad
  178.    do i=1 to DataSizei
  179.       d.i.j=d.i.DataSizej
  180.       end i
  181.    end j
  182.  
  183. /* Replicate the first column as needed. */
  184. do j=1-Pad to DataSizeJ+Pad
  185.    do i=1-Pad to 0
  186.       d.i.j=d.1.j
  187.       end i
  188.    end j
  189.  
  190. /* Replicate the last column as needed. */
  191. do j=1-Pad to DataSizeJ+Pad
  192.    do i=DataSizei+1 to DataSizei+Pad
  193.       d.i.j=d.DataSizei.j
  194.       end i
  195.    end j
  196.  
  197. /* Remove the comments to see the padded array of data values.               */
  198. /*
  199. say 
  200. say 'Data values:'
  201. do j=1-Pad to DataSizej+Pad
  202.    dataout=''
  203.    do i=1-Pad to DataSizei+Pad
  204.       dataout=dataout right(d.i.j,3)', '
  205.       end i
  206.    say dataout
  207.    end j
  208. */
  209. return
  210.  
  211. /* --- end subroutine  - PadData:                               -------------*/
  212. /* --------------------------------------------------------------------------*/
  213.  
  214. /* --------------------------------------------------------------------------*/
  215. /* --- begin subroutine - Convolve:                             -------------*/
  216. Convolve:
  217. /* Loop kernel through the elements (i) then through lines (j) of the data */
  218. /* Build one line of output at a time.                                     */
  219. Output=''
  220. do j=1 to DataSizeJ
  221.    Output=''
  222.    do i=1 to DataSizei
  223.       Sum=0   
  224.       /* for each i,j do a kernel */
  225.       do y=0-Pad to Pad
  226.          kj=j+y
  227.          do x=0-Pad to Pad
  228.             ki=i+x
  229.             Sum=Sum+ (k.x.y*d.ki.kj)
  230.             end x
  231.          end y
  232.       Output=Output format(Sum,3,3)', '
  233.       end i
  234.    rc=lineout(out, Output)
  235.    end j
  236. rc=lineout(out)  /* Closes the output file */
  237.  
  238. /* Show the result */
  239. say 
  240. say 'Filtered array:'
  241. do j=1 to DataSizej
  242.    say linein(out)
  243.    end j
  244. rc=lineout(out)
  245.  
  246. return
  247.  
  248. /* --- end subroutine   - Convolve:                             -------------*/
  249. /* --------------------------------------------------------------------------*/
  250.  
  251. /* --------------------------------------------------------------------------*/
  252. /* --- begin subroutine - Help:                                 -------------*/
  253. Help:
  254. rc= charout(,'1b'x||'[31;7m'||'Convolution:'||'1b'x||'[0m'||'0d0a'x)
  255. say 'Convolve a square kernel with an array of data.  A teaching tool.'
  256.  
  257. say ''
  258. rc= charout(,'1b'x||'[33;1m'||'usage:'||'1b'x||'[0m')
  259. say ' Convolution data kernel OutputFile'
  260. say ''
  261.  
  262. rc= charout(,'1b'x||'[33;1m'||'where:'||'1b'x||'[0m')
  263. say ' data = Array of data values to be filtered.'
  264. say '     kernel = square array used as the kernel.'
  265. say '     Output = file name into which filtered array is written.'
  266. say ''
  267.  
  268. rc= charout(,'1b'x||'[33;1m'||'Exam: '||'1b'x||'[0m')
  269. say ' convolution Conv.data Conv.kernel Conv.out '
  270. say ''
  271.  
  272. rc= charout(,'1b'x||'[33;1m'||'notes:'||'1b'x||'[0m')
  273. say ' The data array is padded by replication of the first & last rows & '
  274. say ' columns.  Values in the data and kernel are blank separated, one row per'
  275. say ' line.  If the sum of kernel weights equals either 0 or 1 the weights are'
  276. say ' used as input, otherwise individual weights are divided by the sum of all'
  277. say ' weights.   The kernel must be square, i.e. equal number of rows and '
  278. say ' columns.'
  279. say ''
  280.  
  281. say ''
  282. say 'Doug Rickman  August 5, 1998'
  283. exit
  284. return
  285.  
  286. /* --- end  subroutine - Help:                                  -------------*/
  287. /* --------------------------------------------------------------------------*/
  288.  
  289. /* --------------------------------------------------------------------------*/
  290. /* --- begin subroutine - Halt:                                 -------------*/
  291. Halt:
  292. say 'This is a graceful exit from a Cntl-C'
  293. exit
  294. /* --- end  subroutine - Halt:                                  -------------*/
  295. /* --------------------------------------------------------------------------*/
  296. /* --- begin subroutine - NotReady:                             -------------*/
  297. NotReady:
  298. say 'It would seem that you are pointing at non-existant data.  Oops.  Bye!'
  299. exit
  300. /* --- end  subroutine - NotReady:                              -------------*/
  301. /* --------------------------------------------------------------------------*/
  302.  
  303.  
  304.  
  305.