sharp/apps/fresize.c

#include <stdio.h>

#include <stdlib.h>

#include <time.h>

#include <string.h>

#include <math.h>

#include <mkl_blas.h>

#include <mkl_service.h>

#include <mkl_lapack.h>

#include <mkl_vml_functions.h>

#include <omp.h>

#include "fresize.h"

#define minval(x,y) (((x) < (y)) ? (x) : (y))

#define maxval(x,y) (((x) < (y)) ? (y) : (x))

#define fresize_sample 1

#define fresize_bin 2

#define fresize_1d 3

#define fresize_2d 4


double sinc(double x)

{

if (fabs(x) < (1.e-10))

  return 1.;

else

  return sin(M_PI*x)/(M_PI*x);

}


int init_fresize_sample(

  struct fresize_struct *pars,

  int nsub

)

{


  pars->method=fresize_sample;

  pars->nsub=nsub;


  return 0;

}


int init_fresize_bin(

  struct fresize_struct *pars,

  int nsub

)

{


  pars->method=fresize_bin;

  pars->nsub=nsub;


  return 0;

}


int init_fresize_boxcar(

  struct fresize_struct *pars,

  int hwidth, // Half width of boxcar. Full is 2*hwidth+1.

  int nsub // Distance between sampled points

)

{

  const int malign=32;

  int fwidth;

  int i,j;


  pars->method=fresize_1d;

  pars->nsub=nsub;

  pars->hwidth=hwidth;

  fwidth=2*hwidth+1;

  pars->kerx=(float *)(MKL_malloc(fwidth*sizeof(float),malign));

  pars->kery=(float *)(MKL_malloc(fwidth*sizeof(float),malign));

  for (i=0;i<fwidth;i++) {

    pars->kerx[i]=1.0f/fwidth;

    pars->kery[i]=1.0f/fwidth;

  }


  return 0;

}


int init_fresize_gaussian(

  struct fresize_struct *pars,

  float sigma, // Shape is exp(-(d/sigma)^2/2)

  int hwidth, // Half (truncation) width of kernel. Full is 2*hwidth+1.

  int nsub // Distance between sampled points

)

{

  const int malign=32;

  int fwidth;

  int i,j;

  float sum,x,y;


  pars->method=fresize_1d;

  pars->nsub=nsub;

  pars->hwidth=hwidth;

  fwidth=2*hwidth+1;

  pars->kerx=(float *)(MKL_malloc(fwidth*sizeof(float),malign));

  pars->kery=(float *)(MKL_malloc(fwidth*sizeof(float),malign));

  sum=0.0f;

  for (i=0;i<fwidth;i++) {

    x=(i-hwidth)/sigma;

    y=exp(-x*x/2);

    sum=sum+y;

  }

  for (i=0;i<fwidth;i++) {

    x=(i-hwidth)/sigma;

    y=exp(-x*x/2);

    pars->kerx[i]=y/sum;

    pars->kery[i]=y/sum;

  }


  return 0;

}


int free_fresize(

  struct fresize_struct *pars

)

{

  if (pars->method==fresize_1d) {

    MKL_free (pars->kerx);

    MKL_free (pars->kery);

  }


  return 0;

}


int fsample(

  float *image_in,

  float *image_out,

  int nxin,

  int nyin,

  int nleadin,

  int nxout,

  int nyout,

  int nleadout,

  int nsub,

  int xoff,

  int yoff,

  float fillval

)


