commit 7cdac8ded95ea4ff0ef488e975b30b262ade5612 Author: Pascal Engélibert Date: Sat Jul 26 22:14:37 2025 +0200 Initial commit diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..6aa15dc --- /dev/null +++ b/.gitignore @@ -0,0 +1,3 @@ +build/* +html +*.exr diff --git a/APSF.cpp b/APSF.cpp new file mode 100644 index 0000000..74b2fbf --- /dev/null +++ b/APSF.cpp @@ -0,0 +1,272 @@ +/////////////////////////////////////////////////////////////////////////////////// +// File : APSF.cpp +/////////////////////////////////////////////////////////////////////////////////// +// +// LumosQuad - A Lightning Generator +// Copyright 2007 +// The University of North Carolina at Chapel Hill +// +/////////////////////////////////////////////////////////////////////////////////// +// +// 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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// The University of North Carolina at Chapel Hill makes no representations +// about the suitability of this software for any purpose. It is provided +// "as is" without express or implied warranty. +// +// Permission to use, copy, modify and distribute this software and its +// documentation for educational, research and non-profit purposes, without +// fee, and without a written agreement is hereby granted, provided that the +// above copyright notice and the following three paragraphs appear in all +// copies. +// +// THE UNIVERSITY OF NORTH CAROLINA SPECIFICALLY DISCLAIM ANY WARRANTIES, +// INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND +// FITNESS FOR A PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN +// "AS IS" BASIS, AND THE UNIVERSITY OF NORTH CAROLINA HAS NO OBLIGATION TO +// PROVIDE MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS. +// +// Please send questions and comments about LumosQuad to kim@cs.unc.edu. +// +/////////////////////////////////////////////////////////////////////////////////// +// +// This program uses OpenEXR, which has the following restrictions: +// +// Copyright (c) 2002, Industrial Light & Magic, a division of Lucas +// Digital Ltd. LLC +// +// All rights reserved. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are +// met: +// * Redistributions of source code must retain the above copyright +// notice, this list of conditions and the following disclaimer. +// * Redistributions in binary form must reproduce the above +// copyright notice, this list of conditions and the following disclaimer +// in the documentation and/or other materials provided with the +// distribution. +// * Neither the name of Industrial Light & Magic nor the names of +// its contributors may be used to endorse or promote products derived +// from this software without specific prior written permission. +// + +#include "APSF.h" +#include + +////////////////////////////////////////////////////////////////////// +// Construction/Destruction +////////////////////////////////////////////////////////////////////// + +APSF::APSF(int res) : + _res(res) +{ + // make sure kernel is odd dimensions + if (!(_res % 2)) _res++; + _kernel = new float[_res * _res]; + + _q = 0.999; + _R = 400.0f; + _D = 2000.0f; + _T = 1.001f; + _sigma = 0.5f; + + _maxTerms = 600; + _I0 = 1.0f; + _retinaSize = 0.01f; + _eyeSize = 0.025f; +} + +APSF::~APSF() +{ + delete[] _kernel; +} + +////////////////////////////////////////////////////////////////////// +// Legendre polymonial +////////////////////////////////////////////////////////////////////// +float APSF::legendreM(int m, float mu) +{ + vector memoized; + memoized.push_back(1.0f); + memoized.push_back(mu); + + for (int x = 2; x <= m; x++) + { + float newMemo = ((2.0f * (float)x - 1.0f) * mu * memoized[x - 1] - + ((float)x - 1.0f) * memoized[x - 2]) / (float)x; + memoized.push_back(newMemo); + } + + return memoized[m]; +} + +////////////////////////////////////////////////////////////////////// +// scattering function at a point +////////////////////////////////////////////////////////////////////// +float APSF::pointAPSF(float mu) +{ + float total = 0.0f; + + for (int m = 0; m < _maxTerms; m++) + total += (gM(_I0, m) + gM(_I0, m + 1)) * legendreM(m, mu); + + return total; +} + +////////////////////////////////////////////////////////////////////// +// generate a convolution kernel +////////////////////////////////////////////////////////////////////// +void APSF::generateKernelFast() +{ + float dx = _retinaSize / (float)_res; + float dy = _retinaSize / (float)_res; + int halfRes = _res / 2; + float* oneD = new float[_res]; + + float max = 0.0f; + float min = 1000.0f; + int x,y = halfRes; + for (x = 0; x < _res; x++) + { + // calc angle + float diffX = (x - halfRes) * dx; + float diffY = (y - halfRes) * dy; + float distance = sqrt(diffX * diffX + diffY * diffY); + if ((distance / _eyeSize) > (_R / _D)) + oneD[x] = 0.0f; + else + { + float i = -distance * distance * _D * _D + _eyeSize * _eyeSize * _R * _R + distance * distance * _R * _R; + i = _eyeSize * _eyeSize * _D - _eyeSize * sqrt(i); + i /= _eyeSize * _eyeSize + distance * distance; + float mu = M_PI - atan(_retinaSize / distance) - asin((_D - i) / _R); + oneD[x] = pointAPSF(cos(mu)); + + min = (oneD[x] < min) ? oneD[x] : min; + } + max = (oneD[x] > max) ? oneD[x] : max; + } + + // floor + if (min > 0.0f) + { + for (int i = 0; i < _res; i++) + if (oneD[i] > 0.0f) + oneD[i] -= min; + max -= min; + } + + // normalize + if (max > 1.0f) + { + float maxInv = 1.0f / max; + for (int i = 0; i < _res; i++) + oneD[i] *= maxInv; + } + + // interpolate the kernel + int index = 0; + for (y = 0; y < _res; y++) + for (x = 0; x < _res; x++, index++) + { + float dx = fabs((float)(x - halfRes)); + float dy = fabs((float)(y - halfRes)); + float magnitude = sqrtf(dx * dx + dy * dy); + + int lower = floor(magnitude); + if (lower > halfRes - 1) + { + _kernel[index] = 0.0f; + continue; + } + float lerp = magnitude - lower; + _kernel[index] = (1.0f - lerp) * oneD[halfRes + lower] + + lerp * oneD[halfRes + lower + 1]; + + } + + delete[] oneD; +} + +////////////////////////////////////////////////////////////////////// +// save the kernel in binary +////////////////////////////////////////////////////////////////////// +void APSF::write(const char* filename) +{ + // open file + FILE* file; + file = fopen(filename, "wb"); + + fwrite((void*)&_res, sizeof(int), 1, file); + fwrite((void*)&_q, sizeof(float), 1, file); + fwrite((void*)&_T, sizeof(float), 1, file); + fwrite((void*)&_I0, sizeof(float), 1, file); + fwrite((void*)&_sigma, sizeof(float), 1, file); + fwrite((void*)&_R, sizeof(float), 1, file); + fwrite((void*)&_D, sizeof(float), 1, file); + fwrite((void*)&_retinaSize, sizeof(float), 1, file); + fwrite((void*)&_eyeSize, sizeof(float), 1, file); + fwrite((void*)&_maxTerms, sizeof(int), 1, file); + fwrite((void*)&_kernel, sizeof(float) * _res * _res, 1, file); + + fclose(file); +} + +////////////////////////////////////////////////////////////////////// +// load a binary kernel +////////////////////////////////////////////////////////////////////// +void APSF::read(const char* filename) +{ + // open file + FILE* file; + file = fopen(filename, "rb"); + + if (_kernel) delete[] _kernel; + + fread((void*)&_res, sizeof(int), 1, file); + fread((void*)&_q, sizeof(float), 1, file); + fread((void*)&_T, sizeof(float), 1, file); + fread((void*)&_I0, sizeof(float), 1, file); + fread((void*)&_sigma, sizeof(float), 1, file); + fread((void*)&_R, sizeof(float), 1, file); + fread((void*)&_D, sizeof(float), 1, file); + fread((void*)&_retinaSize, sizeof(float), 1, file); + fread((void*)&_eyeSize, sizeof(float), 1, file); + fread((void*)&_maxTerms, sizeof(int), 1, file); + + _kernel = new float[_res * _res]; + fread((void*)&_kernel, sizeof(float) * _res * _res, 1, file); + + fclose(file); +} + +////////////////////////////////////////////////////////////////////// +// write the kernel to a PPM file +////////////////////////////////////////////////////////////////////// +void APSF::writePPM(const char* filename) +{ + unsigned char* ppm = new unsigned char[3 * _res * _res]; + + for (int x = 0; x < _res * _res; x++) + { + ppm[3 * x] = 255 * _kernel[x]; + ppm[3 * x + 1] = 255 * _kernel[x]; + ppm[3 * x + 2] = 255 * _kernel[x]; + } + WritePPM(filename, ppm, _res, _res); + + delete[] ppm; +} diff --git a/APSF.h b/APSF.h new file mode 100644 index 0000000..466dd62 --- /dev/null +++ b/APSF.h @@ -0,0 +1,155 @@ +/////////////////////////////////////////////////////////////////////////////////// +// File : APSF.h +/////////////////////////////////////////////////////////////////////////////////// +// +// LumosQuad - A Lightning Generator +// Copyright 2007 +// The University of North Carolina at Chapel Hill +// +/////////////////////////////////////////////////////////////////////////////////// +// +// 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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// The University of North Carolina at Chapel Hill makes no representations +// about the suitability of this software for any purpose. It is provided +// "as is" without express or implied warranty. +// +// Permission to use, copy, modify and distribute this software and its +// documentation for educational, research and non-profit purposes, without +// fee, and without a written agreement is hereby granted, provided that the +// above copyright notice and the following three paragraphs appear in all +// copies. +// +// THE UNIVERSITY OF NORTH CAROLINA SPECIFICALLY DISCLAIM ANY WARRANTIES, +// INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND +// FITNESS FOR A PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN +// "AS IS" BASIS, AND THE UNIVERSITY OF NORTH CAROLINA HAS NO OBLIGATION TO +// PROVIDE MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS. +// +// Please send questions and comments about LumosQuad to kim@cs.unc.edu. +// +/////////////////////////////////////////////////////////////////////////////////// +// +// This program uses OpenEXR, which has the following restrictions: +// +// Copyright (c) 2002, Industrial Light & Magic, a division of Lucas +// Digital Ltd. LLC +// +// All rights reserved. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are +// met: +// * Redistributions of source code must retain the above copyright +// notice, this list of conditions and the following disclaimer. +// * Redistributions in binary form must reproduce the above +// copyright notice, this list of conditions and the following disclaimer +// in the documentation and/or other materials provided with the +// distribution. +// * Neither the name of Industrial Light & Magic nor the names of +// its contributors may be used to endorse or promote products derived +// from this software without specific prior written permission. +// + +#ifndef APSF_H +#define APSF_H + +#include +#include +#include +#include "ppm/ppm.hpp" + +using namespace std; + +#ifndef M_PI +#define M_PI 3.1415926535897931f +#endif + +//////////////////////////////////////////////////////////////////// +/// \brief Generates the rendering filter +//////////////////////////////////////////////////////////////////// +class APSF +{ +public: + //! constructor + APSF(int res = 512); + //! destructor + virtual ~APSF(); + + //! read in an APSF file + void read(const char* filename); + //! write out an APSF file + void write(const char* filename); + //! write out the current kernel to a PPM file + void writePPM(const char* filename); + + //! generate one line of the kernel and spin it radially + void generateKernelFast(); + + //! resolution of current kernel + int res() { return _res; }; + + //! returns the float array for the kernel + float* kernel() { return _kernel; }; + +private: + //! kernel resolution + int _res; + //! convolution kernel + float* _kernel; + + //////////////////////////////////////////////////////////////////// + // APSF components + //////////////////////////////////////////////////////////////////// + + // scattering parameters + float _q; + float _T; + float _I0; + float _sigma; + float _R; + float _D; + + float _retinaSize; + float _eyeSize; + + //! number of coefficients + int _maxTerms; + + //! function value at a point + float pointAPSF(float mu); + + //////////////////////////////////////////////////////////////////// + // auxiliary functions + //////////////////////////////////////////////////////////////////// + float legendreM(int m, float mu); + float gM(float I0, int m) { + return (m == 0) ? 0.0f : exp(-(betaM(m, _q) * _T + alphaM(m) * log(_T))); + }; + float alphaM(float m) { + return m + 1.0f; + }; + float betaM(float m, float q) { + return ((2.0f * m + 1.0f) / m) * (1.0f - pow(q, (int)m - 1)); + }; + float factorial(float x) { + return (x <= 1.0f) ? 1.0f : x * factorial(x - 1.0f); + }; + float choose(float x, float y) { + return factorial(x) / (factorial(y) * factorial(x - y)); + }; +}; + +#endif diff --git a/BlueNoise/BLUE_NOISE.cpp b/BlueNoise/BLUE_NOISE.cpp new file mode 100644 index 0000000..3adb10e --- /dev/null +++ b/BlueNoise/BLUE_NOISE.cpp @@ -0,0 +1,448 @@ +/////////////////////////////////////////////////////////////////////////////////// +// File : BLUE_NOISE.h +/////////////////////////////////////////////////////////////////////////////////// +// +// LumosQuad - A Lightning Generator +// Copyright 2007 +// The University of North Carolina at Chapel Hill +// +/////////////////////////////////////////////////////////////////////////////////// +// +// 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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// The University of North Carolina at Chapel Hill makes no representations +// about the suitability of this software for any purpose. It is provided +// "as is" without express or implied warranty. +// +// Permission to use, copy, modify and distribute this software and its +// documentation for educational, research and non-profit purposes, without +// fee, and without a written agreement is hereby granted, provided that the +// above copyright notice and the following three paragraphs appear in all +// copies. +// +// THE UNIVERSITY OF NORTH CAROLINA SPECIFICALLY DISCLAIM ANY WARRANTIES, +// INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND +// FITNESS FOR A PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN +// "AS IS" BASIS, AND THE UNIVERSITY OF NORTH CAROLINA HAS NO OBLIGATION TO +// PROVIDE MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS. +// +// Please send questions and comments about LumosQuad to kim@cs.unc.edu. +// +/////////////////////////////////////////////////////////////////////////////////// +// +// This program uses OpenEXR, which has the following restrictions: +// +// Copyright (c) 2002, Industrial Light & Magic, a division of Lucas +// Digital Ltd. LLC +// +// All rights reserved. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are +// met: +// * Redistributions of source code must retain the above copyright +// notice, this list of conditions and the following disclaimer. +// * Redistributions in binary form must reproduce the above +// copyright notice, this list of conditions and the following disclaimer +// in the documentation and/or other materials provided with the +// distribution. +// * Neither the name of Industrial Light & Magic nor the names of +// its contributors may be used to endorse or promote products derived +// from this software without specific prior written permission. +// +/////////////////////////////////////////////////////////////////////////////////// +// +// This class is a very thin wrapper to Daniel Dunbar's blue noise generator. +// With the exception of BLUE_NOISE.h and BLUE_NOISE.cpp, the other files +// in this directory are unmodified copies of his code. +// +// For the original, untainted code, see: +// http://www.cs.virginia.edu/~gfx/pubs/antimony/ +// +/////////////////////////////////////////////////////////////////////////////////// + +// $Id: PDSampling.cpp,v 1.12 2006/07/11 16:45:22 zr Exp $ + +#define _USE_MATH_DEFINES +#include + +#include + +#include "BLUE_NOISE.h" +#include "RangeList.h" +#include "ScallopedSector.h" +#include "WeightedDiscretePDF.h" + +typedef std::vector IntVector; + +/// + +BLUE_NOISE::BLUE_NOISE(float _radius, bool _isTiled, bool usesGrid) : + m_rng(123456), + radius(_radius), + isTiled(_isTiled) +{ + if (usesGrid) { + // grid size is chosen so that 4*radius search only + // requires searching adjacent cells, this also + // determines max points per cell + m_gridSize = (int) ceil(2./(4.*_radius)); + if (m_gridSize<2) m_gridSize = 2; + + m_gridCellSize = 2.0f/m_gridSize; + m_grid = new int[m_gridSize*m_gridSize][kMaxPointsPerCell]; + + for (int y=0; y=a.x && 1>=a.y; +} + +Vec2 BLUE_NOISE::randomPoint() +{ + return Vec2(2*m_rng.getFloatL()-1, 2*m_rng.getFloatL()-1); +} + +Vec2 BLUE_NOISE::getTiled(Vec2 v) +{ + float x = v.x, y = v.y; + + if (isTiled) { + if (x<-1) x += 2; + else if (x>1) x -= 2; + + if (y<-1) y += 2; + else if (y>1) y -= 2; + } + + return Vec2(x,y); +} + +void BLUE_NOISE::getGridXY(Vec2 &v, int *gx_out, int *gy_out) +{ + int gx = *gx_out = (int) floor(.5*(v.x + 1)*m_gridSize); + int gy = *gy_out = (int) floor(.5*(v.y + 1)*m_gridSize); + if (gx<0 || gx>=m_gridSize || gy<0 || gy>=m_gridSize) { + printf("Internal error, point outside grid was generated, ignoring.\n"); + } +} + +void BLUE_NOISE::addPoint(Vec2 pt) +{ + int i, gx, gy, *cell; + + points.push_back(pt); + + if (m_grid) { + getGridXY(pt, &gx, &gy); + cell = m_grid[gy*m_gridSize + gx]; + for (i=0; i(m_gridSize>>1)) N = m_gridSize>>1; + + m_neighbors.clear(); + getGridXY(pt, &gx, &gy); + for (j=-N; j<=N; j++) { + for (i=-N; i<=N; i++) { + int cx = (gx+i+m_gridSize)%m_gridSize; + int cy = (gy+j+m_gridSize)%m_gridSize; + int *cell = m_grid[cy*m_gridSize + cx]; + + for (k=0; k(m_gridSize>>1)) N = m_gridSize>>1; + + getGridXY(pt, &gx, &gy); + for (j=-N; j<=N; j++) { + for (i=-N; i<=N; i++) { + int cx = (gx+i+m_gridSize)%m_gridSize; + int cy = (gy+j+m_gridSize)%m_gridSize; + int *cell = m_grid[cy*m_gridSize + cx]; + + for (k=0; k(m_gridSize>>1)) N = m_gridSize>>1; + + getGridXY(candidate, &gx, &gy); + + int xSide = (candidate.x - (-1 + gx*m_gridCellSize))>m_gridCellSize*.5; + int ySide = (candidate.y - (-1 + gy*m_gridCellSize))>m_gridCellSize*.5; + int iy = 1; + for (j=-N; j<=N; j++) { + int ix = 1; + + if (j==0) iy = ySide; + else if (j==1) iy = 0; + + for (i=-N; i<=N; i++) { + if (i==0) ix = xSide; + else if (i==1) ix = 0; + + // offset to closest cell point + float dx = candidate.x - (-1 + (gx+i+ix)*m_gridCellSize); + float dy = candidate.y - (-1 + (gy+j+iy)*m_gridCellSize); + + if (dx*dx+dy*dy relaxTmpOut.txt"); + + tmp = fopen("relaxTmpOut.txt", "r"); + fscanf(tmp, "%d\n%d\n", &dim, &numVerts); + + if (dim!=2) { + printf("Error calling out to qvoronoi, skipping relaxation.\n"); + goto exit; + } + + verts = new Vec2[numVerts]; + for (int i=0; i RegionMap; + +void BLUE_NOISE::complete() +{ + RangeList rl(0,0); + IntVector candidates; + + addPoint(randomPoint()); + candidates.push_back((int) points.size()-1); + + while (candidates.size()) { + int c = m_rng.getInt32()%candidates.size(); + int index = candidates[c]; + Vec2 candidate = points[index]; + candidates[c] = candidates[candidates.size()-1]; + candidates.pop_back(); + + rl.reset(0, (float) M_PI*2); + findNeighborRanges(index, rl); + while (rl.numRanges) { + RangeEntry &re = rl.ranges[m_rng.getInt32()%rl.numRanges]; + float angle = re.min + (re.max-re.min)*m_rng.getFloatL(); + Vec2 pt = getTiled(Vec2(candidate.x + cos(angle)*2*radius, + candidate.y + sin(angle)*2*radius)); + + addPoint(pt); + candidates.push_back((int) points.size()-1); + + rl.subtract(angle - (float) M_PI/3, angle + (float) M_PI/3); + } + } +} + +void BLUE_NOISE::writeToBool(bool* noise, int size) +{ + // wipe + int index = 0; + for (index = 0; index < size * size; index++) + noise[index] = false; + + for (int x = 0; x < points.size(); x++) + { + int i = (points[x].x + 1.0f) * 0.5f * size; + int j = (points[x].y + 1.0f) * 0.5f * size; + index = i + j * size; + noise[index] = true; + } +} diff --git a/BlueNoise/BLUE_NOISE.h b/BlueNoise/BLUE_NOISE.h new file mode 100644 index 0000000..dbeab3e --- /dev/null +++ b/BlueNoise/BLUE_NOISE.h @@ -0,0 +1,183 @@ +/////////////////////////////////////////////////////////////////////////////////// +// File : BLUE_NOISE.h +/////////////////////////////////////////////////////////////////////////////////// +// +// LumosQuad - A Lightning Generator +// Copyright 2007 +// The University of North Carolina at Chapel Hill +// +/////////////////////////////////////////////////////////////////////////////////// +// +// 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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// The University of North Carolina at Chapel Hill makes no representations +// about the suitability of this software for any purpose. It is provided +// "as is" without express or implied warranty. +// +// Permission to use, copy, modify and distribute this software and its +// documentation for educational, research and non-profit purposes, without +// fee, and without a written agreement is hereby granted, provided that the +// above copyright notice and the following three paragraphs appear in all +// copies. +// +// THE UNIVERSITY OF NORTH CAROLINA SPECIFICALLY DISCLAIM ANY WARRANTIES, +// INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND +// FITNESS FOR A PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN +// "AS IS" BASIS, AND THE UNIVERSITY OF NORTH CAROLINA HAS NO OBLIGATION TO +// PROVIDE MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS. +// +// Please send questions and comments about LumosQuad to kim@cs.unc.edu. +// +/////////////////////////////////////////////////////////////////////////////////// +// +// This program uses OpenEXR, which has the following restrictions: +// +// Copyright (c) 2002, Industrial Light & Magic, a division of Lucas +// Digital Ltd. LLC +// +// All rights reserved. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are +// met: +// * Redistributions of source code must retain the above copyright +// notice, this list of conditions and the following disclaimer. +// * Redistributions in binary form must reproduce the above +// copyright notice, this list of conditions and the following disclaimer +// in the documentation and/or other materials provided with the +// distribution. +// * Neither the name of Industrial Light & Magic nor the names of +// its contributors may be used to endorse or promote products derived +// from this software without specific prior written permission. +// +/////////////////////////////////////////////////////////////////////////////////// +// +// This class is a very thin wrapper to Daniel Dunbar's blue noise generator. +// With the exception of BLUE_NOISE.h and BLUE_NOISE.cpp, the other files +// in this directory are unmodified copies of his code. +// +// For the original, untainted code, see: +// http://www.cs.virginia.edu/~gfx/pubs/antimony/ +// +/////////////////////////////////////////////////////////////////////////////////// + +// $Id: PDSampling.h,v 1.6 2006/07/06 23:13:18 zr Exp $ + +#include "RNG.h" +#include +#include + +#define kMaxPointsPerCell 9 + +class RangeList; +class ScallopedRegion; + +class Vec2 { +public: + Vec2() {}; + Vec2(float _x, float _y) : x(_x), y(_y) {}; + + float x,y; + + float length() { return sqrt(x*x + y*y); } + + bool operator ==(const Vec2 &b) const { return x==b.x && y==b.y; } + Vec2 operator +(Vec2 b) { return Vec2(x+b.x, y+b.y); } + Vec2 operator -(Vec2 b) { return Vec2(x-b.x, y-b.y); } + Vec2 operator *(Vec2 b) { return Vec2(x*b.x, y*b.y); } + Vec2 operator /(Vec2 b) { return Vec2(x/b.x, y*b.y); } + + Vec2 operator +(float n) { return Vec2(x+n, y+n); } + Vec2 operator -(float n) { return Vec2(x-n, y-n); } + Vec2 operator *(float n) { return Vec2(x*n, y*n); } + Vec2 operator /(float n) { return Vec2(x/n, y*n); } + + Vec2 &operator +=(Vec2 b) { x+=b.x; y+=b.y; return *this; } + Vec2 &operator -=(Vec2 b) { x-=b.x; y-=b.y; return *this; } + Vec2 &operator *=(Vec2 b) { x*=b.x; y*=b.y; return *this; } + Vec2 &operator /=(Vec2 b) { x/=b.x; y/=b.y; return *this; } + + Vec2 &operator +=(float n) { x+=n; y+=n; return *this; } + Vec2 &operator -=(float n) { x-=n; y-=n; return *this; } + Vec2 &operator *=(float n) { x*=n; y*=n; return *this; } + Vec2 &operator /=(float n) { x/=n; y/=n; return *this; } +}; + +/// \brief Daniel Dunbar's blue noise generator +/// +/// The original code has been modified so that the 'boundary sampling' +/// method is the only one available. +class BLUE_NOISE { +protected: + RNG m_rng; + std::vector m_neighbors; + + int (*m_grid)[kMaxPointsPerCell]; + int m_gridSize; + float m_gridCellSize; + +public: + std::vector points; + float radius; + bool isTiled; + +public: + BLUE_NOISE(float radius, bool isTiled=true, bool usesGrid=true); + virtual ~BLUE_NOISE() { }; + + // + + bool pointInDomain(Vec2 &a); + + // return shortest distance between _a_ + // and _b_ (accounting for tiling) + float getDistanceSquared(Vec2 &a, Vec2 &b) { Vec2 v = getTiled(b-a); return v.x*v.x + v.y*v.y; } + float getDistance(Vec2 &a, Vec2 &b) { return sqrt(getDistanceSquared(a, b)); } + + // generate a random point in square + Vec2 randomPoint(); + + // return tiled coordinates of _v_ + Vec2 getTiled(Vec2 v); + + // return grid x,y for point + void getGridXY(Vec2 &v, int *gx_out, int *gy_out); + + // add _pt_ to point list and grid + void addPoint(Vec2 pt); + + // populate m_neighbors with list of + // all points within _radius_ of _pt_ + // and return number of such points + int findNeighbors(Vec2 &pt, float radius); + + // return distance to closest neighbor within _radius_ + float findClosestNeighbor(Vec2 &pt, float radius); + + // find available angle ranges on boundary for candidate + // by subtracting occluded neighbor ranges from _rl_ + void findNeighborRanges(int index, RangeList &rl); + + // extend point set by boundary sampling until domain is + // full + void maximize(); + + // apply one step of Lloyd relaxation + void relax(); + + void complete(); + + void writeToBool(bool* noise, int size); +}; diff --git a/BlueNoise/LICENSE.txt b/BlueNoise/LICENSE.txt new file mode 100644 index 0000000..07697aa --- /dev/null +++ b/BlueNoise/LICENSE.txt @@ -0,0 +1,7 @@ +This code is released into the public domain. You can do whatever +you want with it. + +I do ask that you respect the authorship and credit myself (Daniel +Dunbar) when referencing the code. Additionally, if you use the +code in an interesting or integral manner I would like to hear +about it. diff --git a/BlueNoise/README.txt b/BlueNoise/README.txt new file mode 100644 index 0000000..af89ac4 --- /dev/null +++ b/BlueNoise/README.txt @@ -0,0 +1,101 @@ +PDSample - Poisson-Disk sample set generation +Daniel Dunbar, daniel@zuster.org +---- + + +Overview +-- +PDSample generates Poisson-disk sampling sets in the domain [-1,1]^2 +using a variety of methods. See "A Spatial Data Structure for Fast +Poisson-Disk Sample Generation", in Proc' of SIGGRAPH 2006 for more +information. + + +Building +--- +The code should be portable to any platform with 32-bit float's and int's. + +Windows: There is an included PDSample.sln for MSVS version 7. +Unix: Type 'make' and hope for the best. + + +Usage +--- +PDSample [-m] [-t] [-r ] [-M ] + [-N ] + +Options +-- + o -t + Uses tiled (toroidal) domain for supporting samplers. The resulting point + set will be suitable for tiling repeatedly in the x and y directions. + + o -m + Maximize the resulting point set. For samplers which do not already produce + a maximal point set then this will use the Boundary sampling method to + ensure the resulting point set is maximal. + + o -r + Apply the specified number of relaxations to the resulting point set. This + requires that qvoronoi be in the path. + + o -M + For DartThrowing and BestCandidate methods this determines the factor to + multiply the current number of points by to determine how many samples to + take before exiting (DartThrowing) or accepting the best candidate + (BestCandidate). + + o -N + This specifies a minimum number of samples that will be taken for the + DartThrowing sampler. See below. + + +Available Samplers (for method argument) +-- + o DartThrowing + Standard dart throwing. On each iteration the DartThrowing sampler will + try min(N*multiplier,minMaxThrows) samples before termination. Note that + for regular dart throwing where simply a maximum number of throws is used + to determine the termination point, the multiplier should be set to 0. + + o BestCandidate + Mitchell's Best Candidate algorithm. Uses the multiplier argument. + + o Boundary + Dart throwing by maximizing boundaries. + + o Pure + Dart throwing using scalloped sectors. + o LinearPure + Dart throwing using scalloped sectors but without sampling regions according + to their probability of being hit. + + o Penrose + Ostromoukhov et al.'s sampling method using their quasisampler_prototype.h + + o Uniform + Random point generation. The number of samples to take is calculated as + .75/radius^2 to approximately match the density of Poisson-disk sampling. + + +Output format +-- +Point sets are output in a trivial binary format. The format is not intended +for distribution and does not encode the endianness of the generating platform. + +The format matches the pseudo-C struct below: +struct { + int N; // number of points + float t; // generation time + float r; // radius used in generation + int isTiled; // flag for if the set is tileable + float points[N][2]; +}; + + +Acknowledgments +-- +Thanks to Ares Lagae for comments on a preliminary release of the code, +Ostromoukhov et al. for making available their quasisampler implementation, +as well as Takuji Nishimura and Makoto Matsumoto for their Mersenne +Twister random number generator. diff --git a/BlueNoise/RNG.cpp b/BlueNoise/RNG.cpp new file mode 100644 index 0000000..b7e3f8a --- /dev/null +++ b/BlueNoise/RNG.cpp @@ -0,0 +1,174 @@ +/* + A C-program for MT19937, with initialization improved 2002/1/26. + Coded by Takuji Nishimura and Makoto Matsumoto. + Modified to be a C++ class by Daniel Dunbar. + + Before using, initialize the state by using init_genrand(seed) + or init_by_array(init_key, key_length). + + Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura, + All rights reserved. + + Redistribution and use in source and binary forms, with or without + modification, are permitted provided that the following conditions + are met: + + 1. Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + 2. Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in the + documentation and/or other materials provided with the distribution. + + 3. The names of its contributors may not be used to endorse or promote + products derived from this software without specific prior written + permission. + + THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS + "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT + LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR + A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR + CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, + EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, + PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR + PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF + LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING + NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS + SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + + + Any feedback is very welcome. + http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html + email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space) +*/ + +#include "RNG.h" + +/* initializes mt[N] with a seed */ +RNG::RNG(unsigned long s) +{ + seed(s); +} + +/* initialize by an array with array-length */ +/* init_key is the array for initializing keys */ +/* key_length is its length */ +/* slight change for C++, 2004/2/26 */ +RNG::RNG(unsigned long init_key[], int key_length) +{ + int i, j, k; + seed(19650218UL); + i=1; j=0; + k = (N>key_length ? N : key_length); + for (; k; k--) { + mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1664525UL)) + + init_key[j] + j; /* non linear */ + mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */ + i++; j++; + if (i>=N) { mt[0] = mt[N-1]; i=1; } + if (j>=key_length) j=0; + } + for (k=N-1; k; k--) { + mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1566083941UL)) + - i; /* non linear */ + mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */ + i++; + if (i>=N) { mt[0] = mt[N-1]; i=1; } + } + + mt[0] = 0x80000000UL; /* MSB is 1; assuring non-zero initial array */ +} + +void RNG::seed(unsigned long s) +{ + mt[0]= s & 0xffffffffUL; + for (mti=1; mti> 30)) + mti); + /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */ + /* In the previous versions, MSBs of the seed affect */ + /* only MSBs of the array mt[]. */ + /* 2002/01/09 modified by Makoto Matsumoto */ + mt[mti] &= 0xffffffffUL; + /* for >32 bit machines */ + } +} + +/* generates a random number on [0,0xffffffff]-interval */ +unsigned long RNG::getInt32() +{ + unsigned long y; + static unsigned long mag01[2]={0x0UL, _MATRIX_A}; + /* mag01[x] = x * _MATRIX_A for x=0,1 */ + + if (mti >= N) { /* generate N words at one time */ + int kk; + + for (kk=0;kk> 1) ^ mag01[y & 0x1UL]; + } + for (;kk> 1) ^ mag01[y & 0x1UL]; + } + y = (mt[N-1]&_UPPER_MASK)|(mt[0]&_LOWER_MASK); + mt[N-1] = mt[_M-1] ^ (y >> 1) ^ mag01[y & 0x1UL]; + + mti = 0; + } + + y = mt[mti++]; + + /* Tempering */ + y ^= (y >> 11); + y ^= (y << 7) & 0x9d2c5680UL; + y ^= (y << 15) & 0xefc60000UL; + y ^= (y >> 18); + + return y; +} + +/* generates a random number on [0,0x7fffffff]-interval */ +long RNG::getInt31() +{ + return (long)(getInt32()>>1); +} + +/* generates a random number on [0,1]-real-interval */ +double RNG::getDoubleLR() +{ + return getInt32()*(1.0/4294967295.0); + /* divided by 2^32-1 */ +} + +/* generates a random number on [0,1)-real-interval */ +double RNG::getDoubleL() +{ + return getInt32()*(1.0/4294967296.0); + /* divided by 2^32 */ +} + +/* generates a random number on (0,1)-real-interval */ +double RNG::getDouble() +{ + return (((double)getInt32()) + 0.5)*(1.0/4294967296.0); + /* divided by 2^32 */ +} + +float RNG::getFloatLR() +{ + return getInt32()*(1.0f/4294967295.0f); + /* divided by 2^32-1 */ +} +float RNG::getFloatL() +{ + return getInt32()*(1.0f/4294967296.0f); + /* divided by 2^32 */ +} +float RNG::getFloat() +{ + return (getInt32() + 0.5f)*(1.0f/4294967296.0f); + /* divided by 2^32 */ +} + diff --git a/BlueNoise/RNG.h b/BlueNoise/RNG.h new file mode 100644 index 0000000..cd1ec2b --- /dev/null +++ b/BlueNoise/RNG.h @@ -0,0 +1,37 @@ +#include + +using namespace std; + +class RNG { +private: + /* Period parameters */ + static const long N = 624; + static const long _M = 397; + static const unsigned long _MATRIX_A = 0x9908b0dfUL; /* constant vector a */ + static const unsigned long _UPPER_MASK = 0x80000000UL; /* most significant w-r bits */ + static const unsigned long _LOWER_MASK = 0x7fffffffUL; /* least significant r bits */ + +private: + unsigned long mt[N]; /* the array for the state vector */ + int mti; + +public: + RNG(unsigned long seed=5489UL); + RNG(unsigned long *init_key, int key_length); + + void seed(unsigned long seed); + + /* generates a random number on [0,0xffffffff]-interval */ + unsigned long getInt32(); + /* generates a random number on [0,0x7fffffff]-interval */ + long getInt31(); + /* generates a random number on [0,1]-real-interval */ + double getDoubleLR(); + float getFloatLR(); + /* generates a random number on [0,1)-real-interval */ + double getDoubleL(); + float getFloatL(); + /* generates a random number on (0,1)-real-interval */ + double getDouble(); + float getFloat(); +}; diff --git a/BlueNoise/RangeList.cpp b/BlueNoise/RangeList.cpp new file mode 100644 index 0000000..741ea11 --- /dev/null +++ b/BlueNoise/RangeList.cpp @@ -0,0 +1,139 @@ +// $Id: RangeList.cpp,v 1.4 2006/01/24 03:22:14 zr Exp $ + +#define _USE_MATH_DEFINES +#include +#include +#include + +#include "RangeList.h" + +/// + +static const float kSmallestRange = .000001f; + +RangeList::RangeList(float min, float max) +{ + numRanges = 0; + rangesSize = 8; + ranges = new RangeEntry[rangesSize]; + reset(min, max); +} + +RangeList::~RangeList() +{ + delete[] ranges; +} + +void RangeList::reset(float min, float max) +{ + numRanges = 1; + ranges[0].min = min; + ranges[0].max = max; +} + +void RangeList::deleteRange(int pos) +{ + if (postwoPi) { + subtract(a-twoPi, b-twoPi); + } else if (b<0) { + subtract(a+twoPi, b+twoPi); + } else if (a<0) { + subtract(0, b); + subtract(a+twoPi,twoPi); + } else if (b>twoPi) { + subtract(a, twoPi); + subtract(0, b-twoPi); + } else if (numRanges==0) { + ; + } else { + int pos; + + if (a>1; + if (ranges[mid].minranges[pos+1].min) { + pos++; + } else { + return; + } + } + + while (posranges[pos].min) { + if (ranges[pos].max-b + +typedef struct _RangeEntry { + float min, max; +} RangeEntry; + +class RangeList { +public: + RangeEntry *ranges; + int numRanges, rangesSize; + +public: + RangeList(float min, float max); + ~RangeList(); + + void reset(float min, float max); + + void print(); + + void subtract(float min, float max); + +private: + void deleteRange(int pos); + void insertRange(int pos, float min, float max); +}; diff --git a/BlueNoise/ScallopedSector.cpp b/BlueNoise/ScallopedSector.cpp new file mode 100644 index 0000000..ac9bbc7 --- /dev/null +++ b/BlueNoise/ScallopedSector.cpp @@ -0,0 +1,289 @@ +#define _USE_MATH_DEFINES +#include +#include + +#include + +#include "BLUE_NOISE.h" +#include "ScallopedSector.h" + +static const float kTwoPi = (float) (M_PI*2); + +static float integralOfDistToCircle(float x, float d, float r, float k) +{ + if (r1) y = 1; + + float theta = asin(y); + + return (r*(r*(x + + k*theta) + + k*cos(theta)*d_sin_x) + + d*cos(x)*d_sin_x)*.5f; +} + +ScallopedSector::ScallopedSector(Vec2 &_Pt, float _a1, float _a2, Vec2 &P1, float r1, float sign1, Vec2 &P2, float r2, float sign2) +{ + Vec2 v1 = Vec2(P1.x - _Pt.x, P1.y - _Pt.y); + Vec2 v2 = Vec2(P2.x - _Pt.x, P2.y - _Pt.y); + + P = _Pt; + a1 = _a1; + a2 = _a2; + + arcs[0].P = P1; + arcs[0].r = r1; + arcs[0].sign = sign1; + arcs[0].d = sqrt(v1.x*v1.x + v1.y*v1.y); + arcs[0].rSqrd = arcs[0].r*arcs[0].r; + arcs[0].dSqrd = arcs[0].d*arcs[0].d; + arcs[0].theta = atan2(v1.y,v1.x); + arcs[0].integralAtStart = integralOfDistToCircle(a1 - arcs[0].theta, arcs[0].d, arcs[0].r, arcs[0].sign); + + arcs[1].P = P2; + arcs[1].r = r2; + arcs[1].sign = sign2; + arcs[1].d = sqrt(v2.x*v2.x + v2.y*v2.y); + arcs[1].rSqrd = arcs[1].r*arcs[1].r; + arcs[1].dSqrd = arcs[1].d*arcs[1].d; + arcs[1].theta = atan2(v2.y,v2.x); + arcs[1].integralAtStart = integralOfDistToCircle(a1 - arcs[1].theta, arcs[1].d, arcs[1].r, arcs[1].sign); + + area = calcAreaToAngle(a2); +} + +float ScallopedSector::calcAreaToAngle(float angle) +{ + float underInner = integralOfDistToCircle(angle - arcs[0].theta, arcs[0].d, arcs[0].r, arcs[0].sign) - arcs[0].integralAtStart; + float underOuter = integralOfDistToCircle(angle - arcs[1].theta, arcs[1].d, arcs[1].r, arcs[1].sign) - arcs[1].integralAtStart; + + return underOuter-underInner; +} + +float ScallopedSector::calcAngleForArea(float area, RNG &rng) +{ + float lo = a1, hi = a2, cur = lo + (hi-lo)*rng.getFloat(); + + for (int i=0; i<10; i++) { + if (calcAreaToAngle(cur) *regions) +{ + std::vector angles; + + Vec2 v(C.x - P.x, C.y-P.y); + float d = sqrt(v.x*v.x + v.y*v.y); + + if (rFLT_EPSILON) { + float invD = 1.0f/d; + float x = (d*d - r*r + R*R)*invD*.5f; + float k = R*R - x*x; + + if (k>0) { + float y = sqrt(k); + float vx = v.x*invD; + float vy = v.y*invD; + float vx_x = vx*x, vy_x = vy*x; + float vx_y = vx*y, vy_y = vy*y; + float angle; + + angle = canonizeAngle(atan2(C2.y + vy_x + vx_y - P.y, + C2.x + vx_x - vy_y - P.x)); + if (a1outer) { + regions->push_back(ScallopedSector(P, a1, a2, arcs[0].P, arcs[0].r, arcs[0].sign, arcs[1].P, arcs[1].r, arcs[1].sign)); + } else { + if (innerpush_back(ScallopedSector(P, a1, a2, arcs[0].P, arcs[0].r, arcs[0].sign, C, r, -1)); + } + if (d2push_back(ScallopedSector(P, a1, a2, C, r, 1, arcs[1].P, arcs[1].r, arcs[1].sign)); + } + } + } +} + +/// + +ScallopedRegion::ScallopedRegion(Vec2 &P, float r1, float r2, float _minArea) : + minArea(_minArea) +{ + regions = new std::vector; + regions->push_back(ScallopedSector(P, 0, kTwoPi, P, r1, 1, P, r2, 1)); + area = (*regions)[0].area; +} + +ScallopedRegion::~ScallopedRegion() +{ + delete regions; +} + +void ScallopedRegion::subtractDisk(Vec2 C, float r) +{ + std::vector *newRegions = new std::vector; + + area = 0; + for (unsigned int i=0; isize(); i++) { + ScallopedSector &ss = (*regions)[i]; + std::vector *tmp = new std::vector; + + ss.subtractDisk(C, r, tmp); + + for (unsigned int j=0; jsize(); j++) { + ScallopedSector &nss = (*tmp)[j]; + + if (nss.area>minArea) { + area += nss.area; + + if (newRegions->size()) { + ScallopedSector &last = (*newRegions)[newRegions->size()-1]; + if (last.a2==nss.a1 && (last.arcs[0].P==nss.arcs[0].P && last.arcs[0].r==nss.arcs[0].r && last.arcs[0].sign==nss.arcs[0].sign) && + (last.arcs[1].P==nss.arcs[1].P && last.arcs[1].r==nss.arcs[1].r && last.arcs[1].sign==nss.arcs[1].sign)) { + last.a2 = nss.a2; + last.area = last.calcAreaToAngle(last.a2); + continue; + } + } + + newRegions->push_back(nss); + } + } + + delete tmp; + } + + delete regions; + regions = newRegions; +} + +Vec2 ScallopedRegion::sample(RNG &rng) +{ + if (!regions->size()) { + printf("Fatal error, sampled from empty region."); + exit(1); + return Vec2(0,0); + } else { + float a = area*rng.getFloatL(); + ScallopedSector &ss = (*regions)[0]; + + for (unsigned int i=0; isize(); i++) { + ss = (*regions)[i]; + if (a + +typedef struct { + Vec2 P; + float r, sign, d, theta, integralAtStart; + float rSqrd, dSqrd; +} ArcData; + +class ScallopedSector +{ +public: + Vec2 P; + float a1, a2, area; + + ArcData arcs[2]; + +public: + ScallopedSector(Vec2 &_Pt, float _a1, float _a2, Vec2 &P1, float r1, float sign1, Vec2 &P2, float r2, float sign2); + + float calcAreaToAngle(float angle); + float calcAngleForArea(float area, RNG &rng); + Vec2 sample(RNG &rng); + + float distToCurve(float angle, int index); + + void subtractDisk(Vec2 &C, float r, std::vector *regions); + +private: + float canonizeAngle(float angle); + + void distToCircle(float angle, Vec2 &C, float r, float *d1_out, float *d2_out); +}; + +class ScallopedRegion +{ +public: + std::vector *regions; + float minArea; + float area; + +public: + ScallopedRegion(Vec2 &P, float r1, float r2, float minArea=.00000001); + ~ScallopedRegion(); + + bool isEmpty() { return regions->size()==0; } + void subtractDisk(Vec2 C, float r); + + Vec2 sample(RNG &rng); +}; diff --git a/BlueNoise/WeightedDiscretePDF.cpp b/BlueNoise/WeightedDiscretePDF.cpp new file mode 100644 index 0000000..f9bde2f --- /dev/null +++ b/BlueNoise/WeightedDiscretePDF.cpp @@ -0,0 +1,313 @@ +// $Id: WeightedDiscretePDF.cpp,v 1.4 2006/07/07 05:54:31 zr Exp $ + +#include +#include "WeightedDiscretePDF.h" + +// + +template +WDPDF_Node::WDPDF_Node(T key_, float weight_, WDPDF_Node *parent_) +{ + m_mark = false; + + key = key_; + weight = weight_; + sumWeights = 0; + left = right = 0; + parent = parent_; +} + +template +WDPDF_Node::~WDPDF_Node() +{ + if (left) delete left; + if (right) delete right; +} + +// + +template +WeightedDiscretePDF::WeightedDiscretePDF() +{ + m_root = 0; +} + +template +WeightedDiscretePDF::~WeightedDiscretePDF() +{ + if (m_root) delete m_root; +} + +template +void WeightedDiscretePDF::insert(T item, float weight) +{ + WDPDF_Node *p=0, *n=m_root; + + while (n) { + if (n->leftIsRed() && n->rightIsRed()) + split(n); + + p = n; + if (n->key==item) { + throw std::domain_error("insert: argument(item) already in tree"); + } else { + n = (itemkey)?n->left:n->right; + } + } + + n = new WDPDF_Node(item, weight, p); + + if (!p) { + m_root = n; + } else { + if (itemkey) { + p->left = n; + } else { + p->right = n; + } + + split(n); + } + + propogateSumsUp(n); +} + +template +void WeightedDiscretePDF::remove(T item) +{ + WDPDF_Node **np = lookup(item, 0); + WDPDF_Node *child, *n = *np; + + if (!n) { + throw std::domain_error("remove: argument(item) not in tree"); + } else { + if (n->left) { + WDPDF_Node **leftMaxp = &n->left; + + while ((*leftMaxp)->right) + leftMaxp = &(*leftMaxp)->right; + + n->key = (*leftMaxp)->key; + n->weight = (*leftMaxp)->weight; + + np = leftMaxp; + n = *np; + } + + // node now has at most one child + + child = n->left?n->left:n->right; + *np = child; + + if (child) { + child->parent = n->parent; + + if (n->isBlack()) { + lengthen(child); + } + } + + propogateSumsUp(n->parent); + + n->left = n->right = 0; + delete n; + } +} + +template +void WeightedDiscretePDF::update(T item, float weight) +{ + WDPDF_Node *n = *lookup(item, 0); + + if (!n) { + throw std::domain_error("update: argument(item) not in tree"); + } else { + float delta = weight - n->weight; + n->weight = weight; + + for (; n; n=n->parent) { + n->sumWeights += delta; + } + } +} + +template +T WeightedDiscretePDF::choose(float p) +{ + if (p<0.0 || p>=1.0) { + throw std::domain_error("choose: argument(p) outside valid range"); + } else if (!m_root) { + throw std::logic_error("choose: choose() called on empty tree"); + } else { + float w = m_root->sumWeights * p; + WDPDF_Node *n = m_root; + + while (1) { + if (n->left) { + if (wleft->sumWeights) { + n = n->left; + continue; + } else { + w -= n->left->sumWeights; + } + } + if (wweight || !n->right) { + break; // !n->right condition shouldn't be necessary, just sanity check + } + w -= n->weight; + n = n->right; + } + + return n->key; + } +} + +template +bool WeightedDiscretePDF::inTree(T item) +{ + WDPDF_Node *n = *lookup(item, 0); + + return !!n; +} + +// + +template +WDPDF_Node **WeightedDiscretePDF::lookup(T item, WDPDF_Node **parent_out) +{ + WDPDF_Node *n, *p=0, **np=&m_root; + + while ((n = *np)) { + if (n->key==item) { + break; + } else { + p = n; + if (itemkey) { + np = &n->left; + } else { + np = &n->right; + } + } + } + + if (parent_out) + *parent_out = p; + return np; +} + +template +void WeightedDiscretePDF::split(WDPDF_Node *n) +{ + if (n->left) n->left->markBlack(); + if (n->right) n->right->markBlack(); + + if (n->parent) { + WDPDF_Node *p = n->parent; + + n->markRed(); + + if (p->isRed()) { + p->parent->markRed(); + + // not same direction + if (!( (n==p->left && p==p->parent->left) || + (n==p->right && p==p->parent->right))) { + rotate(n); + p = n; + } + + rotate(p); + p->markBlack(); + } + } +} + +template +void WeightedDiscretePDF::rotate(WDPDF_Node *n) +{ + WDPDF_Node *p=n->parent, *pp=p->parent; + + n->parent = pp; + p->parent = n; + + if (n==p->left) { + p->left = n->right; + n->right = p; + if (p->left) p->left->parent = p; + } else { + p->right = n->left; + n->left = p; + if (p->right) p->right->parent = p; + } + + n->setSum(); + p->setSum(); + + if (!pp) { + m_root = n; + } else { + if (p==pp->left) { + pp->left = n; + } else { + pp->right = n; + } + } +} + +template +void WeightedDiscretePDF::lengthen(WDPDF_Node *n) +{ + if (n->isRed()) { + n->markBlack(); + } else if (n->parent) { + WDPDF_Node *sibling = n->sibling(); + + if (sibling && sibling->isRed()) { + n->parent->markRed(); + sibling->markBlack(); + + rotate(sibling); // node sibling is now old sibling child, must be black + sibling = n->sibling(); + } + + // sibling is black + + if (!sibling) { + lengthen(n->parent); + } else if (sibling->leftIsBlack() && sibling->rightIsBlack()) { + if (n->parent->isBlack()) { + sibling->markRed(); + lengthen(n->parent); + } else { + sibling->markRed(); + n->parent->markBlack(); + } + } else { + if (n==n->parent->left && sibling->rightIsBlack()) { + rotate(sibling->left); // sibling->left must be red + sibling->markRed(); + sibling->parent->markBlack(); + sibling = sibling->parent; + } else if (n==n->parent->right && sibling->leftIsBlack()) { + rotate(sibling->right); // sibling->right must be red + sibling->markRed(); + sibling->parent->markBlack(); + sibling = sibling->parent; + } + + // sibling is black, and sibling's far child is red + + rotate(sibling); + if (n->parent->isRed()) sibling->markRed(); + sibling->left->markBlack(); + sibling->right->markBlack(); + } + } +} + +template +void WeightedDiscretePDF::propogateSumsUp(WDPDF_Node *n) +{ + for (; n; n=n->parent) + n->setSum(); +} diff --git a/BlueNoise/WeightedDiscretePDF.h b/BlueNoise/WeightedDiscretePDF.h new file mode 100644 index 0000000..8705ece --- /dev/null +++ b/BlueNoise/WeightedDiscretePDF.h @@ -0,0 +1,57 @@ +// $Id: WeightedDiscretePDF.h,v 1.4 2006/07/07 05:54:31 zr Exp $ + +template +class WDPDF_Node +{ +private: + bool m_mark; + +public: + WDPDF_Node *parent, *left, *right; + T key; + float weight, sumWeights; + +public: + WDPDF_Node(T key_, float weight_, WDPDF_Node *parent_); + ~WDPDF_Node(); + + WDPDF_Node *sibling() { return this==parent->left?parent->right:parent->left; } + + void markRed() { m_mark = true; } + void markBlack() { m_mark = false; } + bool isRed() { return m_mark; } + bool isBlack() { return !m_mark; } + bool leftIsBlack() { return !left || left->isBlack(); } + bool rightIsBlack() { return !right || right->isBlack(); } + bool leftIsRed() { return !leftIsBlack(); } + bool rightIsRed() { return !rightIsBlack(); } + void setSum() { sumWeights = weight + (left?left->sumWeights:0) + (right?right->sumWeights:0); } +}; + +template +class WeightedDiscretePDF +{ +private: + WDPDF_Node *m_root; + +public: + WeightedDiscretePDF(); + ~WeightedDiscretePDF(); + + void insert(T item, float weight); + void update(T item, float newWeight); + void remove(T item); + bool inTree(T item); + + /* pick a tree element according to its + * weight. p should be in [0,1). + */ + T choose(float p); + +private: + WDPDF_Node **lookup(T item, WDPDF_Node **parent_out); + void split(WDPDF_Node *node); + void rotate(WDPDF_Node *node); + void lengthen(WDPDF_Node *node); + void propogateSumsUp(WDPDF_Node *n); +}; diff --git a/CELL.cpp b/CELL.cpp new file mode 100644 index 0000000..c8067c2 --- /dev/null +++ b/CELL.cpp @@ -0,0 +1,225 @@ +/////////////////////////////////////////////////////////////////////////////////// +// File : CELL.cpp +/////////////////////////////////////////////////////////////////////////////////// +// +// LumosQuad - A Lightning Generator +// Copyright 2007 +// The University of North Carolina at Chapel Hill +// +/////////////////////////////////////////////////////////////////////////////////// +// +// 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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// The University of North Carolina at Chapel Hill makes no representations +// about the suitability of this software for any purpose. It is provided +// "as is" without express or implied warranty. +// +// Permission to use, copy, modify and distribute this software and its +// documentation for educational, research and non-profit purposes, without +// fee, and without a written agreement is hereby granted, provided that the +// above copyright notice and the following three paragraphs appear in all +// copies. +// +// THE UNIVERSITY OF NORTH CAROLINA SPECIFICALLY DISCLAIM ANY WARRANTIES, +// INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND +// FITNESS FOR A PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN +// "AS IS" BASIS, AND THE UNIVERSITY OF NORTH CAROLINA HAS NO OBLIGATION TO +// PROVIDE MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS. +// +// Please send questions and comments about LumosQuad to kim@cs.unc.edu. +// +/////////////////////////////////////////////////////////////////////////////////// +// +// This program uses OpenEXR, which has the following restrictions: +// +// Copyright (c) 2002, Industrial Light & Magic, a division of Lucas +// Digital Ltd. LLC +// +// All rights reserved. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are +// met: +// * Redistributions of source code must retain the above copyright +// notice, this list of conditions and the following disclaimer. +// * Redistributions in binary form must reproduce the above +// copyright notice, this list of conditions and the following disclaimer +// in the documentation and/or other materials provided with the +// distribution. +// * Neither the name of Industrial Light & Magic nor the names of +// its contributors may be used to endorse or promote products derived +// from this software without specific prior written permission. +// + +#include "CELL.h" + +////////////////////////////////////////////////////////////////////// +// Construction/Destruction +////////////////////////////////////////////////////////////////////// + +// normal cell constructor +CELL::CELL(float north, float east, float south, float west, CELL* parent, int depth) : + parent(parent), depth(depth), index(-1), candidate(false), + boundary(false), potential(0.0f), state(EMPTY) +{ + for (int x = 0; x < 4; x++) + children[x] = NULL; + for (int x = 0; x < 8; x++) + neighbors[x] = NULL; + + bounds[0] = north; bounds[1] = east; bounds[2] = south; bounds[3] = west; + + center[0] = (bounds[1] + bounds[3]) * 0.5f; + center[1] = (bounds[0] + bounds[2]) * 0.5f; +} + +// ghost cell constructor +CELL::CELL(int depth) : parent(NULL), depth(depth), index(-1), candidate(false), + boundary(true), potential(0.0f), state(EMPTY) +{ + for (int x = 0; x < 4; x++) + children[x] = NULL; + for (int x = 0; x < 8; x++) + neighbors[x] = NULL; + + bounds[0] = 0.0f; bounds[1] = 0.0f; bounds[2] = 0.0f; bounds[3] = 0.0f; + center[0] = 0.0f; center[1] = 0.0f; +} + +CELL::~CELL() { + int x; + + for (x = 0; x < 4; x++) + if (children[x] != NULL) + { + delete children[x]; + children[x] = NULL; + } +} + +////////////////////////////////////////////////////////////////////// +// refine current cell +////////////////////////////////////////////////////////////////////// +void CELL::refine() { + if (children[0] != NULL) return; + float center[] = {(bounds[0] + bounds[2]) * 0.5f, (bounds[1] + bounds[3]) * 0.5f}; + + children[0] = new CELL(bounds[0], center[1], center[0], bounds[3], this, depth + 1); + children[1] = new CELL(bounds[0], bounds[1], center[0], center[1], this, depth + 1); + children[2] = new CELL(center[0], bounds[1], bounds[2], center[1], this, depth + 1); + children[3] = new CELL(center[0], center[1], bounds[2], bounds[3], this, depth + 1); + + children[0]->potential = potential; + children[1]->potential = potential; + children[2]->potential = potential; + children[3]->potential = potential; +} + +////////////////////////////////////////////////////////////////////// +// return north neighbor to current cell +////////////////////////////////////////////////////////////////////// +CELL* CELL::northNeighbor() +{ + // if it is the root + if (this->parent == NULL) return NULL; + + // if it is the southern child of the parent + if (parent->children[3] == this) return parent->children[0]; + if (parent->children[2] == this) return parent->children[1]; + + // else look up higher + CELL* mu = parent->northNeighbor(); + + // if there are no more children to look at, + // this is the answer + if (mu == NULL || mu->children[0] == NULL) return mu; + // if it is the NW child of the parent + else if (parent->children[0] == this) return mu->children[3]; + // if it is the NE child of the parent + else return mu->children[2]; +} + +////////////////////////////////////////////////////////////////////// +// return north neighbor to current cell +////////////////////////////////////////////////////////////////////// +CELL* CELL::southNeighbor() +{ + // if it is the root + if (this->parent == NULL) return NULL; + + // if it is the northern child of the parent + if (parent->children[0] == this) return parent->children[3]; + if (parent->children[1] == this) return parent->children[2]; + + // else look up higher + CELL* mu = parent->southNeighbor(); + + // if there are no more children to look at, + // this is the answer + if (mu == NULL || mu->children[0] == NULL) return mu; + // if it is the SW child of the parent + else if (parent->children[3] == this) return mu->children[0]; + // if it is the SE child of the parent + else return mu->children[1]; +} + +////////////////////////////////////////////////////////////////////// +// return north neighbor to current cell +////////////////////////////////////////////////////////////////////// +CELL* CELL::westNeighbor() +{ + // if it is the root + if (this->parent == NULL) return NULL; + + // if it is the eastern child of the parent + if (parent->children[1] == this) return parent->children[0]; + if (parent->children[2] == this) return parent->children[3]; + + // else look up higher + CELL* mu = parent->westNeighbor(); + + // if there are no more children to look at, + // this is the answer + if (mu == NULL || mu->children[0] == NULL) return mu; + // if it is the NW child of the parent + else if (parent->children[0] == this) return mu->children[1]; + // if it is the SW child of the parent + else return mu->children[2]; +} + +////////////////////////////////////////////////////////////////////// +// return north neighbor to current cell +////////////////////////////////////////////////////////////////////// +CELL* CELL::eastNeighbor() +{ + // if it is the root + if (this->parent == NULL) return NULL; + + // if it is the western child of the parent + if (parent->children[0] == this) return parent->children[1]; + if (parent->children[3] == this) return parent->children[2]; + + // else look up higher + CELL* mu = parent->eastNeighbor(); + + // if there are no more children to look at, + // this is the answer + if (mu == NULL || mu->children[0] == NULL) return mu; + // if it is the NE child of the parent + else if (parent->children[1] == this) return mu->children[0]; + // if it is the SE child of the parent + else return mu->children[3]; +} + diff --git a/CELL.h b/CELL.h new file mode 100644 index 0000000..27256ed --- /dev/null +++ b/CELL.h @@ -0,0 +1,194 @@ +/////////////////////////////////////////////////////////////////////////////////// +// File : CELL.h +/////////////////////////////////////////////////////////////////////////////////// +// +// LumosQuad - A Lightning Generator +// Copyright 2007 +// The University of North Carolina at Chapel Hill +// +/////////////////////////////////////////////////////////////////////////////////// +// +// 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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// The University of North Carolina at Chapel Hill makes no representations +// about the suitability of this software for any purpose. It is provided +// "as is" without express or implied warranty. +// +// Permission to use, copy, modify and distribute this software and its +// documentation for educational, research and non-profit purposes, without +// fee, and without a written agreement is hereby granted, provided that the +// above copyright notice and the following three paragraphs appear in all +// copies. +// +// THE UNIVERSITY OF NORTH CAROLINA SPECIFICALLY DISCLAIM ANY WARRANTIES, +// INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND +// FITNESS FOR A PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN +// "AS IS" BASIS, AND THE UNIVERSITY OF NORTH CAROLINA HAS NO OBLIGATION TO +// PROVIDE MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS. +// +// Please send questions and comments about LumosQuad to kim@cs.unc.edu. +// +/////////////////////////////////////////////////////////////////////////////////// +// +// This program uses OpenEXR, which has the following restrictions: +// +// Copyright (c) 2002, Industrial Light & Magic, a division of Lucas +// Digital Ltd. LLC +// +// All rights reserved. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are +// met: +// * Redistributions of source code must retain the above copyright +// notice, this list of conditions and the following disclaimer. +// * Redistributions in binary form must reproduce the above +// copyright notice, this list of conditions and the following disclaimer +// in the documentation and/or other materials provided with the +// distribution. +// * Neither the name of Industrial Light & Magic nor the names of +// its contributors may be used to endorse or promote products derived +// from this software without specific prior written permission. +// + +#ifndef CELL_H +#define CELL_H +#include + +////////////////////////////////////////////////////////////////////// +/// \enum Possible states of the cell in the DBM simulation +////////////////////////////////////////////////////////////////////// +enum CELL_STATE {EMPTY, NEGATIVE, POSITIVE, REPULSOR, ATTRACTOR}; + +////////////////////////////////////////////////////////////////////// +/// \brief Basic cell data structure of the quadtree +////////////////////////////////////////////////////////////////////// +class CELL +{ +public: + //! normal cell constructor + CELL(float north, + float east, + float south, + float west, + CELL* parent = NULL, + int depth = 0); + + //! ghost cell constructor + CELL(int depth = 0); + + //! destructor + ~CELL(); + + //! The children of the node in the quadtree + /*! + Winding order of children is: + + \verbatim + _________ + | | | + | 0 | 1 | + |___|___| + | | | + | 3 | 2 | + |___|___| + \endverbatim */ + CELL* children[4]; + + //! The physical bounds of the current grid cell + /*! + Winding order of bounds is: + + \verbatim + 0 - north + 1 - east + 2 - south + 3 - west + \endverbatim */ + float bounds[4]; + + //! The neighbors in the balanced quadtree + /*! + winding order of the neighbors is: + + \verbatim + | 0 | 1 | + ____|____|____|_____ + | | + 7 | | 2 + ____| |_____ + | | + 6 | | 3 + ____|_________|_____ + | | | + | 5 | 4 | + \endverbatim + + Neighbors 0,2,4,6 should always exist. Depending on + if the neighbor is on a lower refinement level, + neighbors 1,3,5,7 may or may not exist. If they are not + present, the pointer value should ne NULL. */ + CELL* neighbors[8]; + + //! Poisson stencil coefficients + /*! + winding order of the stencil coefficients: + + \verbatim + | 0 | 1 | + ____|____|____|_____ + | | + 7 | | 2 + ____| 8 |_____ + | | + 6 | | 3 + ____|_________|_____ + | | | + | 5 | 4 | + \endverbatim + Stencils 0,2,4,6 should always exist. Depending on + if the neighbor is on a lower refinement level, + stencils 1,3,5,7 may or may not exist. If they are not + present, the pointer value should ne NULL. */ + float stencil[9]; + + float center[2]; ///< center of the cell + int depth; ///< current tree depth + bool candidate; ///< already a member of candidate list? + + CELL* parent; ///< parent node in the quadtree + CELL_STATE state; ///< DBM state of the cell + + void refine(); ///< subdivide the cell + + //////////////////////////////////////////////////////////////// + // solver-related variables + //////////////////////////////////////////////////////////////// + bool boundary; ///< boundary node to include in the solver? + float potential; ///< current electric potential + float b; ///< rhs of the linear system + float residual; ///< residual in the linear solver + int index; ///< lexicographic index for the solver + + //////////////////////////////////////////////////////////////// + // neighbor lookups + //////////////////////////////////////////////////////////////// + CELL* northNeighbor(); ///< lookup northern neighbor + CELL* southNeighbor(); ///< lookup southern neighbor + CELL* westNeighbor(); ///< lookup western neighbor + CELL* eastNeighbor(); ///< lookup eastern neighbor +}; + +#endif diff --git a/CG_SOLVER.cpp b/CG_SOLVER.cpp new file mode 100644 index 0000000..6260c34 --- /dev/null +++ b/CG_SOLVER.cpp @@ -0,0 +1,309 @@ +/////////////////////////////////////////////////////////////////////////////////// +// File : CG_SOLVER.cpp +/////////////////////////////////////////////////////////////////////////////////// +// +// LumosQuad - A Lightning Generator +// Copyright 2007 +// The University of North Carolina at Chapel Hill +// +/////////////////////////////////////////////////////////////////////////////////// +// +// 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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// The University of North Carolina at Chapel Hill makes no representations +// about the suitability of this software for any purpose. It is provided +// "as is" without express or implied warranty. +// +// Permission to use, copy, modify and distribute this software and its +// documentation for educational, research and non-profit purposes, without +// fee, and without a written agreement is hereby granted, provided that the +// above copyright notice and the following three paragraphs appear in all +// copies. +// +// THE UNIVERSITY OF NORTH CAROLINA SPECIFICALLY DISCLAIM ANY WARRANTIES, +// INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND +// FITNESS FOR A PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN +// "AS IS" BASIS, AND THE UNIVERSITY OF NORTH CAROLINA HAS NO OBLIGATION TO +// PROVIDE MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS. +// +// Please send questions and comments about LumosQuad to kim@cs.unc.edu. +// +/////////////////////////////////////////////////////////////////////////////////// +// +// This program uses OpenEXR, which has the following restrictions: +// +// Copyright (c) 2002, Industrial Light & Magic, a division of Lucas +// Digital Ltd. LLC +// +// All rights reserved. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are +// met: +// * Redistributions of source code must retain the above copyright +// notice, this list of conditions and the following disclaimer. +// * Redistributions in binary form must reproduce the above +// copyright notice, this list of conditions and the following disclaimer +// in the documentation and/or other materials provided with the +// distribution. +// * Neither the name of Industrial Light & Magic nor the names of +// its contributors may be used to endorse or promote products derived +// from this software without specific prior written permission. +// + +#include "CG_SOLVER.h" + +////////////////////////////////////////////////////////////////////// +// Construction/Destruction +////////////////////////////////////////////////////////////////////// + +CG_SOLVER::CG_SOLVER(int maxDepth, int iterations, int digits) : + _iterations(iterations), + _arraySize(0), _listSize(0), _digits(digits), + _direction(NULL), _residual(NULL), _q(NULL), _potential(NULL) +{ + // compute the physical size of various grid cells + _dx = new float[maxDepth + 1]; + _dx[0] = 1.0f; + for (int x = 1; x <= maxDepth; x++) + _dx[x] = _dx[x - 1] * 0.5f; +} + +CG_SOLVER::~CG_SOLVER() +{ + if (_direction) delete[] _direction; + if (_residual) delete[] _residual; + if (_potential) delete[] _potential; + if (_q) delete[] _q; +} + +////////////////////////////////////////////////////////////////////// +// reallocate scratch arrays if necessary +////////////////////////////////////////////////////////////////////// +void CG_SOLVER::reallocate() +{ + // if we have enough size already, return + if (_arraySize >= _listSize) return; + + // made sure it SSE aligns okay + _arraySize = _listSize * 2; + if (_arraySize % 4) + _arraySize += 4 - _arraySize % 4; + + // delete the old ones + if (_direction) delete[] _direction; + if (_residual) delete[] _residual; + if (_q) delete[] _q; + + // allocate the new ones + _direction = new float[_arraySize]; + _residual = new float[_arraySize]; + _q = new float[_arraySize]; + + // wipe the new ones + for (int x = 0; x < _arraySize; x++) + _direction[x] = _residual[x] = _q[x] = 0.0f; +} + +////////////////////////////////////////////////////////////////////// +// conjugate gradient solver +////////////////////////////////////////////////////////////////////// +int CG_SOLVER::solve(list cells) +{ + // counters + int x, y, index; + list::iterator cellIterator; + + // i = 0 + int i = 0; + + // precalculate stencils + calcStencils(cells); + + // reallocate scratch arrays if necessary + _listSize = cells.size(); + reallocate(); + + // compute a new lexicographical order + cellIterator = cells.begin(); + for (x = 0; x < _listSize; x++, cellIterator++) + (*cellIterator)->index = x; + + // r = b - Ax + calcResidual(cells); + + // copy residual into easy array + // d = r + cellIterator = cells.begin(); + float deltaNew = 0.0f; + for (x = 0; x < _listSize; x++, cellIterator++) + { + _direction[x] = _residual[x]; + deltaNew += _residual[x] * _residual[x]; + } + + // delta0 = deltaNew + float delta0 = deltaNew; + + // While deltaNew > (eps^2) * delta0 + float eps = pow(10.0f, (float)-_digits); + float maxR = 2.0f * eps; + while ((i < _iterations) && (maxR > eps)) + { + // q = Ad + cellIterator = cells.begin(); + for (y = 0; y < _listSize; y++, cellIterator++) + { + CELL* currentCell = *cellIterator; + CELL** neighbors = currentCell->neighbors; + + float neighborSum = 0.0f; + for (int x = 0; x < 4; x++) + { + int j = x * 2; + neighborSum += _direction[neighbors[j]->index] * currentCell->stencil[j]; + if (neighbors[j+1]) + neighborSum += _direction[neighbors[j+1]->index] * currentCell->stencil[j+1]; + } + _q[y] = -neighborSum + _direction[y] * currentCell->stencil[8]; + } + // alpha = deltaNew / (transpose(d) * q) + float alpha = 0.0f; + for (x = 0; x < _listSize; x++) + alpha += _direction[x] * _q[x]; + if (fabs(alpha) > 0.0f) + alpha = deltaNew / alpha; + + // x = x + alpha * d + cellIterator = cells.begin(); + for (x = 0; x < _listSize; x++, cellIterator++) + (*cellIterator)->potential += alpha * _direction[x]; + + // r = r - alpha * q + maxR = 0.0f; + for (x = 0; x < _listSize; x++) + { + _residual[x] -= _q[x] * alpha; + maxR = (_residual[x] > maxR) ? _residual[x] : maxR; + } + + // deltaOld = deltaNew + float deltaOld = deltaNew; + + // deltaNew = transpose(r) * r + deltaNew = 0.0f; + for (x = 0; x < _listSize; x++) + deltaNew += _residual[x] * _residual[x]; + + // beta = deltaNew / deltaOld + float beta = deltaNew / deltaOld; + + // d = r + beta * d + for (x = 0; x < _listSize; x++) + _direction[x] = _residual[x] + beta * _direction[x]; + + // i = i + 1 + i++; + } + + return i; +} + +////////////////////////////////////////////////////////////////////// +// calculate the residuals +////////////////////////////////////////////////////////////////////// +float CG_SOLVER::calcResidual(list cells) +{ + float maxResidual = 0.0f; + + list::iterator cellIterator = cells.begin(); + for (int i = 0; i < _listSize; i++, cellIterator++) + { + CELL* currentCell = *cellIterator; + float dx = _dx[currentCell->depth]; + float neighborSum = 0.0f; + + for (int x = 0; x < 4; x++) + { + int i = x * 2; + neighborSum += currentCell->neighbors[i]->potential * currentCell->stencil[i]; + if (currentCell->neighbors[i+1]) + neighborSum += currentCell->neighbors[i+1]->potential * currentCell->stencil[i+1]; + } + _residual[i] = currentCell->b - (-neighborSum + currentCell->potential * currentCell->stencil[8]); + + if (fabs(_residual[i]) > maxResidual) + maxResidual = fabs(_residual[i]); + } + return maxResidual; +} + +////////////////////////////////////////////////////////////////////// +// compute stencils once and store +////////////////////////////////////////////////////////////////////// +void CG_SOLVER::calcStencils(list cells) +{ + list::iterator cellIterator = cells.begin(); + for (cellIterator = cells.begin(); cellIterator != cells.end(); cellIterator++) + { + CELL* currentCell = *cellIterator; + float invDx = 1.0f / _dx[currentCell->depth]; + + // sum over faces + float deltaSum = 0.0f; + float bSum = 0.0f; + + for (int x = 0; x < 4; x++) + { + int i = x * 2; + currentCell->stencil[i] = 0.0f; + currentCell->stencil[i+1] = 0.0f; + + if (currentCell->neighbors[i + 1] == NULL) { + // if it is the same refinement level (case 1) + if (currentCell->depth == currentCell->neighbors[i]->depth) { + deltaSum += invDx; + if (!currentCell->neighbors[i]->boundary) + currentCell->stencil[i] = invDx; + else + bSum += (currentCell->neighbors[i]->potential) * invDx; + } + // else it is less refined (case 3) + else { + deltaSum += 0.5f * invDx; + if (!currentCell->neighbors[i]->boundary) + currentCell->stencil[i] = 0.5f * invDx; + else + bSum += currentCell->neighbors[i]->potential * 0.5f * invDx; + } + } + // if the neighbor is at a lower level (case 2) + else { + deltaSum += 2.0f * invDx; + if (!currentCell->neighbors[i]->boundary) + currentCell->stencil[i] = invDx; + else + bSum += currentCell->neighbors[i]->potential * invDx; + if (!currentCell->neighbors[i+1]->boundary) + currentCell->stencil[i+1] = invDx; + else + bSum += currentCell->neighbors[i+1]->potential * invDx; + } + } + + currentCell->stencil[8] = deltaSum; + currentCell->b = bSum; + } +} diff --git a/CG_SOLVER.h b/CG_SOLVER.h new file mode 100644 index 0000000..add6e2f --- /dev/null +++ b/CG_SOLVER.h @@ -0,0 +1,120 @@ +/////////////////////////////////////////////////////////////////////////////////// +// File : CG_SOLVER.h +/////////////////////////////////////////////////////////////////////////////////// +// +// LumosQuad - A Lightning Generator +// Copyright 2007 +// The University of North Carolina at Chapel Hill +// +/////////////////////////////////////////////////////////////////////////////////// +// +// 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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// The University of North Carolina at Chapel Hill makes no representations +// about the suitability of this software for any purpose. It is provided +// "as is" without express or implied warranty. +// +// Permission to use, copy, modify and distribute this software and its +// documentation for educational, research and non-profit purposes, without +// fee, and without a written agreement is hereby granted, provided that the +// above copyright notice and the following three paragraphs appear in all +// copies. +// +// THE UNIVERSITY OF NORTH CAROLINA SPECIFICALLY DISCLAIM ANY WARRANTIES, +// INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND +// FITNESS FOR A PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN +// "AS IS" BASIS, AND THE UNIVERSITY OF NORTH CAROLINA HAS NO OBLIGATION TO +// PROVIDE MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS. +// +// Please send questions and comments about LumosQuad to kim@cs.unc.edu. +// +/////////////////////////////////////////////////////////////////////////////////// +// +// This program uses OpenEXR, which has the following restrictions: +// +// Copyright (c) 2002, Industrial Light & Magic, a division of Lucas +// Digital Ltd. LLC +// +// All rights reserved. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are +// met: +// * Redistributions of source code must retain the above copyright +// notice, this list of conditions and the following disclaimer. +// * Redistributions in binary form must reproduce the above +// copyright notice, this list of conditions and the following disclaimer +// in the documentation and/or other materials provided with the +// distribution. +// * Neither the name of Industrial Light & Magic nor the names of +// its contributors may be used to endorse or promote products derived +// from this software without specific prior written permission. +// + +#ifndef CG_SOLVER_H +#define CG_SOLVER_H + +#include "CELL.h" +#include +#include + +using namespace std; + +//////////////////////////////////////////////////////////////////// +/// \brief Conjugate gradient Poisson solver. +//////////////////////////////////////////////////////////////////// +class CG_SOLVER +{ +public: + //! constructor + CG_SOLVER(int maxDepth, int iterations = 10, int digits = 8); + //! destructor + virtual ~CG_SOLVER(); + + //! solve the Poisson problem + virtual int solve(list cells); + + //! calculate the residual + float calcResidual(list cells); + + //! accessor for the maximum number of iterations + int& iterations() { return _iterations; }; + +protected: + int _iterations; ///< maximum number of iterations + int _digits; ///< desired digits of precision + + //////////////////////////////////////////////////////////////// + // conjugate gradient arrays + //////////////////////////////////////////////////////////////// + float* _direction; ///< conjugate gradient 'd' array + float* _potential; ///< conjugate gradient solution, 'x' array + float* _residual; ///< conjugate gradient residual, 'r' array + float* _q; ///< conjugate gradient 'q' array + + int _arraySize; ///< currently allocated array size + int _listSize; ///< current system size + + //! compute stencils once and store + void calcStencils(list cells); + + //! reallocate the scratch arrays + virtual void reallocate(); + + //! physical lengths of various cell sizes + float* _dx; +}; + +#endif diff --git a/CG_SOLVER_SSE.cpp b/CG_SOLVER_SSE.cpp new file mode 100644 index 0000000..fa34894 --- /dev/null +++ b/CG_SOLVER_SSE.cpp @@ -0,0 +1,416 @@ +/////////////////////////////////////////////////////////////////////////////////// +// File : CG_SOLVER_SSE.cpp +/////////////////////////////////////////////////////////////////////////////////// +// +// LumosQuad - A Lightning Generator +// Copyright 2007 +// The University of North Carolina at Chapel Hill +// +/////////////////////////////////////////////////////////////////////////////////// +// +// 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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// The University of North Carolina at Chapel Hill makes no representations +// about the suitability of this software for any purpose. It is provided +// "as is" without express or implied warranty. +// +// Permission to use, copy, modify and distribute this software and its +// documentation for educational, research and non-profit purposes, without +// fee, and without a written agreement is hereby granted, provided that the +// above copyright notice and the following three paragraphs appear in all +// copies. +// +// THE UNIVERSITY OF NORTH CAROLINA SPECIFICALLY DISCLAIM ANY WARRANTIES, +// INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND +// FITNESS FOR A PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN +// "AS IS" BASIS, AND THE UNIVERSITY OF NORTH CAROLINA HAS NO OBLIGATION TO +// PROVIDE MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS. +// +// Please send questions and comments about LumosQuad to kim@cs.unc.edu. +// +/////////////////////////////////////////////////////////////////////////////////// +// +// This program uses OpenEXR, which has the following restrictions: +// +// Copyright (c) 2002, Industrial Light & Magic, a division of Lucas +// Digital Ltd. LLC +// +// All rights reserved. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are +// met: +// * Redistributions of source code must retain the above copyright +// notice, this list of conditions and the following disclaimer. +// * Redistributions in binary form must reproduce the above +// copyright notice, this list of conditions and the following disclaimer +// in the documentation and/or other materials provided with the +// distribution. +// * Neither the name of Industrial Light & Magic nor the names of +// its contributors may be used to endorse or promote products derived +// from this software without specific prior written permission. +// + +#include "CG_SOLVER_SSE.h" + +////////////////////////////////////////////////////////////////////// +// Construction/Destruction +////////////////////////////////////////////////////////////////////// + +CG_SOLVER_SSE::CG_SOLVER_SSE(int maxDepth, int iterations, int digits) : + CG_SOLVER(maxDepth, iterations, digits) +{ +} + +CG_SOLVER_SSE::~CG_SOLVER_SSE() +{ + if (_direction) free(_direction); + if (_potential) free(_potential); + if (_residual) free(_residual); + if (_q) free(_q); + + _direction = NULL; + _residual = NULL; + _q = NULL; + _potential = NULL; +} + +////////////////////////////////////////////////////////////////////// +// reallocate the sse arrays +////////////////////////////////////////////////////////////////////// +void CG_SOLVER_SSE::reallocate() +{ + if (_arraySize >= _listSize) return; + _arraySize = _listSize * 2; + + if (_arraySize % 4) + _arraySize += 4 - _arraySize % 4; + + if (_direction) free(_direction); + if (_potential) free(_potential); + if (_residual) free(_residual); + if (_q) free(_q); + + _direction = (float*) aligned_alloc(_arraySize * sizeof(float), 16); + _potential = (float*) aligned_alloc(_arraySize * sizeof(float), 16); + _residual = (float*) aligned_alloc(_arraySize * sizeof(float), 16); + _q = (float*) aligned_alloc(_arraySize * sizeof(float), 16); + + return; +} + +////////////////////////////////////////////////////////////////////// +// solve the linear system +////////////////////////////////////////////////////////////////////// +int CG_SOLVER_SSE::solve(list cells) +{ + // counters + int x, y, index; + list::iterator cellIterator; + + // i = 0 + int i = 0; + + // precalculate stencils + calcStencils(cells); + + // reallocate scratch arrays if necessary + _listSize = cells.size(); + reallocate(); + wipeSSE(_potential); + wipeSSE(_direction); + wipeSSE(_residual); + wipeSSE(_q); + + // compute a new lexicographical order + cellIterator = cells.begin(); + for (x = 0; x < _listSize; x++, cellIterator++) + { + CELL* cell = *cellIterator; + cell->index = x; + _potential[x] = cell->potential; + } + + // r = b - Ax + calcResidual(cells); + + // d = r + copySSE(_direction, _residual); + + // deltaNew = r^T r + float deltaNew = dotSSE(_residual, _residual); + + // delta0 = deltaNew + float delta0 = deltaNew; + + // While deltaNew > (eps^2) * delta0 + float eps = pow(10.0f, (float)-_digits); + float maxR = 2.0f * eps; + while ((i < _iterations) && (maxR > eps)) + { + // q = Ad + cellIterator = cells.begin(); + for (y = 0; y < _listSize; y++, cellIterator++) + { + CELL* currentCell = *cellIterator; + CELL** neighbors = currentCell->neighbors; + float* stencil = currentCell->stencil; + + float neighborSum = 0.0f; + for (int x = 0; x < 8; x++) + { + if (neighbors[x]) + neighborSum += _direction[neighbors[x]->index] * stencil[x]; + } + _q[y] = -neighborSum + _direction[y] * currentCell->stencil[8]; + } + + // alpha = deltaNew / (transpose(d) * q) + float alpha = dotSSE(_q, _direction); + if (fabs(alpha) > 0.0f) + alpha = deltaNew / alpha; + + // x = x + alpha * d + saxpySSE(alpha, _direction, _potential); + + // r = r - alpha * q + saxpySSE(-alpha, _q, _residual); + maxR = maxSSE(_residual); + + // deltaOld = deltaNew + float deltaOld = deltaNew; + + // deltaNew = transpose(r) * r + deltaNew = dotSSE(_residual, _residual); + + // beta = deltaNew / deltaOld + float beta = deltaNew / deltaOld; + + // d = r + beta * d + saypxSSE(beta, _residual, _direction); + + // i = i + 1 + i++; + } + + // copy back into the tree + cellIterator = cells.begin(); + for (x = 0; x < _listSize; x++, cellIterator++) + (*cellIterator)->potential = _potential[x]; + + return i; +} + +////////////////////////////////////////////////////////////////////// +// dot product of two vectors +////////////////////////////////////////////////////////////////////// +float CG_SOLVER_SSE::dotSSE(float* x, float* y) +{ + __m128 sum = _mm_set_ps1(0.0f); + __m128* xSSE = (__m128*)x; + __m128* ySSE = (__m128*)y; + __m128 temp; + for (int index = 0; index < _arraySize / 4; index++) + { + temp = _mm_mul_ps(*xSSE, *ySSE); + sum = _mm_add_ps(sum, temp); + xSSE++; + ySSE++; + } + union u { + __m128 m; + float f[4]; + } extract; + extract.m = sum; + return extract.f[0] + extract.f[1] + extract.f[2] + extract.f[3]; +} + +////////////////////////////////////////////////////////////////////// +// scalar 'a' x + y +// Y = aX + Y +////////////////////////////////////////////////////////////////////// +void CG_SOLVER_SSE::saxpySSE(float s, float* x, float* y) +{ + __m128* ySSE = (__m128*)y; + __m128* xSSE = (__m128*)x; + __m128 sSSE = _mm_set_ps1(s); + __m128 temp; + for (int index = 0; index < _arraySize / 4; index++) + { + temp = _mm_mul_ps(*xSSE, sSSE); + *ySSE = _mm_add_ps(*ySSE, temp); + + xSSE++; + ySSE++; + } +} + +////////////////////////////////////////////////////////////////////// +// scalar 'a' y + x +// Y = aY + X +////////////////////////////////////////////////////////////////////// +void CG_SOLVER_SSE::saypxSSE(float s, float* x, float* y) +{ + __m128* ySSE = (__m128*)y; + __m128* xSSE = (__m128*)x; + __m128 sSSE = _mm_set_ps1(s); + __m128 temp; + for (int index = 0; index < _arraySize / 4; index++) + { + temp = _mm_mul_ps(*ySSE, sSSE); + *ySSE = _mm_add_ps(*xSSE, temp); + + xSSE++; + ySSE++; + } +} + +////////////////////////////////////////////////////////////////////// +// scalar 'a' y + x +// Y = aY + X +////////////////////////////////////////////////////////////////////// +float CG_SOLVER_SSE::maxSSE(float* x) +{ + __m128 maxFoundSSE= _mm_set_ps1(0.0f); + __m128* xSSE = (__m128*)x; + for (int index = 0; index < _arraySize / 4; index++) + { + maxFoundSSE = _mm_max_ps(*xSSE, maxFoundSSE); + xSSE++; + } + union u { + __m128 m; + float f[4]; + } extract; + extract.m = maxFoundSSE; + float maxFound = extract.f[0] > extract.f[1] ? extract.f[0] : extract.f[1]; + maxFound = maxFound > extract.f[2] ? maxFound : extract.f[2]; + return maxFound > extract.f[3] ? maxFound : extract.f[3]; +} + +////////////////////////////////////////////////////////////////////// +// SSE add +// Y = X + Y +////////////////////////////////////////////////////////////////////// +void CG_SOLVER_SSE::addSSE(float* x, float* y) +{ + __m128* ySSE = (__m128*)y; + __m128* xSSE = (__m128*)x; + for (int index = 0; index < _arraySize / 4; index++) + { + *ySSE = _mm_add_ps(*ySSE, *xSSE); + xSSE++; + ySSE++; + } +} + +////////////////////////////////////////////////////////////////////// +// SSE multiply +// Y = X * Y +////////////////////////////////////////////////////////////////////// +void CG_SOLVER_SSE::multiplySSE(float* x, float* y) +{ + __m128* ySSE = (__m128*)y; + __m128* xSSE = (__m128*)x; + for (int index = 0; index < _arraySize / 4; index++) + { + *ySSE = _mm_mul_ps(*ySSE, *xSSE); + xSSE++; + ySSE++; + } +} + +////////////////////////////////////////////////////////////////////// +// SSE multiply +// Z = X * Y +////////////////////////////////////////////////////////////////////// +void CG_SOLVER_SSE::multiplySSE(float* x, float* y, float* z) +{ + __m128* zSSE = (__m128*)z; + __m128* ySSE = (__m128*)y; + __m128* xSSE = (__m128*)x; + for (int index = 0; index < _arraySize / 4; index++) + { + *zSSE = _mm_mul_ps(*ySSE, *xSSE); + xSSE++; + ySSE++; + zSSE++; + } +} + +////////////////////////////////////////////////////////////////////// +// SSE multiply +// Z = W - X * Y +////////////////////////////////////////////////////////////////////// +void CG_SOLVER_SSE::multiplySubtractSSE(float* w, float* x, float* y, float* z) +{ + __m128* zSSE = (__m128*)z; + __m128* ySSE = (__m128*)y; + __m128* xSSE = (__m128*)x; + __m128* wSSE = (__m128*)w; + for (int index = 0; index < _arraySize / 4; index++) + { + *zSSE = _mm_mul_ps(*ySSE, *xSSE); + *zSSE = _mm_sub_ps(*wSSE, *zSSE); + + xSSE++; + ySSE++; + zSSE++; + wSSE++; + } +} + +////////////////////////////////////////////////////////////////////// +// SSE set +// X = val +////////////////////////////////////////////////////////////////////// +void CG_SOLVER_SSE::setSSE(float* x, float val) +{ + __m128* xSSE = (__m128*)x; + for (int index = 0; index < _arraySize / 4; index++) + { + *xSSE = _mm_set_ps1(val); + xSSE++; + } +} + +////////////////////////////////////////////////////////////////////// +// SSE set +// X = 0 +////////////////////////////////////////////////////////////////////// +void CG_SOLVER_SSE::wipeSSE(float* x) +{ + __m128* xSSE = (__m128*)x; + for (int index = 0; index < _arraySize / 4; index++) + { + *xSSE = _mm_setzero_ps(); + xSSE++; + } +} + +////////////////////////////////////////////////////////////////////// +// SSE set +// X = Y +////////////////////////////////////////////////////////////////////// +void CG_SOLVER_SSE::copySSE(float* x, float* y) +{ + __m128* ySSE = (__m128*)y; + for (int index = 0; index < _arraySize / 4; index++) + { + _mm_store_ps(x,*ySSE); + x += 4; + ySSE++; + } +} diff --git a/CG_SOLVER_SSE.h b/CG_SOLVER_SSE.h new file mode 100644 index 0000000..4475944 --- /dev/null +++ b/CG_SOLVER_SSE.h @@ -0,0 +1,103 @@ +/////////////////////////////////////////////////////////////////////////////////// +// File : CG_SOLVER_SSE.h +// +// LumosQuad - A Lightning Generator +// Copyright 2007 +// The University of North Carolina at Chapel Hill +// +/////////////////////////////////////////////////////////////////////////////////// +// +// 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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// The University of North Carolina at Chapel Hill makes no representations +// about the suitability of this software for any purpose. It is provided +// "as is" without express or implied warranty. +// +// Permission to use, copy, modify and distribute this software and its +// documentation for educational, research and non-profit purposes, without +// fee, and without a written agreement is hereby granted, provided that the +// above copyright notice and the following three paragraphs appear in all +// copies. +// +// THE UNIVERSITY OF NORTH CAROLINA SPECIFICALLY DISCLAIM ANY WARRANTIES, +// INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND +// FITNESS FOR A PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN +// "AS IS" BASIS, AND THE UNIVERSITY OF NORTH CAROLINA HAS NO OBLIGATION TO +// PROVIDE MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS. +// +// Please send questions and comments about LumosQuad to kim@cs.unc.edu. +// +/////////////////////////////////////////////////////////////////////////////////// +// +// This program uses OpenEXR, which has the following restrictions: +// +// Copyright (c) 2002, Industrial Light & Magic, a division of Lucas +// Digital Ltd. LLC +// +// All rights reserved. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are +// met: +// * Redistributions of source code must retain the above copyright +// notice, this list of conditions and the following disclaimer. +// * Redistributions in binary form must reproduce the above +// copyright notice, this list of conditions and the following disclaimer +// in the documentation and/or other materials provided with the +// distribution. +// * Neither the name of Industrial Light & Magic nor the names of +// its contributors may be used to endorse or promote products derived +// from this software without specific prior written permission. +// + +#ifndef CG_SOLVER_SSE_H +#define CG_SOLVER_SSE_H + +#include "CG_SOLVER.h" +#include + +//////////////////////////////////////////////////////////////////// +/// \brief Conjugate gradient Poisson solver that exploits SSE. +//////////////////////////////////////////////////////////////////// +class CG_SOLVER_SSE : public CG_SOLVER +{ +public: + //! constructor + CG_SOLVER_SSE(int maxDepth, int iterations = 10, int digits = 1); + //! destructor + ~CG_SOLVER_SSE(); + + //! solve the Poisson problem using SSE + virtual int solve(list cells); + +private: + //! reallocate the SSE-friendly scratch arrays + virtual void reallocate(); + + // SSE linear algebra operators + inline float dotSSE(float* x, float* y); + inline void saxpySSE(float a, float* x, float* y); + inline void saypxSSE(float a, float* x, float* y); + inline float maxSSE(float* x); + inline void addSSE(float* x, float* y); + inline void multiplySSE(float* x, float* y); + inline void multiplySSE(float* x, float* y, float* z); + inline void multiplySubtractSSE(float* w, float* x, float* y, float* z); + inline void setSSE(float* x, float val); + inline void wipeSSE(float* x); + inline void copySSE(float* x, float* y); +}; + +#endif diff --git a/DAG.cpp b/DAG.cpp new file mode 100644 index 0000000..a2d7d28 --- /dev/null +++ b/DAG.cpp @@ -0,0 +1,566 @@ +/////////////////////////////////////////////////////////////////////////////////// +// File : DAG.cpp +/////////////////////////////////////////////////////////////////////////////////// +// +// LumosQuad - A Lightning Generator +// Copyright 2007 +// The University of North Carolina at Chapel Hill +// +/////////////////////////////////////////////////////////////////////////////////// +// +// 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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// The University of North Carolina at Chapel Hill makes no representations +// about the suitability of this software for any purpose. It is provided +// "as is" without express or implied warranty. +// +// Permission to use, copy, modify and distribute this software and its +// documentation for educational, research and non-profit purposes, without +// fee, and without a written agreement is hereby granted, provided that the +// above copyright notice and the following three paragraphs appear in all +// copies. +// +// THE UNIVERSITY OF NORTH CAROLINA SPECIFICALLY DISCLAIM ANY WARRANTIES, +// INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND +// FITNESS FOR A PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN +// "AS IS" BASIS, AND THE UNIVERSITY OF NORTH CAROLINA HAS NO OBLIGATION TO +// PROVIDE MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS. +// +// Please send questions and comments about LumosQuad to kim@cs.unc.edu. +// +/////////////////////////////////////////////////////////////////////////////////// +// +// This program uses OpenEXR, which has the following restrictions: +// +// Copyright (c) 2002, Industrial Light & Magic, a division of Lucas +// Digital Ltd. LLC +// +// All rights reserved. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are +// met: +// * Redistributions of source code must retain the above copyright +// notice, this list of conditions and the following disclaimer. +// * Redistributions in binary form must reproduce the above +// copyright notice, this list of conditions and the following disclaimer +// in the documentation and/or other materials provided with the +// distribution. +// * Neither the name of Industrial Light & Magic nor the names of +// its contributors may be used to endorse or promote products derived +// from this software without specific prior written permission. +// + +#include "DAG.h" + +////////////////////////////////////////////////////////////////////// +// Construction/Destruction +////////////////////////////////////////////////////////////////////// + +DAG::DAG(int xRes, int yRes) : + _xRes(xRes), + _yRes(yRes), + _width(xRes), + _height(yRes), + _dx(1.0f / xRes), + _dy(1.0f / yRes), + _root(NULL), + _totalNodes(0), + _bottomHit(-1), + _secondaryIntensity(0.3f), + _leaderIntensity(0.75f) +{ + (_dx < _dy) ? _dy = _dx : _dy = _dx; + _offscreenBuffer = NULL; +} + +DAG::~DAG() +{ + if (_offscreenBuffer) delete[] _offscreenBuffer; + + deleteNode(_root); +} + +////////////////////////////////////////////////////////////////////// +// delete the nodes +////////////////////////////////////////////////////////////////////// +void DAG::deleteNode(NODE* root) +{ + if (root == NULL) return; + + for (int x = 0; x < root->neighbors.size(); x++) + deleteNode(root->neighbors[x]); + + delete root; + root = NULL; +} + +////////////////////////////////////////////////////////////////////// +// add line segment 'index' to segment list +////////////////////////////////////////////////////////////////////// +bool DAG::addSegment(int index, int neighbor) +{ + if (_root) + { + // find corresponding root in DAG + map::iterator iter = _hash.find(neighbor); + NODE* root = (*iter).second; + + if (root == NULL) return false; + + // add to DAG + NODE* newNode = new NODE(index); + newNode->parent = root; + root->neighbors.push_back(newNode); + + // add to hash table + _hash.insert(map::value_type(index, newNode)); + } + else + { + // make the root + _root = new NODE(neighbor); + _hash.insert(map::value_type(neighbor, _root)); + + // then do the add + NODE* newNode = new NODE(index); + newNode->parent = _root; + _root->neighbors.push_back(newNode); + + // add to hash table + _hash.insert(map::value_type(index, newNode)); + } + + _totalNodes++; + + return true; +} + +////////////////////////////////////////////////////////////////////// +// build the leader chain +////////////////////////////////////////////////////////////////////// +void DAG::buildLeader(int bottomHit) +{ + _bottomHit = bottomHit; + + // get pointer to bottommost node + map::iterator iter = _hash.find(bottomHit); + NODE* child = (*iter).second; + + // crawl up the tree + while (child != NULL) + { + // tag segment + child->leader = true; + child->secondary = false; + + // look for side branches + for (int x = 0; x < child->neighbors.size(); x++) + if (!(child->neighbors[x]->leader)) + buildBranch(child->neighbors[x], 1); + + // advance child + child = child->parent; + } + buildIntensity(_root); +} + +////////////////////////////////////////////////////////////////////// +// build the side branch +////////////////////////////////////////////////////////////////////// +void DAG::buildBranch(NODE* node, int depth) +{ + node->depth = depth; + node->leader = false; + + // look for side branches + for (int x = 0; x < node->neighbors.size(); x++) + if (!(node->neighbors[x]->leader)) + buildBranch(node->neighbors[x], depth + 1); +} + +////////////////////////////////////////////////////////////////////// +// draw the DAG segments +////////////////////////////////////////////////////////////////////// +void DAG::drawNode(NODE* root) +{ + if (root == NULL) return; + + // draw segments + int beginIndex = root->index; + int begin[2]; + int end[2]; + float dWidth = _dx; + float dHeight = _dy; + + for (int x = 0; x < root->neighbors.size(); x++) + { + // get end node + NODE* endNode = root->neighbors[x]; + int endIndex = endNode->index; + + // draw segments + begin[1] = beginIndex / _xRes; + begin[0] = beginIndex - begin[1] * _xRes; + end[1] = endIndex / _xRes; + end[0] = endIndex - end[1] * _xRes; + + if (_bottomHit != -1) + glColor4f(endNode->intensity, + endNode->intensity, + endNode->intensity,1.0f); + else + glColor4f(0.0f, 1.0f, 0.0f, 1.0f); + + glBegin(GL_LINES); + glVertex3f(begin[0] * dWidth + dWidth * 0.5f, 1.0f - begin[1] * dHeight + dHeight * 0.5f, 0.1f); + glVertex3f(end[0] * dWidth + dWidth * 0.5f, 1.0f - end[1] * dHeight + dHeight * 0.5f, 0.1f); + glEnd(); + + // call recursively + if (endNode->neighbors.size() > 0) + drawNode(endNode); + } +} + +////////////////////////////////////////////////////////////////////// +// set the intensity of the node +////////////////////////////////////////////////////////////////////// +void DAG::buildIntensity(NODE* root) +{ + if (root == NULL) return; + + // draw segments + int beginIndex = root->index; + int begin[2]; + int end[2]; + float dWidth = _dx; + float dHeight = _dy; + + for (int x = 0; x < root->neighbors.size(); x++) + { + // get end node + NODE* endNode = root->neighbors[x]; + int endIndex = endNode->index; + + // draw segments + begin[1] = beginIndex / _xRes; + begin[0] = beginIndex - begin[1] * _xRes; + end[1] = endIndex / _xRes; + end[0] = endIndex - end[1] * _xRes; + + // set color + if (endNode->leader) + endNode->intensity = _leaderIntensity; + else + { + // find max depth of current channel + if (endNode->maxDepthNode == NULL) + findDeepest(endNode, endNode->maxDepthNode); + + int maxDepth = endNode->maxDepthNode->depth; + + // calc standard deviation + float stdDev = -(float)(maxDepth * maxDepth) / (float)(log(_secondaryIntensity) * 2.0f); + + // calc falloff + float eTerm = -(float)(endNode->depth) * (float)(endNode->depth); + eTerm /= (2.0f * stdDev); + eTerm = exp(eTerm) * 0.5f; + endNode->intensity = eTerm; + } + + // call recursively + if (endNode->neighbors.size() > 0) + buildIntensity(endNode); + } +} + +////////////////////////////////////////////////////////////////////// +// draw the tree offscreen +////////////////////////////////////////////////////////////////////// +float*& DAG::drawOffscreen(int scale) +{ + // allocate buffer + _width = _xRes * scale; + _height = _yRes * scale; + _scale = scale; + if (_offscreenBuffer) delete[] _offscreenBuffer; + _offscreenBuffer = new float[_width * _height]; + + // wipe the buffer + for (int x = 0; x < _width * _height; x++) + _offscreenBuffer[x] = 0.0f; + + // recursively draw the tree + drawOffscreenNode(_root); + + return _offscreenBuffer; +} + +////////////////////////////////////////////////////////////////////// +// draw the DAG segments +////////////////////////////////////////////////////////////////////// +void DAG::drawOffscreenNode(NODE* root) +{ + if (root == NULL) return; + + // draw segments + int beginIndex = root->index; + int begin[2]; + int end[2]; + float dWidth = _dx; + float dHeight = _dy; + + for (int x = 0; x < root->neighbors.size(); x++) + { + // get end node + NODE* endNode = root->neighbors[x]; + int endIndex = endNode->index; + + // get endpoints + begin[0] = beginIndex % _xRes * _scale; + begin[1] = beginIndex / _xRes * _scale; + end[0] = endIndex % _xRes * _scale; + end[1] = endIndex / _xRes * _scale; + + // make sure the one with the smaller x comes first + if (end[0] < begin[0]) + { + int temp[] = {end[0], end[1]}; + end[0] = begin[0]; + end[1] = begin[1]; + + begin[0] = temp[0]; + begin[1] = temp[1]; + } + + // rasterize + drawLine(begin, end, endNode->intensity); + + // call recursively + if (endNode->neighbors.size() > 0) + drawOffscreenNode(endNode); + } +} + +////////////////////////////////////////////////////////////////////// +// rasterize the line +// I assume all the lines are purely horizontal, vertical or diagonal +// to avoid using something like Bresenham +////////////////////////////////////////////////////////////////////// +void DAG::drawLine(int begin[], int end[], float intensity) +{ + int scaledX = _xRes * _scale; + + // if it is a horizontal line + if (begin[1] == end[1]) + { + for (int x = begin[0]; x < end[0]; x++) + { + int index = x + end[1] * scaledX; + if (intensity > _offscreenBuffer[index]) + _offscreenBuffer[index] = intensity; + } + return; + } + + // if it is a vertical line + if (begin[0] == end[0]) + { + int bottom = (begin[1] > end[1]) ? end[1] : begin[1]; + int top = (begin[1] > end[1]) ? begin[1] : end[1]; + + for (int y = bottom; y < top; y++) + { + int index = begin[0] + y * scaledX; + if (intensity > _offscreenBuffer[index]) + _offscreenBuffer[index] = intensity; + } + return; + } + + // else it is diagonal + int slope = (begin[1] < end[1]) ? 1 : -1; + int interval = end[0] - begin[0]; + for (int x = 0; x <= interval; x++) + { + int index = begin[0] + x + (begin[1] + x * slope) * scaledX; + if (intensity > _offscreenBuffer[index]) + _offscreenBuffer[index] = intensity; + } +} + +////////////////////////////////////////////////////////////////////// +// find deepest depth from current root +////////////////////////////////////////////////////////////////////// +void DAG::findDeepest(NODE* root, NODE*& deepest) +{ + deepest = root; + + for (int x = 0; x < root->neighbors.size(); x++) + { + NODE* child = root->neighbors[x]; + NODE* candidate = NULL; + findDeepest(child, candidate); + if (candidate->depth > deepest->depth) + deepest = candidate; + } +} + +////////////////////////////////////////////////////////////////////// +// dump out line segments +////////////////////////////////////////////////////////////////////// +void DAG::write(const char* filename) +{ + // open file + FILE* file; + file = fopen(filename, "wb"); + + // write out total number of DAG nodes + fwrite((void*)&_totalNodes, sizeof(int), 1, file); + fwrite((void*)&_xRes, sizeof(int), 1, file); + fwrite((void*)&_yRes, sizeof(int), 1, file); + fwrite((void*)&_dx, sizeof(float), 1, file); + fwrite((void*)&_dy, sizeof(float), 1, file); + fwrite((void*)&_bottomHit, sizeof(int), 1, file); + fwrite((void*)&_inputWidth, sizeof(int), 1, file); + fwrite((void*)&_inputHeight, sizeof(int), 1, file); + + // write out nodes + writeNode(_root, file); + + fclose(file); +} + +////////////////////////////////////////////////////////////////////// +// read in line segments +////////////////////////////////////////////////////////////////////// +void DAG::read(const char* filename) +{ + + // erase old DAG + deleteNode(_root); + + // open file + FILE* file; + file = fopen(filename, "rb"); + if (file == NULL) + { + cout << "ERROR: " << filename << " is invalid." << endl; + exit(1); + } + + // read in total number of DAG nodes + fread((void*)&_totalNodes, sizeof(int), 1, file); + fread((void*)&_xRes, sizeof(int), 1, file); + fread((void*)&_yRes, sizeof(int), 1, file); + fread((void*)&_dx, sizeof(float), 1, file); + fread((void*)&_dy, sizeof(float), 1, file); + fread((void*)&_bottomHit, sizeof(int), 1, file); + fread((void*)&_inputWidth, sizeof(int), 1, file); + fread((void*)&_inputHeight, sizeof(int), 1, file); + + // clear _hash for use + _hash.clear(); + + // read in all the DAG nodes + for (int x = 0; x <= _totalNodes; x++) + readNode(file); + + fclose(file); + + if (_bottomHit != -1) + buildLeader(_bottomHit); +} + +////////////////////////////////////////////////////////////////////// +// write out a DAG node +////////////////////////////////////////////////////////////////////// +void DAG::writeNode(NODE* root, FILE* file) +{ + int x; + + // write out neighbors + int numNeighbors = root->neighbors.size(); + for (x = 0; x < numNeighbors; x++) + writeNode(root->neighbors[x], file); + + // write out this node + fwrite((void*)&(root->index), sizeof(int), 1, file); + int parent = (root->parent == NULL) ? -1 : root->parent->index; + fwrite((void*)&parent, sizeof(int), 1, file); + fwrite((void*)&(root->leader), sizeof(bool), 1, file); + fwrite((void*)&(root->secondary), sizeof(bool), 1, file); + fwrite((void*)&(root->depth), sizeof(int), 1, file); + fwrite((void*)&numNeighbors, sizeof(int), 1, file); + for (x = 0; x < numNeighbors; x++) + fwrite((void*)&(root->neighbors[x]->index), sizeof(int), 1, file); +} + +////////////////////////////////////////////////////////////////////// +// read in a DAG node +////////////////////////////////////////////////////////////////////// +void DAG::readNode(FILE* file) +{ + // read in DAG data + int index; + int parent; + bool leader; + bool secondary; + int depth; + int numNeighbors; + float potential; + fread((void*)&index, sizeof(int), 1, file); + fread((void*)&parent, sizeof(int), 1, file); + fread((void*)&leader, sizeof(bool), 1, file); + fread((void*)&secondary, sizeof(bool), 1, file); + fread((void*)&depth, sizeof(int), 1, file); + fread((void*)&numNeighbors, sizeof(int), 1, file); + + // create the node + NODE* node = new NODE(index); + node->leader = leader; + node->secondary = secondary; + node->depth = depth; + + // look up neighbors + map::iterator iter; + for (int x = 0; x < numNeighbors; x++) + { + // read in the child index + int neighborIndex; + fread((void*)&neighborIndex, sizeof(int), 1, file); + + // look up in hash table + iter = _hash.find(neighborIndex); + + // push onto the neighbor vector + node->neighbors.push_back((*iter).second); + + // set child's parent + (*iter).second->parent = node; + } + + // add to hash table + _hash.insert(map::value_type(index, node)); + + // search for root node + if (parent == -1) + { + node->parent = NULL; + _root = node; + } +} diff --git a/DAG.h b/DAG.h new file mode 100644 index 0000000..9ba3a9f --- /dev/null +++ b/DAG.h @@ -0,0 +1,208 @@ +/////////////////////////////////////////////////////////////////////////////////// +// File : DAG.h +/////////////////////////////////////////////////////////////////////////////////// +// +// LumosQuad - A Lightning Generator +// Copyright 2007 +// The University of North Carolina at Chapel Hill +// +/////////////////////////////////////////////////////////////////////////////////// +// +// 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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// The University of North Carolina at Chapel Hill makes no representations +// about the suitability of this software for any purpose. It is provided +// "as is" without express or implied warranty. +// +// Permission to use, copy, modify and distribute this software and its +// documentation for educational, research and non-profit purposes, without +// fee, and without a written agreement is hereby granted, provided that the +// above copyright notice and the following three paragraphs appear in all +// copies. +// +// THE UNIVERSITY OF NORTH CAROLINA SPECIFICALLY DISCLAIM ANY WARRANTIES, +// INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND +// FITNESS FOR A PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN +// "AS IS" BASIS, AND THE UNIVERSITY OF NORTH CAROLINA HAS NO OBLIGATION TO +// PROVIDE MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS. +// +// Please send questions and comments about LumosQuad to kim@cs.unc.edu. +// +/////////////////////////////////////////////////////////////////////////////////// +// +// This program uses OpenEXR, which has the following restrictions: +// +// Copyright (c) 2002, Industrial Light & Magic, a division of Lucas +// Digital Ltd. LLC +// +// All rights reserved. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are +// met: +// * Redistributions of source code must retain the above copyright +// notice, this list of conditions and the following disclaimer. +// * Redistributions in binary form must reproduce the above +// copyright notice, this list of conditions and the following disclaimer +// in the documentation and/or other materials provided with the +// distribution. +// * Neither the name of Industrial Light & Magic nor the names of +// its contributors may be used to endorse or promote products derived +// from this software without specific prior written permission. +// + +#ifndef DAG_H +#define DAG_H + +#include +#include +#include +#include +#include + +using namespace std; + +//////////////////////////////////////////////////////////////////// +/// \brief Directed acyclic graph that renders the final lightning +//////////////////////////////////////////////////////////////////// +class DAG +{ +public: + //! constructor + DAG(int xRes, int yRes); + //! destructor + virtual ~DAG(); + + //! build the stepped ladder + void buildLeader(int bottomHit); + + //! add DAG segment + bool addSegment(int index, int neighbor); + + //! draw to OpenGL + void draw() { drawNode(_root); }; + + //! draw to an offscreen buffer + float*& drawOffscreen(int scale = 1); + + //! read in a new DAG + void read(const char* filename); + //! write out the current DAG + void write(const char* filename); + + //! quadtree x resolution accessor + int xRes() { return _xRes; }; + //! quadtree y resolution accessor + int yRes() { return _yRes; }; + + //! input image x resolution accessor + int& inputWidth() { return _inputWidth; }; + //! input image y resolution accessor + int& inputHeight() { return _inputHeight; }; + +private: + //! x resolution of quadtree + int _xRes; + //! y resolution of quadtree + int _yRes; + //! physical length of one grid cell + float _dx; + //! physical length of one grid cell + float _dy; + + //////////////////////////////////////////////////////////////////// + /// \brief node for line segment tree + //////////////////////////////////////////////////////////////////// + struct NODE { + int index; + vector neighbors; + NODE* parent; + bool leader; + bool secondary; + int depth; + NODE* maxDepthNode; + float intensity; + + NODE(int indexIn) { + index = indexIn; + parent = NULL; + leader = false; + depth = 0; + maxDepthNode = NULL; + }; + }; + //! recursive destructor + void deleteNode(NODE* root); + + //! root of DAG + NODE* _root; + + //! hash table of DAG nodes + map _hash; + + //! build side branch + void buildBranch(NODE* node, int depth); + //! draw a node to OpenGL + void drawNode(NODE* root); + //! find the deepest node in a given subtree + void findDeepest(NODE* root, NODE*& deepest); + + //! read in a DAG node + void readNode(FILE* file); + //! write out a DAG node + void writeNode(NODE* root, FILE* file); + + //! total number of nodes in scene + int _totalNodes; + + //! node that finally hit bottom + int _bottomHit; + + //! set the line segment intensities + void buildIntensity(NODE* root); + + //! brightness of secondary branch + float _secondaryIntensity; + //! brightness of primary branch + float _leaderIntensity; + + //////////////////////////////////////////////////////////////// + // offscreen buffer variables + //////////////////////////////////////////////////////////////// + //! offscreen buffer + float* _offscreenBuffer; + + //! width of offscreen buffer + int _width; + + //! height of offscreen buffer + int _height; + + //! scale of offscreen buffer compared to original image + int _scale; + + //! draw a given node in the DAG + void drawOffscreenNode(NODE* root); + + //! rasterize a single line to the offscreen buffer + void drawLine(int begin[], int end[], float intensity); + + //! input image x resolution + int _inputWidth; + //! input image y resolution + int _inputHeight; +}; + +#endif diff --git a/EXR.cpp b/EXR.cpp new file mode 100644 index 0000000..8cdbabd --- /dev/null +++ b/EXR.cpp @@ -0,0 +1,96 @@ +/////////////////////////////////////////////////////////////////////////////////// +// File : EXR.cpp +/////////////////////////////////////////////////////////////////////////////////// +// +// LumosQuad - A Lightning Generator +// Copyright 2007 +// The University of North Carolina at Chapel Hill +// +/////////////////////////////////////////////////////////////////////////////////// +// +// 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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// The University of North Carolina at Chapel Hill makes no representations +// about the suitability of this software for any purpose. It is provided +// "as is" without express or implied warranty. +// +// Permission to use, copy, modify and distribute this software and its +// documentation for educational, research and non-profit purposes, without +// fee, and without a written agreement is hereby granted, provided that the +// above copyright notice and the following three paragraphs appear in all +// copies. +// +// THE UNIVERSITY OF NORTH CAROLINA SPECIFICALLY DISCLAIM ANY WARRANTIES, +// INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND +// FITNESS FOR A PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN +// "AS IS" BASIS, AND THE UNIVERSITY OF NORTH CAROLINA HAS NO OBLIGATION TO +// PROVIDE MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS. +// +// Please send questions and comments about LumosQuad to kim@cs.unc.edu. +// +/////////////////////////////////////////////////////////////////////////////////// +// +// This program uses OpenEXR, which has the following restrictions: +// +// Copyright (c) 2002, Industrial Light & Magic, a division of Lucas +// Digital Ltd. LLC +// +// All rights reserved. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are +// met: +// * Redistributions of source code must retain the above copyright +// notice, this list of conditions and the following disclaimer. +// * Redistributions in binary form must reproduce the above +// copyright notice, this list of conditions and the following disclaimer +// in the documentation and/or other materials provided with the +// distribution. +// * Neither the name of Industrial Light & Magic nor the names of +// its contributors may be used to endorse or promote products derived +// from this software without specific prior written permission. +// + +#include "EXR.h" + +////////////////////////////////////////////////////////////////////// +// Construction/Destruction +////////////////////////////////////////////////////////////////////// + +EXR::EXR() +{ + +} + +EXR::~EXR() +{ + +} + +void EXR::writeEXR(const char filename[], float* image, int width, int height) +{ + Imf::Rgba* exrImage = new Imf::Rgba[width * height]; + for (int x = 0; x < width * height; x++) + { + exrImage[x].r = exrImage[x].g = exrImage[x].b = image[x]; + exrImage[x].a = 1.0f; + } + + Imf::RgbaOutputFile file (filename, width, height); + file.setFrameBuffer (exrImage, 1, width); + file.writePixels (height); + + delete[] exrImage; +} diff --git a/EXR.h b/EXR.h new file mode 100644 index 0000000..5dbdca5 --- /dev/null +++ b/EXR.h @@ -0,0 +1,92 @@ +/////////////////////////////////////////////////////////////////////////////////// +// File : EXR.h +/////////////////////////////////////////////////////////////////////////////////// +// +// LumosQuad - A Lightning Generator +// Copyright 2007 +// The University of North Carolina at Chapel Hill +// +/////////////////////////////////////////////////////////////////////////////////// +// +// 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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// The University of North Carolina at Chapel Hill makes no representations +// about the suitability of this software for any purpose. It is provided +// "as is" without express or implied warranty. +// +// Permission to use, copy, modify and distribute this software and its +// documentation for educational, research and non-profit purposes, without +// fee, and without a written agreement is hereby granted, provided that the +// above copyright notice and the following three paragraphs appear in all +// copies. +// +// THE UNIVERSITY OF NORTH CAROLINA SPECIFICALLY DISCLAIM ANY WARRANTIES, +// INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND +// FITNESS FOR A PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN +// "AS IS" BASIS, AND THE UNIVERSITY OF NORTH CAROLINA HAS NO OBLIGATION TO +// PROVIDE MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS. +// +// Please send questions and comments about LumosQuad to kim@cs.unc.edu. +// +/////////////////////////////////////////////////////////////////////////////////// +// +// This program uses OpenEXR, which has the following restrictions: +// +// Copyright (c) 2002, Industrial Light & Magic, a division of Lucas +// Digital Ltd. LLC +// +// All rights reserved. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are +// met: +// * Redistributions of source code must retain the above copyright +// notice, this list of conditions and the following disclaimer. +// * Redistributions in binary form must reproduce the above +// copyright notice, this list of conditions and the following disclaimer +// in the documentation and/or other materials provided with the +// distribution. +// * Neither the name of Industrial Light & Magic nor the names of +// its contributors may be used to endorse or promote products derived +// from this software without specific prior written permission. +// + +#ifndef EXR_H +#define EXR_H +#include +#include +#include +#include +#include + +//////////////////////////////////////////////////////////////////// +/// \brief Wrapper for the ILM OpenEXR routines +//////////////////////////////////////////////////////////////////// +class EXR +{ +public: + EXR(); + virtual ~EXR(); + + /// \brief write a float array out to an EXR file + /// + /// \param filename name of output file + /// \param image float array to write to an EXR file + /// \param width width of float array + /// \param height height of float array + static void writeEXR(const char* filename, float* image, int width, int height); +}; + +#endif diff --git a/FFT.cpp b/FFT.cpp new file mode 100644 index 0000000..34bfb8a --- /dev/null +++ b/FFT.cpp @@ -0,0 +1,233 @@ +/////////////////////////////////////////////////////////////////////////////////// +// File : FFT.cpp +/////////////////////////////////////////////////////////////////////////////////// +// +// LumosQuad - A Lightning Generator +// Copyright 2007 +// The University of North Carolina at Chapel Hill +// +/////////////////////////////////////////////////////////////////////////////////// +// +// 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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// The University of North Carolina at Chapel Hill makes no representations +// about the suitability of this software for any purpose. It is provided +// "as is" without express or implied warranty. +// +// Permission to use, copy, modify and distribute this software and its +// documentation for educational, research and non-profit purposes, without +// fee, and without a written agreement is hereby granted, provided that the +// above copyright notice and the following three paragraphs appear in all +// copies. +// +// THE UNIVERSITY OF NORTH CAROLINA SPECIFICALLY DISCLAIM ANY WARRANTIES, +// INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND +// FITNESS FOR A PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN +// "AS IS" BASIS, AND THE UNIVERSITY OF NORTH CAROLINA HAS NO OBLIGATION TO +// PROVIDE MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS. +// +// Please send questions and comments about LumosQuad to kim@cs.unc.edu. +// +/////////////////////////////////////////////////////////////////////////////////// +// +// This program uses OpenEXR, which has the following restrictions: +// +// Copyright (c) 2002, Industrial Light & Magic, a division of Lucas +// Digital Ltd. LLC +// +// All rights reserved. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are +// met: +// * Redistributions of source code must retain the above copyright +// notice, this list of conditions and the following disclaimer. +// * Redistributions in binary form must reproduce the above +// copyright notice, this list of conditions and the following disclaimer +// in the documentation and/or other materials provided with the +// distribution. +// * Neither the name of Industrial Light & Magic nor the names of +// its contributors may be used to endorse or promote products derived +// from this software without specific prior written permission. +// + +#include "FFT.h" + +////////////////////////////////////////////////////////////////////// +// Construction/Destruction +////////////////////////////////////////////////////////////////////// + +FFT::FFT() +{ + +} + +FFT::~FFT() +{ + +} + +bool FFT::convolve(float* source, float* kernel, int xSource, int ySource, int xKernel, int yKernel) +{ + int x, y, index; + + // get normalization params + float maxCurrent = 0.0f; + for (x = 0; x < xSource * ySource; x++) + maxCurrent = (maxCurrent < source[x]) ? source[x] : maxCurrent; + float maxKernel = 0.0f; + for (x = 0; x < xKernel * yKernel; x++) + maxKernel = (maxKernel < kernel[x]) ? kernel[x] : maxKernel; + float maxProduct = maxCurrent * maxKernel; + + // retrieve dimensions + int xHalf = xKernel / 2; + int yHalf = yKernel / 2; + int xResPadded = xSource + xKernel; + int yResPadded = ySource + yKernel; + + if (xResPadded != yResPadded) + (xResPadded > yResPadded) ? yResPadded = xResPadded : xResPadded = yResPadded; + + // create padded field + fftw_complex* padded = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * xResPadded * yResPadded); + if (!padded) + { + cout << " IMAGE: Not enough memory! Try a smaller final image size." << endl; + return false; + } + + // init padded field + for (index = 0; index < xResPadded * yResPadded; index++) + padded[index][0] = padded[index][1] = 0.0f; + index = 0; + for (y = 0; y < ySource; y++) + for (x = 0; x < xSource; x++, index++) + { + int paddedIndex = (x + xKernel / 2) + (y + yKernel / 2) * xResPadded; + padded[paddedIndex][0] = source[index]; + } + + // create padded filter + fftw_complex* filter = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * xResPadded * yResPadded); + if (!filter) + { + cout << " FILTER: Not enough memory! Try a smaller final image size." << endl; + return false; + } + + // init padded filter + for (index = 0; index < xResPadded * yResPadded; index++) + filter[index][0] = filter[index][1] = 0.0f; + + // quadrant IV + for (y = 0; y < (yHalf + 1); y++) + for (x = 0; x < (xHalf + 1); x++) + { + int filterIndex = x + xHalf + y * xKernel; + int fieldIndex = x + (y + yResPadded - (yHalf + 1)) * xResPadded; + filter[fieldIndex][0] = kernel[filterIndex]; + } + + // quadrant I + for (y = 0; y < yHalf; y++) + for (x = 0; x < (xHalf + 1); x++) + { + int filterIndex = (x + xHalf) + (y + yHalf + 1) * xKernel; + int fieldIndex = x + y * xResPadded; + filter[fieldIndex][0] = filter[fieldIndex][1] = kernel[filterIndex]; + } + + // quadrant III + for (y = 0; y < (yHalf + 1); y++) + for (x = 0; x < xHalf; x++) + { + int filterIndex = x + y * xKernel; + int fieldIndex = (x + xResPadded - xHalf) + (y + yResPadded - (yHalf + 1)) * xResPadded; + filter[fieldIndex][0] = filter[fieldIndex][1] = kernel[filterIndex]; + } + + // quadrant II + for (y = 0; y < yHalf; y++) + for (x = 0; x < xHalf; x++) + { + int filterIndex = x + (y + yHalf + 1) * xKernel; + int fieldIndex = (x + xResPadded - xHalf) + y * xResPadded; + filter[fieldIndex][0] = filter[fieldIndex][1] = kernel[filterIndex]; + } + + // perform forward FFT on field + fftw_complex* paddedTransformed = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * xResPadded * yResPadded); + if (!paddedTransformed) + { + cout << " T-IMAGE: Not enough memory! Try a smaller final image size." << endl; + return false; + } + fftw_plan forwardField = fftw_plan_dft_2d(xResPadded, yResPadded, padded, paddedTransformed, FFTW_FORWARD, FFTW_ESTIMATE); + fftw_execute(forwardField); + + // perform forward FFT on filter + fftw_complex* filterTransformed = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * xResPadded * yResPadded); + if (!filterTransformed) + { + cout << " T-FILTER: Not enough memory! Try a smaller final image size." << endl; + return false; + } + fftw_plan forwardFilter = fftw_plan_dft_2d(xResPadded, yResPadded, filter, filterTransformed, FFTW_FORWARD, FFTW_ESTIMATE); + fftw_execute(forwardFilter); + + // apply frequency space filter + for (index = 0; index < xResPadded * yResPadded; index++) + { + float newReal = paddedTransformed[index][0] * filterTransformed[index][0] - + paddedTransformed[index][1] * filterTransformed[index][1]; + float newIm = paddedTransformed[index][0] * filterTransformed[index][1] + + paddedTransformed[index][1] * filterTransformed[index][0]; + paddedTransformed[index][0] = newReal; + paddedTransformed[index][1] = newIm; + } + + // transform back + fftw_plan backwardField = fftw_plan_dft_2d(xResPadded, yResPadded, paddedTransformed, padded, FFTW_BACKWARD, FFTW_ESTIMATE); + fftw_execute(backwardField); + + // copy back into padded + index = 0; + for (y = 0; y < ySource; y++) + for (x = 0; x < xSource; x++, index++) + { + int paddedIndex = (x + xKernel / 2) + (y + yKernel / 2) * xResPadded; + source[index] = padded[paddedIndex][0]; + } + + // clean up + fftw_free(padded); + fftw_free(paddedTransformed); + fftw_free(filter); + fftw_free(filterTransformed); + + // if normalization is exceeded, renormalize + float newMax = 0.0f; + for (x = 0; x < xSource * ySource; x++) + newMax = (newMax < source[x]) ? source[x] : newMax; + if (newMax > maxProduct) + { + float scale = maxProduct / newMax; + for (x = 0; x < xSource * ySource; x++) + source[x] *= scale; + } + + return true; +} diff --git a/FFT.h b/FFT.h new file mode 100644 index 0000000..58205f6 --- /dev/null +++ b/FFT.h @@ -0,0 +1,95 @@ +/////////////////////////////////////////////////////////////////////////////////// +// File : FFT.h +/////////////////////////////////////////////////////////////////////////////////// +// +// LumosQuad - A Lightning Generator +// Copyright 2007 +// The University of North Carolina at Chapel Hill +// +/////////////////////////////////////////////////////////////////////////////////// +// +// 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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// The University of North Carolina at Chapel Hill makes no representations +// about the suitability of this software for any purpose. It is provided +// "as is" without express or implied warranty. +// +// Permission to use, copy, modify and distribute this software and its +// documentation for educational, research and non-profit purposes, without +// fee, and without a written agreement is hereby granted, provided that the +// above copyright notice and the following three paragraphs appear in all +// copies. +// +// THE UNIVERSITY OF NORTH CAROLINA SPECIFICALLY DISCLAIM ANY WARRANTIES, +// INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND +// FITNESS FOR A PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN +// "AS IS" BASIS, AND THE UNIVERSITY OF NORTH CAROLINA HAS NO OBLIGATION TO +// PROVIDE MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS. +// +// Please send questions and comments about LumosQuad to kim@cs.unc.edu. +// +/////////////////////////////////////////////////////////////////////////////////// +// +// This program uses OpenEXR, which has the following restrictions: +// +// Copyright (c) 2002, Industrial Light & Magic, a division of Lucas +// Digital Ltd. LLC +// +// All rights reserved. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are +// met: +// * Redistributions of source code must retain the above copyright +// notice, this list of conditions and the following disclaimer. +// * Redistributions in binary form must reproduce the above +// copyright notice, this list of conditions and the following disclaimer +// in the documentation and/or other materials provided with the +// distribution. +// * Neither the name of Industrial Light & Magic nor the names of +// its contributors may be used to endorse or promote products derived +// from this software without specific prior written permission. +// + +#ifndef FFT_H +#define FFT_H + +#include +#include + +using namespace std; + +//////////////////////////////////////////////////////////////////// +/// \brief Wrapper class for FFTW +//////////////////////////////////////////////////////////////////// +class FFT +{ +public: + FFT(); + virtual ~FFT(); + /// \brief convolve image and filter using FFTW + /// + /// \param source source image + /// \param kernel convolution kernel + /// \param xSource width of source image + /// \param ySource height of source image + /// \param xKernel width of kernel + /// \param yKernel height yidth of kernel + /// + /// \return Returns the convolved image in the 'image' array. If the convolve fails, returns false + static bool convolve(float* source, float* kernel, int xSource, int ySource, int xKernel, int yKernel); +}; + +#endif diff --git a/LumosQuad.sln b/LumosQuad.sln new file mode 100644 index 0000000..a461672 --- /dev/null +++ b/LumosQuad.sln @@ -0,0 +1,20 @@ + +Microsoft Visual Studio Solution File, Format Version 9.00 +# Visual C++ Express 2005 +Project("{8BC9CEB8-8B4A-11D0-8D11-00A0C91BC942}") = "LumosQuad", "LumosQuad.vcproj", "{5BC95B61-43C1-46C2-BA59-D95F6C336355}" +EndProject +Global + GlobalSection(SolutionConfigurationPlatforms) = preSolution + Debug|Win32 = Debug|Win32 + Release|Win32 = Release|Win32 + EndGlobalSection + GlobalSection(ProjectConfigurationPlatforms) = postSolution + {5BC95B61-43C1-46C2-BA59-D95F6C336355}.Debug|Win32.ActiveCfg = Debug|Win32 + {5BC95B61-43C1-46C2-BA59-D95F6C336355}.Debug|Win32.Build.0 = Debug|Win32 + {5BC95B61-43C1-46C2-BA59-D95F6C336355}.Release|Win32.ActiveCfg = Release|Win32 + {5BC95B61-43C1-46C2-BA59-D95F6C336355}.Release|Win32.Build.0 = Release|Win32 + EndGlobalSection + GlobalSection(SolutionProperties) = preSolution + HideSolutionNode = FALSE + EndGlobalSection +EndGlobal diff --git a/LumosQuad.vcproj b/LumosQuad.vcproj new file mode 100644 index 0000000..5a62773 --- /dev/null +++ b/LumosQuad.vcproj @@ -0,0 +1,322 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/Makefile b/Makefile new file mode 100644 index 0000000..de3a969 --- /dev/null +++ b/Makefile @@ -0,0 +1,52 @@ +all: build/lumosquad + +build/lumosquad: build/lumosquad.o build/APSF.o build/CELL.o build/CG_SOLVER_SSE.o build/CG_SOLVER.o build/DAG.o build/EXR.o build/FFT.o build/QUAD_DBM_2D.o build/QUAD_POISSON.o build/ppm.o build/BLUE_NOISE.o build/RangeList.o build/RNG.o build/ScallopedSector.o build/WeightedDiscretePDF.o + g++ -g -o build/lumosquad build/*.o -lGL -lglut -lGLU -lfftw3 -lOpenEXR -lImath + +build/lumosquad.o: main.cpp + g++ -g -O3 -o build/lumosquad.o -c main.cpp -I/usr/include/Imath + +build/APSF.o: + g++ -g -O3 -o build/APSF.o -c APSF.cpp + +build/CELL.o: + g++ -g -O3 -o build/CELL.o -c CELL.cpp + +build/CG_SOLVER_SSE.o: + g++ -g -O3 -o build/CG_SOLVER_SSE.o -c CG_SOLVER_SSE.cpp + +build/CG_SOLVER.o: + g++ -g -O3 -o build/CG_SOLVER.o -c CG_SOLVER.cpp + +build/DAG.o: + g++ -g -O3 -o build/DAG.o -c DAG.cpp + +build/EXR.o: + g++ -g -O3 -o build/EXR.o -c EXR.cpp -I/usr/include/Imath + +build/FFT.o: + g++ -g -O3 -o build/FFT.o -c FFT.cpp + +build/QUAD_DBM_2D.o: + g++ -g -O3 -o build/QUAD_DBM_2D.o -c QUAD_DBM_2D.cpp + +build/QUAD_POISSON.o: + g++ -g -O3 -o build/QUAD_POISSON.o -c QUAD_POISSON.cpp + +build/ppm.o: + g++ -g -O3 -o build/ppm.o -c ppm/ppm.cpp + +build/BLUE_NOISE.o: + g++ -g -O3 -o build/BLUE_NOISE.o -c BlueNoise/BLUE_NOISE.cpp + +build/RangeList.o: + g++ -g -O3 -o build/RangeList.o -c BlueNoise/RangeList.cpp + +build/RNG.o: + g++ -g -O3 -o build/RNG.o -c BlueNoise/RNG.cpp + +build/ScallopedSector.o: + g++ -g -O3 -o build/ScallopedSector.o -c BlueNoise/ScallopedSector.cpp + +build/WeightedDiscretePDF.o: + g++ -g -O3 -o build/WeightedDiscretePDF.o -c BlueNoise/WeightedDiscretePDF.cpp diff --git a/QUAD_DBM_2D.cpp b/QUAD_DBM_2D.cpp new file mode 100644 index 0000000..036469c --- /dev/null +++ b/QUAD_DBM_2D.cpp @@ -0,0 +1,461 @@ +/////////////////////////////////////////////////////////////////////////////////// +// File : QUAD_DBM_2D.cpp +/////////////////////////////////////////////////////////////////////////////////// +// +// LumosQuad - A Lightning Generator +// Copyright 2007 +// The University of North Carolina at Chapel Hill +// +/////////////////////////////////////////////////////////////////////////////////// +// +// 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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// The University of North Carolina at Chapel Hill makes no representations +// about the suitability of this software for any purpose. It is provided +// "as is" without express or implied warranty. +// +// Permission to use, copy, modify and distribute this software and its +// documentation for educational, research and non-profit purposes, without +// fee, and without a written agreement is hereby granted, provided that the +// above copyright notice and the following three paragraphs appear in all +// copies. +// +// THE UNIVERSITY OF NORTH CAROLINA SPECIFICALLY DISCLAIM ANY WARRANTIES, +// INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND +// FITNESS FOR A PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN +// "AS IS" BASIS, AND THE UNIVERSITY OF NORTH CAROLINA HAS NO OBLIGATION TO +// PROVIDE MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS. +// +// Please send questions and comments about LumosQuad to kim@cs.unc.edu. +// +/////////////////////////////////////////////////////////////////////////////////// +// +// This program uses OpenEXR, which has the following restrictions: +// +// Copyright (c) 2002, Industrial Light & Magic, a division of Lucas +// Digital Ltd. LLC +// +// All rights reserved. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are +// met: +// * Redistributions of source code must retain the above copyright +// notice, this list of conditions and the following disclaimer. +// * Redistributions in binary form must reproduce the above +// copyright notice, this list of conditions and the following disclaimer +// in the documentation and/or other materials provided with the +// distribution. +// * Neither the name of Industrial Light & Magic nor the names of +// its contributors may be used to endorse or promote products derived +// from this software without specific prior written permission. +// + +#include "QUAD_DBM_2D.h" + +////////////////////////////////////////////////////////////////////// +// Construction/Destruction +////////////////////////////////////////////////////////////////////// + +QUAD_DBM_2D::QUAD_DBM_2D(int xRes, int yRes, int iterations) : + _xRes(xRes), + _yRes(yRes), + _bottomHit(0), + _iterations(iterations), + _quadPoisson(NULL), + _dag(NULL), + _skips(10), + _twister(123456) +{ + allocate(); + _dag = new DAG(_xRes, _yRes); + + // calculate dimensions + _dx = 1.0f / (float)_xRes; + _dy = 1.0f / (float)_yRes; + if (_dx < _dy) _dy = _dx; + else _dx = _dy; + + _maxRes = _xRes * _yRes; +} + +QUAD_DBM_2D::~QUAD_DBM_2D() +{ + deallocate(); +} + +void QUAD_DBM_2D::allocate() +{ + _quadPoisson = new QUAD_POISSON(_xRes, _yRes, _iterations); + _xRes = _yRes = _quadPoisson->maxRes(); +} + +void QUAD_DBM_2D::deallocate() +{ + if (_dag) delete _dag; + if (_quadPoisson) delete _quadPoisson; +} + +////////////////////////////////////////////////////////////////////// +// check neighbors for any candidate nodes +////////////////////////////////////////////////////////////////////// +void QUAD_DBM_2D::checkForCandidates(CELL* cell) +{ + int maxDepth = _quadPoisson->maxDepth(); + + CELL* north = cell->northNeighbor(); + if (north) { + if (north->depth == maxDepth) { + if (!north->candidate) { + _candidates.push_back(north); + north->candidate = true; + } + + CELL* northeast = north->eastNeighbor(); + if (northeast && !northeast->candidate) { + _candidates.push_back(northeast); + northeast->candidate = true; + } + CELL* northwest = north->westNeighbor(); + if (northwest && !northwest->candidate) { + _candidates.push_back(northwest); + northwest->candidate = true; + } + } + } + + CELL* east = cell->eastNeighbor(); + if (east && !east->candidate) { + _candidates.push_back(east); + east->candidate = true; + } + + CELL* south = cell->southNeighbor(); + if (south) { + if (!south->candidate) { + _candidates.push_back(south); + south->candidate = true; + } + + CELL* southeast = south->eastNeighbor(); + if (southeast && !southeast->candidate) { + _candidates.push_back(southeast); + southeast->candidate = true; + } + + CELL* southwest = south->westNeighbor(); + if (southwest && !southwest->candidate) { + _candidates.push_back(southwest); + southwest->candidate = true; + } + } + + CELL* west = cell->westNeighbor(); + if (west && !west->candidate) { + _candidates.push_back(west); + west->candidate = true; + } +} + +////////////////////////////////////////////////////////////////////// +// add particle to the aggregate +////////////////////////////////////////////////////////////////////// +bool QUAD_DBM_2D::addParticle() +{ + static float invSqrtTwo = 1.0f / sqrt(2.0f); + static int totalParticles = 0; + static int skipSolve = 0; + + // compute the potential + int iterations = 0; + if (!skipSolve) + iterations = _quadPoisson->solve(); + skipSolve++; + if (skipSolve == _skips) skipSolve = 0; + + // construct probability distribution + vector probabilities; + float totalPotential = 0.0f; + for (int x = 0; x < _candidates.size(); x++) + { + if (_candidates[x]->candidate) + { + probabilities.push_back(_candidates[x]->potential); + totalPotential += _candidates[x]->potential; + } + else + probabilities.push_back(0.0f); + } + + // get all the candidates + // if none are left, stop + if (_candidates.size() == 0) { + return false; + } + + // if there is not enough potential, go Brownian + int toAddIndex = 0; + if (totalPotential < 1e-8) + toAddIndex = _candidates.size() * _twister.getDoubleLR(); + // else follow DBM algorithm + else + { + // add a neighbor + float random = _twister.getDoubleLR(); + float invTotalPotential = 1.0f / totalPotential; + float potentialSeen = probabilities[0] * invTotalPotential; + while ((potentialSeen < random) && (toAddIndex < _candidates.size())) + { + toAddIndex++; + potentialSeen += probabilities[toAddIndex] * invTotalPotential; + } + } + _candidates[toAddIndex]->boundary = true; + _candidates[toAddIndex]->potential = 0.0f; + _candidates[toAddIndex]->state = NEGATIVE; + + CELL* neighbor = NULL; + CELL* added = _candidates[toAddIndex]; + CELL* north = added->northNeighbor(); + if (north) + { + if (north->state == NEGATIVE) + neighbor = north; + CELL* northeast = north->eastNeighbor(); + if (northeast && northeast->state == NEGATIVE) + neighbor = northeast; + CELL* northwest = north->westNeighbor(); + if (northwest && northwest->state == NEGATIVE) + neighbor = northwest; + } + CELL* east = added->eastNeighbor(); + if (east && east->state == NEGATIVE) + neighbor = east; + + CELL* south = added->southNeighbor(); + if (south) + { + if (south->state == NEGATIVE) + neighbor = south; + CELL* southeast = south->eastNeighbor(); + if (southeast && southeast->state == NEGATIVE) + neighbor = southeast; + CELL* southwest = south->westNeighbor(); + if (southwest && southwest->state == NEGATIVE) + neighbor = southwest; + } + CELL* west = added->westNeighbor(); + if (west && west->state == NEGATIVE) + neighbor = west; + + // insert it as a node for bookkeeping + _quadPoisson->insert(added->center[0], added->center[1]); + checkForCandidates(added); + + // insert into the dag + int newIndex = (int)(added->center[0] * _xRes) + + (int)(added->center[1] * _yRes) * _xRes; + int neighborIndex = (int)(neighbor->center[0] * _xRes) + + (int)(neighbor->center[1] * _yRes) * _xRes; + _dag->addSegment(newIndex, neighborIndex); + + totalParticles++; + if (!(totalParticles % 200)) + cout << " " << totalParticles; + + hitGround(added); + + return true; +} + +////////////////////////////////////////////////////////////////////// +// hit ground yet? +////////////////////////////////////////////////////////////////////// +bool QUAD_DBM_2D::hitGround(CELL* cell) +{ + if (_bottomHit) + return true; + + if (!cell) + return false; + + bool hit = false; + if (cell->northNeighbor()) + { + CELL* north = cell->northNeighbor(); + if (north->state == POSITIVE) + hit = true; + if (north->eastNeighbor()->state == POSITIVE) + hit = true; + if (north->westNeighbor()->state == POSITIVE) + hit = true; + } + if (cell->eastNeighbor()) + if (cell->eastNeighbor()->state == POSITIVE) + hit = true; + if (cell->southNeighbor()) + { + CELL* south = cell->southNeighbor(); + if (south->state == POSITIVE) + hit = true; + if (south->eastNeighbor()->state == POSITIVE) + hit = true; + if (south->westNeighbor()->state == POSITIVE) + hit = true; + } + if (cell->westNeighbor()) + if (cell->westNeighbor()->state == POSITIVE) + hit = true; + + if (hit) + { + _bottomHit = (int)(cell->center[0] * _xRes) + + (int)(cell->center[1] * _yRes) * _xRes; + _dag->buildLeader(_bottomHit); + return true; + } + return false; +} + +//////////////////////////////////////////////////////////////////// +// drawing functions +//////////////////////////////////////////////////////////////////// +void QUAD_DBM_2D::draw() { + glPushMatrix(); + glTranslatef(-0.5, -0.5, 0); + + list leaves; + list::iterator cellIterator = leaves.begin(); + _quadPoisson->getAllLeaves(leaves); + for (cellIterator = leaves.begin(); cellIterator != leaves.end(); cellIterator++) + { + float color = (*cellIterator)->potential; + + if ((*cellIterator)->boundary) { + if (color <= 0.0f) + _quadPoisson->drawCell(*cellIterator, 0,0,0); + else + _quadPoisson->drawCell(*cellIterator, 0,0,color); + } + else + _quadPoisson->drawCell(*cellIterator, color, 0,0); + } + _quadPoisson->draw(NULL); + + glPopMatrix(); +} + +//////////////////////////////////////////////////////////////////// +// read in attractors from an image +//////////////////////////////////////////////////////////////////// +bool QUAD_DBM_2D::readImage(unsigned char* initial, + unsigned char* attractors, + unsigned char* repulsors, + unsigned char* terminators, + int xRes, int yRes) +{ + _dag->inputWidth() = xRes; + _dag->inputHeight() = yRes; + + bool initialFound = false; + bool terminateFound = false; + + int index = 0; + for (int y = 0; y < yRes; y++) + for (int x = 0; x < xRes; x++, index++) + { + // insert initial condition + if (initial[index]) + { + // insert something + CELL* negative = _quadPoisson->insert(x, y); + negative->boundary = true; + negative->potential = 0.0f; + negative->state = NEGATIVE; + negative->candidate = true; + + checkForCandidates(negative); + + initialFound = true; + } + + // insert attractors + if (attractors[index]) + { + // insert something + CELL* positive = _quadPoisson->insert(x, y); + positive->boundary = true; + positive->potential = 1.0f; + positive->state = ATTRACTOR; + positive->candidate = true; + } + + // insert repulsors + if (repulsors[index]) + { + // only insert the repulsor if it is the edge of a repulsor + bool edge = false; + + if (x != 0) + { + if (!repulsors[index - 1]) edge = true; + if (y != 0 && !repulsors[index - xRes - 1]) edge = true; + if (y != yRes - 1 && !repulsors[index + xRes - 1]) edge = true; + } + if (x != _xRes - 1) + { + if (!repulsors[index + 1]) edge = true; + if (y != 0 && !repulsors[index - xRes + 1]) edge = true; + if (y != yRes - 1 && !repulsors[index + xRes + 1]) edge = true; + } + if (y != 0 && !repulsors[index - xRes]) edge = true; + if (y != yRes - 1 && !repulsors[index + xRes]) edge = true; + + if (edge) + { + // insert something + CELL* negative = _quadPoisson->insert(x, y); + negative->boundary = true; + negative->potential = 0.0f; + negative->state = REPULSOR; + negative->candidate = true; + } + } + + // insert terminators + if (terminators[index]) + { + // insert something + CELL* positive = _quadPoisson->insert(x, y); + positive->boundary = true; + positive->potential = 1.0f; + positive->state = POSITIVE; + positive->candidate = true; + + terminateFound = true; + } + } + + if (!initialFound) { + cout << " The lightning does not start anywhere! " << endl; + return false; + } + if (!terminateFound) { + cout << " The lightning does not end anywhere! " << endl; + return false; + } + + return true; +} diff --git a/QUAD_DBM_2D.h b/QUAD_DBM_2D.h new file mode 100644 index 0000000..bdeba3f --- /dev/null +++ b/QUAD_DBM_2D.h @@ -0,0 +1,196 @@ +/////////////////////////////////////////////////////////////////////////////////// +// File : QUAD_DBM_2D.h +/////////////////////////////////////////////////////////////////////////////////// +// +// LumosQuad - A Lightning Generator +// Copyright 2007 +// The University of North Carolina at Chapel Hill +// +/////////////////////////////////////////////////////////////////////////////////// +// +// 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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// The University of North Carolina at Chapel Hill makes no representations +// about the suitability of this software for any purpose. It is provided +// "as is" without express or implied warranty. +// +// Permission to use, copy, modify and distribute this software and its +// documentation for educational, research and non-profit purposes, without +// fee, and without a written agreement is hereby granted, provided that the +// above copyright notice and the following three paragraphs appear in all +// copies. +// +// THE UNIVERSITY OF NORTH CAROLINA SPECIFICALLY DISCLAIM ANY WARRANTIES, +// INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND +// FITNESS FOR A PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN +// "AS IS" BASIS, AND THE UNIVERSITY OF NORTH CAROLINA HAS NO OBLIGATION TO +// PROVIDE MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS. +// +// Please send questions and comments about LumosQuad to kim@cs.unc.edu. +// +/////////////////////////////////////////////////////////////////////////////////// +// +// This program uses OpenEXR, which has the following restrictions: +// +// Copyright (c) 2002, Industrial Light & Magic, a division of Lucas +// Digital Ltd. LLC +// +// All rights reserved. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are +// met: +// * Redistributions of source code must retain the above copyright +// notice, this list of conditions and the following disclaimer. +// * Redistributions in binary form must reproduce the above +// copyright notice, this list of conditions and the following disclaimer +// in the documentation and/or other materials provided with the +// distribution. +// * Neither the name of Industrial Light & Magic nor the names of +// its contributors may be used to endorse or promote products derived +// from this software without specific prior written permission. +// + +#ifndef QUAD_DBM_2D_H +#define QUAD_DBM_2D_H + +#include +#include +#include "DAG.h" +#include "QUAD_POISSON.h" + +//////////////////////////////////////////////////////////////////// +/// \brief Quadtree DBM solver. This is the highest level class. +//////////////////////////////////////////////////////////////////// +class QUAD_DBM_2D +{ +public: + /// \brief DBM constructor + /// + /// \param xRes maximum x resolution + /// \param yRes maximum y resolution + /// \param iterations maximum conjugate gradient iterations + QUAD_DBM_2D(int xRes = 128, int yRes = 128, int iterations = 10); + + //! destructor + virtual ~QUAD_DBM_2D(); + + //! add to aggregate + bool addParticle(); + + /// \brief Hit ground yet? + /// \return returns true if a terminator as already been hit + bool hitGround(CELL* cell = NULL); + + //! draw the quadtree cells to OpenGL + void draw(); + + //! draw the DAG to OpenGL + void drawSegments() { + glLineWidth(1.0f); + glPushMatrix(); + glTranslatef(-0.5f, -0.5f, 0.0f); + _dag->draw(); + glPopMatrix(); + }; + + //////////////////////////////////////////////////////////////// + // file IO + //////////////////////////////////////////////////////////////// + + //! write everything to a file + void writeFields(const char* filename); + + //! read everything from a file + void readFields(const char* filename); + + /// \brief read in control parameters from an input file + /// + /// \param initial initial pixels of lightning + /// \param attractors pixels that attract the lightning + /// \param repulsors pixels that repulse the lightning + /// \param terminators pixels that halt the simulation if hit + /// \param xRes x resolution of the image + /// \param yRes y resolution of the image + /// + /// \return Returns false if it finds something wrong with the images + bool readImage(unsigned char* initial, + unsigned char* attractors, + unsigned char* repulsors, + unsigned char* terminators, + int xRes, int yRes); + + //! read in a new DAG + void readDAG(const char* filename) { _dag->read(filename); }; + + //! write out the current DAG + void writeDAG(const char* filename) { _dag->write(filename); }; + + /// \brief render to a software-only buffer + /// + /// \param scale a (scale * xRes) x (scale * yRes) image is rendered + float*& renderOffscreen(int scale = 1) { return _dag->drawOffscreen(scale); }; + + //! access the DBM x resolution + int xRes() { return _xRes; }; + //! access the DBM y resolution + int yRes() { return _yRes; }; + //! access the DAG x resolution + int xDagRes() { return _dag->xRes(); }; + //! access the DAG y resolution + int yDagRes() { return _dag->yRes(); }; + //! access the x resolution of the input image + int inputWidth() { return _dag->inputWidth(); }; + //! access the y resolution of the input image + int inputHeight() { return _dag->inputHeight(); }; + +private: + void allocate(); + void deallocate(); + + //////////////////////////////////////////////////////////////////// + // dielectric breakdown model components + //////////////////////////////////////////////////////////////////// + + // field dimensions + int _xRes; + int _yRes; + int _maxRes; + float _dx; + float _dy; + int _iterations; + + // which cell did it hit bottom with? + int _bottomHit; + + DAG* _dag; + + QUAD_POISSON* _quadPoisson; + + // current candidate list + vector _candidates; + + // check if any of the neighbors of cell should be added to the + // candidate list + void checkForCandidates(CELL* cell); + + // number of particles to add before doing another Poisson solve + int _skips; + + // Mersenne Twister + RNG _twister; +}; + +#endif diff --git a/QUAD_POISSON.cpp b/QUAD_POISSON.cpp new file mode 100644 index 0000000..d8d12b5 --- /dev/null +++ b/QUAD_POISSON.cpp @@ -0,0 +1,620 @@ +/////////////////////////////////////////////////////////////////////////////////// +// File : QUAD_POISSON.cpp +/////////////////////////////////////////////////////////////////////////////////// +// +// LumosQuad - A Lightning Generator +// Copyright 2007 +// The University of North Carolina at Chapel Hill +// +/////////////////////////////////////////////////////////////////////////////////// +// +// 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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// The University of North Carolina at Chapel Hill makes no representations +// about the suitability of this software for any purpose. It is provided +// "as is" without express or implied warranty. +// +// Permission to use, copy, modify and distribute this software and its +// documentation for educational, research and non-profit purposes, without +// fee, and without a written agreement is hereby granted, provided that the +// above copyright notice and the following three paragraphs appear in all +// copies. +// +// THE UNIVERSITY OF NORTH CAROLINA SPECIFICALLY DISCLAIM ANY WARRANTIES, +// INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND +// FITNESS FOR A PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN +// "AS IS" BASIS, AND THE UNIVERSITY OF NORTH CAROLINA HAS NO OBLIGATION TO +// PROVIDE MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS. +// +// Please send questions and comments about LumosQuad to kim@cs.unc.edu. +// +/////////////////////////////////////////////////////////////////////////////////// +// +// This program uses OpenEXR, which has the following restrictions: +// +// Copyright (c) 2002, Industrial Light & Magic, a division of Lucas +// Digital Ltd. LLC +// +// All rights reserved. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are +// met: +// * Redistributions of source code must retain the above copyright +// notice, this list of conditions and the following disclaimer. +// * Redistributions in binary form must reproduce the above +// copyright notice, this list of conditions and the following disclaimer +// in the documentation and/or other materials provided with the +// distribution. +// * Neither the name of Industrial Light & Magic nor the names of +// its contributors may be used to endorse or promote products derived +// from this software without specific prior written permission. +// + +#include "QUAD_POISSON.h" + +////////////////////////////////////////////////////////////////////// +// Construction/Destruction +////////////////////////////////////////////////////////////////////// + +QUAD_POISSON::QUAD_POISSON(int xRes, int yRes, int iterations) : + _root(new CELL(1.0f, 1.0f, 0.0f, 0.0f)), + _noise() +{ + _root->refine(); + + // figure out the max depth needed + float xMax = log((float)xRes) / log(2.0f); + float yMax = log((float)yRes) / log(2.0f); + + float max = (xMax > yMax) ? xMax : yMax; + if (max - floor(max) > 1e-7) + max = max + 1; + max = floor(max); + + _maxRes = pow(2.0f, (float)max); + _maxDepth = max; + + // create the blue noise + _noiseFunc = new BLUE_NOISE(5.0f / (float)_maxRes); + _noise = new bool[_maxRes * _maxRes]; + _noiseFunc->complete(); + _noiseFunc->maximize(); + _noiseFunc->writeToBool(_noise, _maxRes); + + _solver = new CG_SOLVER(_maxDepth, iterations); +} + +QUAD_POISSON::~QUAD_POISSON() +{ + deleteGhosts(); + delete _root; + delete _solver; + delete _noiseFunc; + delete[] _noise; +} + +////////////////////////////////////////////////////////////////////// +// draw boundaries to OGL +////////////////////////////////////////////////////////////////////// +void QUAD_POISSON::draw(CELL* cell) +{ + // see if it's the root + if (cell == NULL) { + draw(_root); + return; + } + + // draw the current cell + glColor4f(1,1,1,0.1); + + glBegin(GL_LINE_STRIP); + glVertex2f(cell->bounds[1], 1.0f - cell->bounds[0]); + glVertex2f(cell->bounds[1], 1.0f - cell->bounds[2]); + glVertex2f(cell->bounds[3], 1.0f - cell->bounds[2]); + glVertex2f(cell->bounds[3], 1.0f - cell->bounds[0]); + glVertex2f(cell->bounds[1], 1.0f - cell->bounds[0]); + glEnd(); + + // draw the children + for (int x = 0; x < 4; x++) + if (cell->children[x] != NULL) + draw(cell->children[x]); +} + +////////////////////////////////////////////////////////////////////// +// draw the given cell +////////////////////////////////////////////////////////////////////// +void QUAD_POISSON::drawCell(CELL* cell, float r, float g, float b) +{ + // draw the current cell + glColor4f(r,g,b,1.0f); + glBegin(GL_QUADS); + glVertex2f(cell->bounds[1], 1.0f - cell->bounds[0]); + glVertex2f(cell->bounds[1], 1.0f - cell->bounds[2]); + glVertex2f(cell->bounds[3], 1.0f - cell->bounds[2]); + glVertex2f(cell->bounds[3], 1.0f - cell->bounds[0]); + glEnd(); +} + +////////////////////////////////////////////////////////////////////// +// subdivide quadtree to max level for (xPos, yPos) +// returns a pointer to the cell that was created +////////////////////////////////////////////////////////////////////// +CELL* QUAD_POISSON::insert(float xPos, float yPos) +{ + int currentDepth = 0; + CELL* currentCell = _root; + bool existed = true; + + while (currentDepth < _maxDepth) { + // find quadrant of current point + float diff[2]; + diff[0] = xPos - currentCell->center[0]; + diff[1] = yPos - currentCell->center[1]; + int quadrant = 1; + if (diff[0] > 0.0f) { + if (diff[1] < 0.0f) + quadrant = 2; + } + else if (diff[1] < 0.0f) + quadrant = 3; + else + quadrant = 0; + + // check if it exists + if (currentCell->children[quadrant] == NULL) { + existed = false; + currentCell->refine(); + } + + // recurse to next level + currentCell = currentCell->children[quadrant]; + + // increment depth + currentDepth++; + } + // if we had to subdivide to get the cell, add them to the list + if (!existed) + for (int i = 0; i < 4; i++) + { + _smallestLeaves.push_back(currentCell->parent->children[i]); + setNoise(currentCell->parent->children[i]); + } + + /////////////////////////////////////////////////////////////////// + // force orthogonal neighbors to be same depth + // I have commented the first block, the rest follow the same flow + /////////////////////////////////////////////////////////////////// + + // see if neighbor exists + CELL* north = currentCell->northNeighbor(); + if (north && north->depth != _maxDepth) { + // while the neighbor needs to be refined + while (north->depth != _maxDepth) { + + // refine it + north->refine(); + + // set to the newly refined neighbor + north = currentCell->northNeighbor(); + } + // add newly created nodes to the list + for (int i = 0; i < 4; i++) + { + _smallestLeaves.push_back(north->parent->children[i]); + setNoise(north->parent->children[i]); + } + } + CELL* south = currentCell->southNeighbor(); + if (south && south->depth != _maxDepth) { + while (south->depth != _maxDepth) { + south->refine(); + south = currentCell->southNeighbor(); + } + for (int i = 0; i < 4; i++) + { + _smallestLeaves.push_back(south->parent->children[i]); + setNoise(south->parent->children[i]); + } + } + CELL* west = currentCell->westNeighbor(); + if (west && west->depth != _maxDepth) { + while (west->depth != _maxDepth) { + west->refine(); + west = currentCell->westNeighbor(); + } + for (int i = 0; i < 4; i++) + { + _smallestLeaves.push_back(west->parent->children[i]); + setNoise(west->parent->children[i]); + } + } + CELL* east = currentCell->eastNeighbor(); + if (east && east->depth != _maxDepth) { + while (east->depth != _maxDepth) { + east->refine(); + east = currentCell->eastNeighbor(); + } + for (int i = 0; i < 4; i++) + { + _smallestLeaves.push_back(east->parent->children[i]); + setNoise(east->parent->children[i]); + } + } + + /////////////////////////////////////////////////////////////////// + // force diagonal neighbors to be same depth + // The same flow follows as above, except that it makes sure that + // the 'north' and 'south' neighbors already exist + /////////////////////////////////////////////////////////////////// + + if (north) { + CELL* northwest = north->westNeighbor(); + if (northwest && northwest->depth != _maxDepth) { + while (northwest->depth != _maxDepth) { + northwest->refine(); + northwest = northwest->children[2]; + } + for (int i = 0; i < 4; i++) + { + _smallestLeaves.push_back(northwest->parent->children[i]); + setNoise(northwest->parent->children[i]); + } + } + CELL* northeast = north->eastNeighbor(); + if (northeast && northeast->depth != _maxDepth) { + while (northeast->depth != _maxDepth) { + northeast->refine(); + northeast= northeast->children[3]; + } + for (int i = 0; i < 4; i++) + { + _smallestLeaves.push_back(northeast->parent->children[i]); + setNoise(northeast->parent->children[i]); + } + } + } + if (south) { + CELL* southwest = south->westNeighbor(); + if (southwest && southwest->depth != _maxDepth) { + while (southwest->depth != _maxDepth) { + southwest->refine(); + southwest = southwest->children[1]; + } + for (int i = 0; i < 4; i++) + { + _smallestLeaves.push_back(southwest->parent->children[i]); + setNoise(southwest->parent->children[i]); + } + } + CELL* southeast = south->eastNeighbor(); + if (southeast && southeast->depth != _maxDepth) { + while (southeast->depth != _maxDepth) { + southeast->refine(); + southeast= southeast->children[0]; + } + for (int i = 0; i < 4; i++) + { + _smallestLeaves.push_back(southeast->parent->children[i]); + setNoise(southeast->parent->children[i]); + } + } + } + + return currentCell; +} + +////////////////////////////////////////////////////////////////////// +// check if a cell hits a noise node +////////////////////////////////////////////////////////////////////// +void QUAD_POISSON::setNoise(CELL* cell) +{ + if (!(cell->state == EMPTY)) + return; + + int x = cell->center[0] * _maxRes; + int y = cell->center[1] * _maxRes; + + if (_noise[x + y * _maxRes]) + { + cell->boundary = true; + cell->state = ATTRACTOR; + cell->potential = 0.5f; + cell->candidate = true; + } +} + +////////////////////////////////////////////////////////////////////// +// insert all leaves into a list +////////////////////////////////////////////////////////////////////// +void QUAD_POISSON::getAllLeaves(list& leaves, CELL* currentCell) +{ + // if we're at the root + if (currentCell == NULL) + { + getAllLeaves(leaves, _root); + return; + } + + // if we're at a leaf, add it to the list + if (currentCell->children[0] == NULL) + { + leaves.push_back(currentCell); + return; + } + + // if children exist, call recursively + for (int x = 0; x < 4; x++) + getAllLeaves(leaves, currentCell->children[x]); +} + +////////////////////////////////////////////////////////////////////// +// insert all leaves not on the boundary into a list +////////////////////////////////////////////////////////////////////// +void QUAD_POISSON::getEmptyLeaves(list& leaves, CELL* currentCell) +{ + // if we're at the root + if (currentCell == NULL) { + getEmptyLeaves(leaves, _root); + return; + } + + // if we're at a leaf, check if it's a boundary and then + // add it to the list + if (currentCell->children[0] == NULL) { + if (!(currentCell->boundary)) + leaves.push_back(currentCell); + return; + } + + // if children exist, call recursively + for (int x = 0; x < 4; x++) + getEmptyLeaves(leaves, currentCell->children[x]); +} + +////////////////////////////////////////////////////////////////////// +// balance the current tree +////////////////////////////////////////////////////////////////////// +void QUAD_POISSON::balance() +{ + // collect all the leaf nodes + list leaves; + getAllLeaves(leaves); + + // while the list is not empty + list::iterator cellIterator = leaves.begin(); + for (cellIterator = leaves.begin(); cellIterator != leaves.end(); cellIterator++) { + CELL* currentCell = *cellIterator; + + // if a north neighbor exists + CELL* north = currentCell->northNeighbor(); + if (north != NULL) + // while the neighbor is not balanced + while (north->depth < currentCell->depth - 1) { + // refine it + north->refine(); + + // add the newly refined nodes to the list of + // those to be checked + for (int x = 0; x < 4; x++) + leaves.push_back(north->children[x]); + + // set the cell to the newly created one + north = currentCell->northNeighbor(); + } + + // the rest of the blocks flow the same as above + CELL* south = currentCell->southNeighbor(); + if (south!= NULL) + while (south->depth < currentCell->depth - 1) { + south->refine(); + for (int x = 0; x < 4; x++) + leaves.push_back(south->children[x]); + south = currentCell->southNeighbor(); + } + + CELL* west = currentCell->westNeighbor(); + if (west != NULL) + while (west->depth < currentCell->depth - 1) { + west->refine(); + for (int x = 0; x < 4; x++) + leaves.push_back(west->children[x]); + west = currentCell->westNeighbor(); + } + + CELL* east = currentCell->eastNeighbor(); + if (east != NULL) + while (east->depth < currentCell->depth - 1) { + east->refine(); + for (int x = 0; x < 4; x++) + leaves.push_back(east->children[x]); + east = currentCell->eastNeighbor(); + } + } +} + +////////////////////////////////////////////////////////////////////// +// build the neighbor lists of the current quadtree +////////////////////////////////////////////////////////////////////// +void QUAD_POISSON::buildNeighbors() +{ + balance(); + + // collect all the leaf nodes + list leaves; + getAllLeaves(leaves); + + list::iterator cellIterator = leaves.begin(); + for (cellIterator = leaves.begin(); cellIterator != leaves.end(); cellIterator++) + { + CELL* currentCell = *cellIterator; + + // build north neighbors + CELL* north = currentCell->northNeighbor(); + if (north != NULL) { + if (north->children[0] == NULL) { + currentCell->neighbors[0] = north; + currentCell->neighbors[1] = NULL; + } + else { + currentCell->neighbors[0] = north->children[3]; + currentCell->neighbors[1] = north->children[2]; + } + } + // else build a ghost cell + else + currentCell->neighbors[0] = new CELL(currentCell->depth); + + // build east neighbors + CELL* east = currentCell->eastNeighbor(); + if (east != NULL) { + if (east->children[0] == NULL) { + currentCell->neighbors[2] = east; + currentCell->neighbors[3] = NULL; + } + else { + currentCell->neighbors[2] = east->children[0]; + currentCell->neighbors[3] = east->children[3]; + } + } + // else build a ghost cell + else + currentCell->neighbors[2] = new CELL(currentCell->depth); + + // build south neighbors + CELL* south = currentCell->southNeighbor(); + if (south != NULL) { + if (south->children[0] == NULL) { + currentCell->neighbors[4] = south; + currentCell->neighbors[5] = NULL; + } + else { + currentCell->neighbors[4] = south->children[1]; + currentCell->neighbors[5] = south->children[0]; + } + } + // else build a ghost cell + else + currentCell->neighbors[4] = new CELL(currentCell->depth); + + // build west neighbors + CELL* west = currentCell->westNeighbor(); + if (west != NULL) { + if (west->children[0] == NULL) { + currentCell->neighbors[6] = west; + currentCell->neighbors[7] = NULL; + } + else { + currentCell->neighbors[6] = west->children[2]; + currentCell->neighbors[7] = west->children[1]; + } + } + // else build a ghost cell + else + currentCell->neighbors[6] = new CELL(currentCell->depth); + } +} + +////////////////////////////////////////////////////////////////////// +// delete ghost cells +////////////////////////////////////////////////////////////////////// +void QUAD_POISSON::deleteGhosts(CELL* currentCell) +{ + // if at the base, call on the root + if (currentCell == NULL) + { + deleteGhosts(_root); + return; + } + + // if there are children, delete those too + if (currentCell->children[0]) { + // call recursively + for (int x = 0; x < 4; x++) + deleteGhosts(currentCell->children[x]); + return; + } + + // check the neighbors for stuff to delete + for (int x = 0; x < 8; x++) { + // if the neighbor exists + if (currentCell->neighbors[x]) + // and if it is a ghost cell, delete it + if (currentCell->neighbors[x]->parent == NULL) + delete currentCell->neighbors[x]; + } +} + +////////////////////////////////////////////////////////////////////// +// solve the Poisson problem +////////////////////////////////////////////////////////////////////// +int QUAD_POISSON::solve() { + // maintain the quadtree + balance(); + buildNeighbors(); + + // retrieve leaves at the lowest level + _emptyLeaves.clear(); + getEmptyLeaves(_emptyLeaves); + + static bool firstSolve = true; + static int iterations = _solver->iterations(); + + // do a full precision solve the first time + if (firstSolve) + { + iterations = _solver->iterations(); + _solver->iterations() = 10000; + firstSolve = false; + } + else + _solver->iterations() = iterations; + + // return the number of iterations + return _solver->solve(_emptyLeaves); +}; + + +////////////////////////////////////////////////////////////////////// +// get the leafnode that corresponds to the coordinate +////////////////////////////////////////////////////////////////////// +CELL* QUAD_POISSON::getLeaf(float xPos, float yPos) +{ + CELL* currentCell = _root; + + while (currentCell->children[0] != NULL) + { + // find quadrant of current point + float diff[2]; + diff[0] = xPos - currentCell->center[0]; + diff[1] = yPos - currentCell->center[1]; + int quadrant = 1; + if (diff[0] > 0.0f) + { + if (diff[1] < 0.0f) + quadrant = 2; + } + else if (diff[1] < 0.0f) + quadrant = 3; + else + quadrant = 0; + + // check if it exists + if (currentCell->children[quadrant] != NULL) + currentCell = currentCell->children[quadrant]; + } + return currentCell; +} diff --git a/QUAD_POISSON.h b/QUAD_POISSON.h new file mode 100644 index 0000000..6d129ab --- /dev/null +++ b/QUAD_POISSON.h @@ -0,0 +1,190 @@ +/////////////////////////////////////////////////////////////////////////////////// +// File : QUAD_POISSON.h +/////////////////////////////////////////////////////////////////////////////////// +// +// LumosQuad - A Lightning Generator +// Copyright 2007 +// The University of North Carolina at Chapel Hill +// +/////////////////////////////////////////////////////////////////////////////////// +// +// 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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// The University of North Carolina at Chapel Hill makes no representations +// about the suitability of this software for any purpose. It is provided +// "as is" without express or implied warranty. +// +// Permission to use, copy, modify and distribute this software and its +// documentation for educational, research and non-profit purposes, without +// fee, and without a written agreement is hereby granted, provided that the +// above copyright notice and the following three paragraphs appear in all +// copies. +// +// THE UNIVERSITY OF NORTH CAROLINA SPECIFICALLY DISCLAIM ANY WARRANTIES, +// INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND +// FITNESS FOR A PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN +// "AS IS" BASIS, AND THE UNIVERSITY OF NORTH CAROLINA HAS NO OBLIGATION TO +// PROVIDE MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS. +// +// Please send questions and comments about LumosQuad to kim@cs.unc.edu. +// +/////////////////////////////////////////////////////////////////////////////////// +// +// This program uses OpenEXR, which has the following restrictions: +// +// Copyright (c) 2002, Industrial Light & Magic, a division of Lucas +// Digital Ltd. LLC +// +// All rights reserved. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are +// met: +// * Redistributions of source code must retain the above copyright +// notice, this list of conditions and the following disclaimer. +// * Redistributions in binary form must reproduce the above +// copyright notice, this list of conditions and the following disclaimer +// in the documentation and/or other materials provided with the +// distribution. +// * Neither the name of Industrial Light & Magic nor the names of +// its contributors may be used to endorse or promote products derived +// from this software without specific prior written permission. +// + +#ifndef QUAD_POISSON_H +#define QUAD_POISSON_H + +#include +#include +#include "CELL.h" +#include +#include "CG_SOLVER.h" +#include "CG_SOLVER_SSE.h" +#include "BlueNoise/BLUE_NOISE.h" + +#include + +using namespace std; + +////////////////////////////////////////////////////////////////////// +/// \brief Quadtree Poisson solver +////////////////////////////////////////////////////////////////////// +class QUAD_POISSON +{ +public: + /// \brief quadtree constructor + /// + /// \param xRes maximum x resolution + /// \param yRes maximum y resolution + /// \param iterations maximum conjugate gradient iterations + QUAD_POISSON(int xRes, + int yRes, + int iterations = 10); + + //! destructor + virtual ~QUAD_POISSON(); + + /// \brief OpenGL drawing function + /// + /// \param cell internally used param, should always be NULL externally + void draw(CELL* cell = NULL); + + /// \brief Draw a single OpenGL cell + /// + /// \param cell quadtree cell to draw + /// \param r red intensity to draw + /// \param g green intensity to draw + /// \param b blue intensity to draw + void drawCell(CELL* cell, + float r = 1.0f, + float g = 0.0f, + float b = 0.0f); + + //! Solve the Poisson problem + int solve(); + + /// \brief insert point at maximum subdivision level + /// + /// \param xPos x position to insert at + /// \param yPos y position to insert at + CELL* insert(float xPos, float yPos); + + /// \brief insert point at maximum subdivision level + /// + /// \param xPos grid cell x index to insert + /// \param yPos grid cell y index to insert + CELL* insert(int xPos, int yPos) { + return insert((float)xPos / _maxRes, (float)yPos / _maxRes); + }; + + /// \brief get all the leaf nodes + /// \return Leaves are returned in the 'leaves' param + void getAllLeaves(list& leaves, CELL* currentCell = NULL); + + /// \brief get all the leaf nodes at finest subdivision level + /// \return Leaves are returned in the 'leaves' param + list& getSmallestLeaves() { return _smallestLeaves; }; + + //! maximum resolution accessor + int& maxRes() { return _maxRes; } + + //! maximum depth accessor + int& maxDepth() { return _maxDepth; }; + + //! get leaf at coordinate (x,y) + CELL* getLeaf(float xPos, float yPos); + +private: + //! root of the quadtree + CELL* _root; + + //! maximum resolution of quadtree + int _maxRes; + + //! maxmimum depth of quadtree + int _maxDepth; + + //! dependant leaves + list _emptyLeaves; + + //! smallest leaves + list _smallestLeaves; + + //! current Poisson solver + CG_SOLVER* _solver; + + //! balance quadtree + void balance(); + + //! get the leaf nodes not on the boundary + void getEmptyLeaves(list& leaves, CELL* currentCell = NULL); + + //! build the neighbor lists of the cells + void buildNeighbors(); + + //! delete ghost cells + void deleteGhosts(CELL* currentCell = NULL); + + //! Blue noise function + BLUE_NOISE* _noiseFunc; + + //! Blue noise sample locations + bool* _noise; + + //! check if a cell hits a noise node + void setNoise(CELL* cell); +}; + +#endif diff --git a/README.md b/README.md new file mode 100644 index 0000000..240e0a3 --- /dev/null +++ b/README.md @@ -0,0 +1,10 @@ +# LumosQuad + +[Original code source](https://gamma.cs.unc.edu/FAST_LIGHTNING/lumos.html) slightly modified to work on Debian. + +Packages to install: `freeglut3-dev libfftw3-dev libopenexr-dev libopengl-dev` + + make + ./build/lumosquad examples/spine.ppm o.exr + +EXR files are images that can be opened by Gimp. diff --git a/examples/hint.ppm b/examples/hint.ppm new file mode 100644 index 0000000..af55298 Binary files /dev/null and b/examples/hint.ppm differ diff --git a/examples/nopath.ppm b/examples/nopath.ppm new file mode 100644 index 0000000..920df52 Binary files /dev/null and b/examples/nopath.ppm differ diff --git a/examples/spine.lightning b/examples/spine.lightning new file mode 100644 index 0000000..80130ec Binary files /dev/null and b/examples/spine.lightning differ diff --git a/examples/spine.ppm b/examples/spine.ppm new file mode 100644 index 0000000..1c1b749 Binary files /dev/null and b/examples/spine.ppm differ diff --git a/examples/spineless.ppm b/examples/spineless.ppm new file mode 100644 index 0000000..7871fac Binary files /dev/null and b/examples/spineless.ppm differ diff --git a/examples/y-second.ppm b/examples/y-second.ppm new file mode 100644 index 0000000..70c8515 Binary files /dev/null and b/examples/y-second.ppm differ diff --git a/examples/y.ppm b/examples/y.ppm new file mode 100644 index 0000000..238c013 Binary files /dev/null and b/examples/y.ppm differ diff --git a/main.cpp b/main.cpp new file mode 100644 index 0000000..0651087 --- /dev/null +++ b/main.cpp @@ -0,0 +1,408 @@ +/////////////////////////////////////////////////////////////////////////////////// +// File : main.cpp +/////////////////////////////////////////////////////////////////////////////////// +// +// LumosQuad - A Lightning Generator +// Copyright 2007 +// The University of North Carolina at Chapel Hill +// +/////////////////////////////////////////////////////////////////////////////////// +// +// 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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// The University of North Carolina at Chapel Hill makes no representations +// about the suitability of this software for any purpose. It is provided +// "as is" without express or implied warranty. +// +// Permission to use, copy, modify and distribute this software and its +// documentation for educational, research and non-profit purposes, without +// fee, and without a written agreement is hereby granted, provided that the +// above copyright notice and the following three paragraphs appear in all +// copies. +// +// THE UNIVERSITY OF NORTH CAROLINA SPECIFICALLY DISCLAIM ANY WARRANTIES, +// INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND +// FITNESS FOR A PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN +// "AS IS" BASIS, AND THE UNIVERSITY OF NORTH CAROLINA HAS NO OBLIGATION TO +// PROVIDE MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS. +// +// Please send questions and comments about LumosQuad to kim@cs.unc.edu. +// +/////////////////////////////////////////////////////////////////////////////////// +// +// This program uses OpenEXR, which has the following restrictions: +// +// Copyright (c) 2002, Industrial Light & Magic, a division of Lucas +// Digital Ltd. LLC +// +// All rights reserved. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are +// met: +// * Redistributions of source code must retain the above copyright +// notice, this list of conditions and the following disclaimer. +// * Redistributions in binary form must reproduce the above +// copyright notice, this list of conditions and the following disclaimer +// in the documentation and/or other materials provided with the +// distribution. +// * Neither the name of Industrial Light & Magic nor the names of +// its contributors may be used to endorse or promote products derived +// from this software without specific prior written permission. +// +/////////////////////////////////////////////////////////////////////////////////// +/// +/// \mainpage Fast Animation of Lightning Using An Adaptive Mesh +/// \section Introduction +/// +/// This project is an implementation of the paper +/// Fast Animation of Lightning Using An Adaptive Mesh. It +/// includes both the simulation and rendering components described in that paper. +/// +/// Several pieces of software are used in this project that the respective +/// authors were kind enough to make freely available: +/// +///
    +///
  • Mersenne twister +/// (Thanks to Makoto Matsumoto) +///
  • FFTW +/// (Thanks to Matteo Frigo and Steven G. Johnson) +///
  • OpenEXR +/// (Thanks to ILM) +///
  • GLVU +/// (Thanks to Walkthru) +///
  • Antimony +/// (Thanks to Daniel Dunbar and Greg Humphreys) +///
