home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 27 / IOPROG_27.ISO / SOFT / CODICE / COMPR.ZIP / DEC.C next >
Encoding:
C/C++ Source or Header  |  1997-07-16  |  21.6 KB  |  569 lines
  1. /**************************************************************************/
  2. /* Decode an image encoded with a quadtree partition based fractal scheme */
  3. /*                                                                        */
  4. /*       Copyright 1993,1994 Yuval Fisher. All rights reserved.           */
  5. /*                                                                        */
  6. /* Version 0.05 10/10/94                                                  */
  7. /**************************************************************************/
  8.  
  9. /* The following belong in a encdec.h file, but nevermind...              */
  10. /* -----------------------------------------------------------------------*/
  11. #include <stdio.h>
  12. #include <math.h>
  13.  
  14. #define DEBUG 0
  15. #define GREY_LEVELS 255
  16.  
  17. #define IMAGE_TYPE unsigned char /* may be different in some applications */
  18.  
  19. /* The following #define allocates an hsize x vsize  matrix of type type */
  20. #define matrix_allocate(matrix, hsize,vsize, TYPE) {\
  21.     TYPE *imptr; \
  22.     int _i; \
  23.     matrix = (TYPE **)malloc(vsize*sizeof(TYPE *));\
  24.     imptr = (TYPE*)malloc((long)hsize*(long)vsize*sizeof(TYPE));\
  25.     if (imptr == NULL) \
  26.         fatal("\nNo memory in matrix allocate."); \
  27.     for (_i = 0; _i<vsize; ++_i, imptr += hsize) \
  28.         matrix[_i] = imptr; \
  29.  }
  30.  
  31. #define bound(a)   ((a) < 0.0 ? 0 : ((a)>255.0? 255 : a))
  32.  
  33. /* various function declarations to keep compiler warnings away           */
  34. void fatal();
  35. char *malloc();
  36. char *strcpy(); 
  37.  
  38. /* ---------------------------------------------------------------------- */
  39.  
  40. IMAGE_TYPE **image,*imptr,**image1; 
  41.                              /* The input image data and a dummy         */ 
  42.  
  43. double max_scale = 1.0;      /* maximum allowable grey level scale factor */
  44.  
  45. int    s_bits = 5,           /* number of bits used to store scale factor */
  46.        o_bits = 7,           /* number of bits used to store offset       */
  47.        hsize = -1,           /* The horizontal size of the input image    */
  48.        vsize = -1,           /* The vertical size                         */
  49.        scaledhsize,          /* hsize*scalefactor                         */
  50.        scaledvsize,          /* vsize*scalefactor                         */
  51.        size,                 /* largest square image that fits in image   */
  52.        min_part = 3,         /* min and max _part determine a range of    */
  53.        max_part = 4,         /* Range sizes from hsize>>min to hsize>>max */
  54.        dom_step = 4,         /* Density of domains relative to size       */
  55.        dom_step_type = 0,    /* Flag for dom_step a multiplier or divisor */
  56.        dom_type = 0,         /* Method of generating domain pool 0,1,2..  */
  57.        post_process = 1,     /* Flag for postprocessing.                  */
  58.        num_iterations= 10,   /* Number of decoding iterations used.       */
  59.        *no_h_domains,        /* Number of horizontal domains.             */
  60.        *domain_hstep,        /* Domain density step size.                 */
  61.        *domain_vstep,        /* Domain density step size.                 */
  62.        *bits_needed,         /* Number of bits to encode domain position. */
  63.        zero_ialpha,          /* the const ialpha when alpha = 0           */
  64.        output_partition=0,   /* A flag for outputing the partition        */
  65.        max_exponent;         /* The max power of 2 side of square image   */
  66.                              /* that fits in our input image.             */
  67.  
  68. struct transformation_node {
  69.        int rx,ry,            /* The range position and size in a trans.   */
  70.               xsize, ysize,  
  71.               rrx,rry,
  72.               dx,dy;         /* The domain position.                      */
  73.        int sym_op;           /* The symmetry operation used in the trans. */
  74.        int depth;            /* The depth in the quadtree partition.      */
  75.        double scale, offset; /* scalling and offset values.               */
  76.        struct transformation_node *next;       /* The next trans. in list */
  77. } transformations, *trans; 
  78.  
  79. FILE *input,*output,*other_input; 
  80.  
  81.  
  82. /* fatal is for when there is a fatal error... print a message and exit */
  83. void fatal(s)
  84. char *s;
  85. {
  86.      printf("\n%s\n",s);
  87.      exit(-1);
  88. }
  89.  
  90. /* ************************************************************ */
  91. /* unpack value using size bits read from fin.                  */
  92. /* ************************************************************ */
  93. long unpack(size, fin)
  94. int size;
  95. FILE *fin;
  96. {
  97.     int i;
  98.     int value = 0;
  99.     static int ptr = 1; /* how many bits are packed in sum so far */
  100.     static int sum; 
  101.  
  102.  
  103.     /* size == -2 means we initialize things */
  104.     if (size == -2) {
  105.         sum = fgetc(fin);
  106.         sum <<= 1;
  107.         return((long)0);
  108.     }
  109.  
  110.     /* size == -1 means we want to peek at the next bit without */
  111.     /* advancing the pointer */
  112.     if (size == -1)
  113.         return((long)((sum&256)>>8));
  114.  
  115.      for (i=0; i<size; ++i, ++ptr,  sum <<= 1) {
  116.         if (sum & 256) value |= 1<<i;
  117.  
  118.          if (ptr == 8) {
  119.             sum = getc(fin);
  120.             ptr=0;
  121.         }
  122.      }
  123.     return((long)value);
  124. }
  125.  
  126. main(argc,argv)
  127. int argc;
  128. char **argv;
  129. {
  130.     /* Defaults are set initially */
  131.     double         scalefactor = 1.0;           /* Scale factor for output */
  132.     char           inputfilename[200];
  133.     char           outputfilename[200];
  134.     char           other_input_file[200];
  135.     int            i,j, x_exponent, y_exponent;
  136.     int            domain_size, no_domains;
  137.  
  138.  
  139.     inputfilename[0] = 1;    /* We initially set the input to this and */
  140.     outputfilename[0] = 1;   /* then check if the input/output names   */
  141.     other_input_file[0] = 1; /* have been set below.                   */
  142.  
  143.     /* scan through the input line and read in the arguments */
  144.     for (i=1; i<argc; ++i) 
  145.        if (argv[i][0] != '-' ) 
  146.             if (inputfilename[0] == 1)
  147.                 strcpy(inputfilename, argv[i]);
  148.             else if (outputfilename[0] == 1)
  149.                 strcpy(outputfilename, argv[i]);
  150.             else;
  151.        else { /* we have a flag */
  152.             if (strlen(argv[i]) == 1) break;
  153.             switch(argv[i][1]) {
  154.                case 'i': strcpy(other_input_file,argv[++i]); 
  155.                          break;
  156.                case 'n': num_iterations = atoi(argv[++i]);
  157.                          break;
  158.                case 'f': scalefactor = atof(argv[++i]);
  159.                          break;
  160.                case 'P': output_partition = 1;
  161.                          break;
  162.                case 'p': post_process = 0;
  163.                          break;
  164.                case 's': s_bits = atoi(argv[++i]);
  165.                          break;
  166.                case 'o': o_bits = atoi(argv[++i]);
  167.                          break;
  168.                case 'N': max_scale = atof(argv[++i]);
  169.                          break;
  170.                case '?':
  171.                case 'H':
  172.                default:
  173.            printf("\nUsage: dec -[options] [inputfile [outputfile]]");
  174.            printf("\nOptions are: (# = number), defaults show in ()");
  175.            printf("\n  -f # scale factor of output size. (%lf)", scalefactor);
  176.            printf("\n  -i file name. An initial image to iteration from.");
  177.            printf("\n  -n # no. of decoding iterations. (%d)", num_iterations);
  178.            printf("\n  -N # maximum allowed scaling. (%lf)",max_scale);
  179.            printf("\n  -s # number of scaling quantizing bits. (%d)",s_bits);
  180.            printf("\n  -o # number of offset quantizing bits. (%d)",o_bits);
  181.            printf("\n  -P   output the partition into the output file. (off)");
  182.            printf("\n  -p   supress artifact postprocessing. (off)");
  183.            fatal(" ");
  184.             }
  185.        }
  186.  
  187.     if (inputfilename[0] == 1) strcpy(inputfilename, "lenna.trn");
  188.     if (outputfilename[0] == 1) strcpy(outputfilename, "lenna.out");
  189.  
  190.     if ((input = fopen(inputfilename, "r")) == NULL)
  191.         fatal("Can't open input file.");
  192.         
  193.     unpack(-2,input); /* initialize the unpacking routine */
  194.  
  195.     /* read the header data from the input file. This should probably */
  196.     /* be put into one read which reads a structure with the info     */
  197.     min_part = (int)unpack(4,input);
  198.     max_part = (int)unpack(4,input);
  199.     dom_step = (int)unpack(4,input);
  200.     dom_step_type = (int)unpack(1,input);
  201.     dom_type = (int)unpack(2,input);
  202.     hsize = (int)unpack(12,input);
  203.     vsize = (int)unpack(12,input);
  204.  
  205.     /* we now compute size */
  206.     x_exponent = (int)floor(log((double)hsize)/log(2.0));
  207.     y_exponent = (int)floor(log((double)vsize)/log(2.0));
  208.    
  209.     /* exponent is min of x_ and y_ exponent */
  210.     max_exponent = (x_exponent > y_exponent ? y_exponent : x_exponent);
  211.     /* size is the size of the largest square that fits in the image */
  212.     /* It is used to compute the domain and range sizes.             */
  213.     size =  1<<max_exponent; 
  214.  
  215.     /* This is the quantized value of zero scaling */
  216.     zero_ialpha = 0.5 + (max_scale)/(2.0*max_scale)*(1<<s_bits);
  217.  
  218.     /* allocate memory for the output image. Allocating one chunck saves  */
  219.     /* work and time later.                                              */
  220.     scaledhsize = (int)(scalefactor*hsize);
  221.     scaledvsize = (int)(scalefactor*vsize);
  222.     matrix_allocate(image, scaledhsize,scaledvsize, IMAGE_TYPE);
  223.     matrix_allocate(image1, scaledhsize, scaledvsize, IMAGE_TYPE);
  224.  
  225.     if (other_input_file[0] != 1) {
  226.         other_input = fopen(other_input_file, "r");
  227.         i = fread(image[0], sizeof(IMAGE_TYPE), 
  228.                 scaledhsize*scaledvsize, other_input);
  229.         if (i < scaledhsize*scaledvsize) 
  230.            fatal("Couldn't read input... not enough data.");
  231.         else
  232.             printf("\n%d pixels read from %s.\n", i,other_input_file);
  233.         fclose(other_input);
  234.     }
  235.           
  236.     /* since max_ and min_ part are variable, these must be allocated */
  237.     i = max_part - min_part + 1;
  238.     bits_needed = (int *)malloc(sizeof(int)*i);
  239.     no_h_domains = (int *)malloc(sizeof(int)*i);
  240.     domain_hstep = (int *)malloc(sizeof(int)*i);
  241.     domain_vstep = (int *)malloc(sizeof(int)*i);
  242.     
  243.     /* compute bits needed to read each domain type */
  244.     for (i=0; i <= max_part-min_part; ++i) {
  245.        /* first compute how many domains there are horizontally */
  246.        domain_size = size >> (min_part+i-1);
  247.        if (dom_type == 2) 
  248.              domain_hstep[i] = dom_step; 
  249.        else if (dom_type == 1)
  250.              if (dom_step_type ==1)
  251.                 domain_hstep[i] = (size >> (max_part - i-1))*dom_step;
  252.              else 
  253.                 domain_hstep[i] = (size >> (max_part - i-1))/dom_step;
  254.        else 
  255.              if (dom_step_type ==1)
  256.                 domain_hstep[i] = domain_size*dom_step;
  257.              else 
  258.                 domain_hstep[i] = domain_size/dom_step;
  259.  
  260.        no_h_domains[i] = 1+(hsize-domain_size)/domain_hstep[i];
  261.        /* bits_needed[i][0] = ceil(log(no_domains)/log(2));  */
  262.  
  263.        /* now compute how many domains there are vertically */
  264.        if (dom_type == 2) 
  265.              domain_vstep[i] = dom_step; 
  266.        else if (dom_type == 1)
  267.              if (dom_step_type ==1)
  268.              domain_vstep[i] = (size >> (max_part - i-1))*dom_step;
  269.              else 
  270.              domain_vstep[i] = (size >> (max_part - i-1))/dom_step;
  271.        else 
  272.              if (dom_step_type ==1)
  273.                 domain_vstep[i] = domain_size*dom_step;
  274.              else 
  275.                 domain_vstep[i] = domain_size/dom_step;
  276.  
  277.        no_domains = 1+(vsize-domain_size)/domain_vstep[i];
  278.        bits_needed[i] = ceil(log((double)no_domains*(double)no_h_domains[i])/
  279.                            log(2.0)); 
  280.      }
  281.  
  282.     if ((output = fopen(outputfilename, "w")) == NULL)
  283.             fatal("Can't open output file.");
  284.  
  285.     /* Read in the transformation data */
  286.     trans = &transformations;
  287.     printf("\nReading transformations.....");
  288.     fflush(stdout);
  289.     partition_image(0, 0, hsize,vsize );
  290.     fclose(input);
  291.     printf("Done.");
  292.     fflush(stdout);
  293.    
  294.     if (scalefactor != 1.0) {
  295.         printf("\nScaling image to %d x %d.", scaledhsize,scaledvsize);
  296.         scale_transformations(scalefactor); 
  297.     }
  298.  
  299.     /* when we output the partition, we just read the transformations */
  300.     /* in and write them to the outputfile                            */ 
  301.     if (output_partition) {
  302.           fprintf(output,"\n%d %d\n %d %d\n%d %d\n\n", 
  303.                   0, 0, scaledhsize, 0, scaledhsize, scaledvsize);
  304.           printf("\nOutputed partition data in %s\n",outputfilename);
  305.           fclose(output);
  306.           return;
  307.     }
  308.  
  309.     for (i=0; i<num_iterations; ++i) 
  310.       apply_transformations();
  311.  
  312.     if (post_process)     
  313.       smooth_image();
  314.  
  315.     i = fwrite(image[0], sizeof(IMAGE_TYPE), scaledhsize*scaledvsize, output);
  316.     if (i < scaledhsize*scaledvsize) 
  317.         fatal("Couldn't write output... not enough disk space ?.");
  318.     else
  319.         printf("\n%d pixels written to output file.\n", i);
  320.  
  321.     fclose(output);
  322. }
  323.  
  324. /* ************************************************************ */
  325. /* Read in the transformation data from *input.                 */
  326. /* This is a recursive routine whose recursion tree follows the */
  327. /* recursion done by the encoding program.                      */
  328. /* ************************************************************ */
  329. read_transformations(atx,aty,xsize,ysize,depth)
  330. int atx,aty,xsize,ysize,depth; 
  331. {
  332.     /* Having all these locals in a recursive procedure is hard on the  */
  333.     /* stack.. but it is more readable.                                 */
  334.     int i,j,
  335.         sym_op,                   /* A symmetry operation of the square */
  336.         ialpha,                   /* Intgerized scalling factor         */
  337.         ibeta;                    /* Intgerized offset                  */
  338.  
  339.     long domain_ref; 
  340.  
  341.     double alpha, beta;
  342.  
  343.     /* keep breaking it down until we are small enough */
  344.     if (depth < min_part) {
  345.          read_transformations(atx,aty, xsize/2, ysize/2, depth+1);
  346.          read_transformations(atx+xsize/2,aty, xsize/2, ysize/2, depth+1);
  347.          read_transformations(atx,aty+ysize/2,xsize/2, ysize/2,  depth+1);
  348.          read_transformations(atx+xsize/2,aty+ysize/2,xsize/2,ysize/2,depth+1);
  349.          return;
  350.     }
  351.  
  352.     if (depth < max_part && unpack(1,input)) {
  353.          /* A 1 means we subdivided.. so quadtree */
  354.          read_transformations(atx,aty, xsize/2, ysize/2, depth+1);
  355.          read_transformations(atx+xsize/2,aty, xsize/2, ysize/2, depth+1);
  356.          read_transformations(atx,aty+ysize/2, xsize/2, ysize/2, depth+1);
  357.          read_transformations(atx+xsize/2,aty+ysize/2,xsize/2,ysize/2,depth+1);
  358.     } else {  
  359.          /* we have a transformation to read */
  360.          trans->next = (struct transformation_node *)
  361.                   malloc(sizeof(struct transformation_node ));
  362.          trans = trans->next; trans->next = NULL;
  363.          ialpha = (int)unpack(s_bits,  input);
  364.          ibeta = (int)unpack(o_bits,  input);
  365.          alpha = (double)ialpha/(double)(1<<s_bits)*(2.0*max_scale)-max_scale;
  366.  
  367.          beta = (double)ibeta/(double)((1<<o_bits)-1)*
  368.                ((1.0+fabs(alpha))*GREY_LEVELS);
  369.          if (alpha > 0.0) beta -= alpha*GREY_LEVELS;
  370.  
  371.          trans->scale = alpha;
  372.          trans->offset = beta;
  373.          if (ialpha != zero_ialpha) {
  374.             trans-> sym_op = (int)unpack(3, input);
  375.             domain_ref = unpack(bits_needed[depth-min_part], input);
  376.             trans->dx = (double)(domain_ref % no_h_domains[depth-min_part])
  377.                                               * domain_hstep[depth-min_part];
  378.             trans->dy = (double)(domain_ref / no_h_domains[depth-min_part])
  379.                                               * domain_vstep[depth-min_part];
  380.          } else {
  381.             trans-> sym_op = 0;
  382.             trans-> dx  = 0;
  383.             trans-> dy = 0;
  384.          }
  385.          trans->rx = atx;
  386.          trans->ry = aty;
  387.          trans->depth = depth;
  388.         
  389.          trans->rrx = atx + xsize;
  390.          trans->rry = aty + ysize;
  391.  
  392.          if (output_partition) 
  393.                 fprintf(output,"\n%d %d\n %d %d\n%d %d\n\n", 
  394.                   atx,       vsize-aty-ysize, 
  395.                   atx,       vsize-aty, 
  396.                   atx+xsize, vsize-aty);
  397.   
  398.     }
  399. }
  400.      
  401. /* ************************************************************ */
  402. /* Apply the transformations once to an initially black image.  */
  403. /* ************************************************************ */
  404. apply_transformations()
  405. {
  406.     IMAGE_TYPE **tempimage;
  407.     int    i,j,i1,j1,count=0;
  408.     double pixel;
  409.  
  410.     trans = &transformations;
  411.     while (trans->next != NULL) {
  412.         trans = trans->next;
  413.         ++count;
  414.  
  415. /* Since the inner loop is the same in each case of the switch below */
  416. /* we just define it once for easy modification.                     */
  417. #define COMPUTE_LOOP              {                                  \
  418.         pixel = (image[j1][i1]+image[j1][i1+1]+image[j1+1][i1]+      \
  419.                 image[j1+1][i1+1])/4.0;                              \
  420.         image1[j][i] = bound(0.5 + trans->scale*pixel+trans->offset);\
  421.         }
  422.  
  423.         switch(trans->sym_op) {
  424.            case 0: for (i=trans->rx, i1 = trans->dx; 
  425.                          i<trans->rrx; ++i, i1 += 2)
  426.                    for (j=trans->ry, j1 = trans->dy; 
  427.                          j<trans->rry; ++j, j1 += 2) 
  428.                        COMPUTE_LOOP
  429.                    break;
  430.            case 1: for (j=trans->ry, i1 = trans->dx;
  431.                            j<trans->rry; ++j, i1 += 2)
  432.                    for (i=trans->rrx-1,
  433.                            j1 = trans->dy; i>=(int)trans->rx; --i, j1 += 2) 
  434.                        COMPUTE_LOOP
  435.                    break;
  436.            case 2: for (i=trans->rrx-1, 
  437.                           i1 = trans->dx; i>=(int)trans->rx; --i, i1 += 2)
  438.                    for (j=trans->rry-1,
  439.                            j1 = trans->dy; j>=(int)trans->ry; --j, j1 += 2) 
  440.                        COMPUTE_LOOP
  441.                    break;
  442.            case 3: for (j=trans->rry-1,
  443.                            i1 = trans->dx; j>=(int)trans->ry; --j, i1 += 2)
  444.                    for (i=trans->rx, j1 = trans->dy;
  445.                            i<trans->rrx; ++i, j1 += 2) 
  446.                        COMPUTE_LOOP
  447.                    break;
  448.            case 4: for (j=trans->ry, i1 = trans->dx;
  449.                            j<trans->rry; ++j, i1 += 2)
  450.                    for (i=trans->rx, j1 = trans->dy;
  451.                            i<trans->rrx; ++i, j1 += 2) 
  452.                        COMPUTE_LOOP
  453.                    break;
  454.            case 5: for (i=trans->rx, i1 = trans->dx;
  455.                            i<trans->rrx; ++i, i1 += 2)
  456.                    for (j=trans->rry-1,
  457.                            j1 = trans->dy; j>=(int)trans->ry; --j, j1 += 2) 
  458.                        COMPUTE_LOOP
  459.                    break;
  460.            case 6: for (j=trans->rry-1,
  461.                            i1 = trans->dx; j>=(int)trans->ry; --j, i1 += 2)
  462.                    for (i=trans->rrx-1,
  463.                            j1 = trans->dy; i>=(int)trans->rx; --i, j1 += 2) 
  464.                        COMPUTE_LOOP
  465.                    break;
  466.            case 7: for (i=trans->rrx-1, 
  467.                           i1 = trans->dx; i>=(int)trans->rx; --i, i1 += 2)
  468.                    for (j=trans->ry, j1 = trans->dy; 
  469.                           j<trans->rry; ++j, j1 += 2) 
  470.                        COMPUTE_LOOP
  471.                    break;
  472.         }  
  473.   
  474.     }
  475.     tempimage = image;
  476.     image = image1;
  477.     image1 = tempimage;
  478.  
  479.     printf("\n%d transformations applied.",count);
  480. }
  481.  
  482. /*        This should really be done when they are read in.     */
  483. /* ************************************************************ */
  484. scale_transformations(scalefactor)
  485. double scalefactor;
  486. {
  487.     trans = &transformations;
  488.     while (trans->next != NULL) {
  489.         trans = trans->next;
  490.  
  491.         trans->rrx *= scalefactor;
  492.         trans->rry *= scalefactor;
  493.         trans->rx *= scalefactor;
  494.         trans->ry *= scalefactor;
  495.         trans->dx *= scalefactor;
  496.         trans->dy *= scalefactor;
  497.     }
  498. }
  499.  
  500. /* ************************************************************* */
  501. /* Recursively partition an image, finding the largest contained */
  502. /* square and call read_transformations .                        */
  503. /* ************************************************************* */
  504. partition_image(atx, aty, hsize,vsize )
  505. int atx, aty, hsize,vsize;
  506. {
  507.    int x_exponent,    /* the largest power of 2 image size that fits */
  508.        y_exponent,    /* horizontally or vertically the rectangle.   */
  509.        exponent,      /* The actual size of image that's encoded.    */
  510.        size, 
  511.        depth; 
  512.    
  513.    x_exponent = (int)floor(log((double)hsize)/log(2.0));
  514.    y_exponent = (int)floor(log((double)vsize)/log(2.0));
  515.    
  516.    /* exponent is min of x_ and y_ exponent */
  517.    exponent = (x_exponent > y_exponent ? y_exponent : x_exponent);
  518.    size = 1<<exponent; 
  519.    depth = max_exponent - exponent;
  520.  
  521.    read_transformations(atx,aty,size,size,depth);
  522.  
  523.    if (size != hsize) 
  524.       partition_image(atx+size, aty, hsize-size,vsize );
  525.  
  526.    if (size != vsize) 
  527.       partition_image(atx, aty+size, size,vsize-size );
  528. }        
  529.   
  530. /* ************************************************************* */
  531. /* Scan the image and average the transformation boundaries.     */
  532. /* ************************************************************* */
  533. smooth_image()
  534. {
  535.     IMAGE_TYPE pixel1, pixel2;
  536.     int i,j;
  537.     int w1,w2;
  538.  
  539.     printf("\nPostprocessing Image.");
  540.     trans = &transformations;
  541.     while (trans->next != NULL) {
  542.         trans = trans->next;
  543.         if (trans->rx == 0 || trans->ry == 0) 
  544.             continue;
  545.         
  546.         if (trans->depth == max_part) {
  547.             w1 = 5;
  548.             w2 = 1;
  549.         } else {
  550.             w1 = 2;
  551.             w2 = 1;
  552.         }
  553.  
  554.         for (i=trans->rx; i<trans->rrx; ++i) {
  555.              pixel1 = image[(int)trans->ry][i];
  556.              pixel2 = image[(int)trans->ry-1][i];
  557.              image[(int)trans->ry][i] = (w1*pixel1 + w2*pixel2)/(w1+w2);
  558.              image[(int)trans->ry-1][i] = (w2*pixel1 + w1*pixel2)/(w1+w2);
  559.         }
  560.  
  561.         for (j=trans->ry; j<trans->rry; ++j) {
  562.              pixel1 = image[j][(int)trans->rx];
  563.              pixel2 = image[j][(int)trans->rx-1];
  564.              image[j][(int)trans->rx] = (w1*pixel1 + w2*pixel2)/(w1+w2);
  565.              image[j][(int)trans->rx-1] = (w2*pixel1 + w1*pixel2)/(w1+w2);
  566.         }
  567.     }
  568. }
  569.