imgtransform.c

Go to the documentation of this file.

 
/*****************************************************************************/
#include "libtpcimgp.h"
/*****************************************************************************/
 
/*****************************************************************************/
int img2cube(
  IMG *img1,
  int dim,
  IMG *img2
) {
  int ret, n, testf;
  int xi, yi, zi, xj, yj, zj, fi, xp, yp, zp;
  double xsize, ysize, zsize, newsize;
  double xdist, ydist, zdist;
  double xpc, ypc, zpc, xd, yd, zd;
  double w, wsum, d;
 
  if(IMG_TEST>0) printf("img2cube(img1, img2)\n"); 
  if(img1->status!=IMG_STATUS_OCCUPIED) return(1);
  if(img2->status==IMG_STATUS_OCCUPIED) imgEmpty(img2);
 
  if(IMG_TEST>1)
    printf("original dimensions (x,y,z) := %d,%d,%d\n", img1->dimx, img1->dimy, img1->dimz);
  if(dim<1) {
    /* Get the largest dimension of original image volume */
    dim=img1->dimz; n=0;
    if(img1->dimx>dim) {dim=img1->dimx; n++;}
    if(img1->dimy>dim) {dim=img1->dimy; n++;}
  }
  if(IMG_TEST>1) printf("new dimensions (x,y,z) := %d,%d,%d\n", dim, dim, dim);
  testf=dim*dim*dim/10;
 
  /* Allocate memory for the new IMG */
  ret=imgAllocate(img2, dim, dim, dim, img1->dimt);
  if(ret) return(10+ret);
  imgCopyhdr(img1, img2);
 
  /* Calculate the physical extensions of original image volume */
  xsize=img1->dimx*(img1->sizex + img1->gapx); 
  ysize=img1->dimy*(img1->sizey + img1->gapy); 
  zsize=img1->dimz*(img1->sizez + img1->gapz);
  if(IMG_TEST>1) printf("original image size (x,y,z) in mm := %g,%g,%g\n", xsize, ysize, zsize);
  /* New image volume will be a cube of highest size of those */
  newsize=xsize; n=0;
  if(ysize>newsize) {newsize=ysize; n++;}
  if(zsize>newsize) {newsize=zsize; n++;}
  if(IMG_TEST>1) {
    printf("original image size (x,y,z) in mm := %g,%g,%g\n",
      xsize, ysize, zsize);
    if(n>0) printf("new image size (x,y,z) in mm:= %g,%g,%g\n",
      newsize, newsize, newsize);
    else printf("original image size needs not to be changed.\n");
  }
  img2->gapx=img2->gapy=img2->gapz=0.0;
  img2->sizex=img2->sizey=img2->sizez=newsize/(double)dim;
 
  /* Resample the image volume */
  n=0;
  for(zj=0; zj<dim; zj++) for(yj=0; yj<dim; yj++) for(xj=0; xj<dim; xj++) {
    if(IMG_TEST>3 && n%testf==0) printf("zj=%d yj=%d xj=%d\n", zj, yj, xj);
    /* Calculate the distance of this voxel from the mid of new image volume */
    xdist=((double)xj+0.5)*img2->sizex - 0.5*newsize;
    ydist=((double)yj+0.5)*img2->sizey - 0.5*newsize;
    zdist=((double)zj+0.5)*img2->sizez - 0.5*newsize;
    if(IMG_TEST>3 && n%testf==0)
      printf("  xdist=%g ydist=%g zdist=%g\n", xdist, ydist, zdist);
    /* The place in coordinates of original image */
    xpc=(0.5*xsize+xdist)/(img1->sizex + img1->gapx);
    ypc=(0.5*ysize+ydist)/(img1->sizey + img1->gapy);
    zpc=(0.5*zsize+zdist)/(img1->sizez + img1->gapz);
    /* Inside which voxel it resides? */
    xp=(int)(xpc+0.5); yp=(int)(ypc+0.5); zp=(int)(zpc+0.5);
    if(IMG_TEST>3 && n%testf==0)
      printf("  inside pixel %d,%d,%d\n", xp, yp, zp);
    /* Calculate 1/distance-weighted average of 3x3x3 pixels around it */
    for(fi=0; fi<img1->dimt; fi++) img2->m[zj][yj][xj][fi]=0.0;
    wsum=0.0;
    for(zi=zp-1; zi<=zp+1; zi++)
      for(yi=yp-1; yi<=yp+1; yi++) for(xi=xp-1; xi<=xp+1; xi++) {
        /* Distance between voxels */
        xd=(0.5+xi)-xpc; yd=(0.5+yi)-ypc; zd=(0.5+zi)-zpc;
        d=xd*xd+yd*yd+zd*zd; //if(d>0.0) d=sqrt(d); else d=0.0;
        if(IMG_TEST>3 && n%testf==0)
          printf("    distance^2 from (%d,%d,%d) := %g\n", xi, yi, zi, d);
        /* Add weight to the weight sum */
        w=exp(-d); wsum+=w;
        if(IMG_TEST>3 && n%testf==0) printf("    weight := %g\n", w);
        /* If voxel is inside the image, add value to the sum */
        if(zi>=0 && yi>=0 && xi>=0 && zi<img1->dimz && yi<img1->dimy && xi<img1->dimx) {
          for(fi=0; fi<img1->dimt; fi++) {
            img2->m[zj][yj][xj][fi]+=w*img1->m[zi][yi][xi][fi];
          }
        }
      }
    /* Divide the sum with sum weight */
    for(fi=0; fi<img1->dimt; fi++) img2->m[zj][yj][xj][fi]/=wsum;
    if(IMG_TEST>3 && n%testf==0)
      printf("  weighted avg in 1st frame := %g\n", img2->m[zj][yj][xj][0]);
    n++;
  }
 
  return 0;
}
/*****************************************************************************/
 
