From help-request at octave dot org Sat Jul 16 11:41:47 2005 Subject: bicubic spline interpolation in octave From: Hoxide Ma To: help at octave dot org Date: Sun, 17 Jul 2005 00:40:12 +0800 (CST) --0-424385468-1121532012=:8310 Content-Type: text/plain; charset=gb2312 Content-Transfer-Encoding: 8bit Content-Id: Content-Disposition: inline When I try to use interp2.m in octave-forge, i found that in the source code, there is no bicubic interpolation method. "At the moment only 'linear' and 'nearest' methods" ! there is the code: %------------------------------------------------------ function F = cubic(arg1,arg2,arg3,arg4,arg5) %CUBIC 2-D bicubic data interpolation. % CUBIC(...) is the same as LINEAR(....) except that it uses % bicubic interpolation. % % This function needs about 7-8 times SIZE(XI) memory to be available. % % See also LINEAR. % Clay M. Thompson 4-26-91, revised 7-3-91, 3-22-93 by CMT. % Based on "Cubic Convolution Interpolation for Digital Image % Processing", Robert G. Keys, IEEE Trans. on Acoustics, Speech, and % Signal Processing, Vol. 29, No. 6, Dec. 1981, pp. 1153-1160. do_fortran_indexing =1 if nargin==1, % cubic(z), Expand Z [nrows,ncols] = size(arg1); s = 1:.5:ncols; sizs = size(s); t = (1:.5:nrows)'; sizt = size(t); s = s(ones(sizt),:); t = t(:,ones(sizs)); elseif nargin==2, % cubic(z,n), Expand Z n times [nrows,ncols] = size(arg1); ntimes = floor(arg2); s = 1:1/(2^ntimes):ncols; sizs = size(s); t = (1:1/(2^ntimes):nrows)'; sizt = size(t); s = s(ones(sizt),:); t = t(:,ones(sizs)); elseif nargin==3, % cubic(z,s,t), No X or Y specified. [nrows,ncols] = size(arg1); s = arg2; t = arg3; elseif nargin==4, error('Wrong number of input arguments.'); elseif nargin==5, % cubic(x,y,z,s,t), X and Y specified. [nrows,ncols] = size(arg3); mx = prod(size(arg1)); my = prod(size(arg2)); if any([mx my] ~= [ncols nrows]) & ... ~isequal(size(arg1),size(arg2),size(arg3)) error('The lengths of the X and Y vectors must match Z.'); end if any([nrows ncols]<[3 3]), error('Z must be at least 3-by-3.'); end s = 1 + (arg4-arg1(1))/(arg1(mx)-arg1(1))*(ncols-1); t = 1 + (arg5-arg2(1))/(arg2(my)-arg2(1))*(nrows-1); end if any([nrows ncols]<[3 3]), error('Z must be at least 3-by-3.'); end if ~isequal(size(s),size(t)), error('XI and YI must be the same size.'); end % Check for out of range values of s and set to 1 sout = find((s<1)|(s>ncols)); if length(sout)>0, s(sout) = ones(size(sout)); end % Check for out of range values of t and set to 1 tout = find((t<1)|(t>nrows)); if length(tout)>0, t(tout) = ones(size(tout)); end % Matrix element indexing ndx = floor(t)+floor(s-1)*(nrows+2); % Compute intepolation parameters, check for boundary value. if isempty(s), d = s; else d = find(s==ncols); end s(:) = (s - floor(s)); if length(d)>0, s(d) = s(d)+1; ndx(d) = ndx(d)-nrows-2; end % Compute intepolation parameters, check for boundary value. if isempty(t), d = t; else d = find(t==nrows); end t(:) = (t - floor(t)); if length(d)>0, t(d) = t(d)+1; ndx(d) = ndx(d)-1; end d = []; if nargin==5, % Expand z so interpolation is valid at the boundaries. zz = zeros(size(arg3)+2); zz(1,2:ncols+1) = 3*arg3(1,:)-3*arg3(2,:)+arg3(3,:); zz(2:nrows+1,2:ncols+1) = arg3; zz(nrows+2,2:ncols+1) = 3*arg3(nrows,:)-3*arg3(nrows-1,:)+arg3(nrows-2,:); zz(:,1) = 3*zz(:,2)-3*zz(:,3)+zz(:,4); zz(:,ncols+2) = 3*zz(:,ncols+1)-3*zz(:,ncols)+zz(:,ncols-1); nrows = nrows+2; ncols = ncols+2; else % Expand z so interpolation is valid at the boundaries. zz = zeros(size(arg1)+2); zz(1,2:ncols+1) = 3*arg1(1,:)-3*arg1(2,:)+arg1(3,:); zz(2:nrows+1,2:ncols+1) = arg1; zz(nrows+2,2:ncols+1) = 3*arg1(nrows,:)-3*arg1(nrows-1,:)+arg1(nrows-2,:); zz(:,1) = 3*zz(:,2)-3*zz(:,3)+zz(:,4); zz(:,ncols+2) = 3*zz(:,ncols+1)-3*zz(:,ncols)+zz(:,ncols-1); nrows = nrows+2; ncols = ncols+2; end ndx % Now interpolate using computationally efficient algorithm. t0 = ((2-t).*t-1).*t; t1 = (3*t-5).*t.*t+2; t2 = ((4-3*t).*t+1).*t; t(:) = (t-1).*t.*t; F = ( zz(ndx).*t0 + zz(ndx+1).*t1 + zz(ndx+2).*t2 + zz(ndx+3).*t ) ... .* (((2-s).*s-1).*s); ndx(:) = ndx + nrows; F(:) = F + ( zz(ndx).*t0 + zz(ndx+1).*t1 + zz(ndx+2).*t2 + zz(ndx+3).*t ) ... .* ((3*s-5).*s.*s+2); ndx(:) = ndx + nrows; F(:) = F + ( zz(ndx).*t0 + zz(ndx+1).*t1 + zz(ndx+2).*t2 + zz(ndx+3).*t ) ... .* (((4-3*s).*s+1).*s); ndx(:) = ndx + nrows; F(:) = F + ( zz(ndx).*t0 + zz(ndx+1).*t1 + zz(ndx+2).*t2 + zz(ndx+3).*t ) ... .* ((s-1).*s.*s); F(:) = F/4; % Now set out of range values to NaN. if length(sout)>0, F(sout) = NaN; end if length(tout)>0, F(tout) = NaN; end endfunction by changeing the interp2.m code in matlab, I can do bicubic spline interpolation in some case, you can see a demo ploted by gnuplot. But, you know, the work haven't finished, it not work with interp2.m, and it use the code of matlab, there may be some problem in copyright. I plan to rewrite this codes. I write this mail to ask for some suggestions. thx. ___________________________________________________________ 雅虎免费G邮箱－中国第一绝无垃圾邮件骚扰超大邮箱 http://cn.mail.yahoo.com/?id=77071 --0-424385468-1121532012=:8310 Content-Type: image/png; name="cub1.png" Content-Transfer-Encoding: base64 Content-Description: pat397867616 Content-Disposition: inline; filename="cub1.png" iVBORw0KGgoAAAANSUhEUgAAAoAAAAHgCAMAAAACDyzWAAABKVBMVEX///8A AACgoKD/AAAAwAAAgP/AAP8A7u7AQADu7gAgIMD/wCAAgECggP+AQAD/gP8A wGAAwMAAYIDAYIAAgABA/4AwYICAYABAQEBAgAAAAICAYBCAYGCAYIAAAMAA AP8AYADjsMBAwIBgoMBgwABgwKCAAACAAIBgIIBgYGAgICAgQEAgQIBggCBg gGBggICAgEAggCCAgICgoKCg0ODAICAAgIDAYACAwODAYMDAgADAgGD/QAD/ QECAwP//gGD/gIDAoADAwMDA/8D/AAD/AP//gKDAwKD/YGAA/wD/gAD/oACA 4OCg4OCg/yDAAADAAMCgICCgIP+AIACAICCAQCCAQICAYMCAYP+AgADAwAD/ gED/oED/oGD/oHD/wMD//wD//4D//8BUJrxzAAAgAElEQVR4nO2di5KrOAxE SVHz/7+8eycPsK1uScYhE9ynavdOjDFgGlmSHbIsQgghhBBCCCGEEEIIIYQQ QgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQ QkzH7Vb/EdsmxBAkQPFRmLakO/F2bv+4//v8eFfe//+Z24QYyVNv2/9+lbf/ p9z2YPX41AWJ76IV2VLrsNwmxEgkQPFReodgIYZwe4Yd27+vUnubEEIIIYQQ QoxH4Yb4KIUAI2r8edupiOmRAMW72a072JYf7POAi5MHlADFEejcx3MxDLSF P3dOPF9xMQIC5COx1CeOEBEglaAEKI7wcveAAPfiNJEAxRGeyiqWH+zk6AYh QhxB0hIfRQIUQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBC nMDz1ySE+BgSoPgk0p/4KHsBakAWZ1O+Wf1jpyFmRQIUn6QccSVAcTISoPhD SIDio0iA4qNIgOKjSIDio0iA4qNIgOKjSIDio0iA4qNIgOKjSIDiZMr1VxLg Oaz/8+lz+BtUv+YpAZ7BXX0S4T+q5acS4Pv5X3YP4a2LJFgtgdaK6Hezbqpb H5+nRkPwqfyT21Nx665sXiTAM1mXVn/ln9MhAZ4I1N/cClQa5hxWYv9mVmCB BPg21tf/FmP8/afOn7NP6Q8iAb6JR+bv+Wm3YftDApQA30Vh/kD8sUqAEuB7 WAP6+7lz4ln9SSTA8awR/f378KNQRAIcznPC19ffIgFKgKNZE/r7PwiRAiXA kbyWHawR/dVFUyIBjmOb9910xfUnBUqAw/iVH9GfLcD5FFi9I1oCHMNa8Crd 11jMDxMKsPz4odO4GGu56q+xhFB/8ymwXpD6uTO5Dsaqg/06wGUh+ptQgOVy LK2IPkox67sbfmP6m0+Bi9YDjqSc9QCL/rj+5lOgBDiKytkDi/4c+zeXALUi ehivSQ+86K8utsff2RSoIGQAW7aFTXo0i59rrU0owBIJsIddsm8zfzuXb+/9 +fqbWoESYJ5dzi+05gB8WBbbSZyLqwhw9Rl2IMvqYf2hLcXneRX4DQLEq4Zr dVn3cW0rdwty3U34LkBAjcmDMiPSnIdvFCDWEBVgXZaTIp3pja25YgnoaRX4 PQKsBGDfsZQCt09YipXudm+5WqP6M7MxzSkYJzkFf16AP//f7Z+fn5i1CwvQ fD/LJrVadc0827IT6+rpz8jGMH9wLv60AB8SsH3AYybQKDXdw92nsriZBWna LDd4gpMA/xZ7HYDv7oRNIFVgozrbzDXye/4b0l9EbzMp8G9PxTXOWEaAcS8Q ByDG60yzf7f681y+md7Z9ndXRENJmJXtJrzCvXNHzqP8CNsyM9HWJyfmrWtc mb+5JP+hiqPGzk+7uHWbGLeSXN2s1ZTxCYUqr5JZBPi/+v7WglSUXNtXsXcM lpJsX+OaNcEtlt9Leq16zE/kEs08zUX5U8uxDGEMNYH+aFsNztXZUfntsywB /QFbSTdckT/zrbiMu5cvbYyeo0CQhcZHi0e8lQCNE1nrmlfn4wK8Z3vBtnip XQimNOC5tKbv2dBegiQsqSwgNunAFnP344p8VoDbxIK9ub/0Kb2EHwl/PaY8 yXpwbj6tYFtjKemjMY0CN84WYJHbtWuAHWlhIJAB0QnasPu7ydFZagTDaF2V m2YJ8K00xmaICbRG3EATRLBOXAQGY9PD80+M+JJTcJoA435ZvDScXalK6ZPQ yI+uqq/GT09/zYk5O1yfcwSY0klIgdTdY03wCKXNClWBA7PhzQBr2L+0hbw6 JwhwTMRb+Xu8LmoCPAgg1Ni1A+Km1EfDREYHhQvzbgE+b/lxdw/4e1EFspm+ Z+PtPlWDnuCoQ9d2hHE2EuBIui2V3VYicd0O2LS26SQ0o60XnJRFYHQlNfgZ XpT3CTA0ERESIF2n4BevjSFrKxvNe/GCaQ5XtHXnQTinYpZelzcJsCsz0pb2 pfd2xZEh2/QS2vN3lvTVBg4edrUroONenXcI0BwrsybwgLv3LA6N2WZiHEQj 2KQ3Twk5deYNMDf1Krz53TAo0Eyopze9FzqPUgv2cgPoDq51gdWm0YJxFHiB Vxfge5dj4YxLUD3H3L1II1BFz00sW+zbN0tEPFqpD29tuRhvEiBfeucX7wzf gZglEPEW1apNfMLGtF7NR8d1QA+Hd94XoRqCh62I7hPPq9cPBM25iPfhYeJt dUnSvpmCY03W1a4uwOUdFjBy41mw4N5EWrw+zwF7ZkU5ll/k+7u+/rzT4A6u BJgmduNTwUKHiIO1yVefWg/CX0FhPTpmZLMv4ed6bQUOD0KiKT9zdINBSyaU gVO85hFxK+1e5oDM5VW3byqYy/raAhychnnb9FhQxZlFMSynbbkQyPwxedVR it9EW+3iAiw5JMDunF9hs7pNYC5qZkGOmQ6E4+Qm4/ZxqNpHpo6P4zMp8IAA e3N+wRkOrql62HUUSI9p+xDm8FueAzraAnS2AGUGfZhL0i3Avpxf4gtrWNzZ qJkHKPbGQDqQngPSMKyLGro6fQJ0k71m+YBZuo5GeP1E+FJtd2TjhRrUBZlI gT0C7Mn5kYA3bgJZIzm5PnZA8nPTgcZp+DVss2edggRISPlez124wQwJsKMR Zw00lp8fFnQUtCbSroYKr0lSgNElUmXEm6gN23AbMYwuOyScjTOEYbmcvkJd I0pObx4FpgSYXucXnR1z3LSY6EvL7GRokHG0UoW2tPlA0DZvmsjcE3lFEgJM zZkt9uxY1gQmoubCw3e8VGgcjf1QxG06dEWB46sQ/2CiddFRAZIJ/Ez9TCOJ Gbai3LW60Dha7i3W/2Yrkf8XegysDZc2gT1TcaH1nbva2WDVkELukFuxazGh dTTz0aw1GEa0Ak8/jdcVYMdiBO7N165O1xRvOX52Li1EgS2JJ6iNspprJ/QC EgX6o8+Rve0CpAXoqeFosFoUx9al5oZ9NliSYdQ8TCtfWoc9i8jC0+f3IiTe ER3PoRyb4gUBR1jFkUQ3CeKBiwb1+vrsaBRbMzTOeo/7BUj4gKGkc3+wWhwp 4yblhv3H3cYPBDZ/cMBGJcZTSC6MNHhdAZYGL2j97p9RtZStGxy0RIZ95ByS TWurEPMZ87oIXBh6YKhPeg2qARcLMGbTDk/xHlqXGs5Rk7vJLS8Vd2Pe8Bhv F7Ih/roCDL0lP+SUmy+1wNWb4kNBS3zYpxmd6Bc0LMNdl5st0eaJfq8qwApb gG44uRTWp0uAh4KW0HsPnlVZM9amznx0yEZW9Zjk5lCgKUB46XvjFdhhgPU6 OuyzJBIY/I2P9GzZEcxTAh9ovctiCJDc3qc56FkesDV/KGhxhn1kpIk0vVAj OEjDSrSviP2cVID2qPTaCka/oAJ7Zthywz400vWNB3Yo4vrWRWCMt40vsnPR EeFqVAJ0E2r91qt7hm3TXMhJf7mY9GTsTeazZeSj290jY7RxYGY/51BgKUDH vnSnXAa8+SW3FJZb6cgKq/0W78jx359B9WxnxSq8GnsB8lV0ve++OvzmlyXz 4o1HddY8zDBhK+SZXmf8pM81dnkmEyAd4Nzh0+5bFHBkBJgOWnDW+fdeo8bA 8Nscy9Qf8uzKj6gvcsPKtXgKEK8TqTcGxdO/ou9QI6w+to2oMeNPa2hwB2ky zsJHAu9wKR4CJGNI+KVnuz8Pxxt9jXAvgd5q63uW1mMXcuLsSnk7d1EBFrPB v3+T1SDpVS7HVvQ9/+1pxEvRYONopp4so4Z7wzldZseJzC6pwFstQOxipRy4 ZcCbJvuDFppzvlfgWxo5W6dGTmFti6yDWA3MZQJvlQXEPcP6zKyelKvdRlcj boqGjOam3YwsvSqroVOBx91qd7gMX80mwJ+fn9v//1mVUs7/Qze5rjfbSPpJ yE2o9YjNk204g0+B4fgxoaEmyT6XVGBlAY0arjNf1uWxQkBTw7+GRKxZ+cEc dtEKVe4lki5DUqIXzW/BN1P5gHi4Cox8kdkxrqljkyT+qiwyOoPUE5SDM85D g4WkRPuZGsfvpglC9hvD4Wc8VmDW60DQEhqzWRAPJEEsKn/anNddkyvF4ryg AI0V0dhc4HgjozXYSLyJ1qtz7tDDO8QtgW1o+K0PFp8bROMse9D9Qegy3POA v3+Gcn403ogKkBqvwKDqpVzu5Wx49odu+9Csnr1uEOzCJz6dE7oSd1u4wltS evNdM8LJRlyXMVadTXBB+TmpFjoymo8UMnP0IiJO0GV4TsWR2/WoEFqN54xC x1Zl2aM2GjPxrWNPj+OTwtDAzhugLovWs491KV6LEVCF1XbX0ibwaNAy6GtI wMy99oqEuuGYAX2kljJm4S/DKxzBQ1nPSNm04ThXvDiVpOYjPPRgTXVYz0x4 Z2jnmBmPhYEXYrcesNnWE3A0vvbhSZJYwLvVZq3DpuITIk7Zbmdk8FkL3Iu8 IkiAI77BcfSbw0s18PoK5H4qFKe9G0rUkOPTYfbhZLJnMar4C7FLShfS2ddJ q8cYcztMYNOIJ0A+fOH4x97NdiaorHe7IZeDXkPGT7kKxXdCFuDzxdWzIn8v K+Lkt0AW07RAyfGh0g6HF1uZcfXax9p2yznb3wt5RW8i0qxL15W7exkTmA18 8N0L+PZgnI06aqZ6cwvZwBkG9vlK6BtSc53wMDy10cuIx4gGelZlOQ8OtqYZ pVuaNRvGZo4JsGPLV8Jf0Zu4863yeBsxK5pvw1UsERnY4k2IUPOJz4WdJ9Hm 9QR4I6/oDShwp7xDxq4o7foWyHbjcH0ytGWUaUYtQH9svEDnQnR2OQHe6EvK me05NNriyr2rsva7pb1+Yv6MjdZHuHCG2UW8bRYT6L4l3xp9ksNtvDQ7SbIP K0yjVLVOB0QWaNgmb19ATpt0lOnLoC14ny/G/5mGl+de6q7HXWfFXVHzK3YI mEzqJYAjmybVTNSkzBwz1dyhpBu+E/+XkmxzN87fOzJJsoAUoeEd4obQM2W7 F/FnDO3CvBbXi7icAEsSb0gdYgIPT5KE8mXcO0Tb7E/A00MWN9xuvc80Y3AJ eEd0vDRc9fgkSTRHyGfw0D0HuWp0xnyEIDYPb5rSBIK35B/TWmOOBkySuIti Nt+VNQTFaTuE5Imhdhw+DaVzyZrYtTChAIco8NgkSXETw4tieKgO1RBe+YKc S9gwPp3Y0j++uPv7gT9U0+F6vyQ3dJIk6DP+bvBmhNG2xHmZNrNn7V+se+j5 XQHyS0lO2Wrg7J42gTkZO/EJVic4q+gr23BXkcsNecOu2/H9kN+KsxQWF1vn aFsePl6bJ50XJD9iNi1l2o4eEjY5/1j2nSr8ItAfKwyWHXYY256nQYtlKYjB 2bbbm1gpCi32n6FBRocMLr3yJxqvAP251jco0DdqO8sXFKCf1KbyY4P8WtVt 63Wt/QtsCC1KvwDO7wUfKOsxgfWwG7i5AWcADs7Abtoht33d8AyZ6AMCrJ8x tMfXkxfgW0ZmEHC4AvQjZDzVDwytnSnMrdEmDgT1LeDRZhXgMbGF7xlcE8MV aO1UiwmO5iCzmJgPWdCKGrYPdRvx5mkF+IZBuHSrXtrzncOq2A8kyeCMRnrg D+JAJW5ZnU27gzGn9HJ4AtwnY7Yyq565MyqM5g2h6XF9f5a/RiMcmg9h10bc hsym7dRyQ8HX4wpw32ssG+gJ0JskCSvQSdE8R2fcEtqWTjtaW4NWNxCGxMq/ Hl+A3N452epw3hoWl2Oln6LxRme0Dcov+HPBVuPBTftKEqBB0Ck5GJx4xeEU DbuJ2Rts+6focUqMuWa/EMs+tQCDYVlfHOIWP13GYO21S5qZJfbIuSQxuRuu byWzmcCYAM9RIHC24n45DaidfLR1PZujwQ6LopGtMDIMeOsOJhegYT3ePwg/ tReuzU4NShPdebDwJRfB4E3Qqk42BkcFuKyNLQoqsGsQLo4VMQr8zEL5aBol cDPXM++GTlgCBDxGk00abzKBxpgbkCvbAec/YFhj2/tUosbbtFM+OovgMb6a uADL+7O66QlaBgZD5O+xm2LtAm0jNTbOOJtL1Czo3OqjBZ+2qyowIUDTRYt0 n6upXUuJ8GQheRWsJWR4dttS57BgL5UnzGGAMpcJzAjQDgDLXHNuED6WpIbx 8b/6P/casCk4GwNscOrd+rsyLsDEvK8EuAAFbh/Q5Ec9eFvTI/AQ4MA0QP5/ ww+1jnAT24ecl2HH8KatStrfuB45AUYclv0Y1zEjF+j+0ITcD7GOJGi1S4Ex w8kTP62CQ5q5FJgU4MhcTK6fXyNnJGb8+flZwS9xkxuPxtm1+cNoJ/jBOxbe QQK883T79p/bGpGyjoAjNKV1H9t+UDtrMmZAxgyPBYF1Byzi4mb4aqQF2Jii cxSIA466+FHRFmB6QgSZOTw3HZi17ph2Iwbzu+kS4HNsS2QDDwgwMyO3yfQn FRhkzRy+Qv/ac1M82y4S4IN12TtXZjg7SIFl2/7t4u5hKGQg7W3bcPTqx7Ud SSfPYH41HQL8TXEUn59SoVqJxyFPNzMg6125Ux+LE0b26IyhbHwvIeAeNmfn PoFfTV6ARoRZmSzbzY+awPSE3OMBoDtgccK4Bp1DNoez25J8pOo9JMAnP04y xs4ze6JyJ+RgsWd5cPBMBlO7SZrW49LMz3rE7OV30yfA5yj5KrHvcJl/xgI8 MCEHDe6+dbQJCze98MCd9+1ZZRGs/dV0CnDzux5/NHWqh5fSHiHqnTvu+doj MrIPUpl7HtEcEj87CXDH3u9bfAWiIhaI0NJSuNnBGYtsgalqqDIvTOiY9Mg5 wd9MtwAb741VIEUpAW7OUcCdZ4OzMx9iPU9IZe7MJPNqkdJTnfLV9AvQGGTZ dlSUNoEgRm4USUZF6DfAVDVUmRcnpCc9WNJdAixpOn91x8WDJhBPyNVJPzIq huZDgplCJ07wF8U0TzGtfkEFHhFgrYctMMGPfbcCnQm5xikF7ZAnJD3t68UJ oSRe4cl41SXAirW8o/ubtgJPpkOBoQm5Tf+kdeIlIK8SjvfAFiPr6AmQPjZe I1/MMQG+0nj7j9tGK8uSE2A0n/3YwN1DEighYwZV5rkC4SRe1YO8tgTY8Orv dfep2uznoy1XHGQIiS3hyTaizfSsG87hYN1ixYaCKt7I93JUgJXr5TrldgJ6 rTdbO4M2X3vSYY6Ik6ad0QbmCXiDc9VYqjY87LdyWIB1PpBsborcGZHYXfAn RFjWmcjPNozQju+3WxvsyvGnjB72WzkowH8/PlwGe53pwFR/Q83b9xgbLEdl qRzOaztzUtszy126BPjc7cHv37UIzOF1z2Gfp/Q8SX2aEcqrzBc72tKWu1JO lH4xWQHuhLc1UcafT8ce3+DDClzt4b6qzh4DvA2bVCdR4nsCsVOjxfMK0FDe xv3Bbx9pmA48pEAYH1fHxu3gbdiBcKytn0XeykMpQgnwUYEq70E5JNb30BJM lxvoJHN29Vk0RE4Gx0FOhBTJIr+e00BfJEu/GCysiPBelHfXuvlVoJs0ga2I yU1m+UOWqsabeEIyKqkFpGgkwKooo7wnpX9DzQ/O2RnWDiekwQ2io3PPNm/W I5F0HpHzu5oCdzoLDbaQXDKmzgC24H1h8WM3KBcyf0EzhXgvzzbWhxgQb1xR gIeE98JKxuw3tztYjdhN+6Xu8E5nSvA2OusRiYXKmrKASyW1MeL7pfS+Xq7b /nOzQ6jMMxyNvTVsIzkiyRYx2XBJW82gllJXXX8p+9v4p7VSgePaLr38Wh4x E5hV4Lh0YGKTZ9H2G/q8CVIqAWKe0msUB9KBxxRI0oHrvhJuiMTVLORmii7L 35Fy+X4BNhMZ4yif++a2eaEJLCtLA+nA16OA26HpQKwbd7h/lfu2OV/68/Nz gy8+/ApuvwosiwY2X1u9Zmslw5wbaORkiAJ5yg9uS+4Gg5dUDBEx+49rH3m7 PkAzBC+DFViZqmrrtmHFTnylYjM1Y1Uuj0yl2bcpeALk0vpMYHH1EiCnsnA0 GWNnAIHkojf0uR8UIEv5dWxq7RQ4MVYMS9sH78sFWKdh7mUjD/B784vsy1p+ bKsfKKzcM2d4Z2aXmC1m0WqDnThfpxQY/W8XoMlwBW6qW6uPoHp34aZxzz0j jwF9QhyTVl9w5HS9UuhwLBJggCqa3H80kzG5SKTZG94s7I4ST7VRAj+lxwbz FHrcPaq9Xy4pwHeYwMfftSdnxRNdN2/XkOOeEetobCpG1MA54YAnawJd7f1y TQG+T4FLYx4eI/PKFUTuaSthGvGyGw6PYW5FAe8Yd8/X3i8XFeBbFViHxu02 I+617Ca8UTB4ILFwekitDkNjK3ZWVRth6d2RACNYCcCtp40eL32gGqsi2H/f DKjOzZaTQ7TrpBXYobwHVxXgcAWausFj7kEnvs7GoOr7bWAsxMeoH6TESe2P cOtT3oPLCvANobA9SP2WhoJhNLDSA7a2sDgwbt/R5rNKb8Cxs3mHeloCDEJt nunOhURpViTe4YLNFhlS2yOzMd08p620OTcJ0OZN2UBg89adKwQHZjoM1z4i kqvj1IHNwTGdFcPH4khPX1iA71PgYmYampjSCjtwMByQ68tspYfUhYzpuHp1 oqjlQz19ZQG+VYG/MijvzFuD4XIHe9ResDbT/t5utA2EGP09LQHGgTbveZ9C 4WR0FdRW2bFbq6NNz99LjbYACdBmfCDCbJ41loYnVa2omckASrMNVfmBk6Mt oLunry3AN4XCxOYdDkWCsQjTSUibTXGn8p709vTFBfg2BT4+GibP2N54fWtj pUzPkDRZt9LUCJ3//vD2DmEkQJvR12cGGuWwbAQI218tZkWrAUMn9XBfadP4 tg8wswPofTXjsBP4o5yiwPu/zztKouFdUdYTJK4g0mb9S9xvkN2Ovp6+vADP UKA97PKdoCe4lq1QT5DGDOv9rQPvVd2erp6WANNY5qapslR33vTtDmSlm+HW PMl/PwfvXs8wJECb8Qr8/V+hFeL3LY24TODOVakl0J3/WTR39je+e3p6AgEu N/vLc/3Y6b9KGiwY4WV1sWca4cbTXznQ0ccXEWApsPaVWYMPx0ZIZt98CTrG cfvIJPlEAjyL8vvpb31j0Z02xOSRbmQQbtTUao9awr9AvqenEOBbFNiYl+gI a6nHMHRMmn+XdE9fRoD1S7Le98aiO9abFe1wdP93hx38MtK/1PKWszib+iVZ jdc3+jIf7xULK+rrdRVmVgEWNs8IOsZf5zV6bjzJfrlGN74EuP2vqTL8mKMb vAq5jrlINz4SL/d/7BeVj77Si/TceKYU4BOS8ZMAzyLVMxN14+BLnajnsmS6 ZqZuHHutM/VclkTfzNSNEuA5pKY+p+rGoRc7Vc/FeUSD8frvO5U/yMirnavn soR7Z7JuHHi5k/VcEgnQRgI8i2j3zNaN4653tp7LEuyf6bpx2AVP13NZYh00 XTdKgGchAdqMuuL5ei5LqIcm7MZBlzxhz2WJdNGM3TjmmmfsuSQSoI0EeBaB PpqyG4dc9JQ9l8XvpDm7ccRVz9lzWdxemrMbJcCzkABtBlz2pD2XxeumWbvx +HXP2nNZnH6athsPX/i0PZdkYgE2340rvru+JJfuNq337zoXvKMu3I3tO4qG vjPrwj03GNpTF+7GWoC1QTx46RfuucFMLMBb+b6OoW8sunDPjYZ11YW7sXxj UfvONgnwNEhfXbgb65cGtu/rOHTxF+658Zz3woo/RPnGomVpLOCxq79wz41n SgGWbyy6F9Q1jrR+YN/5gL01Qze+541FM/TcQFB3zd6N/dc/e88lkQAB3R0w fc8lAf01fTdKgGdhd5i6sbcH1HMp0G9VqRt7u0A9F4e8MkvdKAGehtVl6sbe PlDP5TH6TN24dHaCei6PBAjo6QX1XAdtp6kb/yEBnkXTa+rGXzq6QT3Xxcjp +CuR7wf1XBcSoI0EeBYj12ReiXRHqOc6Gbks/Upke0I914kECEh2hXqul5Hf TbwSEuBZ3MDfs5PrC/VcPzfzT5HqDPVcPxKgjQR4FjfjL5HrDfXcEW7NH9NQ Ls2tFure8NLdpp3B5zUX8wqwfl/CUr/AKNzQ2POajVv17zS0Aiw3xxsadUaT civ+mYfynVntC2NkAc/itvv/RAx7Z9Z0PTeaaQXIfMDfsmBDI89qSm6v/81E +c4sCfCTwO9qXprynVlW0iXWJ/P13HDmFOCTo+/MmrfnxnFTN5pIgGdxUzea RHpFPTeCg7+VcVkC3aKeG4EECPD7RT13nOi8+4RIgG/lZv1ogdjjdo56rwdT eVKihdcp6rQc2ObJFJpIgIPwhtvmJ9TEHadX1GkeYUdPArTh3aJOg2RDDAkQ QPtFndbSGdxKgAAJ0GH35bYjaRUJEME6Rp223AOIoxk9pQQJpF9m7zLlks9A AmyR8s7k4KLBS1Epj37BWowC9upM3Y3mz/iXa8QQUK9O0dtkuJUAT2JOAfqO XvkFaw3B7wP060W7OzF/5nzBWoziZj7f39/b1RfyO+bPJMBR8FdmmT17gd5+ vJugf/6Mf8FahHGeZtsJf/tZvZ1jCT3/C9YiijecXM0CDs0lS3jHcSM6o4+/ sts1i/E3cSO6rxdgqzxNZPwhel6Z9SX3C9k8BbF/iUBE19ydP3+7+GjrvClW nIsf0X2RAEOOnvemWPEREq/M+oP3KxNiaCLj6/AD43OpQoikCZMP+H24gfHJ 3I6kVTSR8X38KQEeTuhpIuP7cGdHTjmJwblkCe+LcCeI33rwvfIaDUpHU3Az /3z3QVub166XkACn4FwB4tG2FqBcuFm4GX+94Sien9esntAQPAu35o+xzcci DKWRp+VdAswFt3UaWUPwPNyqfw+30yOfMo28LLKAE3Er/jnQzhHD1eSPJcB5 uO3+37X/wCFTupuQfgHKWRMjuL3+F96jlp4WxIsD3JaoAMFwq32QbWIAAADH SURBVMVQ4ggBAXJHTwIUh7hhAXYsiH+UCBHG+qHb7gXxWksg0tQxRE5C1aAr +Ykst21BfId62gXxQjg0K0CP6KaYyVB6UPi8RyLSnYihFXfiw1jpZIlSnEb7 uhjlj8WJSIDio1gC1BgsTsJIldTTGUKcioZg8VEkQPFZlIYRQgghhBBCCCGE EEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGE EEIIIYQQQgghhBBCCCGE+Az/AVyUZ+PdK2rXAAAAAElFTkSuQmCC --0-424385468-1121532012=:8310-- ------------------------------------------------------------- Octave is freely available under the terms of the GNU GPL. Octave's home on the web: http://www.octave.org How to fund new projects: http://www.octave.org/funding.html Subscription information: http://www.octave.org/archive.html -------------------------------------------------------------