function xp=angularpart513to128(x,tipe) % xp=ANGULARPART513TO128(x,tipe) % % INPUT: % % x The three-dimensional array of size 513x513xQ % tipe 'pixel' by numerical integration [fast] % 'wavelets' by wavelet transformation [slow] % % OUTPUT: % % xp The three-dimensional output array of size 128x128xQ % % EXAMPLE: % % angularpart513to128('demo1') % % Written by Ignace Loris (igloris@vub.ac.be) on 23/06/2009 % Last modified by fjsimons-at-alum.mit.edu, 11/07/2011 if ~isstr(x) % Error check on hardwired dimensions [neta,nxi,nrad]=size(x); if neta~=513 || nxi~=513 error('Input should have size 513 in first and second dimension!') end % Prepare output array xp=zeros(128,128,nrad); % Produce the Jacobian and include the sampling interval [J,coordd,dxi,deta]=cubejac(neta,nxi); J=J.*dxi.*deta; % For every depth interval for i=1:nrad % Multiply by Jacobian and sampling interval x(:,:,i)=x(:,:,i).*J; switch tipe case 'pixel' % Use Simpson's rule, with weights [1 4 2 4 1]/3; % First in xi direction y=zeros(128,513); y=(1/3)*x(1:4:509,:,i)+(4/3)*x(2:4:510,:,i)+(2/3)*x(3:4:511,:,i)+... (4/3)*x(4:4:512,:,i)+(1/3)*x(5:4:513,:,i); % Then in eta direction xp(:,:,i)=(1/3)*y(:,1:4:509)+(4/3)*y(:,2:4:510)+(2/3)*y(:,3:4:511)+... (4/3)*y(:,4:4:512)+(1/3)*y(:,5:4:513); % Here is how FJS would have done it... Numerical Recipes 4.1.4 vs 4.1.13 % tst=randi(128); difer([simpson(1:5,x(1:5,tst,i))-y(1,tst)]/y(1,tst)) case 'wavelets' % Change the dimensions from 513x513xQ to 512x512xQ by averaging temp=[x(1:512,1:512,i)+x(2:513,1:512,i)+x(1:512,2:513,i)+x(2:513,2:513,i)]/4; % Compute two levels of the preconditioned inverse transpose wavelet % transform tempw=angularD4WT(temp,[2 2],[1 1],'inversetranspose'); % Throw away the wavelet subbands xp(:,:,i)=tempw(1:128,1:128); otherwise error('The allowable options are ''pixel'' or ''wavelets''') end end elseif strcmp(x,'demo1') % Surface area test kernel=ones(513,513,1); f=angularpart513to128(kernel,'pixel'); err=abs(sum(sum(f))-4*pi/6)/(4*pi/6); disp(['kernel=constant. Model=constant. Integral err is:' num2str(err)]) % test integral of sin(xi)^2*J(xi,eta) x=ones(1,513)'*(sin(linspace(-pi/4,pi/4,513)).^2); kernel(:,:,1)=x; f=angularpart513to128(kernel,'pixel'); int=sum(sum(f)); realint=-sqrt(3)+2*pi/3; err=abs(int-realint)/realint; disp(['kernel=sin(xi)^2. Model=constant. Integral err is:' num2str(err)]) % test integral of xi^2*J(xi,eta) x=ones(1,513)'*(linspace(-pi/4,pi/4,513).^2); kernel(:,:,1)=x; f=angularpart513to128(kernel,'pixel'); int=sum(sum(f)); realint=0.40882059352406189314; err=abs(int-realint)/realint; disp(['kernel=xi^2. Model=constant. Integral err is:' num2str(err)]) % test integral of sin(eta)^2*cos(xi)^2*J(xi,eta) x=exp(linspace(-pi/4,pi/4,513))'*ones(1,513); kernel(:,:,1)=x; f=angularpart513to128(kernel,'pixel'); int=sum(sum(f)); realint=2.3050494965745932318; err=abs(int-realint)/realint; disp(['kernel=sin(eta)^2*cos(xi)^2. Model=constant. Integral err is:' num2str(err)]) elseif strcmp(x,'demo2') % Surface area test % precompute the coarse wavelet representation of the constant model: % constant model modelfull=ones(512,512,1); % reduction temp=angularD4WT(modelfull,[2 2],[1 1],'forward'); modelsmall=temp(1:128,1:128,:); % test integral of constant kernel * constant model kernel=ones(513,513,1); f=angularpart513to128(kernel,'wavelets'); int=sum(sum(f.*modelsmall)); err=abs(int-4*pi/6)/(4*pi/6); disp(['kernel=constant. Constant model. Integral err is:' num2str(err)]) % test integral of sin(xi)^2*J(xi,eta) x=ones(1,513)'*(sin(linspace(-pi/4,pi/4,513)).^2); kernel(:,:,1)=x; f=angularpart513to128(kernel,'wavelets'); int=sum(sum(f.*modelsmall)); realint=-sqrt(3)+2*pi/3; err=abs(int-realint)/realint; disp(['kernel=sin(xi)^2. Constant model. Integral err is:' num2str(err)]) % test integral of xi^2*J(xi,eta) x=ones(1,513)'*(linspace(-pi/4,pi/4,513).^2); kernel(:,:,1)=x; f=angularpart513to128(kernel,'wavelets'); int=sum(sum(f.*modelsmall)); realint=0.40882059352406189314; err=abs(int-realint)/realint; disp(['kernel=xi^2. Constant model. Integral err is:' num2str(err)]) % test integral of sin(eta)^2*cos(xi)^2*J(xi,eta) x=exp(linspace(-pi/4,pi/4,513))'*ones(1,513); kernel(:,:,1)=x; f=angularpart513to128(kernel,'wavelets'); int=sum(sum(f.*modelsmall)); realint=2.3050494965745932318; err=abs(int-realint)/realint; disp(['kernel=sin(eta)^2*cos(xi)^2. Constant model. Integral err is:' num2str(err)]) end