globalhs/apps/imageinterp.c


/*
 *  imageinterp (modified from obs2helio)                       
 *
 *  Purpose:
 *     Interpolate CCD data to estimate value that would
 *     be obtained at specified equal incrememts of heliographic
 *     longitude and sine of latitude.
 *
 *  Description:
 *
 *     Employs bi-linear or cubic convolution interpolation to map a
 *     portion of the solar disk to a rectangular array "cols" by "rows" in
 *     width and height respectively.
 *
 *  Reference:
 *
 *     Green, Robin M. "Spherical Astronomy," Cambridge University
 *     Press, Cambridge, 1985.
 *
 *  Responsible:  Kay Leibrand                   KLeibrand@solar.Stanford.EDU
 *
 *  Restrictions:
 *     The number of output map rows must be even.
 *     The number of output map columns must be less than MAXCOLS.
 *     Assumes orientation is a 4 character string with one of the following
 *     8 legal values: SESW SENE SWSE SWNW NWNE NWSW NENW NESE 
 *
 *  Bugs:
 *     Possible division by 0 
 *
 *  Planned updates:
 *     Optimize.
 *
 *  Revision history is at end of file.
 */

#include <math.h>
#include "projection.h"

char *cvsinfo_imageinterp = "cvsinfo: $Header: imageinterp.c,v $";

const int kMaxCols = 8192;
const double kOmegaCarr = (360 /  27.27527); /* degrees/day - synodic Carrington rotation rate */
const double kRadsPerDegree = M_PI/180.;

static void Ccker(double *u, double s);
static double Ccintd(double *f, int nx, int ny, double x, double y);
static double Linintd(double *f, int nx, int ny, double x, double y);
static void Distort(double x, 
                    double y, 
                    double rsun, 
                    int sign, 
                    double *xd, 
                    double *yd, 
                    const LIBPROJECTION_Dist_t *distpars);

/* Calculate the interpolation kernel. */
void Ccker(double *u, double s)
{
   double s2, s3;
     
   s2= s * s;
   s3= s2 * s;
   u[0] = s2 - 0.5 * (s3 + s);
   u[1] = 1.5*s3 - 2.5*s2 + 1.0;
   u[2] = -1.5*s3 + 2.0*s2 + 0.5*s;
   u[3] = 0.5 * (s3 - s2);
}

/* Cubic convolution interpolation */
/*
 *  Cubic convolution interpolation, based on GONG software, received Apr-88
 *  Interpolation by cubic convolution in 2 dimensions with 32-bit double data
 *
 *  Reference:
 *    R.G. Keys in IEEE Trans. Acoustics, Speech, and Signal Processing, 
 *    Vol. 29, 1981, pp. 1153-1160.
 *
 *  Usage:
 *    double ccintd (double *f, int nx, int ny, double x, double y)
 *
 *  Bugs:
 *    Currently interpolation within one pixel of edge is not
 *      implemented.  If x or y is in this range or off the picture, the
 *      function value returned is zero.
 *    Only works for double data.
 */
double Ccintd(double *f, int nx, int ny, double x, double y)
{
   double  ux[4], uy[4];
   /* Sum changed to double to speed up things */
   double sum;

   int ix, iy, ix1, iy1, i, j;
     
   if (x < 1. || x >= (double)(nx-2) ||
       y < 1. || y >= (double)(ny-2))
     return 0.0;
     
   ix = (int)x;
   iy = (int)y;
   Ccker(ux,  x - (double)ix);
   Ccker(uy,  y - (double)iy);
     
   ix1 = ix - 1;
   iy1 = iy - 1;
   sum = 0.;
   for (i = 0; i < 4; i++)
     for (j = 0; j < 4; j++)
       sum = sum + f[(iy1+i) * nx + ix1 + j] * uy[i] * ux[j];
   return sum;
}

/*
 *  Bilinear interpolation, based on pipeLNU remap.
 *
 *  Reference: 
 *     Abramowitz, Milton, and Stegun, Irene, "Handbook of Mathematical
 *     Functions," Dover Publications, New York, 1964.
 *
 *  Usage:
 *    double linintd (double *f, int nx, int ny, double x, double y)
 *
 *  Bugs:
 *    Only works for double data.
 */

