Wayne,
You have to include the revision numbers in the commit message otherwise
it's difficult to do another merge later should the need arise. It's
dangerous to re-merge already merged changes, so you need to know the last
change merged in during the last merge.
Paul
On 5/19/07, waynezacmp@... <
waynezacmp@...> wrote:
>
> Revision: 1049
> http://svn.sourceforge.net/aqsis/?rev=1049&view=rev
> Author: waynezacmp
> Date: 2007-05-18 19:48:41 -0700 (Fri, 18 May 2007)
>
> Log Message:
> -----------
> Fix build error by doing an svn merge with the trunk. Scons builds
> successfully now
>
> Modified Paths:
> --------------
> branches/soc2007-dsm/aqsis/renderer/render/bucket.cpp
> branches/soc2007-dsm/aqsis/shadercompiler/codegenvm/vmoutput.cpp
> branches/soc2007-dsm/aqsis/shadercompiler/shadervm/shadervm1.cpp
> branches/soc2007-dsm/aqsis/shadercompiler/slparse/parser.yy
>
> Added Paths:
> -----------
> branches/soc2007-dsm/aqsis/testing/prototypes/
> branches/soc2007-dsm/aqsis/testing/prototypes/texfilt/
>
> branches/soc2007-dsm/aqsis/testing/prototypes/texfilt/downsample_test.m
>
> branches/soc2007-dsm/aqsis/testing/prototypes/texfilt/downsample_test.py
> branches/soc2007-dsm/aqsis/testing/prototypes/texfilt/texfilt.py
>
> Removed Paths:
> -------------
> branches/soc2007-dsm/aqsis/testing/misc/
> branches/soc2007-dsm/aqsis/testing/prototypes/texfilt/
>
> branches/soc2007-dsm/aqsis/testing/prototypes/texfilt/downsample_test.m
>
> branches/soc2007-dsm/aqsis/testing/prototypes/texfilt/downsample_test.py
> branches/soc2007-dsm/aqsis/testing/prototypes/texfilt/texfilt.py
>
> Modified: branches/soc2007-dsm/aqsis/renderer/render/bucket.cpp
> ===================================================================
> --- branches/soc2007-dsm/aqsis/renderer/render/bucket.cpp 2007-05-15
> 12:34:19 UTC (rev 1048)
> +++ branches/soc2007-dsm/aqsis/renderer/render/bucket.cpp 2007-05-19
> 02:48:41 UTC (rev 1049)
> @@ -42,8 +42,6 @@
> #include <valarray>
>
>
> -//comment added to bucket.cpp
> -
> START_NAMESPACE( Aqsis )
>
>
>
> Modified: branches/soc2007-dsm/aqsis/shadercompiler/codegenvm/vmoutput.cpp
> ===================================================================
> ---
> branches/soc2007-dsm/aqsis/shadercompiler/codegenvm/vmoutput.cpp 2007-05-15
> 12:34:19 UTC (rev 1048)
> +++
> branches/soc2007-dsm/aqsis/shadercompiler/codegenvm/vmoutput.cpp 2007-05-19
> 02:48:41 UTC (rev 1049)
> @@ -639,25 +639,31 @@
> TqInt iLabelC = m_gcLabels++;
> TqInt iLabelD = m_gcLabels++;
>
> - IqParseNode* pArg = pNode->pChild();
> + IqParseNode* pArgs = pNode->pChild();
> assert( pArg != 0 );
> - IqParseNode* pHitStmt = pArg->pNextSibling();
> + IqParseNode* pHitStmt = pArgs->pNextSibling();
> assert( pHitStmt != 0 );
> IqParseNode* pNoHitStmt = pHitStmt->pNextSibling();
>
> - IqParseNode* pSamples = pArg->pChild();
> + // samples is the 5th argument to the gather function, we need
> that to
> + // initialise the init_gather call.
> + IqParseNode* pSamples = pArgs->pChild();
> TqInt iArg = 5;
> while(--iArg > 0)
> {
> pSamples = pSamples->pNextSibling();
> assert(pSamples);
> }
> + pSamples->Accept(*this);
>
> - pSamples->Accept(*this);
> m_slxFile << "\tinit_gather" << std::endl;
>
> m_slxFile << ":" << iLabelA << std::endl; // loop
> back label
> m_slxFile << "\tS_CLEAR" << std::endl; // clear
> current state
> +
> + // Output the arguments in reverse order so that the gather()
> function
> + // can pop them in the expected order.
> + IqParseNode* pArg = pArgs->pChild();
> while ( pArg->pNextSibling() != 0 )
> pArg = pArg->pNextSibling();
> iArg = 0;
> @@ -668,6 +674,7 @@
> pArg->Accept( *this );
> pArg = pArg->pPrevSibling();
> }
> + // Now output the count of 'additional' arguments.
> iArg -= 5;
> CqParseNodeFloatConst C( static_cast<TqFloat>( abs( iArg ) ) );
> C.Accept( *this );
>
> Modified: branches/soc2007-dsm/aqsis/shadercompiler/shadervm/shadervm1.cpp
> ===================================================================
> ---
> branches/soc2007-dsm/aqsis/shadercompiler/shadervm/shadervm1.cpp 2007-05-15
> 12:34:19 UTC (rev 1048)
> +++
> branches/soc2007-dsm/aqsis/shadercompiler/shadervm/shadervm1.cpp 2007-05-19
> 02:48:41 UTC (rev 1049)
> @@ -451,7 +451,6 @@
> SqLabel lab = ReadNext().m_Label;
> if ( m_pEnv->RunningState().Count() == 0 )
> {
> - std::cout << m_PO << ", " << lab.m_Offset << std::endl;
> m_PO = lab.m_Offset;
> m_PC = lab.m_pAddress;
> }
>
> Modified: branches/soc2007-dsm/aqsis/shadercompiler/slparse/parser.yy
> ===================================================================
> --- branches/soc2007-dsm/aqsis/shadercompiler/slparse/parser.yy 2007-05-15
> 12:34:19 UTC (rev 1048)
> +++ branches/soc2007-dsm/aqsis/shadercompiler/slparse/parser.yy 2007-05-19
> 02:48:41 UTC (rev 1049)
> @@ -1,6 +1,6 @@
> /* -------------- declaration section -------------- */
>
> -%expect 24
> +%expect 26
>
> %{
> #ifdef WIN32
> @@ -1771,6 +1771,7 @@
> }
> ;
>
> +
> unresolvedcall
> : IDENTIFIER '(' proc_arguments ')'
> {
> @@ -1891,12 +1892,6 @@
> $$=new
> CqParseNodeFunctionCall(func);
>
> $$->SetPos(ParseLineNumber,
> ParseStreamName.c_str());
> }
> - | OCCLUSION {
> -
> std::vector<SqFuncRef> func;
> -
> CqFuncDef::FindFunction("occlusion", func);
> - $$=new
> CqParseNodeFunctionCall(func);
> -
> $$->SetPos(ParseLineNumber,ParseStreamName.c_str());
> - }
> ;
>
> texture_filename
>
> Copied: branches/soc2007-dsm/aqsis/testing/prototypes (from rev 1048,
> trunk/aqsis/testing/prototypes)
>
> Copied: branches/soc2007-dsm/aqsis/testing/prototypes/texfilt (from rev
> 1048, trunk/aqsis/testing/prototypes/texfilt)
>
> Deleted:
> branches/soc2007-dsm/aqsis/testing/prototypes/texfilt/downsample_test.m
> ===================================================================
> --- trunk/aqsis/testing/prototypes/texfilt/downsample_test.m 2007-05-15
> 12:34:19 UTC (rev 1048)
> +++
> branches/soc2007-dsm/aqsis/testing/prototypes/texfilt/downsample_test.m
> 2007-05-19 02:48:41 UTC (rev 1049)
> @@ -1,379 +0,0 @@
> -% Matlab script to test image filtering for mipmap creation in Aqsis
> -% Copyright (C) 2007, Christopher J. Foster
> -%
> -% This program is free software; you can redistribute it and/or
> -% modify it under the terms of the GNU General Public License
> -% as published by the Free Software Foundation; either version 2
> -% of the License, or (at your option) any later version.
> -%
> -% This program is distributed in the hope that it will be useful,
> -% but WITHOUT ANY WARRANTY; without even the implied warranty of
> -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
> -% GNU General Public License for more details.
> -%
> -% You should have received a copy of the GNU General Public License
> -% along with this program; if not, write to the Free Software
> -% Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
> MA 02110-1301, USA.
> -
> -%----------------------------------------------------------------------
> -function img_filter()
> -
> -%im = double(imread('lena_std.tif'))/255;
> -%im = double(imread('kitten_tst.png'))/255;
> -im = getGridImg(128,21);
> -
> -%investigateWindowTypes(im, @sincKer);
> -%plotWindowFuncs(10);
> -compareOddEven(im);
> -
> -%windowFunc = makeWinWithoutZeros(@lanczos);
> -%windowFunc = @(N) lanczos(N,1.6);
> -%kerFunc = @sincKer;
> -%
> -% Investigate factor of 2 vs factor of 4 downsampling
> -%mipmap = genMipmapIter(im, 7, 1, windowFunc, kerFunc);
> -%mipmap = genMipmapIterDownsamp4(im, 7, 1, windowFunc, kerFunc);
> -%
> -%imwrite(cellarrayToImage(mipmap(2:end),2), 'mipmap.png','png');
> -
> -%image(cellarrayToImage(mipmap,1));
> -%axis image;
> -%axis off;
> -
> -
> -%----------------------------------------------------------------------
> -% High level driver functions.
> -%----------------------------------------------------------------------
> -function compareOddEven(im)
> -% Compare the difference between sampling an image with an odd or
> even-sized
> -% filter.
> -
> -kerFunc = @sincKer;
> -windowFunc = makeWinWithoutZeros(@lanczos);
> -
> -% Use an odd-width filter
> -mipmap = genMipmapIter(im, 7, 0, windowFunc, kerFunc);
> -imwrite(cellarrayToImage(mipmap,2), 'mipmap_odd_width.png','png');
> -
> -%pause(4);
> -% Use an even-width filter
> -mipmap = genMipmapIter(im, 8, 0, windowFunc, kerFunc);
> -imwrite(cellarrayToImage(mipmap,2), 'mipmap_even_width.png', 'png');
> -
> -
> -%----------------------------------------------------------------------
> -function investigateWindowTypes(im, kerFunc)
> -% Grab the level-three mipmaps arising from using a bunch of windows with
> -% a fixed filter kernel.
> -
> -% Generate a version using the boxcar kernel [1,1] for reference.
> -mipmap = genMipmapIter(im, 1, 1, @rectwin, @boxKer);
> -miplevel3{1} = mipmap{3};
> -
> -winFuncs = {@lanczos, @hamming, @hanning, @blackman, @rectwin, @kaiser};
> -% Generate mipmaps using the other window functions
> -for ii=1:length(winFuncs)
> - mipmap = genMipmapIter(im, 7, 1, winFuncs{ii}, kerFunc);
> - miplevel3{ii+1} = mipmap{3};
> -end
> -
> -imwrite(cellarrayToImage(miplevel3,2), 'mipmap_lvl3_all_wins.png','png');
> -
> -
> -%----------------------------------------------------------------------
> -function plotWindowFuncs(N)
> -% Plot various window functions.
> -winFuncs = {@lanczos, @hamming, @hanning, @blackman, @rectwin, @kaiser};
> -winFuncNames = {'lanczos,' 'hamming,' 'hanning,' 'blackman', 'boxcar',
> 'kaiser'};
> -clf();
> -hold on;
> -for ii=1:length(winFuncs)
> - plot(winFuncs{ii}(N), rnd_ln_style(ii));
> -end
> -legend(winFuncNames);
> -axis([1,N,0,1.05]);
> -set(gca(),'box','on');
> -
> -
> -%----------------------------------------------------------------------
> -% Functions for creating/manipulating mipmaps
> -%----------------------------------------------------------------------
> -function mipmap = genMipmapIter(im, width, selectWidth, windowFunc,
> kerFunc)
> -% Generate a mipmap by downsampling each image in turn to get the next
> level.
> -%
> -% im - input image
> -% width - desired "width" of filter kernel. the kernel will have width+1
> points.
> -% selectWidth - whether the function should try to adjust the filter
> width for
> -% even/odd sized images. To avoid loosing information,
> even
> -% filters should go with even sized images and vice versa
> -% windowFunc - windowing function to use on the filter.
> -%
> -
> -mipmap = {};
> -
> -scale = 2;
> -
> -ker = kerFunc(width, scale);
> -ker = applyWindow(ker, windowFunc);
> -
> -ii = 1;
> -while(size(im,1) > 1)
> - mipmap{ii} = im;
> - ii = ii + 1;
> -
> - if(selectWidth)
> - if(abs(mod(size(im,1) + width + 1, 2)) < eps(10))
> - ker = kerFunc(width,scale);
> - else
> - ker = kerFunc(width+1,scale);
> - end
> - ker = applyWindow(ker, windowFunc);
> - end
> - disp(sprintf('side length = %d, number of points in filter = %d', ...
> - size(im,1), length(ker)));
> -
> - im = imDownsamp(im, ker, scale);
> -end
> -
> -
> -%----------------------------------------------------------------------
> -function mipmap = genMipmapIterDownsamp4(im, width, selectWidth,
> windowFunc, kerFunc)
> -% Same as genMipmapIter, except tries to downsample by a factor of 4
> rather
> -% than two when possible.
> -%
> -% (Beware: ugly hack)
> -
> -mipmap = {};
> -
> -scale = 4;
> -
> -scales = scale*ones(1,20);
> -scales(1) = 2;
> -widths = 2*width*ones(1,20);
> -widths(1) = width;
> -
> -ker = kerFunc(width, scale);
> -ker = applyWindow(ker, windowFunc);
> -prevIm = im;
> -
> -ii = 1;
> -while(size(im,1) > 1)
> - scale = scales(ii);
> - width = widths(ii);
> -
> - mipmap{ii} = im;
> - ii = ii + 1;
> -
> - if(selectWidth)
> - if(abs(mod(size(im,1) + width + 1, 2)) < eps(10))
> - ker = kerFunc(width,scale);
> - else
> - ker = kerFunc(width+1,scale);
> - end
> - ker = applyWindow(ker, windowFunc);
> - end
> - disp(sprintf('side length = %d, number of points in filter = %d', ...
> - size(im,1), length(ker)));
> -
> - im = imDownsamp(mipmap{max(ii-2,1)}, ker, scale);
> -end
> -
> -
> -%----------------------------------------------------------------------
> -function mipmap = genmipmap(im)
> -% Generate a mipmap by applying succesively larger filter kernels to the
> -% original image (Note: impractically slow!)
> -%
> -% Currently broken.
> -
> -mipmap = {};
> -
> -imSmall = im;
> -
> -scale = 2;
> -currScale = 2;
> -xpos = 1;
> -for ii = 1:3
> - mipmap(ii) = imSmall;
> - xpos = xpos + size(imSmall,1);
> - ker = sincKer(2*currScale,currScale);
> -% window = sin(linspace(0,pi,length(ker)));
> -% window = window(:) * window(:).';
> -% ker = ker .* window;
> - disp(ii);
> -
> - imSmall = imDownsamp(im, ker, currScale);
> -
> - currScale = currScale*scale;
> -end
> -
> -
> -%----------------------------------------------------------------------
> -function mipmapImg = cellarrayToImage(mipmap, dim)
> -% "Flatten" mipmap levels into a single image.
> -%
> -% mipmap - cell array of scaled images
> -% dim - dimension number to lay the images out along
> -
> -len = 0;
> -for ii=1:length(mipmap)
> - len = len + size(mipmap{ii},dim);
> - %disp(size(mipmap{ii},dim));
> -end
> -
> -mipSize = size(mipmap{1});
> -mipSize(dim) = len;
> -
> -mipmapImg = ones(mipSize)*0.7;
> -pos = 1;
> -for ii = 1:length(mipmap)
> - if dim == 1
> - mipmapImg(pos:pos+size(mipmap{ii},1)-1, 1:size(mipmap{ii},2), :) =
> mipmap{ii};
> - else
> - mipmapImg(1:size(mipmap{ii},1), pos:pos+size(mipmap{ii},2)-1, :) =
> mipmap{ii};
> - end
> - pos = pos + size(mipmap{ii},1);
> -end
> -
> -
> -%----------------------------------------------------------------------
> -%----------------------------------------------------------------------
> -% Windowing functions
> -%----------------------------------------------------------------------
> -function win = lanczos(N, tau)
> -% return N-point Lanczos window with width tau
> -if nargin == 1
> - tau = 1;
> -end
> -win = sinc(linspace(-1,1,N)/tau);
> -
> -
> -%----------------------------------------------------------------------
> -function win = makeWinWithoutZeros(winFunc)
> -% Return a window function constructed from the given window function,
> but
> -% without the two edge points (see natWin).
> -win = @(N, varargin) natWin(N, winFunc, varargin{:});
> -
> -
> -function win = natWin(N, winFunc, varargin)
> -% Windowing functions often include zeros at the edges. This is a waste
> of
> -% computational effort, since convolution with zeros doesn't achieve
> anything.
> -%
> -% This function removes the two end points in any window.
> -win = winFunc(N+2, varargin{:});
> -win = win(2:end-1);
> -
> -%----------------------------------------------------------------------
> -function kerWin = applyWindow(ker, windowFunc)
> -% Window a filter kernel with the given windowing function.
> -
> -window = windowFunc(length(ker));
> -window = window(:) * window(:).';
> -kerWin = ker .* window;
> -%kerWin = ker
> -
> -
> -%----------------------------------------------------------------------
> -%----------------------------------------------------------------------
> -% Filter Kernels
> -%----------------------------------------------------------------------
> -function k = sincKer(width, scale)
> -% Return a sinc filter kernel with specified width and scale
> -%
> -% width - total number of points in the filter (not radius!)
> -% scale - zeros will be centered around 0 with spacing "scale"
> -
> -x = (0:width)-width*0.5;
> -k = sinc(x/scale)' * sinc(x/scale);
> -
> -
> -%----------------------------------------------------------------------
> -function k = boxKer(width, scale)
> -% Trivial box filter kernel
> -
> -k = ones(width+1);
> -
> -
> -%----------------------------------------------------------------------
> -%----------------------------------------------------------------------
> -% Image manipuation
> -%----------------------------------------------------------------------
> -function res = imConv(im, ker)
> -% Convolve an image with the given filter kernel
> -
> -res = zeros(size(im));
> -normalisation = conv2(ones(size(im,1), size(im,2)), ker, 'full');
> -% compute how much of the full convolution to trim from top left and
> bottom right
> -trim_tl = ceil((size(ker)-1)/2);
> -trim_br = floor((size(ker)-1)/2);
> -for ii=1:size(im,3)
> - fullconv = conv2(im(:,:,ii), ker, 'full')./normalisation;
> - res(:,:,ii) = fullconv(1+trim_tl(1):end-trim_br(1),
> 1+trim_tl(2):end-trim_br(2));
> -end
> -
> -
> -%----------------------------------------------------------------------
> -function res = imDownsamp(im, ker, skip)
> -% Downsample an image by
> -% (1) Applying the filter kernel ker at each pixel
> -% (2) Decimate by taking elements 1:skip:end from the result
> -% (3) Clamp the resulting intensities between 0 and 1
> -
> -res = imConv(im,ker);
> -res = res(1+floor((skip-1)/2):skip:end,1+floor((skip-1)/2):skip:end,:);
> -res = clamp(res);
> -
> -
> -%----------------------------------------------------------------------
> -function trilinear()
> -% Trilinear interpolation (not implemented...)
> -blah();
> -
> -
> -%----------------------------------------------------------------------
> -function res = clamp(im)
> -% Clamp an image between zero and one.
> -res = min(max(im,0),1);
> -
> -
> -%----------------------------------------------------------------------
> -function im = getGridImg(N, gridSize)
> -% Get a test image containig a one pixel wide grid
> -%
> -% N - side length
> -% gridSize - size of the one-pixel wide black grid in the image.
> -
> -im = zeros(N, N, 3);
> -
> -x = linspace(0,1,N);
> -[xx,yy] = ndgrid(x,x);
> -im(:,:,1) = xx;
> -im(:,:,2) = 0.5;
> -im(:,:,3) = yy;
> -
> -im(1:gridSize:end-gridSize+1,:,:) = 0;
> -im(:,1:gridSize:end-gridSize+1,:) = 0;
> -
> -im(:,end,:) = 0;
> -im(end,:,:) = 0;
> -
> -
> -%----------------------------------------------------------------------
> -% Misc
> -%----------------------------------------------------------------------
> -function lineStyle = rnd_ln_style(idx, colOverride)
> -% lsty = rnd_ln_style(idx, colOverride)
> -%
> -% Return a 'random' line style which may be used to usefully distinguish
> -% between a large number of curves.
> -
> -lineDashes = {'-','--','-.',':'};
> -numDashed = length(lineDashes);
> -lineColors = {'b','k','r'};
> -if nargin >= 2
> - lineColors = {colOverride};
> -end
> -numCols = length(lineColors);
> -
> -idx = idx - 1;
> -lineStyle = [lineColors{mod(idx, numCols)+1}, lineDashes{mod(idx,
> numDashed)+1}];
>
> Copied:
> branches/soc2007-dsm/aqsis/testing/prototypes/texfilt/downsample_test.m
> (from rev 1048, trunk/aqsis/testing/prototypes/texfilt/downsample_test.m)
> ===================================================================
> ---
> branches/soc2007-dsm/aqsis/testing/prototypes/texfilt/downsample_test.m
> (rev 0)
> +++
> branches/soc2007-dsm/aqsis/testing/prototypes/texfilt/downsample_test.m
> 2007-05-19 02:48:41 UTC (rev 1049)
> @@ -0,0 +1,379 @@
> +% Matlab script to test image filtering for mipmap creation in Aqsis
> +% Copyright (C) 2007, Christopher J. Foster
> +%
> +% This program is free software; you can redistribute it and/or
> +% modify it under the terms of the GNU General Public License
> +% as published by the Free Software Foundation; either version 2
> +% of the License, or (at your option) any later version.
> +%
> +% This program is distributed in the hope that it will be useful,
> +% but WITHOUT ANY WARRANTY; without even the implied warranty of
> +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
> +% GNU General Public License for more details.
> +%
> +% You should have received a copy of the GNU General Public License
> +% along with this program; if not, write to the Free Software
> +% Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
> MA 02110-1301, USA.
> +
> +%----------------------------------------------------------------------
> +function img_filter()
> +
> +%im = double(imread('lena_std.tif'))/255;
> +%im = double(imread('kitten_tst.png'))/255;
> +im = getGridImg(128,21);
> +
> +%investigateWindowTypes(im, @sincKer);
> +%plotWindowFuncs(10);
> +compareOddEven(im);
> +
> +%windowFunc = makeWinWithoutZeros(@lanczos);
> +%windowFunc = @(N) lanczos(N,1.6);
> +%kerFunc = @sincKer;
> +%
> +% Investigate factor of 2 vs factor of 4 downsampling
> +%mipmap = genMipmapIter(im, 7, 1, windowFunc, kerFunc);
> +%mipmap = genMipmapIterDownsamp4(im, 7, 1, windowFunc, kerFunc);
> +%
> +%imwrite(cellarrayToImage(mipmap(2:end),2), 'mipmap.png','png');
> +
> +%image(cellarrayToImage(mipmap,1));
> +%axis image;
> +%axis off;
> +
> +
> +%----------------------------------------------------------------------
> +% High level driver functions.
> +%----------------------------------------------------------------------
> +function compareOddEven(im)
> +% Compare the difference between sampling an image with an odd or
> even-sized
> +% filter.
> +
> +kerFunc = @sincKer;
> +windowFunc = makeWinWithoutZeros(@lanczos);
> +
> +% Use an odd-width filter
> +mipmap = genMipmapIter(im, 7, 0, windowFunc, kerFunc);
> +imwrite(cellarrayToImage(mipmap,2), 'mipmap_odd_width.png','png');
> +
> +%pause(4);
> +% Use an even-width filter
> +mipmap = genMipmapIter(im, 8, 0, windowFunc, kerFunc);
> +imwrite(cellarrayToImage(mipmap,2), 'mipmap_even_width.png', 'png');
> +
> +
> +%----------------------------------------------------------------------
> +function investigateWindowTypes(im, kerFunc)
> +% Grab the level-three mipmaps arising from using a bunch of windows with
> +% a fixed filter kernel.
> +
> +% Generate a version using the boxcar kernel [1,1] for reference.
> +mipmap = genMipmapIter(im, 1, 1, @rectwin, @boxKer);
> +miplevel3{1} = mipmap{3};
> +
> +winFuncs = {@lanczos, @hamming, @hanning, @blackman, @rectwin, @kaiser};
> +% Generate mipmaps using the other window functions
> +for ii=1:length(winFuncs)
> + mipmap = genMipmapIter(im, 7, 1, winFuncs{ii}, kerFunc);
> + miplevel3{ii+1} = mipmap{3};
> +end
> +
> +imwrite(cellarrayToImage(miplevel3,2), 'mipmap_lvl3_all_wins.png','png');
> +
> +
> +%----------------------------------------------------------------------
> +function plotWindowFuncs(N)
> +% Plot various window functions.
> +winFuncs = {@lanczos, @hamming, @hanning, @blackman, @rectwin, @kaiser};
> +winFuncNames = {'lanczos,' 'hamming,' 'hanning,' 'blackman', 'boxcar',
> 'kaiser'};
> +clf();
> +hold on;
> +for ii=1:length(winFuncs)
> + plot(winFuncs{ii}(N), rnd_ln_style(ii));
> +end
> +legend(winFuncNames);
> +axis([1,N,0,1.05]);
> +set(gca(),'box','on');
> +
> +
> +%----------------------------------------------------------------------
> +% Functions for creating/manipulating mipmaps
> +%----------------------------------------------------------------------
> +function mipmap = genMipmapIter(im, width, selectWidth, windowFunc,
> kerFunc)
> +% Generate a mipmap by downsampling each image in turn to get the next
> level.
> +%
> +% im - input image
> +% width - desired "width" of filter kernel. the kernel will have width+1
> points.
> +% selectWidth - whether the function should try to adjust the filter
> width for
> +% even/odd sized images. To avoid loosing information,
> even
> +% filters should go with even sized images and vice versa
> +% windowFunc - windowing function to use on the filter.
> +%
> +
> +mipmap = {};
> +
> +scale = 2;
> +
> +ker = kerFunc(width, scale);
> +ker = applyWindow(ker, windowFunc);
> +
> +ii = 1;
> +while(size(im,1) > 1)
> + mipmap{ii} = im;
> + ii = ii + 1;
> +
> + if(selectWidth)
> + if(abs(mod(size(im,1) + width + 1, 2)) < eps(10))
> + ker = kerFunc(width,scale);
> + else
> + ker = kerFunc(width+1,scale);
> + end
> + ker = applyWindow(ker, windowFunc);
> + end
> + disp(sprintf('side length = %d, number of points in filter = %d', ...
> + size(im,1), length(ker)));
> +
> + im = imDownsamp(im, ker, scale);
> +end
> +
> +
> +%----------------------------------------------------------------------
> +function mipmap = genMipmapIterDownsamp4(im, width, selectWidth,
> windowFunc, kerFunc)
> +% Same as genMipmapIter, except tries to downsample by a factor of 4
> rather
> +% than two when possible.
> +%
> +% (Beware: ugly hack)
> +
> +mipmap = {};
> +
> +scale = 4;
> +
> +scales = scale*ones(1,20);
> +scales(1) = 2;
> +widths = 2*width*ones(1,20);
> +widths(1) = width;
> +
> +ker = kerFunc(width, scale);
> +ker = applyWindow(ker, windowFunc);
> +prevIm = im;
> +
> +ii = 1;
> +while(size(im,1) > 1)
> + scale = scales(ii);
> + width = widths(ii);
> +
> + mipmap{ii} = im;
> + ii = ii + 1;
> +
> + if(selectWidth)
> + if(abs(mod(size(im,1) + width + 1, 2)) < eps(10))
> + ker = kerFunc(width,scale);
> + else
> + ker = kerFunc(width+1,scale);
> + end
> + ker = applyWindow(ker, windowFunc);
> + end
> + disp(sprintf('side length = %d, number of points in filter = %d', ...
> + size(im,1), length(ker)));
> +
> + im = imDownsamp(mipmap{max(ii-2,1)}, ker, scale);
> +end
> +
> +
> +%----------------------------------------------------------------------
> +function mipmap = genmipmap(im)
> +% Generate a mipmap by applying succesively larger filter kernels to the
> +% original image (Note: impractically slow!)
> +%
> +% Currently broken.
> +
> +mipmap = {};
> +
> +imSmall = im;
> +
> +scale = 2;
> +currScale = 2;
> +xpos = 1;
> +for ii = 1:3
> + mipmap(ii) = imSmall;
> + xpos = xpos + size(imSmall,1);
> + ker = sincKer(2*currScale,currScale);
> +% window = sin(linspace(0,pi,length(ker)));
> +% window = window(:) * window(:).';
> +% ker = ker .* window;
> + disp(ii);
> +
> + imSmall = imDownsamp(im, ker, currScale);
> +
> + currScale = currScale*scale;
> +end
> +
> +
> +%----------------------------------------------------------------------
> +function mipmapImg = cellarrayToImage(mipmap, dim)
> +% "Flatten" mipmap levels into a single image.
> +%
> +% mipmap - cell array of scaled images
> +% dim - dimension number to lay the images out along
> +
> +len = 0;
> +for ii=1:length(mipmap)
> + len = len + size(mipmap{ii},dim);
> + %disp(size(mipmap{ii},dim));
> +end
> +
> +mipSize = size(mipmap{1});
> +mipSize(dim) = len;
> +
> +mipmapImg = ones(mipSize)*0.7;
> +pos = 1;
> +for ii = 1:length(mipmap)
> + if dim == 1
> + mipmapImg(pos:pos+size(mipmap{ii},1)-1, 1:size(mipmap{ii},2), :) =
> mipmap{ii};
> + else
> + mipmapImg(1:size(mipmap{ii},1), pos:pos+size(mipmap{ii},2)-1, :) =
> mipmap{ii};
> + end
> + pos = pos + size(mipmap{ii},1);
> +end
> +
> +
> +%----------------------------------------------------------------------
> +%----------------------------------------------------------------------
> +% Windowing functions
> +%----------------------------------------------------------------------
> +function win = lanczos(N, tau)
> +% return N-point Lanczos window with width tau
> +if nargin == 1
> + tau = 1;
> +end
> +win = sinc(linspace(-1,1,N)/tau);
> +
> +
> +%----------------------------------------------------------------------
> +function win = makeWinWithoutZeros(winFunc)
> +% Return a window function constructed from the given window function,
> but
> +% without the two edge points (see natWin).
> +win = @(N, varargin) natWin(N, winFunc, varargin{:});
> +
> +
> +function win = natWin(N, winFunc, varargin)
> +% Windowing functions often include zeros at the edges. This is a waste
> of
> +% computational effort, since convolution with zeros doesn't achieve
> anything.
> +%
> +% This function removes the two end points in any window.
> +win = winFunc(N+2, varargin{:});
> +win = win(2:end-1);
> +
> +%----------------------------------------------------------------------
> +function kerWin = applyWindow(ker, windowFunc)
> +% Window a filter kernel with the given windowing function.
> +
> +window = windowFunc(length(ker));
> +window = window(:) * window(:).';
> +kerWin = ker .* window;
> +%kerWin = ker
> +
> +
> +%----------------------------------------------------------------------
> +%----------------------------------------------------------------------
> +% Filter Kernels
> +%----------------------------------------------------------------------
> +function k = sincKer(width, scale)
> +% Return a sinc filter kernel with specified width and scale
> +%
> +% width - total number of points in the filter (not radius!)
> +% scale - zeros will be centered around 0 with spacing "scale"
> +
> +x = (0:width)-width*0.5;
> +k = sinc(x/scale)' * sinc(x/scale);
> +
> +
> +%----------------------------------------------------------------------
> +function k = boxKer(width, scale)
> +% Trivial box filter kernel
> +
> +k = ones(width+1);
> +
> +
> +%----------------------------------------------------------------------
> +%----------------------------------------------------------------------
> +% Image manipuation
> +%----------------------------------------------------------------------
> +function res = imConv(im, ker)
> +% Convolve an image with the given filter kernel
> +
> +res = zeros(size(im));
> +normalisation = conv2(ones(size(im,1), size(im,2)), ker, 'full');
> +% compute how much of the full convolution to trim from top left and
> bottom right
> +trim_tl = ceil((size(ker)-1)/2);
> +trim_br = floor((size(ker)-1)/2);
> +for ii=1:size(im,3)
> + fullconv = conv2(im(:,:,ii), ker, 'full')./normalisation;
> + res(:,:,ii) = fullconv(1+trim_tl(1):end-trim_br(1),
> 1+trim_tl(2):end-trim_br(2));
> +end
> +
> +
> +%----------------------------------------------------------------------
> +function res = imDownsamp(im, ker, skip)
> +% Downsample an image by
> +% (1) Applying the filter kernel ker at each pixel
> +% (2) Decimate by taking elements 1:skip:end from the result
> +% (3) Clamp the resulting intensities between 0 and 1
> +
> +res = imConv(im,ker);
> +res = res(1+floor((skip-1)/2):skip:end,1+floor((skip-1)/2):skip:end,:);
> +res = clamp(res);
> +
> +
> +%----------------------------------------------------------------------
> +function trilinear()
> +% Trilinear interpolation (not implemented...)
> +blah();
> +
> +
> +%----------------------------------------------------------------------
> +function res = clamp(im)
> +% Clamp an image between zero and one.
> +res = min(max(im,0),1);
> +
> +
> +%----------------------------------------------------------------------
> +function im = getGridImg(N, gridSize)
> +% Get a test image containig a one pixel wide grid
> +%
> +% N - side length
> +% gridSize - size of the one-pixel wide black grid in the image.
> +
> +im = zeros(N, N, 3);
> +
> +x = linspace(0,1,N);
> +[xx,yy] = ndgrid(x,x);
> +im(:,:,1) = xx;
> +im(:,:,2) = 0.5;
> +im(:,:,3) = yy;
> +
> +im(1:gridSize:end-gridSize+1,:,:) = 0;
> +im(:,1:gridSize:end-gridSize+1,:) = 0;
> +
> +im(:,end,:) = 0;
> +im(end,:,:) = 0;
> +
> +
> +%----------------------------------------------------------------------
> +% Misc
> +%----------------------------------------------------------------------
> +function lineStyle = rnd_ln_style(idx, colOverride)
> +% lsty = rnd_ln_style(idx, colOverride)
> +%
> +% Return a 'random' line style which may be used to usefully distinguish
> +% between a large number of curves.
> +
> +lineDashes = {'-','--','-.',':'};
> +numDashed = length(lineDashes);
> +lineColors = {'b','k','r'};
> +if nargin >= 2
> + lineColors = {colOverride};
> +end
> +numCols = length(lineColors);
> +
> +idx = idx - 1;
> +lineStyle = [lineColors{mod(idx, numCols)+1}, lineDashes{mod(idx,
> numDashed)+1}];
>
> Deleted:
> branches/soc2007-dsm/aqsis/testing/prototypes/texfilt/downsample_test.py
> ===================================================================
> --- trunk/aqsis/testing/prototypes/texfilt/downsample_test.py 2007-05-15
> 12:34:19 UTC (rev 1048)
> +++
> branches/soc2007-dsm/aqsis/testing/prototypes/texfilt/downsample_test.py 2007-05-19
> 02:48:41 UTC (rev 1049)
> @@ -1,337 +0,0 @@
> -#!/usr/bin/python
> -#
> -# Scientific python script to test image filtering for mipmap creation in
> Aqsis
> -# Copyright (C) 2007, Christopher J. Foster
> -#
> -# This program is free software; you can redistribute it and/or
> -# modify it under the terms of the GNU General Public License
> -# as published by the Free Software Foundation; either version 2
> -# of the License, or (at your option) any later version.
> -#
> -# This program is distributed in the hope that it will be useful,
> -# but WITHOUT ANY WARRANTY; without even the implied warranty of
> -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
> -# GNU General Public License for more details.
> -#
> -# You should have received a copy of the GNU General Public License
> -# along with this program; if not, write to the Free Software
> -# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
> MA 02110-1301, USA.
> -
> -# This script investigates some of the factors important when
> downsampling an
> -# image:
> -#
> -# * Filter types
> -# * Filter windows
> -# * Correctly filtering odd and even-sized images
> -
> -from __future__ import division
> -
> -import matplotlib.pylab as pylab
> -
> -import numpy
> -from numpy import r_, size, zeros, ones, ceil, floor, array, linspace,
> meshgrid
> -
> -import scipy
> -from scipy.signal import boxcar, convolve2d as conv2
> -from scipy.misc import imsave, imread
> -
> -#----------------------------------------------------------------------
> -#----------------------------------------------------------------------
> -# Windowing functions
> -#----------------------------------------------------------------------
> -def lanczos(N, tau=1):
> - '''
> - return N-point Lanczos window with width tau
> - '''
> - return scipy.sinc(linspace(-1,1,N)/tau)
> -
> -
> -#----------------------------------------------------------------------
> -def makeTrimmedWindow(winFunc):
> - '''
> - Return a window function constructed from the given window
> function, but
> - without the two edge points (see natWin).
> - '''
> - return lambda N, *args: winFunc(N+2, *args)[1:-1]
> -
> -
> -#----------------------------------------------------------------------
> -def applyWindow(ker, windowFunc):
> - '''
> - Window a filter kernel with the given windowing function.
> - '''
> - window = windowFunc(len(ker))
> - window = numpy.outer(window, window)
> - return ker * window
> -
> -
> -#----------------------------------------------------------------------
> -#----------------------------------------------------------------------
> -# Filter Kernels
> -#----------------------------------------------------------------------
> -def sincKer(width, scale):
> - '''
> - Return a sinc filter kernel with specified width and scale
> -
> - width - total number of points in the filter (not radius!)
> - scale - zeros will be centered around 0 with spacing "scale"
> - '''
> - x = r_[0:width+1]-width*0.5
> - k = numpy.outer(scipy.sinc(x/scale), scipy.sinc(x/scale))
> - return k
> -
> -
> -#----------------------------------------------------------------------
> -def boxKer(width, scale):
> - '''
> - Trivial box filter kernel
> - '''
> - return ones(width+1)
> -
> -
> -#----------------------------------------------------------------------
> -# Image manipuation
> -#----------------------------------------------------------------------
> -def seperableConv2(A, k, *args):
> - '''
> - perform a 2D convolution, assuming that the convolution kernel k
> is
> - square-seperable. In this case we may just do two 1D
> convolutions.
> - '''
> - center = array(k.shape)//2
> - k0 = k[center[0]:center[0]+1,:]
> - k1 = k[:,center[1]:center[1]+1]
> - return conv2(conv2(A,k0, *args),k1, *args)
> -
> -
> -#----------------------------------------------------------------------
> -def imConv(im, ker, conv2Func=conv2):
> - '''
> - Convolve an image with the given filter kernel
> - '''
> - res = zeros(im.shape)
> - normalisation = conv2Func(ones(im.shape[0:2],'d'), ker, 'full')
> - # compute how much of the full convolution to trim from top left
> and bottom right
> - trimTl = numpy.int_(ceil((array(ker.shape)-1)/2))
> - trimBr = numpy.int_(floor((array(ker.shape)-1)/2))
> - for ii in range(0,size(im,2)):
> - fullconv = conv2Func(im[:,:,ii], ker, 'full') /
> normalisation
> - res[:,:,ii] = fullconv[trimTl[0]:-trimBr[0],
> trimTl[1]:-trimBr[1]]
> - return res
> -
> -
> -#----------------------------------------------------------------------
> -def imDownsamp(im, ker, skip):
> - '''
> - Downsample an image by
> - (1) Applying the filter kernel ker at each pixel
> - (2) Decimate by taking elements 0::skip from the result
> - (3) Clamp the resulting intensities between 0 and 1
> - '''
> - res = imConv(im,ker)
> - #res = imConv(im,ker,seperableConv2)
> - res = res[(skip-1)//2::skip, (skip-1)//2::skip, :]
> - res = clamp(res)
> - return res
> -
> -
> -#----------------------------------------------------------------------
> -def trilinear():
> - '''
> - Trilinear interpolation (not implemented...)
> - '''
> - pass
> -
> -
> -#----------------------------------------------------------------------
> -def clamp(im):
> - '''
> - Clamp an image between zero and one.
> - '''
> - return numpy.minimum(numpy.maximum(im,0),1)
> -
> -
> -#----------------------------------------------------------------------
> -def getGridImg(N, gridSize):
> - '''
> - Get a test image containig a one pixel wide grid
> -
> - N - side length
> - gridSize - size of the one-pixel wide black grid in the image.
> - '''
> -
> - im = zeros((N, N, 3))
> -
> - x = linspace(0,1,N)
> - [xx,yy] = meshgrid(x,x)
> - im[:,:,0] = xx
> - im[:,:,1] = 0.5
> - im[:,:,2] = yy
> -
> - im[0:-gridSize-1:gridSize,:,:] = 0
> - im[:,0:-gridSize-1:gridSize,:] = 0
> -
> - im[:,-1,:] = 0
> - im[-1,:,:] = 0
> -
> - return im
> -
> -
> -#----------------------------------------------------------------------
> -#----------------------------------------------------------------------
> -# Functions for creating/manipulating mipmaps
> -#----------------------------------------------------------------------
> -def genMipmapIter(im, width, selectWidth, windowFunc, kerFunc):
> - '''
> - Generate a mipmap by downsampling each image in turn to get the
> next level.
> -
> - im - input image
> - width - desired "width" of filter kernel. the kernel will have
> width+1 points.
> - selectWidth - whether the function should try to adjust the filter
> width for
> - even/odd sized images. To avoid loosing information, even
> - filters should go with even sized images and vice versa
> - windowFunc - windowing function to use on the filter.
> -
> - '''
> -
> - mipmap = []
> -
> - scale = 2
> -
> - ker = kerFunc(width, scale)
> - ker = applyWindow(ker, windowFunc)
> -
> - while size(im,0) > 1:
> - mipmap.append(im)
> -
> - if selectWidth:
> - if (size(im,0) + width + 1) % 2 == 0:
> - ker = kerFunc(width,scale)
> - else:
> - ker = kerFunc(width+1,scale)
> - ker = applyWindow(ker, windowFunc)
> -
> - print 'side length = %d, number of points in filter = %d'
> \
> - % (size(im,0), len(ker))
> - im = imDownsamp(im, ker, scale)
> -
> - return mipmap
> -
> -
> -
> -#----------------------------------------------------------------------
> -def genMipmapIterDownsamp4(im, width, selectWidth, windowFunc, kerFunc):
> - '''
> - Same as genMipmapIter, except tries to downsample by a factor of 4
> rather
> - than two when possible.
> -
> - (Beware: ugly hack)
> - '''
> - pass
> -
> -
> -#----------------------------------------------------------------------
> -def flattenImageList(mipmap, dim):
> - '''
> - "Flatten" mipmap levels into a single image.
> -
> - mipmap - cell array of scaled images
> - dim - dimension number to lay the images out along
> - '''
> -
> - len = 0
> - for mipLevel in mipmap:
> - len += size(mipLevel,dim)
> - mipSize = list(mipmap[0].shape)
> - mipSize[dim] = len
> -
> - mipmapImg = ones(mipSize)*0.7
> - pos = 0
> - for mipLevel in mipmap:
> - if dim == 0:
> - mipmapImg[pos:pos+size(mipLevel,0),
> 0:size(mipLevel,1), :] = mipLevel
> - else:
> - mipmapImg[0:size(mipLevel,0),
> pos:pos+size(mipLevel,1), :] = mipLevel
> - pos += size(mipLevel,0)
> -
> - return mipmapImg
> -
> -
> -#----------------------------------------------------------------------
> -# High level driver functions.
> -#----------------------------------------------------------------------
> -def compareOddEven(im):
> - '''
> - Compare the difference between sampling an image with an odd or
> even-sized
> - filter.
> - '''
> -
> - kerFunc = sincKer
> - windowFunc = makeTrimmedWindow(lanczos)
> -
> - # Use an odd-width filter
> - mipmap = genMipmapIter(im, 7, False, windowFunc, kerFunc)
> - imsave('mipmap_odd_width.png', flattenImageList(mipmap,1))
> -
> - # Use an even-width filter
> - mipmap = genMipmapIter(im, 8, False, windowFunc, kerFunc)
> - imsave('mipmap_even_width.png', flattenImageList(mipmap,1))
> -
> -
> -#----------------------------------------------------------------------
> -def investigateWindowTypes(im, kerFunc):
> - '''
> - grab the level-three mipmaps arising from using a bunch of windows
> with
> - a fixed filter kernel.
> - '''
> -
> - # generate a version using the boxcar kernel [1,1] for reference.
> - mipmap = genMipmapIter(im, 1, True, boxcar, boxker)
> - miplevel3 = [mipmap[3]]
> -
> - winfuncs = [lanczos, hamming, hanning, blackman, boxcar]
> - # generate mipmaps using the other window functions
> - for func in winfuncs:
> - mipmap = genMipmapIter(im, 7, True, func, kerfunc)
> - miplevel3.append(mipmap[3])
> -
> - imsave('mipmap_lvl3_all_wins.png', flattenImageList(miplevel3,1))
> -
> -
> -#----------------------------------------------------------------------
> -def plotWindowFuncs(N):
> - '''
> - Plot various window functions.
> - '''
> - # Try kaiser window - just need to choose a width...
> - winFuncs = [lanczos, hamming, hanning, blackman, boxcar]
> - winFuncNames = ['lanczos', 'hamming', 'hanning', 'blackman',
> 'boxcar']
> - clf()
> - hold(True)
> - for func in winFuncs:
> - plot(func(N))
> -
> - legend(winFuncNames)
> - axis([1,N,0,1.05])
> -
> -
> -#----------------------------------------------------------------------
> -# Main program
> -#----------------------------------------------------------------------
> -
> -if __name__ == '__main__':
> - # Choose image
> - im = imread('lena_std.png')[:,:,:3]/256.0
> - #im = getGridImg(246, 20)
> -
> - # Compare odd and even filter sizes on the same image
> - compareOddEven(im)
> -
> - # Compare convolution strategies for square seperable filters.
> -
> - # Get a mipmap
> - #mipmap = genMipmapIter(im, 7, True, boxcar, sincKer)
> - # show/save mipmap
> - #mipAll = flattenImageList(mipmap, 1)
> - #pylab.imshow(mipAll, interpolation='nearest')
> - #pylab.imshow(mipAll)
> - #imsave('mipmap.png', mipAll)
>
> Copied:
> branches/soc2007-dsm/aqsis/testing/prototypes/texfilt/downsample_test.py
> (from rev 1048, trunk/aqsis/testing/prototypes/texfilt/downsample_test.py)
> ===================================================================
> ---
> branches/soc2007-dsm/aqsis/testing/prototypes/texfilt/downsample_test.py (rev
> 0)
> +++
> branches/soc2007-dsm/aqsis/testing/prototypes/texfilt/downsample_test.py 2007-05-19
> 02:48:41 UTC (rev 1049)
> @@ -0,0 +1,337 @@
> +#!/usr/bin/python
> +#
> +# Scientific python script to test image filtering for mipmap creation in
> Aqsis
> +# Copyright (C) 2007, Christopher J. Foster
> +#
> +# This program is free software; you can redistribute it and/or
> +# modify it under the terms of the GNU General Public License
> +# as published by the Free Software Foundation; either version 2
> +# of the License, or (at your option) any later version.
> +#
> +# This program is distributed in the hope that it will be useful,
> +# but WITHOUT ANY WARRANTY; without even the implied warranty of
> +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
> +# GNU General Public License for more details.
> +#
> +# You should have received a copy of the GNU General Public License
> +# along with this program; if not, write to the Free Software
> +# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
> MA 02110-1301, USA.
> +
> +# This script investigates some of the factors important when
> downsampling an
> +# image:
> +#
> +# * Filter types
> +# * Filter windows
> +# * Correctly filtering odd and even-sized images
> +
> +from __future__ import division
> +
> +import matplotlib.pylab as pylab
> +
> +import numpy
> +from numpy import r_, size, zeros, ones, ceil, floor, array, linspace,
> meshgrid
> +
> +import scipy
> +from scipy.signal import boxcar, convolve2d as conv2
> +from scipy.misc import imsave, imread
> +
> +#----------------------------------------------------------------------
> +#----------------------------------------------------------------------
> +# Windowing functions
> +#----------------------------------------------------------------------
> +def lanczos(N, tau=1):
> + '''
> + return N-point Lanczos window with width tau
> + '''
> + return scipy.sinc(linspace(-1,1,N)/tau)
> +
> +
> +#----------------------------------------------------------------------
> +def makeTrimmedWindow(winFunc):
> + '''
> + Return a window function constructed from the given window
> function, but
> + without the two edge points (see natWin).
> + '''
> + return lambda N, *args: winFunc(N+2, *args)[1:-1]
> +
> +
> +#----------------------------------------------------------------------
> +def applyWindow(ker, windowFunc):
> + '''
> + Window a filter kernel with the given windowing function.
> + '''
> + window = windowFunc(len(ker))
> + window = numpy.outer(window, window)
> + return ker * window
> +
> +
> +#----------------------------------------------------------------------
> +#----------------------------------------------------------------------
> +# Filter Kernels
> +#----------------------------------------------------------------------
> +def sincKer(width, scale):
> + '''
> + Return a sinc filter kernel with specified width and scale
> +
> + width - total number of points in the filter (not radius!)
> + scale - zeros will be centered around 0 with spacing "scale"
> + '''
> + x = r_[0:width+1]-width*0.5
> + k = numpy.outer(scipy.sinc(x/scale), scipy.sinc(x/scale))
> + return k
> +
> +
> +#----------------------------------------------------------------------
> +def boxKer(width, scale):
> + '''
> + Trivial box filter kernel
> + '''
> + return ones(width+1)
> +
> +
> +#----------------------------------------------------------------------
> +# Image manipuation
> +#----------------------------------------------------------------------
> +def seperableConv2(A, k, *args):
> + '''
> + perform a 2D convolution, assuming that the convolution kernel k
> is
> + square-seperable. In this case we may just do two 1D
> convolutions.
> + '''
> + center = array(k.shape)//2
> + k0 = k[center[0]:center[0]+1,:]
> + k1 = k[:,center[1]:center[1]+1]
> + return conv2(conv2(A,k0, *args),k1, *args)
> +
> +
> +#----------------------------------------------------------------------
> +def imConv(im, ker, conv2Func=conv2):
> + '''
> + Convolve an image with the given filter kernel
> + '''
> + res = zeros(im.shape)
> + normalisation = conv2Func(ones(im.shape[0:2],'d'), ker, 'full')
> + # compute how much of the full convolution to trim from top left
> and bottom right
> + trimTl = numpy.int_(ceil((array(ker.shape)-1)/2))
> + trimBr = numpy.int_(floor((array(ker.shape)-1)/2))
> + for ii in range(0,size(im,2)):
> + fullconv = conv2Func(im[:,:,ii], ker, 'full') /
> normalisation
> + res[:,:,ii] = fullconv[trimTl[0]:-trimBr[0],
> trimTl[1]:-trimBr[1]]
> + return res
> +
> +
> +#----------------------------------------------------------------------
> +def imDownsamp(im, ker, skip):
> + '''
> + Downsample an image by
> + (1) Applying the filter kernel ker at each pixel
> + (2) Decimate by taking elements 0::skip from the result
> + (3) Clamp the resulting intensities between 0 and 1
> + '''
> + res = imConv(im,ker)
> + #res = imConv(im,ker,seperableConv2)
> + res = res[(skip-1)//2::skip, (skip-1)//2::skip, :]
> + res = clamp(res)
> + return res
> +
> +
> +#----------------------------------------------------------------------
> +def trilinear():
> + '''
> + Trilinear interpolation (not implemented...)
> + '''
> + pass
> +
> +
> +#----------------------------------------------------------------------
> +def clamp(im):
> + '''
> + Clamp an image between zero and one.
> + '''
> + return numpy.minimum(numpy.maximum(im,0),1)
> +
> +
> +#----------------------------------------------------------------------
> +def getGridImg(N, gridSize):
> + '''
> + Get a test image containig a one pixel wide grid
> +
> + N - side length
> + gridSize - size of the one-pixel wide black grid in the image.
> + '''
> +
> + im = zeros((N, N, 3))
> +
> + x = linspace(0,1,N)
> + [xx,yy] = meshgrid(x,x)
> + im[:,:,0] = xx
> + im[:,:,1] = 0.5
> + im[:,:,2] = yy
> +
> + im[0:-gridSize-1:gridSize,:,:] = 0
> + im[:,0:-gridSize-1:gridSize,:] = 0
> +
> + im[:,-1,:] = 0
> + im[-1,:,:] = 0
> +
> + return im
> +
> +
> +#----------------------------------------------------------------------
> +#----------------------------------------------------------------------
> +# Functions for creating/manipulating mipmaps
> +#----------------------------------------------------------------------
> +def genMipmapIter(im, width, selectWidth, windowFunc, kerFunc):
> + '''
> + Generate a mipmap by downsampling each image in turn to get the
> next level.
> +
> + im - input image
> + width - desired "width" of filter kernel. the kernel will have
> width+1 points.
> + selectWidth - whether the function should try to adjust the filter
> width for
> + even/odd sized images. To avoid loosing information, even
> + filters should go with even sized images and vice versa
> + windowFunc - windowing function to use on the filter.
> +
> + '''
> +
> + mipmap = []
> +
> + scale = 2
> +
> + ker = kerFunc(width, scale)
> + ker = applyWindow(ker, windowFunc)
> +
> + while size(im,0) > 1:
> + mipmap.append(im)
> +
> + if selectWidth:
> + if (size(im,0) + width + 1) % 2 == 0:
> + ker = kerFunc(width,scale)
> + else:
> + ker = kerFunc(width+1,scale)
> + ker = applyWindow(ker, windowFunc)
> +
> + print 'side length = %d, number of points in filter = %d'
> \
> + % (size(im,0), len(ker))
> + im = imDownsamp(im, ker, scale)
> +
> + return mipmap
> +
> +
> +
> +#----------------------------------------------------------------------
> +def genMipmapIterDownsamp4(im, width, selectWidth, windowFunc, kerFunc):
> + '''
> + Same as genMipmapIter, except tries to downsample by a factor of 4
> rather
> + than two when possible.
> +
> + (Beware: ugly hack)
> + '''
> + pass
> +
> +
> +#----------------------------------------------------------------------
> +def flattenImageList(mipmap, dim):
> + '''
> + "Flatten" mipmap levels into a single image.
> +
> + mipmap - cell array of scaled images
> + dim - dimension number to lay the images out along
> + '''
> +
> + len = 0
> + for mipLevel in mipmap:
> + len += size(mipLevel,dim)
> + mipSize = list(mipmap[0].shape)
> + mipSize[dim] = len
> +
> + mipmapImg = ones(mipSize)*0.7
> + pos = 0
> + for mipLevel in mipmap:
> + if dim == 0:
> + mipmapImg[pos:pos+size(mipLevel,0),
> 0:size(mipLevel,1), :] = mipLevel
> + else:
> + mipmapImg[0:size(mipLevel,0),
> pos:pos+size(mipLevel,1), :] = mipLevel
> + pos += size(mipLevel,0)
> +
> + return mipmapImg
> +
> +
> +#----------------------------------------------------------------------
> +# High level driver functions.
> +#----------------------------------------------------------------------
> +def compareOddEven(im):
> + '''
> + Compare the difference between sampling an image with an odd or
> even-sized
> + filter.
> + '''
> +
> + kerFunc = sincKer
> + windowFunc = makeTrimmedWindow(lanczos)
> +
> + # Use an odd-width filter
> + mipmap = genMipmapIter(im, 7, False, windowFunc, kerFunc)
> + imsave('mipmap_odd_width.png', flattenImageList(mipmap,1))
> +
> + # Use an even-width filter
> + mipmap = genMipmapIter(im, 8, False, windowFunc, kerFunc)
> + imsave('mipmap_even_width.png', flattenImageList(mipmap,1))
> +
> +
> +#----------------------------------------------------------------------
> +def investigateWindowTypes(im, kerFunc):
> + '''
> + grab the level-three mipmaps arising from using a bunch of windows
> with
> + a fixed filter kernel.
> + '''
> +
> + # generate a version using the boxcar kernel [1,1] for reference.
> + mipmap = genMipmapIter(im, 1, True, boxcar, boxker)
> + miplevel3 = [mipmap[3]]
> +
> + winfuncs = [lanczos, hamming, hanning, blackman, boxcar]
> + # generate mipmaps using the other window functions
> + for func in winfuncs:
> + mipmap = genMipmapIter(im, 7, True, func, kerfunc)
> + miplevel3.append(mipmap[3])
> +
> + imsave('mipmap_lvl3_all_wins.png', flattenImageList(miplevel3,1))
> +
> +
> +#----------------------------------------------------------------------
> +def plotWindowFuncs(N):
> + '''
> + Plot various window functions.
> + '''
> + # Try kaiser window - just need to choose a width...
> + winFuncs = [lanczos, hamming, hanning, blackman, boxcar]
> + winFuncNames = ['lanczos', 'hamming', 'hanning', 'blackman',
> 'boxcar']
> + clf()
> + hold(True)
> + for func in winFuncs:
> + plot(func(N))
> +
> + legend(winFuncNames)
> + axis([1,N,0,1.05])
> +
> +
> +#----------------------------------------------------------------------
> +# Main program
> +#----------------------------------------------------------------------
> +
> +if __name__ == '__main__':
> + # Choose image
> + im = imread('lena_std.png')[:,:,:3]/256.0
> + #im = getGridImg(246, 20)
> +
> + # Compare odd and even filter sizes on the same image
> + compareOddEven(im)
> +
> + # Compare convolution strategies for square seperable filters.
> +
> + # Get a mipmap
> + #mipmap = genMipmapIter(im, 7, True, boxcar, sincKer)
> + # show/save mipmap
> + #mipAll = flattenImageList(mipmap, 1)
> + #pylab.imshow(mipAll, interpolation='nearest')
> + #pylab.imshow(mipAll)
> + #imsave('mipmap.png', mipAll)
>
> Deleted: branches/soc2007-dsm/aqsis/testing/prototypes/texfilt/texfilt.py
> ===================================================================
> --- trunk/aqsis/testing/prototypes/texfilt/texfilt.py 2007-05-15
> 12:34:19 UTC (rev 1048)
> +++
> branches/soc2007-dsm/aqsis/testing/prototypes/texfilt/texfilt.py 2007-05-19
> 02:48:41 UTC (rev 1049)
> @@ -1,893 +0,0 @@
> -#!/usr/bin/python
> -#
> -# Scientific python script to test image filtering for mipmap creation in
> Aqsis
> -# Copyright (C) 2007, Christopher J. Foster
> -#
> -# This program is free software; you can redistribute it and/or
> -# modify it under the terms of the GNU General Public License
> -# as published by the Free Software Foundation; either version 2
> -# of the License, or (at your option) any later version.
> -#
> -# This program is distributed in the hope that it will be useful,
> -# but WITHOUT ANY WARRANTY; without even the implied warranty of
> -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
> -# GNU General Public License for more details.
> -#
> -# You should have received a copy of the GNU General Public License
> -# along with this program; if not, write to the Free Software
> -# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
> MA 02110-1301, USA.
> -
> -'''
> -This script investigates texture access. There are two main parts:
> - 1) Mipmap generation & access
> - 2) Texture filtering using anisotropic filtering. This includes
> - * stochastic filtering with a weight function over a
> quadrilateral
> - * elliptical gaussian filtering (EWA)
> -
> -To run the script, you'll need scientific python; specifically, the
> packages
> - * scipy
> - * matplotlib
> -On a unix machine, ipython is highly recommended for a nice interactive
> -interface which meshes very well with matplotlib.
> -'''
> -
> -from __future__ import division
> -
> -import matplotlib.pylab as pylab
> -from pylab import imshow
> -
> -import numpy
> -from numpy import pi, r_, size, zeros, ones, eye, ceil, floor, array, \
> - linspace, meshgrid, random, arange, dot, diag, asarray
> -
> -import scipy
> -from scipy import cos, sin, exp, log, log2, sqrt, ogrid, mgrid, ndimage
> -import scipy.linalg as linalg
> -from scipy.misc import imsave, imread
> -
> -# The purpose of this script is to investigate mipmap generation and
> indexing
> -# in more-or-less the same way which it will be used internally by Aqsis.
> -
> -#----------------------------------------------------------------------
> -#----------------------------------------------------------------------
> -# Filter Kernels, including any windowing.
> -#----------------------------------------------------------------------
> -def sincKer(width, even, scale = 2):
> - '''
> - Return a sinc filter kernel with specified width and scale
> -
> - width - width of the filter (used by the window)
> - even - true if the generated filter should have an even number of
> points.
> - scale - zeros will be centred around 0 with spacing "scale"
> - '''
> - if even:
> - numPts = max(2*int((width+1)/2), 2)
> - else:
> - numPts = max(2*int(width/2) + 1, 3)
> - x = r_[0:numPts]-(numPts-1)*0.5
> - k = scipy.sinc(x/scale)
> - k *= cos(linspace(-pi/4,pi/4,len(k)))
> - ker = numpy.outer(k,k)
> - ker /= sum(ker.ravel())
> - return ker
> -
> -
> -#----------------------------------------------------------------------
> -# Image manipuation
> -#----------------------------------------------------------------------
> -def applyToChannels(img, fxn, *args, **kwargs):
> - '''
> - Apply the given function to each of the channels in img in turn.
> - Additional positional and keyword arguments are passed to the
> called fxn.
> - '''
> - numChannels = img.shape[-1]
> - outList = []
> - for i in range(0,numChannels):
> - outList.append(fxn(img[:,:,i], *args, **kwargs))
> - outArr = array(outList)
> - return outArr.transpose(range(1,len(outArr.shape)) + [0])
> -
> -
> -#----------------------------------------------------------------------
> -#----------------------------------------------------------------------
> -# Functions for creating test images
> -#----------------------------------------------------------------------
> -def getGridImg(N, gridSize, gridWidth=1):
> - '''
> - Get a test image containig a one pixel wide grid
> -
> - N - side length
> - gridSize - size of the one-pixel wide black grid in the image.
> - '''
> -
> - im = zeros((N, N, 3))
> -
> - x = linspace(0,1,N)
> - [xx,yy] = meshgrid(x,x)
> - im[:,:,0] = xx
> - im[:,:,1] = 0.5
> - im[:,:,2] = yy
> -
> - for i in range(0,gridWidth):
> - im[i:-gridSize//2:gridSize,:,:] = 0
> - im[:,i:-gridSize//2:gridSize,:] = 0
> -
> - im[:,-1,:] = 0
> - im[-1,:,:] = 0
> -
> - return im
> -
> -
> -def getCentreTest(N=512):
> - '''
> - Get a test image consisting of a 2x2 checkered plane.
> -
> - N - side length
> - '''
> - im = zeros((N, N, 3))
> -
> - im[:N/2,N/2:,:] = 1
> - im[N/2:,:N/2,:] = 1
> -
> - return im
> -
> -
> -def getEdgeTest(N=512, borderwidth=20):
> - '''
> - Get a white test image with a black bordering strip
> -
> - N - side length
> - '''
> - im = ones((N, N, 3))
> -
> - im[borderwidth:2*borderwidth,:,:] = 0
> - im[-2*borderwidth:-borderwidth,:,:] = 0
> - im[:,borderwidth:2*borderwidth,:] = 0
> - im[:,-2*borderwidth:-borderwidth,:] = 0
> -
> - return im
> -
> -#----------------------------------------------------------------------
> -#----------------------------------------------------------------------
> -# Utils for creating/manipulating mipmaps
> -#----------------------------------------------------------------------
> -
> -class ImageDownsampler:
> - '''
> - Class which holds a filter kernel (function & width) and is able
> to
> - correctly use it to downsample odd or even-sized images.
> - '''
> - def __init__(self, kerWidth, kerFunc=sincKer):
> - self.kerFunc = kerFunc
> - self.kerWidth = kerWidth
> - def downsample(self, img):
> - '''
> - Downsample an image by
> - (1) Applying the filter kernel at each pixel
> - (2) Decimate by taking every second element from the
> result
> - (3) Clamp the resulting intensities between 0 and 1
> - '''
> - ker = self.kerFunc(self.kerWidth, (size(img,0) % 2) == 0)
> - res = applyToChannels(img, ndimage.correlate, ker,
> mode='constant', cval=0)
> - res = res[(1-(size(img,0) % 2))::2, (1-(size(img,1) %
> 2))::2, :]
> - res = res.clip(min=0, max=1)
> - return res
> -
> -#----------------------------------------------------------------------
> -class MipLevel:
> - '''
> - Class to hold a level of mipmap, which knows how to bilinear
> interpolation
> - on itself.
> - '''
> - def __init__(self, image, offset, mult):
> - '''
> - image - pixel data; a NxMx3 array
> - offset, mult - length two vectors which specify an affine
> - transformation from texture coordinates to this
> mipmap's raster
> - coordinates:
> - i = mult[0] * s + offset[0]
> - j = mult[1] * t + offset[1]
> - '''
> - # store the image with some ghost points to work around a
> bug in
> - # scipy.ndimage, where incorrect interpolation happens for
> negative
> - # indices.
> - self.image = zeros((image.shape[0]+2,image.shape[1]+2,
> image.shape[2]))
> - self.image[1:-1,1:-1,:] = image
> - self.offset = array(offset,'d')
> - self.mult = array(mult,'d')
> - # Variances for the reconstruction & prefilters when using
> EWA. A
> - # variance of 1/(2*pi) gives a filter with centeral weight
> 1, but in
> - # practise this is slightly too small (resulting in a
> little bit of
> - # aliasing). Therefore it's adjusted up slightly.
> - self.varianceXY = 1.3/(2*pi)
> - self.varianceST = 1.3/(2*pi)
> -
> - def rasterCoords(self,s,t):
> - '''
> - Return the raster coordinates for this mipmap from the
> texture
> - coordinates (s,t)
> - '''
> - i = scipy.atleast_1d(self.mult[0]*s + self.offset[0] + 1)
> - j = scipy.atleast_1d(self.mult[1]*t + self.offset[1] + 1)
> - return (i,j)
> -
> - def sampleBilinear(self, s, t):
> - '''
> - Sample the mipmap level at the coordinates specified by s
> and t. For
> - performance, s & t can be 1D arrays.
> - '''
> - i,j = self.rasterCoords(s,t)
> - return applyToChannels(self.image, ndimage.map_coordinates,
> \
> - array([i,j]), order=1, mode='constant',
> cval=0)
> -
> - def getRandPointsInQuad(self, s, t, numPoints):
> - '''
> - Get a random set of points inside the quadrilatiral given
> by the
> - texture-space vertices
> - [s[1], t[1]], ... , [s[4], t[4]].
> -
> - The sample should be uniform if the quadrialteral can be
> made from a
> - square with a _linear_ transformation. Otherwise it will
> be biased.
> - The form of biasing is probably exactly what we want if
> we're taking
> - pixel samples...
> - '''
> - lerp1 = numpy.random.rand(numPoints)
> - lerp1m = 1-lerp1
> - lerp2 = numpy.random.rand(numPoints)
> - st = array([s,t])
> - op = numpy.outer
> - randPts = lerp2*(op(st[:,0],lerp1) + op(st[:,1],lerp1m)) \
> - + (1-lerp2)*(op(st[:,3],lerp1) +
> op(st[:,2],lerp1m))
> - #weights = exp(-8*((lerp1-0.5)**2 + (lerp2-0.5)**2))
> - weights = ones(numPoints)
> - return (randPts[0,:], randPts[1,:], weights)
> -
> - def filterQuadAF(self, s, t, numSamples):
> - '''
> - Filter the texture by sampling inside the given
> quadriateral specified
> - by the texture-space vertices
> - [s[1], t[1]], ... , [s[4], t[4]],
> - and taking an average. Sample points are chosen
> stochastically over
> - the quadrilateral. This is essentially Monte Carlo
> integration, with a
> - reconstruction filter equal to bilinear interpolation.
> - '''
> - sRnd, tRnd, w = self.getRandPointsInQuad(s, t, numSamples)
> - samples = self.sampleBilinear(sRnd, tRnd)
> - return ( samples*w.reshape(w.shape+(1,)) ).sum(0) / w.sum
> ()
> -
> - def estimateJacobianInverse(self, s, t):
> - '''
> - Estimate the inverse Jacobian from a small (s,t) box using
> finite
> - differences. The assumption is that the (s,t) box has
> come from a 1x1
> - box in the final coordinate system.
> - '''
> - return 0.5 * array([[s[1]-s[0] + s[2]-s[3], s[3]-s[0] +
> s[2]-s[1]],
> - [t[1]-t[0] +
> t[2]-t[3], t[3]-t[0] + t[2]-t[1]]])
> -
> - def filterEWA(self, sBox, tBox, deltaX=None, deltaY=None, J=None,
> returnWeights=False):
> - '''
> - Filter a portion of the texture using an elliptical
> gaussian filter, as
> - described in Heckbert's 1989 thesis, "Fundamentals of
> Texture Mapping
> - and Image Warping". This kind of filter is also called an
> EWA filter for
> - "Elliptical Weighted Average".
> -
> - s, t - texture coordinates of a box to be filtered over
> - deltaX, deltaY - pixel side length in the final discrete
> image.
> - J - Jacobian of the coordinate transformation.
> - returnWeights - true if you just want the filter weights
> returned for
> - debugging rather than the
> filtered image.
> - '''
> - # convert to raster coords
> - sBox,tBox = self.rasterCoords(sBox,tBox)
> - # centre of box
> - s0,t0 = sBox.mean(), tBox.mean()
> - if not returnWeights and (s0 < 0 or s0 > self.image.shape[0]
> \
> - or t0 < 0
> or t0 > self.image.shape[1]):
> - return 0
> - if J is None:
> - # Estimate the Jacobian from the texture box.
> - Ji = self.estimateJacobianInverse(sBox, tBox)
> - else:
> - # Else we have an analytical Jacobian specified.
> - # Adjust J for raster->texture coords
> - J = dot(dot(diag((1/deltaX, 1/deltaY)), J),
> diag((1/self.image.shape[0], 1/self.image.shape[1])))
> - # get the inverse of J
> - Ji = linalg.inv(J)
> - # width parameter is log(n) where the filter falls to 1/n
> of it's
> - # centeral value at the edge.
> - width = 4
> - # V = covariance matrix
> - V = self.varianceXY*dot(Ji, Ji.transpose()) +
> self.varianceST*eye(2)
> - # Q = (+ve definite) quadratic form matrix.
> - Q = 0.5*linalg.inv(V)
> - det = Q[0,0]*Q[1,1] - Q[0,1]*Q[1,0]
> - # compute bounding radii
> - sRad = sqrt(Q[1,1]*width/det)
> - tRad = sqrt(Q[0,0]*width/det)
> - # integer filter width
> - sWidth = int(floor(s0+sRad)) - int(ceil(s0-sRad)) + 1
> - tWidth = int(floor(t0+tRad)) - int(ceil(t0-tRad)) + 1
> - # compute edge of filter region in level0 image:
> - sStart = ceil(s0-sRad)
> - tStart = ceil(t0-tRad)
> - # compute filter weights. I'm being sloppy here, and not
> bothering to
> - # truncate outside the relevant ellipse
> - s,t = ogrid[0:sWidth, 0:tWidth]
> - s = s - (s0-sStart)
> - t = t - (t0-tStart)
> - if not returnWeights and (sStart < 0 or sStart + sWidth-1
> >= self.image.shape[0] \
> - or tStart < 0 or tStart + tWidth-1 >=
> self.image.shape[1]):
> - return 0.5
> - weights = exp(-(Q[0,0]*s*s + (Q[0,1]+Q[1,0])*s*t +
> Q[1,1]*t*t))
> - weights /= weights.sum()
> - if returnWeights:
> - return weights
> - imSec = self.image
> [sStart:sStart+sWidth,tStart:tStart+tWidth]
> - return applyToChannels(imSec, lambda x: (x*weights).sum())
> -
> - def display(self, s=None, t=None):
> - if s is None or t is None:
> - imshow(self.image[1:-1,1:-1,:],
> interpolation='nearest')
> - else:
> - imshow(self.sampleBilinear(s,t),
> interpolation='nearest')
> -
> -
> -#----------------------------------------------------------------------
> -class TextureMap:
> - '''
> - Class to encapsulate a (power-of-2) mipmapped texture in much the
> same way
> - mipmapping will be done internally in Aqsis.
> - '''
> - def __init__(self, img, filterWidth, kerFunc=sincKer, \
> - levelCalcMethod='minQuadWidth', numSamples=100):
> - '''
> - img - texture image
> - filterWidth - filter width for mipmap downsampling
> - kerFunc - kernel function for mipmap downsampling
> - levelCalcMethod - method used for calculating required
> mipmap level
> - from a quadrialteral. filtering. May be one of:
> 'minSideLen',
> - 'minQuadWidth', 'trilinear', 'minDiag', 'level0'.
> - numSamples - number of samples to used when performing
> stochastic
> - filtering
> - '''
> - self.downsampler = ImageDownsampler(filterWidth, kerFunc)
> - self.levels = []
> - self.levelCalcMethod = levelCalcMethod
> - self.genMipMap(img)
> - self.numSamples = numSamples
> -
> - def genMipMap(self, img):
> - '''
> - Generate a mipmap, by downsampling by successive powers of
> two, using
> - the held kernel function and width.
> -
> - In this function we make the requirement that the filter
> coefficients
> - for filtering between levels are constant across the image
> for a
> - particular mipmap level. This restriction means that
> even-sized images
> - have to be filtered with even- sized filter kernels, and
> vice versa.
> - Even sizes also require that an offset be taken into
> account, since the
> - edges of the downsampled levels don't lie on the edges of
> the original
> - image.
> -
> - Examples showing mipmap sample points (represented as x's)
> follow for
> - the 1D case. Several variables are also shown as a
> function of the
> - mipmap level, l on the rhs:
> -
> - i_l(t) = raster coordinates for level l as a function of
> texture coordinate t.
> - o_l = offset for level l in level 0 units
> - s_l = span of level l in level 0 units
> - d_l = interval between pixels of level l in level 0 raster
> units.
> - w_l = width of level l in pixels.
> -
> -
> - (mostly) odd numbers of points:
> - 0 1 2 3 4 <-- i_0
> - l | | | |
> | o_l s_l w_l d_l
> - 0 [x x x x x] 0 4 5 1
> - 1 [x x x] 0 2 3 2
> - 2 [x x] 0 1 2 4
> - 3 [ x ] 2 0 1 ?
> -
> - Another example (powers of two)
> - 0 1 2 3 4 5 6 7 <-- i_0
> - l | | | | | | | | o_l s_l w_l
> - 0 [x x x x x x x x] 0 7 8
> - 1 [ x x x x ] 0.5 6 4
> - 2 [ x x ] 1.5 4 2
> - 3 [ x ] 3.5 0 1
> -
> - Notice that the higher level mipmaps for the power of
> two-sized image
> - do *not* have points which lie on the edge of the
> image. This means
> - that we've got to incorporate an offset when we sample the
> level.
> -
> - Formulas for texture coordinates -> raster space
> -
> - s_l = (s_0 - 2*o_l)
> -
> - Raster units for level 0 are:
> - i_0(t) = t * s_0
> - Raster units for level l are:
> - i_l(t) = (i_0(t) - o_0)/d_l
> -
> - Or, what we really need, which is raster units for level l
> in terms of
> - the texture coordinate t:
> - i_l(t) = t * s_0/d_l - o_0/d_l
> -
> - ((
> - Note that we could abandon the
> constant-coefficients approach
> - (which would make it very slow to generate
> mipmaps), we would have:
> -
> - 0 1 2 3 4 5 6 7 <-- i_0
> - l | | | | | | | | o_l s_l w_l
> - 0 [x x x x x x x x] 0 7 8
> - 1 [x x x x] 0 7 4
> - 2 [x x] 0 7 2
> - 3 [ x ] 0 7 1
> -
> - The the raster units for level l in terms of the
> texture coordinates t
> - are then very simple:
> -
> - i_l(t) = t*(w_l - 1)
> - ))
> - '''
> - scale = 2
> - level0Span = array(img.shape[0:2])-1
> - offset = array([0,0],'d')
> - self.levels = [MipLevel(img, offset, level0Span)]
> - # delta is distance between pixels; level 0 units
> - delta = 1
> - # offsets are between 1st pixel of 0th level to 1st pixel
> of current level
> - #
> - # eg, in 1D
> - # [x x x x..]
> - # <---> <- level 0 delta
> - # [ x x ...]
> - # <-> <- level 1 offset
> - while (array(img.shape[0:2]) > 1).all():
> - # calculate start offsets for this level.
> - offset += (1-(array(img.shape[0:2]) % 2))*delta/2
> - delta *= 2
> - # downsample
> - img = self.downsampler.downsample(img)
> - # create a new MipLevel to hold the image.
> - self.levels.append(MipLevel(img, -offset/delta,
> level0Span/delta))
> - # The below is a theoretically broken, but very
> usable version, -
> - # mipmap offsets are set to zero.
> - #self.levels.append(MipLevel(img, 0*offset, array(
> img.shape[0:2])-1))
> - print 'generated level size', img.shape[0:2]
> -
> - def filterTrilinear(self, sBox, tBox, level=None):
> - '''
> - Trilinear filtering over the mipmaps
> -
> - s, t - texture coordinates of points to sample at
> - ds,dt - size of a "box" on the screen.
> - level - mipmap level to use. If equal to None, calculate
> from ds & dt
> - '''
> - if level is None:
> - level = min(len(self.levels)-1,
> self.calculateLevel(sBox, tBox))
> - level = max(level,0)
> - level1 = int(floor(level))
> - interp = level - level1
> - level2 = min(len(self.levels)-1, int(level1 + 1))
> - s0 = sBox.mean()
> - t0 = tBox.mean()
> - samp1 = self.levels[level1].sampleBilinear(s0,t0)
> - samp2 = self.levels[level2].sampleBilinear(s0,t0)
> - return (1-interp)*samp1 + interp*samp2
> -
> - def calculateLevel(self, s, t):
> - '''
> - Calculate the appropriate mipmap level for texture
> filtering over a
> - quadrilateral given by the texture-space vertices
> - [s[1], t[1]], ... , [s[4], t[4]].
> -
> - There are many ways to do this; the instance variable
> levelCalcMethod
> - selects the desired one. The most correct way is to
> choose the
> - minSideLen method, as long as the quadrilateral is vaguely
> - rectangular-shaped. This only works if you're happy to
> use lots of
> - samples however, otherwise you get aliasing.
> - '''
> - s = s.copy()*self.levels[0].image.shape[0]
> - t = t.copy()*self.levels[0].image.shape[1]
> - if self.levelCalcMethod == 'minSideLen':
> - # Get mipmap level with minimum feature size equal
> to the shortest
> - # quadrilateral side
> - s1 = pylab.concatenate((s, s[0:1]))
> - t1 = pylab.concatenate((t, t[0:1]))
> - minSideLen2 = (numpy.diff(s1)**2 + numpy.diff
> (t1)**2).min()
> - level = log2(minSideLen2)/2
> - elif self.levelCalcMethod == 'minQuadWidth':
> - # Get mipmap level with minimum feature size equal
> to the width of
> - # the quadrilateral. This one is kinda tricky.
> - # v1,v2 = vectors along edges
> - v1 = array([0.5*(s[1]-s[0] + s[2]-s[3]), 0.5*(t[1]-t[0]
> + t[2]-t[3]),])
> - v2 = array([0.5*(s[3]-s[0] + s[2]-s[1]), 0.5*(t[3]-t[0]
> + t[2]-t[1]),])
> - v1Sq = dot(v1,v1)
> - v2Sq = dot(v2,v2)
> - level = 0.5*log2(min(v1Sq,v2Sq) * (1 -
> dot(v1,v2)**2/(v1Sq*v2Sq)))
> - elif self.levelCalcMethod == 'minDiag':
> - # Get mipmap level with minimum feature size equal
> to the minimum
> - # distance between the centre of the quad and the
> vertices. Sort
> - # of a "quad radius"
> - #
> - # This is more-or-less the algorithm used in
> Pixie...
> - minDiag2 = ((s - s.mean())**2 + (t - t.mean
> ())**2).min()
> - level = log2(minDiag2)/2
> - #elif self.levelCalcMethod == 'sqrtArea':
> - # Get mipmap level with minimum feature size
> estimated as the
> - # square root of the area of the box.
> - elif self.levelCalcMethod == 'trilinear':
> - # Get mipmap level which will result in no
> aliasing when plain
> - # trilinear filtering is used (no integration)
> - maxDiag2 = ((s - s.mean())**2 + (t - t.mean
> ())**2).max()
> - level = log2(maxDiag2)/2
> - elif self.levelCalcMethod == 'level0':
> - # Else just use level 0. Correct texture
> filtering will take care
> - # of any aliasing...
> - level = 0
> - else:
> - raise "Invalid mipmap level calculation type: %s"
> % self.levelCalcMethod
> - return max(level,0)
> -
> - def filterQuadAF(self, s, t):
> - '''
> - Filter the texture by sampling inside the given
> quadriateral specified
> - by the texture-space vertices
> - [s[1], t[1]], ... , [s[4], t[4]],
> - and taking an average. Sample points are chosen
> stochastically over
> - the quadrilateral at the mipmap level returned by
> calculateLevel.
> - This is essentially Monte Carlo integration, with a
> reconstruction
> - filter equal to bilinear interpolation.
> - '''
> - level = int(self.calculateLevel(s,t))
> - return self.levels[level].filterQuadAF(s,t,
> self.numSamples)
> -
> - def filterEWA(self, sBox, tBox, deltaX=None, deltaY=None, J=None):
> - '''
> - Filter a portion of the texture using an elliptical
> gaussian filter, as
> - described in Heckbert's 1989 thesis, "Fundamentals of
> Texture Mapping
> - and Image Warping". This kind of filter is also called an
> EWA filter for
> - "Elliptical Weighted Average".
> -
> - s, t - Texture box over which to filter
> - deltaX, deltaY - pixel side length in the final discrete
> image.
> - J - Jacobian of the coordinate transformation.
> - '''
> - level = int(self.calculateLevel(sBox,tBox))
> - col = self.levels[level].filterEWA(sBox, tBox, deltaX,
> deltaY, J)
> - return col #* 0.9**level
> -
> - def displayLevel(self, level, s=None, t=None):
> - '''
> - Display a mipmap level
> - '''
> - self.levels[level].display(s,t)
> -
> -
> -#----------------------------------------------------------------------
> -#----------------------------------------------------------------------
> -# 2D texture transformation classes ("warps")
> -#
> -# These are useful as example texture warps to test anisotropic filtering
> -# techniques.
> -#
> -# Each warp class defines three functions
> -# fwd - the forward transform
> -# inv - the inverse transform
> -# jacobian - the jacobian of the forward transform
> -#----------------------------------------------------------------------
> -class PerspectiveWarp:
> - '''
> - Simplified perspective transformation, between "texture
> coordinates", (s,t)
> - on an infinite plane with parametric form P(s,t) = [s,-1,t],
> - and "image coordinates", (x,y) in the x,y plane.
> -
> - We use the LH coordinate system, with x=right, y=up, z=into the
> screen.
> - '''
> - def __init__(self, d=1):
> - '''
> - d - distance from view window to viewer
> - '''
> - self.d = d
> -
> - def fwd(self, s, t):
> - '''
> - Forward perspective transformation
> -
> - Params:
> - s,t - texture coordinates along x & z directions
> -
> - Returns:
> - (x,y) - coordinates on the view window.
> - '''
> - x = 1/(1+t/self.d) * s
> - y = -1/(1+t/self.d)
> - return (x,y)
> -
> - def jacobian(self, s, t):
> - '''
> - return the Jacobian matrix of the forward transformation
> - '''
> - D = 1/self.d*(1+t/self.d)**(-2)
> - return D*array([[self.d+t, -s], [0, 1]])
> -
> - def inv(self, x, y):
> - '''
> - Inverse perspective transformation
> -
> - Params:
> - x,y - coordinates in the view window
> -
> - Returns:
> - (s,t) - texture coordinates along x & z directions
> - '''
> - s = -x/y
> - t = -self.d*(1+1/y)
> - return (s,t)
> -
> -#----------------------------------------------------------------------
> -class AffineWarp:
> - '''
> - A class representing an affine transformation
> - '''
> - def __init__(self, A=eye(2), b=array([0,0]), theta=None,
> scale=None):
> - '''
> - Affine transformation is x' = A*x + b
> - A - linear part of the affine trans.
> - b - translation
> - theta - if this is set to a number, make A into a rotation
> with angle theta.
> - scale - if this is set to a number, scale A by this
> factor.
> - '''
> - self.A = asarray(A)
> - self.b = asarray(b)
> - if theta is not None:
> - self.theta = theta
> - theta *= pi/180
> - self.A = array([[cos(theta), sin(theta)],
> [-sin(theta), cos(theta)]])
> - if scale is not None:
> - self.A *= scale
> - self.Ai = linalg.inv(self.A)
> -
> - def fwd(self, s, t):
> - return (self.A[0,0]*s + self.A[0,1]*t + self.b[0], \
> - self.A[1,0]*s + self.A[1,1]*t + self.b[1])
> -
> - def inv(self, x, y):
> - return (self.Ai[0,0]*(x-self.b[0]) + self.Ai[0,1]*(y -
> self.b[1]), \
> - self.Ai[1,0]*(x-self.b[0]) + self.Ai[1,1]*(y
> - self.b[1]))
> -
> - def jacobian(self, s, t):
> - return self.A
> -
> - def __mul__(self, rhs):
> - '''
> - The * operator composes affine transformations
> - '''
> - return AffineWarp(A=dot(self.A,rhs.A), b=dot(self.A,rhs.b
> )+self.b)
> -
> -#----------------------------------------------------------------------
> -class CompositionWarp:
> - '''
> - A class for holding the composition of several transformations
> - '''
> - def __init__(self, *args):
> - '''
> - Takes a list of transformations to be composed.
> - The transformation list is evaluated from right to left in
> the same
> - manner as it is usually written mathematically.
> - '''
> - self.transList = list(args)
> - self.transList.reverse()
> -
> - def fwd(self, s,t):
> - x,y = s,t
> - for trans in self.transList:
> - x,y = trans.fwd(x,y)
> - return (x,y)
> -
> - def inv(self, x,y):
> - s,t = x,y
> - for trans in reversed(self.transList):
> - s,t = trans.inv(s,t)
> - return (s,t)
> -
> - def jacobian(self, s,t):
> - J = eye(2)
> - for trans in self.transList:
> - J = dot(trans.jacobian(s,t),J)
> - s,t = trans.fwd(s,t)
> - return J
> -
> -#----------------------------------------------------------------------
> -def plotPerspectiveWarp(u, v, direc='f'):
> - '''
> - Plot u,v on an axis, and the things they transform into under an
> (possibly
> - inverse) perspective transformation on another axis.
> -
> - u,v - either image coords if direc='b' or texture coords if
> direc='f'
> - direc - transformation direction; 'f' = forward, 'b'=backward
> - '''
> - perspec = PerspectiveWarp()
> - if direc == 'f':
> - s,t = u,v
> - x,y = perspec.fwd(s,t)
> - else:
> - x,y = u,v
> - s,t = perspec.inv(x,y)
> -
> - # plot texture coords
> - pylab.subplot(211)
> - pylab.plot(s,t)
> - pylab.axis('equal')
> -
> - # plot image coords
> - pylab.subplot(212)
> - pylab.plot(x,y)
> - pylab.axis('equal')
> - pylab.axis([-1,1,-1,0])
> -
> -#----------------------------------------------------------------------
> -# High level driver functions.
> -#----------------------------------------------------------------------
> -def genTrilinearSeq(mipmap, s, t, ds, dt):
> - '''
> - Generate and save a sequence of images by trilinear interpolation
> at
> - various levels, but with the same texture coordinates (s,t)
> - '''
> - for i in range(0,31):
> - level = i/31.0*8
> - image = mipmap.sampleTrilinear(s,t,ds,dt,level)
> - imsave('downsamp%0.2d.png' % i, image)
> -
> -#----------------------------------------------------------------------
> -def plotPixelBoxes():
> - '''
> - Plot the outlines of some "pixels" in the image plane, as they
> transform to
> - the texture coordinates. This allows us to see what kinds of
> filter shapes
> - we *should* have to do proper texture filtering... Naive
> trilinear
> - mipmapping only gets square filters correct...
> - '''
> - x = array([1,1,-1,-1,1])
> - y = array([-1,1,1,-1,-1])
> - plotPerspectiveWarp(x,y+1,'f')
> - plotPerspectiveWarp(x,y+7,'f')
> - plotPerspectiveWarp(x,y+4,'f')
> - plotPerspectiveWarp(x+3,y+4,'f')
> - plotPerspectiveWarp(x*0.05,y*0.05-0.1,'b')
> - plotPerspectiveWarp(x*0.01,y*0.01-0.1,'b')
> - plotPerspectiveWarp(x*0.01-0.5,y*0.01-0.1,'b')
> - plotPerspectiveWarp(x*0.01+0.5,y*0.01-0.1,'b')
> - plotPerspectiveWarp(x*0.01+0.5,y*0.01-0.05,'b')
> -
> -#----------------------------------------------------------------------
> -def textureWarp(x, y, tex, trans=PerspectiveWarp(), method='EWA'):
> - '''
> - Warp the given texture.
> -
> - x,y - image space coordinates for output.
> - '''
> - im = zeros((x.shape[1]-1, x.shape[0]-1, 3))
> - deltaX = x[1,0]-x[0,0]
> - deltaY = y[0,1]-y[0,0]
> - print "Warping the texture..."
> - for i in range(0,x.shape[0]-1):
> - for j in range(0,x.shape[1]-1):
> - xBox = array((x[i,j], x[i+1,j], x[i+1,j+1],
> x[i,j+1]))
> - yBox = array((y[i,j], y[i+1,j], y[i+1,j+1],
> y[i,j+1]))
> - sBox, tBox = trans.inv(xBox, yBox)
> - if method == 'EWA':
> - # use Elliptical Weighted Average filter,
> with Jacobian
> - # estimation via finite differences
> - col = tex.filterEWA(sBox, tBox)
> - elif method == 'EWA_analytical_J':
> - # Use EWA with an analytical
> Jacobian. This can be used as a
> - # check against the numerical
> estimate. Last time I looked, it
> - # was spot-on!
> - J = trans.jacobian(sBox.mean(),tBox.mean
> ())
> - col = tex.filterEWA(sBox, tBox, deltaX,
> deltaY, J)
> - elif method == 'MC_integration':
> - # Use Monte Carlo integration over a box.
> - col = tex.filterQuadAF(sBox, tBox)
> - elif method == 'trilinear':
> - # Use plain old trilinear filtering
> - col = tex.filterTrilinear(sBox, tBox)
> - else:
> - raise "Invalid method specified."
> - im[-1-j,i,:] = col
> - if i % 20 == 0:
> - print "done col %d, " % i
> - return im
> -
> -#----------------------------------------------------------------------
> -def showMipmapOffsetDifferences():
> - '''
> - Show the differences in warped images generated when the mipmap
> offsets are
> - turned off (very small).
> - '''
> - im = getEdgeTest(1000, 100)
> - texture = TextureMap(im, 8)
> -
> - N = 1000
> - x,y = mgrid[-0.5:0.5:N*1j,-0.22:-0.005:N*2//5*1j]
> - rot = AffineWarp(b=(0.5,0.5)) * AffineWarp(theta=45) *
> AffineWarp(b=(-0.5,-0.5))
> - trans = CompositionWarp(PerspectiveWarp(0.4),AffineWarp(A=diag((10,30)),
> b = array((-5,4)))*rot)
> - imWarp = textureWarp(x,y, texture, trans)
> - imsave('tmp_offset_bad.png',imWarp)
> -
> -#----------------------------------------------------------------------
> -def showLevelCalcDifferences():
> - '''
> - Compare several of the best mipmap level calculation methods for
> EWA filters.
> - '''
> - im = getGridImg(1000, 30, 2)
> - texture = TextureMap(im, 8, levelCalcMethod='minQuadWidth')
> -
> - N = 500
> - x,y = mgrid[-0.5:0.5:N*1j,-0.22:-0.005:N*2//5*1j]
> - trans =
> CompositionWarp(PerspectiveWarp(),AffineWarp(A=diag((10,30)), b =
> array((-5,4))))
> - imWarp = textureWarp(x,y, texture, trans)
> - imshow(imWarp)
> - imsave('level_minQuadWidth.png', imWarp)
> -
> - texture.levelCalcMethod='minSideLen'
> - imWarp = textureWarp(x,y, texture, trans)
> - imshow(imWarp)
> - imsave('level_minSideLen.png', imWarp)
> -
> - texture.levelCalcMethod='level0'
> - imWarp = textureWarp(x,y, texture, trans)
> - imshow(imWarp)
> - imsave('level_level0.png', imWarp)
> -
> -#----------------------------------------------------------------------
> -def showFilteringMethodDifferences():
> - '''
> - Compare several different filtering methods
> - '''
> - im = getGridImg(1000, 30, 2)
> - texture = TextureMap(im, 8, levelCalcMethod='minQuadWidth')
> -
> - N = 700
> - x,y = mgrid[-0.5:0.5:N*1j,-0.22:-0.005:N*2//5*1j]
> - trans =
> CompositionWarp(PerspectiveWarp(),AffineWarp(A=diag((10,30)), b =
> array((-5,4))))
> - imWarp = textureWarp(x,y, texture, trans, method='EWA')
> - imshow(imWarp)
> - imsave('method_EWA.png', imWarp)
> - imWarp = textureWarp(x,y, texture, trans,
> method='EWA_analytical_J')
> - imshow(imWarp)
> - imsave('method_EWA_analytical_J.png', imWarp)
> - imWarp = textureWarp(x,y, texture, trans, method='MC_integration')
> - imshow(imWarp)
> - imsave('method_MC_integration.png', imWarp)
> - texture.levelCalcMethod='trilinear' # need a special level
> calculation for trilin.
> - imWarp = textureWarp(x,y, texture, trans, method='trilinear')
> - imshow(imWarp)
> - imsave('method_trilinear.png', imWarp)
> -
> -#----------------------------------------------------------------------
> -# Main program
> -#----------------------------------------------------------------------
> -
> -if __name__ == '__main__':
> - # a usage example
> - if True:
> - # Choose image & create mipmap
> - #im = getCentreTest(1000)
> - im = getGridImg(1000, 30, 2)
> - #im = getEdgeTest(1000, 100)
> - #im = imread('lena_std.png')[:,:,:3]/256.0
> - texture = TextureMap(im, 8,
> levelCalcMethod='minQuadWidth')
> - # Warp image to new coordinate system.
> - N = 500
> - x,y = mgrid[-0.5:0.5:N*1j,-0.22:-0.005:N*2//5*1j]
> - trans =
> CompositionWarp(PerspectiveWarp(),AffineWarp(A=diag((10,30)), b =
> array((-5,4))))
> - imWarp = textureWarp(x,y, texture, trans)
> - pylab.hold(False)
> - imshow(imWarp)
> - imsave('imWarp.png',imWarp)
> - pylab.draw()
> - else:
> - # Area for hacking out solutions.
> - #showLevelCalcDifferences()
> - showFilteringMethodDifferences()
>
> Copied: branches/soc2007-dsm/aqsis/testing/prototypes/texfilt/texfilt.py
> (from rev 1048, trunk/aqsis/testing/prototypes/texfilt/texfilt.py)
> ===================================================================
> ---
> branches/soc2007-dsm/aqsis/testing/prototypes/texfilt/texfilt.py (rev
> 0)
> +++
> branches/soc2007-dsm/aqsis/testing/prototypes/texfilt/texfilt.py 2007-05-19
> 02:48:41 UTC (rev 1049)
> @@ -0,0 +1,893 @@
> +#!/usr/bin/python
> +#
> +# Scientific python script to test image filtering for mipmap creation in
> Aqsis
> +# Copyright (C) 2007, Christopher J. Foster
> +#
> +# This program is free software; you can redistribute it and/or
> +# modify it under the terms of the GNU General Public License
> +# as published by the Free Software Foundation; either version 2
> +# of the License, or (at your option) any later version.
> +#
> +# This program is distributed in the hope that it will be useful,
> +# but WITHOUT ANY WARRANTY; without even the implied warranty of
> +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
> +# GNU General Public License for more details.
> +#
> +# You should have received a copy of the GNU General Public License
> +# along with this program; if not, write to the Free Software
> +# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
> MA 02110-1301, USA.
> +
> +'''
> +This script investigates texture access. There are two main parts:
> + 1) Mipmap generation & access
> + 2) Texture filtering using anisotropic filtering. This includes
> + * stochastic filtering with a weight function over a
> quadrilateral
> + * elliptical gaussian filtering (EWA)
> +
> +To run the script, you'll need scientific python; specifically, the
> packages
> + * scipy
> + * matplotlib
> +On a unix machine, ipython is highly recommended for a nice interactive
> +interface which meshes very well with matplotlib.
> +'''
> +
> +from __future__ import division
> +
> +import matplotlib.pylab as pylab
> +from pylab import imshow
> +
> +import numpy
> +from numpy import pi, r_, size, zeros, ones, eye, ceil, floor, array, \
> + linspace, meshgrid, random, arange, dot, diag, asarray
> +
> +import scipy
> +from scipy import cos, sin, exp, log, log2, sqrt, ogrid, mgrid, ndimage
> +import scipy.linalg as linalg
> +from scipy.misc import imsave, imread
> +
> +# The purpose of this script is to investigate mipmap generation and
> indexing
> +# in more-or-less the same way which it will be used internally by Aqsis.
> +
> +#----------------------------------------------------------------------
> +#----------------------------------------------------------------------
> +# Filter Kernels, including any windowing.
> +#----------------------------------------------------------------------
> +def sincKer(width, even, scale = 2):
> + '''
> + Return a sinc filter kernel with specified width and scale
> +
> + width - width of the filter (used by the window)
> + even - true if the generated filter should have an even number of
> points.
> + scale - zeros will be centred around 0 with spacing "scale"
> + '''
> + if even:
> + numPts = max(2*int((width+1)/2), 2)
> + else:
> + numPts = max(2*int(width/2) + 1, 3)
> + x = r_[0:numPts]-(numPts-1)*0.5
> + k = scipy.sinc(x/scale)
> + k *= cos(linspace(-pi/4,pi/4,len(k)))
> + ker = numpy.outer(k,k)
> + ker /= sum(ker.ravel())
> + return ker
> +
> +
> +#----------------------------------------------------------------------
> +# Image manipuation
> +#----------------------------------------------------------------------
> +def applyToChannels(img, fxn, *args, **kwargs):
> + '''
> + Apply the given function to each of the channels in img in turn.
> + Additional positional and keyword arguments are passed to the
> called fxn.
> + '''
> + numChannels = img.shape[-1]
> + outList = []
> + for i in range(0,numChannels):
> + outList.append(fxn(img[:,:,i], *args, **kwargs))
> + outArr = array(outList)
> + return outArr.transpose(range(1,len(outArr.shape)) + [0])
> +
> +
> +#----------------------------------------------------------------------
> +#----------------------------------------------------------------------
> +# Functions for creating test images
> +#----------------------------------------------------------------------
> +def getGridImg(N, gridSize, gridWidth=1):
> + '''
> + Get a test image containig a one pixel wide grid
> +
> + N - side length
> + gridSize - size of the one-pixel wide black grid in the image.
> + '''
> +
> + im = zeros((N, N, 3))
> +
> + x = linspace(0,1,N)
> + [xx,yy] = meshgrid(x,x)
> + im[:,:,0] = xx
> + im[:,:,1] = 0.5
> + im[:,:,2] = yy
> +
> + for i in range(0,gridWidth):
> + im[i:-gridSize//2:gridSize,:,:] = 0
> + im[:,i:-gridSize//2:gridSize,:] = 0
> +
> + im[:,-1,:] = 0
> + im[-1,:,:] = 0
> +
> + return im
> +
> +
> +def getCentreTest(N=512):
> + '''
> + Get a test image consisting of a 2x2 checkered plane.
> +
> + N - side length
> + '''
> + im = zeros((N, N, 3))
> +
> + im[:N/2,N/2:,:] = 1
> + im[N/2:,:N/2,:] = 1
> +
> + return im
> +
> +
> +def getEdgeTest(N=512, borderwidth=20):
> + '''
> + Get a white test image with a black bordering strip
> +
> + N - side length
> + '''
> + im = ones((N, N, 3))
> +
> + im[borderwidth:2*borderwidth,:,:] = 0
> + im[-2*borderwidth:-borderwidth,:,:] = 0
> + im[:,borderwidth:2*borderwidth,:] = 0
> + im[:,-2*borderwidth:-borderwidth,:] = 0
> +
> + return im
> +
> +#----------------------------------------------------------------------
> +#----------------------------------------------------------------------
> +# Utils for creating/manipulating mipmaps
> +#----------------------------------------------------------------------
> +
> +class ImageDownsampler:
> + '''
> + Class which holds a filter kernel (function & width) and is able
> to
> + correctly use it to downsample odd or even-sized images.
> + '''
> + def __init__(self, kerWidth, kerFunc=sincKer):
> + self.kerFunc = kerFunc
> + self.kerWidth = kerWidth
> + def downsample(self, img):
> + '''
> + Downsample an image by
> + (1) Applying the filter kernel at each pixel
> + (2) Decimate by taking every second element from the
> result
> + (3) Clamp the resulting intensities between 0 and 1
> + '''
> + ker = self.kerFunc(self.kerWidth, (size(img,0) % 2) == 0)
> + res = applyToChannels(img, ndimage.correlate, ker,
> mode='constant', cval=0)
> + res = res[(1-(size(img,0) % 2))::2, (1-(size(img,1) %
> 2))::2, :]
> + res = res.clip(min=0, max=1)
> + return res
> +
> +#----------------------------------------------------------------------
> +class MipLevel:
> + '''
> + Class to hold a level of mipmap, which knows how to bilinear
> interpolation
> + on itself.
> + '''
> + def __init__(self, image, offset, mult):
> + '''
> + image - pixel data; a NxMx3 array
> + offset, mult - length two vectors which specify an affine
> + transformation from texture coordinates to this
> mipmap's raster
> + coordinates:
> + i = mult[0] * s + offset[0]
> + j = mult[1] * t + offset[1]
> + '''
> + # store the image with some ghost points to work around a
> bug in
> + # scipy.ndimage, where incorrect interpolation happens for
> negative
> + # indices.
> + self.image = zeros((image.shape[0]+2,image.shape[1]+2,
> image.shape[2]))
> + self.image[1:-1,1:-1,:] = image
> + self.offset = array(offset,'d')
> + self.mult = array(mult,'d')
> + # Variances for the reconstruction & prefilters when using
> EWA. A
> + # variance of 1/(2*pi) gives a filter with centeral weight
> 1, but in
> + # practise this is slightly too small (resulting in a
> little bit of
> + # aliasing). Therefore it's adjusted up slightly.
> + self.varianceXY = 1.3/(2*pi)
> + self.varianceST = 1.3/(2*pi)
> +
> + def rasterCoords(self,s,t):
> + '''
> + Return the raster coordinates for this mipmap from the
> texture
> + coordinates (s,t)
> + '''
> + i = scipy.atleast_1d(self.mult[0]*s + self.offset[0] + 1)
> + j = scipy.atleast_1d(self.mult[1]*t + self.offset[1] + 1)
> + return (i,j)
> +
> + def sampleBilinear(self, s, t):
> + '''
> + Sample the mipmap level at the coordinates specified by s
> and t. For
> + performance, s & t can be 1D arrays.
> + '''
> + i,j = self.rasterCoords(s,t)
> + return applyToChannels(self.image, ndimage.map_coordinates,
> \
> + array([i,j]), order=1, mode='constant',
> cval=0)
> +
> + def getRandPointsInQuad(self, s, t, numPoints):
> + '''
> + Get a random set of points inside the quadrilatiral given
> by the
> + texture-space vertices
> + [s[1], t[1]], ... , [s[4], t[4]].
> +
> + The sample should be uniform if the quadrialteral can be
> made from a
> + square with a _linear_ transformation. Otherwise it will
> be biased.
> + The form of biasing is probably exactly what we want if
> we're taking
> + pixel samples...
> + '''
> + lerp1 = numpy.random.rand(numPoints)
> + lerp1m = 1-lerp1
> + lerp2 = numpy.random.rand(numPoints)
> + st = array([s,t])
> + op = numpy.outer
> + randPts = lerp2*(op(st[:,0],lerp1) + op(st[:,1],lerp1m)) \
> + + (1-lerp2)*(op(st[:,3],lerp1) +
> op(st[:,2],lerp1m))
> + #weights = exp(-8*((lerp1-0.5)**2 + (lerp2-0.5)**2))
> + weights = ones(numPoints)
> + return (randPts[0,:], randPts[1,:], weights)
> +
> + def filterQuadAF(self, s, t, numSamples):
> + '''
> + Filter the texture by sampling inside the given
> quadriateral specified
> + by the texture-space vertices
> + [s[1], t[1]], ... , [s[4], t[4]],
> + and taking an average. Sample points are chosen
> stochastically over
> + the quadrilateral. This is essentially Monte Carlo
> integration, with a
> + reconstruction filter equal to bilinear interpolation.
> + '''
> + sRnd, tRnd, w = self.getRandPointsInQuad(s, t, numSamples)
> + samples = self.sampleBilinear(sRnd, tRnd)
> + return ( samples*w.reshape(w.shape+(1,)) ).sum(0) / w.sum
> ()
> +
> + def estimateJacobianInverse(self, s, t):
> + '''
> + Estimate the inverse Jacobian from a small (s,t) box using
> finite
> + differences. The assumption is that the (s,t) box has
> come from a 1x1
> + box in the final coordinate system.
> + '''
> + return 0.5 * array([[s[1]-s[0] + s[2]-s[3], s[3]-s[0] +
> s[2]-s[1]],
> + [t[1]-t[0] +
> t[2]-t[3], t[3]-t[0] + t[2]-t[1]]])
> +
> + def filterEWA(self, sBox, tBox, deltaX=None, deltaY=None, J=None,
> returnWeights=False):
> + '''
> + Filter a portion of the texture using an elliptical
> gaussian filter, as
> + described in Heckbert's 1989 thesis, "Fundamentals of
> Texture Mapping
> + and Image Warping". This kind of filter is also called an
> EWA filter for
> + "Elliptical Weighted Average".
> +
> + s, t - texture coordinates of a box to be filtered over
> + deltaX, deltaY - pixel side length in the final discrete
> image.
> + J - Jacobian of the coordinate transformation.
> + returnWeights - true if you just want the filter weights
> returned for
> + debugging rather than the
> filtered image.
> + '''
> + # convert to raster coords
> + sBox,tBox = self.rasterCoords(sBox,tBox)
> + # centre of box
> + s0,t0 = sBox.mean(), tBox.mean()
> + if not returnWeights and (s0 < 0 or s0 > self.image.shape[0]
> \
> + or t0 < 0
> or t0 > self.image.shape[1]):
> + return 0
> + if J is None:
> + # Estimate the Jacobian from the texture box.
> + Ji = self.estimateJacobianInverse(sBox, tBox)
> + else:
> + # Else we have an analytical Jacobian specified.
> + # Adjust J for raster->texture coords
> + J = dot(dot(diag((1/deltaX, 1/deltaY)), J),
> diag((1/self.image.shape[0], 1/self.image.shape[1])))
> + # get the inverse of J
> + Ji = linalg.inv(J)
> + # width parameter is log(n) where the filter falls to 1/n
> of it's
> + # centeral value at the edge.
> + width = 4
> + # V = covariance matrix
> + V = self.varianceXY*dot(Ji, Ji.transpose()) +
> self.varianceST*eye(2)
> + # Q = (+ve definite) quadratic form matrix.
> + Q = 0.5*linalg.inv(V)
> + det = Q[0,0]*Q[1,1] - Q[0,1]*Q[1,0]
> + # compute bounding radii
> + sRad = sqrt(Q[1,1]*width/det)
> + tRad = sqrt(Q[0,0]*width/det)
> + # integer filter width
> + sWidth = int(floor(s0+sRad)) - int(ceil(s0-sRad)) + 1
> + tWidth = int(floor(t0+tRad)) - int(ceil(t0-tRad)) + 1
> + # compute edge of filter region in level0 image:
> + sStart = ceil(s0-sRad)
> + tStart = ceil(t0-tRad)
> + # compute filter weights. I'm being sloppy here, and not
> bothering to
> + # truncate outside the relevant ellipse
> + s,t = ogrid[0:sWidth, 0:tWidth]
> + s = s - (s0-sStart)
> + t = t - (t0-tStart)
> + if not returnWeights and (sStart < 0 or sStart + sWidth-1
> >= self.image.shape[0] \
> + or tStart < 0 or tStart + tWidth-1 >=
> self.image.shape[1]):
> + return 0.5
> + weights = exp(-(Q[0,0]*s*s + (Q[0,1]+Q[1,0])*s*t +
> Q[1,1]*t*t))
> + weights /= weights.sum()
> + if returnWeights:
> + return weights
> + imSec = self.image
> [sStart:sStart+sWidth,tStart:tStart+tWidth]
> + return applyToChannels(imSec, lambda x: (x*weights).sum())
> +
> + def display(self, s=None, t=None):
> + if s is None or t is None:
> + imshow(self.image[1:-1,1:-1,:],
> interpolation='nearest')
> + else:
> + imshow(self.sampleBilinear(s,t),
> interpolation='nearest')
> +
> +
> +#----------------------------------------------------------------------
> +class TextureMap:
> + '''
> + Class to encapsulate a (power-of-2) mipmapped texture in much the
> same way
> + mipmapping will be done internally in Aqsis.
> + '''
> + def __init__(self, img, filterWidth, kerFunc=sincKer, \
> + levelCalcMethod='minQuadWidth', numSamples=100):
> + '''
> + img - texture image
> + filterWidth - filter width for mipmap downsampling
> + kerFunc - kernel function for mipmap downsampling
> + levelCalcMethod - method used for calculating required
> mipmap level
> + from a quadrialteral. filtering. May be one of:
> 'minSideLen',
> + 'minQuadWidth', 'trilinear', 'minDiag', 'level0'.
> + numSamples - number of samples to used when performing
> stochastic
> + filtering
> + '''
> + self.downsampler = ImageDownsampler(filterWidth, kerFunc)
> + self.levels = []
> + self.levelCalcMethod = levelCalcMethod
> + self.genMipMap(img)
> + self.numSamples = numSamples
> +
> + def genMipMap(self, img):
> + '''
> + Generate a mipmap, by downsampling by successive powers of
> two, using
> + the held kernel function and width.
> +
> + In this function we make the requirement that the filter
> coefficients
> + for filtering between levels are constant across the image
> for a
> + particular mipmap level. This restriction means that
> even-sized images
> + have to be filtered with even- sized filter kernels, and
> vice versa.
> + Even sizes also require that an offset be taken into
> account, since the
> + edges of the downsampled levels don't lie on the edges of
> the original
> + image.
> +
> + Examples showing mipmap sample points (represented as x's)
> follow for
> + the 1D case. Several variables are also shown as a
> function of the
> + mipmap level, l on the rhs:
> +
> + i_l(t) = raster coordinates for level l as a function of
> texture coordinate t.
> + o_l = offset for level l in level 0 units
> + s_l = span of level l in level 0 units
> + d_l = interval between pixels of level l in level 0 raster
> units.
> + w_l = width of level l in pixels.
> +
> +
> + (mostly) odd numbers of points:
> + 0 1 2 3 4 <-- i_0
> + l | | | |
> | o_l s_l w_l d_l
> + 0 [x x x x x] 0 4 5 1
> + 1 [x x x] 0 2 3 2
> + 2 [x x] 0 1 2 4
> + 3 [ x ] 2 0 1 ?
> +
> + Another example (powers of two)
> + 0 1 2 3 4 5 6 7 <-- i_0
> + l | | | | | | | | o_l s_l w_l
> + 0 [x x x x x x x x] 0 7 8
> + 1 [ x x x x ] 0.5 6 4
> + 2 [ x x ] 1.5 4 2
> + 3 [ x ] 3.5 0 1
> +
> + Notice that the higher level mipmaps for the power of
> two-sized image
> + do *not* have points which lie on the edge of the
> image. This means
> + that we've got to incorporate an offset when we sample the
> level.
> +
> + Formulas for texture coordinates -> raster space
> +
> + s_l = (s_0 - 2*o_l)
> +
> + Raster units for level 0 are:
> + i_0(t) = t * s_0
> + Raster units for level l are:
> + i_l(t) = (i_0(t) - o_0)/d_l
> +
> + Or, what we really need, which is raster units for level l
> in terms of
> + the texture coordinate t:
> + i_l(t) = t * s_0/d_l - o_0/d_l
> +
> + ((
> + Note that we could abandon the
> constant-coefficients approach
> + (which would make it very slow to generate
> mipmaps), we would have:
> +
> + 0 1 2 3 4 5 6 7 <-- i_0
> + l | | | | | | | | o_l s_l w_l
> + 0 [x x x x x x x x] 0 7 8
> + 1 [x x x x] 0 7 4
> + 2 [x x] 0 7 2
> + 3 [ x ] 0 7 1
> +
> + The the raster units for level l in terms of the
> texture coordinates t
> + are then very simple:
> +
> + i_l(t) = t*(w_l - 1)
> + ))
> + '''
> + scale = 2
> + level0Span = array(img.shape[0:2])-1
> + offset = array([0,0],'d')
> + self.levels = [MipLevel(img, offset, level0Span)]
> + # delta is distance between pixels; level 0 units
> + delta = 1
> + # offsets are between 1st pixel of 0th level to 1st pixel
> of current level
> + #
> + # eg, in 1D
> + # [x x x x..]
> + # <---> <- level 0 delta
> + # [ x x ...]
> + # <-> <- level 1 offset
> + while (array(img.shape[0:2]) > 1).all():
> + # calculate start offsets for this level.
> + offset += (1-(array(img.shape[0:2]) % 2))*delta/2
> + delta *= 2
> + # downsample
> + img = self.downsampler.downsample(img)
> + # create a new MipLevel to hold the image.
> + self.levels.append(MipLevel(img, -offset/delta,
> level0Span/delta))
> + # The below is a theoretically broken, but very
> usable version, -
> + # mipmap offsets are set to zero.
> + #self.levels.append(MipLevel(img, 0*offset, array(
> img.shape[0:2])-1))
> + print 'generated level size', img.shape[0:2]
> +
> + def filterTrilinear(self, sBox, tBox, level=None):
> + '''
> + Trilinear filtering over the mipmaps
> +
> + s, t - texture coordinates of points to sample at
> + ds,dt - size of a "box" on the screen.
> + level - mipmap level to use. If equal to None, calculate
> from ds & dt
> + '''
> + if level is None:
> + level = min(len(self.levels)-1,
> self.calculateLevel(sBox, tBox))
> + level = max(level,0)
> + level1 = int(floor(level))
> + interp = level - level1
> + level2 = min(len(self.levels)-1, int(level1 + 1))
> + s0 = sBox.mean()
> + t0 = tBox.mean()
> + samp1 = self.levels[level1].sampleBilinear(s0,t0)
> + samp2 = self.levels[level2].sampleBilinear(s0,t0)
> + return (1-interp)*samp1 + interp*samp2
> +
> + def calculateLevel(self, s, t):
> + '''
> + Calculate the appropriate mipmap level for texture
> filtering over a
> + quadrilateral given by the texture-space vertices
> + [s[1], t[1]], ... , [s[4], t[4]].
> +
> + There are many ways to do this; the instance variable
> levelCalcMethod
> + selects the desired one. The most correct way is to
> choose the
> + minSideLen method, as long as the quadrilateral is vaguely
> + rectangular-shaped. This only works if you're happy to
> use lots of
> + samples however, otherwise you get aliasing.
> + '''
> + s = s.copy()*self.levels[0].image.shape[0]
> + t = t.copy()*self.levels[0].image.shape[1]
> + if self.levelCalcMethod == 'minSideLen':
> + # Get mipmap level with minimum feature size equal
> to the shortest
> + # quadrilateral side
> + s1 = pylab.concatenate((s, s[0:1]))
> + t1 = pylab.concatenate((t, t[0:1]))
> + minSideLen2 = (numpy.diff(s1)**2 + numpy.diff
> (t1)**2).min()
> + level = log2(minSideLen2)/2
> + elif self.levelCalcMethod == 'minQuadWidth':
> + # Get mipmap level with minimum feature size equal
> to the width of
> + # the quadrilateral. This one is kinda tricky.
> + # v1,v2 = vectors along edges
> + v1 = array([0.5*(s[1]-s[0] + s[2]-s[3]), 0.5*(t[1]-t[0]
> + t[2]-t[3]),])
> + v2 = array([0.5*(s[3]-s[0] + s[2]-s[1]), 0.5*(t[3]-t[0]
> + t[2]-t[1]),])
> + v1Sq = dot(v1,v1)
> + v2Sq = dot(v2,...
>
> [Message clipped]
--
Paul Gregory
http://www.aqsis.org
|