diff --git a/src/plugins_experimental/filter_screened_poisson/Src/CmdLineParser.cpp b/src/plugins_experimental/filter_screened_poisson/Src/CmdLineParser.cpp
new file mode 100755
index 000000000..9061b4a85
--- /dev/null
+++ b/src/plugins_experimental/filter_screened_poisson/Src/CmdLineParser.cpp
@@ -0,0 +1,279 @@
+/*
+Copyright (c) 2006, Michael Kazhdan and Matthew Bolitho
+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 the Johns Hopkins University nor the names of its contributors
+may 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.
+*/
+
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <assert.h>
+#include "CmdLineParser.h"
+
+
+#ifdef WIN32
+int strcasecmp(char* c1,char* c2){return _stricmp(c1,c2);}
+#endif
+
+cmdLineReadable::cmdLineReadable(const char* name)
+{
+	set=false;
+	this->name=new char[strlen(name)+1];
+	strcpy(this->name,name);
+}
+cmdLineReadable::~cmdLineReadable(void)
+{
+	if(name) delete[] name;
+	name=NULL;
+}
+int cmdLineReadable::read(char**,int){
+	set=true;
+	return 0;
+}
+void cmdLineReadable::writeValue(char* str)
+{
+	str[0] = 0;
+}
+
+////////////////
+// cmdLineInt //
+////////////////
+cmdLineInt::cmdLineInt(const char* name) : cmdLineReadable(name) {value=0;}
+cmdLineInt::cmdLineInt(const char* name,const int& v) : cmdLineReadable(name) {value=v;}
+int cmdLineInt::read(char** argv,int argc){
+	if(argc>0){
+		value=atoi(argv[0]);
+		set=true;
+		return 1;
+	}
+	else{return 0;}
+}
+void cmdLineInt::writeValue(char* str)
+{
+	sprintf(str,"%d",value);
+}
+
+//////////////////
+// cmdLineFloat //
+//////////////////
+cmdLineFloat::cmdLineFloat(const char* name) : cmdLineReadable(name) {value=0;}
+cmdLineFloat::cmdLineFloat(const char* name, const float& v) : cmdLineReadable(name) {value=v;}
+int cmdLineFloat::read(char** argv,int argc){
+	if(argc>0){
+		value=(float)atof(argv[0]);
+		set=true;
+		return 1;
+	}
+	else{return 0;}
+}
+void cmdLineFloat::writeValue(char* str)
+{
+	sprintf(str,"%f",value);
+}
+
+///////////////////
+// cmdLineString //
+///////////////////
+cmdLineString::cmdLineString(const char* name) : cmdLineReadable(name) {value=NULL;}
+cmdLineString::~cmdLineString(void)
+{
+	if(value)	delete[] value;
+	value=NULL;
+}
+int cmdLineString::read(char** argv,int argc){
+	if(argc>0)
+	{
+		value=new char[strlen(argv[0])+1];
+		strcpy(value,argv[0]);
+		set=true;
+		return 1;
+	}
+	else{return 0;}
+}
+void cmdLineString::writeValue(char* str)
+{
+	sprintf(str,"%s",value);
+}
+
+////////////////////
+// cmdLineStrings //
+////////////////////
+cmdLineStrings::cmdLineStrings(const char* name,int Dim) : cmdLineReadable(name)
+{
+	this->Dim=Dim;
+	values=new char*[Dim];
+	for(int i=0;i<Dim;i++)	values[i]=NULL;
+}
+cmdLineStrings::~cmdLineStrings(void)
+{
+	for(int i=0;i<Dim;i++)
+	{
+		if(values[i])	delete[] values[i];
+		values[i]=NULL;
+	}
+	delete[] values;
+	values=NULL;
+}
+int cmdLineStrings::read(char** argv,int argc)
+{
+	if(argc>=Dim)
+	{
+		for(int i=0;i<Dim;i++)
+		{
+			values[i]=new char[strlen(argv[i])+1];
+			strcpy(values[i],argv[i]);
+		}
+		set=true;
+		return Dim;
+	}
+	else	return 0;
+}
+void cmdLineStrings::writeValue(char* str)
+{
+	char* temp=str;
+	for(int i=0;i<Dim;i++)
+	{
+		sprintf(temp,"%s ",values[i]);
+		temp=str+strlen(str);
+	}
+}
+
+
+char* GetFileExtension(char* fileName){
+	char* fileNameCopy;
+	char* ext=NULL;
+	char* temp;
+
+	fileNameCopy=new char[strlen(fileName)+1];
+	assert(fileNameCopy);
+	strcpy(fileNameCopy,fileName);
+	temp=strtok(fileNameCopy,".");
+	while(temp!=NULL)
+	{
+		if(ext!=NULL){delete[] ext;}
+		ext=new char[strlen(temp)+1];
+		assert(ext);
+		strcpy(ext,temp);
+		temp=strtok(NULL,".");
+	}
+	delete[] fileNameCopy;
+	return ext;
+}
+char* GetLocalFileName(char* fileName){
+	char* fileNameCopy;
+	char* name=NULL;
+	char* temp;
+
+	fileNameCopy=new char[strlen(fileName)+1];
+	assert(fileNameCopy);
+	strcpy(fileNameCopy,fileName);
+	temp=strtok(fileNameCopy,"\\");
+	while(temp!=NULL){
+		if(name!=NULL){delete[] name;}
+		name=new char[strlen(temp)+1];
+		assert(name);
+		strcpy(name,temp);
+		temp=strtok(NULL,"\\");
+	}
+	delete[] fileNameCopy;
+	return name;
+}
+
+void cmdLineParse(int argc, char **argv,int num,cmdLineReadable** readable,int dumpError)
+{
+	int i,j;
+	while (argc > 0)
+	{
+		if (argv[0][0] == '-' && argv[0][1]=='-')
+		{
+			for(i=0;i<num;i++)
+			{
+				if (!strcmp(&argv[0][2],readable[i]->name))
+				{
+					argv++, argc--;
+					j=readable[i]->read(argv,argc);
+					argv+=j,argc-=j;
+					break;
+				}
+			}
+			if(i==num){
+				if(dumpError)
+				{
+					fprintf(stderr, "invalid option: %s\n",*argv);
+					fprintf(stderr, "possible options are:\n");
+					for(i=0;i<num;i++)	fprintf(stderr, "  %s\n",readable[i]->name);
+				}
+				argv++, argc--;
+			}
+		}
+		else
+		{
+			if(dumpError)
+			{
+				fprintf(stderr, "invalid option: %s\n", *argv);
+				fprintf(stderr, "  options must start with a \'--\'\n");
+			}
+			argv++, argc--;
+		}
+	}
+}
+char** ReadWords(const char* fileName,int& cnt)
+{
+	char** names;
+	char temp[500];
+	FILE* fp;
+
+	fp=fopen(fileName,"r");
+	if(!fp){return NULL;}
+	cnt=0;
+	while(fscanf(fp," %s ",temp)==1){cnt++;}
+	fclose(fp);
+
+	names=new char*[cnt];
+	if(!names){return NULL;}
+
+	fp=fopen(fileName,"r");
+	if(!fp){
+		delete[] names;
+		cnt=0;
+		return NULL;
+	}
+	cnt=0;
+	while(fscanf(fp," %s ",temp)==1){
+		names[cnt]=new char[strlen(temp)+1];
+		if(!names){
+			for(int j=0;j<cnt;j++){delete[] names[j];}
+			delete[] names;
+			cnt=0;
+			fclose(fp);
+			return NULL;
+		}
+		strcpy(names[cnt],temp);
+		cnt++;
+	}
+	fclose(fp);
+	return names;
+}
diff --git a/src/plugins_experimental/filter_screened_poisson/Src/CmdLineParser.h b/src/plugins_experimental/filter_screened_poisson/Src/CmdLineParser.h
new file mode 100755
index 000000000..6d8b1759f
--- /dev/null
+++ b/src/plugins_experimental/filter_screened_poisson/Src/CmdLineParser.h
@@ -0,0 +1,120 @@
+/*
+Copyright (c) 2006, Michael Kazhdan and Matthew Bolitho
+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 the Johns Hopkins University nor the names of its contributors
+may 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.
+*/
+
+#ifndef CMD_LINE_PARSER_INCLUDED
+#define CMD_LINE_PARSER_INCLUDED
+#include <stdarg.h>
+#include <string.h>
+
+
+#ifdef WIN32
+int strcasecmp(char* c1,char* c2);
+#endif
+
+class cmdLineReadable{
+public:
+	bool set;
+	char* name;
+	cmdLineReadable(const char* name);
+	virtual ~cmdLineReadable(void);
+	virtual int read(char** argv,int argc);
+	virtual void writeValue(char* str);
+};
+
+class cmdLineInt : public cmdLineReadable {
+public:
+	int value;
+	cmdLineInt(const char* name);
+	cmdLineInt(const char* name,const int& v);
+	int read(char** argv,int argc);
+	void writeValue(char* str);
+};
+template<int Dim>
+class cmdLineIntArray : public cmdLineReadable {
+public:
+	int values[Dim];
+	cmdLineIntArray(const char* name);
+	cmdLineIntArray(const char* name,const int v[Dim]);
+	int read(char** argv,int argc);
+	void writeValue(char* str);
+};
+
+class cmdLineFloat : public cmdLineReadable {
+public:
+	float value;
+	cmdLineFloat(const char* name);
+	cmdLineFloat(const char* name,const float& f);
+	int read(char** argv,int argc);
+	void writeValue(char* str);
+};
+template<int Dim>
+class cmdLineFloatArray : public cmdLineReadable {
+public:
+	float values[Dim];
+	cmdLineFloatArray(const char* name);
+	cmdLineFloatArray(const char* name,const float f[Dim]);
+	int read(char** argv,int argc);
+	void writeValue(char* str);
+};
+class cmdLineString : public cmdLineReadable {
+public:
+	char* value;
+	cmdLineString(const char* name);
+	~cmdLineString();
+	int read(char** argv,int argc);
+	void writeValue(char* str);
+};
+class cmdLineStrings : public cmdLineReadable {
+	int Dim;
+public:
+	char** values;
+	cmdLineStrings(const char* name,int Dim);
+	~cmdLineStrings(void);
+	int read(char** argv,int argc);
+	void writeValue(char* str);
+};
+template<int Dim>
+class cmdLineStringArray : public cmdLineReadable {
+public:
+	char* values[Dim];
+	cmdLineStringArray(const char* name);
+	~cmdLineStringArray();
+	int read(char** argv,int argc);
+	void writeValue(char* str);
+};
+
+// This reads the arguments in argc, matches them against "names" and sets
+// the values of "r" appropriately. Parameters start with "--"
+void cmdLineParse(int argc, char **argv,int num,cmdLineReadable** r,int dumpError=1);
+
+char* GetFileExtension(char* fileName);
+char* GetLocalFileName(char* fileName);
+char** ReadWords(const char* fileName,int& cnt);
+
+#include "CmdLineParser.inl"
+#endif // CMD_LINE_PARSER_INCLUDED
diff --git a/src/plugins_experimental/filter_screened_poisson/Src/CmdLineParser.inl b/src/plugins_experimental/filter_screened_poisson/Src/CmdLineParser.inl
new file mode 100755
index 000000000..eeded6805
--- /dev/null
+++ b/src/plugins_experimental/filter_screened_poisson/Src/CmdLineParser.inl
@@ -0,0 +1,141 @@
+/* -*- C++ -*-
+Copyright (c) 2006, Michael Kazhdan and Matthew Bolitho
+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 the Johns Hopkins University nor the names of its contributors
+may 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.
+*/
+
+/////////////////////
+// cmdLineIntArray //
+/////////////////////
+template<int Dim>
+cmdLineIntArray<Dim>::cmdLineIntArray(const char* name) : cmdLineReadable(name)
+{
+	for(int i=0;i<Dim;i++)	values[i]=0;
+}
+template<int Dim>
+cmdLineIntArray<Dim>::cmdLineIntArray(const char* name,const int v[Dim]) : cmdLineReadable(name)
+{
+	for(int i=0;i<Dim;i++)	values[i]=v[i];
+}
+template<int Dim>
+int cmdLineIntArray<Dim>::read(char** argv,int argc)
+{
+	if(argc>=Dim)
+	{
+		for(int i=0;i<Dim;i++)	values[i]=atoi(argv[i]);
+		set=true;
+		return Dim;
+	}
+	else{return 0;}
+}
+template<int Dim>
+void cmdLineIntArray<Dim>::writeValue(char* str)
+{
+	char* temp=str;
+	for(int i=0;i<Dim;i++)
+	{
+		sprintf(temp,"%d ",values[i]);
+		temp=str+strlen(str);
+	}
+}
+
+///////////////////////
+// cmdLineFloatArray //
+///////////////////////
+template<int Dim>
+cmdLineFloatArray<Dim>::cmdLineFloatArray(const char* name) : cmdLineReadable(name)
+{
+	for(int i=0;i<Dim;i++)	values[i]=0;
+}
+template<int Dim>
+cmdLineFloatArray<Dim>::cmdLineFloatArray(const char* name,const float f[Dim]) : cmdLineReadable(name)
+{
+	for(int i=0;i<Dim;i++)	values[i]=f[i];
+}
+template<int Dim>
+int cmdLineFloatArray<Dim>::read(char** argv,int argc)
+{
+	if(argc>=Dim)
+	{
+		for(int i=0;i<Dim;i++)	values[i]=(float)atof(argv[i]);
+		set=true;
+		return Dim;
+	}
+	else{return 0;}
+}
+template<int Dim>
+void cmdLineFloatArray<Dim>::writeValue(char* str)
+{
+	char* temp=str;
+	for(int i=0;i<Dim;i++)
+	{
+		sprintf(temp,"%f ",values[i]);
+		temp=str+strlen(str);
+	}
+}
+
+
+////////////////////////
+// cmdLineStringArray //
+////////////////////////
+template<int Dim>
+cmdLineStringArray<Dim>::cmdLineStringArray(const char* name) : cmdLineReadable(name)
+{
+	for(int i=0;i<Dim;i++)	values[i]=NULL;
+}
+template<int Dim>
+cmdLineStringArray<Dim>::~cmdLineStringArray(void)
+{
+	for(int i=0;i<Dim;i++)
+	{
+		if(values[i])	delete[] values[i];
+		values[i]=NULL;
+	}
+}
+template<int Dim>
+int cmdLineStringArray<Dim>::read(char** argv,int argc)
+{
+	if(argc>=Dim)
+	{
+		for(int i=0;i<Dim;i++)
+		{
+			values[i]=new char[strlen(argv[i])+1];
+			strcpy(values[i],argv[i]);
+		}
+		set=true;
+		return Dim;
+	}
+	else{return 0;}
+}
+template<int Dim>
+void cmdLineStringArray<Dim>::writeValue(char* str)
+{
+	char* temp=str;
+	for(int i=0;i<Dim;i++)
+	{
+		sprintf(temp,"%s ",values[i]);
+		temp=str+strlen(str);
+	}
+}
diff --git a/src/plugins_experimental/filter_screened_poisson/Src/Factor.cpp b/src/plugins_experimental/filter_screened_poisson/Src/Factor.cpp
new file mode 100755
index 000000000..145c100dc
--- /dev/null
+++ b/src/plugins_experimental/filter_screened_poisson/Src/Factor.cpp
@@ -0,0 +1,264 @@
+/*
+Copyright (c) 2006, Michael Kazhdan and Matthew Bolitho
+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 the Johns Hopkins University nor the names of its contributors
+may 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.
+*/
+
+//////////////////////
+// Polynomial Roots //
+//////////////////////
+#include <math.h>
+#include "Factor.h"
+int Factor(double a1,double a0,double roots[1][2],double EPS){
+	if(fabs(a1)<=EPS){return 0;}
+	roots[0][0]=-a0/a1;
+	roots[0][1]=0;
+	return 1;
+}
+int Factor(double a2,double a1,double a0,double roots[2][2],double EPS){
+	double d;
+	if(fabs(a2)<=EPS){return Factor(a1,a0,roots,EPS);}
+
+	d=a1*a1-4*a0*a2;
+	a1/=(2*a2);
+	if(d<0){
+		d=sqrt(-d)/(2*a2);
+		roots[0][0]=roots[1][0]=-a1;
+		roots[0][1]=-d;
+		roots[1][1]= d;
+	}
+	else{
+		d=sqrt(d)/(2*a2);
+		roots[0][1]=roots[1][1]=0;
+		roots[0][0]=-a1-d;
+		roots[1][0]=-a1+d;
+	}
+	return 2;
+}
+// Solution taken from: http://mathworld.wolfram.com/CubicFormula.html
+// and http://www.csit.fsu.edu/~burkardt/f_src/subpak/subpak.f90
+int Factor(double a3,double a2,double a1,double a0,double roots[3][2],double EPS){
+	double q,r,r2,q3;
+
+	if(fabs(a3)<=EPS){return Factor(a2,a1,a0,roots,EPS);}
+	a2/=a3;
+	a1/=a3;
+	a0/=a3;
+
+	q=-(3*a1-a2*a2)/9;
+	r=-(9*a2*a1-27*a0-2*a2*a2*a2)/54;
+	r2=r*r;
+	q3=q*q*q;
+
+	if(r2<q3){
+		double sqrQ=sqrt(q);
+		double theta = acos ( r / (sqrQ*q) );
+		double cTheta=cos(theta/3)*sqrQ;
+		double sTheta=sin(theta/3)*sqrQ*SQRT_3/2;
+		roots[0][1]=roots[1][1]=roots[2][1]=0;
+		roots[0][0]=-2*cTheta;
+		roots[1][0]=-2*(-cTheta*0.5-sTheta);
+		roots[2][0]=-2*(-cTheta*0.5+sTheta);
+	}
+	else{
+		double s1,s2,sqr=sqrt(r2-q3);
+		double t;
+		t=-r+sqr;
+		if(t<0){s1=-pow(-t,1.0/3);}
+		else{s1=pow(t,1.0/3);}
+		t=-r-sqr;
+		if(t<0){s2=-pow(-t,1.0/3);}
+		else{s2=pow(t,1.0/3);}
+		roots[0][1]=0;
+		roots[0][0]=s1+s2;
+		s1/=2;
+		s2/=2;
+		roots[1][0]= roots[2][0]=-s1-s2;
+		roots[1][1]= SQRT_3*(s1-s2);
+		roots[2][1]=-roots[1][1];
+	}
+	roots[0][0]-=a2/3;
+	roots[1][0]-=a2/3;
+	roots[2][0]-=a2/3;
+	return 3;
+}
+double ArcTan2(double y,double x){
+	/* This first case should never happen */
+	if(y==0 && x==0){return 0;}
+	if(x==0){
+		if(y>0){return PI/2.0;}
+		else{return -PI/2.0;}
+	}
+	if(x>=0){return atan(y/x);}
+	else{
+		if(y>=0){return atan(y/x)+PI;}
+		else{return atan(y/x)-PI;}
+	}
+}
+double Angle(const double in[2]){
+	if((in[0]*in[0]+in[1]*in[1])==0.0){return 0;}
+	else{return ArcTan2(in[1],in[0]);}
+}
+void Sqrt(const double in[2],double out[2]){
+	double r=sqrt(sqrt(in[0]*in[0]+in[1]*in[1]));
+	double a=Angle(in)*0.5;
+	out[0]=r*cos(a);
+	out[1]=r*sin(a);
+}
+void Add(const double in1[2],const double in2[2],double out[2]){
+	out[0]=in1[0]+in2[0];
+	out[1]=in1[1]+in2[1];
+}
+void Subtract(const double in1[2],const double in2[2],double out[2]){
+	out[0]=in1[0]-in2[0];
+	out[1]=in1[1]-in2[1];
+}
+void Multiply(const double in1[2],const double in2[2],double out[2]){
+	out[0]=in1[0]*in2[0]-in1[1]*in2[1];
+	out[1]=in1[0]*in2[1]+in1[1]*in2[0];
+}
+void Divide(const double in1[2],const double in2[2],double out[2]){
+	double temp[2];
+	double l=in2[0]*in2[0]+in2[1]*in2[1];
+	temp[0]= in2[0]/l;
+	temp[1]=-in2[1]/l;
+	Multiply(in1,temp,out);
+}
+// Solution taken from: http://mathworld.wolfram.com/QuarticEquation.html
+// and http://www.csit.fsu.edu/~burkardt/f_src/subpak/subpak.f90
+int Factor(double a4,double a3,double a2,double a1,double a0,double roots[4][2],double EPS){
+	double R[2],D[2],E[2],R2[2];
+
+	if(fabs(a4)<EPS){return Factor(a3,a2,a1,a0,roots,EPS);}
+	a3/=a4;
+	a2/=a4;
+	a1/=a4;
+	a0/=a4;
+
+	Factor(1.0,-a2,a3*a1-4.0*a0,-a3*a3*a0+4.0*a2*a0-a1*a1,roots,EPS);
+
+	R2[0]=a3*a3/4.0-a2+roots[0][0];
+	R2[1]=0;
+	Sqrt(R2,R);
+	if(fabs(R[0])>10e-8){
+		double temp1[2],temp2[2];
+		double p1[2],p2[2];
+
+		p1[0]=a3*a3*0.75-2.0*a2-R2[0];
+		p1[1]=0;
+
+		temp2[0]=((4.0*a3*a2-8.0*a1-a3*a3*a3)/4.0);
+		temp2[1]=0;
+		Divide(temp2,R,p2);
+
+		Add     (p1,p2,temp1);
+		Subtract(p1,p2,temp2);
+
+		Sqrt(temp1,D);
+		Sqrt(temp2,E);
+	}
+	else{
+		R[0]=R[1]=0;
+		double temp1[2],temp2[2];
+		temp1[0]=roots[0][0]*roots[0][0]-4.0*a0;
+		temp1[1]=0;
+		Sqrt(temp1,temp2);
+		temp1[0]=a3*a3*0.75-2.0*a2+2.0*temp2[0];
+		temp1[1]=                  2.0*temp2[1];
+		Sqrt(temp1,D);
+		temp1[0]=a3*a3*0.75-2.0*a2-2.0*temp2[0];
+		temp1[1]=                 -2.0*temp2[1];
+		Sqrt(temp1,E);
+	}
+
+	roots[0][0]=-a3/4.0+R[0]/2.0+D[0]/2.0;
+	roots[0][1]=        R[1]/2.0+D[1]/2.0;
+
+	roots[1][0]=-a3/4.0+R[0]/2.0-D[0]/2.0;
+	roots[1][1]=        R[1]/2.0-D[1]/2.0;
+
+	roots[2][0]=-a3/4.0-R[0]/2.0+E[0]/2.0;
+	roots[2][1]=       -R[1]/2.0+E[1]/2.0;
+
+	roots[3][0]=-a3/4.0-R[0]/2.0-E[0]/2.0;
+	roots[3][1]=       -R[1]/2.0-E[1]/2.0;
+	return 4;
+}
+
+int Solve(const double* eqns,const double* values,double* solutions,int dim){
+	int i,j,eIndex;
+	double v,m;
+	int *index=new int[dim];
+	int *set=new int[dim];
+	double* myEqns=new double[dim*dim];
+	double* myValues=new double[dim];
+
+	for(i=0;i<dim*dim;i++){myEqns[i]=eqns[i];}
+	for(i=0;i<dim;i++){
+		myValues[i]=values[i];
+		set[i]=0;
+	}
+	for(i=0;i<dim;i++){
+		// Find the largest equation that has a non-zero entry in the i-th index
+		m=-1;
+		eIndex=-1;
+		for(j=0;j<dim;j++){
+			if(set[j]){continue;}
+			if(myEqns[j*dim+i]!=0 && fabs(myEqns[j*dim+i])>m){
+				m=fabs(myEqns[j*dim+i]);
+				eIndex=j;
+			}
+		}
+		if(eIndex==-1){
+			delete[] index;
+			delete[] myValues;
+			delete[] myEqns;
+			delete[] set;
+			return 0;
+		}
+		// The position in which the solution for the i-th variable can be found
+		index[i]=eIndex;
+		set[eIndex]=1;
+
+		// Normalize the equation
+		v=myEqns[eIndex*dim+i];
+		for(j=0;j<dim;j++){myEqns[eIndex*dim+j]/=v;}
+		myValues[eIndex]/=v;
+
+		// Subtract it off from everything else
+		for(j=0;j<dim;j++){
+			if(j==eIndex){continue;}
+			double vv=myEqns[j*dim+i];
+			for(int k=0;k<dim;k++){myEqns[j*dim+k]-=myEqns[eIndex*dim+k]*vv;}
+			myValues[j]-=myValues[eIndex]*vv;
+		}
+	}
+	for(i=0;i<dim;i++){solutions[i]=myValues[index[i]];}
+	delete[] index;
+	delete[] myValues;
+	delete[] myEqns;
+	delete[] set;
+	return 1;
+}
diff --git a/src/plugins_experimental/filter_screened_poisson/Src/Factor.h b/src/plugins_experimental/filter_screened_poisson/Src/Factor.h
new file mode 100755
index 000000000..d6ee4b200
--- /dev/null
+++ b/src/plugins_experimental/filter_screened_poisson/Src/Factor.h
@@ -0,0 +1,50 @@
+/*
+Copyright (c) 2006, Michael Kazhdan and Matthew Bolitho
+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 the Johns Hopkins University nor the names of its contributors
+may 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.
+*/
+
+#ifndef FACTOR_INCLUDED
+#define FACTOR_INCLUDED
+
+#define PI 3.1415926535897932384
+#define SQRT_3 1.7320508075688772935
+
+double ArcTan2(double y,double x);
+double Angle(const double in[2]);
+void Sqrt(const double in[2],double out[2]);
+void Add(const double in1[2],const double in2[2],double out[2]);
+void Subtract(const double in1[2],const double in2[2],double out[2]);
+void Multiply(const double in1[2],const double in2[2],double out[2]);
+void Divide(const double in1[2],const double in2[2],double out[2]);
+
+int Factor(double a1,double a0,double roots[1][2],double EPS);
+int Factor(double a2,double a1,double a0,double roots[2][2],double EPS);
+int Factor(double a3,double a2,double a1,double a0,double roots[3][2],double EPS);
+int Factor(double a4,double a3,double a2,double a1,double a0,double roots[4][2],double EPS);
+
+int Solve(const double* eqns,const double* values,double* solutions, int dim);
+
+#endif // FACTOR_INCLUDED
\ No newline at end of file
diff --git a/src/plugins_experimental/filter_screened_poisson/Src/FunctionData.inl b/src/plugins_experimental/filter_screened_poisson/Src/FunctionData.inl
new file mode 100755
index 000000000..c8d5293ab
--- /dev/null
+++ b/src/plugins_experimental/filter_screened_poisson/Src/FunctionData.inl
@@ -0,0 +1,415 @@
+/*
+Copyright (c) 2006, Michael Kazhdan and Matthew Bolitho
+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 the Johns Hopkins University nor the names of its contributors
+may 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.
+*/
+
+//////////////////
+// FunctionData //
+//////////////////
+template<int Degree,class Real>
+FunctionData<Degree,Real>::FunctionData(void)
+{
+	dotTable=dDotTable=d2DotTable=NULL;
+	valueTables=dValueTables=NULL;
+	res=0;
+}
+
+template<int Degree,class Real>
+FunctionData<Degree,Real>::~FunctionData(void)
+{
+	if(res)
+	{
+		if(  dotTable) delete[]   dotTable;
+		if( dDotTable) delete[]  dDotTable;
+		if(d2DotTable) delete[] d2DotTable;
+		if( valueTables) delete[]  valueTables;
+		if(dValueTables) delete[] dValueTables;
+	}
+	dotTable=dDotTable=d2DotTable=NULL;
+	valueTables=dValueTables=NULL;
+	res=0;
+}
+
+template<int Degree,class Real>
+#if BOUNDARY_CONDITIONS
+void FunctionData<Degree,Real>::set( const int& maxDepth , const PPolynomial<Degree>& F , const int& normalize , bool useDotRatios , bool reflectBoundary )
+#else // !BOUNDARY_CONDITIONS
+void FunctionData<Degree,Real>::set(const int& maxDepth,const PPolynomial<Degree>& F,const int& normalize , bool useDotRatios )
+#endif // BOUNDARY_CONDITIONS
+{
+	this->normalize       = normalize;
+	this->useDotRatios    = useDotRatios;
+#if BOUNDARY_CONDITIONS
+	this->reflectBoundary = reflectBoundary;
+#endif // BOUNDARY_CONDITIONS
+
+	depth = maxDepth;
+	res = BinaryNode<double>::CumulativeCenterCount( depth );
+	res2 = (1<<(depth+1))+1;
+	baseFunctions = new PPolynomial<Degree+1>[res];
+	// Scale the function so that it has:
+	// 0] Value 1 at 0
+	// 1] Integral equal to 1
+	// 2] Square integral equal to 1
+	switch( normalize )
+	{
+		case 2:
+			baseFunction=F/sqrt((F*F).integral(F.polys[0].start,F.polys[F.polyCount-1].start));
+			break;
+		case 1:
+			baseFunction=F/F.integral(F.polys[0].start,F.polys[F.polyCount-1].start);
+			break;
+		default:
+			baseFunction=F/F(0);
+	}
+	dBaseFunction = baseFunction.derivative();
+#if BOUNDARY_CONDITIONS
+	leftBaseFunction   = baseFunction + baseFunction.shift( -1 );
+	rightBaseFunction  = baseFunction + baseFunction.shift(  1 );
+	dLeftBaseFunction  =  leftBaseFunction.derivative();
+	dRightBaseFunction = rightBaseFunction.derivative();
+#endif // BOUNDARY_CONDITIONS
+	double c1,w1;
+	for( int i=0 ; i<res ; i++ )
+	{
+		BinaryNode< double >::CenterAndWidth( i , c1 , w1 );
+#if BOUNDARY_CONDITIONS
+		if( reflectBoundary )
+		{
+			int d , off;
+			BinaryNode< double >::DepthAndOffset( i , d , off );
+			if     ( off==0          ) baseFunctions[i] =  leftBaseFunction.scale( w1 ).shift( c1 );
+			else if( off==((1<<d)-1) ) baseFunctions[i] = rightBaseFunction.scale( w1 ).shift( c1 );
+			else                       baseFunctions[i] =      baseFunction.scale( w1 ).shift( c1 );
+		}
+		else baseFunctions[i] = baseFunction.scale(w1).shift(c1);
+#else // !BOUNDARY_CONDITIONS
+		baseFunctions[i] = baseFunction.scale(w1).shift(c1);
+#endif // BOUNDARY_CONDITIONS
+		// Scale the function so that it has L2-norm equal to one
+		switch( normalize )
+		{
+			case 2:
+				baseFunctions[i]/=sqrt(w1);
+				break;
+			case 1:
+				baseFunctions[i]/=w1;
+				break;
+		}
+	}
+}
+template<int Degree,class Real>
+void FunctionData<Degree,Real>::setDotTables( const int& flags )
+{
+	clearDotTables( flags );
+	int size;
+	size = ( res*res + res )>>1;
+	if( flags & DOT_FLAG )
+	{
+		dotTable = new Real[size];
+		memset( dotTable , 0 , sizeof(Real)*size );
+	}
+	if( flags & D_DOT_FLAG )
+	{
+		dDotTable = new Real[size];
+		memset( dDotTable , 0 , sizeof(Real)*size );
+	}
+	if( flags & D2_DOT_FLAG )
+	{
+		d2DotTable = new Real[size];
+		memset( d2DotTable , 0 , sizeof(Real)*size );
+	}
+	double t1 , t2;
+	t1 = baseFunction.polys[0].start;
+	t2 = baseFunction.polys[baseFunction.polyCount-1].start;
+	for( int i=0 ; i<res ; i++ )
+	{
+		double c1 , c2 , w1 , w2;
+		BinaryNode<double>::CenterAndWidth( i , c1 , w1 );
+#if BOUNDARY_CONDITIONS
+		int d1 , d2 , off1 , off2;
+		BinaryNode< double >::DepthAndOffset( i , d1 , off1 );
+		int boundary1 = 0;
+		if     ( reflectBoundary && off1==0             ) boundary1 = -1;
+		else if( reflectBoundary && off1==( (1<<d1)-1 ) ) boundary1 =  1;
+#endif // BOUNDARY_CONDITIONS
+		double start1 = t1 * w1 + c1;
+		double end1   = t2 * w1 + c1;
+		for( int j=0 ; j<=i ; j++ )
+		{
+			BinaryNode<double>::CenterAndWidth( j , c2 , w2 );
+#if BOUNDARY_CONDITIONS
+			BinaryNode< double >::DepthAndOffset( j , d2 , off2 );
+			int boundary2 = 0;
+			if     ( reflectBoundary && off2==0             ) boundary2 = -1;
+			else if( reflectBoundary && off2==( (1<<d2)-1 ) ) boundary2 =  1;
+#endif // BOUNDARY_CONDITIONS
+			int idx = SymmetricIndex( i , j );
+
+			double start = t1 * w2 + c2;
+			double end   = t2 * w2 + c2;
+#if BOUNDARY_CONDITIONS
+			if( reflectBoundary )
+			{
+				if( start<0 ) start = 0;
+				if( start>1 ) start = 1;
+				if( end  <0 )   end = 0;
+				if( end  >1 )   end = 1;
+			}
+#endif // BOUNDARY_CONDITIONS
+
+			if( start<  start1 ) start = start1;
+			if( end  >  end1   )   end = end1;
+			if( start>= end    ) continue;
+
+#if BOUNDARY_CONDITIONS
+			Real dot = dotProduct( c1 , w1 , c2 , w2 , boundary1 , boundary2 );
+#else // !BOUNDARY_CONDITIONS
+			Real dot = dotProduct( c1 , w1 , c2 , w2 );
+#endif // BOUNDARY_CONDITIONS
+			if( fabs(dot)<1e-15 ) continue;
+			if( flags & DOT_FLAG ) dotTable[idx]=dot;
+			if( useDotRatios )
+			{
+#if BOUNDARY_CONDITIONS
+				if( flags &  D_DOT_FLAG )  dDotTable[idx] = -dDotProduct( c1 , w1 , c2 , w2 , boundary1 , boundary2 ) / dot;
+				if( flags & D2_DOT_FLAG ) d2DotTable[idx] = d2DotProduct( c1 , w1 , c2 , w2 , boundary1 , boundary2 ) / dot;
+#else // !BOUNDARY_CONDITIONS
+				if( flags &  D_DOT_FLAG )  dDotTable[idx] = -dDotProduct(c1,w1,c2,w2)/dot;
+				if( flags & D2_DOT_FLAG ) d2DotTable[idx] = d2DotProduct(c1,w1,c2,w2)/dot;
+#endif // BOUNDARY_CONDITIONS
+			}
+			else
+			{
+#if BOUNDARY_CONDITIONS
+				if( flags &  D_DOT_FLAG )  dDotTable[idx] =  dDotProduct( c1 , w1 , c2 , w2 , boundary1 , boundary2 );
+				if( flags & D2_DOT_FLAG ) d2DotTable[idx] = d2DotProduct( c1 , w1 , c2 , w2 , boundary1 , boundary2 );
+#else // !BOUNDARY_CONDTIONS
+				if( flags &  D_DOT_FLAG )  dDotTable[idx] =  dDotProduct(c1,w1,c2,w2);
+				if( flags & D2_DOT_FLAG ) d2DotTable[idx] = d2DotProduct(c1,w1,c2,w2);
+#endif // BOUNDARY_CONDITIONS
+			}
+		}
+	}
+}
+template<int Degree,class Real>
+void FunctionData<Degree,Real>::clearDotTables( const int& flags )
+{
+	if((flags & DOT_FLAG) && dotTable)
+	{
+		delete[] dotTable;
+		dotTable=NULL;
+	}
+	if((flags & D_DOT_FLAG) && dDotTable)
+	{
+		delete[] dDotTable;
+		dDotTable=NULL;
+	}
+	if((flags & D2_DOT_FLAG) && d2DotTable)
+	{
+		delete[] d2DotTable;
+		d2DotTable=NULL;
+	}
+}
+template<int Degree,class Real>
+void FunctionData<Degree,Real>::setValueTables( const int& flags , const double& smooth )
+{
+	clearValueTables();
+	if( flags &   VALUE_FLAG )  valueTables = new Real[res*res2];
+	if( flags & D_VALUE_FLAG ) dValueTables = new Real[res*res2];
+	PPolynomial<Degree+1> function;
+	PPolynomial<Degree>  dFunction;
+	for( int i=0 ; i<res ; i++ )
+	{
+		if(smooth>0)
+		{
+			function=baseFunctions[i].MovingAverage(smooth);
+			dFunction=baseFunctions[i].derivative().MovingAverage(smooth);
+		}
+		else
+		{
+			function=baseFunctions[i];
+			dFunction=baseFunctions[i].derivative();
+		}
+		for( int j=0 ; j<res2 ; j++ )
+		{
+			double x=double(j)/(res2-1);
+			if(flags &   VALUE_FLAG){ valueTables[j*res+i]=Real( function(x));}
+			if(flags & D_VALUE_FLAG){dValueTables[j*res+i]=Real(dFunction(x));}
+		}
+	}
+}
+template<int Degree,class Real>
+void FunctionData<Degree,Real>::setValueTables(const int& flags,const double& valueSmooth,const double& normalSmooth){
+	clearValueTables();
+	if(flags &   VALUE_FLAG){ valueTables=new Real[res*res2];}
+	if(flags & D_VALUE_FLAG){dValueTables=new Real[res*res2];}
+	PPolynomial<Degree+1> function;
+	PPolynomial<Degree>  dFunction;
+	for(int i=0;i<res;i++){
+		if(valueSmooth>0)	{ function=baseFunctions[i].MovingAverage(valueSmooth);}
+		else				{ function=baseFunctions[i];}
+		if(normalSmooth>0)	{dFunction=baseFunctions[i].derivative().MovingAverage(normalSmooth);}
+		else				{dFunction=baseFunctions[i].derivative();}
+
+		for(int j=0;j<res2;j++){
+			double x=double(j)/(res2-1);
+			if(flags &   VALUE_FLAG){ valueTables[j*res+i]=Real( function(x));}
+			if(flags & D_VALUE_FLAG){dValueTables[j*res+i]=Real(dFunction(x));}
+		}
+	}
+}
+
+
+template<int Degree,class Real>
+void FunctionData<Degree,Real>::clearValueTables(void){
+	if( valueTables){delete[]  valueTables;}
+	if(dValueTables){delete[] dValueTables;}
+	valueTables=dValueTables=NULL;
+}
+
+#if BOUNDARY_CONDITIONS
+template<int Degree,class Real>
+Real FunctionData<Degree,Real>::dotProduct( const double& center1 , const double& width1 , const double& center2 , const double& width2 , int boundary1 , int boundary2 ) const
+{
+	const PPolynomial< Degree > *b1 , *b2;
+	if     ( boundary1==-1 ) b1 = & leftBaseFunction;
+	else if( boundary1== 0 ) b1 = &     baseFunction;
+	else if( boundary1== 1 ) b1 = &rightBaseFunction;
+	if     ( boundary2==-1 ) b2 = & leftBaseFunction;
+	else if( boundary2== 0 ) b2 = &     baseFunction;
+	else if( boundary2== 1 ) b2 = &rightBaseFunction;
+	double r=fabs( baseFunction.polys[0].start );
+	switch( normalize )
+	{
+	case 2:
+		return Real(((*b1)*b2->scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r)*width1/sqrt(width1*width2));
+	case 1:
+		return Real(((*b1)*b2->scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r)*width1/(width1*width2));
+	default:
+		return Real(((*b1)*b2->scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r)*width1);
+	}
+}
+template<int Degree,class Real>
+Real FunctionData<Degree,Real>::dDotProduct( const double& center1 , const double& width1 , const double& center2 , const double& width2 , int boundary1 , int boundary2 ) const
+{
+	const PPolynomial< Degree-1 > *b1;
+	const PPolynomial< Degree   > *b2;
+	if     ( boundary1==-1 ) b1 = & dLeftBaseFunction;
+	else if( boundary1== 0 ) b1 = &     dBaseFunction;
+	else if( boundary1== 1 ) b1 = &dRightBaseFunction;
+	if     ( boundary2==-1 ) b2 = &  leftBaseFunction;
+	else if( boundary2== 0 ) b2 = &      baseFunction;
+	else if( boundary2== 1 ) b2 = & rightBaseFunction;
+	double r=fabs(baseFunction.polys[0].start);
+	switch(normalize){
+		case 2:
+			return Real(((*b1)*b2->scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r)/sqrt(width1*width2));
+		case 1:
+			return Real(((*b1)*b2->scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r)/(width1*width2));
+		default:
+			return Real(((*b1)*b2->scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r));
+	}
+}
+template<int Degree,class Real>
+Real FunctionData<Degree,Real>::d2DotProduct( const double& center1 , const double& width1 , const double& center2 , const double& width2 , int boundary1 , int boundary2 ) const
+{
+	const PPolynomial< Degree-1 > *b1 , *b2;
+	if     ( boundary1==-1 ) b1 = & dLeftBaseFunction;
+	else if( boundary1== 0 ) b1 = &     dBaseFunction;
+	else if( boundary1== 1 ) b1 = &dRightBaseFunction;
+	if     ( boundary2==-1 ) b2 = & dLeftBaseFunction;
+	else if( boundary2== 0 ) b2 = &     dBaseFunction;
+	else if( boundary2== 1 ) b2 = &dRightBaseFunction;
+	double r=fabs(baseFunction.polys[0].start);
+	switch( normalize )
+	{
+		case 2:
+			return Real(((*b1)*b2->scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r)/width2/sqrt(width1*width2));
+		case 1:
+			return Real(((*b1)*b2->scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r)/width2/(width1*width2));
+		default:
+			return Real(((*b1)*b2->scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r)/width2);
+	}
+}
+#else // !BOUNDARY_CONDITIONS
+template<int Degree,class Real>
+Real FunctionData<Degree,Real>::dotProduct(const double& center1,const double& width1,const double& center2,const double& width2) const{
+	double r=fabs(baseFunction.polys[0].start);
+	switch( normalize )
+	{
+		case 2:
+			return Real((baseFunction*baseFunction.scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r)*width1/sqrt(width1*width2));
+		case 1:
+			return Real((baseFunction*baseFunction.scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r)*width1/(width1*width2));
+		default:
+			return Real((baseFunction*baseFunction.scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r)*width1);
+	}
+}
+template<int Degree,class Real>
+Real FunctionData<Degree,Real>::dDotProduct(const double& center1,const double& width1,const double& center2,const double& width2) const{
+	double r=fabs(baseFunction.polys[0].start);
+	switch(normalize){
+		case 2:
+			return Real((dBaseFunction*baseFunction.scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r)/sqrt(width1*width2));
+		case 1:
+			return Real((dBaseFunction*baseFunction.scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r)/(width1*width2));
+		default:
+			return Real((dBaseFunction*baseFunction.scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r));
+	}
+}
+template<int Degree,class Real>
+Real FunctionData<Degree,Real>::d2DotProduct(const double& center1,const double& width1,const double& center2,const double& width2) const{
+	double r=fabs(baseFunction.polys[0].start);
+	switch(normalize){
+		case 2:
+			return Real((dBaseFunction*dBaseFunction.scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r)/width2/sqrt(width1*width2));
+		case 1:
+			return Real((dBaseFunction*dBaseFunction.scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r)/width2/(width1*width2));
+		default:
+			return Real((dBaseFunction*dBaseFunction.scale(width2/width1).shift((center2-center1)/width1)).integral(-2*r,2*r)/width2);
+	}
+}
+#endif // BOUNDARY_CONDITIONS
+template<int Degree,class Real>
+inline int FunctionData<Degree,Real>::SymmetricIndex( const int& i1 , const int& i2 )
+{
+	if( i1>i2 ) return ((i1*i1+i1)>>1)+i2;
+	else        return ((i2*i2+i2)>>1)+i1;
+}
+template<int Degree,class Real>
+inline int FunctionData<Degree,Real>::SymmetricIndex( const int& i1 , const int& i2 , int& index )
+{
+	if( i1<i2 )
+	{
+		index = ((i2*i2+i2)>>1)+i1;
+		return 1;
+	}
+	else{
+		index = ((i1*i1+i1)>>1)+i2;
+		return 0;
+	}
+}
diff --git a/src/plugins_experimental/filter_screened_poisson/Src/Hash.h b/src/plugins_experimental/filter_screened_poisson/Src/Hash.h
new file mode 100755
index 000000000..b66c07362
--- /dev/null
+++ b/src/plugins_experimental/filter_screened_poisson/Src/Hash.h
@@ -0,0 +1,29 @@
+#ifndef HASH_INCLUDED
+#define HASH_INCLUDED
+#ifdef WIN32
+#include <hash_map>
+using stdext::hash_map;
+#else // !WIN32
+#include <ext/hash_map>
+using namespace __gnu_cxx;
+
+namespace __gnu_cxx
+{
+  template<> struct hash<long long> {
+    size_t operator()(long long __x) const { return __x; }
+  };
+  template<> struct hash<const long long> {
+    size_t operator()(const long long __x) const { return __x; }
+  };
+  
+  
+  template<> struct hash<unsigned long long> {
+    size_t operator()(unsigned long long __x) const { return __x; }
+  };
+  template<> struct hash<const unsigned long long> {
+    size_t operator()(const unsigned long long __x) const { return __x; }
+  };
+}
+#endif // WIN32
+#endif // HASH_INCLUDED
+
diff --git a/src/plugins_experimental/filter_screened_poisson/Src/Time.cpp b/src/plugins_experimental/filter_screened_poisson/Src/Time.cpp
new file mode 100755
index 000000000..b6b6869bc
--- /dev/null
+++ b/src/plugins_experimental/filter_screened_poisson/Src/Time.cpp
@@ -0,0 +1,46 @@
+/*
+Copyright (c) 2006, Michael Kazhdan and Matthew Bolitho
+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 the Johns Hopkins University nor the names of its contributors
+may 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.
+*/
+
+#include <string.h>
+#include <sys/timeb.h>
+#ifndef WIN32
+#include <sys/time.h>
+#endif // WIN32
+
+double Time( void )
+{
+#ifdef WIN32
+	struct _timeb t;
+	_ftime( &t );
+	return double( t.time ) + double( t.millitm ) / 1000.0;
+#else // WIN32
+	struct timeval t;
+	gettimeofday( &t , NULL );
+	return t.tv_sec + double( t.tv_usec ) / 1000000;
+#endif // WIN32
+}
diff --git a/src/plugins_experimental/filter_screened_poisson/Src/Time.h b/src/plugins_experimental/filter_screened_poisson/Src/Time.h
new file mode 100755
index 000000000..2b3cd59f3
--- /dev/null
+++ b/src/plugins_experimental/filter_screened_poisson/Src/Time.h
@@ -0,0 +1,32 @@
+/*
+Copyright (c) 2006, Michael Kazhdan and Matthew Bolitho
+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 the Johns Hopkins University nor the names of its contributors
+may 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.
+*/
+
+#ifndef TIME_INCLUDED
+#define TIME_INCLUDED
+double Time(void);
+#endif // TIME_INCLUDED