double Linintd(double *f, int nx, int ny, double x, double y) 
{
   double p0, p1, q0, q1;          /* weights                              */
   int ilow, jlow;   /* selected CCD pixels                  */
   double *fptr;                    /* input array temp                     */
   double u;                        
     
   ilow = (int)floor(x);
   jlow = (int)floor(y);
   if (ilow < 0) ilow = 0;
   if (ilow+1 >= nx) ilow -= 1;
   if (jlow < 0) jlow = 0;
   if (jlow+1 >= ny) jlow -= 1;
   p1 = x - ilow;
   p0 = 1.0 - p1;
   q1 = y - jlow;
   q0 = 1.0 - q1;
     
   /* Bilinear interpolation from four data points. */
   u = 0.0;
   fptr = f + jlow*nx + ilow;
   u += p0 * q0 * *fptr++;
   u += p1 * q0 * *fptr;
   fptr += nx - 1;
   u += p0 * q1 * *fptr++;
   u += p1 * q1 * *fptr;
   return u;
}


int SetDistort(int dist, 
               double cubic, 
               double alpha, 
               double beta, 
               double feff,
               LIBPROJECTION_Dist_t *dOut)
{
   int error = 0;

   if (dOut != NULL)
   {
      int status = 0;

      dOut->disttype = dist;

      /* disttype == 0 is 'no distortion' and is not an error */
      if (dOut->disttype != 0)
      {
         if (dOut->disttype == 1) 
         { 
            /* FD */
            dOut->xc = 511.5;
            dOut->yc = 511.5;
            dOut->scale = 1.0;
         }
         else if (dOut->disttype == 2)
         { 
            /* vw */
            dOut->xc = 95.1;
            dOut->yc = 95.1;
            dOut->scale = 5.0;
         }
         else 
         {
            error = 1;
         }

         if (!error)
         {
            dOut->feff = feff;
            dOut->alpha = alpha * kRadsPerDegree;
            beta *= kRadsPerDegree;
            dOut->cosbeta = cos(beta);
            dOut->sinbeta = sin(beta);
            dOut->cdist = cubic;

            if (status != CMDPARAMS_SUCCESS)
            {
               error = 1;
            }

         }

      } /* disttype != 0*/
   }
   else
   {
      error = 1;
   }

   return error;
}


void Distort(double x, 
             double y, 
             double rsun, 
             int sign, 
             double *xd, 
             double *yd, 
             const LIBPROJECTION_Dist_t *distpars)
{
   /*
     For sign=+1 transforms ideal coordinates (x,y) to distorted
     coordinates (xd,yd). rsun needed to correct for scale error.
     For sign=-1 ideal and distorted coordinates are swapped.
     Units are pixels from lower left corner.
     Adapted from /home/schou/calib/dist/distort.pro
     Note that rsun must be passed in FD pixels!
   */
     
   double xx,yy,xxl,yyl,x0,y0,x00,y00,x0l,y0l,epsx,epsy,rdist;
     
   if (distpars->disttype == 0) 
   {
      *xd=x;
      *yd=y;
      return;
   }
   /* Convert fo FD pixel scale*/

     rsun=rsun*distpars->scale;

   /* Don't do anyway since rsun already corrected. Don't get me started... */
   /* above comment no longer applies, rsun is now untouched in o2helio (v2helio) */  
     
   /* Shift zero point to nominal image center and convert to FD pixel scale */
     
   x00=distpars->scale*(x-distpars->xc);
   y00=distpars->scale*(y-distpars->yc);
     
   /* Radial distortion */
     
   rdist=distpars->cdist*(x00*x00+y00*y00-rsun*rsun);
     
   x0=x00+rdist*x00;
   y0=y00+rdist*y00;
     
   /* Distortion due to CCD tilt */
     
   /* First rotate coordinate system */
   x0l=x0*distpars->cosbeta+y0*distpars->sinbeta;
   y0l=-x0*distpars->sinbeta+y0*distpars->cosbeta;
     
   /* Apply distortion */
   xxl=x0l*(1.-distpars->alpha*distpars->alpha/4.+y0l/distpars->feff*distpars->alpha);
   yyl=y0l*(1.+distpars->alpha*distpars->alpha/4.+y0l/distpars->feff*distpars->alpha) - 
     rsun*rsun/distpars->feff*distpars->alpha;
     
   /* Rotate coordinates back */
   xx=xxl*distpars->cosbeta-yyl*distpars->sinbeta;
   yy=xxl*distpars->sinbeta+yyl*distpars->cosbeta;
     
   /* Total distortion */
     
   epsx=xx-x00;
   epsy=yy-y00;
     
   /* Return distorted values. Cheating with inverse transform. */
     
   *xd=x+sign*epsx/distpars->scale;
   *yd=y+sign*epsy/distpars->scale;
}


