Index: /branches/eam_branches/ipp-20110710/psphot/src/psphotSersicModelClass.c =================================================================== --- /branches/eam_branches/ipp-20110710/psphot/src/psphotSersicModelClass.c (revision 31966) +++ /branches/eam_branches/ipp-20110710/psphot/src/psphotSersicModelClass.c (revision 31966) @@ -0,0 +1,720 @@ +# include "psphotInternal.h" + +psF64 psPolynomial1DSolve(const psPolynomial1D* poly, psF64 y, bool upper); +void psphotSersicModelNorm (pmPCMdata *pcm, const pmSource *source); + +static psArray *classArray = NULL; +static float PSFratioMin = 1.0; +static float PSFratioMax = 4.0; +static float ASPratioMin = 0.2; +static float ASPratioMax = 1.0; + +// structure to define the sersic model guess information +typedef struct { + float PSFratioMin; + float PSFratioMax; + float ASPratioMin; + float ASPratioMax; + + float indexMin; + float indexMax; + float kronMin; + float kronMax; + + bool KronUpper; + psPolynomial1D *KronToIndex; + psPolynomial1D *IndexToTotalMag; + psPolynomial1D *IndexToMajor; + psPolynomial1D *IndexToMinor; +} psphotSersicModelClass; + +static void psphotSersicModelClassFree (psphotSersicModelClass *class) { + psFree(class->KronToIndex); + psFree(class->IndexToTotalMag); + psFree(class->IndexToMajor); + psFree(class->IndexToMinor); + return; +} + +psphotSersicModelClass *psphotSersicModelClassAlloc() { + + psphotSersicModelClass *class = (psphotSersicModelClass *) psAlloc(sizeof(psphotSersicModelClass)); + psMemSetDeallocator(class, (psFreeFunc) psphotSersicModelClassFree); + + class->PSFratioMin = NAN; + class->PSFratioMax = NAN; + class->ASPratioMin = NAN; + class->ASPratioMax = NAN; + + class->KronUpper = false; + class->KronToIndex = NULL; + class->IndexToTotalMag = NULL; + class->IndexToMajor = NULL; + class->IndexToMinor = NULL; + + return class; +} + +void psphotSersicModelClassSetIndexRange(psphotSersicModelClass *class, float indexMin, float indexMax) { + + class->indexMin = indexMin; + class->indexMax = indexMax; + float modelMinX = -0.5 * class->KronToIndex->coeff[1] / class->KronToIndex->coeff[2]; + float modelMinY = psPolynomial1DEval (class->KronToIndex, modelMinX); + float kronIndexMin = psPolynomial1DEval (class->KronToIndex, class->indexMin); + float kronIndexMax = psPolynomial1DEval (class->KronToIndex, class->indexMax); + + // if modelMinX is between indexMin and indexMax, use the modelMinY as the min/max kron value + // else, use the kron values at the min and max positions + + if ((modelMinX > class->indexMin) && (modelMinX < class->indexMax)) { + if (class->KronToIndex->coeff[2] < 0.0) { + class->kronMax = modelMinY - 0.001; // pad slightly to avoid falling off the curve + class->kronMin = PS_MIN(kronIndexMin, kronIndexMax); + } else { + class->kronMin = modelMinY + 0.001; // pad slightly to avoid falling off the curve + class->kronMax = PS_MAX(kronIndexMin, kronIndexMax); + } + } else { + class->kronMin = PS_MIN(kronIndexMin, kronIndexMax); + class->kronMax = PS_MAX(kronIndexMin, kronIndexMax); + } + + // KronUpper specifies which of the 2 quadratic solutions to accept: + if (modelMinX < class->indexMin) { + if (class->KronToIndex->coeff[2] < 0.0) { + class->KronUpper = false; + } else { + class->KronUpper = true; + } + return; + } + if (modelMinX > class->indexMax) { + if (class->KronToIndex->coeff[2] < 0.0) { + class->KronUpper = true; + } else { + class->KronUpper = false; + } + return; + } + + if (fabs(modelMinX - class->indexMin) < fabs(modelMinX - class->indexMax)) { + if (class->KronToIndex->coeff[2] < 0.0) { + class->KronUpper = false; + } else { + class->KronUpper = true; + } + return; + } + if (class->KronToIndex->coeff[2] < 0.0) { + class->KronUpper = true; + } else { + class->KronUpper = false; + } + return; +} + +void psphotSersicModelClassInit () { + + psphotSersicModelClass *class; + + if (classArray) return; + + // hardwired trends for now (move into the recipe?) + classArray = psArrayAllocEmpty (4); + + // image.00.01.fit.dat: PSFratio : 1.0, ASPratio : 1.0 + class = psphotSersicModelClassAlloc(); + + class->PSFratioMin = 1.0; + class->PSFratioMax = 1.5; + class->ASPratioMin = 0.71; + class->ASPratioMax = 1.00; + class->KronToIndex = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2); + class->IndexToTotalMag = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 1); + class->IndexToMajor = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 1); + class->IndexToMinor = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 1); + + // major: C0: -0.430070611162, C1: 0.0511725134843 + class->IndexToMajor->coeff[0] = -0.430070611162; + class->IndexToMajor->coeff[1] = 0.0511725134843; + + // minor: C0: -0.250221103738, C1: 0.0136204422354 + class->IndexToMinor->coeff[0] = -0.250221103738; + class->IndexToMinor->coeff[1] = 0.0136204422354; + + // MagOffset: C0: -0.0718507604436, C1: -0.0518624470415 + class->IndexToTotalMag->coeff[0] = -0.0718507604436; + class->IndexToTotalMag->coeff[1] = -0.0518624470415; + + // KronMag: C0: -1.80936771396, C1: 0.344933296013, C2: -0.0314000083805 + class->KronToIndex->coeff[0] = -1.80936771396; + class->KronToIndex->coeff[1] = 0.344933296013; + class->KronToIndex->coeff[2] = -0.0314000083805; + + psphotSersicModelClassSetIndexRange(class, 0.5, 6.0); + + psArrayAdd(classArray, 4, class); + psFree (class); + + // image.00.02.fit.dat: PSFratio : 2.0, ASPratio : 1.0 + class = psphotSersicModelClassAlloc(); + + class->PSFratioMin = 1.5; + class->PSFratioMax = 4.0; // note : missing 4.0, 1.0 + class->ASPratioMin = 0.71; + class->ASPratioMax = 1.00; + class->KronToIndex = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2); + class->IndexToTotalMag = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 1); + class->IndexToMajor = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 1); + class->IndexToMinor = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 1); + + // major: C0: -0.0611883230251, C1: 0.0498266924052 + class->IndexToMajor->coeff[0] = -0.0611883230251; + class->IndexToMajor->coeff[1] = 0.0498266924052; + + // minor: C0: -0.0471497475843, C1: 0.0498093464447 + class->IndexToMinor->coeff[0] = -0.0471497475843; + class->IndexToMinor->coeff[1] = 0.0498093464447; + + // MagOffset: C0: 0.133350843774, C1: -0.147490621217 + class->IndexToTotalMag->coeff[0] = 0.133350843774; + class->IndexToTotalMag->coeff[1] = -0.147490621217; + + // KronMag: C0: -3.81043970737, C1: 1.11738875295, C2: -0.121059914305 + class->KronToIndex->coeff[0] = -3.81043970737; + class->KronToIndex->coeff[1] = 1.11738875295; + class->KronToIndex->coeff[2] = -0.121059914305; + + psphotSersicModelClassSetIndexRange(class, 0.5, 6.0); + + psArrayAdd(classArray, 4, class); + psFree (class); + + // image.01.01.fit.dat: PSFratio : 1.0, ASPratio : 0.5 + class = psphotSersicModelClassAlloc(); + + class->PSFratioMin = 1.0; + class->PSFratioMax = 1.5; + class->ASPratioMin = 0.41; + class->ASPratioMax = 0.71; + class->KronToIndex = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2); + class->IndexToTotalMag = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 1); + class->IndexToMajor = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 1); + class->IndexToMinor = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 1); + + // major: C0: -0.315072669277, C1: 0.00877759328113 + class->IndexToMajor->coeff[0] = -0.315072669277; + class->IndexToMajor->coeff[1] = 0.00877759328113; + + // minor: C0: -0.315072669277, C1: 0.00877759328113 + class->IndexToMinor->coeff[0] = -0.315072669277; + class->IndexToMinor->coeff[1] = 0.00877759328113; + + // MagOffset: C0: -0.441283222989, C1: 0.0531287339075 + class->IndexToTotalMag->coeff[0] = -0.441283222989; + class->IndexToTotalMag->coeff[1] = 0.0531287339075; + + // KronMag: C0: -0.819564188427, C1: -0.0333590360661, C2: 0.00814022994725 + // XXX use a linear fit? + class->KronToIndex->coeff[0] = -0.819564188427; + class->KronToIndex->coeff[1] = -0.0333590360661; + class->KronToIndex->coeff[2] = 0.00814022994725; + + psphotSersicModelClassSetIndexRange(class, 0.5, 6.0); + + psArrayAdd(classArray, 4, class); + psFree (class); + + // image.01.02.fit.dat: PSFratio : 2.0, ASPratio : 0.5 + class = psphotSersicModelClassAlloc(); + + class->PSFratioMin = 1.5; + class->PSFratioMax = 3.0; + class->ASPratioMin = 0.41; + class->ASPratioMax = 0.71; + class->KronToIndex = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2); + class->IndexToTotalMag = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 1); + class->IndexToMajor = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 1); + class->IndexToMinor = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 1); + + // major: C0: 0.0394644381614, C1: -0.00261362634784 + class->IndexToMajor->coeff[0] = 0.0394644381614; + class->IndexToMajor->coeff[1] = -0.00261362634784; + + // minor: C0: -0.314123499959, C1: 0.0215627071372 + class->IndexToMinor->coeff[0] = -0.314123499959; + class->IndexToMinor->coeff[1] = 0.0215627071372; + + // MagOffset: C0: 0.00475524952411, C1: -0.0960576544787 + class->IndexToTotalMag->coeff[0] = 0.00475524952411; + class->IndexToTotalMag->coeff[1] = -0.0960576544787; + + // KronMag: C0: -2.77961654208, C1: 0.719625939973, C2: -0.0772103821458 + class->KronToIndex->coeff[0] = -2.77961654208; + class->KronToIndex->coeff[1] = 0.719625939973; + class->KronToIndex->coeff[2] = -0.0772103821458; + + psphotSersicModelClassSetIndexRange(class, 0.5, 6.0); + + psArrayAdd(classArray, 4, class); + psFree (class); + + // image.01.03.fit.dat: PSFratio : 4.0, ASPratio : 0.5 + class = psphotSersicModelClassAlloc(); + + class->PSFratioMin = 3.0; + class->PSFratioMax = 4.0; + class->ASPratioMin = 0.41; + class->ASPratioMax = 0.71; + class->KronToIndex = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2); + class->IndexToTotalMag = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 1); + class->IndexToMajor = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 1); + class->IndexToMinor = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 1); + + // major: C0: 0.0966828791109, C1: 0.0508749927221 + class->IndexToMajor->coeff[0] = 0.0966828791109; + class->IndexToMajor->coeff[1] = 0.0508749927221; + + // minor: C0: -0.163878057365, C1: 0.0621817242826 + class->IndexToMinor->coeff[0] = -0.163878057365; + class->IndexToMinor->coeff[1] = 0.0621817242826; + + // MagOffset: C0: -0.0584285249407, C1: -0.115758644152 + class->IndexToTotalMag->coeff[0] = -0.0584285249407; + class->IndexToTotalMag->coeff[1] = -0.115758644152; + + // KronMag: C0: -4.41677343695, C1: 1.23569951267, C2: -0.131733958476 + class->KronToIndex->coeff[0] = -4.41677343695; + class->KronToIndex->coeff[1] = 1.23569951267; + class->KronToIndex->coeff[2] = -0.131733958476; + + psphotSersicModelClassSetIndexRange(class, 0.5, 6.0); + + psArrayAdd(classArray, 4, class); + psFree (class); + + // image.02.01.fit.dat: PSFratio : 1.0, ASPratio : 0.33 + class = psphotSersicModelClassAlloc(); + + class->PSFratioMin = 1.0; + class->PSFratioMax = 1.5; + class->ASPratioMin = 0.25; + class->ASPratioMax = 0.41; + class->KronToIndex = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2); + class->IndexToTotalMag = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 1); + class->IndexToMajor = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 1); + class->IndexToMinor = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 1); + + // major: C0: -0.154899753946, C1: -0.0193502412043 + class->IndexToMajor->coeff[0] = -0.154899753946; + class->IndexToMajor->coeff[1] = -0.0193502412043; + + // minor: C0: -0.154899753946, C1: -0.0193502412043 + class->IndexToMinor->coeff[0] = -0.154899753946; + class->IndexToMinor->coeff[1] = -0.0193502412043; + + // MagOffset: C0: -0.696799097791, C1: 0.098950200942 + class->IndexToTotalMag->coeff[0] = -0.696799097791; + class->IndexToTotalMag->coeff[1] = 0.098950200942; + + // KronMag: C0: -0.567489477137, C1: 0.00400934329328, C2: -0.00641637764069 + // XXX use a linear fit? + class->KronToIndex->coeff[0] = -0.567489477137; + class->KronToIndex->coeff[1] = 0.00400934329328; + class->KronToIndex->coeff[2] = -0.00641637764069; + + psphotSersicModelClassSetIndexRange(class, 0.5, 6.0); + + psArrayAdd(classArray, 4, class); + psFree (class); + + // image.02.02.fit.dat: PSFratio : 2.0, ASPratio : 0.33 + class = psphotSersicModelClassAlloc(); + + class->PSFratioMin = 1.5; + class->PSFratioMax = 3.0; + class->ASPratioMin = 0.25; + class->ASPratioMax = 0.41; + class->KronToIndex = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2); + class->IndexToTotalMag = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 1); + class->IndexToMajor = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 1); + class->IndexToMinor = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 1); + + // major: C0: -0.227923466363, C1: 0.078789082243 + class->IndexToMajor->coeff[0] = -0.227923466363; + class->IndexToMajor->coeff[1] = 0.078789082243; + + // minor: C0: -0.39010838322, C1: -0.0191311211209 + class->IndexToMinor->coeff[0] = -0.39010838322; + class->IndexToMinor->coeff[1] = -0.0191311211209; + + // MagOffset: C0: -0.196058611899, C1: -0.0364610432252 + class->IndexToTotalMag->coeff[0] = -0.196058611899; + class->IndexToTotalMag->coeff[1] = -0.0364610432252; + + // KronMag: C0: -1.94707725307, C1: 0.296502645791, C2: -0.0225310468023 + class->KronToIndex->coeff[0] = -1.94707725307; + class->KronToIndex->coeff[1] = 0.296502645791; + class->KronToIndex->coeff[2] = -0.0225310468023; + + psphotSersicModelClassSetIndexRange(class, 0.5, 6.0); + + psArrayAdd(classArray, 4, class); + psFree (class); + + // image.02.03.fit.dat: PSFratio : 4.0, ASPratio : 0.33 + class = psphotSersicModelClassAlloc(); + + class->PSFratioMin = 3.0; + class->PSFratioMax = 4.0; + class->ASPratioMin = 0.25; + class->ASPratioMax = 0.41; + class->KronToIndex = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2); + class->IndexToTotalMag = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 1); + class->IndexToMajor = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 1); + class->IndexToMinor = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 1); + + // major: C0: -0.00875235291975, C1: 0.0644897947208 + class->IndexToMajor->coeff[0] = -0.00875235291975; + class->IndexToMajor->coeff[1] = 0.0644897947208; + + // minor: C0: -0.250701690848, C1: 0.0431820606974 + class->IndexToMinor->coeff[0] = -0.250701690848; + class->IndexToMinor->coeff[1] = 0.0431820606974; + + // MagOffset: C0: -0.158069871256, C1: -0.0811671351817 + class->IndexToTotalMag->coeff[0] = -0.158069871256; + class->IndexToTotalMag->coeff[1] = -0.0811671351817; + + // KronMag: C0: -3.6119244064, C1: 0.8631113861, C2: -0.0853296107645 + class->KronToIndex->coeff[0] = -3.6119244064; + class->KronToIndex->coeff[1] = 0.8631113861; + class->KronToIndex->coeff[2] = -0.0853296107645; + + psphotSersicModelClassSetIndexRange(class, 0.5, 6.0); + + psArrayAdd(classArray, 4, class); + psFree (class); + + // image.03.01.fit.dat: PSFratio : 1.0, ASPratio : 0.2 + class = psphotSersicModelClassAlloc(); + + class->PSFratioMin = 1.0; + class->PSFratioMax = 1.5; + class->ASPratioMin = 0.20; + class->ASPratioMax = 0.25; + class->KronToIndex = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2); + class->IndexToTotalMag = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 1); + class->IndexToMajor = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 1); + class->IndexToMinor = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 1); + + // major: C0: -0.212431708382, C1: 0.00426679037981 + class->IndexToMajor->coeff[0] = -0.212431708382; + class->IndexToMajor->coeff[1] = 0.00426679037981; + + // minor: C0: -0.212431708382, C1: 0.00426679037981 + class->IndexToMinor->coeff[0] = -0.212431708382; + class->IndexToMinor->coeff[1] = 0.00426679037981; + + // MagOffset: C0: -0.62070857795, C1: 0.000773428325368 + class->IndexToTotalMag->coeff[0] = -0.62070857795; + class->IndexToTotalMag->coeff[1] = 0.000773428325368; + + // KronMag: C0: -0.460478728787, C1: 0.00587256401903, C2: -0.002344873666 + class->KronToIndex->coeff[0] = -0.460478728787; + class->KronToIndex->coeff[1] = 0.00587256401903; + class->KronToIndex->coeff[2] = -0.002344873666; + + psphotSersicModelClassSetIndexRange(class, 0.5, 6.0); + + psArrayAdd(classArray, 4, class); + psFree (class); + + // image.03.02.fit.dat: PSFratio : 2.0, ASPratio : 0.2 + class = psphotSersicModelClassAlloc(); + + class->PSFratioMin = 1.5; + class->PSFratioMax = 3.0; + class->ASPratioMin = 0.20; + class->ASPratioMax = 0.25; + class->KronToIndex = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2); + class->IndexToTotalMag = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 1); + class->IndexToMajor = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 1); + class->IndexToMinor = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 1); + + // major: C0: -0.385105744329, C1: 0.066302247629 + class->IndexToMajor->coeff[0] = -0.385105744329; + class->IndexToMajor->coeff[1] = 0.066302247629; + + // minor: C0: -0.385105744329, C1: 0.066302247629 + class->IndexToMinor->coeff[0] = -0.385105744329; + class->IndexToMinor->coeff[1] = 0.066302247629; + + // MagOffset: C0: -0.359723538992, C1: 0.00907435268225 + class->IndexToTotalMag->coeff[0] = -0.359723538992; + class->IndexToTotalMag->coeff[1] = 0.00907435268225; + + // KronMag: C0: -1.42088663998, C1: 0.0988810020043, C2: -0.00428297294884 + class->KronToIndex->coeff[0] = -1.42088663998; + class->KronToIndex->coeff[1] = 0.0988810020043; + class->KronToIndex->coeff[2] = -0.00428297294884; + + psphotSersicModelClassSetIndexRange(class, 0.5, 6.0); + + psArrayAdd(classArray, 4, class); + psFree (class); + + // image.03.03.fit.dat: PSFratio : 4.0, ASPratio : 0.2 + class = psphotSersicModelClassAlloc(); + + class->PSFratioMin = 3.0; + class->PSFratioMax = 4.0; + class->ASPratioMin = 0.20; + class->ASPratioMax = 0.25; + class->KronToIndex = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2); + class->IndexToTotalMag = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 1); + class->IndexToMajor = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 1); + class->IndexToMinor = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 1); + + // major: C0: -0.132561970043, C1: 0.064707352393 + class->IndexToMajor->coeff[0] = -0.132561970043; + class->IndexToMajor->coeff[1] = 0.064707352393; + + // minor: C0: -0.195782196124, C1: -0.0457728222883 + class->IndexToMinor->coeff[0] = -0.195782196124; + class->IndexToMinor->coeff[1] = -0.0457728222883; + + // MagOffset: C0: -0.208808020613, C1: -0.0605412708629 + class->IndexToTotalMag->coeff[0] = -0.208808020613; + class->IndexToTotalMag->coeff[1] = -0.0605412708629; + + // KronMag: C0: -3.01650012413, C1: 0.618382890037, C2: -0.0615504992544 + class->KronToIndex->coeff[0] = -3.01650012413; + class->KronToIndex->coeff[1] = 0.618382890037; + class->KronToIndex->coeff[2] = -0.0615504992544; + + psphotSersicModelClassSetIndexRange(class, 0.5, 6.0); + + psArrayAdd(classArray, 4, class); + psFree (class); +} + +void psphotSersicModelClassCleanup () { + psFree (classArray); +} + +void psphotWriteGuess(pmModel *model, float dKronMag, float PSFratio, float ASPratio); + +bool psphotSersicModelClassGuessPCM (pmPCMdata *pcm, pmSource *source) { + + psF32 *PAR = pcm->modelConv->params->data.F32; + + // XXX require: moments, psfMag, moments->kronFlux + + // we attempt to guess the Sersic model parameters based on a number of non-parametric + // measurements: moments, kron mag, etc + + // * first choose the model class based on (a) ratio of Moments Major axis to PSF major + // * axis and (b) moments axial ratio + + // convert the moments to Major,Minor,Theta + psEllipseMoments moments; + + moments.x2 = source->moments->Mxx; + moments.y2 = source->moments->Myy; + moments.xy = source->moments->Mxy; + + // limit axis ratio < 20.0 + psEllipseAxes momentAxes = psEllipseMomentsToAxes (moments, 20.0); + + // convert the PSF shape to Major,Minor,Theta + psEllipseShape shape; + + // XXX make sure this is consistent with the re-definition of PM_PAR_SXX + shape.sx = source->modelPSF->params->data.F32[PM_PAR_SXX] / M_SQRT2; + shape.sy = source->modelPSF->params->data.F32[PM_PAR_SYY] / M_SQRT2; + shape.sxy = source->modelPSF->params->data.F32[PM_PAR_SXY]; + psEllipseAxes psfAxes = psEllipseShapeToAxes (shape, 20.0); + + // get the PSFratio and the ASPratio : we use these to choose the model class + float PSFratio = momentAxes.major / (2.35 * psfAxes.major) ; + float ASPratio = momentAxes.minor / momentAxes.major; + + // saturate PSFratio and ASPratio to min / max values + PSFratio = PS_MAX(PSFratio, PSFratioMin); + PSFratio = PS_MIN(PSFratio, PSFratioMax); + ASPratio = PS_MAX(ASPratio, ASPratioMin); + ASPratio = PS_MIN(ASPratio, ASPratioMax); + + // find the containing model class: + + psphotSersicModelClass *class = NULL; + for (int i = 0; i < classArray->n; i++) { + psphotSersicModelClass *thisClass = classArray->data[i]; + if (PSFratio < thisClass->PSFratioMin) continue; + if (PSFratio > thisClass->PSFratioMax) continue; + if (ASPratio < thisClass->ASPratioMin) continue; + if (ASPratio > thisClass->ASPratioMax) continue; + class = thisClass; + break; + } + + psAssert (class, "PSFratio and ASPratio must be in range"); + + // get the index guess from the KronMag - psfMag: + + // get dKronMag & saturate dKronMag at limits of valid range + float dKronMag = -2.5*log10(source->moments->KronFlux) - source->psfMag; + dKronMag = PS_MIN(dKronMag, class->kronMax); + dKronMag = PS_MAX(dKronMag, class->kronMin); + + // get index (saturate at valid ends) + float index = psPolynomial1DSolve (class->KronToIndex, dKronMag, class->KronUpper); + index = PS_MIN (index, class->indexMax); + index = PS_MAX (index, class->indexMin); + + // float totalMagOffset = psPolynomial1DEval (class->IndexToTotalMag, index); + // float totalMag = -2.5*log10(source->moments->KronFlux) + totalMagOffset; + + // need to go from totalMag to Io in the sersic model + + float majorAxisFactor = psPolynomial1DEval (class->IndexToMajor, index); + float minorAxisFactor = psPolynomial1DEval (class->IndexToMinor, index); + + momentAxes.major *= (1.0 + majorAxisFactor); + momentAxes.minor *= (1.0 + minorAxisFactor); + + psEllipseShape extShapeGuess = psEllipseAxesToShape (momentAxes); + if (!isfinite(extShapeGuess.sx)) return false; + if (!isfinite(extShapeGuess.sy)) return false; + if (!isfinite(extShapeGuess.sxy)) return false; + + // set the actual model parameters: + + // index is standard sersic index (n = 1-4), but PAR7 = 1/2n + PAR[PM_PAR_7] = 0.5 / index; + + // sky is zero (no longer fitted, but not yet deprecated) + PAR[PM_PAR_SKY] = 0.0; + + // set the model position + if (!pmModelSetPosition(&PAR[PM_PAR_XPOS], &PAR[PM_PAR_YPOS], source)) { + return false; + } + PAR[PM_PAR_SXX] = extShapeGuess.sx * M_SQRT2; + PAR[PM_PAR_SXY] = extShapeGuess.sxy; + PAR[PM_PAR_SYY] = extShapeGuess.sy * M_SQRT2; + + // XXX this is a bit of a waste: calculate the flux with Io = 1.0, then renormalize. + PAR[PM_PAR_I0] = 1.0; + + // set the normalization by linear fit between model and data + psphotSersicModelNorm (pcm, source); + + // float flux = pcm->modelConv->modelFlux(pcm->modelConv->params); + // float normMag = -2.5*log10(flux); + // float Io = pow(10.0, -0.4*(totalMag - normMag)); + // PAR[PM_PAR_I0] = Io; + + psphotWriteGuess(pcm->modelConv, dKronMag, PSFratio, ASPratio); + + return true; +} + +void psphotSersicModelNorm (pmPCMdata *pcm, const pmSource *source) { + + psVector *params = pcm->modelConv->params; + // XXX : not needed? psAssert (params->data.F32[PM_PAR_I0] == 1.0, "not normalized?"); + + // generate the convolved model image + // working vector to store local coordinate + psVector *coord = psVectorAlloc(2, PS_TYPE_F32); + + // create the convolved model in situ + psAssert (pcm->modelConvFlux, "not already allocated?"); + psImageInit (pcm->modelConvFlux, 0.0); + + // fill in the coordinate and value entries + for (psS32 i = 0; i < source->pixels->numRows; i++) { + for (psS32 j = 0; j < source->pixels->numCols; j++) { + + // Convert i/j to image space: + coord->data.F32[0] = (psF32) (j + source->pixels->col0); + coord->data.F32[1] = (psF32) (i + source->pixels->row0); + + pcm->modelConvFlux->data.F32[i][j] = pcm->modelConv->modelFunc (NULL, params, coord); + } + } + psFree(coord); + + psImageSmooth (pcm->modelConvFlux, pcm->sigma, pcm->nsigma); + + float YYmod = 0.0; + float Ymod2 = 0.0; + bool usePoisson = false; + + for (psS32 i = 0; i < source->pixels->numRows; i++) { + for (psS32 j = 0; j < source->pixels->numCols; j++) { + // XXX are we doing the right thing with the mask? + // skip masked points + if (source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[i][j]) { + continue; + } + // skip zero-variance points + if (source->variance->data.F32[i][j] == 0) { + continue; + } + // skip nan value points + if (!isfinite(source->pixels->data.F32[i][j])) { + continue; + } + + float ymodel = pcm->modelConvFlux->data.F32[i][j]; + float yweight = (usePoisson) ? 1.0 / source->variance->data.F32[i][j] : 1.0; + + YYmod += ymodel * source->pixels->data.F32[i][j] * yweight; + Ymod2 += PS_SQR(ymodel) * yweight; + } + } + + float Io = YYmod / Ymod2; + + params->data.F32[PM_PAR_I0] *= Io; + + return; +} + +// Given a polynomial y = C0 + C1 x + C2 x^2 etc, solve for x given y +// actually: this only solves 1st and 2nd order polynomials +// for 2nd order, it always returns the (positive/negative) term + +psF64 psPolynomial1DSolve( + const psPolynomial1D* poly, + psF64 y, + bool upper) +{ + psF64 x, C; + switch (poly->nX) { + case 1: + // y = coeff[0] + coeff[1]*x + // x = (y - coeff[0]) / coeff[1] + x = (y - poly->coeff[0]) / poly->coeff[1]; + return x; + + case 2: + // y = coeff[0] + coeff[1]*x + coeff[2]*x^2 + // x = -coeff[1] +/- sqrt(coeff[1]^2 - 4*coeff[0]*coeff[2]) / (2 coeff[0]) + C = poly->coeff[0] - y; + if (upper) { + x = (-poly->coeff[1] + sqrt(PS_SQR(poly->coeff[1]) - 4*poly->coeff[2]*C)) / (2.0 * poly->coeff[2]); + } else { + x = (-poly->coeff[1] - sqrt(PS_SQR(poly->coeff[1]) - 4*poly->coeff[2]*C)) / (2.0 * poly->coeff[2]); + } + return x; + + default: + psAbort("invalid polynomial for 1D solver"); + } + psAbort("invalid polynomial for 1D solver"); +}