src/Tools/mrcImage/mrcImageUnbentROI/src/mrcImageUnbentROI.c

   1 #include <stdio.h>
   2 #include <stdlib.h>
   3 #include <string.h>
   4 #include <math.h>
   5 #define GLOBAL_DECLARATION
   6 #include "../inc/config.h"
   7 #define DEBUG
   8 #include "genUtil.h"
   9 #include "mrcImage.h"
  10
  11 typedef struct lmrcImageUnbentROIInfo {
  12         int    nUnbentPoint;
  13         float* unbentPoint;
  14         int Nx;
  15         int Ny;
  16         float deltax;
  17         float deltay;
  18         float correlationWindow;
  19 } lmrcImageUnbentROIInfo;
  20
  21 extern void lmrcImageUnbentROI(mrcImage* out, mrcImage* in, lmrcImageUnbentROIInfo* linfo, int mode);
  22
  23 int
  24 main(int argc, char* argv[])
  25 {
  26         mrcImageUnbentROIInfo  info;
  27         lmrcImageUnbentROIInfo linfo;
  28         mrcImage in;
  29         mrcImage out;
  30         int i;
  31
  32         init0(&info);
  33     argCheck(&info, argc, argv);
  34     init1(&info);
  35         mrcFileRead(&in, info.In, "in main", 0);
  36
  37         linfo.nUnbentPoint = info.flagroiLine/2;
  38         linfo.unbentPoint  = info.roiLine;
  39
  40         if(info.flagWidth) {
  41                 linfo.Nx     = info.Width;
  42                 linfo.deltax = info.x/info.Width;
  43         } else {
  44                 linfo.Nx     = info.x;
  45                 linfo.deltax = 1.0;
  46         }
  47         if(info.flagHeight) {
  48                 linfo.Ny = info.Height;
  49                 if(!info.flagy) {
  50                         info.y = 0;
  51                 }
  52                 linfo.deltay = info.y/info.Height;
  53         } else {
  54                 if(!info.flagy) {
  55                         info.y = 0;
  56                 }
  57                 linfo.deltay = 1.0;
  58                 linfo.Ny     = info.y;
  59         }
  60         if(info.flagShrink) {
  61                 for(i=0; i<linfo.nUnbentPoint*2; i++) {
  62                         linfo.unbentPoint[i] *= info.Shrink;
  63                 }
  64         }
  65
  66         if(info.flagCor) {
  67                 linfo.correlationWindow = info.Cor;
  68         } else {
  69                 linfo.correlationWindow = 0;
  70         }
  71
  72         lmrcImageUnbentROI(&out, &in, &linfo, info.mode);
  73
  74         mrcFileWrite(&out, info.Out, "in main", 0);
  75         exit(EXIT_SUCCESS);
  76 }
  77
  78 void
  79 lmrcImageUnbentROI(mrcImage* out, mrcImage* in, lmrcImageUnbentROIInfo* linfo, int mode)
  80 {
  81         long i, n;
  82         floatVector x, y, p, a, b;
  83         floatVector newx, newy, newp, newa, newb;
  84         floatVector oldx, oldy, oldp, olda, oldb;
  85         double dx, dy;
  86         double px[3];
  87         double py[3];
  88         double vx, vy, l, sx, sy;
  89         double data, data2, s;
  90         long num;
  91
  92         /* Preparation of Spline Table */
  93         n = linfo->nUnbentPoint;
  94         floatVectorInit(&x, n);
  95         floatVectorInit(&y, n);
  96         for(i=0; i<n; i++) {
  97                 x.data[i] = linfo->unbentPoint[2*i];
  98                 y.data[i] = linfo->unbentPoint[2*i + 1];
  99         }
 100         floatVectorInit(&p, n); /* Parameter along a bent line */
 101         floatVectorInit(&a, n);
 102         floatVectorInit(&b, n);
 103         lVectorSplineTable2DMake(&p, &x, &y, &a, &b);
 104
 105         /* Output File */
 106         out->Header = in->Header;
 107         out->HeaderMode = mrcFloatImage;
 108         out->HeaderN.x = linfo->Nx;
 109         if(0<linfo->Ny) {
 110                 out->HeaderN.y = linfo->Ny;
 111         } else {
 112                 out->HeaderN.y = p.data[n-1];
 113         }
 114
 115         out->HeaderLength.x *= linfo->deltax;
 116         out->HeaderLength.y *= linfo->deltay;
 117
 118         mrcInit(out, NULL);
 119
 120         /* With Correlation */
 121         if(0<linfo->correlationWindow) {
 122                 int m;
 123                 int flagIter;
 124                 int flagOK;
 125                 mrcImage ref;
 126                 mrcImage tmp;
 127                 mrcImage tmp2;
 128                 mrcImage avg;
 129                 mrcImageParaTypeRealCoord shift;
 130                 mrcImage cor;
 131                 mrcImageInformation info;
 132
 133                 /* Reference */
 134                 ref.Header    = out->Header;
 135                 ref.HeaderN.x = out->HeaderN.x;
 136                 ref.HeaderN.y = ((int)(linfo->correlationWindow/out->HeaderLength.y/2))*2;
 137                 mrcInit(&ref, NULL);
 138                 /* Temporary for fitting */
 139                 tmp.Header = ref.Header;
 140                 mrcInit(&tmp, NULL);
 141                 tmp2.Header = ref.Header;
 142                 mrcInit(&tmp2, NULL);
 143                 avg.Header = ref.Header;
 144                 mrcInit(&avg, NULL);
 145
 146                 /* Number of new points */
 147                 m = (int)(out->HeaderN.y*out->HeaderLength.y/linfo->correlationWindow) + 1;
 148                 floatVectorInit(&newx, m);
 149                 floatVectorInit(&newy, m);
 150                 floatVectorInit(&newp, m);
 151                 floatVectorInit(&newa, m);
 152                 floatVectorInit(&newb, m);
 153                 flagIter = 0;
 154                 while(0<=flagIter && flagIter<1) {
 155                         if(0==flagIter) { /* Attention: first time, the size of oldx is n, not m */
 156                                 oldx = x;
 157                                 oldy = y;
 158                                 oldp = p;
 159                                 olda = a;
 160                                 oldb = b;
 161                         } else {
 162                                 if(flagIter == 1) {
 163                                         floatVectorInit(&oldx, m);
 164                                         floatVectorInit(&oldy, m);
 165                                         floatVectorInit(&oldp, m);
 166                                         floatVectorInit(&olda, m);
 167                                         floatVectorInit(&oldb, m);
 168                                 }
 169                                 for(i=0; i<m; i++) { oldx.data[i] = newx.data[i]; }
 170                                 for(i=0; i<m; i++) { oldy.data[i] = newy.data[i]; }
 171                                 for(i=0; i<m; i++) { oldp.data[i] = newp.data[i]; }
 172                                 for(i=0; i<m; i++) { olda.data[i] = newa.data[i]; }
 173                                 for(i=0; i<m; i++) { oldb.data[i] = newb.data[i]; }
 174                         }
 175                         DEBUGPRINT1("Correlation  : %d\n", flagIter);
 176                         DEBUGPRINT1("Reference Set: %d\n", flagIter);
 177                         for(dy=0; dy<ref.HeaderN.y; dy++) {
 178                                 lVectorSpline2D((dy-1.0)*linfo->deltay, &(px[0]), &(py[0]), &oldp, &oldx, &oldy, &olda, &oldb);
 179                                 lVectorSpline2D( dy     *linfo->deltay, &(px[1]), &(py[1]), &oldp, &oldx, &oldy, &olda, &oldb);
 180                                 lVectorSpline2D((dy+1.0)*linfo->deltay, &(px[2]), &(py[2]), &oldp, &oldx, &oldy, &olda, &oldb);
 181                                 vy = -(px[2] - px[0])/2.0;
 182                                 vx =  (py[2] - py[0])/2.0;
 183                                 l  = sqrt(SQR(vx) + SQR(vy));
 184                                 vx /= l;
 185                                 vy /= l;
 186                                 for(dx=0; dx<ref.HeaderN.x; dx++) {
 187                                         s  = (dx - ref.HeaderN.x/2.0)*linfo->deltax;
 188                                         sx = px[1] + vx*s;
 189                                         sy = py[1] + vy*s;
 190                                         mrcPixelDataGet(in,   sx, sy, 0.0, &data, mrcPixelMag, mrcPixelHowNearest);
 191                                         if(flagIter==0) {
 192                                                 mrcPixelDataSet(&ref, dx, dy, 0.0,  data, mrcPixelMag);
 193                                         } else {
 194                                                 mrcPixelDataGet(&avg, dx, dy, 0.0, &data2, mrcPixelMag, mrcPixelHowNearest);
 195                                                 mrcPixelDataSet(&ref, dx, dy, 0.0,  data2, mrcPixelMag);
 196                                         }
 197                                         mrcPixelDataSet(&avg, dx, dy, 0.0,  data, mrcPixelMag);
 198                                 }
 199                         }
 200                         DEBUGPRINT("StartingPoint\n")
 201                         newx.data[0] = x.data[0];
 202                         newy.data[0] = y.data[0];
 203
 204                         flagOK = 0;
 205                         for(i=1; i<m; i++) {
 206                                 DEBUGPRINT("Temporary Image Get\n");
 207                                 for(dy=0; dy<tmp.HeaderN.y; dy++) {
 208                                         lVectorSpline2D((dy-1.0)*linfo->deltay+i*tmp.HeaderN.y,
 209                                                 &(px[0]), &(py[0]), &oldp, &oldx, &oldy, &olda, &oldb);
 210                                         lVectorSpline2D( dy     *linfo->deltay+i*tmp.HeaderN.y,
 211                                                 &(px[1]), &(py[1]), &oldp, &oldx, &oldy, &olda, &oldb);
 212                                         lVectorSpline2D((dy+1.0)*linfo->deltay+i*tmp.HeaderN.y,
 213                                                 &(px[2]), &(py[2]), &oldp, &oldx, &oldy, &olda, &oldb);
 214                                         vy = -(px[2] - px[0])/2.0;
 215                                         vx =  (py[2] - py[0])/2.0;
 216                                         l  = sqrt(SQR(vx) + SQR(vy));
 217                                         vx /= l;
 218                                         vy /= l;
 219                                         for(dx=0; dx<tmp.HeaderN.x; dx++) {
 220                                                 s  = (dx - tmp.HeaderN.x/2.0)*linfo->deltax;
 221                                                 sx = px[1] + vx*s;
 222                                                 sy = py[1] + vy*s;
 223                                                 mrcPixelDataGet(in,   sx, sy, 0.0, &data, mrcPixelMag, mrcPixelHowNearest);
 224                                                 mrcPixelDataSet(&tmp, dx, dy, 0.0,  data, mrcPixelMag);
 225                                         }
 226                                 }
 227
 228                                 /* newx etc. */
 229                                 lmrcImageCorrelation(&cor, &tmp, &ref, 16);
 230                                 info.mode = meanOfAll;
 231                                 lmrcImageInformation(&info, &cor);
 232                                 DEBUGPRINT3("%d: %f %f\n", i, info.maxCoord.x, info.maxCoord.y);
 233                                 if(info.maxCoord.y < cor.HeaderN.y - info.maxCoord.y) {
 234                                         dy =                   info.maxCoord.y;
 235                                 } else {
 236                                         dy = -(cor.HeaderN.y - info.maxCoord.y);
 237                                 }
 238                                 lVectorSpline2D((dy-1.0)*linfo->deltay+i*tmp.HeaderN.y, &(px[0]), &(py[0]), &oldp, &oldx, &oldy, &olda, &oldb);
 239                                 lVectorSpline2D((dy    )*linfo->deltay+i*tmp.HeaderN.y, &(px[1]), &(py[1]), &oldp, &oldx, &oldy, &olda, &oldb);
 240                                 lVectorSpline2D((dy+1.0)*linfo->deltay+i*tmp.HeaderN.y, &(px[2]), &(py[2]), &oldp, &oldx, &oldy, &olda, &oldb);
 241                                 vy = -(px[2] - px[0])/2.0;
 242                                 vx =  (py[2] - py[0])/2.0;
 243                                 l  = sqrt(SQR(vx) + SQR(vy));
 244                                 vx /= l;
 245                                 vy /= l;
 246                                 DEBUGPRINT3("%d: %f %f\n", i, info.maxCoord.x, info.maxCoord.y);
 247                                 if(info.maxCoord.x < cor.HeaderN.x - info.maxCoord.x) {
 248                                         s =  (info.maxCoord.x                )*linfo->deltax;
 249                                 } else {
 250                                         s = -(cor.HeaderN.x - info.maxCoord.x)*linfo->deltax;
 251                                 }
 252                                 newx.data[i] = px[1] + vx*s;
 253                                 newy.data[i] = py[1] + vy*s;
 254
 255                                 DEBUGPRINT3("%d: %f %f\n", i, info.maxCoord.x, info.maxCoord.y);
 256                                 DEBUGPRINT3("%d: old %f %f\n", i, px[1], py[1]);
 257                                 DEBUGPRINT3("%d: new %f %f\n", i, newx.data[i], newy.data[i]);
 258                                 DEBUGPRINT3("%d: v   %f %f\n", i, vx, vy);
 259                                 DEBUGPRINT2("%d: s   %f \n", i, s);
 260                                 flagOK += (dy + s);
 261
 262                                 shift.x = -info.maxCoord.x;
 263                                 shift.y = -info.maxCoord.y;
 264                                 shift.z = 0;
 265                                 lmrcImageShift(&tmp2, &tmp, shift, 0);
 266                                 lmrcImageAdd(&avg, &tmp, &num);
 267                                 {
 268                                         char s[1024];
 269                                         sprintf(s, "/tmp/tmp.test.%d", i);
 270                                         mrcFileWrite(&tmp, s, "", 0);
 271                                 }
 272                         }
 273                         if(0 == flagOK) {
 274                                 flagIter=-1;
 275                         } else {
 276                                 flagIter++;
 277                         }
 278                         lVectorSplineTable2DMake(&newp, &newx, &newy, &newa, &newb);
 279                         fprintf(stdout, "mrcImageUnbentROI -r ");
 280                         for(i=0; i<m; i++) {
 281                                 fprintf(stdout, "%f %f ", newx.data[i], newy.data[i]);
 282                         }
 283                         fprintf(stdout, "\n");
 284                         lmrcImageDevidedByReal(&avg, (double)m);
 285                         mrcFileWrite(&avg, "/tmp/avg.test", "", 0);
 286                 }
 287                 p = newp;
 288                 x = newx;
 289                 y = newy;
 290                 a = newa;
 291                 b = newb;
 292         }
 293
 294
 295         for(dy=0; dy<out->HeaderN.y; dy++) {
 296                 lVectorSpline2D((dy-1.0)*linfo->deltay, &(px[0]), &(py[0]), &p, &x, &y, &a, &b);
 297                 lVectorSpline2D( dy     *linfo->deltay, &(px[1]), &(py[1]), &p, &x, &y, &a, &b);
 298                 lVectorSpline2D((dy+1.0)*linfo->deltay, &(px[2]), &(py[2]), &p, &x, &y, &a, &b);
 299
 300                 vy = -(px[2] - px[0])/2.0;
 301                 vx =  (py[2] - py[0])/2.0;
 302                 l  = sqrt(SQR(vx) + SQR(vy));
 303                 vx /= l;
 304                 vy /= l;
 305                 for(dx=0; dx<out->HeaderN.x; dx++) {
 306                         s = (dx - out->HeaderN.x/2.0)*linfo->deltax;
 307                         sx = px[1] + vx*s;
 308                         sy = py[1] + vy*s;
 309                         mrcPixelDataGet(in,  sx, sy, 0.0, &data, mrcPixelMag, mrcPixelHowNearest);
 310                         mrcPixelDataSet(out, dx, dy, 0.0,  data, mrcPixelMag);
 311                 }
 312         }
 313 }
 314
 315 void
 316 additionalUsage()
 317 {
 318 }