Working cleanset 5 version

src/do_itos.c

@@ -145,7 +145,8 @@ void    newdo_itos(double **im, double dlsds, double zs,
     register    int ii, jj, k, is, js;
     const extern  struct  g_mode  M;
     const extern  struct g_pixel  ps;
-
+    int     offscl=0;  // 0 if OK, 1 if offscale
+    double xpix,ypix;
 
     double ech=1.0;
 
@@ -153,6 +154,10 @@ void    newdo_itos(double **im, double dlsds, double zs,
     int npcont[10];
 
 
+    /*
+     * Read and store all contours temporarily
+     */
+
     for (k = 0; k < imFrame.ncont; k++)
     {
         npcont[k] = readlist(I[k], imFrame.contfile[k]);
@@ -168,6 +173,8 @@ void    newdo_itos(double **im, double dlsds, double zs,
     /*
      * For each pixel in the source plane
      */
+
+   #pragma omp parallel for schedule(dynamic,1)
     for (is = 0; is < ps.nx; is++)
     {
 	    for(js = 0; js < ps.ny; js++)
@@ -177,8 +184,8 @@ void    newdo_itos(double **im, double dlsds, double zs,
 		    struct point spos;
                     struct point pImage[NIMAX]; // list of image positions for a given source position
 
-		    spos.x=ps.xmin+((double)is+1)*ps.pixelx;
-		    spos.y=ps.ymin+((double)js+1)*ps.pixely;
+		    spos.x=ps.xmin+((double)is)*ps.pixelx;
+		    spos.y=ps.ymin+((double)js)*ps.pixely;
 
 		    /*
 		     * Compute all image positions from this source location
@@ -198,18 +205,17 @@ void    newdo_itos(double **im, double dlsds, double zs,
                          {
                               if(inconvexe(pImage[k],npcont[kcont],I[kcont])>0)
 		              {
-				   ii=((pImage[k].x-imFrame.xmin))/imFrame.pixelx+0.5;
-				   jj=((pImage[k].y-imFrame.ymin))/imFrame.pixely+0.5;
+                                   double xw = pImage[k].x / (-3600.) / cos(M.ref_dec * DTR) + M.ref_ra;
+                                   double yw = pImage[k].y / 3600. + M.ref_dec;
+   #pragma omp critical
+				   {
+				   wcs2pix(imFrame.wcsinfo, xw , yw, &xpix, &ypix, &offscl);
+				   }
+				   imult[js][is] += 1;
+				   //nearest pixel
+				   //source[js][is]=(source[js][is]*(imult[js][is]-1)+ech*im[(int)(ypix-0.5)][(int)(xpix-0.5)])/((double)imult[js][is]);
 				   //do bilinear interpolation instead
-                                   imult[js][is] += 1;
-		                   source[js][is]=(source[js][is]*(imult[js][is]-1)+ech*im[(int)jj][(int)ii])/((double)imult[js][is]);
-				   if((is==31)&&(js==35))
-			           {
-					   printf("%d %d %f\n",ii,jj,im[(int)jj][(int)ii]);
-					   printf("x=%f y=%f\n",pImage[k].x,pImage[k].y);
-					   printf("xmin=%f ymin=%f\n",imFrame.xmin,imFrame.ymin);
-					   printf("pixelx=%f pixely=%f\n",imFrame.pixelx,imFrame.pixely);
-			           }
+				   source[js][is]=(source[js][is]*(imult[js][is]-1)+ech*bilinear(im,xpix,ypix,(int)imFrame.wcsinfo->nxpix,(int)imFrame.wcsinfo->nypix))/((double)imult[js][is]);
 			      }
 			 }
 		    }