int imageinterp(
        float *V,  /* input velocity array    */
        /* assumed declared V[ypixels][xpixels]  */
        float *U,  /* output projected array   */
        
        int xpixels, /* x width of input array */
        int ypixels, /* y width of input array */
        double   x0, /* x pixel address of disk center          */
        double   y0, /* y pixel address of disk center          */
        double   P,  /* angle between CCD y-axis and solar vertical */
                     /* positive to the east (CCW) */
        double   rsun,  /* pixels */
        double   Rmax, /* maximum disk radius to use (e.g. 0.95)        */
        int      NaN_beyond_rmax,
        int interpolation, /* option */
        int cols,  /* width of output array   */
        int rows,  /* height of output array   */
        double   vsign,
        const LIBPROJECTION_Dist_t *distpars)
{
     
   float u;   /* output temp    */
   double x, y;   /* CCD location of desired point */
   int row, col;   /* index into output array  */
   int rowoffset;
   double xratio, yratio, Rmax2;

   long i;
   double *vdp;
   float *vp;
   double xtmp, ytmp;
     
   double *VD;
   VD=(double *)(malloc(xpixels*ypixels*sizeof(double)));
   vdp=VD;
   vp=V;
   for (i=0;i<xpixels*ypixels;i++) *vdp++=(double)*vp++;

   Rmax2 = Rmax*Rmax*rsun*rsun;
     
   if (cols > kMaxCols) 
   {
      return kLIBPROJECTION_DimensionMismatch;
   }

     
   if (interpolation != kLIBPROJECTION_InterCubic && interpolation != kLIBPROJECTION_InterBilinear) 
   {     
      return kLIBPROJECTION_UnsupportedInterp;
   }

   xratio=(double)xpixels/cols;
   yratio=(double)ypixels/rows;
     
   for (row = 0; row < rows; row++) {
 
      rowoffset = row*cols;
  
      for (col = 0; col < cols; col++) {

         xtmp = (col + 0.5)*xratio - 0.5 - x0;
         ytmp = (row + 0.5)*yratio - 0.5 - y0;
         if ((xtmp*xtmp + ytmp*ytmp) >= Rmax2)
         {
           *(U + rowoffset + col) = (NaN_beyond_rmax) ? DRMS_MISSING_FLOAT : (float)0.0;
           continue;
         }
         x= x0 + xtmp*cos(P) - ytmp*sin(P);
         y= y0 + xtmp*sin(P) + ytmp*cos(P);

         Distort(x, y, rsun, +1, &x, &y, distpars);

         if (interpolation == kLIBPROJECTION_InterCubic) 
         {
            u = (float)(vsign * Ccintd (VD,xpixels,ypixels,x,y));
         }
         else if (interpolation == kLIBPROJECTION_InterBilinear) 
         {
            u = (float)(vsign * Linintd (VD,xpixels,ypixels,x,y));
         }

         *(U + rowoffset + col) = u;
      }
   }
   free(VD);
   return kLIBPROJECTION_Success;
}
Revision:	1.4
Committed:	Fri May 3 20:47:15 2013 UTC (10 years, 4 months ago) by tplarson
Content type:	text/plain
Branch:	MAIN
CVS Tags:	globalhs_version_5, Ver_8-7, Ver_8-5, globalhs_version_23, globalhs_version_22, globalhs_version_21, globalhs_version_20, Ver_LATEST, globalhs_version_24, Ver_8-3, globalhs_version_8, globalhs_version_9, globalhs_version_0, globalhs_version_1, globalhs_version_2, globalhs_version_3, globalhs_version_4, Ver_9-41, globalhs_version_6, globalhs_version_7, Ver_9-5, Ver_8-8, globalhs_version_19, Ver_8-2, Ver_8-10, Ver_8-1, Ver_8-6, Ver_9-1, Ver_8-4, Ver_9-2, globalhs_version_12, globalhs_version_13, globalhs_version_10, globalhs_version_11, globalhs_version_16, globalhs_version_17, globalhs_version_14, globalhs_version_15, globalhs_version_18, Ver_9-4, Ver_9-3, Ver_8-11, Ver_8-12, Ver_9-0, HEAD
Changes since 1.3:	+12 -42 lines
Log Message:	start using projection.h
#	Content
1
2	/*
3	* imageinterp (modified from obs2helio)
4	*
5	* Purpose:
6	* Interpolate CCD data to estimate value that would
7	* be obtained at specified equal incrememts of heliographic
8	* longitude and sine of latitude.
9	*
10	* Description:
11	*
12	* Employs bi-linear or cubic convolution interpolation to map a
13	* portion of the solar disk to a rectangular array "cols" by "rows" in
14	* width and height respectively.
15	*
16	* Reference:
17	*
18	* Green, Robin M. "Spherical Astronomy," Cambridge University
19	* Press, Cambridge, 1985.
20	*
21	* Responsible: Kay Leibrand KLeibrand@solar.Stanford.EDU
22	*
23	* Restrictions:
24	* The number of output map rows must be even.
25	* The number of output map columns must be less than MAXCOLS.
26	* Assumes orientation is a 4 character string with one of the following
27	* 8 legal values: SESW SENE SWSE SWNW NWNE NWSW NENW NESE
28	*
29	* Bugs:
30	* Possible division by 0
31	*
32	* Planned updates:
33	* Optimize.
34	*
35	* Revision history is at end of file.
36	*/
37
38	#include <math.h>
39	#include "projection.h"
40
41	char *cvsinfo_imageinterp = "cvsinfo: $Header: imageinterp.c,v $";
42
43	const int kMaxCols = 8192;
44	const double kOmegaCarr = (360 / 27.27527); /* degrees/day - synodic Carrington rotation rate */
45	const double kRadsPerDegree = M_PI/180.;
46
47	static void Ccker(double *u, double s);
48	static double Ccintd(double *f, int nx, int ny, double x, double y);
49	static double Linintd(double *f, int nx, int ny, double x, double y);
50	static void Distort(double x,
51	double y,
52	double rsun,
53	int sign,
54	double *xd,
55	double *yd,
56	const LIBPROJECTION_Dist_t *distpars);
57
58	/* Calculate the interpolation kernel. */
59	void Ccker(double *u, double s)
60	{
61	double s2, s3;
62
63	s2= s * s;
64	s3= s2 * s;
65	u[0] = s2 - 0.5 * (s3 + s);
66	u[1] = 1.5s3 - 2.5s2 + 1.0;
67	u[2] = -1.5s3 + 2.0s2 + 0.5*s;
68	u[3] = 0.5 * (s3 - s2);
69	}
70
71	/* Cubic convolution interpolation */
72	/*
73	* Cubic convolution interpolation, based on GONG software, received Apr-88
74	* Interpolation by cubic convolution in 2 dimensions with 32-bit double data
75	*
76	* Reference:
77	* R.G. Keys in IEEE Trans. Acoustics, Speech, and Signal Processing,
78	* Vol. 29, 1981, pp. 1153-1160.
79	*
80	* Usage:
81	* double ccintd (double *f, int nx, int ny, double x, double y)
82	*
83	* Bugs:
84	* Currently interpolation within one pixel of edge is not
85	* implemented. If x or y is in this range or off the picture, the
86	* function value returned is zero.
87	* Only works for double data.
88	*/
89	double Ccintd(double *f, int nx, int ny, double x, double y)
90	{
91	double ux[4], uy[4];
92	/* Sum changed to double to speed up things */
93	double sum;
94
95	int ix, iy, ix1, iy1, i, j;
96
97	if (x < 1. \|\| x >= (double)(nx-2) \|\|
98	y < 1. \|\| y >= (double)(ny-2))
99	return 0.0;
100
101	ix = (int)x;
102	iy = (int)y;
103	Ccker(ux, x - (double)ix);
104	Ccker(uy, y - (double)iy);
105
106	ix1 = ix - 1;
107	iy1 = iy - 1;
108	sum = 0.;
109	for (i = 0; i < 4; i++)
110	for (j = 0; j < 4; j++)
111	sum = sum + f[(iy1+i) * nx + ix1 + j] * uy[i] * ux[j];
112	return sum;
113	}
114
115	/*
116	* Bilinear interpolation, based on pipeLNU remap.
117	*
118	* Reference:
119	* Abramowitz, Milton, and Stegun, Irene, "Handbook of Mathematical
120	* Functions," Dover Publications, New York, 1964.
121	*
122	* Usage:
123	* double linintd (double *f, int nx, int ny, double x, double y)
124	*
125	* Bugs:
126	* Only works for double data.
127	*/
128
129	double Linintd(double *f, int nx, int ny, double x, double y)
130	{
131	double p0, p1, q0, q1; /* weights */
132	int ilow, jlow; /* selected CCD pixels */
133	double fptr; / input array temp */
134	double u;
135
136	ilow = (int)floor(x);
137	jlow = (int)floor(y);
138	if (ilow < 0) ilow = 0;
139	if (ilow+1 >= nx) ilow -= 1;
140	if (jlow < 0) jlow = 0;
141	if (jlow+1 >= ny) jlow -= 1;
142	p1 = x - ilow;
143	p0 = 1.0 - p1;
144	q1 = y - jlow;
145	q0 = 1.0 - q1;
146
147	/* Bilinear interpolation from four data points. */
148	u = 0.0;
149	fptr = f + jlow*nx + ilow;
150	u += p0 * q0 * *fptr++;
151	u += p1 * q0 * *fptr;
152	fptr += nx - 1;
153	u += p0 * q1 * *fptr++;
154	u += p1 * q1 * *fptr;
155	return u;
156	}
157
158
159	int SetDistort(int dist,
160	double cubic,
161	double alpha,
162	double beta,
163	double feff,
164	LIBPROJECTION_Dist_t *dOut)
165	{
166	int error = 0;
167
168	if (dOut != NULL)
169	{
170	int status = 0;
171
172	dOut->disttype = dist;
173
174	/* disttype == 0 is 'no distortion' and is not an error */
175	if (dOut->disttype != 0)
176	{
177	if (dOut->disttype == 1)
178	{
179	/* FD */
180	dOut->xc = 511.5;
181	dOut->yc = 511.5;
182	dOut->scale = 1.0;
183	}
184	else if (dOut->disttype == 2)
185	{
186	/* vw */
187	dOut->xc = 95.1;
188	dOut->yc = 95.1;
189	dOut->scale = 5.0;
190	}
191	else
192	{
193	error = 1;
194	}
195
196	if (!error)
197	{
198	dOut->feff = feff;
199	dOut->alpha = alpha * kRadsPerDegree;
200	beta *= kRadsPerDegree;
201	dOut->cosbeta = cos(beta);
202	dOut->sinbeta = sin(beta);
203	dOut->cdist = cubic;
204
205	if (status != CMDPARAMS_SUCCESS)
206	{
207	error = 1;
208	}
209
210	}
211
212	} /* disttype != 0*/
213	}
214	else
215	{
216	error = 1;
217	}
218
219	return error;
220	}
221
222
223	void Distort(double x,
224	double y,
225	double rsun,
226	int sign,
227	double *xd,
228	double *yd,
229	const LIBPROJECTION_Dist_t *distpars)
230	{
231	/*
232	For sign=+1 transforms ideal coordinates (x,y) to distorted
233	coordinates (xd,yd). rsun needed to correct for scale error.
234	For sign=-1 ideal and distorted coordinates are swapped.
235	Units are pixels from lower left corner.
236	Adapted from /home/schou/calib/dist/distort.pro
237	Note that rsun must be passed in FD pixels!
238	*/
239
240	double xx,yy,xxl,yyl,x0,y0,x00,y00,x0l,y0l,epsx,epsy,rdist;
241
242	if (distpars->disttype == 0)
243	{
244	*xd=x;
245	*yd=y;
246	return;
247	}
248	/* Convert fo FD pixel scale*/
249
250	rsun=rsun*distpars->scale;
251
252	/* Don't do anyway since rsun already corrected. Don't get me started... */
253	/* above comment no longer applies, rsun is now untouched in o2helio (v2helio) */
254
255	/* Shift zero point to nominal image center and convert to FD pixel scale */
256
257	x00=distpars->scale*(x-distpars->xc);
258	y00=distpars->scale*(y-distpars->yc);
259
260	/* Radial distortion */
261
262	rdist=distpars->cdist(x00x00+y00y00-rsunrsun);
263
264	x0=x00+rdist*x00;
265	y0=y00+rdist*y00;
266
267	/* Distortion due to CCD tilt */
268
269	/* First rotate coordinate system */
270	x0l=x0distpars->cosbeta+y0distpars->sinbeta;
271	y0l=-x0distpars->sinbeta+y0distpars->cosbeta;
272
273	/* Apply distortion */
274	xxl=x0l(1.-distpars->alphadistpars->alpha/4.+y0l/distpars->feff*distpars->alpha);
275	yyl=y0l(1.+distpars->alphadistpars->alpha/4.+y0l/distpars->feff*distpars->alpha) -
276	rsunrsun/distpars->feffdistpars->alpha;
277
278	/* Rotate coordinates back */
279	xx=xxldistpars->cosbeta-yyldistpars->sinbeta;
280	yy=xxldistpars->sinbeta+yyldistpars->cosbeta;
281
282	/* Total distortion */
283
284	epsx=xx-x00;
285	epsy=yy-y00;
286
287	/* Return distorted values. Cheating with inverse transform. */
288
289	xd=x+signepsx/distpars->scale;
290	yd=y+signepsy/distpars->scale;
291	}
292
293
294	int imageinterp(
295	float V, / input velocity array */
296	/* assumed declared V[ypixels][xpixels] */
297	float U, / output projected array */
298
299	int xpixels, /* x width of input array */
300	int ypixels, /* y width of input array */
301	double x0, /* x pixel address of disk center */
302	double y0, /* y pixel address of disk center */
303	double P, /* angle between CCD y-axis and solar vertical */
304	/* positive to the east (CCW) */
305	double rsun, /* pixels */
306	double Rmax, /* maximum disk radius to use (e.g. 0.95) */
307	int NaN_beyond_rmax,
308	int interpolation, /* option */
309	int cols, /* width of output array */
310	int rows, /* height of output array */
311	double vsign,
312	const LIBPROJECTION_Dist_t *distpars)
313	{
314
315	float u; /* output temp */
316	double x, y; /* CCD location of desired point */
317	int row, col; /* index into output array */
318	int rowoffset;
319	double xratio, yratio, Rmax2;
320
321	long i;
322	double *vdp;
323	float *vp;
324	double xtmp, ytmp;
325
326	double *VD;
327	VD=(double )(malloc(xpixelsypixels*sizeof(double)));
328	vdp=VD;
329	vp=V;
330	for (i=0;i<xpixelsypixels;i++) vdp++=(double)*vp++;
331
332	Rmax2 = RmaxRmaxrsun*rsun;
333
334	if (cols > kMaxCols)
335	{
336	return kLIBPROJECTION_DimensionMismatch;
337	}
338
339
340	if (interpolation != kLIBPROJECTION_InterCubic && interpolation != kLIBPROJECTION_InterBilinear)
341	{
342	return kLIBPROJECTION_UnsupportedInterp;
343	}
344
345	xratio=(double)xpixels/cols;
346	yratio=(double)ypixels/rows;
347
348	for (row = 0; row < rows; row++) {
349
350	rowoffset = row*cols;
351
352	for (col = 0; col < cols; col++) {
353
354	xtmp = (col + 0.5)*xratio - 0.5 - x0;
355	ytmp = (row + 0.5)*yratio - 0.5 - y0;
356	if ((xtmpxtmp + ytmpytmp) >= Rmax2)
357	{
358	*(U + rowoffset + col) = (NaN_beyond_rmax) ? DRMS_MISSING_FLOAT : (float)0.0;
359	continue;
360	}
361	x= x0 + xtmpcos(P) - ytmpsin(P);
362	y= y0 + xtmpsin(P) + ytmpcos(P);
363
364	Distort(x, y, rsun, +1, &x, &y, distpars);
365
366	if (interpolation == kLIBPROJECTION_InterCubic)
367	{
368	u = (float)(vsign * Ccintd (VD,xpixels,ypixels,x,y));
369	}
370	else if (interpolation == kLIBPROJECTION_InterBilinear)
371	{
372	u = (float)(vsign * Linintd (VD,xpixels,ypixels,x,y));
373	}
374
375	*(U + rowoffset + col) = u;
376	}
377	}
378	free(VD);
379	return kLIBPROJECTION_Success;
380	}