{

  int i,j;

  int imin,imax,jmin,jmax;

  float *imp;


// Find first and last indices for which resulting input point is within image

  if (xoff >= 0) imin=0; else imin=((-xoff+nsub-1)/nsub);

  if (xoff <= 0) imax=nxout-1; else imax=minval(((nxin-xoff-1)/nsub),nxout-1);

  if (yoff >= 0) jmin=0; else jmin=((-yoff+nsub-1)/nsub);

  if (yoff <= 0) jmax=nyout-1; else jmax=minval(((nyin-yoff-1)/nsub),nyout-1);


  imp=image_in+yoff*nleadin+xoff;


#pragma omp parallel default(none) private(i,j) shared(image_in,image_out,imp,fillval) shared(imin,imax,jmin,jmax) shared(nxin,nyin,nleadin,nxout,nyout,nleadout,nsub)

  { // Needed to define parallel region


// Fill below valid region

#pragma omp for

    for (j=0; j<jmin; j++) {

      for (i=0; i<nxout; i++) {

        image_out[j*nleadout+i]=fillval;

      } // i=

    } //j=


// Valid heights

#pragma omp for

    for (j=jmin; j<=jmax; j++) {

// Fill to the left

      for (i=0; i<imin; i++) {

        image_out[j*nleadout+i]=fillval;

      } // i=

// Valid region

      for (i=imin; i<=imax; i++) {

//      image_out[j*nleadout+i]=image_in[j*nsub*nleadin+i*nsub];

        image_out[j*nleadout+i]=imp[j*nsub*nleadin+i*nsub];

      } // i=

// Fill to the right

      for (i=imax+1; i<nxout; i++) {

        image_out[j*nleadout+i]=fillval;

      } // i=

    } //j=


// Fill above valid region

#pragma omp for

    for (j=jmax+1; j<nyout; j++) {

      for (i=0; i<nxout; i++) {

        image_out[j*nleadout+i]=fillval;

      } // i=

    } //j=


  } // End of parallel region


  return 0;

}


int fbin(

  float *image_in,

  float *image_out,

  int nxin,

  int nyin,

  int nleadin,

  int nxout,

  int nyout,

  int nleadout,

  int nsub,

  int xoff,

  int yoff,

  float fillval

)


{

  int i,j,i1,j1;

  int imin,imax,jmin,jmax;

  float *imp,*impi;

  float sum;


// Find first and last indices for which resulting input point is within image

  if (xoff >= 0) imin=0; else imin=((-xoff+nsub-1)/nsub);

  if (xoff <= 0) imax=nxout-1; else imax=minval(((nxin-xoff-nsub)/nsub),nxout-1);

  if (yoff >= 0) jmin=0; else jmin=((-yoff+nsub-1)/nsub);

  if (yoff <= 0) jmax=nyout-1; else jmax=minval(((nyin-yoff-nsub)/nsub),nyout-1);


  imp=image_in+yoff*nleadin+xoff;


#pragma omp parallel default(none) private(i,j,i1,j1,impi,sum) shared(image_in,image_out,imp,fillval) shared(imin,imax,jmin,jmax) shared(nxin,nyin,nleadin,nxout,nyout,nleadout,nsub)

  { // Needed to define parallel region


// Fill below valid region

#pragma omp for

    for (j=0; j<jmin; j++) {

      for (i=0; i<nxout; i++) {

        image_out[j*nleadout+i]=fillval;

      } // i=

    } //j=


// Valid heights

#pragma omp for

    for (j=jmin; j<=jmax; j++) {

// Fill to the left

      for (i=0; i<imin; i++) {

        image_out[j*nleadout+i]=fillval;

      } // i=

// Valid region

      for (i=imin; i<=imax; i++) {

        impi=imp+j*nleadin*nsub+i*nsub;

        sum=0.0f;

        for (j1=0; j1<nsub; j1++) {

          for (i1=0; i1<nsub; i1++) {

            sum=sum+impi[i1];

          }

          impi=impi+nleadin;

        }

        image_out[j*nleadout+i]=sum/(nsub*nsub);

      } // i=

// Fill to the right

      for (i=imax+1; i<nxout; i++) {

        image_out[j*nleadout+i]=fillval;

      } // i=

    } //j=


// Fill above valid region

#pragma omp for

    for (j=jmax+1; j<nyout; j++) {

      for (i=0; i<nxout; i++) {

        image_out[j*nleadout+i]=fillval;

      } // i=

    } //j=


  } // End of parallel region


  return 0;

}