+/// +/// Theodore Kim, kim@cs.unc.edu, October 2006 +/// +/////////////////////////////////////////////////////////////////////////////////// + +#define COMMAND_LINE_VERSION 1 + +#include +#include "ppm/ppm.hpp" +#include "APSF.h" +#include "FFT.h" +#include "QUAD_DBM_2D.h" +#include "EXR.h" + +using namespace std; + +//////////////////////////////////////////////////////////////////////////// +// globals +//////////////////////////////////////////////////////////////////////////// +int iterations = 10; +static QUAD_DBM_2D* potential = new QUAD_DBM_2D(256, 256, iterations); +APSF apsf(512); + +// input image info +int inputWidth = -1; +int inputHeight = -1; + +// input params +string inputFile; +string outputFile; + +// image scale +float scale = 5; + +// pause the simulation? +bool pause = false; + +//////////////////////////////////////////////////////////////////////////// +// render the glow +//////////////////////////////////////////////////////////////////////////// +void renderGlow(string filename, int scale = 1) +{ + int w = potential->xDagRes() * scale; + int h = potential->yDagRes() * scale; + + // draw the DAG + float*& source = potential->renderOffscreen(scale); + + // if there is no input dimensions specified, else there were input + // image dimensions, so crop it + if (inputWidth == -1) + { + inputWidth = potential->inputWidth(); + inputHeight = potential->inputHeight(); + } + + // copy out the cropped version + int wCropped = inputWidth * scale; + int hCropped = inputHeight * scale; + float* cropped = new float[wCropped * hCropped]; + cout << endl << " Generating EXR image width: " << wCropped << " height: " << hCropped << endl; + for (int y = 0; y < hCropped; y++) + for (int x = 0; x < wCropped; x++) + { + int uncroppedIndex = x + y * w; + int croppedIndex = x + y * wCropped; + cropped[croppedIndex] = source[uncroppedIndex]; + } + + // create the filter + apsf.generateKernelFast(); + + // convolve with FFT + bool success = FFT::convolve(cropped, apsf.kernel(), wCropped, hCropped, apsf.res(), apsf.res()); + + if (success) { + EXR::writeEXR(filename.c_str(), cropped, wCropped, hCropped); + cout << " " << filename << " written." << endl; + } + else + cout << " Final image generation failed." << endl; + + delete[] cropped; +} + +//////////////////////////////////////////////////////////////////////////// +// load image file into the DBM simulation +//////////////////////////////////////////////////////////////////////////// +bool loadImages(string inputFile) +{ + // load the files + unsigned char* input = NULL; + LoadPPM(inputFile.c_str(), input, inputWidth, inputHeight); + + unsigned char* start = new unsigned char[inputWidth * inputHeight]; + unsigned char* repulsor = new unsigned char[inputWidth * inputHeight]; + unsigned char* attractor = new unsigned char[inputWidth * inputHeight]; + unsigned char* terminators = new unsigned char[inputWidth * inputHeight]; + + // composite RGB channels into one + for (int x = 0; x < inputWidth * inputHeight; x++) + { + start[x] = (input[3 * x] == 255) ? 255 : 0; + repulsor[x] = (input[3 * x + 1] == 255) ? 255 : 0; + attractor[x] = (input[3 * x + 2] == 255) ? 255 : 0; + terminators[x] = 0; + + if (input[3 * x] + input[3 * x + 1] + input[3 * x + 2] == 255 * 3) + { + terminators[x] = 255; + start[x] = repulsor[x] = attractor[x] = 0; + } + } + + if (potential) delete potential; + potential = new QUAD_DBM_2D(inputWidth, inputHeight, iterations); + bool success = potential->readImage(start, attractor, repulsor, terminators, inputWidth, inputHeight); + + // delete the memory + delete[] input; + delete[] start; + delete[] repulsor; + delete[] attractor; + delete[] terminators; + + return success; +} + +int width = 600; +int height = 600; +bool animate = false; +float camera[2]; +float translate[2]; + +//////////////////////////////////////////////////////////////////////////// +// window Reshape function +//////////////////////////////////////////////////////////////////////////// +void Reshape(int w, int h) +{ + if (h == 0) h = 1; + + glViewport(0, 0, w, h); + + glMatrixMode(GL_PROJECTION); + glLoadIdentity(); + gluOrtho2D(-camera[0] - translate[0], camera[0] + translate[0], -camera[1] - translate[1], camera[1] + translate[1]); + + glMatrixMode(GL_MODELVIEW); + glLoadIdentity(); +} + +//////////////////////////////////////////////////////////////////////////// +// GLUT Display callback +//////////////////////////////////////////////////////////////////////////// +void Display() +{ + glMatrixMode(GL_PROJECTION); + glLoadIdentity(); + gluOrtho2D(-camera[0] + translate[0], camera[0] + translate[0], -camera[1] + translate[1], camera[1] + translate[1]); + + glMatrixMode(GL_MODELVIEW); + glLoadIdentity(); + + glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); + glPolygonMode(GL_FRONT_AND_BACK, GL_FILL); + + potential->draw(); + potential->drawSegments(); + + glutSwapBuffers(); +} + +//////////////////////////////////////////////////////////////////////////// +// GLUT Keyboard callback +//////////////////////////////////////////////////////////////////////////// +void Keyboard(unsigned char key, int x, int y) +{ + switch(key) + { + case 'p': + pause = !pause; + break; + + case 'q': + cout << " You terminated the simulation prematurely." << endl; + exit(0); + break; + } + + glutPostRedisplay(); +} + +//////////////////////////////////////////////////////////////////////////// +// window Reshape function +//////////////////////////////////////////////////////////////////////////// +void Idle() +{ + if (!pause) + for (int x = 0; x < 100; x++) + { + bool success = potential->addParticle(); + + if (!success) + { + cout << " No nodes left to add! Is your terminator reachable?" << endl; + //exit(1); + return; + } + + if (potential->hitGround()) + { + glutPostRedisplay(); + cout << endl << endl; + + // write out the DAG file + string lightningFile = inputFile.substr(0, inputFile.size() - 3) + string("lightning"); + cout << " Intermediate file " << lightningFile << " written." << endl; + potential->writeDAG(lightningFile.c_str()); + + // render the final EXR file + renderGlow(outputFile, scale); + delete potential; + exit(0); + } + } + glutPostRedisplay(); +} + +//////////////////////////////////////////////////////////////////////////// +// GLUT Main +//////////////////////////////////////////////////////////////////////////// +int glutMain(int *argc, char **argv) +{ + float smaller = 1.0f; + camera[0] = smaller * 0.5f; + camera[1] = smaller * 0.5f; + translate[0] = 0.0f; + translate[1] = 0.0f; + + glutInit(argc, argv); + glutInitDisplayMode(GLUT_DOUBLE | GLUT_DEPTH | GLUT_RGBA); + glutInitWindowPosition(50, 50); + glutInitWindowSize(width, height); + glutCreateWindow("Lumos: A Lightning Generator v0.1"); + + glutDisplayFunc(Display); + glutKeyboardFunc(Keyboard); + glutIdleFunc(Idle); + glEnable(GL_BLEND); + glBlendFunc(GL_SRC_ALPHA,GL_ONE); + + Reshape(width, height); + + glClearColor(0.0, 0.0, 0.0, 1.0); + glShadeModel(GL_SMOOTH); + + // Go! + glutMainLoop(); + + return 0; +} + +//////////////////////////////////////////////////////////////////////////// +// Main +//////////////////////////////////////////////////////////////////////////// +int main(int argc, char **argv) +{ + if (argc < 3) + { + cout << endl; + cout << " LumosQuad " << endl; + cout << " =========================================================" << endl; + cout << " - *.ppm file with input colors" << endl; + cout << " --OR--" << endl; + cout << " *.lightning file from a previous run" << endl; + cout << " - The OpenEXR file to output" << endl; + cout << " - Scaling constant for final image." << endl; + cout << " Press 'q' to terminate the simulation prematurely." << endl; + cout << " Send questions and comments to kim@cs.unc.edu" << endl; + return 1; + } + + cout << endl << "Lumos: A lightning generator v0.1" << endl; + cout << "------------------------------------------------------" << endl; + + // store the input params + inputFile = string(argv[1]); + outputFile = string(argv[2]); + if (argc > 3) scale = atoi(argv[3]); + + // see if the input is a *.lightning file + if (inputFile.size() > 10) + { + string postfix = inputFile.substr(inputFile.size() - 9, inputFile.size()); + + cout << " Using intermediate file " << inputFile << endl; + if (postfix == string("lightning")) + { + potential->readDAG(inputFile.c_str()); + renderGlow(outputFile, scale); + delete potential; + return 0; + } + } + + // read in the *.ppm input file + if (!loadImages(inputFile)) + { + cout << " ERROR: " << inputFile.c_str() << " is not a valid PPM file." << endl; + return 1; + } + cout << " " << inputFile << " read." << endl << endl; + + + // loop simulation until it hits a terminator + cout << " Total particles added: "; + glutMain(&argc, argv); + + return 0; +} diff --git a/ppm/ppm.cpp b/ppm/ppm.cpp new file mode 100644 index 0000000..c95b50d --- /dev/null +++ b/ppm/ppm.cpp @@ -0,0 +1,60 @@ +//------------------------------------------------------------------------------ +// File : ppm.cpp +//------------------------------------------------------------------------------ +// GLVU : Copyright 1997 - 2002 +// The University of North Carolina at Chapel Hill +//------------------------------------------------------------------------------ +// Permission to use, copy, modify, distribute and sell this software and its +// documentation for any purpose is hereby granted without fee, provided that +// the above copyright notice appear in all copies and that both that copyright +// notice and this permission notice appear in supporting documentation. +// Binaries may be compiled with this software without any royalties or +// restrictions. +// +// The University of North Carolina at Chapel Hill makes no representations +// about the suitability of this software for any purpose. It is provided +// "as is" without express or implied warranty. + +//============================================================================ +// ppm.cpp : Portable Pixel Map image format module +//============================================================================ + +#include + +#define PPM_VERBOSE 0 + +//---------------------------------------------------------------------------- +// READS AN IMAGE IN FROM A PPM FILE. RETURNS THE COLOR RGB ARRAY AND DIMENSIONS +// PERFORMS AUTO-ALLOCATION OF Color ARRAY IF SET TO NULL BEFORE CALLING; OTHERWISE +// ASSUMES THAT COLOR HAS BEEN PRE-ALLOCED. +//---------------------------------------------------------------------------- +void LoadPPM(const char *FileName, unsigned char* &Color, int &Width, int &Height) +{ + FILE* fp = fopen(FileName, "rb"); + if (fp==NULL) + { printf("PPM ERROR (ReadPPM) : unable to open %s!\n",FileName); + Color=NULL; Width=0; Height=0; return; } + int c,s; + do{ do { s=fgetc(fp); } while (s!='\n'); } while ((c=fgetc(fp))=='#'); + ungetc(c,fp); + fscanf(fp, "%d %d\n255\n", &Width, &Height); +#if PPM_VERBOSE + printf("Reading %dx%d Texture [%s]. . .\n", Width, Height, FileName); +#endif + int NumComponents = Width*Height*3; + if (Color==NULL) Color = new unsigned char[NumComponents]; + fread(Color,NumComponents,1,fp); + fclose(fp); +} + +//---------------------------------------------------------------------------- +// Writes an unsigned byte RGB color array out to a PPM file. +//---------------------------------------------------------------------------- +void WritePPM(const char *FileName, unsigned char* Color, int Width, int Height) +{ + FILE* fp = fopen(FileName, "wb"); + if (fp==NULL) { printf("PPM ERROR (WritePPM) : unable to open %s!\n",FileName); return; } + fprintf(fp, "P6\n%d %d\n255\n", Width, Height); + fwrite(Color,1,Width*Height*3,fp); + fclose(fp); +} diff --git a/ppm/ppm.hpp b/ppm/ppm.hpp new file mode 100644 index 0000000..c887b0f --- /dev/null +++ b/ppm/ppm.hpp @@ -0,0 +1,23 @@ +//------------------------------------------------------------------------------ +// File : ppm.hpp +//------------------------------------------------------------------------------ +// GLVU : Copyright 1997 - 2002 +// The University of North Carolina at Chapel Hill +//------------------------------------------------------------------------------ +// Permission to use, copy, modify, distribute and sell this software and its +// documentation for any purpose is hereby granted without fee, provided that +// the above copyright notice appear in all copies and that both that copyright +// notice and this permission notice appear in supporting documentation. +// Binaries may be compiled with this software without any royalties or +// restrictions. +// +// The University of North Carolina at Chapel Hill makes no representations +// about the suitability of this software for any purpose. It is provided +// "as is" without express or implied warranty. + +//============================================================================ +// ppm.hpp : Portable Pixel Map image format module +//============================================================================ + +void LoadPPM(const char *FileName, unsigned char* &Color, int &Width, int &Height); +void WritePPM(const char *FileName, unsigned char* Color, int Width, int Height);