NDEVR
API Documentation
FourierPoisson.h
1#include <cmath>
2#include <iostream>
3#include <DGtal/base/Common.h>
4#include <DGtal/images/CImage.h>
5#include <DGtal/io/Color.h>
6#include <DGtal/io/readers/GenericReader.h>
7#include <DGtal/io/writers/PPMWriter.h>
8#include <DGtal/io/writers/PGMWriter.h>
9//#include <DGtal/math/RealFFT.h>
10
11namespace DGtal {
12
14 namespace image {
16 struct ColorToRedFunctor {
17 int operator()( int color ) const
18 { return (color >> 16) & 0xff; }
19 };
21 struct ColorToGreenFunctor {
22 int operator()( int color ) const
23 { return (color >> 8) & 0xff; }
24 };
26 struct ColorToBlueFunctor {
27 int operator()( int color ) const
28 { return color & 0xff; }
29 };
31 struct ColorToGrayFunctor {
32 int operator()( int color ) const
33 {
34 double val = 0.2126 * ( (double) ( ( color >> 16 ) & 0xff ) )
35 + 0.7152 * ( (double) ( ( color >> 8 ) & 0xff ) )
36 + 0.0722 * ( (double) ( ( color >> 0 ) & 0xff ) );
37 return ( (int) round( val ) ) & 0xff;
38 }
39 };
41 struct DoubleToGrayFunctor {
42 unsigned char operator()( double scalar ) const
43 { return (unsigned char) std::min( 255.0, std::max( 0.0, round( scalar ) ) ); }
44 };
46 struct SignedDoubleToGrayFunctor {
47 unsigned char operator()( double scalar ) const
48 { return (unsigned char) std::min( 255.0, std::max( 0.0, 128.0 + round( scalar ) ) ); }
49 };
51 struct GrayToGrayFunctor {
52 int operator()( int color ) const
53 { return color & 0xff; }
54 };
56 struct GrayToRedGreen {
57 DGtal::Color operator()( int value ) const
58 {
59 if ( value >= 128 ) return DGtal::Color( 0, (value-128)*2+1, 0 );
60 else return DGtal::Color( 0, 0, 255-value*2 );
61 }
62 };
63
65 struct Zoom {
66 typedef Z2i::Integer Integer;
67 typedef Z2i::Point Point;
68 typedef Z2i::RealPoint RealPoint;
69 typedef Z2i::Domain Domain;
70 typedef RealPoint::Component Scalar;
71
72 Zoom( Domain aDomain, int aZoom )
73 : _domain( aDomain ), _zoom( aZoom ) {
74 _zDomain = getZoomedDomain();
75 _stride = _domain.upperBound()[ 0 ] - _domain.lowerBound()[ 0 ] + 1;
76 _zStride = _zDomain.upperBound()[ 0 ] - _zDomain.lowerBound()[ 0 ] + 1;
77 }
78
79 const Domain& domain() const { return _domain; }
80 const Domain& zoomedDomain() const { return _zDomain; }
81 Domain getZoomedDomain() const {
82 return Domain( zoom( _domain.lowerBound() ), zoom( _domain.upperBound() ) );
83 }
84
85 Integer index( Point p ) const
86 {
87 return p[ 1 ] * _stride + p[ 0 ];
88 }
89
90 Integer zoomedIndex( Point p ) const
91 {
92 return p[ 1 ] * _zStride + p[ 0 ];
93 }
94
95 Point zoom( Point p ) const
96 {
97 return p * _zoom;
98 }
99
100 RealPoint zoom( RealPoint p ) const
101 {
102 return p * ( Scalar ) _zoom;
103 }
104
105 Point unzoom( Point p ) const
106 {
107 return p / _zoom;
108 }
109
110 RealPoint real( Point p ) const
111 { return RealPoint( p[ 0 ], p[ 1 ] ); }
112
113 Point closest( RealPoint p ) const
114 { return Point( (Integer) round( p[ 0 ] ), (Integer) round( p[ 1 ] ) ); }
115
116 RealPoint project( Point p ) const
117 {
118 return real( p ) / (Scalar) _zoom;
119 }
120
121 RealPoint unzoom( RealPoint p ) const
122 {
123 return p / ( Scalar ) _zoom;
124 }
125
126 Domain _domain;
127 Domain _zDomain;
128 Integer _zoom;
129 Integer _stride;
130 Integer _zStride;
131 };
132
134 struct Utils {
135 typedef double Scalar;
136 typedef PointVector< 3, Scalar > Value;
137 // typedef std::array< Value, 2 > VectorValue;
138 struct VectorValue {
139 Value x;
140 Value y;
141 };
142 typedef std::vector<Scalar> ScalarForm;
143 typedef std::vector<Value> ValueForm;
144 typedef std::vector<VectorValue> VectorValueForm;
145
146 template <typename Image>
147 ScalarForm getScalarForm( const Image& image, bool color )
148 {
149 unsigned int v = 0;
150 ScalarForm I( image.size() );
151 if ( color ) {
152 auto F = ColorToGrayFunctor();
153 for ( unsigned int val : image )
154 I[ v++ ] = (Scalar) F( val );
155 } else {
156 auto F = GrayToGrayFunctor();
157 for ( unsigned int val : image )
158 I[ v++ ] = (Scalar) F( val );
159 }
160 return I;
161 }
162 template <typename Image>
163 ValueForm getValueForm( const Image& image, bool color )
164 {
165 unsigned int v = 0;
166 ValueForm I( image.size() );
167 if ( color ) {
168 auto R = ColorToRedFunctor();
169 auto G = ColorToGreenFunctor();
170 auto B = ColorToBlueFunctor();
171 for ( unsigned int val : image )
172 I[ v++ ] = Value( (Scalar) R( val ), (Scalar) G( val ), (Scalar) B( val ) );
173 } else {
174 auto F = GrayToGrayFunctor();
175 for ( unsigned int val : image )
176 I[ v++ ] = Value( (Scalar) F( val ), (Scalar) F( val ), (Scalar) F( val ) );
177 }
178 return I;
179 }
180 };
181
182 // ----------------------------------------------------------------------
183 namespace functions {
184 template <typename TImageDst, typename TImageSrc,
185 typename TFunctor>
186 void transform( TImageDst& dst, const TImageSrc& src,
187 TFunctor F )
188 {
189 std::transform( src.begin(), src.end(), dst.begin(), F );
190 }
191
192 template <typename ScalarImage>
193 bool exportPGM( const char* name, const ScalarImage& image )
194 {
195 return PGMWriter<ScalarImage,DoubleToGrayFunctor>::exportPGM( name, image );
196 }
197 template <typename ScalarImage>
198 bool exportGradPGM( const char* name, const ScalarImage& image )
199 {
200 return PGMWriter<ScalarImage,SignedDoubleToGrayFunctor>::exportPGM( name, image );
201 }
202
203 // template <typename TImage>
204 // void bandFilter( TImage& image,
205 // double low_freq = 0.0, double high_freq = 0.5,
206 // double shift = 0.0,
207 // double output_min = 0.0,
208 // double output_max = 255.0 )
209 // {
210 // typedef TImage Image;
211 // BOOST_CONCEPT_ASSERT(( DGtal::concepts::CImage< Image > ));
212
213 // typedef typename Image::Domain Domain;
214 // typedef typename Domain::Space Space;
215 // typedef typename Space::Point Point;
216 // typedef typename Space::RealVector RealVector;
217 // typedef typename RealVector::Component Scalar;
218 // typedef typename Image::Value Value;
219
220 // Domain domain( image.domain() );
221 // RealFFT< Domain, Scalar > F ( domain );
222 // auto IF = F.getSpatialImage();
223 // auto FF = F.getFreqImage();
224 // // trace.info() << "spatial domain = " << F.getSpatialDomain() << std::endl;
225 // // trace.info() << "frequency domain = " << F.getFreqDomain() << std::endl;
226 // auto IF_it = IF.begin();
227 // for ( auto I_it = image.begin(), I_itE = image.end(); I_it != I_itE; I_it++ )
228 // *IF_it++ = (Scalar) *I_it;
229 // F.forwardFFT();
230 // Scalar l2 = (Scalar) ( low_freq * low_freq );
231 // Scalar h2 = (Scalar) ( high_freq * high_freq );
232 // for ( auto&& fp : F.getFreqDomain() )
233 // {
234 // auto freq = F.calcScaledFreqCoords( fp );
235 // auto nfreq = freq.dot( freq );
236 // if ( nfreq < l2 || nfreq > h2 )
237 // FF.setValue( fp, (Scalar) 0 );
238 // }
239 // F.backwardFFT();
240 // auto I_it = image.begin();
241 // for ( auto IF_it = IF.begin(), IF_itE = IF.end(); IF_it != IF_itE; IF_it++ ) {
242 // Value val = (Value) std::max( output_min,
243 // std::min( output_max, shift + (*IF_it) ) );
244 // *I_it++ = val;
245 // }
246 // }
247
248 // template <typename ScalarImage>
249 // void poisson( ScalarImage& I, const ScalarImage& Gx, const ScalarImage& Gy )
250 // {
251 // typedef typename ScalarImage::Point Point;
252 // typedef typename ScalarImage::Domain Domain;
253 // typedef typename ScalarImage::Value Scalar;
254 // Domain domain( I.domain() );
255 // Domain sym_domain( Domain( -domain.upperBound() - Point::diagonal( 1 ),
256 // domain.upperBound() ) );
257 // RealFFT< Domain, double > U ( sym_domain );
258 // // Calcule l'image gradient.
259 // RealFFT< Domain, double > Vx( sym_domain );
260 // RealFFT< Domain, double > Vy( sym_domain );
261 // auto IU = U.getSpatialImage();
262 // auto IVx = Vx.getSpatialImage();
263 // auto IVy = Vy.getSpatialImage();
264 // for ( auto&& p : domain )
265 // {
266 // // Extend by symetry
267 // Point sym_px ( -1 - p[ 0 ], p[ 1 ] );
268 // Point sym_py ( p[ 0 ], -1 - p[ 1 ] );
269 // Point sym_pxy( -1 - p[ 0 ], -1 - p[ 1 ] );
270 // std::array<Point,4> all_p = { p, sym_px, sym_py, sym_pxy };
271 // double vI = (double) I( p );
272 // double vGx = (double) Gx( p );
273 // double vGy = (double) Gy( p );
274 // for ( auto&&q : all_p ) IU.setValue( q, vI );
275 // IVx.setValue( p, vGx );
276 // IVx.setValue( sym_px, -vGx );
277 // IVx.setValue( sym_py, vGx );
278 // IVx.setValue( sym_pxy, -vGx );
279 // IVy.setValue( p, vGy );
280 // IVy.setValue( sym_px, vGy );
281 // IVy.setValue( sym_py, -vGy );
282 // IVy.setValue( sym_pxy, -vGy );
283 // }
284 // trace.info() << "Compute FFT[V]" << std::endl;
285 // trace.info() << "spatial domain = " << Vx.getSpatialDomain() << std::endl;
286 // trace.info() << "frequency domain = " << Vx.getFreqDomain() << std::endl;
287 // U.forwardFFT();
288 // Vx.forwardFFT();
289 // Vy.forwardFFT();
290 // Domain freq_domain = Vx.getFreqDomain();
291 // auto FU = U.getFreqImage();
292 // auto FVx = Vx.getFreqImage();
293 // auto FVy = Vy.getFreqImage();
294 // typedef typename RealFFT< Domain, double >::Complex Complex;
295 // const double two_pi = 2.0 * M_PI;
296 // const Complex i = Complex( 0.0, 1.0 );
297 // const double J = ( sym_domain.upperBound()[0] - sym_domain.lowerBound()[0] ) + 1;
298 // const double L = ( sym_domain.upperBound()[1] - sym_domain.lowerBound()[1] ) + 1;
299 // for ( auto&& p : freq_domain )
300 // {
301 // auto freq = U.calcScaledFreqCoords( p );
302 // const double mj = freq[ 0 ];
303 // const double nl = freq[ 1 ];
304 // Complex val = i * ( mj * FVx( p ) + nl * FVy( p ) )
305 // / ( two_pi * ( mj * mj + nl * nl ) );
306 // if ( p == Point::zero ) continue;
307 // else FU.setValue( p, -val );
308 // }
309 // U.backwardFFT();
310 // for ( auto&& p : domain )
311 // I.setValue( p, (Scalar) std::max( 0.0, std::min( 255.0, IU( p ) ) ) );
312 // } // poisson
313 } // namespace functions
314 } // namespace image
315} // namespace DGtal