00001
00002
00003
00004
00005
00006
00007
00008
00009
00010 #ifndef HxGaussIIR_h
00011 #define HxGaussIIR_h
00012
00013
00014 #include "HxImageRep.h"
00015 #include "HxImageGenerator.h"
00016
00017 #define MAXRECUR 5
00018 #define MAXDERI 2
00019
00020
00021 class L_HXIMAGEREP GaussIIR: public HxImageGenerator
00022 {
00023 public:
00024 GaussIIR(double sigma, int derivativeOrder = 0,
00025 int filterOrder = 3);
00026 virtual HxSizes domainSize() const
00027 {
00028 return HxSizes(2*_filterOrder+1, 1, 1);
00029 }
00030 virtual HxVec3Double get(int x, int y, int z)
00031 {
00032 return _filter[x];
00033 }
00034 double filterOrder() const
00035 {
00036 return _filterOrder;
00037 }
00038
00039 double sigma() const
00040 {
00041 return _sigma;
00042 }
00043
00044 double derivativeOrder() const
00045 {
00046 return _derivativeOrder;
00047 }
00048
00049 STD_OSTREAM& put(STD_OSTREAM&) const;
00050
00051
00052 template<class TypePtr> void forward3(TypePtr ptr, int n)
00053 {
00054 double val, p1, p2, p3;
00055 double a0 = _filter[3].x();
00056 double a1 = _filter[2].x();
00057 double a2 = _filter[1].x();
00058 double a3 = _filter[0].x();
00059
00060 p1 = *ptr++; p2 = p1; p3 = p2;
00061 n--;
00062 while (n--) {
00063 val = a0 * *ptr + a1 * p1 + a2 * p2 + a3 * p3;
00064 *ptr++ = val;
00065 p3 = p2; p2 = p1; p1 = val;
00066 }
00067 }
00068
00069 template<class TypePtr> void backward3(TypePtr ptr, int n)
00070 {
00071 double val, p1, p2, p3;
00072 double a0 = _filter[3].x();
00073 double a1 = _filter[4].x();
00074 double a2 = _filter[5].x();
00075 double a3 = _filter[6].x();
00076
00077 p1 = *(--ptr); p2 = p1; p3 = p2;
00078 n--;
00079 while (n--) {
00080 val = a0 * *(--ptr) + a1 * p1 + a2 * p2 + a3 * p3;
00081 *ptr = val;
00082 p3 = p2; p2 = p1; p1 = val;
00083 }
00084 }
00085
00086 #ifdef NOMSDEVBUG
00087 template<class TypePtr> void backward3(TypePtr ptr, int n, int deri)
00088 {
00089 double val, p1, p2, p3;
00090 double a0 = _filter[3].x();
00091 double a1 = _filter[4].x();
00092 double a2 = _filter[5].x();
00093 double a3 = _filter[6].x();
00094
00095 p1 = *(--ptr); p2 = p1; p3 = p2;
00096 n--;
00097 switch(deri) {
00098 case 1:
00099 while (n--) {
00100 val = a0 * *(--ptr) + a1 * p1 + a2 * p2 + a3 * p3;
00101 *(ptr+1) = p2-val;
00102 p3 = p2; p2 = p1; p1 = val;
00103 }
00104 break;
00105 case 2:
00106 while (n--) {
00107 val = a0 * *(--ptr) + a1 * p1 + a2 * p2 + a3 * p3;
00108 *(ptr+1) = val-p1-p1+p2;
00109 p3 = p2; p2 = p1; p1 = val;
00110 }
00111 break;
00112 default:
00113 while (n--) {
00114 val = a0 * *(--ptr) + a1 * p1 + a2 * p2 + a3 * p3;
00115 *ptr = val;
00116 p3 = p2; p2 = p1; p1 = val;
00117 }
00118 break;
00119 }
00120 }
00121 #endif
00122
00123 template<class InPtrT, class OutPtrT>
00124 void forward3(InPtrT in, OutPtrT out, int n)
00125 {
00126 double val, p1, p2, p3;
00127 double a0 = _filter[3].x();
00128 double a1 = _filter[2].x();
00129 double a2 = _filter[1].x();
00130 double a3 = _filter[0].x();
00131
00132 p1 = *in++; p2 = p1; p3 = p2;
00133 *out++ = p1;
00134 n--;
00135 while (n--) {
00136 val = a0 * *in++ + a1 * p1 + a2 * p2 + a3 * p3;
00137 *out++ = val;
00138 p3 = p2; p2 = p1; p1 = val;
00139 }
00140 }
00141
00142 template<class TypePtr> void lineFilter(TypePtr ptr, int n, int deri=0)
00143 {
00144 switch(_filterOrder) {
00145 case 3:
00146 forward3(ptr, n);
00147 backward3(ptr+n, n);
00148 break;
00149 }
00150 }
00151
00152 template<class InPtrT, class OutPtrT>
00153 void lineFilter(InPtrT in, OutPtrT out, int n, int deri=0)
00154 {
00155 switch(_filterOrder) {
00156 case 3:
00157 forward3(in, out, n);
00158 backward3(out+n, n);
00159 break;
00160 }
00161 }
00162
00163
00164 private:
00165 HxVec3Double _filter[2*MAXRECUR+1];
00166
00167 double _sigma;
00168 int _derivativeOrder;
00169 int _filterOrder;
00170 };
00171
00172
00173 inline STD_OSTREAM&
00174 operator<<(STD_OSTREAM& os, const GaussIIR& filter)
00175 {
00176 return filter.put(os);
00177 }
00178
00179 #endif