Skip to content

GitLab

Commit 83060a1a authored by Johan Richard's avatar Johan Richard
Browse files

Added extended contours in shapemodel BC

parent 0975ff3a
Hide whitespace changes
Inline Side-by-side
src/keep_cl.c
View file @ 83060a1a
......@@ -210,3 +210,83 @@ static int readlist(struct point I[NPOINT], char file[NIMAX])
return(k);
}
/****************************************************************
* In the cleanlens mode 2, grow the contour inside which the profile
* will be computed according to the psf
* Input
* - the pixlist containing the points inside the contours
* provided by the user
* Return
* - pl_ext the pixlist of points in the extended contour
* - npl_ext the number of elements in pl_ext
*/
void keep_cl_extended(struct pixlist *pl, int npl, struct pixlist *pl_ext, int *npl_ext)
{
extern struct g_mode M;
extern struct g_pixel imFrame;
int n, j, k;
double **tmp;
// allocate temporary array to compute chi2 in o_chi()
tmp = (double **) alloc_square_double(imFrame.ny, imFrame.nx);
// set to 1 the pixels that are inside the contour
for(n = 0; n < npl; n++)
{
j = pl[n].ii;
k = pl[n].jj;
tmp[j][k] = 1.0;
}
// convolve array with psf
d_seeing_omp(tmp, imFrame.nx, imFrame.ny, imFrame.pixelx);
// debug output image
//char *outfile_ext = "tmp.fits";
// wrf_fits(outfile_ext, tmp, imFrame.nx, imFrame.ny, imFrame.xmin, imFrame.xmax, imFrame.ymin, imFrame.ymax);
//exit(-1);
// get pixels that are above threshold in the convolved array
const double ZERO_TH = 1e-10;
int kk = 0;
double xpixel, ypixel, xw, yw;
struct point P;
if (imFrame.wcsinfo != NULL)
{
for (j = 0; j < imFrame.ny; j++)
{
for (k = 0; k < imFrame.nx; k++)
{
if(tmp[j][k] > ZERO_TH)
{
xpixel = (double)(k + 1);
ypixel = (double)(j + 1);
pix2wcs(imFrame.wcsinfo, xpixel, ypixel, &xw, &yw);
P.x = xw - M.ref_ra;
P.x *= -3600.*cos(M.ref_dec * DTR);
P.y = yw - M.ref_dec;
P.y *= 3600.;
pl_ext[kk].i = P.y;
pl_ext[kk].j = P.x;
pl_ext[kk].ii = j;
pl_ext[kk].jj = k;
kk++;
}
}
}
}
else /* TODO: take this case into account if necessary */
{
FPRINTF(stderr, "ERROR: no wcs in input image %s\n",imFrame.pixfile);
exit(-1);
}
NPRINTF(stderr, "KEEP: %d pixels in extended contour\n", kk);
*npl_ext = kk;
free_square_double(tmp, imFrame.ny);
}
src/o_chi.c
View file @ 83060a1a
......@@ -93,8 +93,11 @@ int o_chi_lhood0(double *chi_ext, double *lhood0_ext, double *np_b0)
// mod bclement shapemodel_contour
extern struct pixlist *imo_pl;
extern struct pixlist *imo_pl_dpl;
extern int imo_npl;
// mod bclement shapemodel_contour_ext
extern struct pixlist *imo_pl_ext;
extern struct pixlist *imo_pl_dpl_ext;
extern int imo_npl_ext;
int n;
/* local variables */
......@@ -192,7 +195,7 @@ int o_chi_lhood0(double *chi_ext, double *lhood0_ext, double *np_b0)
// mod bclement shapemodel_contour
//o_pixel(ero, imFrame.nx, imFrame.ny, imFrame.pixelx, imFrame.xmin, imFrame.xmax, imFrame.ymin, imFrame.ymax, source, map_axx, map_ayy);
o_pixel_contour(ero, imFrame.nx, imFrame.ny, imFrame.pixelx, imFrame.xmin, imFrame.xmax, imFrame.ymin, imFrame.ymax, source, imo_pl, imo_pl_dpl, imo_npl); // mod bclement shapemodel_contour
o_pixel_contour(ero, imFrame.nx, imFrame.ny, imFrame.pixelx, imFrame.xmin, imFrame.xmax, imFrame.ymin, imFrame.ymax, source, imo_pl_ext, imo_pl_dpl_ext, imo_npl_ext); // mod bclement shapemodel_contour_ext
if (O.setseeing)
d_seeing_omp(ero, imFrame.nx, imFrame.ny, imFrame.pixelx);
......
src/o_chires.c
View file @ 83060a1a
......@@ -93,8 +93,11 @@ int o_chi_lhood0(double *chi_ext, double *lhood0_ext, double *np_b0)
// mod bclement shapemodel_contour
extern struct pixlist *imo_pl;
extern struct pixlist *imo_pl_dpl;
extern int imo_npl;
// mod bclement shapemodel_contour_ext
extern struct pixlist *imo_pl_ext;
extern struct pixlist *imo_pl_dpl_ext;
extern int imo_npl_ext;
int n;
/* local variables */
......@@ -192,7 +195,7 @@ int o_chi_lhood0(double *chi_ext, double *lhood0_ext, double *np_b0)
// mod bclement shapemodel_contour
//o_pixel(ero, imFrame.nx, imFrame.ny, imFrame.pixelx, imFrame.xmin, imFrame.xmax, imFrame.ymin, imFrame.ymax, source, map_axx, map_ayy);
o_pixel_contour(ero, imFrame.nx, imFrame.ny, imFrame.pixelx, imFrame.xmin, imFrame.xmax, imFrame.ymin, imFrame.ymax, source, imo_pl, imo_pl_dpl, imo_npl); // mod bclement shapemodel_contour
o_pixel_contour(ero, imFrame.nx, imFrame.ny, imFrame.pixelx, imFrame.xmin, imFrame.xmax, imFrame.ymin, imFrame.ymax, source, imo_pl_ext, imo_pl_dpl_ext, imo_npl_ext); // mod bclement shapemodel_contour_ext
if (O.setseeing)
d_seeing_omp(ero, imFrame.nx, imFrame.ny, imFrame.pixelx);
......
src/o_global.c
View file @ 83060a1a
......@@ -54,8 +54,11 @@ double *np_b0; // non parametric accelerator
// mod bclement shapemodel_contour
struct pixlist *imo_pl; // list of points inside the polygons for shapemodel
struct pixlist *imo_pl_dpl; // list of dpl inside the polygons for shapemodel
int imo_npl;
// mod bclement shapemodel_contour_ext
struct pixlist *imo_pl_ext; // list of points inside the psf-extended polygons for shapemodel
struct pixlist *imo_pl_dpl_ext; // list of dpl inside the psf_extended polygons for shapemodel
int imo_npl_ext;
/****************************************************************/
/* nom: global */
......@@ -278,7 +281,8 @@ void o_global_free()
{
free_square_double(ero, imFrame.ny);
free(imo_pl); // mod bclement shapemodel_contour
free(imo_pl_dpl); // mod bclement shapemodel_contour
free(imo_pl_ext); // mod bclement shapemodel_contour_ext
free(imo_pl_dpl_ext); // mod bclement shapemodel_contour_ext
}
else if (M.iclean == 3 )
{
......@@ -472,6 +476,8 @@ void readConstraints()
// Initialize and build the pixlist defining the contours
if(M.iclean == 2)
{
const extern struct g_observ O; // mod bclement shapemodel_contour_ext
/* Allocate the maximum number of points to send to source plane */
imo_pl = (struct pixlist*)malloc((unsigned int)
imFrame.nx * imFrame.ny * sizeof(struct pixlist));
......@@ -485,6 +491,13 @@ void readConstraints()
else
{
keep_cl(imo, imo_pl, &imo_npl);
// mod bclement shapemodel_contour_ext
if (O.setseeing)
{
imo_pl_ext = (struct pixlist*)malloc((unsigned int)
imFrame.nx * imFrame.ny * sizeof(struct pixlist));
keep_cl_extended(imo_pl, imo_npl, imo_pl_ext, &imo_npl_ext);
}
}
// TODO : remove duplicate pixels in list when overlapping contours occur
......@@ -494,6 +507,9 @@ void readConstraints()
imo_pl = (struct pixlist*)realloc(imo_pl, (unsigned int)
imo_npl * sizeof(struct pixlist));
// mod bclement shapemodel_contour_ext
imo_pl_ext = (struct pixlist*)realloc(imo_pl_ext, (unsigned int)
imo_npl_ext * sizeof(struct pixlist));
// redo the pre-computation of the loglikelihood normalization factor
lhood_wo = 0.;
......@@ -530,22 +546,22 @@ void readConstraints()
if( G.no_lens == 0 && pot_nopt == 0 )
{
//mod bclement shapemodel_contour
//mod bclement shapemodel_contour_ext
struct point pi, ps;
double zs_ref = 100.;
double dlsds_ref = dratio(lens[0].z, zs_ref);
if(M.iclean == 2) // compute normalized dpl only inside contour
if(M.iclean == 2) // compute normalized dpl only inside extended contour
{
imo_pl_dpl = (struct pixlist*)malloc((unsigned int)
imo_npl * sizeof(struct pixlist));
imo_pl_dpl_ext = (struct pixlist*)malloc((unsigned int)
imo_npl_ext * sizeof(struct pixlist));
int n;
for(n = 0; n < imo_npl; n++)
for(n = 0; n < imo_npl_ext; n++)
{
pi.y = imo_pl[n].i;
pi.x = imo_pl[n].j;
pi.y = imo_pl_ext[n].i;
pi.x = imo_pl_ext[n].j;
e_dpl(&pi, dlsds_ref, zs_ref, &ps);
imo_pl_dpl[n].i = (pi.x - ps.x)/dlsds_ref;
imo_pl_dpl[n].j = (pi.y - ps.y)/dlsds_ref;
imo_pl_dpl_ext[n].i = (pi.x - ps.x)/dlsds_ref;
imo_pl_dpl_ext[n].j = (pi.y - ps.y)/dlsds_ref;
}
}
else
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment