home *** CD-ROM | disk | FTP | other *** search
/ CD/PC Actual 11 / CDACTUAL11.iso / cdactual / demobin / share / W95 / FRACCOMP / FRAC.CPP < prev    next >
Encoding:
C/C++ Source or Header  |  1995-04-30  |  59.4 KB  |  1,687 lines
  1. #include "stdafx.h"           // Access to Microsoft fundation classes
  2. #include "io.h"
  3. #include "math.h" 
  4. #include "stdlib.h"
  5. #include "frac.h"
  6. #include <fcntl.h>
  7. #include <errno.h>
  8. #include <sys/types.h>
  9. #include <sys/stat.h>
  10.  
  11. //////////////////////////////////////////////////////////////////////
  12. // CFracComp
  13. CFracComp::CFracComp(HWND hwnd,short id,BOOL output)
  14. {
  15.        
  16.          //Set display output device
  17.        m_hwnd = hwnd;
  18.        m_id = id;
  19.        m_output = output;
  20. }    
  21.  
  22.  
  23. CFracComp::~CFracComp()
  24. {
  25.      // Delete DC context
  26.      DeleteDC(m_MemDc);
  27. }                     
  28.    
  29. void CFracComp::Reset()
  30. {
  31.  
  32.        m_hsize = -1;         
  33.        m_vsize = -1;
  34.        m_subclass_search = 0,      
  35.        m_fullclass_search = 0,      
  36.        m_only_positive = 0;       
  37.        m_dom_step_type = 0;   
  38.        m_dom_type = 0;  
  39.        m_max_scale = 1.0;
  40.        m_no_h_domains = NULL;        
  41.        m_domain_hstep = NULL;        
  42.        m_domain_vstep = NULL;        
  43.        m_bits_needed  = NULL;      
  44.        m_output_partition=FALSE;
  45.        m_image = NULL;
  46.        m_imageDummy = NULL;
  47.        
  48.         // The class_transform gives the transforms between classification numbers for      
  49.         // negative scaling values, when brightest becomes darkest, etc...                   
  50.        m_class_transform[0][0] = 23;
  51.        m_class_transform[0][1] = 17;
  52.        m_class_transform[0][2] = 21;
  53.        m_class_transform[0][3] = 11;
  54.        m_class_transform[0][4] = 15;
  55.        m_class_transform[0][5] =  9;
  56.        m_class_transform[0][6] = 22;
  57.        m_class_transform[0][7] = 16;
  58.        m_class_transform[0][8] = 19;
  59.        m_class_transform[0][9] =  5;
  60.        m_class_transform[0][10] = 13;
  61.        m_class_transform[0][11] =  3;
  62.        m_class_transform[0][12] = 20;
  63.        m_class_transform[0][13] = 10;
  64.        m_class_transform[0][14] = 18;
  65.        m_class_transform[0][15] =  4;
  66.        m_class_transform[0][16] =  7;
  67.        m_class_transform[0][17] =  1;
  68.        m_class_transform[0][18] = 14;
  69.        m_class_transform[0][19] =  8;
  70.        m_class_transform[0][20] = 12;
  71.        m_class_transform[0][21] =  2;
  72.        m_class_transform[0][22] =  6;
  73.        m_class_transform[0][23] =  0;
  74.        m_class_transform[1][0] =  16;
  75.        m_class_transform[1][1] =  22;
  76.        m_class_transform[1][2] =  10;
  77.        m_class_transform[1][3] =  20;
  78.        m_class_transform[1][4] =   8;
  79.        m_class_transform[1][5] =  14;
  80.        m_class_transform[1][6] =  17;
  81.        m_class_transform[1][7] =  23;
  82.        m_class_transform[1][8] =   4;
  83.        m_class_transform[1][9] =  18;
  84.        m_class_transform[1][10] =  2;
  85.        m_class_transform[1][11] = 12;
  86.        m_class_transform[1][12] = 11;
  87.        m_class_transform[1][13] = 21;
  88.        m_class_transform[1][14] =  5;
  89.        m_class_transform[1][15] = 19;
  90.        m_class_transform[1][16] =  0;
  91.        m_class_transform[1][17] =  6;
  92.        m_class_transform[1][18] =  9;
  93.        m_class_transform[1][19] = 15;
  94.        m_class_transform[1][20] =  3;
  95.        m_class_transform[1][21] = 13;
  96.        m_class_transform[1][22] =  1;
  97.        m_class_transform[1][23] =  7;
  98.  
  99.        // rot_transform gives the rotations for domains with negative scalings.           
  100.        m_rot_transform[0][0] = 7;
  101.        m_rot_transform[0][1] = 4;
  102.        m_rot_transform[0][2] = 5;
  103.        m_rot_transform[0][3] = 6;
  104.        m_rot_transform[0][4] = 1;
  105.        m_rot_transform[0][5] = 2;
  106.        m_rot_transform[0][6] = 3;
  107.        m_rot_transform[0][7] = 0;
  108.        m_rot_transform[1][0] = 2;    
  109.        m_rot_transform[1][1] = 3;
  110.        m_rot_transform[1][2] = 0;
  111.        m_rot_transform[1][3] = 1;
  112.        m_rot_transform[1][4] = 6;
  113.        m_rot_transform[1][5] = 7;
  114.        m_rot_transform[1][6] = 4;
  115.        m_rot_transform[1][7] = 5;
  116.  
  117. }
  118.  
  119. void CFracComp::SetFileName(CString szfileIn, CString szfileOut)
  120. {
  121.      //Set input file name
  122.     m_szNameIn = szfileIn; 
  123.     //Set output file name    
  124.     m_szNameOut= szfileOut;    
  125. }
  126.  
  127. void CFracComp::Message(CString message)
  128. {
  129.   if(!m_output)
  130.       return;
  131.   if(m_hwnd != NULL)
  132.   {
  133.           if(m_id < 0)
  134.            {
  135.             HDC nDC = GetDC(m_hwnd);
  136.             TEXTMETRIC lpMetrics; 
  137.             GetTextMetrics(nDC,&lpMetrics);
  138.             TextOut(nDC,0,abs(m_id--)*lpMetrics.tmHeight,message,message.GetLength());
  139.         }
  140.     else
  141.         SetDlgItemText(m_hwnd,m_id,message);
  142.   }
  143.   printf("%s\n",message);
  144. }
  145.  
  146. BOOL CFracComp::WriteRawFile(int hsz, int vsz, CString szOut)
  147. {
  148.     int i,j;
  149.     CString msg;
  150.     // Open output file (decompression file)
  151.     if ((m_FileOut = fopen(szOut, "wb")) == NULL)
  152.     {
  153.     
  154.         Message("Error: Can't open output file");
  155.         return FALSE;
  156.  
  157.     }
  158.     // Write file to disk, a byte at a time 
  159.     // (Bug in _write function in VC 2.0) 
  160.     // force to write a byte at a time
  161.      for(i =0; i<hsz;i++)
  162.         for(j=0;j<vsz;j++)
  163.             if((fputc(m_image[i][j],m_FileOut)) != m_image[i][j])
  164.             {
  165.                 msg = "Error:  Can't write byte to disk (Disk may be full)";
  166.                   Message(msg);
  167.                 return FALSE;
  168.             
  169.             }
  170.  
  171.  
  172.     msg.Format("%d pixels written to %s.", hsz*vsz, szOut);
  173.     Message(msg);
  174.     fclose(m_FileOut);
  175.     return TRUE;
  176. }
  177.  
  178. BOOL CFracComp::ReadRawFile(int hsz, int vsz,CString szIn)
  179. {
  180.      int i,handle;
  181.      CString msg;
  182.     if((handle= _open(szIn,_O_BINARY)) == -1)
  183.     {
  184.         if(handle == EACCES)
  185.             msg = "Error: Tried to open read-only file for writing, or file's sharing mode does not allow specified operations, or given path is directory";
  186.         if(handle == EMFILE)
  187.             msg = "Error: No more file handles available (too many open files)";
  188.         if(handle == ENOENT)
  189.             msg = "Error:  File or path not found";
  190.         if(msg.IsEmpty())
  191.             msg = "Error: Unknown Error has occured while open file";
  192.           Message(msg);
  193.         return FALSE;
  194.     }
  195.     for(i=0;i<hsz; i++)
  196.         if((_read(handle,m_image[i],vsz)) != vsz)
  197.         {
  198.             msg = "Error: File is corrupted, can't read blocks";
  199.               Message(msg);
  200.             return FALSE;    
  201.         } 
  202.  
  203.     
  204.     msg.Format("%dx%d = %d pixels read from %s.", hsz,vsz,hsz*vsz,szIn);
  205.     Message(msg);
  206.     _close(handle);
  207.     return TRUE;
  208. }
  209.  
  210. HDC CFracComp::HandleDC()
  211. {
  212.     // Handle to  current DC context. If 
  213.     // Context is given this will have final
  214.     // image
  215.     return(m_MemDc);    
  216. }
  217.  
  218. HBITMAP CFracComp::HandleBitmap()
  219. {
  220.     return(m_bitmap);
  221. }
  222.  
  223. void CFracComp::Display(HWND pDest,int hsz,int vsz,CString stitle)
  224. {
  225.      int color;
  226.     SetWindowText(pDest,stitle);
  227.     for (int dx = 0 ; dx < hsz; dx++) 
  228.         for (int cx = 0 ; cx < vsz ;cx++ ) 
  229.           {
  230.            color = m_image[dx][cx];
  231.            SetPixelV(m_MemDc,cx,dx,m_colorTable[color]); 
  232.           }
  233.  
  234.     BitBlt(GetDC(pDest),0,0,hsz,vsz,m_MemDc,0,0,SRCCOPY);
  235. }
  236.  
  237. ///////////////////////////////////////////////////////////////////////////////////////
  238. // Public Encoder and Decoder
  239. ///////////////////////////////////////////////////////////////////////////////////////
  240. void CFracComp::Compress(HWND hwndIn, double tolerance, int hsize, int vsize, int max_part, int min_part,
  241. int dom_step, int dom_step_type,int scaling_bit, int offset_bit, double max_scale, int divisor, BOOL positive,
  242. BOOL d24, BOOL d3)
  243. {
  244.        
  245.     int  i,j,k;
  246.     // Reset all values
  247.     Reset();
  248.     m_hsize = hsize;
  249.     m_vsize = vsize;
  250.     m_min_part = min_part;         
  251.     m_max_part = max_part; 
  252.     m_dom_step = dom_step;
  253.     m_s_bits = scaling_bit;
  254.     m_o_bits = offset_bit;
  255.     m_max_scale = max_scale;
  256.     m_dom_step_type = dom_step_type;
  257.     m_dom_step_type = divisor;
  258.     m_only_positive = positive;
  259.     m_subclass_search = d24;
  260.     m_fullclass_search = d3; 
  261.     CString msg;
  262.     
  263.     //Check for empty fileString string;
  264.     if((m_szNameIn.IsEmpty()) || (m_szNameOut.IsEmpty()))
  265.     {
  266.          
  267.         Message("Error: Input/Output name not given");
  268.         return;
  269.     }
  270.     // Start wait cursor
  271.     if(hwndIn != NULL)
  272.         BeginWaitCursor(); 
  273.      
  274.     // allocate memory for the input image. Allocating one chunck saves  
  275.     // work and time later.                                              
  276.     unsigned char *row_image; 
  277.     m_image = (unsigned char **)malloc((m_vsize)*sizeof(unsigned char *));
  278.     row_image = (unsigned char *)malloc((long)(m_hsize)*(long)(m_vsize)*sizeof(unsigned char));
  279.     if (row_image == NULL)
  280.     {
  281.     
  282.           Message("Error: Out of memory from image buffer");
  283.         return;
  284.     }
  285.     for (i = 0; i<m_vsize; ++i, row_image += m_hsize) 
  286.         m_image[i] = row_image;
  287.  
  288.     double *row_domimage0;
  289.        m_domimage[0] = (double **)malloc((m_vsize/2)*sizeof(double *));
  290.     row_domimage0 = (double *)malloc((long)(m_hsize/2)*(long)(m_vsize/2)*sizeof(double));
  291.     if (row_domimage0 == NULL)
  292.     {
  293.           Message("Error: Out of memory for domain buffer");
  294.         return;
  295.     }
  296.     for (i = 0; i<m_vsize/2; ++i, row_domimage0 += m_hsize/2) 
  297.         m_domimage[0][i] = row_domimage0;
  298.  
  299.     double *row_domimage1;
  300.        m_domimage[1] = (double **)malloc((m_vsize/2)*sizeof(double *));
  301.     row_domimage1 = (double *)malloc((long)(m_hsize/2)*(long)(m_vsize/2)*sizeof(double));
  302.     if (row_domimage1 == NULL)
  303.     {
  304.         
  305.         msg = "Error: Out of memory for domain buffer";
  306.           Message(msg);
  307.         return;
  308.     }
  309.     for (i = 0; i<m_vsize/2; ++i, row_domimage1 += m_hsize/2) 
  310.         m_domimage[1][i] = row_domimage1;
  311.  
  312.     double *row_domimage2;
  313.        m_domimage[2] = (double **)malloc((m_vsize/2)*sizeof(double *));
  314.     row_domimage2 = (double *)malloc((long)(m_hsize/2)*(long)(m_vsize/2)*sizeof(double));
  315.     if (row_domimage2 == NULL)
  316.     {
  317.         msg = "Error: Out of memory for domain buffer";
  318.           Message(msg);
  319.         return;
  320.     }
  321.     for (i = 0; i<m_vsize/2; ++i, row_domimage2 += m_hsize/2) 
  322.         m_domimage[2][i] = row_domimage2;
  323.  
  324.     double *row_domimage3;
  325.        m_domimage[3] = (double **)malloc((m_vsize/2)*sizeof(double *));
  326.     row_domimage3 = (double *)malloc((long)(m_hsize/2)*(long)(m_vsize/2)*sizeof(double));
  327.     if (row_domimage3 == NULL)
  328.     {
  329.         msg = "Error: Out of memory for domain buffer";
  330.           Message(msg);
  331.         return;
  332.     }
  333.     for (i = 0; i<m_vsize/2; ++i, row_domimage3 += m_hsize/2) 
  334.         m_domimage[3][i] = row_domimage3;
  335.  
  336.  
  337.     // max_ & min_ part are variable, so this must be run time allocated
  338.     m_bits_needed = (int *)malloc(sizeof(int)*(1+m_max_part-m_min_part));
  339.  
  340.    
  341.  
  342.     if(!ReadRawFile(hsize,vsize,m_szNameIn))
  343.         return;        
  344.  
  345.     // allcate memory for domain data and initialize it
  346.     if(!Compute_Sums(m_hsize,m_vsize))
  347.         return;
  348.  
  349.     if ((m_FileOut = fopen(m_szNameOut, "wb")) == NULL)
  350.     {
  351.            msg = "Error: Can't open output file";
  352.           Message(msg);
  353.         return;
  354.     }
  355.  
  356.     // output some data into the outputfile.                      
  357.     Pack(4,(long)m_min_part);
  358.     Pack(4,(long)m_max_part);
  359.     Pack(4,(long)m_dom_step);
  360.     Pack(1,(long)m_dom_step_type);
  361.     Pack(2,(long)m_dom_type);
  362.     Pack(12,(long)m_hsize);
  363.     Pack(12,(long)m_vsize);
  364.  
  365.     // This is the quantized value of zero scaling.. needed later 
  366.     m_zero_ialpha = 0.5 + (m_max_scale)/(2.0*m_max_scale)*(1<<m_s_bits);
  367.  
  368.     // The following routine takes a rectangular image and calls the
  369.     // quadtree routine to encode square sum-images in it.          
  370.     // the tolerance is a parameter since in some applications different 
  371.     // regions of the image may need to be compressed to different tol's 
  372.         
  373.     Message("Encoding image.....");
  374.     if(hwndIn != NULL)
  375.     {
  376.         // Delete old context device (just in case)
  377.         DeleteDC(m_MemDc);                        
  378.         // Create compatible with display device         
  379.         m_MemDc =  CreateCompatibleDC(GetDC(hwndIn));
  380.         CRect rect;
  381.         GetWindowRect(hwndIn,rect);    
  382.         // Set Window size to picture size;
  383.         SetWindowPos(hwndIn,HWND_TOP,rect.left,rect.top,m_vsize,m_hsize,SWP_SHOWWINDOW);
  384.         m_bitmap = CreateCompatibleBitmap(GetDC(hwndIn),m_vsize,m_hsize); 
  385.         SelectObject(m_MemDc,m_bitmap);
  386.         for (i= 0; i<255; i++)
  387.             m_colorTable[i] = GetNearestColor(m_MemDc,RGB(i,i,i));
  388.         ShowWindow(hwndIn,SW_SHOW);
  389.         Display(hwndIn,m_hsize,m_hsize,CString("Loading image....")); 
  390.         MoveToEx(m_MemDc,(0+m_hsize/2),0,NULL);
  391.         LineTo(m_MemDc,(0+m_hsize/2),m_vsize);
  392.         MoveToEx(m_MemDc,0,(0+m_hsize/2),NULL);
  393.         LineTo(m_MemDc,m_hsize,(0+m_vsize/2));    
  394.         SetWindowText(hwndIn,"Quadtree starting....");
  395.         BitBlt(GetDC(hwndIn),0,0,m_hsize,m_vsize,m_MemDc,0,0,SRCCOPY); 
  396.     } 
  397.  
  398.     Partition_Image(0, 0, m_hsize,m_vsize, tolerance, hwndIn);
  399.     
  400.     Message("Done");
  401.  
  402.     // stuff the last byte if needed 
  403.     Pack(-1,(long)0);
  404.     
  405.     fclose(m_FileOut);
  406.     i = Pack(-2,(long)0);
  407.     msg.Format("Compression = %lf from %d bytes written in %s.",
  408.             (double)(m_hsize*m_vsize)/(double)i, i, m_szNameOut);
  409.           Message(msg);
  410.  
  411.     // Free allocated memory
  412.     free(m_bits_needed);
  413.     free(m_domimage[0]);
  414.     free(m_domimage[1]);
  415.     free(m_domimage[2]);
  416.     free(m_domimage[3]);
  417.     free(m_domain.no_h_domains);
  418.     free(m_domain.no_v_domains);
  419.     free(m_domain.domain_hsize);
  420.     free(m_domain.domain_vsize);
  421.     free(m_domain.domain_hstep);
  422.     free(m_domain.domain_vstep);
  423.     for (i=0; i <= m_max_part-m_min_part; ++i) 
  424.         free(m_domain.pixel[i]);
  425.     free(m_domain.pixel);
  426.     free(m_image[0]);
  427.     // free memory allocated in the list structure the_domain
  428.     struct classified_domain *node;
  429.     struct classified_domain *tempnode;
  430.     for (i=0; i <= m_max_part-m_min_part; ++i) 
  431.         for (k=0; k<3; ++k) 
  432.              for (j=0; j<24; ++j) 
  433.              {
  434.                 node =m_dd_domain[k][j][i]; 
  435.                 while(node->next != NULL)
  436.                {
  437.                        tempnode = node;
  438.                     node = node->next;
  439.                     free(tempnode);
  440.                 }
  441.                  free(node);
  442.               }
  443.     SetWindowText(hwndIn,"Done.");
  444.     // End wait cursor
  445.     if(hwndIn != NULL)
  446.         EndWaitCursor();  
  447.    
  448. }  
  449.  
  450. void CFracComp::Decompress(HWND hwndOut,int scale, int iterations, BOOL process,
  451.     double maximum_scale, int scaling_bit, int offset_bit, BOOL output_p)
  452. {
  453.      
  454.    
  455.     //Counter variables
  456.     int            i;                       
  457.     int               x_exponent, y_exponent;
  458.     int            domain_size,no_domains;
  459.     int               scaledvsize,scaledhsize; 
  460.     CString msg;
  461.     // Reset all values
  462.     Reset();
  463.     //Scaling factor
  464.     m_scaledfactor = scale;
  465.     m_min_part = 3;         
  466.     m_max_part = 4; 
  467.     m_dom_step = 4;  
  468.     m_output_partition = output_p;
  469.     m_s_bits = scaling_bit;
  470.     m_o_bits = offset_bit;
  471.     m_max_scale = maximum_scale;
  472.     //Check for empty fileString string;
  473.     if((m_szNameIn.IsEmpty()) || (m_szNameOut.IsEmpty()))
  474.     {
  475.         
  476.           Message("Error: Input/Output name not given");
  477.         return;
  478.     }               
  479.     // Open compression file
  480.     if ((m_FileIn = fopen(m_szNameIn, "rb")) == NULL)
  481.     {
  482.    
  483.           Message("Error: Can't open Input file");
  484.         return;
  485.     }
  486.     if(hwndOut != NULL)
  487.         BeginWaitCursor();      
  488.     Unpack(-2); // initialize the unpacking routine 
  489.  
  490.     // read the header data from the input file. This should probably 
  491.     // be put into one read which reads a structure with the info     
  492.     m_min_part = (int)Unpack(4);
  493.     m_max_part = (int)Unpack(4);
  494.     m_dom_step = (int)Unpack(4);
  495.     m_dom_step_type = (int)Unpack(1);
  496.     m_dom_type = (int)Unpack(2);
  497.     m_hsize = (int)Unpack(12);
  498.     m_vsize = (int)Unpack(12);
  499.  
  500.     // we now compute size 
  501.     x_exponent = (int)floor(log((double)m_hsize)/log(2.0));
  502.     y_exponent = (int)floor(log((double)m_vsize)/log(2.0));
  503.    
  504.     // exponent is min of x_ and y_ exponent 
  505.     m_max_exponent = (x_exponent > y_exponent ? y_exponent : x_exponent);
  506.     // size is the size of the largest square that fits in the image 
  507.     // It is used to compute the domain and range sizes.             
  508.     m_size =  1<<m_max_exponent; 
  509.  
  510.     // This is the quantized value of zero scaling
  511.     m_zero_ialpha = 0.5 + (m_max_scale)/(2.0*m_max_scale)*(1<<m_s_bits);
  512.  
  513.     // allocate memory for the output image. Allocating one chunck saves  
  514.     // work and time later.                                              
  515.     scaledhsize = (int)(m_scaledfactor*m_hsize);
  516.     scaledvsize = (int)(m_scaledfactor*m_vsize);
  517.     if(hwndOut != NULL)
  518.     {
  519.         // Delete old context device (just in case)
  520.         DeleteDC(m_MemDc);                        
  521.         // Create compatible with display device         
  522.         m_MemDc =  CreateCompatibleDC(GetDC(hwndOut));
  523.         CRect rect;
  524.         GetWindowRect(hwndOut,rect);
  525.         // Set Window size to picture size;
  526.         SetWindowPos(hwndOut,HWND_TOP,rect.left,rect.top,scaledvsize,scaledhsize,SWP_SHOWWINDOW);
  527.         m_bitmap = CreateCompatibleBitmap(GetDC(hwndOut),scaledvsize,scaledhsize); 
  528.         SelectObject(m_MemDc,m_bitmap);
  529.         // Create a valid color table for this device
  530.         for (i= 0; i<255; i++)
  531.             m_colorTable[i] = GetNearestColor(m_MemDc,RGB(i,i,i));
  532.     }
  533.     // Allocate memory for image and dummy image.
  534.     unsigned char *row_image;
  535.     m_image = (unsigned char **)malloc(scaledvsize*sizeof(unsigned char *));
  536.     row_image = (unsigned char*)malloc((long)scaledhsize*(long)scaledvsize*sizeof(unsigned char));
  537.     if (row_image == NULL) 
  538.     {
  539.     
  540.           Message("Error: Out of memory for image buffer");
  541.         return;
  542.     }
  543.     for (i = 0; i<scaledvsize; ++i, row_image += scaledhsize) 
  544.         m_image[i] = row_image; 
  545.      
  546.     unsigned char *row_imageDummy;
  547.     m_imageDummy = (unsigned char **)malloc(scaledvsize*sizeof(unsigned char *));
  548.     row_imageDummy = (unsigned char *)malloc((long)scaledhsize*(long)scaledvsize*sizeof(unsigned char));
  549.     if (row_imageDummy == NULL) 
  550.     {
  551.     
  552.           Message("Error: Out of memory for image secondary buffer");
  553.         return;
  554.     }
  555.     for (i = 0; i<scaledvsize; ++i, row_imageDummy += scaledhsize) 
  556.         m_imageDummy[i] = row_imageDummy; 
  557.  
  558.  
  559.           
  560.     // since max_ and min_ part are variable, these must be allocated 
  561.     i = m_max_part - m_min_part + 1;
  562.     m_bits_needed = (int *)malloc(sizeof(int)*i);
  563.     m_no_h_domains = (int *)malloc(sizeof(int)*i);
  564.     m_domain_hstep = (int *)malloc(sizeof(int)*i);
  565.     m_domain_vstep = (int *)malloc(sizeof(int)*i);
  566.     
  567.     // compute bits needed to read each domain type 
  568.     for (i=0; i <= m_max_part-m_min_part; ++i) {
  569.        // first compute how many domains there are horizontally 
  570.        domain_size = m_size >> (m_min_part+i-1);
  571.        if (m_dom_type == 2) 
  572.              m_domain_hstep[i] = m_dom_step; 
  573.        else if (m_dom_type == 1)
  574.              if (m_dom_step_type ==1)
  575.                 m_domain_hstep[i] = (m_size >> (m_max_part - i-1))*m_dom_step;
  576.              else 
  577.                 m_domain_hstep[i] = (m_size >> (m_max_part - i-1))/m_dom_step;
  578.        else 
  579.              if (m_dom_step_type ==1)
  580.                 m_domain_hstep[i] = domain_size*m_dom_step;
  581.              else 
  582.                 m_domain_hstep[i] = domain_size/m_dom_step;
  583.  
  584.        m_no_h_domains[i] = 1+(m_hsize-domain_size)/m_domain_hstep[i];
  585.       
  586.        // now compute how many domains there are vertically 
  587.        if (m_dom_type == 2) 
  588.              m_domain_vstep[i] = m_dom_step; 
  589.        else if (m_dom_type == 1)
  590.              if (m_dom_step_type ==1)
  591.              m_domain_vstep[i] = (m_size >> (m_max_part - i-1))*m_dom_step;
  592.              else 
  593.              m_domain_vstep[i] = (m_size >> (m_max_part - i-1))/m_dom_step;
  594.        else 
  595.              if (m_dom_step_type ==1)
  596.                 m_domain_vstep[i] = domain_size*m_dom_step;
  597.              else 
  598.                 m_domain_vstep[i] = domain_size/m_dom_step;
  599.  
  600.        no_domains = 1+(m_vsize-domain_size)/m_domain_vstep[i];
  601.        m_bits_needed[i] = ceil(log((double)no_domains*(double)m_no_h_domains[i])/log(2.0)); 
  602.     }
  603.     // Open output file (decompression file)
  604.     if ((m_FileOut = fopen(m_szNameOut, "wb")) == NULL)
  605.     {
  606.     
  607.         Message("Error: Can't open output file");
  608.         return;
  609.  
  610.     }
  611.  
  612.     // Read in the transformation data 
  613.     m_trans = &m_transformations;
  614.  
  615.     Message("Reading transformation....");
  616.  
  617.     Partition_Image(0, 0, m_hsize,m_vsize );
  618.     fclose(m_FileIn);
  619.    
  620.     Message("Done");
  621.  
  622.    
  623.  
  624.     // when we output the partition, we just read the transformations 
  625.     // in and write them to the outputfile                             
  626.     if (m_output_partition) 
  627.     {
  628.               // Open output file (decompression file)
  629.             if ((m_FileOut = fopen(m_szNameOut, "wb")) == NULL)
  630.             {
  631.     
  632.                     Message("Error: Can't open output file");
  633.                     return;
  634.  
  635.             }
  636.               fprintf(m_FileOut,"\n%d %d\n %d %d\n%d %d\n\n", 
  637.                   0, 0, scaledhsize, 0, scaledhsize, scaledvsize);
  638.          
  639.             msg.Format("Outputed partition data in %s.",m_szNameOut);
  640.               Message(msg);
  641.         
  642.               fclose(m_FileOut);
  643.               return;
  644.     }
  645.  
  646.     if(hwndOut != NULL)
  647.         ShowWindow(hwndOut,SW_SHOW);
  648.     // Apply transformation 
  649.     for (i=0; i<iterations; ++i) 
  650.         if(hwndOut != NULL)
  651.         {
  652.             msg.Format("Processing ... [%d%%]",(i*100)/iterations);
  653.             Apply_Transformations(hwndOut,msg);
  654.         }
  655.     else
  656.         Apply_Transformations(NULL,"");
  657.     // Fix, edges for a smoother image
  658.     if (process)     
  659.     if(hwndOut != NULL)
  660.         Smooth_Image(hwndOut);
  661.     else
  662.         Smooth_Image(NULL);
  663.     // Write Data to file.
  664.     if(!WriteRawFile(scaledhsize,scaledvsize,m_szNameOut))
  665.         return;
  666.     SetWindowText(hwndOut,"Done.");
  667.     free(m_image);
  668.     free(m_imageDummy);
  669.     if(hwndOut != NULL)
  670.         EndWaitCursor();  
  671. }
  672.  
  673. ///////////////////////////////////////////////////////////////////////////////////////
  674. // Private Decoding routines
  675. ///////////////////////////////////////////////////////////////////////////////////////
  676. void CFracComp::Smooth_Image(HWND hdcOut)
  677. {
  678.       unsigned char pixel1, pixel2;
  679.     int i,j;
  680.     int w1,w2;
  681.  
  682.         
  683.     Message("Postprocessing Image.");
  684.  
  685.     m_trans = &m_transformations;
  686.     while (m_trans->next != NULL) {
  687.         m_trans = m_trans->next;
  688.         if (m_trans->rx == 0 || m_trans->ry == 0) 
  689.             continue;
  690.         
  691.         if (m_trans->depth == m_max_part) 
  692.         {
  693.             w1 = 5;
  694.             w2 = 1;
  695.         } else {
  696.             w1 = 2;
  697.             w2 = 1;
  698.         }
  699.  
  700.         for (i=m_trans->rx; i<m_trans->rrx; ++i) {
  701.              pixel1 = m_image[(int)m_trans->ry][i];
  702.              pixel2 = m_image[(int)m_trans->ry-1][i];
  703.              m_image[(int)m_trans->ry][i] = (w1*pixel1 + w2*pixel2)/(w1+w2);
  704.              m_image[(int)m_trans->ry-1][i] = (w2*pixel1 + w1*pixel2)/(w1+w2);
  705.         }
  706.  
  707.         for (j=m_trans->ry; j<m_trans->rry; ++j) {
  708.              pixel1 = m_image[j][(int)m_trans->rx];
  709.              pixel2 = m_image[j][(int)m_trans->rx-1];
  710.              m_image[j][(int)m_trans->rx] = (w1*pixel1 + w2*pixel2)/(w1+w2);
  711.              m_image[j][(int)m_trans->rx-1] = (w2*pixel1 + w1*pixel2)/(w1+w2);
  712.         }
  713.     }
  714.     if(hdcOut != NULL)  
  715.         Display(hdcOut,m_hsize*m_scaledfactor,m_vsize*m_scaledfactor,CString("Post Processing..."));
  716. }       
  717.   
  718. void CFracComp::Partition_Image(int atx, int aty,  int hsize, int vsize )
  719. {
  720.    int x_exponent,    // the largest power of 2 image size that fits 
  721.        y_exponent,    // horizontally or vertically the rectangle.   
  722.        exponent,      // The actual size of image that's encoded.    
  723.        size, 
  724.        depth; 
  725.    
  726.    x_exponent = (int)floor(log((double)hsize)/log(2.0));
  727.    y_exponent = (int)floor(log((double)vsize)/log(2.0));
  728.    
  729.    // exponent is min of x_ and y_ exponent 
  730.    exponent = (x_exponent > y_exponent ? y_exponent : x_exponent);
  731.    size = 1<<exponent; 
  732.    depth = m_max_exponent - exponent;
  733.  
  734.    Read_Transformations(atx,aty,size,size,depth);
  735.  
  736.    if (size != hsize) 
  737.       Partition_Image(atx+size, aty, hsize-size,vsize );
  738.  
  739.    if (size != vsize) 
  740.       Partition_Image(atx, aty+size, size,vsize-size );
  741. }        
  742.  
  743.  
  744. void CFracComp::Apply_Transformations(HWND hdcOut,CString stitle)
  745. {
  746.      unsigned char **tempimage;
  747.     int    i,j,i1,j1,count=0;
  748.     double pixel;
  749.     CString msg;
  750.     m_trans = &m_transformations;
  751.     while (m_trans->next != NULL) 
  752.     {
  753.         m_trans = m_trans->next;
  754.         ++count;
  755.  
  756. // Since the inner loop is the same in each case of the switch below
  757. // we just define it once for easy modification.                     
  758.  
  759.         switch(m_trans->sym_op) 
  760.         {
  761.            case 0: for (i=m_trans->rx, i1 = m_trans->dx;i<m_trans->rrx; ++i, i1 += 2)
  762.                            for (j=m_trans->ry, j1 = m_trans->dy;j<m_trans->rry; ++j, j1 += 2) 
  763.                            {                                 
  764.                             pixel = (m_image[j1][i1]+m_image[j1][i1+1]+m_image[j1+1][i1]+m_image[j1+1][i1+1])/4.0;                              
  765.                             m_imageDummy[j][i] = bound(0.5 + m_trans->scale*pixel+m_trans->offset);
  766.                         }
  767.  
  768.                    break;
  769.  
  770.            case 1: for (j=m_trans->ry, i1 = m_trans->dx;j<m_trans->rry; ++j, i1 += 2)
  771.                            for (i=m_trans->rrx-1,j1 = m_trans->dy; i>=(int)m_trans->rx; --i, j1 += 2) 
  772.                            {                                 
  773.                             pixel = (m_image[j1][i1]+m_image[j1][i1+1]+m_image[j1+1][i1]+m_image[j1+1][i1+1])/4.0;                              
  774.                             m_imageDummy[j][i] = bound(0.5 + m_trans->scale*pixel+m_trans->offset);
  775.                         }
  776.  
  777.                    break;
  778.            case 2: for (i=m_trans->rrx-1,i1 = m_trans->dx; i>=(int)m_trans->rx; --i, i1 += 2)
  779.                            for (j=m_trans->rry-1,j1 = m_trans->dy; j>=(int)m_trans->ry; --j, j1 += 2) 
  780.                            {                                 
  781.                             pixel = (m_image[j1][i1]+m_image[j1][i1+1]+m_image[j1+1][i1]+m_image[j1+1][i1+1])/4.0;                              
  782.                             m_imageDummy[j][i] = bound(0.5 + m_trans->scale*pixel+m_trans->offset);
  783.                         }
  784.  
  785.                    break;
  786.  
  787.            case 3: for (j=m_trans->rry-1,i1 = m_trans->dx; j>=(int)m_trans->ry; --j, i1 += 2)
  788.                            for (i=m_trans->rx, j1 = m_trans->dy;i<m_trans->rrx; ++i, j1 += 2) 
  789.                            {                                 
  790.                             pixel = (m_image[j1][i1]+m_image[j1][i1+1]+m_image[j1+1][i1]+m_image[j1+1][i1+1])/4.0;                              
  791.                             m_imageDummy[j][i] = bound(0.5 + m_trans->scale*pixel+m_trans->offset);
  792.                         }
  793.  
  794.                    break;
  795.            
  796.            case 4: for (j=m_trans->ry, i1 = m_trans->dx;j<m_trans->rry; ++j, i1 += 2)
  797.                            for (i=m_trans->rx, j1 = m_trans->dy;i<m_trans->rrx; ++i, j1 += 2) 
  798.                           {                                 
  799.                             pixel = (m_image[j1][i1]+m_image[j1][i1+1]+m_image[j1+1][i1]+m_image[j1+1][i1+1])/4.0;                              
  800.                             m_imageDummy[j][i] = bound(0.5 + m_trans->scale*pixel+m_trans->offset);
  801.                         }
  802.  
  803.                    break;
  804.  
  805.            case 5: for (i=m_trans->rx, i1 = m_trans->dx;i<m_trans->rrx; ++i, i1 += 2)
  806.                            for (j=m_trans->rry-1,j1 = m_trans->dy; j>=(int)m_trans->ry; --j, j1 += 2) 
  807.                           {                                 
  808.                             pixel = (m_image[j1][i1]+m_image[j1][i1+1]+m_image[j1+1][i1]+m_image[j1+1][i1+1])/4.0;                              
  809.                             m_imageDummy[j][i] = bound(0.5 + m_trans->scale*pixel+m_trans->offset);
  810.                         }
  811.  
  812.                    break;
  813.  
  814.            case 6: for (j=m_trans->rry-1,i1 = m_trans->dx; j>=(int)m_trans->ry; --j, i1 += 2)
  815.                            for (i=m_trans->rrx-1,j1 = m_trans->dy; i>=(int)m_trans->rx; --i, j1 += 2) 
  816.                            {                                 
  817.                             pixel = (m_image[j1][i1]+m_image[j1][i1+1]+m_image[j1+1][i1]+m_image[j1+1][i1+1])/4.0;                              
  818.                             m_imageDummy[j][i] = bound(0.5 + m_trans->scale*pixel+m_trans->offset);
  819.                         }
  820.  
  821.                    break;
  822.            case 7: for (i=m_trans->rrx-1,i1 = m_trans->dx; i>=(int)m_trans->rx; --i, i1 += 2)
  823.                            for (j=m_trans->ry, j1 = m_trans->dy; j<m_trans->rry; ++j, j1 += 2) 
  824.                            {                                 
  825.                             pixel = (m_image[j1][i1]+m_image[j1][i1+1]+m_image[j1+1][i1]+m_image[j1+1][i1+1])/4.0;                              
  826.                             m_imageDummy[j][i] = bound(0.5 + m_trans->scale*pixel+m_trans->offset);
  827.                         }            
  828.  
  829.                    break;
  830.         }  
  831.   
  832.     }
  833.     tempimage = m_image;
  834.     m_image = m_imageDummy;
  835.     m_imageDummy = tempimage;
  836.  
  837.     msg.Format("%d transformations applied.",count);
  838.     Message(msg);
  839.     if(hdcOut != NULL)  
  840.         Display(hdcOut,m_hsize*m_scaledfactor,m_vsize*m_scaledfactor,stitle);
  841.     
  842. }
  843.  
  844. void CFracComp::Read_Transformations(int atx, int aty, int xsize, int ysize, int depth)
  845. {
  846.    
  847.     // Locals in a recursive procedure             
  848.     
  849.     int ialpha,                   // Intgerized scaling factor         
  850.         ibeta;                    // Intgerized offset                  
  851.  
  852.     long domain_ref; 
  853.  
  854.     double alpha, beta;
  855.  
  856.     // keep breaking it down until we are small enough 
  857.     if (depth < m_min_part) {
  858.          Read_Transformations(atx,aty, xsize/2, ysize/2, depth+1);
  859.          Read_Transformations(atx+xsize/2,aty, xsize/2, ysize/2, depth+1);
  860.          Read_Transformations(atx,aty+ysize/2,xsize/2, ysize/2,  depth+1);
  861.          Read_Transformations(atx+xsize/2,aty+ysize/2,xsize/2,ysize/2,depth+1);
  862.          return;
  863.     }
  864.  
  865.     if (depth < m_max_part && Unpack(1)) {
  866.          // A 1 means we subdivided.. so quadtree 
  867.          Read_Transformations(atx,aty, xsize/2, ysize/2, depth+1);
  868.          Read_Transformations(atx+xsize/2,aty, xsize/2, ysize/2, depth+1);
  869.          Read_Transformations(atx,aty+ysize/2, xsize/2, ysize/2, depth+1);
  870.          Read_Transformations(atx+xsize/2,aty+ysize/2,xsize/2,ysize/2,depth+1);
  871.     } else {  
  872.          // we have a transformation to read 
  873.          m_trans->next = (struct transformation_node *)
  874.                   malloc(sizeof(struct transformation_node ));
  875.          m_trans = m_trans->next;
  876.          m_trans->next = NULL;
  877.          ialpha = (int)Unpack(m_s_bits);
  878.          ibeta = (int)Unpack(m_o_bits);
  879.          alpha = (double)ialpha/(double)(1<<m_s_bits)*(2.0*m_max_scale)-m_max_scale;
  880.  
  881.          beta = (double)ibeta/(double)((1<<m_o_bits)-1)*
  882.                ((1.0+fabs(alpha))*GREY_LEVELS);
  883.          if (alpha > 0.0) beta -= alpha*GREY_LEVELS;
  884.  
  885.          m_trans->scale = alpha;
  886.          m_trans->offset = beta;
  887.          if (ialpha != m_zero_ialpha) {
  888.             m_trans-> sym_op = (int)Unpack(3);
  889.             domain_ref = Unpack(m_bits_needed[depth-m_min_part]);
  890.             m_trans->dx = (double)(domain_ref % m_no_h_domains[depth-m_min_part])
  891.                                               * m_domain_hstep[depth-m_min_part];
  892.             m_trans->dy = (double)(domain_ref / m_no_h_domains[depth-m_min_part])
  893.                                               * m_domain_vstep[depth-m_min_part];
  894.          } else {
  895.             m_trans-> sym_op = 0;
  896.             m_trans-> dx  = 0;
  897.             m_trans-> dy = 0;
  898.          }
  899.          m_trans->rx = atx;
  900.          m_trans->ry = aty;
  901.          m_trans->depth = depth;
  902.         
  903.          m_trans->rrx = atx + xsize;
  904.          m_trans->rry = aty + ysize;
  905.          //Apply scaling factor
  906.         m_trans->rrx *= m_scaledfactor;
  907.         m_trans->rry *= m_scaledfactor;
  908.         m_trans->rx *= m_scaledfactor;
  909.         m_trans->ry *= m_scaledfactor;
  910.         m_trans->dx *= m_scaledfactor;
  911.         m_trans->dy *= m_scaledfactor;
  912.  
  913.          if (m_output_partition) 
  914.                 fprintf(m_FileOut,"\n%d %d\n %d %d\n%d %d\n\n", 
  915.                   atx,       m_vsize-aty-ysize, 
  916.                   atx,       m_vsize-aty, 
  917.                   atx+xsize, m_vsize-aty);
  918.   
  919.     }
  920.  
  921. }
  922.  
  923. long CFracComp::Unpack(int size)
  924. {
  925.     int i;
  926.     int value = 0;
  927.     static int ptr = 1; // how many bits are packed in sum so far 
  928.     static int sum; 
  929.  
  930.  
  931.     // size == -2 means we initialize things 
  932.     if (size == -2) {
  933.         sum = fgetc(m_FileIn);
  934.         sum <<= 1;
  935.         return((long)0);
  936.     }
  937.  
  938.     // size == -1 means we want to peek at the next bit without
  939.     // advancing the pointer 
  940.     if (size == -1)
  941.         return((long)((sum&256)>>8));
  942.  
  943.      for (i=0; i<size; ++i, ++ptr,  sum <<= 1) {
  944.         if (sum & 256) value |= 1<<i;
  945.  
  946.          if (ptr == 8) {
  947.             sum = getc(m_FileIn);
  948.             ptr=0;
  949.         }
  950.      }
  951.     return((long)value);
  952. }
  953. ///////////////////////////////////////////////////////////////////////////////////////
  954. // Private Encoding routines
  955. ///////////////////////////////////////////////////////////////////////////////////////
  956. void CFracComp::Classify(int x, int y, int xsize, int ysize, int *pfirst, int *psecond, 
  957. int *psym, double *psum, double *psum2, int type)
  958. {
  959.  
  960.     int order[4], i,j; 
  961.     double a[4],a2[4];
  962.      
  963.      // Get the average values of each quadrant   
  964.     if (type == 2) 
  965.     {
  966.         Average(x,y,xsize/2,ysize/2,&a[0], &a2[0]);
  967.         Average(x,y+ysize/2,xsize/2,ysize/2,&a[1], &a2[1]);
  968.         Average(x+xsize/2,y+ysize/2, xsize/2,ysize/2, &a[2],&a2[2]);
  969.         Average(x+xsize/2,y,xsize/2,ysize/2,&a[3],&a2[3]);
  970.     
  971.     }
  972.     else 
  973.     {
  974.         Average1(x,y,xsize/2,ysize/2,&a[0], &a2[0]);
  975.         Average1(x,y+ysize/2,xsize/2,ysize/2,&a[1], &a2[1]);
  976.         Average1(x+xsize/2,y+ysize/2, xsize/2,ysize/2, &a[2],&a2[2]);
  977.         Average1(x+xsize/2,y,xsize/2,ysize/2,&a[3],&a2[3]);
  978.  
  979.     }
  980.                        
  981.  
  982.     *psum = a[0] + a[1] + a[2] + a[3];
  983.     *psum2 =  a2[0] + a2[1] + a2[2] + a2[3];
  984.  
  985.     for (i=0; i<4; ++i) {
  986.          // after the sorting below order[i] is the i-th brightest quadrant.                                                
  987.          order[i] = i; 
  988.          // convert a2[] to store the variance of each quadrant      
  989.          a2[i] -= (double)(1<<(2*type))*a[i]*a[i]/(double)(xsize*ysize);
  990.     }
  991.  
  992.     // Now order the average value and also in order[], which will
  993.     // then tell us the indices (in a[]) of the brightest to darkest 
  994.     for (i=2; i>=0; --i)
  995.     for (j=0; j<=i; ++j)
  996.         if (a[j]<a[j+1]) {
  997.             swap(order[j], order[j+1],int)
  998.             swap(a[j], a[j+1],double)
  999.     }
  1000.  
  1001.     // because of the way we ordered the a[] the rotation can be
  1002.     // read right off of order[]. That will make the brightest   
  1003.     // quadrant be in the upper left corner. But we must still   
  1004.     // decide which cannonical class the image portion belogs    
  1005.     // to and whether to do a flip or just a rotation. This is   
  1006.     // the following table summarizes the horrid lines below     
  1007.     // order      class            do a rotation                 
  1008.     // 0,2,1,3      0                   0                        
  1009.     // 0,2,3,1      0                   1                        
  1010.     // 0,1,2,3      1                   0                        
  1011.     // 0,3,2,1      1                   1                        
  1012.     // 0,1,3,2      2                   0                        
  1013.     // 0,3,1,2      2                   1                        
  1014.  
  1015.     *psym = order[0]; 
  1016.     // rotate the values 
  1017.     for (i=0; i<4; ++i)
  1018.         order[i] = (order[i] - (*psym) + 4)%4; 
  1019.  
  1020.     for (i=0; order[i] != 2; ++i); 
  1021.     *pfirst = i-1;
  1022.     if (order[3] == 1 || (*pfirst == 2 && order[2] == 1)) *psym += 4;
  1023.  
  1024.     //Now further classify the sub-image by the variance of its   
  1025.     // quadrants. This give 24 subclasses for each of the 3 classes 
  1026.     for (i=0; i<4; ++i) order[i] = i; 
  1027.  
  1028.     for (i=2; i>=0; --i)
  1029.     for (j=0; j<=i; ++j)
  1030.         if (a2[j]<a2[j+1]) {
  1031.             swap(order[j], order[j+1],int)
  1032.             swap(a2[j], a2[j+1],double)
  1033.     }
  1034.  
  1035.     // Now do the symmetry operation
  1036.     for (i=0; i<4; ++i)
  1037.         order[i] = (order[i] - (*psym%4) + 4)%4; 
  1038.     if (*psym > 3) 
  1039.         for (i=0; i<4; ++i)
  1040.            if (order[i]%2) order[i] = (2 + order[i])%4;
  1041.  
  1042.     // We want to return a class number from 0 to 23 depending on 
  1043.     // the ordering of the quadrants according to their variance  
  1044.     *psecond = 0;
  1045.     for (i=2; i>=0; --i)
  1046.     for (j=0; j<=i; ++j)
  1047.         if (order[j] > order[j+1]) {
  1048.             swap(order[j],order[j+1], int);
  1049.             if (order[j] == 0 || order [j+1] == 0) 
  1050.                  *psecond += 6;
  1051.             else if (order[j] == 1 || order [j+1] == 1) 
  1052.                  *psecond += 2;
  1053.             else if (order[j] == 2 || order [j+1] == 2) 
  1054.                  *psecond += 1;
  1055.         }
  1056. }
  1057.  
  1058. BOOL CFracComp::Compute_Sums(int hsize, int vsize)
  1059. {
  1060.     int i,j,k,l,
  1061.          domain_x, 
  1062.          domain_y, 
  1063.          first_class, 
  1064.          second_class,
  1065.          size,  
  1066.          x_exponent, 
  1067.          y_exponent;
  1068.  
  1069.     struct classified_domain *node;
  1070.     
  1071.     Message("Computing domain sums... ");
  1072.  
  1073.    
  1074.     fflush(stdout); 
  1075.  
  1076.     // pre-decimate the image into domimage to avoid having to    
  1077.     // do repeated averaging of 2x2 pixel groups.                   
  1078.     // There are 4 ways to decimate the image, depending on the    
  1079.     // location of the domain, odd or even address.                 
  1080.     for (i=0; i<2; ++i)
  1081.     for (j=0; j<2; ++j)
  1082.     for (k=i; k<hsize-i; k += 2)
  1083.     for (l=j; l<vsize-j; l += 2) 
  1084.         m_domimage[(i<<1)+j][l>>1][k>>1] = 
  1085.                      ((double)m_image[l][k] + (double)m_image[l+1][k+1] +
  1086.                      (double)m_image[l][k+1] + (double)m_image[l+1][k])*0.25;
  1087.  
  1088.  
  1089.     // Allocate memory for the sum and sum^2 of domain pixels       
  1090.     // We first compute the size of the largest square that fits in
  1091.     // the image.                                                   
  1092.     x_exponent = (int)floor(log((double)hsize)/log(2.0));
  1093.     y_exponent = (int)floor(log((double)vsize)/log(2.0));
  1094.    
  1095.     // exponent is min of x_ and y_ exponent 
  1096.    m_max_exponent = (x_exponent > y_exponent ? y_exponent : x_exponent);
  1097.  
  1098.     // size is the size of the largest square that fits in the image 
  1099.     // It is used to compute the domain and range sizes.             
  1100.     size =  1<<m_max_exponent; 
  1101.  
  1102.     if (m_max_exponent < m_max_part)
  1103.     {
  1104.       
  1105.           Message("Error: Partition is absurd range");
  1106.         return FALSE;
  1107.     }
  1108.     if (m_max_exponent-2 < m_max_part)
  1109.           Message("Warning: so many quadtree partitions yield absurd ranges.");
  1110.  
  1111.       
  1112.  
  1113.     i = m_max_part - m_min_part + 1;
  1114.     m_domain.no_h_domains = (int *)malloc(sizeof(int)*i);
  1115.     m_domain.no_v_domains = (int *)malloc(sizeof(int)*i);
  1116.     m_domain.domain_hsize = (int *)malloc(sizeof(int)*i);
  1117.     m_domain.domain_vsize = (int *)malloc(sizeof(int)*i);
  1118.     m_domain.domain_hstep = (int *)malloc(sizeof(int)*i);
  1119.     m_domain.domain_vstep = (int *)malloc(sizeof(int)*i);
  1120.  
  1121.     m_domain.pixel= (struct domain_pixels ***)
  1122.               malloc(i*sizeof(struct domain_pixels **));
  1123.     
  1124.     if (m_domain.pixel == NULL) 
  1125.     {
  1126.         Message("Error: Input/Output name not given");
  1127.         return FALSE;
  1128.     }
  1129.        
  1130.     for (i=0; i <= m_max_part-m_min_part; ++i) {
  1131.        // first compute how many domains there are horizontally 
  1132.        m_domain.domain_hsize[i] = size >> (m_min_part+i-1);
  1133.        if (m_dom_type == 2) 
  1134.              m_domain.domain_hstep[i] = m_dom_step; 
  1135.        else if (m_dom_type == 1)
  1136.              if (m_dom_step_type == 1) 
  1137.                m_domain.domain_hstep[i] = (size >> (m_max_part - i-1))*m_dom_step;
  1138.              else 
  1139.                m_domain.domain_hstep[i] = (size >> (m_max_part - i-1))/m_dom_step;
  1140.        else 
  1141.              if (m_dom_step_type == 1) 
  1142.                m_domain.domain_hstep[i] = m_domain.domain_hsize[i]*m_dom_step;
  1143.              else 
  1144.                m_domain.domain_hstep[i] = m_domain.domain_hsize[i]/m_dom_step;
  1145.  
  1146.        m_domain.no_h_domains[i] = 1+(hsize-m_domain.domain_hsize[i])/m_domain.domain_hstep[i];
  1147.        
  1148.        m_domain.domain_vsize[i] = size >> (m_min_part+i-1);
  1149.        if (m_dom_type == 2) 
  1150.              m_domain.domain_vstep[i] = m_dom_step; 
  1151.        else if (m_dom_type == 1)
  1152.              if (m_dom_step_type == 1) 
  1153.                m_domain.domain_vstep[i] = (size >> (m_max_part - i-1))*m_dom_step;
  1154.              else
  1155.                m_domain.domain_vstep[i] = (size >> (m_max_part - i-1))/m_dom_step;
  1156.        else 
  1157.              if (m_dom_step_type == 1) 
  1158.                m_domain.domain_vstep[i] = m_domain.domain_vsize[i]*m_dom_step;
  1159.              else
  1160.                m_domain.domain_vstep[i] = m_domain.domain_vsize[i]/m_dom_step;
  1161.  
  1162.        m_domain.no_v_domains[i] = 1+(vsize-m_domain.domain_vsize[i])/m_domain.domain_vstep[i];
  1163.  
  1164.        // Now compute the number of bits needed to store the domain data 
  1165.        m_bits_needed[i] = ceil(log((double)m_domain.no_h_domains[i]*
  1166.                                    (double)m_domain.no_v_domains[i])/log(2.0)); 
  1167.  
  1168.        struct domain_pixels *row_domain_pixel;
  1169.        m_domain.pixel[i]=  (struct domain_pixels **)malloc((m_domain.no_v_domains[i])*sizeof(struct domain_pixels *));
  1170.        row_domain_pixel = (struct domain_pixels*)malloc((long)(m_domain.no_h_domains[i])*(long)(m_domain.no_v_domains[i])*sizeof(struct domain_pixels));
  1171.        if(row_domain_pixel == NULL)
  1172.        {
  1173.               Message("Error: Out of memory for image buffer");
  1174.             return FALSE;
  1175.         }
  1176.         for (int _i = 0; _i<m_domain.no_v_domains[i]; ++_i, row_domain_pixel += m_domain.no_h_domains[i])
  1177.             m_domain.pixel[i][_i] = row_domain_pixel; 
  1178.  
  1179.  
  1180.     }
  1181.  
  1182.     // allocate and zero the list containing the classified domain data 
  1183.     i = m_max_part - m_min_part + 1;
  1184.     for (first_class = 0; first_class < 3; ++first_class)
  1185.     for (second_class = 0; second_class < 24; ++second_class) {
  1186.        m_dd_domain[first_class][second_class] = 
  1187.                            (struct classified_domain **) 
  1188.                            malloc(i*sizeof(struct classified_domain *));
  1189.        for (j=0; j<i; ++j)
  1190.                     m_dd_domain[first_class][second_class][j] = NULL;
  1191.     }
  1192.  
  1193.     // Precompute sum and sum of squares for domains                 
  1194.     // This part can get faster for overlapping domains if repeated  
  1195.     // sums are avoided                                              
  1196.     for (i=0; i <= m_max_part-m_min_part; ++i) {
  1197.       for (j=0,domain_x=0; j<m_domain.no_h_domains[i]; ++j,
  1198.              domain_x+=m_domain.domain_hstep[i]) 
  1199.       for (k=0,domain_y=0; k<m_domain.no_v_domains[i]; ++k,
  1200.              domain_y+=m_domain.domain_vstep[i]) {
  1201.         Classify(domain_x, domain_y, 
  1202.                  m_domain.domain_hsize[i], 
  1203.                  m_domain.domain_vsize[i], 
  1204.                      &first_class, &second_class,
  1205.                      &m_domain.pixel[i][k][j].sym, 
  1206.                      &m_domain.pixel[i][k][j].sum, 
  1207.                      &m_domain.pixel[i][k][j].sum2, 2);
  1208.  
  1209.         // When the domain data is referenced from the list, we need to 
  1210.         // know where the domain is.. so we have to store the position  
  1211.         m_domain.pixel[i][k][j].dom_x = j;
  1212.         m_domain.pixel[i][k][j].dom_y = k;
  1213.         node = (struct classified_domain *)
  1214.                                 malloc(sizeof(struct classified_domain));
  1215.         
  1216.         // put this domain in the classified list structure 
  1217.         node->dd = &m_domain.pixel[i][k][j];
  1218.         node->next = m_dd_domain[first_class][second_class][i];
  1219.         m_dd_domain[first_class][second_class][i] = node;
  1220.       }
  1221.     }
  1222.  
  1223.     // Now we make sure no domain class is actually empty.  
  1224.     for (i=0; i <= m_max_part-m_min_part; ++i) 
  1225.     for (first_class = 0; first_class < 3; ++first_class)
  1226.     for (second_class = 0; second_class < 24; ++second_class) 
  1227.        if (m_dd_domain[first_class][second_class][i] == NULL) {
  1228.            node = (struct classified_domain *)
  1229.                                 malloc(sizeof(struct classified_domain));
  1230.            node->dd = &m_domain.pixel[i][0][0];
  1231.            node->next = NULL;
  1232.            m_dd_domain[first_class][second_class][i] = node;
  1233.        }
  1234.  
  1235.    
  1236.     
  1237.     Message("Done ");
  1238.     return TRUE;
  1239.     
  1240. }
  1241.  
  1242. int CFracComp::Pack(int size, long int value)
  1243. {
  1244.      int i;
  1245.      static int ptr = 1, // how many bits are packed in sum so far 
  1246.                 sum = 0, // packed bits 
  1247.                 num_of_packed_bytes = 0; // total bytes written out 
  1248.  
  1249.     // size == -1 means we are at the end, so write out what is left 
  1250.     if (size == -1 && ptr != 1) {
  1251.         fputc(sum<<(8-ptr), m_FileOut);
  1252.         ++num_of_packed_bytes;
  1253.         return(0);
  1254.     }
  1255.     
  1256.     // size == -2 means we want to know how many bytes we have written 
  1257.     if (size == -2) 
  1258.             return(num_of_packed_bytes);
  1259.     
  1260.     for (i=0; i<size; ++i, ++ptr, value = value>>1, sum = sum<<1) {
  1261.         if (value & 1) sum |= 1;
  1262.  
  1263.          if (ptr == 8) 
  1264.          {
  1265.             fputc(sum,m_FileOut);
  1266.             ++num_of_packed_bytes;
  1267.             sum=0;
  1268.             ptr=0;
  1269.         }
  1270.      }
  1271.     return(0);
  1272. }
  1273.  
  1274. double CFracComp::Compare(int atx, int aty, int xsize, int ysize, int depth, double rsum, double rsum2, int dom_x, int dom_y, 
  1275.                                   int sym_op, int *pialpha, int *pibeta)
  1276.  
  1277. {
  1278.     int i, j, i1, j1, k,
  1279.         domain_x, 
  1280.         domain_y;        // The domain position                
  1281.  
  1282.     double pixel,
  1283.            det,          // determinant of solution 
  1284.            dsum,         // sum of domain values 
  1285.            rdsum = 0,    // sum of range*domain values   
  1286.            dsum2,        // sum of domain^2 values   
  1287.            w2 = 0,       // total number of values tested 
  1288.            rms = 0,      // root means square difference 
  1289.            alpha,        // the scale factor 
  1290.            beta;         // The offset 
  1291.  
  1292.  
  1293.     
  1294.    
  1295.     w2 = xsize * ysize;
  1296.     dsum = m_domain.pixel[depth-m_min_part][dom_y][dom_x].sum;
  1297.     dsum2 = m_domain.pixel[depth-m_min_part][dom_y][dom_x].sum2;
  1298.     domain_x = (dom_x * m_domain.domain_hstep[depth-m_min_part]);
  1299.     domain_y = (dom_y * m_domain.domain_vstep[depth-m_min_part]);
  1300.     k = ((domain_x%2)<<1) + domain_y%2;
  1301.     domain_x >>= 1;
  1302.     domain_y >>= 1;
  1303.  
  1304.     // The main statement in this routine is a switch statement which 
  1305.     // scans through the domain and range to compare them.  
  1306.      switch(sym_op) {
  1307.          case 0: for (i=atx, i1 = domain_x; i<atx+xsize; ++i, ++i1)
  1308.                      for (j=aty, j1 = domain_y; j<aty+ysize; ++j, ++j1) 
  1309.                      {                                 
  1310.                         pixel = m_domimage[k][j1][i1];                                
  1311.                         rdsum += m_image[j][i]*pixel;                                 
  1312.                      }
  1313.                  break;
  1314.  
  1315.          case 1: for (j=aty, i1 = domain_x; j<aty+ysize; ++j, ++i1)
  1316.                      for (i=atx+xsize-1, j1 = domain_y; i>=atx; --i, ++j1) 
  1317.                      {                                 
  1318.                         pixel = m_domimage[k][j1][i1];                                
  1319.                         rdsum += m_image[j][i]*pixel;                                 
  1320.                      }
  1321.                  break;
  1322.  
  1323.          case 2: for (i=atx+xsize-1, i1 = domain_x; i>=atx; --i, ++i1)
  1324.                      for (j=aty+ysize-1, j1 = domain_y; j>=aty; --j, ++j1) 
  1325.                      {                                 
  1326.                         pixel = m_domimage[k][j1][i1];                                
  1327.                         rdsum += m_image[j][i]*pixel;                                 
  1328.                      }
  1329.                  break;
  1330.  
  1331.          case 3: for (j=aty+ysize-1, i1 = domain_x; j>=aty; --j, ++i1)
  1332.                      for (i=atx, j1 = domain_y; i<atx+xsize; ++i, ++j1) 
  1333.                     {                                 
  1334.                         pixel = m_domimage[k][j1][i1];                                
  1335.                         rdsum += m_image[j][i]*pixel;                                 
  1336.                     }
  1337.                  break;
  1338.  
  1339.          case 4: for (j=aty, i1 = domain_x; j<aty+ysize; ++j, ++i1)
  1340.                      for (i=atx, j1 = domain_y; i<atx+xsize; ++i, ++j1) 
  1341.                      {                                 
  1342.                         pixel = m_domimage[k][j1][i1];                                
  1343.                         rdsum += m_image[j][i]*pixel;                                 
  1344.                      }
  1345.                  break;
  1346.  
  1347.          case 5: for (i=atx, i1 = domain_x; i<atx+xsize; ++i, ++i1)
  1348.                      for (j=aty+ysize-1, j1 = domain_y; j>=aty; --j, ++j1) 
  1349.                      {                                 
  1350.                         pixel = m_domimage[k][j1][i1];                                
  1351.                         rdsum += m_image[j][i]*pixel;                                 
  1352.                      }
  1353.                  break;
  1354.  
  1355.          case 6: for (j=aty+ysize-1, i1 = domain_x; j>=aty; --j, ++i1)
  1356.                      for (i=atx+xsize-1, j1 = domain_y; i>=atx; --i, ++j1) 
  1357.                      {                                 
  1358.                         pixel = m_domimage[k][j1][i1];                                
  1359.                         rdsum += m_image[j][i]*pixel;                                 
  1360.                      }
  1361.                  break;
  1362.  
  1363.          case 7: for (i=atx+xsize-1, i1 = domain_x; i>=atx; --i, ++i1)
  1364.                      for (j=aty, j1 = domain_y; j<aty+ysize; ++j, ++j1) 
  1365.                      {                                 
  1366.                         pixel = m_domimage[k][j1][i1];                                
  1367.                         rdsum += m_image[j][i]*pixel;                                 
  1368.                      }
  1369.                  break;
  1370.      }
  1371.     
  1372.      det = (xsize*ysize)*dsum2 - dsum*dsum;
  1373.     
  1374.      if (det == 0.0) 
  1375.          alpha = 0.0; 
  1376.       else 
  1377.          alpha = (w2*rdsum - rsum*dsum)/det;
  1378.         
  1379.      if  (m_only_positive && alpha < 0.0) alpha = 0.0;
  1380.      *pialpha = 0.5 + (alpha + m_max_scale)/(2.0*m_max_scale)*(1<<m_s_bits);
  1381.      if (*pialpha < 0) *pialpha = 0;
  1382.      if (*pialpha >= 1<<m_s_bits) *pialpha = (1<<m_s_bits)-1;
  1383.  
  1384.      // Now recompute alpha back 
  1385.      alpha = (double)*pialpha/(double)(1<<m_s_bits)*(2.0*m_max_scale)-m_max_scale;
  1386.         
  1387.      //compute the offset 
  1388.      beta = (rsum - alpha*dsum)/w2;
  1389.  
  1390.      // Convert beta to an integer 
  1391.      // we use the sign information of alpha to pack efficiently 
  1392.      if (alpha > 0.0)  beta += alpha*GREY_LEVELS;
  1393.      *pibeta = 0.5 + beta/
  1394.              ((1.0+fabs(alpha))*GREY_LEVELS)*((1<<m_o_bits)-1);
  1395.      if (*pibeta< 0) *pibeta = 0;
  1396.      if (*pibeta>= 1<<m_o_bits) *pibeta = (1<<m_o_bits)-1;
  1397.  
  1398.      // Recompute beta from the integer 
  1399.      beta = (double)*pibeta/(double)((1<<m_o_bits)-1)*
  1400.                ((1.0+fabs(alpha))*GREY_LEVELS);
  1401.      if (alpha > 0.0) beta  -= alpha*GREY_LEVELS;
  1402.  
  1403.      // Compute the rms based on the quantized alpha and beta
  1404.      rms= sqrt((rsum2 + alpha*(alpha*dsum2 - 2.0*rdsum + 2.0*beta*dsum) +
  1405.                  beta*(beta*w2 - 2.0*rsum))/w2);
  1406.     
  1407.      return(rms);
  1408.  
  1409. }
  1410.  
  1411. void CFracComp::Quadtree(int atx, int aty, int xsize, int ysize,double tol,int depth,HWND hdcOut)
  1412. {
  1413.     int i,
  1414.         sym_op,                   // A symmetry operation of the square 
  1415.         ialpha,                   // Intgerized scalling factor         
  1416.         ibeta,                    // Intgerized offset                  
  1417.         best_ialpha,              // best ialpha found                  
  1418.         best_ibeta,
  1419.         best_sym_op,
  1420.         best_domain_x,
  1421.         best_domain_y, 
  1422.         first_class, 
  1423.         the_first_class, 
  1424.         first_class_start,        // loop beginning and ending values   
  1425.         first_class_end, 
  1426.         second_class[2],
  1427.         the_second_class, 
  1428.         second_class_start,       // loop beginning and ending values   
  1429.         second_class_end, 
  1430.         range_sym_op[2],          // the operations to bring square to  
  1431.         domain_sym_op;            // cannonical position.               
  1432.  
  1433.     struct classified_domain *node;  // var for domains we scan through 
  1434.  
  1435.     double rms, best_rms,          // rms value and min rms found so far 
  1436.            sum=0, sum2=0;          // sum and sum^2 of range pixels   
  1437.     HDC tempDc = GetDC(hdcOut); 
  1438.     CString msg; 
  1439.     if (depth < m_min_part) {
  1440.          if(hdcOut != NULL)
  1441.          {               
  1442.              msg.Format("Quadtree is at depth %d",depth+1);
  1443.             SetWindowText(hdcOut,msg);
  1444.              MoveToEx(m_MemDc,(atx+xsize/4),aty,NULL);
  1445.             LineTo(m_MemDc,(atx+xsize/4),ysize/2);
  1446.             MoveToEx(m_MemDc,atx,(aty+ysize/4),NULL);
  1447.             LineTo(m_MemDc,xsize/2,(aty+ysize/4));
  1448.             BitBlt(tempDc,0,0,m_hsize,m_vsize,m_MemDc,0,0,SRCCOPY);
  1449.          } 
  1450.          Quadtree(atx,aty, xsize/2, ysize/2, tol,depth+1,hdcOut);
  1451.          if(hdcOut != NULL)
  1452.          {               
  1453.              msg.Format("Quadtree is at depth %d",depth+1);
  1454.             SetWindowText(hdcOut,msg);
  1455.              MoveToEx(m_MemDc,((atx+xsize/2)+xsize/4),aty,NULL);
  1456.             LineTo(m_MemDc,((atx+xsize/2)+xsize/4),ysize/2);
  1457.             MoveToEx(m_MemDc,(atx+xsize/2),(aty+ysize/4),NULL);
  1458.             LineTo(m_MemDc,xsize/2,(aty+ysize/4));
  1459.             BitBlt(tempDc,0,0,m_hsize,m_vsize,m_MemDc,0,0,SRCCOPY);
  1460.          } 
  1461.          Quadtree(atx+xsize/2,aty, xsize/2, ysize/2,tol,depth+1,hdcOut);
  1462.          if(hdcOut != NULL)
  1463.          {               
  1464.              msg.Format("Quadtree is at depth %d",depth+1);
  1465.             SetWindowText(hdcOut,msg);
  1466.              MoveToEx(m_MemDc,(atx+xsize/4),aty+ysize/2,NULL);
  1467.             LineTo(m_MemDc,(atx+xsize/4),ysize/2);
  1468.             MoveToEx(m_MemDc,atx,((aty+ysize/2)+ysize/4),NULL);
  1469.             LineTo(m_MemDc,xsize/2,((aty+ysize/2)+ysize/4));
  1470.             BitBlt(tempDc,0,0,m_hsize,m_vsize,m_MemDc,0,0,SRCCOPY);
  1471.          } 
  1472.          Quadtree(atx,aty+ysize/2, xsize/2, ysize/2,tol,depth+1,hdcOut);
  1473.          if(hdcOut != NULL)
  1474.          {               
  1475.              msg.Format("Quadtree is at depth %d",depth+1);
  1476.             SetWindowText(hdcOut,msg);
  1477.              MoveToEx(m_MemDc,((atx+xsize/2)+xsize/4),aty+ysize/2,NULL);
  1478.             LineTo(m_MemDc,((atx+xsize/2)+xsize/4),ysize/2);
  1479.             MoveToEx(m_MemDc,atx+xsize/2,((aty+ysize/2)+ysize/4),NULL);
  1480.             LineTo(m_MemDc,xsize/2,((aty+ysize/2)+ysize/4));
  1481.             BitBlt(tempDc,0,0,m_hsize,m_vsize,m_MemDc,0,0,SRCCOPY);
  1482.          } 
  1483.          Quadtree(atx+xsize/2,aty+ysize/2, xsize/2, ysize/2,tol,depth+1,hdcOut);
  1484.          return;
  1485.     }
  1486.  
  1487.     // now search for the best domain-range match and write it out
  1488.     best_rms = 10000000000.0;     
  1489.  
  1490.     Classify(atx, aty, xsize,ysize, 
  1491.               &the_first_class, &the_second_class, 
  1492.               &range_sym_op[0], &sum, &sum2, 1);
  1493.  
  1494.  
  1495.     // sort out how much to search based on -f and -F input flags 
  1496.     if (m_fullclass_search) {
  1497.          first_class_start = 0; 
  1498.          first_class_end = 3;
  1499.     } else {
  1500.          first_class_start = the_first_class; 
  1501.          first_class_end = the_first_class+1;
  1502.     } 
  1503.  
  1504.     if (m_subclass_search) {
  1505.          second_class_start = 0; 
  1506.          second_class_end = 24;
  1507.     } else {
  1508.          second_class_start = the_second_class; 
  1509.          second_class_end = the_second_class+1;
  1510.     } 
  1511.  
  1512.     // these for loops vary according to the optimization flags we set 
  1513.     // for subclass_search and fullclass_search==1, we search all the  
  1514.     // domains (except not all rotations).                             
  1515.     for (first_class = first_class_start; 
  1516.          first_class < first_class_end; ++first_class)
  1517.     for (second_class[0] = second_class_start; 
  1518.          second_class[0] < second_class_end; ++second_class[0]) {
  1519.  
  1520.        // We must check each domain twice. Once for positive scaling,  
  1521.        // once for negative scaling. Each has its own class and sym_op 
  1522.        if (!m_only_positive) {
  1523.           second_class[1] = 
  1524.             m_class_transform[(first_class == 2 ? 1 : 0)][second_class[0]];
  1525.           range_sym_op[1] = 
  1526.              m_rot_transform[(the_first_class == 2 ? 1 : 0)][range_sym_op[0]];
  1527.        } 
  1528.  
  1529.        // only_positive is 0 or 1, so we may or not scan                
  1530.        for (i=0; i<(2-m_only_positive); ++i) 
  1531.        for (node = m_dd_domain[first_class][second_class[i]][depth-m_min_part]; 
  1532.            node != NULL; 
  1533.            node = node->next) {
  1534.            domain_sym_op = node->dd->sym;
  1535.            // The following if statement figures out how to transform  
  1536.            // the domain onto the range, given that we know how to get 
  1537.            // each into cannonical position.                           
  1538.            if (((domain_sym_op>3 ? 4: 0) + (range_sym_op[i]>3 ? 4: 0))%8 == 0)
  1539.              sym_op = (4 + domain_sym_op%4 - range_sym_op[i]%4)%4;
  1540.            else 
  1541.              sym_op = (4 + (domain_sym_op%4 + 3*(4-range_sym_op[i]%4))%4)%8;
  1542.             
  1543.            rms = Compare(atx,aty, xsize, ysize,  depth, sum,sum2, 
  1544.                                   node->dd->dom_x, 
  1545.                                   node->dd->dom_y, 
  1546.                                   sym_op, &ialpha,&ibeta); 
  1547.     
  1548.            if (rms < best_rms) {
  1549.                    best_ialpha = ialpha;
  1550.                    best_ibeta = ibeta;
  1551.                    best_rms = rms;
  1552.                    best_sym_op = sym_op;
  1553.                    best_domain_x = node->dd->dom_x;
  1554.                    best_domain_y = node->dd->dom_y;
  1555.            }
  1556.        }
  1557.     }
  1558.         
  1559.     if (best_rms > tol && depth < m_max_part)
  1560.     {
  1561.        // We didn't find a good enough fit so quadtree down 
  1562.        Pack(1,(long)1);  // This bit means we quadtree'd down 
  1563.             if(hdcOut != NULL)
  1564.          {               
  1565.              msg.Format("Quadtree is at depth %d",depth+1);
  1566.             SetWindowText(hdcOut,msg);
  1567.              MoveToEx(m_MemDc,(atx+xsize/4),aty,NULL);
  1568.             LineTo(m_MemDc,(atx+xsize/4),ysize/2);
  1569.             MoveToEx(m_MemDc,atx,(aty+ysize/4),NULL);
  1570.             LineTo(m_MemDc,xsize/2,(aty+ysize/4));
  1571.             BitBlt(tempDc,0,0,m_hsize,m_vsize,m_MemDc,0,0,SRCCOPY);
  1572.          } 
  1573.          Quadtree(atx,aty, xsize/2, ysize/2, tol,depth+1,hdcOut);
  1574.          if(hdcOut != NULL)
  1575.          {               
  1576.              msg.Format("Quadtree is at depth %d",depth+1);
  1577.             SetWindowText(hdcOut,msg);
  1578.              MoveToEx(m_MemDc,((atx+xsize/2)+xsize/4),aty,NULL);
  1579.             LineTo(m_MemDc,((atx+xsize/2)+xsize/4),ysize/2);
  1580.             MoveToEx(m_MemDc,(atx+xsize/2),(aty+ysize/4),NULL);
  1581.             LineTo(m_MemDc,xsize/2,(aty+ysize/4));
  1582.             BitBlt(tempDc,0,0,m_hsize,m_vsize,m_MemDc,0,0,SRCCOPY);
  1583.          } 
  1584.          Quadtree(atx+xsize/2,aty, xsize/2, ysize/2,tol,depth+1,hdcOut);
  1585.          if(hdcOut != NULL)
  1586.          {               
  1587.              msg.Format("Quadtree is at depth %d",depth+1);
  1588.             SetWindowText(hdcOut,msg);
  1589.              msg.Format("Quadtree is at depth %d",depth+1);
  1590.             SetWindowText(hdcOut,msg);
  1591.              MoveToEx(m_MemDc,(atx+xsize/4),aty+ysize/2,NULL);
  1592.             LineTo(m_MemDc,(atx+xsize/4),ysize/2);
  1593.             MoveToEx(m_MemDc,atx,((aty+ysize/2)+ysize/4),NULL);
  1594.             LineTo(m_MemDc,xsize/2,((aty+ysize/2)+ysize/4));
  1595.             SetWindowText(hdcOut,msg);
  1596.             BitBlt(tempDc,0,0,m_hsize,m_vsize,m_MemDc,0,0,SRCCOPY);
  1597.  
  1598.          } 
  1599.          Quadtree(atx,aty+ysize/2, xsize/2, ysize/2,tol,depth+1,hdcOut);
  1600.          if(hdcOut != NULL)
  1601.          {               
  1602.              msg.Format("Quadtree is at depth %d",depth+1);
  1603.             SetWindowText(hdcOut,msg);
  1604.              MoveToEx(m_MemDc,((atx+xsize/2)+xsize/4),aty+ysize/2,NULL);
  1605.             LineTo(m_MemDc,((atx+xsize/2)+xsize/4),ysize/2);
  1606.             MoveToEx(m_MemDc,atx+xsize/2,((aty+ysize/2)+ysize/4),NULL);
  1607.             LineTo(m_MemDc,xsize/2,((aty+ysize/2)+ysize/4));
  1608.             BitBlt(tempDc,0,0,m_hsize,m_vsize,m_MemDc,0,0,SRCCOPY);
  1609.          } 
  1610.          Quadtree(atx+xsize/2,aty+ysize/2, xsize/2, ysize/2,tol,depth+1,hdcOut);
  1611.     } 
  1612.     else 
  1613.     {
  1614.        // The fit was good enough or we just can't get smaller ranges 
  1615.        // So write out the data 
  1616.        if (depth < m_max_part)              // if we are not at the smallest range 
  1617.                Pack(1,(long)0);             // then we must tell the decoder we    
  1618.                                           // stopped quadtreeing                 
  1619.        Pack(m_s_bits, (long)best_ialpha);
  1620.        Pack(m_o_bits, (long)best_ibeta);
  1621.        // When the scaling is zero, there is no need to store the domain
  1622.        if (best_ialpha != m_zero_ialpha) 
  1623.        {
  1624.           Pack(3, (long)best_sym_op);
  1625.           Pack(m_bits_needed[depth-m_min_part], (long)(best_domain_y*
  1626.             m_domain.no_h_domains[depth-m_min_part]+best_domain_x));
  1627.       }
  1628.     }
  1629.  
  1630. }
  1631.  
  1632. void CFracComp::Partition_Image(int atx, int aty,int hsize, int vsize,double tol,HWND hdcOut)
  1633. {
  1634.    int x_exponent,    // the largest power of 2 image size that fits 
  1635.        y_exponent,    // horizontally or vertically the rectangle.   
  1636.        exponent,      // The actual size of image that's encoded.    
  1637.        size, 
  1638.        depth; 
  1639.    
  1640.    x_exponent = (int)floor(log((double)hsize)/log(2.0));
  1641.    y_exponent = (int)floor(log((double)vsize)/log(2.0));
  1642.    
  1643.    // exponent is min of x_ and y_ exponent
  1644.    exponent = (x_exponent > y_exponent ? y_exponent : x_exponent);
  1645.    size = 1<<exponent; 
  1646.    depth = m_max_exponent - exponent;
  1647.    Quadtree(atx,aty,size,size,tol,depth,hdcOut);
  1648.    if (size != hsize) 
  1649.       Partition_Image(atx+size, aty, hsize-size,vsize, tol,hdcOut);
  1650.  
  1651.    if (size != vsize) 
  1652.       Partition_Image(atx, aty+size, size,vsize-size, tol,hdcOut);
  1653. }
  1654.  
  1655. void CFracComp::Average(int x, int y, int xsize, int ysize, double *psum, double *psum2)
  1656. {
  1657.     register int i,j,k;
  1658.     register double pixel;
  1659.     *psum = *psum2 = 0.0;
  1660.     k = ((x%2)<<1) + y%2;
  1661.     x >>= 1; y >>= 1;
  1662.     xsize >>= 1; ysize >>= 1;
  1663.     for (i=x; i<x+xsize; ++i)
  1664.     for (j=y; j<y+ysize; ++j) 
  1665.     {
  1666.         pixel = m_domimage[k][j][i];
  1667.        *psum += pixel;
  1668.        *psum2 += pixel*pixel;
  1669.     }
  1670. }
  1671.  
  1672. void CFracComp::Average1(int x, int y, int xsize, int ysize, double *psum, double *psum2)
  1673. {
  1674.     register int i,j;
  1675.     register double pixel;
  1676.     *psum = *psum2 = 0.0;
  1677.  
  1678.     for (i=x; i<x+xsize; ++i)
  1679.     for (j=y; j<y+ysize; ++j) {
  1680.        pixel = (double)m_image[j][i];
  1681.        *psum += pixel;
  1682.        *psum2 += pixel*pixel;
  1683.     }
  1684. }
  1685.  
  1686.  
  1687.