/*****************************************************************************/
void imgScale(
  IMG *src,
  IMG *targ,
  float zoom,
  int method
) {
  int frame, plane, x, y;
  float **tmp_targ;
 
  if(method==0) {} // prevent compiler warning
 
  // Allocate memory for temporary target matrix:
  tmp_targ = malloc(targ->dimy * sizeof(float*));
  for(y=0; y<targ->dimy; y++)
    tmp_targ[y] = malloc(targ->dimx * sizeof(float));
 
  for(plane=0; plane<src->dimz; plane++)
    for(frame=0; frame<src->dimt; frame++){
      integerScale(frame, src->m[plane], tmp_targ, src->dimx, src->dimy, (int)roundf(zoom));
      // Copy scaled image matrix to target image buffer:
      for(y=0; y<targ->dimy; y++)
        for(x=0; x<targ->dimx; x++)
          targ->m[plane][y][x][frame] = tmp_targ[y][x]; 
    }
}
/*****************************************************************************/
 
/*****************************************************************************/
void integerScale(
  int frame,
  float ***src, 
  float **targ, 
  int width,
  int height,
  int zoom
) {
  if(frame<0) frame=0;
  if(zoom<1) zoom=1;
 
  int w, h, z;
 
  for(h=0; h<height; h++){
    for(w=0; w<width; w++){
      float val=src[h][w][frame];
      for(z=0; z<zoom; z++) {
        targ[h*zoom][w*zoom + z]=val;
      }
    }
    /* Copy to next rows */
    for(z=1; z<zoom; z++) {
      for(int i=0; i<zoom*width; i++) targ[h*zoom+z][i]=targ[h*zoom+z-1][i];
      //memcpy(targ[h*zoom+z], targ[h*zoom], width*zoom*sizeof(float));
    }
  }
}
/*****************************************************************************/
 
