/* This is a C++ source code fragment for the LCIS function described in the SIGGRAPH'99 paper "LCIS: A Boundary Hierarchy For Detail-Preserving Contrast Reduction" ---------------------------------------------------------------------- by Jack Tumblin and Greg Turk. I offer this code only as an aid to people who may want to implement LCIS themselves. This code was written to test and explore LCIS easily, not for speed or efficiency. You may notice that the 'leakfix' and 'diffusion' passes through the image could be combined to cut code size almost in half, and that the code is full of expensive sqrt() calls and seemingly pointless casts between PIXTYPE and double etc.; these are vestiges of earlier experiments. Though it is not elegant, it is readable and easy to debug. Please use, share and modify this code in any way you wish. A few items in the code may need some explanation: --'Raw2D' and 'Raw3D' are my own floating-point image classes whose pixel datatype is set by '#DEFINE PIXTYPE float'. Both hold resizeable arrays of pixels; Raw2D holds a 2D array indexed by (x,y), and Raw3D holds a 1D array of Raw2D objects, indexed by 'z'. In the code below, Raw3D objects all hold three Raw2D objects. member functions you'll see here include: get(x,y) -- returns the PIXTYPE pixel value at (x,y) put(x,y,val) -- write PIXTYPE value sizer(arg) -- set array size to match 'arg' wipecopy(src) -- copy contents of 'src', resize as needed. ones() -- set every pixel value to 1.0 zeros() -- set every pixel value to 0.0 --LCIS inputs: The input image read from: m_pA --a pointer to a Raw3D image object; it holds 3 Raw2D objects, and we'll use the middle one of three. To read input pixel value at (x,y), do this: val = m_pA->get(x,y,1); Several user-settable parameters are held in a dialog box class; these are: --# of timesteps to run: found at m_pParamDlg->m_diffusCount --Conductance Threshold 'K'; found at m_pParamDlg->m_diffusKval --Timestep length 'T'; found at m_pParamDlg->m_diffusTstep --Leakfix enable/disable found at M_pParamDlg->m_diffusLeakfix --LCIS outputs: the output image is written to: m_pOut --a pointer to a Raw3D object; it holds 3 Raw2D objects; two hold the current value for 'leakfix multipliers' (D_E,D_N in the paper), the other holds the result of the LCIS operation. m_pOut->z[0] holds the EW leakfix multiplier, m_pOut->z[1] holds the LCIS output image, m_pOut->z[2] holds the NS leakfix multipler. --Quick overview: We begin by copying the input image to the output image object; we never touch the input object again. After initializing a few pesky details, the LCIS code below uses just a few image objects: 'src','bound' and 'm_pOut' (the output image object). The 'src' object holds the result of the most recently completed timestep; we use 'src' as the source of all pixel information we need as we compute results of the current timestep and write them to m_pOut. (m_pOut and src also hold the values of the leakfix multipliers for each pixel). Once we finish finding all new pixel values for m_pOut, we copy m_pOut to src and start on the next timestep. For each timestep there are two passes through the entire image; The first pass computes the 'leakfix multiplier' value for each link, and determines whether or not the link has become a 'boundary link'. If it has, we mark the pixels on either end of the link as a 'boundary pixel' by setting its value in the Raw2D image object called 'bound' to the value IS_BOUND_PIX. The second pass computes the flux through all links and applies the flux to the output image m_pOut. However, it does not compute or apply any flux at all through any link with a boundary pixel on either end. This accomplishes two goals: --no flux will ever transform a ridge-like shock into a step-like shock (prohibited by the LCIS PDEs--see dissertation for proof) --no curvature estimate will use pixels that straddle a shock, where curvature is (theoretically) infinite. If I've still failed to make this code clear to you and you're still interested, send me a question by e-mail: jet@cs.northwestern.edu Regards, --Jack Tumblin */ //3456789_123456789_123456789_123456789_123456789_123456789_123456789_123456789_ // BOOL CFcnGroup::fcnLCIS4(CString & label) //------------------------------------------------------------------------------ // Low Curvature Image Simplifier(LCIS) // // K, timestep, leakfix enable, # iterations, // from SetFcnParam dialog class. // //------------------------------------------------------------------------------ // --prevents any curvature or motive force est. from using a set of pixels // that straddles a 'shock', and // --ensures that a 'shock' always forms a ridge discontinuity that never // (due to smoothing within a region) evolves into a step discontinuity, // mimicking the behavior of the PDE form of LCIS. // Heres how: // --DEFINE: any pixel on either end of a link that contains a shock is a // 'Boundary Pixel', (shock==leakfix multiplier is >10.0) // --APPLY: any link CONNECTED to a boundary pixel must have ZERO FLUX to // ensure the boundary pixel will not change. // Note that this scheme ADIABATIC and fully consistent with the // continuous PDE form of LCIS. //------------------------------------------------------------------------------ // (Return 1 if successful; retn 0 if text 'label' does not match our internal // name for this function; this is a double-check that index# used to select // functions matches the desired function) { Raw2D src; // source img for each diffusion timestep. Raw2D bound; // boundary-pixel map. int i,iE1,iE2,iW,j,jN1,jN2,jS,k,imax,jmax,kmax; // image indices PIXTYPE N2val,N1val,Pval,S1val,W1val,E1val,E2val,NEval,NWval,SEval; // neighborhood pixels: north,south,.. // (See LCIS paper, Figure 6) PIXTYPE Pxx2,Pyy2,Exx2,Eyy2,Nxx2,Nyy2,Nxy2,Sxy2,Wxy2; // 2nd partial derivative ests. (squared) double Nmag2,Emag2; // squared 'edginess' estimates m_EW, m_NS double cXXX,cXXY,cXYY,cYYY; // 3rd partial derivative ests. double fN,fE; // Motive force for NS, EW links. double condN,condE; // conductances for N,E links double gainN,gainE,threshN,threshE; // leakfix temp. vars double gcN,gcE; // temp: gain*conductance. double fluxN,fluxE; // flow for this link, this timestep; double lambtime; // actual timestep length double k2inv; // 1.0/K^2 (K==diffusion constant) #define BOUND_THRESH 10.0 // Boundary threshold for leakfix mpy. #define IS_BOUND_PIX 0.0f // flag value use to mark boundary pixels // in Raw2D 'bound' //******MAGIC CONSTANTS****** double gainLimit = 1E+15; // upper bound on leakfix gain boost double tstep = 1.0/32.0; // default timestep length if(0!=label.CompareNoCase(JT_LCIS4)) // if label mismatch, { BEEP; // COMPLAIN return(FALSE); } //------------------------------------------------------initialize; m_pOut->wipecopy(*m_pA); // copy input image to output image. // use z=0,z=1 planes of output image // to store leakfix multiplier values: m_pOut->z[0].ones(); // set leakfix EWgain to 1.0. m_pOut->z[2].ones(); // set leakfix NSgain to 1.0. imax = m_pOut->getXsize(); // get image size, jmax = m_pOut->getYsize(); // calc a few quick constants; k2inv = 1.0/(m_pParamDlg->m_diffusKval * // 1/(Conductance Threshold K)^2 m_pParamDlg->m_diffusKval); lambtime = tstep * m_pParamDlg->m_diffusTstep; // timestep length kmax = m_pParamDlg->m_diffusCount; // number of timesteps to run. bound.sizer(&(m_pA->z[1])); // boundary-pixel image (0,1 @all pix) bound.ones(); // set to 'no boundary pixels'. //---------------------------------------------------Timestep loop; // use 'src' for input data, write output data to'm_pOut'. for(k=0; kz[1]));// make local image copy. if(m_pParamDlg->m_diffusLeakfix==TRUE) //--------------------------leakfix; { // gain-setting; for(j=0;j=jmax) jN1=jmax-1; // (stay in-bounds) if(jN2>=jmax) jN2=jmax-1; if(jS<0) jS=0; for(i=0; i=imax) iE1=imax-1; // (stay in-bounds) if(iE2>=imax) iE2=imax-1; if(iW<0) iW=0; N2val = src.get(i, jN2); // find neighborhood pixels NWval = src.get(iW, jN1); N1val = src.get(i, jN1); NEval = src.get(iE1,jN1); W1val = src.get(iW, j ); Pval = src.get(i, j ); E1val = src.get(iE1,j ); E2val = src.get(iE2,j ); S1val = src.get(i, jS ); SEval = src.get(iE1,jS ); // find squared xx,yy curvatures through pixels P,E and N: Pxx2 = E1val+W1val-(PIXTYPE)2.0* Pval; Pxx2 = Pxx2*Pxx2; Pyy2 = N1val+S1val-(PIXTYPE)2.0* Pval; Pyy2 = Pyy2*Pyy2; Exx2 = E2val+ Pval-(PIXTYPE)2.0*E1val; Exx2 = Exx2*Exx2; Eyy2 = NEval+SEval-(PIXTYPE)2.0*E1val; Eyy2 = Eyy2*Eyy2; Nxx2 = NEval+NWval-(PIXTYPE)2.0*N1val; Nxx2 = Nxx2*Nxx2; Nyy2 = N2val+ Pval-(PIXTYPE)2.0*N1val; Nyy2 = Nyy2*Nyy2; // find squared xy curvatures through pixel midpoints N,S,W: Nxy2 = (NEval-E1val)-(N1val-Pval); Nxy2 = Nxy2*Nxy2; Sxy2 = (E1val-SEval)-(Pval-S1val); Sxy2 = Sxy2*Sxy2; Wxy2 = (N1val-Pval)-(NWval-W1val); Wxy2 = Wxy2*Wxy2; // Use wt'd. sums of these to est. squared link curv. mag; Emag2 = 0.25*(double)(Pxx2+Pyy2+Exx2+Eyy2) + 0.5*(double)(Nxy2+Sxy2); Nmag2 = 0.25*(double)(Pxx2+Pyy2+Nxx2+Nyy2) + 0.5*(double)(Wxy2+Nxy2); gainE = (double)m_pOut->z[0].get(i,j); // current E leakfix gain; gainN = (double)m_pOut->z[2].get(i,j); // current N leakfix gain; threshE = Emag2*k2inv; // test; threshN = Nmag2*k2inv; // test: // is edginess m > cond. thresh K? // if not, shrink leakfix multiplier // 1/5th of the way towards 1.0; if(threshE<1.0) // gainE = gainE - 0.2*(gainE - 1.0); { // or put more succinctly, gainE = 0.8*gainE + 0.2; } else // but if so, BOOST leakfix multiplier { // (see LCIS paper, equations 21,22) gainE = gainE*(1.0+ sqrt(Emag2)); } if(threshN<1.0) { gainN = 0.8*gainN + 0.2; } else { gainN = gainN*(1.0+ sqrt(Nmag2)); } // write min(gain,gainLimit) if(gainE>gainLimit) gainE = gainLimit; if(gainN>gainLimit) gainN = gainLimit; m_pOut->z[0].put(i,j,(PIXTYPE)gainE); m_pOut->z[2].put(i,j,(PIXTYPE)gainN); if(gainE > BOUND_THRESH) // mark all boundary pixels { // (pixels on either end of // any link holding a shock) bound.put(i , j,IS_BOUND_PIX); bound.put(iE1,j,IS_BOUND_PIX); } if(gainN > BOUND_THRESH) { bound.put(i,j ,IS_BOUND_PIX); bound.put(i,jN1,IS_BOUND_PIX); } }//for(i...) }//for(j...) }//if(leakfix... //--------------------------diffuse for(j=0; j=jmax) jN1=jmax-1; // (stay in-bounds) if(jN2>=jmax) jN2=jmax-1; if(jS<0) jS=0; for(i=0; i=imax) iE1=imax-1; if(iE2>=imax) iE2=imax-1; if(iW<0) iW=0; fluxE = 0.0; // (clear flux values) fluxN = 0.0; // If 'P' is a boundary pixel then NEITHER // the EW or NS link have any flux; if(IS_BOUND_PIX != bound.get(i,j)) { // get common neighbor values; N1val = src.get(i, jN1); NEval = src.get(iE1,jN1); W1val = src.get(iW, j ); Pval = src.get(i, j ); E1val = src.get(iE1,j ); S1val = src.get(i, jS ); // find curv. at pixel P // (for link conductance calcs) Pxx2 = E1val+W1val-(PIXTYPE)2.0*Pval; Pxx2 = Pxx2*Pxx2; Pyy2 = N1val+S1val-(PIXTYPE)2.0*Pval; Pyy2 = Pyy2*Pyy2; // find cross curv. between // pixels P,N1,NE,E1; Nxy2 = (NEval-E1val)-(N1val-Pval); Nxy2 = Nxy2*Nxy2; if(IS_BOUND_PIX !=bound.get(iE1,j)) { // if EW link doesn't connect // to a boundary pixel, then // find its flux; E2val = src.get(iE2,j ); // get its pixel values SEval = src.get(iE1,jS ); // compute change-of-curvature // for motive force fE; cXXX = (E2val - W1val) + 3.0*(Pval - E1val); cXYY = (NEval - N1val) + (SEval - S1val) + 2.0*(Pval - E1val); fE = cXXX + cXYY; // For EW link conductance calcs, find curv. mag.^2 // estimate Emag2: Exx2 = E2val+Pval-(PIXTYPE)2.0*E1val; Exx2 = Exx2*Exx2; Eyy2 = NEval+SEval-(PIXTYPE)2.0*E1val; Eyy2 = Eyy2*Eyy2; // cross-curvature below link; Sxy2 = (E1val-SEval)-(Pval-S1val); Sxy2 = Sxy2*Sxy2; // Make weighted sum of squares to find a rotationally // invariant curvatures mag^2 estimate for EW link; Emag2 = 0.25*(double)(Pxx2+Pyy2+Exx2+Eyy2) + 0.5*(double)(Nxy2+Sxy2); if(m_pParamDlg->m_diffusLeakfix==TRUE) { // leakfix: apply gain to curvature^2 ests. gainE = (double)m_pOut->z[0].get(i,j); gcE = Emag2 * gainE; condE = 1.0 / (1.0 + (k2inv*gcE)); } else { condE = 1.0 / (1.0 + (k2inv*Emag2)); } // flux=lambtime*conductance*motive force; fluxE = lambtime * condE * fE; } if(IS_BOUND_PIX !=bound.get(i,jN1)) { // if NS link doesn't connect // to a boundary pixel, then // find its flux; N2val = src.get(i, jN2); // get its pixel values NWval = src.get(iW, jN1); // compute change-of-curvature // for motive force fN; cYYY = (N2val - S1val) + 3.0*(Pval - N1val); cXXY = (NEval - E1val) + (NWval - W1val) + 2.0*(Pval - N1val); fN = cXXY + cYYY; // For NS link conductance calcs, find curv. mag.^2 Nxx2 = NEval+NWval-(PIXTYPE)2.0*N1val; Nxx2 = Nxx2*Nxx2; Nyy2 = N2val+Pval-(PIXTYPE)2.0*N1val; Nyy2 = Nyy2*Nyy2; // cross-curvature left of link Wxy2 = (N1val-Pval)-(NWval-W1val); Wxy2 = Wxy2*Wxy2; // Make weighted sum of squares to find a rotationally // invariant curvatures mag^2 estimate for NS link; Nmag2 = 0.25*(double)(Pxx2+Pyy2+Nxx2+Nyy2) +0.5*(double)(Wxy2+Nxy2); if(m_pParamDlg->m_diffusLeakfix==TRUE) { // leakfix: apply gain to curvature^2 ests. gainN = (double)m_pOut->z[2].get(i,j); gcN = Nmag2 * gainN; condN = 1.0 / (1.0 + (k2inv*gcN)); } else { condN = 1.0 / (1.0 + (k2inv * Nmag2)); } // flux=lambtime*conductance*motive force; fluxN = lambtime * condN * fN; } // Now apply flux to output image; m_pOut->addTo((PIXTYPE)(-fluxN-fluxE),i,j,1);// center pix m_pOut->addTo((PIXTYPE)fluxN,i, jN1,1); // North m_pOut->addTo((PIXTYPE)fluxE,iE1,j ,1); // East }//if(IS_BOUND_PIX... }//for(i... }// for(j... }//for(k...) #undef BOUND_THRESH // locals only! #undef IS_BOUND_PIX return(TRUE); }