//static
  void GradientsBase::binarizeImage(imageType_t const & currentImage_p, realType_t dBlackPoint_p,
					   weightImageType_t & rMaskImage_p)
  {
#if 0
    imageType_t image(currentImage_p);
    image.ResetSelections();
    image.Binarize(dBlackPoint_p);
    // first channel
    image.ResetSelections();
    image.SelectChannel(0);
    rMaskImage_p.ResetSelections();
    rMaskImage_p.Assign(image);
    // compress into one channel
    for(int i=1;i<image.NumberOfChannels();++i){
      image.SelectChannel(i);
      rMaskImage_p.Max(image);
    }
#else
    int const nCols=currentImage_p.Width();
    int const nRows=currentImage_p.Height();
    int const nChannels=currentImage_p.NumberOfChannels();

    rMaskImage_p.AllocateData(nCols,nRows);
    rMaskImage_p.Black();
    for(int channel=0;channel<nChannels;++channel){
      for(int row=0;row<nRows;++row){
	for(int col=0;col<nCols;++col){
	  if(currentImage_p.Pixel(col,row,channel)>dBlackPoint_p) {
	    rMaskImage_p.Pixel(col,row)=1.0;
	  }
	}
      }
    }
#endif
  }

//static
void
GradientsHdrCompression::minMaxValImage(imageType_t const &rImage_p, realType_t &rdMin_p,realType_t &rdMax_p)
{
  AssertImage(rImage_p);
  rImage_p.ResetSelections();
  rImage_p.GetExtremePixelValues(rdMin_p, rdMax_p);
}

  //static
  void GradientsMergeMosaic::averageImage(pcl_enum eType_p,
					  weightImageType_t const &rCountImage_p,
					  imageType_t &rImage_p)
  {
    int const nCols=rImage_p.Width();
    int const nRows=rImage_p.Height();
    int const nChannels=rImage_p.NumberOfChannels();
    switch(eType_p) {
      // FIXME sure there must be a faster way to do this via masked ops?!?
    case GradientsMergeMosaicType::Overlay:
      // do nothing
      break;
    case GradientsMergeMosaicType::Average:
      // FIXME use builtin PCL for this
      for(int channel=0;channel<nChannels;++channel){
	for(int row=0;row<nRows;++row){
	  for(int col=0;col<nCols;++col){
       	    if(rCountImage_p.Pixel(col,row)>0.0){
	      rImage_p.Pixel(col,row,channel)/=static_cast<realType_t>(rCountImage_p.Pixel(col,row));
	    }
	  }
	}
      }
      break;
    default:
      Assert(false);
    }
  }

void
GradientsBase::AssertColImage(imageType_t const &rImage_p)
{
  Assert(rImage_p.Height()>2);
  Assert(rImage_p.Width()>2);
  Assert(rImage_p.NumberOfChannels()>=1);
}

  //static
  void GradientsHdrComposition::clipGradients(imageType_t &rImage_p)
  {
    for ( realType_t* f = rImage_p.PixelData(), * f1 = f + rImage_p.NumberOfPixels(); f != f1; ++f )
      {
	if ( std::abs(*f) > -MINLOGVALUE) {
	  *f = 0.0;
	}
      }
  }

  //static
  void GradientsHdrComposition::logImage(imageType_t const & rCurrentImage_p,
					 imageType_t & rResultImage_p)
  {
    int const nCols=rCurrentImage_p.Width();
    int const nRows=rCurrentImage_p.Height();
    int const nChannels=rCurrentImage_p.NumberOfChannels();
    imageType_t::color_space  colorSpace=rCurrentImage_p.ColorSpace();

    if (&rCurrentImage_p != &rResultImage_p){
      rResultImage_p.AllocateData(nCols,nRows, nChannels, colorSpace);
    }
    for(int channel=0;channel<nChannels;++channel){
      for(int row=0;row<nRows;++row){
	for(int col=0;col<nCols;++col){
	  realType_t val=rCurrentImage_p.Pixel(col,row,channel);
	  if(val<=0.0) {
	    rResultImage_p.Pixel(col,row,channel)= -1000; //this is just some value to identify extremes
	  } else {
	    val=std::log(val);
	    rResultImage_p.Pixel(col,row,channel)=val;
	  }
	}
      }
    }
  }