int f1d(

  struct fresize_struct *pars,

  float *image_in,

  float *image_out,

  int nxin,

  int nyin,

  int nleadin,

  int nxout,

  int nyout,

  int nleadout,

  int xoff,

  int yoff,

  float fillval

)


{

  const int malign=32;

  char transpose[] = "t";

  char normal[] = "n";

  const int ione = 1;

  const float fone = 1.0f;

  const float fzero = 0.0f;

  int i,j,i1,j1;

  int imin,imax,jmin,jmax;

  float *inp,*inpi,*inpj;

  float sum;

  int nsub,hwidth;

  float *kerx,*kery,*work;

  int xoffl,xoffh,yoffl,yoffh;

  int nwork;

  double t1,t2,t3;

  int n1,n2,nchunk;


  nsub=pars->nsub;

  hwidth=pars->hwidth;

  kerx=pars->kerx;

  kery=pars->kery;


// Find first and last indices for which resulting input point is within image

  xoffl=xoff-hwidth;

  xoffh=xoff+hwidth;

  if (xoffl >= 0) imin=0; else imin=((-xoffl+nsub-1)/nsub);

  if (xoffh <= 0) imax=nxout-1; else imax=minval(((nxin-xoffh-1)/nsub),nxout-1);


  yoffl=yoff-hwidth;

  yoffh=yoff+hwidth;

  if (yoffl >= 0) jmin=0; else jmin=((-yoffl+nsub-1)/nsub);

  if (yoffh <= 0) jmax=nyout-1; else jmax=minval(((nyin-yoffh-1)/nsub),nyout-1);


// Get work array big enough to hold convolution in x

  nwork=nxout*nyin;

  work=(float *)(MKL_malloc(nwork*sizeof(float),malign));


nchunk=64;

n1=(imax-imin+1); // Size of matrix for sgemv

n2=2*hwidth+1; // Size of matrix for sgemv

//t1=dsecnd();

 
#pragma omp parallel default(none) private(i,j,i1,j1,inpi,inpj,sum) shared(image_in,image_out,fillval) shared(imin,imax,jmin,jmax,hwidth,kerx,kery,work,xoffl,yoffl) shared(normal,transpose,n1,n2,nchunk) shared(nxin,nyin,nleadin,nxout,nyout,nleadout,nsub)

  { // Needed to define parallel region


// First convolve in x

// Brute force loop takes longer than it ought to, but have not

// found an obvious solution.

#pragma omp for

    for (j=0;j<nyin;j++) {

      inpj=image_in+j*nleadin+xoffl; // Note offset to match kernels.

      for (i=imin; i<=imax; i++) {

        sum=0.0f;

        inpi=inpj+i*nsub;

        for (i1=0; i1<=2*hwidth; i1++) {

          sum=sum+kerx[i1]*inpi[i1];

        }

        work[j*nleadout+i]=sum;

// This works but is slower

//      work[j*nleadout+i]=sdot(&n2,kerx,&ione,inpj+i*nsub,&ione);

      }

// This clever trick gives error message

//sgemv(transpose,&n2,&n1,&fone,inpj+imin*nsub,&nsub,kerx,&ione,&fzero,work+j*nleadout+imin,&ione);

    }

/*

// This also works but is slower

#pragma omp for

for (j=0;j<nyin;j=j+nchunk) {

  for (i=imin; i<=imax; i++) {

    sgemv(transpose,&n2,&nchunk,&fone,image_in+j*nleadin+xoffl+i*nsub,&nleadin,kerx,&ione,&fzero,work+j*nleadout+i,&nleadout);

}

*/


// Then convolve in y


#pragma omp for

    for (j=jmin; j<=jmax; j++) {

      inpj=work+(yoffl+j*nsub)*nleadout;

/* Old brute force code

      for (i=imin; i<=imax; i++) {

        sum=0.0f;

        for (j1=0; j1<=2*hwidth; j1++) {

          sum=sum+inpj[j1*nleadout+i]*kery[j1];

        }

        image_out[j*nleadout+i]=sum;

      }

*/

      sgemv(normal,&n1,&n2,&fone,inpj+imin,&nleadout,kery,&ione,&fzero,image_out+j*nleadout+imin,&ione);

    }


// Fill below valid region

#pragma omp for

    for (j=0; j<jmin; j++) {

      for (i=0; i<nxout; i++) {

        image_out[j*nleadout+i]=fillval;

      } // i=

    } //j=


// Valid heights

#pragma omp for

    for (j=jmin; j<=jmax; j++) {

// Fill to the left

      for (i=0; i<imin; i++) {

        image_out[j*nleadout+i]=fillval;

      } // i=

// Fill to the right

      for (i=imax+1; i<nxout; i++) {

        image_out[j*nleadout+i]=fillval;

      } // i=

    } //j=


// Fill above valid region

#pragma omp for

    for (j=jmax+1; j<nyout; j++) {

      for (i=0; i<nxout; i++) {

        image_out[j*nleadout+i]=fillval;

      } // i=

    } //j=


  } // End of parallel region


  return 0;

}


