/* * mfeatures.c 7/3/84 Jim Piper * * compute features of chromosome objects */ #include #include "../h/wstruct.h" #include "../h/chromanal.h" #define RADIANS 57.296 #define AXISSMOOTH 5 #define PROFSMOOTH 1 /* * select which profiles should be used, one for band parameters, * one for centromere location. We already know that profile "0" (density) * is better than profile "-1" (mean density) or profile "-2" (max * density) for wdd features. */ int whichprof[2] = {0,2}; extern double xscale; extern double yscale; struct chromosome *chromfeatures(obj,basic) struct chromosome *obj; { struct chromosome *robj, *vertical(); struct object *mpol, *chromaxis2(), *prof[2]; struct object *cvh, *convhull(); struct histogramdomain *h; if (basic != 0) { basicfeatures(obj); robj = vertical(obj); } else { if (obj->plist->area == 0) basicfeatures(obj); cvh = convhull(obj); obj->plist->hullarea = hullarea(cvh, &obj->plist->hullperim); freeobj(cvh); robj = obj; } /* * Small objects can cause trouble e.g. in axis finding. * Don't do any further processing of anything that is * too small here. Other feature values will be zero, * but since this cannot be a chromosome this does not matter. */ if (robj->idom->lastln - robj->idom->line1 < 8) return(robj); mpol = chromaxis2(robj); /* * find density and moment profiles, and shorten mid-points polygon * and profile if tips lie outside object. Then chromosome length * is resulting profile length. */ multiprofiles(robj,mpol,prof,whichprof,2); shorten(mpol,prof[0],prof[1]); robj->plist->length = ((struct polygondomain *)(mpol->idom))->nvertices; /* set tip values to zero */ h = (struct histogramdomain *)prof[1]->idom; h->hv[0] = h->hv[h->npoints-1] = 0; /* * Smoothing possibly improves the power of wdd features, * so perhaps we should smooth more strongly here. */ histosmooth(prof[0], PROFSMOOTH); histosmooth(prof[1], PROFSMOOTH); /* * find centromere, measure features */ centromerefeatures(robj,mpol,prof[1]); /* * mass difference between top and bottom half of chromosome, * divided by chromosome area */ robj->plist->mdra = massdiff(robj); /* * profile density distribution features */ wddfeatures(robj,prof[0],prof[1]); freeobj(mpol); freeobj(prof[0]); freeobj(prof[1]); return (robj); } struct chromosome *vertical(obj) struct chromosome *obj; { struct object *cvh, *rect, *convhull(), *minwrect(); struct chromosome *robj, *spinsqueeze(); struct chromplist *plist; struct frect *rdom; float atan(),sin(),cos(); int i,j,k; cvh = convhull(obj); obj->plist->hullarea = hullarea(cvh, &obj->plist->hullperim); rect = minwrect(cvh); rdom = (struct frect *) rect->idom; /* * if a FIP digitisation then the aspect ratio, as * reflected in the values xscale and yscale, will have * lead to an incorrect angle (ignoring for the preesent the * effect on choice of minimum width direction) which * must be corrected prior to calling spinsqueeze, since * the angle given to spinsqueeze is the notional angle * of rotation after coordinate squeezing */ rdom->rangle = atan((xscale*sin(rdom->rangle)) / (yscale*cos(rdom->rangle))); robj = spinsqueeze(obj,rdom->rangle,xscale,yscale); /* * swap property list from original to rotated object * clear original pointer to permit correct "freeobj()" */ robj->plist = obj->plist; obj->plist = NULL; /* * fill in the property list */ plist = robj->plist; plist->rangle = (int) RADIANS*rdom->rangle; /* angle through which I have been rotated */ freeobj(cvh); freeobj(rect); return(robj); } basicfeatures (obj) struct chromosome *obj; { struct chromplist *plist; plist = obj->plist; plist->area = area(obj); plist->mass = mass(obj); } centromerefeatures(obj,mpol,prof) struct object *mpol, *prof; struct chromosome *obj; { struct ipoint *point, *ncentromere(); struct chromplist *plist; struct ivector v; /* * find centromere from second moment profile */ point = ncentromere(mpol,prof); plist = obj->plist; if (point != NULL) { plist->cx = point->k; plist->cy = point->l; plist->cangle = 0; /* since chromosome vertical */ plist->cindexa = cindexa(obj,point->l); plist->cindexm = cindexm(obj,point->l); plist->cindexl = cindexl(mpol,point->l); free(point); } else { plist->cindexa = -1; /* failed flag */ } } wddfeatures(obj,prof0,prof2) struct chromosome *obj; struct object *prof0, *prof2; { int i; struct chromplist *plist = obj->plist; struct object *profg, *gradprof(); profg = gradprof(prof0); wdd(prof0,plist->wdd); wdd(prof2,plist->mwdd); wdd(profg,plist->gwdd); freeobj(profg); plist->cvdd = cvdd(prof0); plist->nssd = nssd(prof0); /* * get wdd features the right way round w.r.t. area centromeric index */ if (plist->cindexa > 50) { for (i=0; i<=4; i += 2) { plist->wdd[i] = - plist->wdd[i]; plist->mwdd[i] = - plist->mwdd[i]; plist->gwdd[i] = - plist->gwdd[i]; } plist->mdra = - plist->mdra; } } /* * using density profile, shorten tips of mid points polygon, * density profile and moment profile * until within chromosome body (density at least 5% of mean). */ shorten(mpol,prof0,prof2) struct object *mpol, *prof0, *prof2; { register int i,l,*g,maxh; int nnz; register struct histogramdomain *hdom0, *hdom2; register struct polygondomain *p; hdom0 = (struct histogramdomain *)prof0->idom; hdom2 = (struct histogramdomain *)prof2->idom; p = (struct polygondomain *)mpol->idom; l = hdom0->npoints; g = hdom0->hv; maxh = nnz = 0; for (i=0; i 0) nnz++; g++; } maxh /= nnz; maxh /= 20; while (*--g < maxh) { hdom0->npoints--; hdom2->npoints--; p->nvertices--; } while (*hdom0->hv < maxh) { hdom0->hv++; hdom0->npoints--; hdom2->hv++; hdom2->npoints--; p->vtx++; p->nvertices--; } } xrange(pol) struct object *pol; { register struct polygondomain *pdom; register struct ivertex *vtx; register int i, minx, maxx; pdom = (struct polygondomain *)pol->idom; vtx = pdom->vtx; minx = vtx->vtX; maxx = minx; for (i=1; invertices; i++) { if ((++vtx)->vtX < minx) minx = vtx->vtX; else if (vtx->vtX > maxx) maxx = vtx->vtX; } return(maxx-minx); }