/*****************************************************************************/
int imgSwell(
  IMG *img1,
  IMG *img2,
  int doz
) {
  if(img1==NULL || img2==NULL) return(1);
  if(img1->dimz<1 || img1->dimy<1 || img1->dimx<1 || img1->dimt<1) return(2);
 
  /* Set correct dimensions for the swollen image */
  int dimz, dimy, dimx;
  if(doz!=0) dimz=2*img1->dimz; else dimz=img1->dimz;
  dimy=2*img1->dimy; 
  dimx=2*img1->dimx; 
 
  /* If output image does not yet have these dimensions, then re-allocate it and set headers */
  if(img2->dimz!=dimz || img2->dimy!=dimy || img2->dimx!=dimx || img2->dimt!=img1->dimt) {
    /* Allocate memory for the image */
    int ret=imgAllocateWithHeader(img2, dimz, dimy, dimx, img1->dimt, img1); 
    if(ret) return(10+ret);
    /* Correct the pixel sizes so that the actual size (cm x cm) of the image won't change */
    if(doz!=0) {img2->sizez=0.5*img1->sizez; img2->gapz=0.5*img1->gapz;}
    img2->sizey=0.5*img1->sizey; img2->gapy=0.5*img1->gapy;
    img2->sizex=0.5*img1->sizex; img2->gapx=0.5*img1->gapx;
  }
 
  /* Fill the output image matrix */
  if(doz!=0) {
    int x, y, z;
    for(z=0; z<img1->dimz; z++) for(y=0; y<img1->dimy; y++) for(x=0; x<img1->dimx; x++) {
      int nz, ny, nx;
      for(nz=2*z; nz<2+2*z; nz++) for(ny=2*y; ny<2+2*y; ny++) for(nx=2*x; nx<2+2*x; nx++) {
        for(int t=0; t<img1->dimt; t++) img2->m[nz][ny][nx][t]=img1->m[z][y][x][t];
      }
    }
  } else {
    int x, y, z;
    for(z=0; z<img1->dimz; z++) for(y=0; y<img1->dimy; y++) for(x=0; x<img1->dimx; x++) {
      int ny, nx;
      for(ny=2*y; ny<2+2*y; ny++) for(nx=2*x; nx<2+2*x; nx++) {
        for(int t=0; t<img1->dimt; t++) img2->m[z][ny][nx][t]=img1->m[z][y][x][t];
      }
    }
  }
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/
int imgShrink(
  IMG *img1,
  IMG *img2,
  int doz
) {
  if(img1==NULL || img2==NULL) return(1);
  if(img1->dimz<1 || img1->dimy<2 || img1->dimx<2 || img1->dimt<1) return(2);
 
  /* Set correct dimensions for the shrunken image */
  int dimz, dimy, dimx;
  if(doz!=0) dimz=0.5*img1->dimz; else dimz=img1->dimz;
  dimy=0.5*img1->dimy; 
  dimx=0.5*img1->dimx; 
  /* Correct dimensions (add 1) if the original ones were odd */
  if((img1->dimz % 2!=0.0) && doz!=0) dimz+=1;
  if(img1->dimy % 2!=0.0) dimy+=1;
  if(img1->dimx % 2!=0.0) dimx+=1;
 
  /* If output image does not yet have these dimensions, then re-allocate it and set headers */
  if(img2->dimz!=dimz || img2->dimy!=dimy || img2->dimx!=dimx || img2->dimt!=img1->dimt) {
    /* Allocate memory for the image */
    int ret=imgAllocateWithHeader(img2, dimz, dimy, dimx, img1->dimt, img1); 
    if(ret) return(10+ret);
    /* Correct the pixel sizes so that the actual size (cm x cm) of the image won't change */
    if(doz!=0) {img2->sizez=2.0*img1->sizez; img2->gapz=2.0*img1->gapz;}
    img2->sizey=2.0*img1->sizey; img2->gapy=2.0*img1->gapy;
    img2->sizex=2.0*img1->sizex; img2->gapx=2.0*img1->gapx;
  }
 
  /* Calculate the averages of 2x2x2 or 2x2 pixel areas */
  if(doz!=0) {
    int fi, pi, yi, xi, pc, yc, xc;
    for(fi=0; fi<img2->dimt; fi++){
      for(pi=0, pc=1; pi<img2->dimz; pi++) {
        for(yi=0, yc=1; yi<img2->dimy; yi++) for(xi=0, xc=1; xi<img2->dimx; xi++) {
          img2->m[pi][yi][xi][fi]=0.0;
          if(2*pi+1==img1->dimz) pc=0;
          if(2*yi+1==img1->dimy) yc=0;
          if(2*xi+1==img1->dimx) xc=0;
          img2->m[pi][yi][xi][fi]+=img1->m[2*pi][2*yi][2*xi][fi];
          img2->m[pi][yi][xi][fi]+=img1->m[2*pi][2*yi][2*xi+xc][fi];
          img2->m[pi][yi][xi][fi]+=img1->m[2*pi][2*yi+yc][2*xi][fi];
          img2->m[pi][yi][xi][fi]+=img1->m[2*pi][2*yi+yc][2*xi+xc][fi];
          img2->m[pi][yi][xi][fi]+=img1->m[2*pi+pc][2*yi][2*xi][fi];
          img2->m[pi][yi][xi][fi]+=img1->m[2*pi+pc][2*yi][2*xi+xc][fi];
          img2->m[pi][yi][xi][fi]+=img1->m[2*pi+pc][2*yi+yc][2*xi][fi];
          img2->m[pi][yi][xi][fi]+=img1->m[2*pi+pc][2*yi+yc][2*xi+xc][fi];
          img2->m[pi][yi][xi][fi]/=8.0;
        } /* next row and column */
      } /* next plane */
    } /*next frame */
  } else {
    int fi, pi, yi, xi, yc, xc;
    for(fi=0; fi<img2->dimt; fi++){
      for(pi=0; pi<img2->dimz; pi++) {
        for(yi=0, yc=1; yi<img2->dimy; yi++) for(xi=0, xc=1; xi<img2->dimx; xi++) {
          img2->m[pi][yi][xi][fi]=0.0;
          if(2*yi+1==img1->dimy) yc=0;
          if(2*xi+1==img1->dimx) xc=0;
          img2->m[pi][yi][xi][fi]+=img1->m[pi][2*yi][2*xi][fi];
          img2->m[pi][yi][xi][fi]+=img1->m[pi][2*yi][2*xi+xc][fi];
          img2->m[pi][yi][xi][fi]+=img1->m[pi][2*yi+yc][2*xi][fi];
          img2->m[pi][yi][xi][fi]+=img1->m[pi][2*yi+yc][2*xi+xc][fi];
          img2->m[pi][yi][xi][fi]/=4.0;
        } /* next row and column */
      } /* next plane */
    } /*next frame */
  }
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/