int fresize(

  struct fresize_struct *pars,

  float *image_in,

  float *image_out,

  int nxin,

  int nyin,

  int nleadin,

  int nx,

  int ny,

  int nlead,

  int xoff,

  int yoff,

  float fillval

)

{

  int status;


  switch (pars->method) {


  case fresize_sample:

    status=fsample(image_in,image_out,nxin,nyin,nleadin,nx,ny,nlead,pars->nsub,xoff,yoff,fillval);

    break;


  case fresize_bin:

    status=fbin(image_in,image_out,nxin,nyin,nleadin,nx,ny,nlead,pars->nsub,xoff,yoff,fillval);

    break;


  case fresize_1d:

    status=f1d(pars,image_in,image_out,nxin,nyin,nleadin,nx,ny,nlead,xoff,yoff,fillval);

    break;


  default:

    status=1;

    break;

  }


return status;


}


Revision:	1.1
Committed:	Thu Aug 23 08:38:33 2012 UTC (11 years, 1 month ago) by xudong
Content type:	text/plain
Branch:	MAIN
CVS Tags:	Ver_6-4, Ver_9-1, Ver_LATEST, Ver_9-3, Ver_9-41, Ver_9-2, Ver_8-8, Ver_8-2, Ver_8-3, Ver_8-0, Ver_8-1, Ver_8-6, Ver_8-7, Ver_8-4, Ver_8-5, Ver_7-1, Ver_7-0, Ver_9-5, Ver_9-4, Ver_8-10, Ver_8-11, Ver_8-12, Ver_9-0, HEAD
Log Message:	Update functioning code sharp.c per Monica's request
#	Content
1	#include <stdio.h>
2
3	#include <stdlib.h>
4
5	#include <time.h>
6
7	#include <string.h>
8
9	#include <math.h>
10
11	#include <mkl_blas.h>
12
13	#include <mkl_service.h>
14
15	#include <mkl_lapack.h>
16
17	#include <mkl_vml_functions.h>
18
19	#include <omp.h>
20
21	#include "fresize.h"
22
23	#define minval(x,y) (((x) < (y)) ? (x) : (y))
24
25	#define maxval(x,y) (((x) < (y)) ? (y) : (x))
26
27	#define fresize_sample 1
28
29	#define fresize_bin 2
30
31	#define fresize_1d 3
32
33	#define fresize_2d 4
34
35
36
37	double sinc(double x)
38
39	{
40
41	if (fabs(x) < (1.e-10))
42
43	return 1.;
44
45	else
46
47	return sin(M_PIx)/(M_PIx);
48
49	}
50
51
52
53	int init_fresize_sample(
54
55	struct fresize_struct *pars,
56
57	int nsub
58
59	)
60
61	{
62
63
64
65	pars->method=fresize_sample;
66
67	pars->nsub=nsub;
68
69
70
71	return 0;
72
73	}
74
75
76
77	int init_fresize_bin(
78
79	struct fresize_struct *pars,
80
81	int nsub
82
83	)
84
85	{
86
87
88
89	pars->method=fresize_bin;
90
91	pars->nsub=nsub;
92
93
94
95	return 0;
96
97	}
98
99
100
101	int init_fresize_boxcar(
102
103	struct fresize_struct *pars,
104
105	int hwidth, // Half width of boxcar. Full is 2*hwidth+1.
106
107	int nsub // Distance between sampled points
108
109	)
110
111	{
112
113	const int malign=32;
114
115	int fwidth;
116
117	int i,j;
118
119
120
121	pars->method=fresize_1d;
122
123	pars->nsub=nsub;
124
125	pars->hwidth=hwidth;
126
127	fwidth=2*hwidth+1;
128
129	pars->kerx=(float )(MKL_malloc(fwidthsizeof(float),malign));
130
131	pars->kery=(float )(MKL_malloc(fwidthsizeof(float),malign));
132
133	for (i=0;i<fwidth;i++) {
134
135	pars->kerx[i]=1.0f/fwidth;
136
137	pars->kery[i]=1.0f/fwidth;
138
139	}
140
141
142
143	return 0;
144
145	}
146
147
148
149	int init_fresize_gaussian(
150
151	struct fresize_struct *pars,
152
153	float sigma, // Shape is exp(-(d/sigma)^2/2)
154
155	int hwidth, // Half (truncation) width of kernel. Full is 2*hwidth+1.
156
157	int nsub // Distance between sampled points
158
159	)
160
161	{
162
163	const int malign=32;
164
165	int fwidth;
166
167	int i,j;
168
169	float sum,x,y;
170
171
172
173	pars->method=fresize_1d;
174
175	pars->nsub=nsub;
176
177	pars->hwidth=hwidth;
178
179	fwidth=2*hwidth+1;
180
181	pars->kerx=(float )(MKL_malloc(fwidthsizeof(float),malign));
182
183	pars->kery=(float )(MKL_malloc(fwidthsizeof(float),malign));
184
185	sum=0.0f;
186
187	for (i=0;i<fwidth;i++) {
188
189	x=(i-hwidth)/sigma;
190
191	y=exp(-x*x/2);
192
193	sum=sum+y;
194
195	}
196
197	for (i=0;i<fwidth;i++) {
198
199	x=(i-hwidth)/sigma;
200
201	y=exp(-x*x/2);
202
203	pars->kerx[i]=y/sum;
204
205	pars->kery[i]=y/sum;
206
207	}
208
209
210
211	return 0;
212
213	}
214
215
216
217	int free_fresize(
218
219	struct fresize_struct *pars
220
221	)
222
223	{
224
225	if (pars->method==fresize_1d) {
226
227	MKL_free (pars->kerx);
228
229	MKL_free (pars->kery);
230
231	}
232
233
234
235	return 0;
236
237	}
238
239
240
241	int fsample(
242
243	float *image_in,
244
245	float *image_out,
246
247	int nxin,
248
249	int nyin,
250
251	int nleadin,
252
253	int nxout,
254
255	int nyout,
256
257	int nleadout,
258
259	int nsub,
260
261	int xoff,
262
263	int yoff,
264
265	float fillval
266
267	)
268
269
270
271	{
272
273	int i,j;
274
275	int imin,imax,jmin,jmax;
276
277	float *imp;
278
279
280
281	// Find first and last indices for which resulting input point is within image
282
283	if (xoff >= 0) imin=0; else imin=((-xoff+nsub-1)/nsub);
284
285	if (xoff <= 0) imax=nxout-1; else imax=minval(((nxin-xoff-1)/nsub),nxout-1);
286
287	if (yoff >= 0) jmin=0; else jmin=((-yoff+nsub-1)/nsub);
288
289	if (yoff <= 0) jmax=nyout-1; else jmax=minval(((nyin-yoff-1)/nsub),nyout-1);
290
291
292
293	imp=image_in+yoff*nleadin+xoff;
294
295
296
297	#pragma omp parallel default(none) private(i,j) shared(image_in,image_out,imp,fillval) shared(imin,imax,jmin,jmax) shared(nxin,nyin,nleadin,nxout,nyout,nleadout,nsub)
298
299	{ // Needed to define parallel region
300
301
302
303	// Fill below valid region
304
305	#pragma omp for
306
307	for (j=0; j<jmin; j++) {
308
309	for (i=0; i<nxout; i++) {
310
311	image_out[j*nleadout+i]=fillval;
312
313	} // i=
314
315	} //j=
316
317
318
319	// Valid heights
320
321	#pragma omp for
322
323	for (j=jmin; j<=jmax; j++) {
324
325	// Fill to the left
326
327	for (i=0; i<imin; i++) {
328
329	image_out[j*nleadout+i]=fillval;
330
331	} // i=
332
333	// Valid region
334
335	for (i=imin; i<=imax; i++) {
336
337	// image_out[jnleadout+i]=image_in[jnsubnleadin+insub];
338
339	image_out[jnleadout+i]=imp[jnsubnleadin+insub];
340
341	} // i=
342
343	// Fill to the right
344
345	for (i=imax+1; i<nxout; i++) {
346
347	image_out[j*nleadout+i]=fillval;
348
349	} // i=
350
351	} //j=
352
353
354
355	// Fill above valid region
356
357	#pragma omp for
358
359	for (j=jmax+1; j<nyout; j++) {
360
361	for (i=0; i<nxout; i++) {
362
363	image_out[j*nleadout+i]=fillval;
364
365	} // i=
366
367	} //j=
368
369
370
371	} // End of parallel region
372
373
374
375	return 0;
376
377	}
378
379
380
381	int fbin(
382
383	float *image_in,
384
385	float *image_out,
386
387	int nxin,
388
389	int nyin,
390
391	int nleadin,
392
393	int nxout,
394
395	int nyout,
396
397	int nleadout,
398
399	int nsub,
400
401	int xoff,
402
403	int yoff,
404
405	float fillval
406
407	)
408
409
410
411	{
412
413	int i,j,i1,j1;
414
415	int imin,imax,jmin,jmax;
416
417	float imp,impi;
418
419	float sum;
420
421
422
423	// Find first and last indices for which resulting input point is within image
424
425	if (xoff >= 0) imin=0; else imin=((-xoff+nsub-1)/nsub);
426
427	if (xoff <= 0) imax=nxout-1; else imax=minval(((nxin-xoff-nsub)/nsub),nxout-1);
428
429	if (yoff >= 0) jmin=0; else jmin=((-yoff+nsub-1)/nsub);
430
431	if (yoff <= 0) jmax=nyout-1; else jmax=minval(((nyin-yoff-nsub)/nsub),nyout-1);
432
433
434
435	imp=image_in+yoff*nleadin+xoff;
436
437
438
439	#pragma omp parallel default(none) private(i,j,i1,j1,impi,sum) shared(image_in,image_out,imp,fillval) shared(imin,imax,jmin,jmax) shared(nxin,nyin,nleadin,nxout,nyout,nleadout,nsub)
440
441	{ // Needed to define parallel region
442
443
444
445	// Fill below valid region
446
447	#pragma omp for
448
449	for (j=0; j<jmin; j++) {
450
451	for (i=0; i<nxout; i++) {
452
453	image_out[j*nleadout+i]=fillval;
454
455	} // i=
456
457	} //j=
458
459
460
461	// Valid heights
462
463	#pragma omp for
464
465	for (j=jmin; j<=jmax; j++) {
466
467	// Fill to the left
468
469	for (i=0; i<imin; i++) {
470
471	image_out[j*nleadout+i]=fillval;
472
473	} // i=
474
475	// Valid region
476
477	for (i=imin; i<=imax; i++) {
478
479	impi=imp+jnleadinnsub+i*nsub;
480
481	sum=0.0f;
482
483	for (j1=0; j1<nsub; j1++) {
484
485	for (i1=0; i1<nsub; i1++) {
486
487	sum=sum+impi[i1];
488
489	}
490
491	impi=impi+nleadin;
492
493	}
494
495	image_out[jnleadout+i]=sum/(nsubnsub);
496
497	} // i=
498
499	// Fill to the right
500
501	for (i=imax+1; i<nxout; i++) {
502
503	image_out[j*nleadout+i]=fillval;
504
505	} // i=
506
507	} //j=
508
509
510
511	// Fill above valid region
512
513	#pragma omp for
514
515	for (j=jmax+1; j<nyout; j++) {
516
517	for (i=0; i<nxout; i++) {
518
519	image_out[j*nleadout+i]=fillval;
520
521	} // i=
522
523	} //j=
524
525
526
527	} // End of parallel region
528
529
530
531	return 0;
532
533	}
534
535
536
537	int f1d(
538
539	struct fresize_struct *pars,
540
541	float *image_in,
542
543	float *image_out,
544
545	int nxin,
546
547	int nyin,
548
549	int nleadin,
550
551	int nxout,
552
553	int nyout,
554
555	int nleadout,
556
557	int xoff,
558
559	int yoff,
560
561	float fillval
562
563	)
564
565
566
567	{
568
569	const int malign=32;
570
571	char transpose[] = "t";
572
573	char normal[] = "n";
574
575	const int ione = 1;
576
577	const float fone = 1.0f;
578
579	const float fzero = 0.0f;
580
581	int i,j,i1,j1;
582
583	int imin,imax,jmin,jmax;
584
585	float inp,inpi,*inpj;
586
587	float sum;
588
589	int nsub,hwidth;
590
591	float kerx,kery,*work;
592
593	int xoffl,xoffh,yoffl,yoffh;
594
595	int nwork;
596
597	double t1,t2,t3;
598
599	int n1,n2,nchunk;
600
601
602
603	nsub=pars->nsub;
604
605	hwidth=pars->hwidth;
606
607	kerx=pars->kerx;
608
609	kery=pars->kery;
610
611
612
613	// Find first and last indices for which resulting input point is within image
614
615	xoffl=xoff-hwidth;
616
617	xoffh=xoff+hwidth;
618
619	if (xoffl >= 0) imin=0; else imin=((-xoffl+nsub-1)/nsub);
620
621	if (xoffh <= 0) imax=nxout-1; else imax=minval(((nxin-xoffh-1)/nsub),nxout-1);
622
623
624
625	yoffl=yoff-hwidth;
626
627	yoffh=yoff+hwidth;
628
629	if (yoffl >= 0) jmin=0; else jmin=((-yoffl+nsub-1)/nsub);
630
631	if (yoffh <= 0) jmax=nyout-1; else jmax=minval(((nyin-yoffh-1)/nsub),nyout-1);
632
633
634
635	// Get work array big enough to hold convolution in x
636
637	nwork=nxout*nyin;
638
639	work=(float )(MKL_malloc(nworksizeof(float),malign));
640
641
642
643	nchunk=64;
644
645	n1=(imax-imin+1); // Size of matrix for sgemv
646
647	n2=2*hwidth+1; // Size of matrix for sgemv
648
649	//t1=dsecnd();
650
651
652
653	#pragma omp parallel default(none) private(i,j,i1,j1,inpi,inpj,sum) shared(image_in,image_out,fillval) shared(imin,imax,jmin,jmax,hwidth,kerx,kery,work,xoffl,yoffl) shared(normal,transpose,n1,n2,nchunk) shared(nxin,nyin,nleadin,nxout,nyout,nleadout,nsub)
654
655	{ // Needed to define parallel region
656
657
658
659	// First convolve in x
660
661	// Brute force loop takes longer than it ought to, but have not
662
663	// found an obvious solution.
664
665	#pragma omp for
666
667	for (j=0;j<nyin;j++) {
668
669	inpj=image_in+j*nleadin+xoffl; // Note offset to match kernels.
670
671	for (i=imin; i<=imax; i++) {
672
673	sum=0.0f;
674
675	inpi=inpj+i*nsub;
676
677	for (i1=0; i1<=2*hwidth; i1++) {
678
679	sum=sum+kerx[i1]*inpi[i1];
680
681	}
682
683	work[j*nleadout+i]=sum;
684
685	// This works but is slower
686
687	// work[jnleadout+i]=sdot(&n2,kerx,&ione,inpj+insub,&ione);
688
689	}
690
691	// This clever trick gives error message
692
693	//sgemv(transpose,&n2,&n1,&fone,inpj+iminnsub,&nsub,kerx,&ione,&fzero,work+jnleadout+imin,&ione);
694
695	}
696
697	/*
698
699	// This also works but is slower
700
701	#pragma omp for
702
703	for (j=0;j<nyin;j=j+nchunk) {
704
705	for (i=imin; i<=imax; i++) {
706
707	sgemv(transpose,&n2,&nchunk,&fone,image_in+jnleadin+xoffl+insub,&nleadin,kerx,&ione,&fzero,work+j*nleadout+i,&nleadout);
708
709	}
710
711	*/
712
713
714
715	// Then convolve in y
716
717
718
719	#pragma omp for
720
721	for (j=jmin; j<=jmax; j++) {
722
723	inpj=work+(yoffl+jnsub)nleadout;
724
725	/* Old brute force code
726
727	for (i=imin; i<=imax; i++) {
728
729	sum=0.0f;
730
731	for (j1=0; j1<=2*hwidth; j1++) {
732
733	sum=sum+inpj[j1nleadout+i]kery[j1];
734
735	}
736
737	image_out[j*nleadout+i]=sum;
738
739	}
740
741	*/
742
743	sgemv(normal,&n1,&n2,&fone,inpj+imin,&nleadout,kery,&ione,&fzero,image_out+j*nleadout+imin,&ione);
744
745	}
746
747
748
749	// Fill below valid region
750
751	#pragma omp for
752
753	for (j=0; j<jmin; j++) {
754
755	for (i=0; i<nxout; i++) {
756
757	image_out[j*nleadout+i]=fillval;
758
759	} // i=
760
761	} //j=
762
763
764
765	// Valid heights
766
767	#pragma omp for
768
769	for (j=jmin; j<=jmax; j++) {
770
771	// Fill to the left
772
773	for (i=0; i<imin; i++) {
774
775	image_out[j*nleadout+i]=fillval;
776
777	} // i=
778
779	// Fill to the right
780
781	for (i=imax+1; i<nxout; i++) {
782
783	image_out[j*nleadout+i]=fillval;
784
785	} // i=
786
787	} //j=
788
789
790
791	// Fill above valid region
792
793	#pragma omp for
794
795	for (j=jmax+1; j<nyout; j++) {
796
797	for (i=0; i<nxout; i++) {
798
799	image_out[j*nleadout+i]=fillval;
800
801	} // i=
802
803	} //j=
804
805
806
807	} // End of parallel region
808
809
810
811	return 0;
812
813	}
814
815
816
817
818
819	int fresize(
820
821	struct fresize_struct *pars,
822
823	float *image_in,
824
825	float *image_out,
826
827	int nxin,
828
829	int nyin,
830
831	int nleadin,
832
833	int nx,
834
835	int ny,
836
837	int nlead,
838
839	int xoff,
840
841	int yoff,
842
843	float fillval
844
845	)
846
847	{
848
849	int status;
850
851
852
853	switch (pars->method) {
854
855
856
857	case fresize_sample:
858
859	status=fsample(image_in,image_out,nxin,nyin,nleadin,nx,ny,nlead,pars->nsub,xoff,yoff,fillval);
860
861	break;
862
863
864
865	case fresize_bin:
866
867	status=fbin(image_in,image_out,nxin,nyin,nleadin,nx,ny,nlead,pars->nsub,xoff,yoff,fillval);
868
869	break;
870
871
872
873	case fresize_1d:
874
875	status=f1d(pars,image_in,image_out,nxin,nyin,nleadin,nx,ny,nlead,xoff,yoff,fillval);
876
877	break;
878
879
880
881	default:
882
883	status=1;
884
885	break;
886
887	}
888
889
890
891	return status;
892
893
894
895	}
896
897
898