void
GradientsHdrCompression::hdrCompression(realType_t maxGradient_p, realType_t minGradient_p, realType_t expGradient_p,
			   bool bRescale01_p, imageType_t &rResultImage_p) const
{
  TimeMessage startHdrCompression("Gradient Domain HDR Compression");
  imageType_t li;
  {
    TimeMessage startClipValue("Determine clip values for gradients");
    imageType_t dx(m_imageDx);
    imageType_t dy(m_imageDy);
    // now clip values
    realType_t minGradientX,maxGradientX;
    minMaxValImage(dx,minGradientX,maxGradientX);
    realType_t minGradientY,maxGradientY;
    minMaxValImage(dy,minGradientY,maxGradientY);
    realType_t minValue=Min(minGradientX, minGradientY);
    realType_t maxValue=Max(maxGradientX, maxGradientY);
    double rangeValue=Max(-minValue, maxValue);
    realType_t const clipRange=rangeValue*maxGradient_p;
    realType_t const zeroRange=rangeValue*minGradient_p;
    startClipValue.Stop();

    TimeMessage startClip("Clipping Gradients");
    clipImage(dx,-clipRange,clipRange);
    zeroRangeImage(dx,-zeroRange,zeroRange);
    clipImage(dy,-clipRange,clipRange);
    zeroRangeImage(dy,-zeroRange,zeroRange);
    startClip.Stop();

    if(expGradient_p!=1.0){
      TimeMessage start("Pow() transformation of gradients");
      absPowImage(dx,expGradient_p);
      absPowImage(dy,expGradient_p);
    }

    TimeMessage startLaplace("Computing 2nd derivative");
    createLaplaceVonNeumannImage(dx,dy,li);
  }
  TimeMessage startSolve("Solving image");
  solveImage(li,rResultImage_p);
  startSolve.Stop();

  TimeMessage startRescale("Rescaling image");
  if(bRescale01_p){
    rResultImage_p.Rescale(0.0,1.0);
  } else {
    rResultImage_p.Rescale(m_dMinVal,m_dMaxVal);
  }
}

