22 """Support utilities for Measuring sources"""
24 __all__ = [
"splitId",
"showSourceSet",
"showPsfSpatialCells",
25 "showPsfCandidates",
"makeSubplots",
"plotPsfSpatialModel",
26 "showPsf",
"showPsfMosaic",
"showPsfResiduals",
"saveSpatialCellSet"]
43 from .
import subtractPsf, fitKernelParamsToImage
47 afwDisplay.setDefaultMaskTransparency(75)
52 objId = int((oid & 0xffff) - 1)
55 return dict(objId=objId)
60 def showSourceSet(sSet, xy0=(0, 0), display=
None, ctype=afwDisplay.GREEN, symb=
"+", size=2):
61 """Draw the (XAstrom, YAstrom) positions of a set of Sources. Image has the given XY0"""
64 display = afwDisplay.Display()
65 with display.Buffering():
67 xc, yc = s.getXAstrom() - xy0[0], s.getYAstrom() - xy0[1]
70 display.dot(str(
splitId(s.getId(),
True)[
"objId"]), xc, yc, ctype=ctype, size=size)
72 display.dot(symb, xc, yc, ctype=ctype, size=size)
80 symb=None, ctype=None, ctypeUnused=None, ctypeBad=None, size=2, display=None):
81 """Show the SpatialCells.
83 If symb is something that afwDisplay.Display.dot() understands (e.g. "o"),
84 the top nMaxPerCell candidates will be indicated with that symbol, using
89 display = afwDisplay.Display()
90 with display.Buffering():
91 origin = [-exposure.getMaskedImage().getX0(), -exposure.getMaskedImage().getY0()]
92 for cell
in psfCellSet.getCellList():
93 displayUtils.drawBBox(cell.getBBox(), origin=origin, display=display)
99 goodies = ctypeBad
is None
100 for cand
in cell.begin(goodies):
104 xc, yc = cand.getXCenter() + origin[0], cand.getYCenter() + origin[1]
110 color = ctypeBad
if cand.isBad()
else ctype
118 display.dot(symb, xc, yc, ctype=ct, size=size)
120 source = cand.getSource()
123 rchi2 = cand.getChi2()
126 display.dot(
"%d %.1f" % (
splitId(source.getId(),
True)[
"objId"], rchi2),
127 xc - size, yc - size - 4, ctype=color, size=2)
130 display.dot(
"%.2f %.2f %.2f" % (source.getIxx(), source.getIxy(), source.getIyy()),
131 xc-size, yc + size + 4, ctype=color, size=size)
135 def showPsfCandidates(exposure, psfCellSet, psf=None, display=None, normalize=True, showBadCandidates=True,
136 fitBasisComponents=False, variance=None, chi=None):
137 """Display the PSF candidates.
139 If psf is provided include PSF model and residuals; if normalize is true normalize the PSFs
142 If chi is True, generate a plot of residuals/sqrt(variance), i.e. chi
144 If fitBasisComponents is true, also find the best linear combination of the PSF's components
148 display = afwDisplay.Display()
151 if variance
is not None:
156 mos = displayUtils.Mosaic()
158 candidateCenters = []
159 candidateCentersBad = []
162 for cell
in psfCellSet.getCellList():
163 for cand
in cell.begin(
False):
164 rchi2 = cand.getChi2()
168 if not showBadCandidates
and cand.isBad():
172 im_resid = displayUtils.Mosaic(gutter=0, background=-5, mode=
"x")
175 im = cand.getMaskedImage()
176 xc, yc = cand.getXCenter(), cand.getYCenter()
178 margin = 0
if True else 5
179 w, h = im.getDimensions()
183 bim = im.Factory(w + 2*margin, h + 2*margin)
187 bim.getVariance().
set(stdev**2)
194 im = im.Factory(im,
True)
195 im.setXY0(cand.getMaskedImage().getXY0())
200 im_resid.append(im.Factory(im,
True))
204 psfIm = mi.getImage()
205 config = measBase.SingleFrameMeasurementTask.ConfigClass()
206 config.slots.centroid =
"base_SdssCentroid"
208 schema = afwTable.SourceTable.makeMinimalSchema()
209 measureSources = measBase.SingleFrameMeasurementTask(schema, config=config)
213 miBig = mi.Factory(im.getWidth() + 2*extra, im.getHeight() + 2*extra)
214 miBig[extra:-extra, extra:-extra, afwImage.LOCAL] = mi
215 miBig.setXY0(mi.getX0() - extra, mi.getY0() - extra)
225 footprintSet.makeSources(catalog)
227 if len(catalog) == 0:
228 raise RuntimeError(
"Failed to detect any objects")
230 measureSources.run(catalog, exp)
231 if len(catalog) == 1:
235 for i, s
in enumerate(catalog):
236 d = numpy.hypot(xc - s.getX(), yc - s.getY())
237 if i == 0
or d < dmin:
239 xc, yc = source.getCentroid()
249 resid = resid.getImage()
250 var = im.getVariance()
251 var = var.Factory(var,
True)
252 numpy.sqrt(var.getArray(), var.getArray())
255 im_resid.append(resid)
258 if fitBasisComponents:
259 im = cand.getMaskedImage()
261 im = im.Factory(im,
True)
262 im.setXY0(cand.getMaskedImage().getXY0())
265 noSpatialKernel = psf.getKernel()
267 noSpatialKernel =
None
276 outImage = afwImage.ImageD(outputKernel.getDimensions())
277 outputKernel.computeImage(outImage,
False)
279 im -= outImage.convertF()
283 bim = im.Factory(w + 2*margin, h + 2*margin)
287 bim.assign(resid, bbox)
291 resid = resid.getImage()
294 im_resid.append(resid)
296 im = im_resid.makeMosaic()
298 im = cand.getMaskedImage()
303 objId =
splitId(cand.getSource().getId(),
True)[
"objId"]
305 lab =
"%d chi^2 %.1f" % (objId, rchi2)
306 ctype = afwDisplay.RED
if cand.isBad()
else afwDisplay.GREEN
308 lab =
"%d flux %8.3g" % (objId, cand.getSource().getPsfInstFlux())
309 ctype = afwDisplay.GREEN
311 mos.append(im, lab, ctype)
313 if False and numpy.isnan(rchi2):
314 display.mtv(cand.getMaskedImage().getImage(), title=
"showPsfCandidates: candidate")
315 print(
"amp", cand.getAmplitude())
317 im = cand.getMaskedImage()
318 center = (candidateIndex, xc - im.getX0(), yc - im.getY0())
321 candidateCentersBad.append(center)
323 candidateCenters.append(center)
326 title =
"chi(Psf fit)"
328 title =
"Stars & residuals"
329 mosaicImage = mos.makeMosaic(display=display, title=title)
331 with display.Buffering():
332 for centers, color
in ((candidateCenters, afwDisplay.GREEN), (candidateCentersBad, afwDisplay.RED)):
334 bbox = mos.getBBox(cen[0])
335 display.dot(
"+", cen[1] + bbox.getMinX(), cen[2] + bbox.getMinY(), ctype=color)
340 def makeSubplots(fig, nx=2, ny=2, Nx=1, Ny=1, plottingArea=(0.1, 0.1, 0.85, 0.80),
341 pxgutter=0.05, pygutter=0.05, xgutter=0.04, ygutter=0.04,
342 headroom=0.0, panelBorderWeight=0, panelColor=
'black'):
343 """Return a generator of a set of subplots, a set of Nx*Ny panels of nx*ny plots. Each panel is fully
344 filled by row (starting in the bottom left) before the next panel is started. If panelBorderWidth is
345 greater than zero a border is drawn around each panel, adjusted to enclose the axis labels.
348 subplots = makeSubplots(fig, 2, 2, Nx=1, Ny=1, panelColor='k')
349 ax = subplots.next(); ax.text(0.3, 0.5, '[0, 0] (0,0)')
350 ax = subplots.next(); ax.text(0.3, 0.5, '[0, 0] (1,0)')
351 ax = subplots.next(); ax.text(0.3, 0.5, '[0, 0] (0,1)')
352 ax = subplots.next(); ax.text(0.3, 0.5, '[0, 0] (1,1)')
357 fig : `matplotlib.pyplot.figure`
358 The matplotlib figure to draw
360 The number of plots in each row of each panel
362 The number of plots in each column of each panel
364 The number of panels in each row of the figure
366 The number of panels in each column of the figure
367 plottingArea : `tuple`
368 (x0, y0, x1, y1) for the part of the figure containing all the panels
370 Spacing between columns of panels in units of (x1 - x0)
372 Spacing between rows of panels in units of (y1 - y0)
374 Spacing between columns of plots within a panel in units of (x1 - x0)
376 Spacing between rows of plots within a panel in units of (y1 - y0)
378 Extra spacing above each plot for e.g. a title
379 panelBorderWeight : `int`
380 Width of border drawn around panels
382 Colour of border around panels
387 import matplotlib.pyplot
as plt
388 except ImportError
as e:
389 log.warning(
"Unable to import matplotlib: %s", e)
395 except AttributeError:
396 fig.__show = fig.show
403 fig.show = types.MethodType(myShow, fig)
412 Callback to draw the panel borders when the plots are drawn to the canvas
414 if panelBorderWeight <= 0:
417 for p
in axes.keys():
420 bboxes.append(ax.bbox.union([label.get_window_extent()
for label
in
421 ax.get_xticklabels() + ax.get_yticklabels()]))
428 bbox = ax.bbox.union(bboxes)
430 xy0, xy1 = ax.transData.inverted().
transform(bbox)
433 w, h = x1 - x0, y1 - y0
442 rec = ax.add_patch(plt.Rectangle((x0, y0), w, h, fill=
False,
443 lw=panelBorderWeight, edgecolor=panelColor))
444 rec.set_clip_on(
False)
448 fig.canvas.mpl_connect(
'draw_event', on_draw)
452 x0, y0 = plottingArea[0:2]
453 W, H = plottingArea[2:4]
454 w = (W - (Nx - 1)*pxgutter - (nx*Nx - 1)*xgutter)/float(nx*Nx)
455 h = (H - (Ny - 1)*pygutter - (ny*Ny - 1)*ygutter)/float(ny*Ny)
459 for panel
in range(Nx*Ny):
463 for window
in range(nx*ny):
464 x = nx*px + window%nx
465 y = ny*py + window//nx
466 ax = fig.add_axes((x0 + xgutter + pxgutter + x*w + (px - 1)*pxgutter + (x - 1)*xgutter,
467 y0 + ygutter + pygutter + y*h + (py - 1)*pygutter + (y - 1)*ygutter,
468 w, h), frame_on=
True, facecolor=
'w')
474 matchKernelAmplitudes=False, keepPlots=True):
475 """Plot the PSF spatial model."""
479 import matplotlib.pyplot
as plt
480 import matplotlib
as mpl
481 except ImportError
as e:
482 log.warning(
"Unable to import matplotlib: %s", e)
485 noSpatialKernel = psf.getKernel()
492 for cell
in psfCellSet.getCellList():
493 for cand
in cell.begin(
False):
494 if not showBadCandidates
and cand.isBad():
498 im = cand.getMaskedImage()
506 for p, k
in zip(params, kernels):
507 amp += p * k.getSum()
509 targetFits = badFits
if cand.isBad()
else candFits
510 targetPos = badPos
if cand.isBad()
else candPos
511 targetAmps = badAmps
if cand.isBad()
else candAmps
513 targetFits.append([x / amp
for x
in params])
514 targetPos.append(candCenter)
515 targetAmps.append(amp)
517 xGood = numpy.array([pos.getX()
for pos
in candPos]) - exposure.getX0()
518 yGood = numpy.array([pos.getY()
for pos
in candPos]) - exposure.getY0()
519 zGood = numpy.array(candFits)
521 xBad = numpy.array([pos.getX()
for pos
in badPos]) - exposure.getX0()
522 yBad = numpy.array([pos.getY()
for pos
in badPos]) - exposure.getY0()
523 zBad = numpy.array(badFits)
526 xRange = numpy.linspace(0, exposure.getWidth(), num=numSample)
527 yRange = numpy.linspace(0, exposure.getHeight(), num=numSample)
529 kernel = psf.getKernel()
530 nKernelComponents = kernel.getNKernelParameters()
534 nPanelX = int(math.sqrt(nKernelComponents))
535 nPanelY = nKernelComponents//nPanelX
536 while nPanelY*nPanelX < nKernelComponents:
542 fig.canvas._tkcanvas._root().lift()
548 mpl.rcParams[
"figure.titlesize"] =
"x-small"
549 subplots =
makeSubplots(fig, 2, 2, Nx=nPanelX, Ny=nPanelY, xgutter=0.06, ygutter=0.06, pygutter=0.04)
551 for k
in range(nKernelComponents):
552 func = kernel.getSpatialFunction(k)
553 dfGood = zGood[:, k] - numpy.array([func(pos.getX(), pos.getY())
for pos
in candPos])
557 dfBad = zBad[:, k] - numpy.array([func(pos.getX(), pos.getY())
for pos
in badPos])
558 yMin =
min([yMin, dfBad.min()])
559 yMax =
max([yMax, dfBad.max()])
560 yMin -= 0.05 * (yMax - yMin)
561 yMax += 0.05 * (yMax - yMin)
566 fRange = numpy.ndarray((len(xRange), len(yRange)))
567 for j, yVal
in enumerate(yRange):
568 for i, xVal
in enumerate(xRange):
569 fRange[j][i] = func(xVal, yVal)
573 ax.set_autoscale_on(
False)
574 ax.set_xbound(lower=0, upper=exposure.getHeight())
575 ax.set_ybound(lower=yMin, upper=yMax)
576 ax.plot(yGood, dfGood,
'b+')
578 ax.plot(yBad, dfBad,
'r+')
580 ax.set_title(
'Residuals(y)')
584 if matchKernelAmplitudes
and k == 0:
591 norm = mpl.colors.Normalize(vmin=vmin, vmax=vmax)
592 im = ax.imshow(fRange, aspect=
'auto', origin=
"lower", norm=norm,
593 extent=[0, exposure.getWidth()-1, 0, exposure.getHeight()-1])
594 ax.set_title(
'Spatial poly')
595 plt.colorbar(im, orientation=
'horizontal', ticks=[vmin, vmax])
598 ax.set_autoscale_on(
False)
599 ax.set_xbound(lower=0, upper=exposure.getWidth())
600 ax.set_ybound(lower=yMin, upper=yMax)
601 ax.plot(xGood, dfGood,
'b+')
603 ax.plot(xBad, dfBad,
'r+')
605 ax.set_title(
'K%d Residuals(x)' % k)
609 photoCalib = exposure.getPhotoCalib()
611 if photoCalib.getCalibrationMean() <= 0:
614 ampMag = [photoCalib.instFluxToMagnitude(candAmp)
for candAmp
in candAmps]
615 ax.plot(ampMag, zGood[:, k],
'b+')
617 badAmpMag = [photoCalib.instFluxToMagnitude(badAmp)
for badAmp
in badAmps]
618 ax.plot(badAmpMag, zBad[:, k],
'r+')
620 ax.set_title(
'Flux variation')
625 if keepPlots
and not keptPlots:
628 print(
"%s: Please close plots when done." % __name__)
633 print(
"Plots closed, exiting...")
635 atexit.register(show)
639 def showPsf(psf, eigenValues=None, XY=None, normalize=True, display=None):
640 """Display a PSF's eigen images
642 If normalize is True, set the largest absolute value of each eigenimage to 1.0 (n.b. sum == 0.0 for i > 0)
648 coeffs = psf.getLocalKernel(
lsst.geom.PointD(XY[0], XY[1])).getKernelParameters()
652 mos = displayUtils.Mosaic(gutter=2, background=-0.1)
653 for i, k
in enumerate(psf.getKernel().getKernelList()):
654 im = afwImage.ImageD(k.getDimensions())
655 k.computeImage(im,
False)
657 im /= numpy.max(numpy.abs(im.getArray()))
660 mos.append(im,
"%g" % (coeffs[i]/coeffs[0]))
665 display = afwDisplay.Display()
666 mos.makeMosaic(display=display, title=
"Kernel Basis Functions")
671 def showPsfMosaic(exposure, psf=None, nx=7, ny=None, showCenter=True, showEllipticity=False,
672 showFwhm=False, stampSize=0, display=None, title=None):
673 """Show a mosaic of Psf images. exposure may be an Exposure (optionally with PSF),
674 or a tuple (width, height)
676 If stampSize is > 0, the psf images will be trimmed to stampSize*stampSize
681 showEllipticity =
True
682 scale = 2*math.log(2)
684 mos = displayUtils.Mosaic()
687 width, height = exposure.getWidth(), exposure.getHeight()
688 x0, y0 = exposure.getXY0()
690 psf = exposure.getPsf()
691 except AttributeError:
693 width, height = exposure[0], exposure[1]
696 raise RuntimeError(
"Unable to extract width/height from object of type %s" %
type(exposure))
699 ny = int(nx*float(height)/width + 0.5)
703 centroidName =
"SdssCentroid"
704 shapeName =
"base_SdssShape"
706 schema = afwTable.SourceTable.makeMinimalSchema()
707 schema.getAliasMap().
set(
"slot_Centroid", centroidName)
708 schema.getAliasMap().
set(
"slot_Centroid_flag", centroidName+
"_flag")
710 control = measBase.SdssCentroidControl()
711 centroider = measBase.SdssCentroidAlgorithm(control, centroidName, schema)
713 sdssShape = measBase.SdssShapeControl()
714 shaper = measBase.SdssShapeAlgorithm(sdssShape, shapeName, schema)
715 table = afwTable.SourceTable.make(schema)
717 table.defineCentroid(centroidName)
718 table.defineShape(shapeName)
723 if stampSize <= w
and stampSize <= h:
731 x = int(ix*(width-1)/(nx-1)) + x0
732 y = int(iy*(height-1)/(ny-1)) + y0
738 im = im.Factory(im, bbox)
739 lab =
"PSF(%d,%d)" % (x, y)
if False else ""
744 w, h = im.getWidth(), im.getHeight()
745 centerX = im.getX0() + w//2
746 centerY = im.getY0() + h//2
747 src = table.makeRecord()
750 foot.addPeak(centerX, centerY, 1)
751 src.setFootprint(foot)
754 centroider.measure(src, exp)
755 centers.append((src.getX() - im.getX0(), src.getY() - im.getY0()))
757 shaper.measure(src, exp)
758 shapes.append((src.getIxx(), src.getIxy(), src.getIyy()))
763 display = afwDisplay.Display()
764 mos.makeMosaic(display=display, title=title
if title
else "Model Psf", mode=nx)
766 if centers
and display:
767 with display.Buffering():
768 for i, (cen, shape)
in enumerate(zip(centers, shapes)):
769 bbox = mos.getBBox(i)
770 xc, yc = cen[0] + bbox.getMinX(), cen[1] + bbox.getMinY()
772 display.dot(
"+", xc, yc, ctype=afwDisplay.BLUE)
775 ixx, ixy, iyy = shape
779 display.dot(
"@:%g,%g,%g" % (ixx, ixy, iyy), xc, yc, ctype=afwDisplay.RED)
785 mimIn = exposure.getMaskedImage()
786 mimIn = mimIn.Factory(mimIn,
True)
788 psf = exposure.getPsf()
789 psfWidth, psfHeight = psf.getLocalKernel().getDimensions()
793 w, h = int(mimIn.getWidth()/scale), int(mimIn.getHeight()/scale)
795 im = mimIn.Factory(w + psfWidth, h + psfHeight)
799 x, y = s.getX(), s.getY()
801 sx, sy = int(x/scale + 0.5), int(y/scale + 0.5)
805 sim = smim.getImage()
809 flux = s.getApInstFlux()
810 elif magType ==
"model":
811 flux = s.getModelInstFlux()
812 elif magType ==
"psf":
813 flux = s.getPsfInstFlux()
815 raise RuntimeError(
"Unknown flux type %s" % magType)
818 except Exception
as e:
822 expIm = mimIn.getImage().
Factory(mimIn.getImage(),
824 int(y) - psfHeight//2),
830 cenPos.append([x - expIm.getX0() + sx, y - expIm.getY0() + sy])
835 display = afwDisplay.Display()
836 display.mtv(im, title=
"showPsfResiduals: image")
837 with display.Buffering():
839 display.dot(
"+", x, y)
845 """Write the contents of a SpatialCellSet to a many-MEF fits file"""
848 for cell
in psfCellSet.getCellList():
849 for cand
in cell.begin(
False):
855 md.set(
"CELL", cell.getLabel())
856 md.set(
"ID", cand.getId())
857 md.set(
"XCENTER", cand.getXCenter())
858 md.set(
"YCENTER", cand.getYCenter())
859 md.set(
"BAD", cand.isBad())
860 md.set(
"AMPL", cand.getAmplitude())
861 md.set(
"FLUX", cand.getSource().getPsfInstFlux())
862 md.set(
"CHI2", cand.getSource().getChi2())
864 im.writeFits(fileName, md, mode)
868 display.mtv(im, title=
"saveSpatialCellSet: image")
A Threshold is used to pass a threshold value to detection algorithms.
A compact representation of a collection of pixels.
The photometric calibration of an exposure.
A kernel that is a linear combination of fixed basis kernels.
A class that can be used to generate sequences of random numbers according to a number of different a...
Class for storing generic metadata.
An integer coordinate rectangle.
static Log getLogger(Log const &logger)
Provides consistent interface for LSST exceptions.
daf::base::PropertyList * list
daf::base::PropertySet * set
std::shared_ptr< FrameSet > append(FrameSet const &first, FrameSet const &second)
Construct a FrameSet that performs two transformations in series.
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
std::shared_ptr< Exposure< ImagePixelT, MaskPixelT, VariancePixelT > > makeExposure(MaskedImage< ImagePixelT, MaskPixelT, VariancePixelT > &mimage, std::shared_ptr< geom::SkyWcs const > wcs=std::shared_ptr< geom::SkyWcs const >())
A function to return an Exposure of the correct type (cf.
MaskedImage< ImagePixelT, MaskPixelT, VariancePixelT > * makeMaskedImage(typename std::shared_ptr< Image< ImagePixelT >> image, typename std::shared_ptr< Mask< MaskPixelT >> mask=Mask< MaskPixelT >(), typename std::shared_ptr< Image< VariancePixelT >> variance=Image< VariancePixelT >())
A function to return a MaskedImage of the correct type (cf.
int positionToIndex(double pos)
Convert image position to nearest integer index.
void randomGaussianImage(ImageT *image, Random &rand)
Set image to random numbers with a gaussian N(0, 1) distribution.
Statistics makeStatistics(lsst::afw::image::Image< Pixel > const &img, lsst::afw::image::Mask< image::MaskPixel > const &msk, int const flags, StatisticsControl const &sctrl=StatisticsControl())
Handle a watered-down front-end to the constructor (no variance)
std::shared_ptr< ImageT > offsetImage(ImageT const &image, float dx, float dy, std::string const &algorithmName="lanczos5", unsigned int buffer=0)
Return an image offset by (dx, dy) using the specified algorithm.
def makeSubplots(fig, nx=2, ny=2, Nx=1, Ny=1, plottingArea=(0.1, 0.1, 0.85, 0.80), pxgutter=0.05, pygutter=0.05, xgutter=0.04, ygutter=0.04, headroom=0.0, panelBorderWeight=0, panelColor='black')
def showPsfSpatialCells(exposure, psfCellSet, nMaxPerCell=-1, showChi2=False, showMoments=False, symb=None, ctype=None, ctypeUnused=None, ctypeBad=None, size=2, display=None)
def showPsf(psf, eigenValues=None, XY=None, normalize=True, display=None)
def plotPsfSpatialModel(exposure, psf, psfCellSet, showBadCandidates=True, numSample=128, matchKernelAmplitudes=False, keepPlots=True)
def splitId(oid, asDict=True)
def showPsfResiduals(exposure, sourceSet, magType="psf", scale=10, display=None)
def saveSpatialCellSet(psfCellSet, fileName="foo.fits", display=None)
def showSourceSet(sSet, xy0=(0, 0), display=None, ctype=afwDisplay.GREEN, symb="+", size=2)
def showPsfCandidates(exposure, psfCellSet, psf=None, display=None, normalize=True, showBadCandidates=True, fitBasisComponents=False, variance=None, chi=None)
def showPsfMosaic(exposure, psf=None, nx=7, ny=None, showCenter=True, showEllipticity=False, showFwhm=False, stampSize=0, display=None, title=None)
std::pair< std::vector< double >, afw::math::KernelList > fitKernelParamsToImage(afw::math::LinearCombinationKernel const &kernel, Image const &image, geom::Point2D const &pos)
Fit a LinearCombinationKernel to an Image, allowing the coefficients of the components to vary.
double subtractPsf(afw::detection::Psf const &psf, ImageT *data, double x, double y, double psfFlux=std::numeric_limits< double >::quiet_NaN())