Imagine++
Algos.h
1// ===========================================================================
2// Imagine++ Libraries
3// Copyright (C) Imagine
4// For detailed information: http://imagine.enpc.fr/software
5// ===========================================================================
6
7namespace Imagine {
10
11 // ===============================================
12 // Scalings
13
22 template <typename T, int dim>
24 {
25 Image<T,dim> I1(I.sizes() * fact);
26 // coordinates aware iterator
27 for (CoordsIterator<dim> r = I.coordsBegin() ; r != I.coordsEnd() ; ++r) {
28 Coords<dim> o=*r;
29 T b=I(o);
30 for (CoordsIterator<dim> zr(o * fact, o * fact + Coords<dim>(fact-1));zr!=I.coordsEnd();++zr){
31 I1(*zr)=b;
32 }
33 }
34 return I1;
35 }
44 template <typename T, int dim>
46 {
47 Image<T,dim> I1(I.sizes()/fact);
48 // coordinates aware iterator
49 for (CoordsIterator<dim> r = I1.coordsBegin() ; r != I1.coordsEnd() ; ++r)
50 I1(*r)=I((*r)*fact);
51 return I1;
52 }
61 template <typename T, int dim>
62 Image<T,dim> reduce(const Image<T,dim>&I, int fact)
63 {
64 typedef typename PixelTraits<T>::template CastPixel<double>::value_type doubleT;
65 Coords<dim> d=I.sizes()/fact;
66 double nb=pow(double(fact),dim);
67 Image<T,dim> Ir(d);
68 for (CoordsIterator<dim> r = Ir.coordsBegin() ; r != Ir.coordsEnd() ; ++r) {
69 doubleT p(0.);
70 for (CoordsIterator<dim> r2((*r)*fact,((*r)*fact+Coords<dim>(fact-1))); r2!=CoordsIterator<dim>(); ++r2) {
71 p+=doubleT(I(*r2));
72 }
73 Ir(*r)=T(p/nb);
74 }
75 return Ir;
76 }
87 template <typename T, int dim>
88 Image<T,dim> reduce(const Image<T,dim>&I,Coords<dim> nd,bool keepRatio=false)
89 {
90 typedef typename PixelTraits<T>::template CastPixel<double>::value_type doubleT;
92 Coords<dim> od=I.sizes();
94 std::pair<double,double> mM = range(f);
95 assert(mM.first>=1);
96 if (keepRatio) {
97 f.fill(mM.second);
98 nd=Coords<dim>(FVector<double,dim>(od)/mM.second);
99 }
100 inPlaceBlur(oI,1.5*(sqrt(f)-.99)); // Todo: better filtering
101 Image<T,dim> nI(nd);
102 for (CoordsIterator<dim> r = nI.coordsBegin() ; r != nI.coordsEnd() ; ++r) {
104 nI(*r)=T(oI.interpolate(x));
105 }
106 return nI;
107 }
109 template <typename T> Image<T,2> inline reduce(const Image<T,2>&I,int w,int h,bool keepRatio=false) { return reduce(I,Coords<2>(w,h),keepRatio); } // Reduce image (given dimensions), 3D alias.
111 template <typename T> Image<T,3> inline reduce(const Image<T,3>&I,int w,int h,int d,bool keepRatio=false) { return reduce(I,Coords<3>(w,h,d),keepRatio); }
120 template <typename T,int dim>
121 inline Image<T,dim> reduce(const Image<T,dim>&I,double fact)
122 {
124 return reduce(I,nd);
125 }
136 template <typename T,int dim>
137 Image<T,dim> enlarge(const Image<T,dim>&I,Coords<dim> nd,bool keepRatio=false)
138 {
139 Coords<dim> od=I.sizes();
141 std::pair<double,double> mM = range(f);
142 assert(mM.second<=1);
143 if (keepRatio) {
144 f.fill(mM.second);
145 nd=Coords<dim>(FVector<double,dim>(od)/mM.second);
146 }
147 Image<T,dim> nI(nd);
148 for (CoordsIterator<dim> r = nI.coordsBegin() ; r != nI.coordsEnd() ; ++r) {
150 nI(*r)=T(I.interpolate(x));
151 }
152 return nI;
153 }
155 template <typename T> Image<T,2> inline enlarge(const Image<T,2>&I,int w,int h,bool keepRatio=false) { return enlarge(I,Coords<2>(w,h),keepRatio); }
157 template <typename T> Image<T,3> inline enlarge(const Image<T,3>&I,int w,int h,int d,bool keepRatio=false) { return enlarge(I,Coords<3>(w,h,d),keepRatio); }
166 template <typename T,int dim>
167 inline Image<T,dim> enlarge(const Image<T,dim>&I,double fact)
168 {
170 return enlarge(I,nd);
171 }
172
173 // ===============================================
174 // Deriche
175
187 template <typename T,int dim>
188 void inPlaceDeriche(Image<T,dim>&I,typename PixelTraits<T>::scalar_type sigma, int order, int d, bool neumann = true) {
189 // Checks parameter values
190 assert(sigma>0 && order>=0 && order<3 && d>=0 && d<dim);
191
192 // Computes coefficients of the recursive filter
193 const typename PixelTraits<T>::scalar_type
194 alpha = 1.695f/sigma,
195 ea = std::exp(alpha),
196 ema = std::exp(-alpha),
197 em2a = ema*ema,
198 b1 = 2*ema,
199 b2 = -em2a;
200
201 typename PixelTraits<T>::scalar_type ek,ekn,parity,a1,a2,a3,a4,g0,sumg1,sumg0;
202
203 switch(order) {
204
205 // first-order derivative
206 case 1:
207 ek = -(1-ema)*(1-ema)*(1-ema)/(2*(ema+1)*ema);
208 a1 = a4 = 0;
209 a2 = ek*ema;
210 a3 = -ek*ema;
211 parity = -1;
212 if (neumann) {
213 sumg1 = (ek*ea) / ((ea-1)*(ea-1));
214 g0 = 0;
215 sumg0 = g0 + sumg1;
216 }
217 else
218 g0 = sumg0 = sumg1 = 0;
219 break;
220
221 // second-order derivative
222 case 2:
223 ekn = ( -2*(-1+3*ea-3*ea*ea+ea*ea*ea)/(3*ea+1+3*ea*ea+ea*ea*ea) );
224 ek = -(em2a-1)/(2*alpha*ema);
225 a1 = ekn;
226 a2 = -ekn*(1+ek*alpha)*ema;
227 a3 = ekn*(1-ek*alpha)*ema;
228 a4 = -ekn*em2a;
229 parity = 1;
230 if (neumann) {
231 sumg1 = ekn/2;
232 g0 = ekn;
233 sumg0 = g0 + sumg1;
234 }
235 else
236 g0 = sumg0 = sumg1 = 0;
237 break;
238
239 // smoothing
240 default:
241 ek = (1-ema)*(1-ema) / (1+2*alpha*ema - em2a);
242 a1 = ek;
243 a2 = ek*ema*(alpha-1);
244 a3 = ek*ema*(alpha+1);
245 a4 = -ek*em2a;
246 parity = 1;
247 if (neumann) {
248 sumg1 = ek*(alpha*ea+ea-1) / ((ea-1)*(ea-1));
249 g0 = ek;
250 sumg0 = g0 + sumg1;
251 }
252 else
253 g0 = sumg0 = sumg1 = 0;
254 break;
255 }
256
257 // filter init
258 T *Y = new T[I.size(d)];
259 const size_t offset = I.stride(d);
260 const size_t nb = I.size(d);
261
262 // Iterates on dimensions other than d
263 Coords<dim> beg(0), end = I.sizes() - Coords<dim>(1);
264 end[d]=0;
265 for (CoordsIterator<dim> p(beg,end); p != CoordsIterator<dim>(); ++p) {
266 T *ima = &(I(*p));
267 T I2 = *ima; ima += offset;
268 T I1 = *ima; ima += offset;
269 T Y2 = *(Y++) = sumg0*I2;
270 T Y1 = *(Y++) = g0*I1 + sumg1*I2;
271 for (size_t i=2; i<nb; i++) {
272 I1 = *ima; ima+=offset;
273 T Y0 = *(Y++) = a1*I1 + a2*I2 + b1*Y1 + b2*Y2;
274 I2=I1; Y2=Y1; Y1=Y0;
275 }
276 ima -= offset;
277 I2 = *ima;
278 Y2 = Y1 = (parity*sumg1)*I2;
279 *ima = *(--Y)+Y2;
280 ima-=offset;
281 I1 = *ima;
282 *ima = *(--Y)+Y1;
283 for (size_t i=nb-3; ; i--) {
284 T Y0 = a3*I1+a4*I2+b1*Y1+b2*Y2;
285 ima-=offset;
286 I2=I1;
287 I1=*ima;
288 *ima=*(--Y)+Y0;
289 Y2=Y1;
290 Y1=Y0;
291 if (i==0)
292 break;
293 }
294 }
295 delete [] Y;
296
297 }
310 template <typename T,int dim>
311 Image<T,dim> deriche(const Image<T,dim>&I,typename PixelTraits<T>::scalar_type sigma, int order, int d, bool neumann = true) {
312 Image<T,dim> J=I.clone();
313 inPlaceDeriche(J,sigma,order,d,neumann);
314 return J;
315 }
324 template <typename T,int dim>
325 void inPlaceBlur(Image<T,dim>&I,const FVector<typename PixelTraits<T>::scalar_type,dim>& sigmas, bool neumann = true) {
326 for (int i=0;i<dim;i++) {
327 inPlaceDeriche(I,sigmas[i], 0, i, neumann);
328 }
329 }
338 template <typename T,int dim>
339 void inPlaceBlur(Image<T,dim>&I,typename PixelTraits<T>::scalar_type sigma, bool neumann = true) {
340 inPlaceBlur(I,FVector<typename PixelTraits<T>::scalar_type,dim>(sigma),neumann);
341 }
351 template <typename T,int dim>
352 Image<T,dim> blur(const Image<T,dim>&I,typename PixelTraits<T>::scalar_type sigma, bool neumann = true) {
353 Image<T,dim> J=I.clone();
354 inPlaceBlur(J,sigma,neumann);
355 return J;
356 }
366 template <typename T,int dim>
367 Image<T,dim> blur(const Image<T,dim>&I,const FVector<typename PixelTraits<T>::scalar_type,dim>& sigmas, bool neumann = true) {
368 Image<T,dim> J=I.clone();
369 inPlaceBlur(J,sigmas,neumann);
370 return J;
371 }
372
374}
Coordinates.
Definition: Coords.h:16
Iterator on Coordinates.
Definition: Coords.h:83
Vector of fixed size.
Definition: FVector.h:17
FVector & fill(const T &v)
Filling.
Definition: FVector.h:75
Image.
Definition: Image.h:24
Image clone() const
Cloning.
Definition: Image.h:146
PixelTraits< T >::template CastPixel< V >::value_type interpolate(const FVector< V, dim > &c) const
Interpolation.
Definition: Image.h:198
FArray< size_t, dim > stride() const
Stride.
Definition: MultiArray.h:253
CoordsIterator< dim > coordsBegin() const
Begin coords iterator.
Definition: MultiArray.h:376
CoordsIterator< dim > coordsEnd() const
End coords iterator.
Definition: MultiArray.h:384
int size(int i) const
ith size.
Definition: MultiArray.h:225
Coords< dim > sizes() const
Sizes.
Definition: MultiArray.h:211
Image< T, dim > enlarge(const Image< T, dim > &I, Coords< dim > nd, bool keepRatio=false)
Enlarge image (given dimensions).
Definition: Algos.h:137
Image< T, dim > scaleDown(const Image< T, dim > &I, int fact)
Down scaling: fast naive version.
Definition: Algos.h:45
Image< T, dim > reduce(const Image< T, dim > &I, int fact)
Reduce image (integer factor).
Definition: Algos.h:62
Image< T, dim > scaleUp(const Image< T, dim > &I, int fact)
Up scaling: fast naive version.
Definition: Algos.h:23
Image< T, dim > blur(const Image< T, dim > &I, typename PixelTraits< T >::scalar_type sigma, bool neumann=true)
Blur.
Definition: Algos.h:352
void inPlaceBlur(Image< T, dim > &I, const FVector< typename PixelTraits< T >::scalar_type, dim > &sigmas, bool neumann=true)
In Place Blur (anisotropic).
Definition: Algos.h:325
Image< T, dim > deriche(const Image< T, dim > &I, typename PixelTraits< T >::scalar_type sigma, int order, int d, bool neumann=true)
Deriche filter.
Definition: Algos.h:311
void inPlaceDeriche(Image< T, dim > &I, typename PixelTraits< T >::scalar_type sigma, int order, int d, bool neumann=true)
In place Deriche filter.
Definition: Algos.h:188
Imagine++ namespace.