//static
void
GradientsBase::createDyImage(imageType_t const &rImage_p, imageType_t &rResultImage_p)
{
  AssertColImage(rImage_p);
  int const nRows=rImage_p.Height();
  int const nCols=rImage_p.Width();
  int const nChannels=rImage_p.NumberOfChannels();
  imageType_t::color_space colorSpace=rImage_p.ColorSpace();

  rResultImage_p.AllocateData(nCols,nRows-1,nChannels,colorSpace);
  for(int chan=0;chan<nChannels;++chan){
    for (int row=0;row<nRows-1;++row){
      for(int col=0;col<nCols;++col){
	rResultImage_p.Pixel(col,row,chan)=rImage_p.Pixel(col,row+1,chan)-rImage_p.Pixel(col,row,chan);
      }
    }
  }
}

  //static
  void GradientsHdrComposition::expImage(imageType_t const & rCurrentImage_p,
					 imageType_t & rResultImage_p)
  {
    int const nCols=rCurrentImage_p.Width();
    int const nRows=rCurrentImage_p.Height();
    int const nChannels=rCurrentImage_p.NumberOfChannels();
    imageType_t::color_space  colorSpace=rCurrentImage_p.ColorSpace();

    if (&rCurrentImage_p != &rResultImage_p){
      rResultImage_p.AllocateData(nCols,nRows, nChannels, colorSpace);
    }
    for(int channel=0;channel<nChannels;++channel){
      for(int row=0;row<nRows;++row){
	for(int col=0;col<nCols;++col){
	  realType_t val=rCurrentImage_p.Pixel(col,row,channel);
	  rResultImage_p.Pixel(col,row,channel)=std::exp(val);
	}
      }
    }
  }

  //static
  void GradientsHdrComposition::addToImage(imageType_t const & rImage_p,
					   int imageNum_p,
					   weightImageType_t const & rMask_p,
					   imageType_t &rSumImageDx_p,imageType_t &rSumImageDy_p,
					   numImageType_t &rDxImage_p,
					   numImageType_t &rDyImage_p)
  {
    int const nCols=rImage_p.Width();
    int const nRows=rImage_p.Height();
    int const nChannels=rImage_p.NumberOfChannels();

    Assert(nCols==rSumImageDx_p.Width()+1);
    Assert(nRows==rSumImageDx_p.Height());
    Assert(nChannels==rSumImageDx_p.NumberOfChannels());

    Assert(nCols==rSumImageDy_p.Width());
    Assert(nRows==rSumImageDy_p.Height()+1);
    Assert(nChannels==rSumImageDy_p.NumberOfChannels());

    Assert(nCols==rDxImage_p.Width());
    Assert(nRows==rDxImage_p.Height());
    Assert(nChannels==rDxImage_p.NumberOfChannels());

    Assert(nCols==rDyImage_p.Width());
    Assert(nRows==rDyImage_p.Height());
    Assert(nChannels==rDyImage_p.NumberOfChannels());

    imageType_t dxImage, dyImage;
    const double zeroLimit=0.0; /// limit for weight that is considered zero
    // handle Dx
    TimeMessage startDx("Creating Dx");
    createDxImage(rImage_p,dxImage);
    startDx.Stop();
    TimeMessage startAddDx("Adding Dx to gradients");
    // transfer useful dx pixels
    for(int row=0;row<nRows;++row){
      for(int col=0;col<nCols-1;++col){
	if(rMask_p.Pixel(col,row)>zeroLimit && rMask_p.Pixel(col+1,row)>zeroLimit) {
	  // we are inside of image, and have useful gradient there
	  for(int channel=0;channel<nChannels;++channel){
	    realType_t currentVal=dxImage.Pixel(col,row,channel);
	    realType_t sumVal=rSumImageDx_p.Pixel(col,row,channel);
	    if(std::abs(currentVal)>std::abs(sumVal)){
	      rSumImageDx_p.Pixel(col,row,channel)=currentVal;
	      rDxImage_p.Pixel(col,row,channel)=imageNum_p;
	    } //if abs
	  } // for chan
	  // FIXME may need to add border handling
	}  //if inside
      }  //for col
    }  //for row
    dxImage.AllocateData(0,0); //save some memory
    startAddDx.Stop();

    // handle Dy just as dx
    TimeMessage startDy("Creating Dy");
    createDyImage(rImage_p,dyImage);
    startDy.Stop();
    TimeMessage startAddDy("Adding Dxy to gradients");
    // transfer useful dy pixels
    for(int row=0;row<nRows-1;++row){
      for(int col=0;col<nCols;++col){
	if(rMask_p.Pixel(col,row)>zeroLimit && rMask_p.Pixel(col,row+1)>zeroLimit) {
	  for(int channel=0;channel<nChannels;++channel){
	    // we are inside of image, and have useful gradient there
	    realType_t currentVal=dyImage.Pixel(col,row,channel);
	    realType_t sumVal=rSumImageDy_p.Pixel(col,row,channel);
	    if(std::abs(currentVal)>std::abs(sumVal)){
	      rSumImageDy_p.Pixel(col,row,channel)=currentVal;
	      rDyImage_p.Pixel(col,row,channel)=imageNum_p;
	    } //if abs()
	  } // for chan
	  // FIXME may need to add border handling
	}  // if inside
      }  //for col
    }  //for row
    startAddDy.Stop();
  }

  //static
  void GradientsMergeMosaic::addToImage(imageType_t const & rImage_p,pcl_enum eType_p,
					weightImageType_t const & rMask_p,
					imageType_t &rSumImageDx_p,imageType_t &rSumImageDy_p,
					weightImageType_t &rCountImageDx_p, weightImageType_t &rCountImageDy_p)
  {
    int const nCols=rImage_p.Width();
    int const nRows=rImage_p.Height();
    int const nChannels=rImage_p.NumberOfChannels();
    imageType_t dxImage, dyImage;

    const double zeroLimit=0.0; /// limit for weight that is considered zero


    TimeMessage startDx("Creating Dx");
    createDxImage(rImage_p,dxImage);
    startDx.Stop();
    TimeMessage startAddDx("Adding Dx");
    // transfer useful dx pixels
    // FIXME this is a relatively slow loop. Think about making it faster
    for(int row=0;row<nRows;++row){
      for(int col=0;col<nCols-1;++col){
	if(rMask_p.Pixel(col,row)>zeroLimit && rMask_p.Pixel(col+1,row)>zeroLimit) {
	  // we are inside of image
	  realType_t weight=(rMask_p.Pixel(col,row)+rMask_p.Pixel(col+1,row))/2.0;
	  if(eType_p==GradientsMergeMosaicType::Average){
	    if(rCountImageDx_p.Pixel(col,row)<=0.0){
	      // first foreground pixel on this location
	      for(int channel=0;channel<nChannels;++channel){
		rSumImageDx_p.Pixel(col,row,channel)=dxImage.Pixel(col,row,channel)*weight;
	      }
	    } else {
	      // there have been other pixels. Create average
	      for(int channel=0;channel<nChannels;++channel){
		rSumImageDx_p.Pixel(col,row,channel)+=dxImage.Pixel(col,row,channel)*weight;
	      }
	    }
	    rCountImageDx_p.Pixel(col,row)+=weight;
	  } else {
	    // type overlay, last gradient wins
	    if(rCountImageDx_p.Pixel(col,row)<=0.0){
		// first foreground pixel on this location
	      for(int channel=0;channel<nChannels;++channel){
		rSumImageDx_p.Pixel(col,row,channel)=dxImage.Pixel(col,row,channel);
	      }
	      rCountImageDx_p.Pixel(col,row)=1.0; //mark as used
	    } else {
	      // there have been other pixels. Blend
	      for(int channel=0;channel<nChannels;++channel){
		rSumImageDx_p.Pixel(col,row,channel)=dxImage.Pixel(col,row,channel)*weight+rSumImageDx_p.Pixel(col,row,channel)*(1.0-weight);
	      }
	    }
	  } //if type
	} else if(rCountImageDx_p.Pixel(col,row)==0.0 &&
		  (rMask_p.Pixel(col,row)<0.0 || rMask_p.Pixel(col+1,row)<0.0)) {
	  // we are at border of image and dont have values there. Just copy in gradient so if nothing else comes in, we have at least the border
	  for(int channel=0;channel<nChannels;++channel){
	    rSumImageDx_p.Pixel(col,row,channel)=dxImage.Pixel(col,row,channel);
	  }
	  // add if first should win. Otherwise last will win
	  //rCountImageDx_p.Pixel(col,row)=-1.0; // mark as background already occupied
	}
      }  //for col
    }  //for row
    dxImage.AllocateData(0,0); //save some memory
    startAddDx.Stop();

    TimeMessage startDy("Creating Dy");
    // transfer useful dy pixels
    createDyImage(rImage_p,dyImage);
    startDy.Stop();

    TimeMessage startAddDy("Adding Dy");
    for(int row=0;row<nRows-1;++row){
      for(int col=0;col<nCols;++col){
	if(rMask_p.Pixel(col,row)>zeroLimit && rMask_p.Pixel(col,row+1)>zeroLimit) {
	  realType_t weight=(rMask_p.Pixel(col,row)+rMask_p.Pixel(col,row+1))/2.0;
	  if(eType_p==GradientsMergeMosaicType::Average){
	    // type average and inside
	    if(rCountImageDy_p.Pixel(col,row)<=0.0){
	      // first foreground pixel on this location
	      for(int channel=0;channel<nChannels;++channel){
		rSumImageDy_p.Pixel(col,row,channel)=dyImage.Pixel(col,row,channel)*weight;
	      }
	    } else {
	      // we already were there. Blend
	      for(int channel=0;channel<nChannels;++channel){
		rSumImageDy_p.Pixel(col,row,channel)+=dyImage.Pixel(col,row,channel)*weight;
	      }
	    }
	    rCountImageDy_p.Pixel(col,row)+=weight;
	  } else {
	    // type overlay and inside, last gradient wins
	    if(rCountImageDy_p.Pixel(col,row)<=0.0){
	      // first foreground pixel on this location
	      for(int channel=0;channel<nChannels;++channel){
		rSumImageDy_p.Pixel(col,row,channel)=dyImage.Pixel(col,row,channel);
	      }
	      rCountImageDy_p.Pixel(col,row)=1.0; //mark as used
	    } else {
//.........这里部分代码省略.........

  //static
  void GradientsMergeMosaic::mergeMosaic(imageListType_t const & rImageList_p,
					 realType_t dBlackPoint_p,
					 pcl_enum eType_p,
					 int32 shrinkCount_p,
					 int32 featherRadius_p,
					 imageType_t &rResultImage_p,
					 sumMaskImageType_t &rSumMaskImage_p)
  {
    Assert(rImageList_p.Length()>=1);
    bool firstImage=true;
    int nCols=0,nRows=0,nChannels=0; /// size and color space of first image
    imageType_t::color_space colorSpace;
    weightImageType_t countImageDx, countImageDy; /// number of pixels that contributed to sumImageDx,Dy in average mode
    imageType_t sumImageDx, sumImageDy; /// combined gradients in x and y direction. Note: these gradients are *between* the pixels
                                        /// of the original image, so size is one less then original image is direction of derivative
    int nImages=0; /// number of images read
    const int enlargeSize=1; // number of pixels added at the border

    TimeMessage startMergeMosaic("Gradient Domain Merge Mosaic");

    TimeMessage startLoadImages("Loading images");
    for(std::size_t i=0;i<rImageList_p.Length();++i){
      imageType_t currentImage;
      int imageIndex=0;
      // allow for multi-image files
      while(loadFile(rImageList_p[i],imageIndex,currentImage)){
	++nImages;
	++imageIndex;
	// expand image dimensions so I have sufficient border for morpological transform and convolution
	TimeMessage startEnlarge("creating border");
	currentImage.CropBy(enlargeSize,enlargeSize,enlargeSize,enlargeSize);
	startEnlarge.Stop();

	if(firstImage){
	  firstImage=false;
	  // determine those parameters that must be shared by all images
	  nCols=currentImage.Width();
	  nRows=currentImage.Height();
	  nChannels=currentImage.NumberOfChannels();
	  colorSpace=currentImage.ColorSpace();

	  //allocate necessary helper images
	  rSumMaskImage_p.AllocateData(nCols,nRows);
	  rSumMaskImage_p.ResetSelections();
	  rSumMaskImage_p.Black();
	  sumImageDx.AllocateData(nCols-1,nRows,nChannels,colorSpace);
	  sumImageDx.ResetSelections();
	  sumImageDx.Black();
	  sumImageDy.AllocateData(nCols,nRows-1,nChannels,colorSpace);
	  sumImageDy.ResetSelections();
	  sumImageDy.Black();

	  countImageDx.AllocateData(nCols-1,nRows);
	  countImageDx.Black();
	  countImageDy.AllocateData(nCols,nRows-1);
	  countImageDy.Black();
	} else {
	  // FIXME I wonder if I should check color space etc as well...
	  // check if properties of this image are identical to those of the first image
	  if(nCols!=currentImage.Width()) {
	    throw Error("Current image width differs from first image width.");
	  } else if(nRows!=currentImage.Height()) {
	    throw Error("Current image height differs from first image height.");
	  }	else if(nChannels!=currentImage.NumberOfChannels()) {
	    throw Error("Current image number of channels differs from first image number of channels.");
	  }
	}
	TimeMessage startProcessImage("Processing Image"+String(nImages));
	mergeMosaicProcessImage(currentImage,dBlackPoint_p,eType_p,shrinkCount_p,featherRadius_p,sumImageDx, sumImageDy ,rSumMaskImage_p,countImageDx, countImageDy);
      }
    }
    startLoadImages.Stop();
    if (eType_p==GradientsMergeMosaicType::Average) {
      TimeMessage startAverage("Averaging images");
      averageImage(eType_p,countImageDx,sumImageDx);
      averageImage(eType_p,countImageDy,sumImageDy);
      // we do not need count images any longer
      countImageDx.AllocateData(0,0);
      countImageDy.AllocateData(0,0);
    }
    // at this point:
    // sumImageDx: Average or overlay of gradients of images read in x direction
    // sumImageDy: Average or overlay of gradients of images read in y direction
    // rSumMaskImage_p: mask with different values for the different sources of images. 0 is background.
    //              We use this later for information of the user, but it is not needed in the following process
    TimeMessage startMerge("Merging Images");
    imageType_t laplaceImage;
    TimeMessage startLaplace("Creating Laplace image");
    createLaplaceVonNeumannImage(sumImageDx,sumImageDy,laplaceImage);
    startLaplace.Stop();
    TimeMessage startSolve("Solving Laplace");
    solveImage(laplaceImage,rResultImage_p);
    startSolve.Stop();
    startMerge.Stop();
    rResultImage_p.ResetSelections();
#if 0
    // for debugging laplaceImage
    //    rResultImage_p.Assign(laplaceImage);
    rResultImage_p.Assign(sumImageDx);
//.........这里部分代码省略.........

  //static
  void GradientsMergeMosaic::mergeMosaicProcessImage(imageType_t const & rImage_p, realType_t dBlackPoint_p,
						     pcl_enum eType_p,
						     int32 shrinkCount_p,
						     int32 featherRadius_p,
						     imageType_t &rSumImageDx_p,
						     imageType_t &rSumImageDy_p,
						     sumMaskImageType_t &rSumMaskImage_p,
						     weightImageType_t &rCountImageDx_p,
						     weightImageType_t &rCountImageDy_p)
  {
#ifdef DEBUG
    int const nCols=rImage_p.Width();
    int const nRows=rImage_p.Height();
    int const nChannels=rImage_p.NumberOfChannels();
#endif
    Assert(nCols==rSumImageDx_p.Width()+1);
    Assert(nRows==rSumImageDx_p.Height());
    Assert(nChannels==rSumImageDx_p.NumberOfChannels());

    Assert(nCols==rSumImageDy_p.Width());
    Assert(nRows==rSumImageDy_p.Height()+1);
    Assert(nChannels==rSumImageDy_p.NumberOfChannels());

    Assert(nCols==rSumMaskImage_p.Width());
    Assert(nRows==rSumMaskImage_p.Height());
    Assert(1==rSumMaskImage_p.NumberOfChannels());

    Assert(eType_p!=GradientsMergeMosaicType::Average || nCols==rCountImageDx_p.Width()+1);
    Assert(eType_p!=GradientsMergeMosaicType::Average || nRows==rCountImageDx_p.Height());
    Assert(eType_p!=GradientsMergeMosaicType::Average || 1==rCountImageDx_p.NumberOfChannels());
    Assert(eType_p!=GradientsMergeMosaicType::Average || nCols==rCountImageDy_p.Width());
    Assert(eType_p!=GradientsMergeMosaicType::Average || nRows==rCountImageDy_p.Height()+1);
    Assert(eType_p!=GradientsMergeMosaicType::Average || 1==rCountImageDy_p.NumberOfChannels());

    Assert(shrinkCount_p>=0);
    Assert(featherRadius_p>=0);
    Assert(dBlackPoint_p>=0.0);
    weightImageType_t maskImage;

    TimeMessage startAddImage("Adding image to data");

    TimeMessage startBinarize("Binarize Image");
    binarizeImage(rImage_p,dBlackPoint_p,maskImage);
    // save this for border computation later
    weightImageType_t fullMask(maskImage);
    startBinarize.Stop();

    // we are doing this because image after StarAlign usually contain aliased pixels.
    // These must not to be used during merge.
    TimeMessage startShrink("Shrinking mask");
    erodeMask(maskImage,shrinkCount_p);
    startShrink.Stop();
    TimeMessage startFeather("Feathering mask");
    featherMask(maskImage,featherRadius_p);
    startFeather.Stop();

    TimeMessage startBorder("Computing border");
    addBorder(fullMask,maskImage);
    fullMask.AllocateData(0,0); // save memory
    startBorder.Stop();

    TimeMessage startSumMask("Creating combined mask");
    addToMask(maskImage,eType_p,rSumMaskImage_p);
    startSumMask.Stop();
    TimeMessage startAddGradients("Adding gradients data");
    addToImage(rImage_p,eType_p,maskImage,rSumImageDx_p,rSumImageDy_p,rCountImageDx_p, rCountImageDy_p);
    startAddGradients.Stop();
  }